On a Possible Solution to the Tidal Realignment Problem for Hot Jupiters

We consider a star–planet system with masses and radii ${M}_{\star }$, ${R}_{\star }$, and ${M}_{{\rm{p}}}$, ${R}_{{\rm{p}}}$. We assume the star is rotating uniformly with frequency ${{\rm{\Omega }}}_{\star }$, and has a moment of inertia ${I}_{\star }={k}_{\star }{M}_{\star }{R}_{\star }^{2}$ with ${k}_{\star }\simeq 0.06$. The stellar angular momentum is denoted ${\boldsymbol{S}}=S\hat{{\boldsymbol{s}}}$ and the orbital angular momentum as ${\boldsymbol{L}}=L\hat{{\boldsymbol{l}}}$, with $\hat{{\boldsymbol{s}}}$ and $\hat{{\boldsymbol{l}}}$ unit vectors. We restrict our attention to circular, possibly inclined orbits, having in mind that the planet's orbital eccentricity has already been extensively damped, and the planet's rotation is synchronized with its orbit. We consider the tides raised by the planet on the host star, in concert with the magnetic braking torque on the host star.

The spin and orbital angular momenta each evolve due to the tidal torque and energy transfer. We adopt the tidal model of Lai (2012), which parameterizes the dissipation in terms of a separate phase lag for each tidal component. Below, we summarize the key features of this model.

Placing $\hat{{\boldsymbol{s}}}$ along the $\hat{{\boldsymbol{z}}}$-axis and $\hat{{\boldsymbol{l}}}$ in the x − z plane, the as-yet unspecified tidal torque and energy dissipation rate are ${\boldsymbol{T}}={T}_{x}\hat{{\boldsymbol{x}}}+{T}_{y}\hat{{\boldsymbol{y}}}+{T}_{z}\hat{{\boldsymbol{z}}}$ and $\dot{E}$, respectively. The $\hat{{\boldsymbol{x}}}$ and $\hat{{\boldsymbol{z}}}$ components introduce dissipation, leading to changes in orbital distance, spin rate, and obliquity. The $\hat{{\boldsymbol{y}}}$ component contributes only to precession, and does not otherwise affect the orbital and spin parameters; for this reason, we did not include the $\hat{{\boldsymbol{y}}}$ component in our calculations. The spin and orbit thus evolve under tides as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{\boldsymbol{S}}}{{dt}} & = & {T}_{x}\hat{{\boldsymbol{x}}}+{T}_{z}\hat{{\boldsymbol{z}}}\\ \displaystyle \frac{d{\boldsymbol{L}}}{{dt}} & = & -\displaystyle \frac{d{\boldsymbol{S}}}{{dt}}.\end{array}\end{eqnarray} \tag{ 1 }$$

In terms of semimajor axis a, spin frequency ${{\rm{\Omega }}}_{\star }$, and obliquity $\theta \equiv {\cos }^{-1}(\hat{{\boldsymbol{s}}}\cdot \hat{{\boldsymbol{l}}})$, the evolution equations take the form

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{\mathrm{Tide}}=-\displaystyle \frac{2{a}^{2}\dot{E}}{{{GM}}_{\star }{M}_{{\rm{p}}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{{\rm{\Omega }}}_{\star }}{{dt}}\right)}_{\mathrm{Tide}}=\displaystyle \frac{{T}_{z}}{{I}_{\star }},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d\theta }{{dt}}\right)}_{\mathrm{Tide}}=-\displaystyle \frac{{T}_{x}}{S}-\displaystyle \frac{{T}_{x}\cos \theta }{L}+\displaystyle \frac{{T}_{z}\sin \theta }{L}.\end{eqnarray} \tag{ 4 }$$

The tidal torque components and energy dissipation rate are obtained by expanding the tidal potential in spherical harmonics and transforming to the rotating frame of the host star to analyze the fluid response (see Sections 2.1–2.3 of Lai 2012). For a circular, inclined orbit, there are seven independent Fourier components of the tidal potential, operating at frequencies

$$\begin{eqnarray}&&{\omega }_{{mm}^{\prime} }=m^{\prime} {\rm{\Omega }}-m{{\rm{\Omega }}}_{\star }.\end{eqnarray} \tag{ 5 }$$

For a derivation of the tidal components for circular, spin–orbit misaligned systems, see Barker & Ogilvie (2009). The physically distinct tidal components that lead to dissipation are $(m,m^{\prime} )=(0,2)$, $(\pm 1,2)$, $(\pm 2,2)$, $(1,0)$, and $(2,0)$. For each component, the dissipation is specified by a phase shift ${{\rm{\Delta }}}_{{mm}^{\prime} }$ and a Love coefficient ${\kappa }_{{mm}^{\prime} }$. Since only the product of these two parameters is relevant to our calculations, we define ${\rm{\Delta }}{{\prime} }_{{mm}^{\prime} }\equiv {{\rm{\Delta }}}_{{mm}^{\prime} }{\kappa }_{{mm}^{\prime} }$. The torque components T_x and T_z and the energy dissipation rate $\dot{E}$ are calculated as a sum over contributions from all of the tidal components, and can be written explicitly as (see Lai 2012, Equations (27), (28), (35))

$$\begin{eqnarray}&&\begin{array}{l}{T}_{x}=\displaystyle \frac{3\pi }{20}{T}_{0}\left[\displaystyle \frac{1}{2}{s}_{\theta }{\left(1+{c}_{\theta }\right)}^{3}{\rm{\Delta }}{{\prime} }_{22}+{s}_{\theta }{\left(1+{c}_{\theta }\right)}^{2}(2-{c}_{\theta }){\rm{\Delta }}{{\prime} }_{12}\right.\\ \,+\,3{s}_{\theta }^{3}{\rm{\Delta }}{{\prime} }_{02}+{s}_{\theta }{\left(1-{c}_{\theta }\right)}^{2}(2+{c}_{\theta }){\rm{\Delta }}{{\prime} }_{-12}\\ \,+\,\left.\displaystyle \frac{1}{2}{s}_{\theta }{\left(1-{c}_{\theta }\right)}^{3}{\rm{\Delta }}{{\prime} }_{-22}\right]\\ \,-\,\displaystyle \frac{3\pi }{5}{T}_{0}\left[\displaystyle \frac{1}{2}{s}_{\theta }^{3}{c}_{\theta }{\rm{\Delta }}{{\prime} }_{20}+{s}_{\theta }{c}_{\theta }^{3}{\rm{\Delta }}{{\prime} }_{10}\right],\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\begin{array}{l}{T}_{z}=\displaystyle \frac{3\pi }{20}{T}_{0}\left[\displaystyle \frac{1}{2}{\left(1+{c}_{\theta }\right)}^{4}{\rm{\Delta }}{{\prime} }_{22}+{s}_{\theta }^{2}{\left(1+{c}_{\theta }\right)}^{2}{\rm{\Delta }}{{\prime} }_{12}\right.\\ \,\left.-\,{s}_{\theta }^{2}{\left(1-{c}_{\theta }\right)}^{2}{\rm{\Delta }}{{\prime} }_{-12}-\displaystyle \frac{1}{2}{\left(1-{c}_{\theta }\right)}^{4}{\rm{\Delta }}{{\prime} }_{-22}\right]\\ \,+\,\displaystyle \frac{3\pi }{5}{T}_{0}\left[\displaystyle \frac{1}{2}{s}_{\theta }^{4}{\rm{\Delta }}{{\prime} }_{20}+{s}_{\theta }^{2}{c}_{\theta }^{2}{\rm{\Delta }}{{\prime} }_{10}\right],\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\begin{array}{l}\dot{E}=\displaystyle \frac{3\pi }{20}{T}_{0}{\rm{\Omega }}\left[\displaystyle \frac{1}{2}{\left(1+{c}_{\theta }\right)}^{4}{\rm{\Delta }}{{\prime} }_{22}+2{s}_{\theta }^{2}{\left(1+{c}_{\theta }\right)}^{2}{\rm{\Delta }}{{\prime} }_{12}\right.\\ \,\left.+\,3{s}_{\theta }^{4}{\rm{\Delta }}{{\prime} }_{02}+2{s}_{\theta }^{2}{\left(1-{c}_{\theta }\right)}^{2}{\rm{\Delta }}{{\prime} }_{-12}+\displaystyle \frac{1}{2}{\left(1-{c}_{\theta }\right)}^{4}{\rm{\Delta }}{{\prime} }_{-22}\right],\end{array}\end{eqnarray} \tag{ 8 }$$

where ${c}_{\theta }=\cos \theta$, ${s}_{\theta }=\sin \theta$, and

$$\begin{eqnarray}&&{T}_{0}=\displaystyle \frac{{{GM}}_{{\rm{p}}}^{2}{R}_{\star }^{5}}{{a}^{6}}.\end{eqnarray} \tag{ 9 }$$

This tidal model thus allows for general specification of each phase lag ${\rm{\Delta }}{{\prime} }_{{mm}^{\prime} }$, which we relate to the modified quality factor $Q{{\prime} }_{{mm}^{\prime} }$ via $Q{{\prime} }_{{mm}^{\prime} }=1/| {\rm{\Delta }}{{\prime} }_{{mm}^{\prime} }|$. ⁴ By modeling the tidal spin-up of hot Jupiter hosts and comparing to observations, Penev et al. (2018) found evidence that the $m=m^{\prime} =2$ component of the stellar modified tidal quality factor exhibits the frequency dependence

$$\begin{eqnarray}&&{Q}_{22}^{{\prime} }\propto {\omega }_{22}^{3}.\end{eqnarray} \tag{ 10 }$$

We decided to investigate a more general power-law dependence, ${Q}_{{mm}^{\prime} }^{{\prime} }\propto {\omega }_{{mm}^{\prime} }^{p}$, with a single value of p that applies to all of the tidal components. We thus adopt phase shifts of the form

$$\begin{eqnarray}&&{\rm{\Delta }}{{\prime} }_{{mm}^{\prime} }\propto {\omega }_{{mm}^{\prime} }^{-p}.\end{eqnarray} \tag{ 11 }$$

Most of the results we will describe are for the case p = 3, motivated by the work of Penev et al. (2018), but we also investigated some other possible values of p. We note that $p=-1$ corresponds to a constant time lag, while p = 0 corresponds to a constant phase lag. As the frequency decreases, the scaling in Equation (11) must break down at some point, and when the tidal frequency approaches zero, the dissipation should vanish. We therefore introduce a softening parameter and impose a maximum (in magnitude) phase shift ${\rm{\Delta }}{{\prime} }_{\max }$, so that the phase shift for the ${mm}^{\prime}$ component takes the form

$$\begin{eqnarray}&&{\rm{\Delta }}{{\prime} }_{{mm}^{\prime} }=\mathrm{sign}({\tilde{\omega }}_{{mm}^{\prime} })\min \left[\displaystyle \frac{{\rm{\Delta }}{{\prime} }_{0}| {\tilde{\omega }}_{{mm}^{\prime} }| }{{\left(| {\tilde{\omega }}_{{mm}^{\prime} }| +\epsilon \right)}^{p+1}},{\rm{\Delta }}{{\prime} }_{\max }\right],\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{\tilde{\omega }}_{{mm}^{\prime} }\equiv \displaystyle \frac{{\omega }_{{mm}^{\prime} }}{2\pi \ {\mathrm{days}}^{-1}},\end{eqnarray} \tag{ 13 }$$

and where the scaling constants ${\rm{\Delta }}{{\prime} }_{0},{\rm{\Delta }}{{\prime} }_{\max }$ are positive quantities. For $| {\omega }_{{mm}^{\prime} }| \gg \epsilon$, Equation (12) yields ${\rm{\Delta }}{{\prime} }_{{mm}^{\prime} }\propto {\omega }_{{mm}^{\prime} }^{-p}$, and for $| {\omega }_{{mm}^{\prime} }| \ll \epsilon$, ${\rm{\Delta }}{{\prime} }_{{mm}^{\prime} }\to 0$. The scaling found by Penev et al. (2018) was calibrated using hot Jupiter systems with tidal forcing periods as short as 0.5 days. This implies that the softening parameter should satisfy $\epsilon \ll 1$. Somewhat arbitrarily, we chose $\epsilon =0.05$. We experimented with values of a factor of two larger and smaller than 0.05, and found no qualitative differences in the results of individual integrations, and no difference at all in the statistical results.

With these choices, there is a minimum possible value for the modified quality factor $Q{{\prime} }_{\min }=1/{\rm{\Delta }}{{\prime} }_{\max }$ and an overall scale $Q{{\prime} }_{0}=1/{\rm{\Delta }}{{\prime} }_{0}$. We chose $Q{{\prime} }_{\min }={10}^{5}$, motivated by observations of eclipsing binaries (Milliman et al. 2014). The adopted scaling of ${Q}_{{mm}^{\prime} }$ (from Equation (12)) is depicted in Figure 1. The results of Penev et al. (2018) indicate that $Q{{\prime} }_{0}\sim {10}^{6}$. In Section 3 we explore the implications of varying ${Q}_{0}^{\prime}$, and in Section 4 we adjust the exact value of $Q{{\prime} }_{0}$ system by system in order to reproduce the observed spin and orbital periods.

Zoom In Zoom Out Reset image size

Note that in this model, $Q{{\prime} }_{{mm}^{\prime} }$ loses its frequency dependence and becomes a constant for low enough tidal frequencies:

$$\begin{eqnarray}&&| {\omega }_{{mm}^{\prime} }| \leqslant \displaystyle \frac{2\pi }{\mathrm{days}}{\left(\displaystyle \frac{Q{{\prime} }_{\min }}{Q{{\prime} }_{0}}\right)}^{1/p}.\end{eqnarray} \tag{ 14 }$$

For the $(1,0)$ component, this corresponds to a spin frequency

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\star }\leqslant \displaystyle \frac{2\pi }{\mathrm{days}}{\left(\displaystyle \frac{Q{{\prime} }_{\min }}{Q{{\prime} }_{0}}\right)}^{1/p}.\end{eqnarray} \tag{ 15 }$$

Along with the tidal torques, we include the spin-down torque on the host star due to magnetic braking (${T}_{\mathrm{MB}}$), via

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{{\rm{\Omega }}}_{\star }}{{dt}}\right)}_{\mathrm{MB}}=-\alpha \min {\left({{\rm{\Omega }}}_{\mathrm{sat}},{{\rm{\Omega }}}_{\star }\right)}^{2}\,{{\rm{\Omega }}}_{\star }=-\displaystyle \frac{{T}_{\mathrm{MB}}}{{I}_{\star }},\end{eqnarray} \tag{ 16 }$$

(see, e.g., Skumanich 1972; Kawaler 1988; Gallet & Bouvier 2013), where ${{\rm{\Omega }}}_{\mathrm{sat}}$ is the frequency at which the stellar dynamo saturates, and

$$\begin{eqnarray}&&\alpha \simeq 1.5\times {10}^{-14}\mathrm{yr}\,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-3/2}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-3/2}.\end{eqnarray} \tag{ 17 }$$

The mass and radius scaling in Equation (17) is motivated by Equation (3) of Amard et al. (2016). The coefficient, taken from Barker & Ogilvie (2009), reproduces the spin period of the Sun at its current age. Although the dependence on stellar mass and radius is taken from spin-down models that may not be completely realistic, we will ultimately restrict our attention to a relatively narrow range of stellar effective temperatures to help alleviate this issue. We set ${{\rm{\Omega }}}_{\mathrm{sat}}=10\,{{\rm{\Omega }}}_{\odot }$ throughout this paper.

In summary, the spin and orbit evolve under the influence of both tides and magnetic braking according to Equations (2)–(4) and (16), using Equation (12) to specify the dissipation for each tidal component. In practice, we evolve the system in terms of the Cartesian components of the angular momentum vectors in an inertial frame, rather than in terms of orbital elements.

If at any point the orbital distance becomes less than the tidal disruption radius ${R}_{\mathrm{tide}}$, with

$$\begin{eqnarray}&&{R}_{\mathrm{tide}}=\eta {R}_{{\rm{p}}}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{{\rm{p}}}}\right)}^{1/3},\end{eqnarray} \tag{ 18 }$$

we terminate the integration and consider the planet to be destroyed. The value of η depends on the internal constitution of the planet, which is of course unknown. Estimates based on numerical simulations are typically between 2 and 3 (e.g., Guillochon et al. 2011). Unless otherwise specified, we set $\eta =2.5$ for the remainder of this paper.

Here we discuss the general features of the tidal model, and present several illustrative examples of tidal evolution in spin–orbit misaligned systems.

As discussed by Lai (2012), the $(m,m^{\prime} )=(1,0)$ and $(2,0)$ tidal components do not transfer energy, because they are static in the inertial frame. These components thereby realign obliquities without inducing orbital decay. For typical hot Jupiter systems, the orbital frequency Ω usually exceeds the spin frequency ${{\rm{\Omega }}}_{\star }$ by a factor of at least a few. Consequently, the $(m,m^{\prime} )=(1,0)$ component typically has the lowest forcing frequency of all of the tidal components. Under the assumption that ${Q}_{{mm}^{\prime} }^{{\prime} }\propto {\omega }_{{mm}^{\prime} }^{p}$ (with $p\gt 0$), Q₁₀ can be much smaller than the other ${Q}_{{mm}^{\prime} }^{{\prime} }$. This is why the model is capable of reducing the obliquity without destroying the planet.

It is well known that the $(m,m^{\prime} )=(1,0)$ component alone may drive the obliquity toward 0°, 90°, or 180° depending on the initial obliquity (Lai 2012; Rogers & Lin 2013). However, considering the combined effects of all of the tidal components, the perpendicular and antiparallel obliquity states are at most temporarily stable (Xue et al. 2014; Li & Winn 2016), although the system may reside in such states for long periods of time.

To illustrate this phenomenon using the specific tidal model adopted in this paper, we calculate $\dot{\theta }$ as a function of θ, using Equation (4). We denote by ${\theta }^{* }$ the values of θ for which $\dot{\theta }=0$ (i.e., the "fixed points" ⁵ ). The sign of $\dot{\theta }$ in the neighborhood of ${\theta }^{* }$ is related to the stability of the fixed point. We denote by ${\theta }_{-}^{* }$ and ${\theta }_{+}^{* }$ the values of θ just less than ${\theta }^{* }$, and just greater than ${\theta }^{* }$, respectively. If $\dot{\theta }({\theta }_{-}^{* })\gt 0$ and $\dot{\theta }({\theta }_{+}^{* })\lt 0$, the system may be driven toward ${\theta }^{* }$. In contrast, if $\dot{\theta }({\theta }_{-}^{* })\lt 0$ and $\dot{\theta }({\theta }_{+}^{* })\gt 0$, the system will be repelled from ${\theta }^{* }$. Finally, if the sign of $\dot{\theta }$ does not change from ${\theta }_{-}^{* }$ to ${\theta }_{+}^{* }$, ${\theta }^{* }$ is a saddle point; it is stable only when approaching from one direction.

Figure 2 shows $\dot{\theta }$ as a function of θ, with the tidal model parameters set to p = 3 and ${Q}_{0}^{\prime} ={10}^{6}$. We fixed the orbital period and vary the spin period. We compare results for the $(1,0)$ component in isolation, and the combined effects of all components. Although 90° is a fixed point when considering the $(1,0)$ component alone, it is not a fixed point when including all of the other tidal components. Figure 3 illustrates the behavior of ${\theta }^{* }$ as a function of ${P}_{\star }$. The number of ${\theta }^{* }$ and their stability depends upon the stellar spin period (with all other system parameters held fixed). The spin–orbit aligned state is always stable, while the anti-aligned state may be stable or unstable depending upon the spin period. This indicates that the obliquity may spend periods of time locked at 180° (as we show later via numerical integrations).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even if the obliquity reaches one of the nonzero ${\theta }^{* }$ states, this does not necessarily imply that the system will remain in such a state for long, since ${\theta }^{* }$ and its stability depend on the spin and orbital period, which both continue to evolve when $\theta ={\theta }^{* }$. For sufficiently large obliquities, ${\dot{{\rm{\Omega }}}}_{\star }\lt 0$, causing the star to spin down even as the orbital period shrinks. When $\theta =180^\circ$, the stellar spin magnitude evolves according to

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Omega }}}_{\star }}{{dt}}=-\displaystyle \frac{6\pi {T}_{0}}{5{I}_{\star }}{\rm{\Delta }}{{\prime} }_{-22}-\displaystyle \frac{{T}_{\mathrm{MB}}}{{I}_{\star }}\lt 0.\end{eqnarray} \tag{ 19 }$$

Since ${{\rm{\Delta }}}_{-22}\gt 0$, the star's rotation must slow down in the anti-aligned state.

Figure 4 depicts the spin and orbital evolution as a function of time for a range of different initial obliquities, spanning both initially prograde and retrograde orientations. The initial orbital period and arrival time for the planet are held fixed at 3 days and 2 Gyr, respectively, which is consistent with a hot Jupiter that previously underwent high-eccentricity migration. We assume an initial stellar spin period of 3 days, and evolve the spin magnitude due to magnetic braking alone up to the time of the planet's arrival, and then evolve the system under the influence of both tides and magnetic braking. In these examples, the orbital evolution is similar in all cases, but the spin and obliquity evolution differs for the initially prograde and retrograde systems. For the initially prograde systems, the spin frequency briefly decreases due to magnetic braking, before the (positive) tidal torque overcomes the negative spin-down torque, thus spinning up the host star. Meanwhile, the obliquity is driven toward alignment. For the initially retrograde systems, the star is de-spun while the obliquity is temporarily driven toward anti-alignment. Then, the star halts momentarily, before spinning up in the opposite direction. As a result, the obliquity flips from 180° to 0°. Systems with initially retrograde obliquities experience slower orbital decay than those with initially prograde obliquities. This is in contrast with previous studies employing a constant time lag (in which $Q^{\prime}$ decreases with forcing frequency; e.g., Barker & Ogilvie 2009). This behavior in our model is a consequence of the rise in $Q^{\prime}$ with forcing frequency.

Zoom In Zoom Out Reset image size

In all of the examples depicted in Figure 4, the star spends the majority of the time in a spin–orbit aligned state. The same initial conditions, evolved under a constant phase-lag formalism (rather than frequency-dependent), fail to realign before the planet reaches the tidal disruption radius.

In this section we calculate possible tidal histories for real hot Jupiters orbiting cool host stars. Since the rotational evolution of the star depends on the extent of the surface convective region, and because magnetic braking is best understood for solar-like stars, we restrict our attention to systems with temperatures in the range 5500–6000 K. In order to construct as homogeneous a sample as possible, we adopt the effective temperatures and metallicities from SWEET-Cat (Santos et al. 2013), which is regularly updated with new planet discoveries. We then cross-match with planets from the NASA Exoplanet Archive and select transiting planets with masses larger than $0.5{M}_{{\rm{J}}}$ and orbital periods less than 5 days. The reason we focus on transiting planets is that their true masses are known, as opposed to Doppler planets for which only the minimum mass is known. The upper limit on the orbital period was chosen to ensure that the sample largely contains systems susceptible to tidal evolution. We reject any systems with an eccentricity above 0.05, to avoid introducing additional parameters in the model to account for tides raised on the planet by the star. We note that many of the planets in our sample have orbital eccentricities that are poorly constrained by the available data. However, observations have shown that highly eccentric systems around cool stars are rare when the period is shorter than 5 days. For each system, we perform a literature search to determine the stellar rotation period, obtained either from photometric variability (when available), or from the stellar $v\sin i$ and radius (assuming spin–orbit alignment). We discard systems lacking any information on the spin period.

We use the Isochrones package to calculate stellar ages (Morton 2015). Isochrones provides Bayesian stellar age and mass estimates by interpolating over stellar isochrones obtained from the MESA (Modules for Experiments in Stellar Astrophysics) Isochrones and Stellar Tracks (MIST) models (Choi et al. 2016; Dotter 2016). Isochrones accepts a variety of user inputs, including any combination of effective temperature, metallicity, surface gravity, stellar density, V-band extinction, parallax, and photometry. We input the stellar effective temperature and metallicity from SWEET-Cat, along with the stellar density obtained from the transit light-curve parameters ${P}_{\mathrm{orb}}$ and ${\text{}}a/{\text{}}{R}_{\star }$, via

$$\begin{eqnarray}&&{\rho }_{\star }=\displaystyle \frac{3\pi }{{{GP}}_{\mathrm{orb}}^{2}}{\left(\displaystyle \frac{a}{{R}_{\star }}\right)}^{3}.\end{eqnarray} \tag{ 20 }$$

We also input the Gaia parallax and photometry ($G,{BP},{RP}$), and the Two Micron All Sky Survey (2MASS) ${\mathrm{JHK}}_{s}$ and Wide-field Infrared Survey Explorer (WISE) 123 bands. We use the 3D galactic extinction model MWDUST (Bovy et al. 2016) to estimate the extinction at the distance implied by the Gaia parallax, and impose a flat prior with maximum ${A}_{{\rm{V}},\max }+0.1\ \mathrm{mag}$, where ${A}_{{\rm{V}},\max }$ is the maximum extinction along the entire line of sight. We allow for additional extinction of 0.1 mag beyond the maximum value to account for uncertainties in the galactic extinction model. See Appendix for additional details.

After discarding several systems with unreliable parameters or output ages (see Appendix), our final sample consists of 46 hot Jupiters orbiting cool host stars. This sample is a subset of the observed planets depicted in Figure 5, since here we restrict our attention to transiting planets.

Next, we present tidal histories for the sample of hot Jupiters orbiting cool stars, adopting the fiducial scaling of the quality factor with forcing frequency (p = 3). Figure 8 depicts possible tidal histories for the HATS-18 system (Penev et al. 2016), an illustrative case. We assumed the age to be 4.2 Gyr, and sampled a range of initial obliquities and arrival times for the planet. We then iteratively determined the values of the initial orbital period and ${Q}_{0}^{\prime}$ that reproduce the observed spin and orbital periods at the adopted system age. For all of the initial conditions explored, the obliquity has realigned at the adopted system age.

Zoom In Zoom Out Reset image size

Similar to the example depicted in Figure 8 for HATS-18, we calculate possible tidal histories for each of the 46 observed hot Jupiters, sampling the stellar age directly from the Isochrones age posteriors. We conduct two separate experiments, corresponding to early arrival of the planet due to disk migration, and late arrival due to high-eccentricity migration. We refer to the two experiments as EARLY and LATE, respectively. The procedure is as follows:

• 1.  
  Draw a random sample ${t}_{\mathrm{age}}$ from the Isochrones age posterior.
• 2.  
  Sample an arrival time ${t}_{\mathrm{arr}}$ for the planet (3 Myr for the EARLY experiment, and uniformly in the range between $50\mathrm{Myr}$ and $0.9{t}_{\mathrm{age}}$ for the LATE experiment).
• 3.  
  Sample the initial obliquity from an isotropic distribution (uniform in $\cos \theta$).
• 4.  
  Sample an initial stellar T-Tauri spin period uniformly between 1–10 days, and evolve the spin period due to magnetic braking up to the time of the planet's arrival.
• 5.  
  Evolve the system due to both tides and magnetic braking up to ${t}_{\mathrm{age}}$, iteratively adjusting the initial orbital period and the value of ${Q}_{0}^{\prime}$ to reproduce the observed spin and orbital periods at ${t}_{\mathrm{age}}$ (within the observational uncertainties). We impose a minimum uncertainty of $10 \%$ for the spin period, regardless of the quoted observational uncertainty.

For each system, this procedure results in a distribution of stellar obliquities at the current system age. Figure 9 illustrates the final obliquity distributions from our simulations of each of the 46 systems, as a function of the observed orbital period. The black error bars indicate the 16th and 84th percentiles of the obliquity distribution, while red error bars indicate the observed sky-projected obliquities for the systems for which obliquity measurements are available. Our calculations generally imply that systems with orbital periods less than 2.5 days have undergone substantial obliquity damping, and have final obliquities smaller than 1°.

Zoom In Zoom Out Reset image size

For the majority of systems, the current obliquity is predicted to be close to zero, indicating efficient realignment. However, for several stars that are observed to be slowly rotating, our calculations predict retrograde obliquities. To demonstrate this, we compare the observed spin period of the star to the expected spin period of an isolated star at the same age (for which the spin period evolves due to magnetic braking alone). The spin period as a function of age for an isolated star is approximately

$$\begin{eqnarray}&&{P}_{\star ,\mathrm{mb}}\simeq {P}_{\star ,0}{\left(1+\displaystyle \frac{8{\pi }^{2}\alpha {t}_{\mathrm{age}}}{{P}_{\star ,0}^{2}}\right)}^{1/2},\end{eqnarray} \tag{ 21 }$$

where ${P}_{\star ,0}$ is the spin period during the T-Tauri phase, and α is given by Equation (17). Thus, for each system we can compare ${P}_{\star ,\mathrm{obs}}$ and ${P}_{\star ,\mathrm{mb}}$. Figure 10 shows the final obliquity distribution versus ${P}_{\star ,\mathrm{obs}}/{P}_{\star ,\mathrm{mb}}$, where we have used the median value of ${P}_{\star ,\mathrm{mb}}$ obtained from our initial spin periods and age posteriors. For the slowly rotating stars (${P}_{\star ,\mathrm{obs}}/{P}_{\star ,\mathrm{mb}}\gtrsim 1$), there is often a preference for initially retrograde obliquities. This preference may be understood by noticing that tides in sufficiently misaligned configurations act to de-spin the star (see, e.g., Equation (19)). However, we stop short of predicting that such stars are truly retrograde, because all of the rotation periods for the anomalously slowly rotating stars were derived from the stellar $v\sin i$ under the assumption $\sin i=1$, rather than from the less ambiguous method of detecting a periodicity in the photometric variability. As a result, another possible interpretation is that the spin period has been overestimated because $\sin i\ne 1$, i.e., the obliquity is currently large. In this scenario, repeating the tidal calculations with the (unknown) true spin period would possibly allow initially prograde configurations to reproduce the system properties, by decreasing or eliminating the need for tidal de-spinning. Given our present knowledge, we cannot distinguish between these two possibilities.

Zoom In Zoom Out Reset image size

We also calculate posteriors for the remaining lifetime, defined as the time until the planet crosses the tidal disruption radius ${R}_{\mathrm{tide}}$, given by Equation (18) (with $\eta =2.5$). Each set of initial conditions (with ${t}_{\mathrm{arr}}$, ${\theta }_{0}$, and ${t}_{\mathrm{age}}$ randomly sampled) is integrated forward in time up to 5 Gyr, using the values of ${Q}_{0}^{\prime}$ and the initial orbital period that were previously determined to reproduce the observed spin and orbital periods. If the planet crosses ${R}_{\mathrm{tide}}$, we terminate the integration and record the remaining lifetime. If the planet remains exterior to ${R}_{\mathrm{tide}}$ over the entire 5 Gyr integration, we set the remaining lifetime to 5 Gyr. Figure 11 depicts the results as a function of orbital period. Most of the planets are not in immediate danger of being devoured by their host stars. The typical remaining lifetime is at least hundreds of millions years. An exception to this trend is WASP-19b (the shortest-period planet in our sample), which our calculations predict will cross ${R}_{\mathrm{tide}}$ within 10–20 Myr. In contrast, the only other planet in our sample with a similar orbital period (HATS-18b) has a remaining lifetime of ∼500 Myr. The key difference between these two planets is that WASP-19b is located extremely close to ${R}_{\mathrm{tide}}$, according to our definition. This illustrates that the planet's fate depends sensitively on the chosen disruption criterion, through the value of η (see Equation (18)).

Zoom In Zoom Out Reset image size

To explore this issue further, we perform additional calculations for the WASP-19 and HATS-18 systems, varying η between 2 and a maximum possible value

$$\begin{eqnarray}&&{\eta }_{\max }=\displaystyle \frac{a}{{R}_{{\rm{p}}}}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\right)}^{1/3}.\end{eqnarray} \tag{ 22 }$$

For each value of η, we integrate the system forward in time and calculate the remaining lifetime, assuming spin–orbit alignment and varying ${Q}_{0}^{\prime}$. The results are shown in Figure 12. For a given value of ${Q}_{0}^{\prime}$, reducing η by a small amount prolongs the remaining lifetime for both systems. We also compare the results for our fiducial frequency-dependent quality factor (with p = 3, solid curves) to a constant quality factor (p = 0, dashed curves). Figure 11 reveals that the frequency-dependent quality factor prolongs the remaining lifetime by more than a factor of 10.

Zoom In Zoom Out Reset image size

Tidal realignment remains a promising explanation for at least some aspects of the well-known pattern relating measurements of stellar obliquity and the effective temperature of the host star (Winn et al. 2010). This paper describes our investigations of the tidal evolution of spin–orbit misaligned systems with a frequency-dependent stellar tidal quality factor. Our model is motivated by the recent study of stellar rotation in hot Jupiter systems by Penev et al. (2018), who found evidence that the tidal quality factor increases strongly with forcing frequency (assuming spin–orbit alignment). For spin–orbit misaligned systems, there are multiple independent tidal Fourier components, each operating at a different forcing frequency (a linear combination of the orbital and spin frequencies) and each associated with a potentially different quality factor (Lai 2012; Ogilvie 2014). We have implemented a frequency-dependent quality factor for each tidal component. The tidal component operating at the stellar spin frequency affects the obliquity without introducing orbital decay. Since the strong frequency dependence causes the components with the lowest forcing frequencies to have the smallest quality factors, this model efficiently re-aligns obliquities before the planetary orbit catastrophically decays.

Including the effects of both tides and magnetic braking, we have explored this model for simulated hot Jupiters with a range of properties, as well as for a sample of real hot Jupiters orbiting cool host stars. Using Gaia DR2 data, we have derived stellar ages for the hot Jupiter host stars, and explored possible tidal histories, starting with a broad range of initial stellar obliquities and two very different possibilities for the arrival times of the planet. Our calculations confirm that large obliquities are efficiently damped. The model predicts extremely low obliquities (typically much less than 1°) for hot Jupiters orbiting cool stars within 2.4 days (Figure 9).

The most important assumption in our model is that the tidal quality factor rises with forcing frequency over the range of frequencies typical in hot Jupiter systems. This assumption was inspired by the earlier findings of Penev et al. (2018), who investigated the possible tidal evolution histories of a large collection of hot Jupiter systems. A nearly contemporaneous population-level study of tidal dissipation in hot Jupiter systems was performed by Collier Cameron & Jardine (2018), who did not report evidence for a frequency-dependent tidal quality factor. This makes it important to compare the two studies to see if there is any contradiction.

The Bayesian hierarchical analysis of Collier Cameron & Jardine (2018) was based on the observed distribution of planetary orbital periods, and assumed an initially log-normal semimajor axis distribution. Those authors chose to analyze all of the planets in an online catalog of transiting planets called TEP-Cat. ⁶ They divided the sample into two groups: systems with $0.5\lt {P}_{\mathrm{orb}}/{P}_{\star }\lt 2$ (for which inertial waves can be excited), and systems with relative spin and orbital periods outside of the aforementioned range (for which they assumed the dissipation of the equilibrium tide is the dominant mechanism). They inferred a quality factor $Q^{\prime} \sim {10}^{7}$ for the inertial-wave systems and $Q^{\prime} \sim {10}^{8}$ for the equilibrium-tide systems. In contrast, Penev et al. (2018) restricted their sample to have orbital periods shorter than 3 days, and host stars cooler than 6200 K, narrower ranges than in the sample of Collier Cameron & Jardine (2018). As a result of the short orbital period cutoff, none of the Penev et al. (2018) systems are unambiguously in the inertial-wave regime. The inertial-wave sample of Collier Cameron & Jardine (2018) consists of planets with longer tidal forcing periods compared to the equilibrium-tide sample (so that the condition $0.5\lt {P}_{\mathrm{orb}}/{P}_{\star }\lt 2$ is satisfied). As a result, the findings of Collier Cameron & Jardine (2018) are broadly consistent with the finding of Penev et al. (2018) that longer-period systems have lower values of $Q{{\prime} }_{\star }$.

The results of this paper are subject to some uncertainties and caveats. Notably, we adopted the empirically motivated scaling for the stellar quality factor for each tidal component, as found by Penev et al. (2018) for spin–orbit aligned systems. For systems with spin–orbit misalignment, it is far from obvious whether the dissipation introduced by the various tidal components should be expected to scale in the same way. In the face of this uncertainty, we adopted the same quality factor for each tidal component, representing the simplest possible assumption. The theoretically expected dependence of $Q$ on tidal frequency remains highly uncertain and depends on the dissipation mechanism under consideration.

Tidal dissipation in solar-type stars may arise from the equilibrium tide or from dynamical tides (the excitation and damping of various kinds of waves inside the star). In the case of equilibrium tides, dissipation is thought to arise from damping of turbulent viscosity in convection zones (Zahn 1966). In this theory, the scaling of the quality factor with tidal frequency depends on the eddy viscosity. If the eddy viscosity is independent of forcing frequency, the quality factor scales as $Q\propto {\omega }^{-1}$ (constant lag time). However, for tidal frequencies larger than the convective turnover frequency ${\omega }_{{\rm{c}}}$, the eddy viscosity is thought to be reduced below this expectation. The exact form of the reduction factor has been disputed for decades (Zahn 1966, 1977; Goldreich & Nicholson 1977; Goodman & Oh 1997). Zahn (1966) argued for a reduction factor of $\omega /{\omega }_{c}$, which implies that the quality factor should scale as $Q\propto {\omega }^{0}$ (i.e., constant Q). In contrast, Goldreich & Nicholson (1977) argued for a reduction factor of ${(\omega /{\omega }_{{\rm{c}}})}^{2}$, giving $Q\propto \omega$. Since then, numerical simulations have been undertaken to determine the frequency dependence of the quality factor (Penev et al. 2009; Ogilvie & Lesur 2012), with a recent study finding support for the quadratic reduction factor (Duguid et al. 2020). This is qualitatively consistent with the models explored in this paper, especially the case of our p = 1 calculations. However, although the scaling with frequency is in general agreement with our models, the numerical simulations lead to quality factors that are larger in absolute terms, which would lead to negligible obliquity and spin evolution. Various uncertainties remain in the models, however, and further studies are needed.

Regarding dynamical tides, a leading contender for the orbital evolution in hot Jupiter systems is the excitation and damping of internal gravity waves. For the solar-type stars considered in this paper, these waves are excited at the core-envelope boundary, and propagate inwards in the radiative zone. The degree of dissipation depends strongly on whether or not these waves break at the center (Goodman & Dickson 1998; Barker & Ogilvie 2010), which itself depends on the planetary and stellar masses, as well as the stellar age (see Figures 8 and 9 of Barker (2020)). If the waves are fully damped, this mechanism predicts the scaling $Q\propto {\omega }^{-8/3}$, which is very different from the form adopted in this paper. Since the critical mass for wave breaking decreases strongly with age, the results of Barker (2020) predict that the quality factor should decrease with age, eventually leading to planetary engulfment. However, it is not clear what fraction of the observed hot Jupiters are currently in this wave-breaking regime, especially given the large uncertainties in stellar ages. Further studies and more precise stellar ages are needed to clarify this issue.

The excitation and damping of inertial waves is another contender for driving the orbital and obliquity evolution in hot Jupiter systems, particularly during the pre-main-sequence stage of stellar evolution. The theory of inertial waves predicts a strong and complicated frequency dependence of the quality factor. Averaged over frequency, the theory predicts $Q\propto {{\rm{\Omega }}}_{\star }^{-2}$ (e.g., Ogilvie & Lin 2007). Mathis (2015) has argued for greatly enhanced levels of dissipation due to inertial waves during the pre-main-sequence phase, based on a frequency-averaged two-layer model derived by Ogilvie (2013). Assuming spin–orbit alignment, and coupling the Ogilvie (2013) model to stellar and orbital evolution, Bolmont & Mathis (2016) investigated the tidal evolution of various star–planet systems, with subsequent calculations performed by Heller (2019) in concert with disk torques. These studies have recently been improved upon by Barker (2020), accounting for the realistic structure (density profile) of the star. For spin–orbit misaligned systems, the complexity of the calculation increases dramatically due to the multiple Fourier components at work. Accounting for the changes in dissipation due to stellar evolution for spin–orbit misaligned systems is therefore beyond the scope of this paper (but see Lin & Ogilvie 2017, for an analysis of the dissipation introduced by the $(m,m^{\prime} )=(1,0)$ component in isolation). In any case, the previous works, while intriguing, are subject to their own uncertainties and caveats. If we had taken into account the possible increase in dissipation during the pre-main-sequence phase, the results of our simulations involving early arriving planets would probably have been very different. Even without inertial waves enhancing the dissipation, pre-main-sequence evolution could still be important simply because of the larger stellar radius. This issue therefore depends critically on the timing of hot Jupiter emplacement. If the majority of hot Jupiters arrive at their observed orbital locations early, due to disk migration or in situ formation, some of the results of this paper would need to be modified. If most hot Jupiters arrive after the star has settled onto the main sequence, possibly due to some form of high-eccentricity migration, our results are valid as they are.

Since the tidal dissipation is expected to increase as the star leaves the main sequence (due to the increased stellar radius, as well as the excitation and damping of inertial waves and internal gravity waves), the remaining lifetimes derived in Section 4.3 should be interpreted as maximum lifetimes. Our calculations predict that most hot Jupiters are not in immediate danger of being tidally destructed, with remaining lifetimes ranging from hundreds of megayears to several gigayears or more (see Figure 11). Since evidence for hot Jupiters to be tidally destroyed has been found (Hamer & Schlaufman 2019), this may indicate tidal destruction is often initiated by stellar evolution. Our calculations suggest that WASP-19b (the shortest-period planet in our sample) may be on the verge of tidal disruption, although this conclusion depends strongly on the choice of the tidal disruption criterion (see Figure 12). So far, the WASP-12 system remains the only hot Jupiter for which the period has been observed to be decreasing with time (Maciejewski et al. 2016; Yee et al. 2020). We note, however, that WASP-12 is a hot star (${T}_{\mathrm{eff}}\approx 6300\,{\rm{K}}$) and might be a subgiant (Weinberg et al. 2017; Bailey & Goodman 2019). As a result, the tidal model adopted in this paper cannot be applied to this system; the frequency dependence of the quality factor inferred by Penev et al. (2018) was based on observations of cooler main-sequence stars that spin down with age.

In addition, the magnetic braking model adopted in this paper is relatively simplistic. For example, while we have assumed the star is rotating uniformly, some degree of core-envelope decoupling may be necessary to reproduce the observed rotation periods in open clusters of different ages (e.g., Irwin et al. 2007). Allowing for the core-envelope decoupling would decrease the effective moment of inertia of the star participating in the tidal interaction, leading to more efficient obliquity realignment. Since our tidal model already predicts extremely efficient obliquity realignment, accounting for core-envelope decoupling would make our conclusion even stronger.

As a result, regardless of the possibilities for enhanced tidal dissipation due to stellar evolution, and modifications to the magnetic braking model, the main result of this paper is likely to persist even after these physical effects are taken into account. Our tidal model implies that most hot Jupiters with orbital periods less than several days should have obliquities much smaller than 1°. These values are smaller than those predicted by a tidal model with a constant stellar quality factor (see the left-hand panel of Figure 6), and they are smaller than the 6° obliquity of the Sun. Our prediction may be testable through future high-precision Rossiter–McLaughlin measurements.

This paper has been motivated by the observed obliquity-temperature for hot Jupiters. While we have shown that our tidal model efficiently damps obliquities at short orbital periods ($\lt 3$ days), at longer periods, the tidal damping is not efficient enough to erase large obliquities (see Figure 5). The typically low obliquities of planets orbiting cool stars even at longer orbital periods may instead reflect the formation conditions of these planets, or tidal effects not captured in our model.

We find it intriguing that some hot Jupiter hosts may have been tidally de-spun, rather than spun up, if the initial obliquities were retrograde. In such cases, the obliquity may be temporarily driven toward 180°, at which point the spin rate continues to decrease due to tides and magnetic braking. Eventually the star comes to rest, before spinning up in the opposite direction (see Figure 4). Thus, stars with anomalously slow rotation periods for their ages may be a clue that the obliquity was initially large. There are a few cases in our sample that appear to be in this category, subject to an ambiguity in the current spin periods. It would be interesting to measure their rotation periods more directly, i.e., based on periodic photometric variations, and test the notion that the planets have tidally de-spun their stars.

K.R.A. thanks Joel Hartman for assistance with the stellar age calculations, and Roberto Tejada Arevalo for sharing his compilation of stellar rotation periods for use in Figure 5. We also thank the referee, Adrian Barker, for constructive criticism that significantly improved the paper. This research has been supported in part by NASA grant ATP 80NSSC18K1009. K.R.A. is supported by a Lyman Spitzer, Jr. Postdoctoral Fellowship at Princeton University. The simulations presented in this paper were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.