TIDAL EVOLUTION OF THE SPIN–ORBIT ANGLE IN EXOPLANETARY SYSTEMS

More than 1000 exoplanets discovered so far have exhibited a surprising diversity in their physical properties; this diversity provides important observational clues that can be used to unveil their formation and evolution processes. In particular, unexpectedly large fractions of hot Jupiters (giant gaseous planets orbiting the central star with periods less than one week) and planets with highly eccentric and/or misaligned orbits (their orbital axes are not aligned with the spin axes of the central star) are regarded as serious challenges to conventional models of planet formation that have been proposed to explain the properties of our solar system. According to the standard core-accretion scenario, such gas giants are supposed to form beyond the ice line in nearly circular orbits. Thus it is widely believed that the discovered hot Jupiters should have formed at a large distance beyond the ice line and then somehow migrated toward the orbits close to the central star (e.g., Ida & Lin 2004a, 2004b).

A widely accepted scenario is Type II migration in which gas giants migrate inward by planet–disk interaction and halt at the inner edge of the disk (e.g., Lin et al. 1996). Other scenarios require dynamical processes after the depletion of the gas disk, including planet–planet scattering (e.g., Rasio & Ford 1996; Nagasawa et al. 2008; Nagasawa & Ida 2011), secular evolution (Wu & Lithwick 2011), and the Kozai mechanism (Kozai 1962; Wu & Murray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2012). While the former normally predicts alignment of the stellar spin axis and the planetary orbital axis, the latter dynamical processes are found to enhance the eccentricity and inclination of planets. Thus, the spin–orbit angle distribution may be a useful discriminator of the planetary migration scenarios.

The projected angle (λ) between the stellar spin and the planetary orbit of transiting planets can be determined through the Rossiter–McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924; Queloz et al. 2000; Ohta et al. 2005; Winn et al. 2005; Hirano et al. 2011). The spin–orbit angle distribution may be a useful discriminator of the planetary migration scenarios. Indeed, the RM effect has been measured successfully for around 70 transiting exoplanets as illustrated in Figure 7 below. Among them, 31 systems exhibit significant misalignment(λ > π/8), including 12 polar orbits and 7 retrograde orbits (see also Addison 2013). These unexpected and counterintuitive discoveries imply that close-in giant planets should have experienced violent dynamical processes, for instance, planet–planet scattering.

There is also a possibility that the spin–orbit misalignment may be imprinted even at the initial condition of the protoplanetary disk (e.g., Lai et al. 2011; Batygin 2013), but a reliable estimate of the distribution function of the initial spin–orbit angles is challenging. Thus in this article we focus on the dynamical process between the star and planet in the later stage of planet formation.

If the dynamical process is the major path to form hot Jupiters, one may expect that the spin–orbit angle distribution just after the formation of hot Jupiters should be very broad and even close to random. In order to be consistent with the observed distribution with some overabundance around alignment configuration, a fairly efficient physical process of the spin–orbit (re)alignment is required. While the subsequent tidal interaction between the central star and the innermost planet is generally believed to be responsible for the alignment, its detailed model is still unknown.

A conventional equilibrium tide model could realign the system but inevitably accompanies the orbital decay of the planet within a similar timescale (Barker & Ogilvie 2009; Levrard et al. 2009; Matsumura et al. 2010); see Equation (4) below. Therefore a simple equilibrium tide model does not explain the majority of the realigned systems with finite semimajor axes.

In order to solve this problem, Winn et al. (2010) proposed a decoupling model in which the stellar convective envelope is weakly coupled to its radiative core, thus preferentially reducing the realignment timescale of the stellar envelope relative to that of the orbital decay. In reality, however, it is not clear to what extent the decoupling model can explain the observed distribution since there is no systematic and quantitative study of this model. In addition, the Sun is known to have a strong coupling between the convective and radiative zones (e.g., Howe 2009), and theoretical support for this model seems weak.

Recently, Lai (2012) proposed a new model in which the damping timescale of the spin–orbit angle could be significantly smaller than that of the planetary orbit. When the stellar spin is misaligned with respect to the planetary orbital axis, one component in the tidal potential may excite inertial waves in the convective zone and provide a dynamical tidal response (e.g., Goodman & Lackner 2009). This effect increases the efficiency of the realignment process without contributing to the orbital decay.

The model was studied in detail by Rogers & Lin (2013), who numerically integrated a set of simplified equations for the semimajor axis a, the spin angular frequency Ω_s, and the spin–orbit angle Θ, originally derived by Lai (2012). They found that planetary systems with initially arbitrary spin–orbit angles have three stable configurations; parallel, antiparallel, and polar orbits.

In the present article, we use a full set of equations to trace the three-dimensional orbit of a planet (we assume that this represents the innermost planet in the case of multiplanetary systems) and the spin vector of the central star simultaneously. We also examine the tidal evolution of exoplanetary systems on a longer time-scale. As a result, we find that both antiparallel and polar orbits eventually approach the parallel orbit. We present a detailed comparison with the previous results by Rogers & Lin (2013) and argue that the full set of equations is important in understanding the longer-term evolution of the tidal model.

We briefly review the two tidal models, the equilibrium tide, and the inertial wave dissipation in Section 2 and describe how to combine the two models in our treatment. The basic set of equations is summarized explicitly in the Appendix. After a comparison with the previous result in Section 3.1, we present the evolution in the combined tidal model in Section 3.2. We also consider the dependence on the fairly uncertain parameters of the model. Section 4 is devoted to a summary and a discussion of the implications of the present result.

The entire dynamical evolution of planetary systems has many unresolved aspects including the initial structure of protoplanetary disks, formation processes of protoplanets, planetary migration, planet–planet gravitational scattering, and the tidal interaction between the central star and orbiting planets. It is beyond the scope of the present article to consider such complicated processes in a self-consistent fashion. Thus we consider a very simple system comprising a star and a hot Jupiter and focus on their tidal interaction in order to examine the dynamical behavior of the stellar spin and the planetary orbit. For definiteness, we fix the mass and radius of the star and the planet as m_s = 1 M_☉, m_p = 10⁻³ M_☉, and R_s = 1 R_☉. The initial semimajor axis of the planetary orbit is set as a_ini = 0.02 AU.

We do not assume any specific formation mechanism, but the above configuration is expected generically from any successful models for hot Jupiter formation. The eccentricity and inclination of the hot Jupiter would depend on the details of the formation mechanism. The present study, however, focuses on the tidal evolution between the star and the hot Jupiter after the orbit circularization. We then set the eccentricity of the planet to zero and consider a wide range of the initial inclination of its orbit with respect to the spin axis of the star.

We numerically solve a set of equations (Correia et al. 2011) for the stellar spin angular momentum S = I_sΩ_s (I_s is the inertia moment of the star and Ω_s is the spin angular frequency) and the orbital angular momentum L = m_pm_s/(m_p + m_s)a²Ω_p (Ω_p is the orbital mean motion) in order to explore the tidal evolution of the system. Correia et al. (2011) derived those equations for the equilibrium tide model (ET model, hereafter), which is based on the weak-friction model with a constant delay time Δt_L (Mignard 1979).

We summarize the full equations in the Appendix that directly trace the three-dimensional orbit of planets. If one neglects the eccentricity, the general relativity effect, and the tidal deformation of the planet, however, they lead to the following simplified equations for the semimajor axis a, the spin angular frequency Ω_s, and the spin–orbit angle Θ in the ET model (Lai 2012):

$$\begin{eqnarray} \frac{(\dot{a})_{\rm e}}{a} &=& -\frac{1}{\tau _{\rm e}} \left(1-\frac{\Omega _{\rm s}}{\Omega _{\rm p}}\cos \Theta \right),\qquad\quad\qquad\qquad\quad \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{(\dot{\Omega }_{\rm s})_{\rm e}}{\Omega _{\rm s}} &=& \frac{1}{\tau _{\rm e}} \left(\frac{L}{2S}\right) \left[\cos \Theta -\frac{\Omega _{\rm s}}{2\Omega _{\rm p}} (1+\cos ^{2}\Theta)\right],\qquad\quad\;\; \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} (\dot{\Theta })_{\rm e} &=& -\frac{1}{\tau _{\rm e}} \left(\frac{L}{2S}\right) \sin \Theta \left[1-\frac{\Omega _{\rm s}}{2\Omega _{\rm p}} \left(\cos \Theta -\frac{S}{L}\right)\right].\quad \end{eqnarray} \tag{ 3 }$$

The subscript e refers to the term due to the ET model.

As is clear from Equations (1)–(3), a, Ω_s, and Θ have similar damping timescales characterized by

$$\begin{eqnarray} \tau _{\rm e} &=& \left(\frac{Q_{\rm e}}{3k_{2}}\right) \left(\frac{m_{\rm s}}{m_{\rm p}}\right) \left(\frac{a}{R_{\rm s}}\right)^5 \Omega _{\rm p}^{-1} \nonumber\\ \ms &\approx & 1.3 \left(\frac{Q_{\rm e}/3k_{2}}{2\times 10^6}\right) \left(\frac{m_{\rm s}}{10^3m_{\rm p}}\right) \left(\frac{\bar{\rho }_{\rm s}}{\bar{\rho }_\odot }\right)^{5/3} \left(\frac{P_{\rm p}}{1 {\rm \,day}}\right)^{13/3} {\rm \,Gyr},\nonumber\\ \end{eqnarray} \tag{ 4 }$$

where $\bar{\rho }_{\rm s}$ is the mean density of the star, P_p is the orbital period of the planet, Q_e ≡ (2Ω_pΔt_L)⁻¹ is the tidal quality factor of the star, and k₂ is the Love number of the star that represents gravity field deformation at the stellar surface in response to an external perturbing potential of spherical harmonic degree 2. We note that k₂ depends only on the internal density distribution of the star. The specific values of Δt_L and k₂ that we employ are shown in Table 1.

The fact that the ET model is basically governed by the single timescale is inconsistent with the presence of a number of well-aligned hot Jupiters orbiting at a finite distance from the star, as long as the ET model is the dominant mechanism for the spin–orbit (re)alignment. This is why Winn et al. (2010) proposed a decoupling model in which the star's convective envelope is weakly coupled to its radiative core and Θ is damped efficiently without significant orbital decay of the planet.

Another model that we consider in detail below is proposed by Lai (2012), who pointed out the importance of the inertial waves of the star that are driven by the Coriolis force and dissipated by the tidal interaction of the misaligned S and L. More specifically, he expanded the tidal potential due to the planet orbiting at r in the spherical coordinate system (r, θ, ϕ) centered at the star with the z-axis along the stellar spin S:

$$\begin{eqnarray} &&U_2(\mathbf {r},t) = - \sum _{m, m^{\prime }}U_{m m^{\prime }}(m_{\rm p}, \Theta)r^2 Y_{2m}(\theta,\phi)e^{-im^{\prime }\Omega _{\rm p}t}, \end{eqnarray} \tag{ 5 }$$

where the index m' refers to that of spherical harmonics $Y_{2m^{\prime }}(\theta _{\rm L}, \phi _{\rm L})$ defined in the coordinate with its z-axis along the planetary orbital angular momentum L. Then he found that the only additional tidal torque in this model comes from (m, m') = (1, 0) component.

The corresponding tidal torque components due to such dynamical tides are given in Lai (2012) as

$$\begin{equation} T_{10,x} = \frac{3k_{10}}{4Q_{10}}T_{0}\sin {\Theta }\cos ^{3}{\Theta },\quad\; \end{equation} \tag{ 6 }$$

$$\begin{equation} T_{10,y} = -\frac{3k_{10}}{16}T_{0}\sin {4\Theta },\qquad\;\; \end{equation} \tag{ 7 }$$

$$\begin{equation} T_{10,z} = -\frac{3k_{10}}{4Q_{10}}T_{0}(\sin {\Theta }\cos {\Theta })^{2}, \end{equation} \tag{ 8 }$$

where k₁₀ and Q₁₀ are the dimensionless tidal Love number and tidal quality factor corresponding to the (1, 0) component, and

$$\begin{eqnarray} &&T_{0}= G \left(\frac{m_{\rm p}}{a^3}\right)^2 R_{\rm s}^5 \end{eqnarray} \tag{ 9 }$$

(cf. Murray & Dermott 1999).

We note here that Lai (2012) chooses the y-direction along the direction S × L. Thus T_{10, y} does not change the spin–orbit elements, contributing to the spin precession, but in any case we compute T_{10, y} on the basis of Equation (20) of Lai (2012). Up to the leading-order of the delay time, T_{10, y} does not depend on Q₁₀ unlike T_{10, x} nor T_{10, z}.

From the corresponding tidal torques due to the inertial wave dissipation described above, Lai (2012) derived the following equations for a, Ω_s, and Θ:

$$\begin{equation} (\dot{a})_{10}=0,\qquad\qquad\qquad\qquad\qquad \end{equation} \tag{ 10 }$$

$$\begin{equation} \frac{(\dot{\Omega }_{\rm s})_{10}}{\Omega _{\rm s}} = -\frac{1}{\tau _{10}} \left(\sin \Theta \cos \Theta \right)^{2},\qquad\qquad\; \end{equation} \tag{ 11 }$$

$$\begin{equation} (\dot{\Theta })_{10} = -\frac{1}{\tau _{10}}\sin \Theta \cos ^{2}\Theta \left(\cos \Theta +\frac{S}{L}\right). \end{equation} \tag{ 12 }$$

In the above equations, τ₁₀ is the characteristic tidal damping timescale corresponding to the (1, 0) component⁷:

$$\begin{eqnarray} \tau _{10} &=& \left(\frac{4Q_{10}}{3k_{10}}\right) \left(\frac{m_{\rm s}}{m_{\rm p}}\right) \left(\frac{a}{R_{\rm s}}\right)^{5} \left(\frac{S}{L}\right)\Omega _{\rm p}^{-1} \nonumber\\ \ms &=& 4 \left(\frac{Q_{10}}{Q_{\rm e}}\right) \left(\frac{k_{2}}{k_{10}}\right) \left(\frac{S}{L}\right) \tau _{\rm e}. \end{eqnarray} \tag{ 13 }$$

Thus the inertial wave dissipation adds the damping terms with the timescale of τ₁₀ for Θ and Ω_s, while the planetary orbit is unaffected since the (1, 0) component of the tidal potential is static in the inertial frame (Lai 2012). Thus in this model the spin–orbit angle aligns faster before the planet falls into the central star, or more strictly, closer to its Roche limit ${\approx } 2R_{\rm s}(\bar{\rho }_{\rm s}/\bar{\rho }_{\rm p})^{1/3}$, if τ₁₀ is much smaller than τ_e. The value of τ₁₀/τ_e, however, is difficult to estimate in a reliable fashion. Thus we assume a fiducial value of 10⁻³ at the start of our simulation following Rogers & Lin (2013), and examine the dependence in Section 3.2.

We numerically integrate the set of equations shown in the Appendix for the ET model. For the Lai model, we modify the equations by adding the tidal torque T₁₀ as

$$\begin{equation} \dot{\mathbf {S}} = \dot{\mathbf {S}}_{\rm (e)} + {\mathbf {T}}_{10} - {\mathbf {T}}_{10,\rm e}, \end{equation} \tag{ 14 }$$

$$\begin{equation} \dot{\mathbf {L}} = \dot{\mathbf {L}}_{\rm (e)} - {\mathbf {T}}_{10} + {\mathbf {T}}_{10,\rm e}, \end{equation} \tag{ 15 }$$

where the subscript (e) indicates the terms for the ET model, and T_{10, e} is the term that is introduced to avoid the double counting in the above equations:

$$\begin{eqnarray} \frac{{\mathbf {T}}_{10,\rm e}}{{\mathbf {T}}_{10}} &=& \left(\frac{Q_{10}}{Q_{\rm e}}\right) \left(\frac{k_{2}}{k_{10}}\right). \end{eqnarray} \tag{ 16 }$$

Table 1 summarizes the fiducial values that we employ in the simulations. We use the same values for the stellar principal moment of inertia, the Love number, and the tidal delay time that were adopted by Correia et al. (2011) for ∼1 M_☉ star.

3.1. Comparison with Previous Results

Before presenting a detailed analysis of the Lai model, we compare typical results arising from different tidal models and approximations. We consider two different tidal models: the ET model and the Lai model. Unlike in Rogers & Lin (2013), we refer to the Lai model that incorporates both the equilibrium tide and the inertial wave dissipation effects.

We numerically integrate the full set of equations in the Appendix throughout this article. If one focuses on the evolution for a, Ω_s, and Θ, the evolution equations in the Lai model can be reduced as follows:

$$\begin{equation} \dot{a} = (\dot{a})_{\rm e},\;\;\qquad\qquad\qquad \end{equation} \tag{ 17 }$$

$$\begin{equation} \dot{\Omega }_{\rm s} = (\dot{\Omega }_{\rm s})_{\rm e} + (\dot{\Omega }_{\rm s})_{10} - (\dot{\Omega }_{\rm s})_{10,\rm e}, \end{equation} \tag{ 18 }$$

$$\begin{equation} \dot{\Theta }= (\dot{\Theta })_{\rm e} + (\dot{\Theta })_{10} - (\dot{\Theta })_{10,\rm e},\; \end{equation} \tag{ 19 }$$

where

$$\begin{equation} \frac{(\dot{\Omega }_{\rm s})_{10,\rm e}}{(\dot{\Omega }_{\rm s})_{10}} =\frac{(\dot{\Theta })_{10,\rm e}}{(\dot{\Theta })_{10}} = \left(\frac{Q_{10}}{Q_{\rm e}}\right) \left(\frac{k_{2}}{k_{10}}\right). \end{equation} \tag{ 20 }$$

In order to compare with the results from the full set of equations, we also integrate the above set of equations assuming the constant value for (Q₁₀/Q_e)(k₂/k₁₀). In practice, we use Equation (13), and fix its value from the initial values of τ₁₀/τ_e and S/L:

$$\begin{equation} \frac{Q_{10}}{Q_{\rm e}}\frac{k_{2}}{k_{10}}= \frac{1}{4}\left(\frac{\tau _{10}}{\tau _{\rm e}}\right)_{\rm ini} \left(\frac{S}{L}\right)^{-1}_{\rm ini}. \end{equation} \tag{ 21 }$$

In the above simplified set of equations, m_s and m_p show up implicitly, but not explicitly, in τ₁₀, τ_e, and I_s. For definiteness, we fix $I_s=0.08m_{\rm s}R_{\rm s}^2$ (Wu & Murray 2003; Correia et al. 2011). Thus S/L and Ω_s/Ω_p are related to each other as

$$\begin{eqnarray} &&\frac{\Omega _{\rm s}}{\Omega _{\rm p}} \approx 0.23 \left(\frac{M_\odot }{m_{\rm s}}\right) \left(\frac{m_{\rm p}}{10^{-3}M_\odot }\right) \left(\frac{R_\odot }{R_{\rm s}}\right)^2 \left(\frac{a}{0.02{\rm \,AU}}\right)^2 \left(\frac{S}{L}\right).\nonumber\\ \end{eqnarray} \tag{ 22 }$$

Figure 1 shows the evolution of a (upper panels), Θ (center), and Ω_s/Ω_{s, ini} (lower) for the ET model. The left and right panels plot the results on the basis of the simplified set of equations for a, Ω_s and Θ (Equations (1)–(3)), and the full equations in the Appendix, respectively. Note that τ_e completely specifies the units of time in the simulations. Thus the time evolution is scaled with respect to the initial value of τ_e. This is why those panels are plotted against t/τ_{e, ini}.

Zoom In Zoom Out Reset image size

In each panel, we start the models with three different initial spin–orbit angles; prograde (Θ_ini = 58°), polar (Θ_ini = 90°), and retrograde (Θ_ini = 116°). They are shown in red, blue, and black lines, respectively. For each initial condition, we plot two cases with the (S/L)_ini = 0.1 (solid) and 2 (dashed) roughly corresponding to the upper and lower limits of S/L for planets observed via the RM effect (Rogers & Lin 2013). We note here that the system evolved toward Θ = 0 regardless of Θ_ini.

As pointed out earlier by several authors, the spin–orbit alignment occurs almost simultaneously with the planetary orbital decay. More strictly, Equations (1) and (3) indicate that the damping timescales of the orbit and the spin–orbit angle are given by τ_e and (S/L)τ_e, respectively. This explains the dependence on (S/L)_ini in the plots (compare solid and dashed lines).

The situation changes drastically in the Lai model, which introduces an additional timescale τ₁₀. Figure 2 shows an example of (τ₁₀/τ_e)_ini = 10⁻³. In this case, the orbital decay proceeds according to the timescale of τ_e, while the alignment timescale is controlled by τ₁₀ as well as (S/L)τ_e. Thus if τ₁₀ ≪ τ_e, one can neglect the change of a during the spin–orbit evolution. This is the approximation adopted by Rogers & Lin (2013), which corresponds to the left panels of Figure 2; we numerically integrate Equation (19) neglecting the time evolution of a and Ω_s.

Zoom In Zoom Out Reset image size

Rogers & Lin (2013) found that the system has three distinct stable configurations, i.e., antiparallel (Θ = π), polar (π/2), and parallel (0) orbits, which is easily expected from the right-hand side of Equation (19) as well. Nevertheless the damping on the order of τ_e shows up in the later stage, and the polar orbit evolves toward Θ = 0. This can be hardly recognized in Figure 2 of Rogers & Lin (2013) because we suspect they stop the integration at t = 100τ₁₀ = 0.1τ_e, before their approximation $\dot{a}=0$ becomes invalid.

In reality, the evolution beyond the epoch should be computed taking account of the change of a and Ω_s properly. The middle panels show the result and confirm that the polar orbit is a metastable configuration. Since the right-hand sides of both Equations (3) and (19) are proportional to sin Θ, the system should eventually evolve toward either Θ = 0 or π, but not π/2.

We found that the Θ = π configuration finally approaches Θ = 0 by integrating the full set of equations for the three-dimensional planetary orbit. In the antiparallel case, the spin–orbit angle has a sharp change between π and 0. This is because as the orbit damps, the stellar spin S continuously decreases according to the total angular momentum conservation. At some point, therefore, the stellar spin becomes 0, changes the direction, and then starts to increase (aligned to the orbital axis). Thus the really stable configuration is Θ = 0 alone. Nevertheless the duration of such metastable configurations is also sensitive to the choice of τ₁₀/τ_e and/or (S/L)_ini; see Figure 4. Thus the retrograde and polar–orbit systems can be observed in the real systems depending on their age.

This behavior cannot be traced properly by the simplified approach, which is based on the differential equation for Θ, combining Equations (3) and (12); the right-hand side of Equation (3) diverges at S = 0 or Ω_s = 0 and cannot be numerically integrated beyond the point. In contrast, the full set of equations in the Appendix computes S and L first, and Θ later. Thus one can compute the evolution continuously beyond S = 0. This is one of the advantages of using the full set of equations, even in the case of the simple star–planet system.

In any case, we confirm the original claim of Lai (2012) that one can have an aligned system with a finite semimajor axis as long as τ₁₀/τ_e ≪ 1 is satisfied.

3.2. Spin–Orbit Angle Evolution in the Lai Model

Now we examine the Lai model more systematically using the full set of equations. We run 30 models of a planet with a regularly spaced Θ_ini between 0° and 180°. We plot the results in Figure 3 for (S/L)_ini = 2 (left), 0.5 (middle), and 0.1 (right). All the simulations adopt (τ₁₀/τ_e)_ini = 10⁻³ so as to compare the middle panels of Figure 2 of Rogers & Lin (2013).

Zoom In Zoom Out Reset image size

As explained in Section 3.1, the system first approaches parallel, polar, or antiparallel orbits within a timescale of τ₁₀. They are plotted in black, blue, and red, respectively, so that they are easily distinguished visually. Next the polar, and subsequently antiparallel, orbits approach toward the parallel orbits in a timescale of τ_e, eventually falling below the Roche limit of the star (0.012 AU in the present case), where we stop the simulations.

The transition from Θ = 180° to 0° (red curves) happens through a state of Ω_s = 0, implying that the stellar spin starts counterrotating due to the tidal effect of the orbiting planet. As mentioned in Section 3.1, the evolution beyond Ω_s = 0 is difficult to trace with the simplified set of equations. Thus our simulations on the basis of the full set of equations are essential. Note that the transition to the three metastable configurations is fairly rapid. Therefore if the age of the system is larger than τ₁₀ and smaller than τ_e, one may expect basically three distinct spin–orbit angles, but their broad distribution as observed (cf. Figure 7 below) is not likely to be explained even taking into account the projection effect; see Figure 3 of Rogers & Lin (2013).

While Figure 3 is the main result of this article, it remains to consider the dependence on the initial ratios of (τ₁₀/τ_e)_ini and (Ω_s/Ω_p)_ini. As we will show below, the behavior presented in Figure 3 is indeed robust against the choice of those parameters.

Figure 4 shows results of initially prograde (top), polar (center), and retrograde (bottom) orbits with (S/L)_ini = 2.0, 0.5, and 0.1 for (τ₁₀/τ_e)_ini = 10⁻⁴ (black), 10⁻³ (blue), 10⁻² (green), and 10⁻¹ (red). The ratio τ₁₀/τ_e reflects the property of the stellar fluid itself, and thus is not easy to predict in a reliable fashion. Therefore we consider a fairly wide range of its possible value. The initially prograde cases approach Θ = 0 with a timescale of τ₁₀ fairly clearly. The initially polar–orbit cases stay the configuration for a significantly longer period than τ₁₀ but eventually approach Θ = 0. The initially retrograde cases are somewhat intermediate between the two. In any case, the behavior changes systematically with the value of (τ₁₀/τ_e)_ini and can be interpolated/extrapolated relatively easily from Figure 4.

Zoom In Zoom Out Reset image size

Figure 5 shows the dependence on the planet mass m_p or equivalently on (Ω_s/Ω_p)_ini through Equation (22). Since our simulations adopt m_s = 1 M_☉, R_s = 1 R_☉, and m_p = 10⁻³ M_☉, Equation (22) determines the value of Ω_s/Ω_p uniquely through (S/L)_ini. In contrast, Rogers & Lin (2013) vary the value randomly in the range of 0.1 < (Ω_s/Ω_p)_ini < 10, while they do not describe exactly how. Equation (22) implies that the corresponding values of (Ω_s/Ω_p)_ini in our model with m_p = 10⁻³ M_☉ are 0.46, 0.11, and 0.02 for (S/L)_ini = 2, 0.5, and 0.1, respectively. In order to see the dependence, we simply repeat the simulations with m_p = 10⁻² M_☉, keeping the other parameters exactly the same. Thus the m_p = 10⁻² M_☉ case corresponds to an order of magnitude increase of (Ω_s/Ω_p)_ini relative to the m_p = 10⁻³ M_☉ case. We do not run the case with m_p = 10⁻² M_☉ and (S/L)_ini = 2 since it does not satisfy the criterion Ω_p ≫ Ω_s under which the (m, m') = (1, 0) is the only excitation mode (Lai 2012). Figure 5 basically shows that the result depends on (Ω_s/Ω_p)_ini very weakly.

Zoom In Zoom Out Reset image size

We have considered tidal evolution of star–hot Jupiter systems paying particular attention to their spin–orbit alignment. We focused on the inertial wave dissipation model proposed by Lai (2012) and examined the extent to which the model reproduces the observed distribution of spin–orbit angles for transiting exoplanets.

Basically we confirmed the conclusion of Rogers & Lin (2013) that the Lai model has three distinct stable configurations—antiparallel, polar, and parallel orbits. In reality, however, the former two turn out to be metastable and approach the parallel orbits over a longer timescale as the equilibrium tide effect exceeds that of the inertial wave dissipation. We also found that the later evolution stage needs to be examined using direct three-dimensional integrations for the spin vector and orbital angular momentum vector of the system, rather than the simplified differential equations for the semimajor axis, spin angular frequency, and spin–orbit angle, even when a simple star–planet system is considered.

The relative importance of the two tidal effects is determined by the ratio of τ_e and τ₁₀. Unfortunately there is a huge uncertainty in predicting the value of each parameter. Nevertheless in order to achieve a spin–orbit alignment at a finite planetary orbit, τ₁₀ ≪ τ_e is required. In such cases, however, the alignment due to the inertial wave dissipation works too efficiently, and the observed broad distribution of Θ for transiting planets (more precisely, the projected angle λ of Θ onto the sky plane) is difficult to reproduce.

To illustrate this point, we simulate 50 systems with a planet located initially at 0.02 AU and a randomly chosen Θ_ini for (S/L)_ini = 2, 0.5, and 0.1. Figure 6 plots resulting Θ against a at four different epochs: t/τ_{e, ini} = 0, 0.03, 0.07, and 0.1. While the ET model (upper) predicts a relatively continuous correlation between Θ and a, the Lai model (lower) has the three distinct (meta-)stable states, but they subsequently become completely aligned by t = 0.1τ_e. The evolution is so rapid that even the different value of the initial semimajor axis cannot broaden the distribution significantly.

Zoom In Zoom Out Reset image size

The observed distribution of the projected spin–orbit angles λ is plotted in Figure 7 on the basis of the Holt–Rossiter–McLaughlin Encyclopedia compiled by René Heller⁸; the radial coordinate of each symbol corresponds to the logarithm of the orbital period of the planet (P∝a^3/2), and its angular coordinate represents the observed value of λ. The black and red circles indicate the innermost planets in single transiting systems and the largest planets in the multitransiting systems, respectively. The range of the solar-system planets is plotted in a blue region.

Zoom In Zoom Out Reset image size

There is a clear tendency of clumping around 0 < λ < 30° in addition to the dominance of the prograde orbits relative to the retrograde ones. Nevertheless, the distribution is rather broad and does not seem to be consistent with that expected from Figure 6. In order to make more quantitative comparison, one needs to consider the effect of projection since λ is different from Θ but corresponds to its projected angle on the sky. The effect, however, does not change the main conclusion, and we will leave the detailed analysis to our future study.

In this sense, we agree with the overall conclusion of Rogers & Lin (2013) that the Lai model, at least in the current simple star–planet system, cannot explain the wide range of observed λ. Furthermore, the current model is unlikely to explain the empirical trend that the realigned systems are preferentially found in the host stars with T_eff < 6250 K (Winn et al. 2010). Nevertheless, we have to recognize that the planetary system considered in this article is oversimplified; we ignore the outer planets that may influence the dynamics of the innermost planet significantly and the host-star dependence of the tidal parameters. In addition, we entirely neglect the dependence on the initial conditions before the tidal realignment. Thus it is premature to make a negative conclusion, and we should explore the wider range of system configurations. For this purpose, the simplified set of equations is inappropriate, and we need to integrate three-dimensional orbits of multi-planets using the full set of equations (Correia et al. 2011). Also it is important and interesting to examine three-dimensional evolution of the spin and orbital angular momenta directly instead of that of their mutual angle alone.

We thank Shoya Kamiaka for his careful reading of the manuscript and a referee for very valuable and pertinent comments, which improved the earlier manuscript. Y.S. gratefully acknowledges the supports from the Global Collaborative Research FundA Worldwide Investigation of Other Worldsgrant, the Global Scholars Program of Princeton University, and the Grant-in Aid for Scientific Research by JSPS (No. 24340035). K.M. is supported by the Leading Graduate Course for Frontiers of Mathematical Sciences and Physics. T.H. is supported by the Japan Society for Promotion of Science (JSPS) Fellowship for Research (No. 25-3183). A.T. acknowledges the support from the Grant-in-Aid for Scientific Research by JSPS (No. 24540257).

The present article considers a hybrid tidal model that combines a conventional equilibrium tide model with the inertial tidal dissipation model proposed by Lai (2012). For completeness, we explicitly give a full set of basic equations that we numerically integrate. The equations are derived by Correia et al. (2011) for the equilibrium tide model for a system of a central star and two planets with arbitrary orbital eccentricities. The tidal interaction between the stellar spin and the planetary orbit is based on the weak-friction model with a constant delay time Δt (Mignard 1979). Furthermore, the general relativity correction and the stellar oblateness are taken into account.

In the present article, we focus on the tidal interaction between the star and the innermost planet. Thus we neglect the distant planet for simplicity. The tide on the planet is not considered either because it should have a very negligible effect on the dynamics of the star and planet. Furthermore since we fix the initial eccentricity as 0 in the current simulations, the general relativity effect is also neglected.

Let us denote their mass by m_s and m_p, the semimajor axis by a, the spin angular velocity by Ω_s, the orbital angular velocity by Ω_p, radius by R_s, and gravity coefficient by $J_{2_{s}}$. The equations are written in terms of the Jacobi coordinates with r₁ being the relative position from m_s to m_p under the quadrupole approximation for the conservative and tidal potentials of the system (e.g., Smart 1953; Kaula 1964; Correia et al. 2011).

The evolution of spin of star and orbit of planet can be specified in the quadrupole approximation by two parameters;

• 1.  
  star rotational angular momentum:

  $$\begin{equation} \mathbf {S}=C_{\rm s}\Omega _{\rm s}\hat{\mathbf {s}}, \end{equation} \tag{ A1 }$$

  where $\hat{\mathbf {s}}$ is the unit vector of $\hat{\mathbf {S}}$ and C_s is the principal moment of inertia.
• 2.  
  orbital angular momentum:

  $$\begin{equation} \mathbf {L}=\beta \sqrt{{\mu } a}\hat{\mathbf {k}}, \end{equation} \tag{ A2 }$$

  where $\hat{\mathbf {k}}$ is the unit vector of $\hat{\mathbf {L}}$, μ = G(m_s + m_p), and β = m_sm_p/(m_s + m_p).We define Θ, the angle between the spin of the star and the orbit of the innermost planet via

  $$\begin{equation} \cos \Theta =\hat{\mathbf {s}}\cdot \hat{\mathbf {k}}. \end{equation} \tag{ A3 }$$

As Correia et al. (2011) described in detail, the evolution of the system is governed by the conservative motion and the tidal effects. First, the equation for conservative motion is obtained by averaging the equations of motion over the mean anomaly:

$$\begin{equation} \dot{\mathbf {L}}_{(\rm c)} = \alpha \cos \Theta \hat{\mathbf {s}}\times \hat{\mathbf {k}},\quad \end{equation} \tag{ A4 }$$

$$\begin{equation} \dot{\mathbf {S}}_{(\rm c)} = -\alpha \cos \Theta \hat{\mathbf {s}}\times \hat{\mathbf {k}}, \end{equation} \tag{ A5 }$$

where L is the norm of L and

$$\begin{equation} \alpha = \frac{3Gm_{\rm s}m_{\rm p}J_{2_{s}}R_{s}^{2}}{2a^{3}}. \end{equation} \tag{ A6 }$$

Second, the equilibrium tidal effect is considered under the quadrupole approximation of the tidal potential assuming the constant delay time Δt_L (Mignard 1979). After averaging equations of motion over the mean anomaly, one obtains

$$\begin{eqnarray} \dot{\mathbf {L}}_{(\rm t)} &=& -\dot{\mathbf {S}}_{(\rm t)}, \nonumber\\ \dot{\mathbf {S}}_{(\rm t)} &=& K \left [ \frac{\Omega _{\rm s}}{2} (\hat{\mathbf {s}}- \cos \Theta \hat{\mathbf {k}}) - \Omega _{\rm s} \hat{\mathbf {s}} + \Omega _{\rm p} \hat{\mathbf {k}}\right ], \end{eqnarray} \tag{ A7 }$$

where

$$\begin{equation} K = \Delta {t_{\rm L}} \frac{3k_{2}Gm_{\rm p}^{2}R_{\rm s}^{5}}{a^{6}}, \end{equation} \tag{ A8 }$$

Thus the total rates of change of L and S under the equilibrium tidal model are given by the sum of the above terms corresponding to the conservative motion and tidal effect:

$$\begin{equation} \dot{\mathbf {L}}_{(\rm e)} = \dot{\mathbf {L}}_{(\rm c)} +\dot{\mathbf {L}}_{(\rm t)}, \end{equation} \tag{ A9 }$$

$$\begin{equation} \dot{\mathbf {S}}_{(\rm e)} = \dot{\mathbf {S}}_{(\rm c)} +\dot{\mathbf {S}}_{(\rm t)}. \end{equation} \tag{ A10 }$$

Finally, we add the tidal torque due to the inertial wave dissipation, Equations (6)–(8) to the above equations as in Equations (14) and (15).