EARLY EXCITATION OF SPIN–ORBIT MISALIGNMENTS IN CLOSE-IN PLANETARY SYSTEMS

Typically quoted lifetimes of protoplanetary disks fall in the range ∼1–10 Myr and almost certainly depend on various parameters such as the host stellar mass (Williams & Cieza 2011). We adopt several approximations for the physical evolution of the star and disk, which are specific to Sun-like stars, which host the best observationally characterized hot Jupiters. While generally difficult to accurately parameterize, the disk mass can be taken to evolve as (Laughlin et al. 2004)

$$\begin{equation} M_{\rm {disk}} = \frac{M_{\rm {disk}}^0}{1 + t/\tau _{\rm {disk}}}. \end{equation} \tag{ 1 }$$

Interpreting the time derivative of M_disk to represent the accretionary flow, following Batygin & Adams (2013) we find that the initial disk mass, $M_{\rm {disk}}^0 = 5 \times 10^{-2} \ M_{\odot }$ and evaporation timescale τ_disk = 5 × 10⁻¹ Myr provide an acceptable match to the observations (Hartmann 2008; Herczeg & Hillenbrand 2008; Hillenbrand 2008).

For simplicity, we model the interior structure of the central star with a polytrope of index ξ = 3/2 (appropriate for a fully convective object; Chandrasekhar 1939). A polytropic body of this index is characterized by a specific moment of inertia I = 0.21 and a Love number (twice the apsidal motion constant) of k₂ = 0.14. Because T-Tauri stars derive a dominant fraction of their luminosity from gravitational contraction, we adopt the following expression for the radiative loss of binding energy (Hansen & Kawaler 1994):

$$\begin{equation} - 4 \pi R_\star ^2 \sigma T_{\rm {eff}}^4 = \left(\frac{3}{5-\xi } \right) {G M_\star ^2 \over 2 R_\star ^2} {d R_\star \over dt}. \end{equation} \tag{ 2 }$$

Equation (2) effectively dictates the process of Kelvin–Helmholtz contraction, and is satisfied by the following solution:

$$\begin{equation} R_{\star } = (R_{\star }^0) \left[1+ \left(\frac{5-\xi }{3} \right) \frac{24 \pi \sigma T_{\rm {eff}}^4}{G M_{\star } (R_{\star }^0)^3} t \right]^{-1/3}. \end{equation} \tag{ 3 }$$

A good match to the numerical evolutionary track of Siess et al. (2000) for a M_⋆ = 1 M_☉ star can be obtained by assuming an initial radius of $R_{\star }^0\, {\simeq}\, 4 \ R_{\odot }$ and an effective temperature of T_eff = 4100 K.

In order to model the magnetic disk–star interactions, we consider a T-Tauri star possessing a pure dipole magnetic field, whose north pole is aligned with the stellar spin axis. In the region of interest (i.e., in the domain of the disk), the field is current-free and can be expressed as a gradient of a scalar potential:

$$\begin{equation} \boldsymbol {B}_{\rm {dip}}=-\boldsymbol {\nabla } V. \end{equation} \tag{ 4 }$$

To retain generality, we take the field to be tilted at an angle β_i with respect to the disk plane into a direction specified by an azimuthal angle $\tilde{\phi }_i$:

$$\begin{eqnarray} V&=&B_{\star }R_\star \left (\frac{R_\star }{r}\right)^2 \left[ P_0^1(\cos (\tilde{\theta })) \cos (\beta _i)\right.\nonumber \\ &&\left. -\sin (\beta _i) (\sin (\tilde{\phi }_i)\sin (\tilde{\phi })\right. \nonumber \\ &&\left. + \cos (\tilde{\phi }_i){\rm cos}(\tilde{\phi })) P_1^1({\rm cos}(\tilde{\theta })) \right], \end{eqnarray} \tag{ 5 }$$

where B_⋆ is the stellar surface field and $P_l^m$ are associated Legendre polynomials.

If we assume the disk to be circular and Keplerian, there exists a radius $a^{\prime }_{\rm {co}}=(G\,M_{\star }/\omega ^2)^{1/3}$ at which the mean motion of the disk material, n', equals the spin rate of the stellar magnetic field, ω. At larger and smaller radii, Keplerian shear will give rise to relative fluid velocity with respect to the stellar rotation. Accordingly, as a result of the thermal ionization of alkali metals in the disk (Draine et al. 1983), the magnetic field will be dragged azimuthally by differential rotation, while slipping backward diffusively (Livio & Pringle 1992). Following Armitage & Clarke (1996), we parameterize the magnitude of the azimuthally induced field B_ϕ as a fraction, γ = B_ϕ/B_z, of the vertical component of the dipole field B_z.

As shown by Agapitou & Papaloizou (2000) and Uzdensky et al. (2002), beyond a critical value of γ ≃ 1, field lines are stretched to a sufficient degree to reconnect and transition from a closed to an open topology. Thus, the condition |γ| ≲ 1 defines a magnetically connected region within the disk with $\hat{a}^{\prime }_{\rm {in}}<a^{\prime }<\hat{a}^{\prime }_{\rm {out}}$. Outside of this region, we assume there to be no appreciable magnetic coupling to the disk.

The radial profile of γ is determined by the magnetic diffusivity of the disk, which in turn may be represented by the following dimensionless parameter (Matt & Pudritz 2004):

$$\begin{equation} \zeta \,=\,\frac{\alpha }{P_m}\frac{h}{a^{\prime }}, \end{equation} \tag{ 6 }$$

where α is the disk viscosity parameter introduced by Shakura & Sunyaev (1973), P_m is the magnetic Prandtl number, and h is the scale height of the disk. As argued by Matt & Pudritz (2004, 2005), any realistic choice of parameters yields ζ ⩽ 10⁻², which is the value we adopt throughout our work here. In terms of ζ, a steady state magnetic twist angle may be expressed as (Uzdensky et al. 2002)

$$\begin{equation} \gamma \,=\,\frac{(a^{\prime }/a^{\prime }_{\rm {co}})^{3/2}-1}{\zeta }. \end{equation} \tag{ 7 }$$

For our adopted value of ζ, this so-called magnetically connected region does not diverge from the corotation radius by more than ∼1%.

The above discussion highlights a crucial aspect of the magnetic star–disk interaction, which is discussed in detail by Matt & Pudritz (2004, 2005). If the disk is truncated at $a^{\prime }_{\rm {in}}\,>\,\hat{a}^{\prime }_{\rm {out}}$, then there is no magnetically connected region within the disk. The picture is slightly more complicated for the case where $a^{\prime }_{\rm {in}}\,<\hat{a}^{\prime }_{\rm {in}}$, as one may speculate that magnetic effects arising from differential rotation outside $a^{\prime }_{\rm {co}}$ may cancel those associated with differential rotation inside $a^{\prime }_{\rm {co}}$ to first order. In all of the following work, we circumvent these issues by assuming a disk-locked condition (Shu et al. 1994; Mohanty & Shu 2008) where $a^{\prime }_{\rm {in}}\,=\,a^{\prime }_{\rm {co}}$, but add a cautionary note that this assumption may lead to somewhat overly favorable results.

In order to derive the analytical form of the torques, we take a similar approach to that of Lai et al. (2011), and assume the disk to be razor thin. The disk current loops are envisioned to follow the magnetic field lines in a force-free fashion (see, e.g., Krasnopolsky et al. 2009; Zanni & Ferreira 2013) and connect onto the stellar surface. Accordingly, the induced azimuthal magnetic field arises from a radial current within the disk (Lai 1999).

The magnitude of the radial current is calculated using Amp$\grave{\rm {e}}$re's law (Jackson 1998) in the form

$$\begin{equation} \int _C \boldsymbol {B}\cdot d\boldsymbol {l} = \mu _0 \int \int _A \boldsymbol {J}\cdot d\boldsymbol {S}, \end{equation} \tag{ 8 }$$

where $d\boldsymbol {l}$ is a vector path length, $d\boldsymbol {S}$ is a vector area element, and C is a loop encompassing the surface A. $\boldsymbol {J}$ is the (induced) current density within the disk. Because the induced field is entirely toroidal, we can integrate the left-hand side along an azimuthal loop around the disk. This yields

$$\begin{equation} 4\pi a^{\prime } B_\phi = 2 \pi \mu _0 a^{\prime } K_{r}, \end{equation} \tag{ 9 }$$

where K_r = ∫J_rdz = 2B_ϕ/μ₀ is the inward radial surface current.

Combining this expression with Equation (7), we obtain

$$\begin{equation} K_r\,=\frac{2 B_z}{\mu _0} \left[ \frac{(a^{\prime }/a^{\prime }_{\rm {co}})^{3/2}-1}{\zeta }\right]. \end{equation} \tag{ 10 }$$

With an expression for the induced current at hand, we immediately arrive at an expression for the associated Lorentz torque, considering the induced current to interact only with the stellar dipole field:

$$\begin{equation} \boldsymbol {\tau }_L=(a^{\prime }\, {\hat{\boldsymbol { \rho }}}\,)\times (K_r\,{\hat {\boldsymbol {\rho }}} \times \boldsymbol {B}_{\rm {dip}}). \end{equation} \tag{ 11 }$$

In the above expression, $\hat{\boldsymbol \rho }$ is the radial unit vector in the plane of the disk.

At this point, in order to cast the magnetic torques into a usable form, we project $\boldsymbol {\tau }_L$ onto each of the Cartesian axes in the disk's frame and subsequently integrate over the entire magnetically connected region:

$$\begin{equation} \tau _{i^{\prime }}=\int _0^{2\pi } \int _{a^{\prime }_{\rm {co}}}^{\hat{a}^{\prime }_{\rm {out}}}{\boldsymbol \tau }_L \cdot {{\hat{\boldsymbol x}}}_{i^{\prime }}\, \rho \, d\rho \, d\phi, \end{equation} \tag{ 12 }$$

where the subscript i' represents the Cartesian axes in the disk's frame. With the variables and parameters given above, these torques evaluate to

$$\begin{eqnarray} &&\tau _{x^{\prime }}= \left(\frac{2 \pi B_{\star }^2\,R_{\star }^6\, \zeta \, \sin (\beta _i)\, \cos (\beta _i)}{3\mu _0 (1+\zeta)^2 (a^{\prime }_{\rm {co}})^3} \right) \,\cos (\tilde{\phi }_i) \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&\tau _{y^{\prime }} = \left(\frac{2 \pi B_{\star }^2\,R_\star ^6\, \zeta \, \sin (\beta _i)\, \cos (\beta _i)}{3\mu _0 (1+\zeta)^2 (a^{\prime }_{\rm {co}})^3} \right) \,\sin (\tilde{\phi }_i) \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\tau _{z^{\prime }} = \left(\frac{4\pi B_{\star }^2R_\star ^6\, \zeta \, \cos ^2(\beta _i)}{3\mu _0(1+\zeta)^2 (a^{\prime }_{\rm {co}})^3} \right). \end{eqnarray} \tag{ 15 }$$

Note that for β_i > π/2, the star and disk spin in opposite directions and so the concept of a corotation radius loses meaning. Accordingly, under the assumptions of our model, no magnetically connected region exists as presented above. Indeed, it is unclear how the magnetic field would interact with the disk in such a case and consequently, additional torques or accretional processes may have been omitted. For the purposes of our idealized model, we incorporate the loss of a magnetically connected region by artificially forcing the torques to equal zero for angles of β_i > π/2. Mathematically, this is done by multiplying the torques by an approximation to a step function $\mathcal {S}(\beta _i)$ given by

$$\begin{equation} \mathcal {S}(\beta _i)=1- \left(\frac{\pi /2+ \arctan \left(\frac{\beta _i-\pi /2}{\ell }\right)}{\pi }\right), \end{equation} \tag{ 16 }$$

where ℓ = 10⁻⁴. The importance of such a term becomes apparent once the torques are coupled to gravity and inclinations above 90 deg are naturally attained.

Angular momentum transport among neighboring annuli of the disk is facilitated by propagation of bending waves (Foucart & Lai 2011) as well as disk self-gravity (Batygin et al. 2011) and generally occurs on a much shorter timescale than magnetic tilting of the host star. Taking advantage of this, in our analyses, we assume that the effective moment of inertia of the disk around all axes is much greater than that of the star, allowing us to ignore any torques from the star on the disk and to consider $-\tau _{i^{\prime }}$ as a back-reaction on the star's dynamics.

As will become apparent below, it is beneficial to carry out all calculations in the frame of a distant, binary companion to the central star. As such, we follow the approach of Peale et al. (2014) and define Euler angles within the binary frame related to the nutation, precession and rotation of the rigid body while assuming exclusively principal axis rotation (this is an excellent approximation for a T-Tauri star, spinning at a period of 1–10 days). Specifically, $\tilde{\beta }$ is the angle between the central star's spin axis and the binary orbit normal; $\tilde{\Omega }$ is the longitude of ascending node of the star in the binary frame where $\tilde{\Omega } = 0$ implies collinear disk and stellar lines of nodes; and the third Euler angle φ is the angle through which the star rotates as it spins (φ only enters the equations as a rate of change: $\dot{\varphi } = \omega$).

The equations for the evolution of $\tilde{\beta }$ and $\tilde{\Omega }$, adapted from Peale et al. (2014) are

$$\begin{eqnarray} \frac{d\tilde{\beta }}{dt}&=&-\frac{1}{\omega } [\cos (\tilde{\beta })(-N_{\bar{x}} \sin (\tilde{\Omega })+N_{\bar{y}} \cos (\tilde{\Omega })) + N_{\bar{z}} \sin (\tilde{\beta }))], \nonumber\\ \end{eqnarray} \tag{ 17 }$$

$$\begin{equation} \frac{d\tilde{\Omega }}{dt}=-\frac{1}{\omega \, \sin (\tilde{\beta })} [N_{\bar{x}} \cos (\tilde{\Omega })+N_{\bar{y}} \sin (\tilde{\Omega })], \end{equation} \tag{ 18 }$$

where $N_{\bar{i}}$ are projected torques. Note that by fixing the disk's longitude of the ascending node at Ω' = 0, we have implicitly placed ourselves into a frame coprecessing with the disk's angular momentum vector. The effect of precession shall be included within the gravitational part of the equations and we need not retain it here.

Throughout the pre-main-sequence, we assume a constant rotation rate of the host star. Although stellar rotation is almost certainly modulated by the presence of the disk (Herbst et al. 2007; Affer et al. 2013; Bouvier 2013), the present lack of detailed understanding of the physical processes behind rotational breaking render this assumption reasonable (see Gallet & Bouvier 2013).

The projected quantities $N_{\bar{i}}$ are directly related to the torques calculated above, although the components of the torques in the disk frame, $-\tau _{i^{\prime }}$, must first be projected onto the Cartesian axes in the binary frame. Such a projection constitutes a simple rotation of co-ordinates because, as discussed below, the disk–binary inclination is a constant of motion. The rotation angle is fixed at some prescribed angle, β', counterclockwise about the x-axis. As such, the components, $N_{\bar{i}}$ are given in terms of τ_i by

$$\begin{equation} N_{\bar{x}}=-\tau _{x^{\prime }}/(I M_{\star } R_{\star }^2), \end{equation} \tag{ 19 }$$

$$\begin{equation} N_{\bar{y}}=-(\cos (\beta ^{\prime })\tau _{y^{\prime }}-\tau _{z^{\prime }} \sin (\beta ^{\prime }))/(I M_{\star } R_{\star }^2), \end{equation} \tag{ 20 }$$

$$\begin{equation} N_{\bar{z}}=-(\cos (\beta ^{\prime })\tau _{x^{\prime }}+\tau _{y^{\prime }}\sin (\beta ^{\prime }))/(I M_{\star } R_{\star }^2). \end{equation} \tag{ 21 }$$

The above equations can be used to analyze the dynamics of the central star owing to its magnetic field interacting with its protoplanetary disk. It is noteworthy that we have made no mathematical assumption of small angles, but physically we have not taken into account the changes in the parameterized geometry of the problem that arise when mutual disk–star inclinations approach β_i → π/2. Without a doubt, future calculations should consider this effect, as we expect changes in magnetic field geometry to affect the detailed nature of the relevant torques. At the same time, it is plausible that despite quantitative differences, the overall qualitative picture envisioned within the context of more precise calculations will not be in stark contrast to the parameterized model employed here.

2.3.1. Binary Star–Disk Interactions

In this section, we calculate the gravitational response of a disk to a distant, binary companion, whose orbit is taken to lie in the reference plane. We derive a Hamiltonian for the disk subject to perturbations from the companion, working under the secular approximation. In other words, we assume that the disk's outer radius $a_{\rm {out}}^\prime$ is sufficiently small compared to the binary orbit's semi-major axis $\bar{a}$ that no meaningful commensurabilities between the orbital motions exist.

The Gaussian averaging method (see chapter 7 of Murray & Dermott 2000) dictates that in the aforementioned regime, the (orbit-averaged) treatment of the gravitational interactions of the disk-companion system is mathematically equivalent to considering the companion to be a circular ring with line density $\lambda =\bar{M}/(2\pi \bar{a})$ and the disk as an infinite series of annular wires at every radius between $a^{\prime }_{\rm {in}}$ and $a^{\prime }_{\rm {out}}$ (Murray & Dermott 2000; Morbidelli et al. 2012).

It is well known that within the secular framework, the semi-major axes are constants of motion (Morbidelli 2002). Consequently, the Keplerian contribution to the Hamiltonian can be dropped, rendering the Hamiltonian of this set up, simply the total gravitational potential energy $\mathcal {U}$ possessed by the disk in the field of the companion ring:

$$\begin{equation} \mathcal {U} = {-} \int _{\rm {disk}} \int _{\rm {ring}}\frac{G}{r}\,dM_{\rm {disk}}\,dM_{\rm {ring}}, \end{equation} \tag{ 22 }$$

where r is the separation between two mass elements dM_ring (companion) and dM_disk (disk). The integral is carried out over all angles ($\bar{\phi }$) within the ring and over all radii (a') and angles (ϕ') in the disk. The evaluation of $r(\phi ^{\prime },a^{\prime },\bar{\phi })$ is a purely geometric problem and can be simplified by approximating $a_{\rm {out}}^\prime / \bar{a}\ll 1$. Under such an approximation, we expand r to second order in Equation (22) (first-order terms are axisymmetric and therefore cancel out).

In order to compute the integral (Equation (22)), we must specify the disk surface density profile, Σ. For definitiveness, in this work, we shall follow (Mestel 1963; Batygin 2012) and adopt

$$\begin{eqnarray} &&\Sigma = \Sigma _0 \left (\frac{a^{\prime }_0}{a^{\prime }} \right), \end{eqnarray} \tag{ 23 }$$

where Σ₀ is the surface density at semi-major axis $a^{\prime }_0$. We note that the passive disk model of Chiang & Goldreich (1997) is characterized by a very similar power-law index: Σ∝(a')^−15/14 (Rafikov 2012).

With this prescription, Equation (22) becomes

$$\begin{equation} \mathcal {U}=-\int _{a^{\prime }_{\rm {in}}}^{a^{\prime }_{\rm {out}}}\int _0^{2\pi }\int _0^{2\pi }\frac{G}{r}\Sigma _0\,a^{\prime }_0\,\frac{\bar{M}}{2\pi }\,d\bar{\phi }\,d\phi ^{\prime }\,da^{\prime }. \end{equation} \tag{ 24 }$$

Noting that $a_{\rm {in}}^\prime \ll a_{\rm {out}}^\prime$, we arrive at the expression for the Hamiltonian.

Switching to canonically conjugated variables, we introduce the scaled Poincar$\acute{\rm {e}}$ action-angle coordinates:

$$\begin{equation} Z^{\prime }=1-\cos (\beta ^{\prime }) \ \ \ \ \ \ \ z^{\prime } = -\Omega ^{\prime }. \end{equation} \tag{ 25 }$$

This definition of the coordinates differs from the standard definition (see chapter 2 of Murray & Dermott 2000) in that at each disk annulus, the standard definition multiplies Z' by $d\Lambda = dm^{\prime } \sqrt{G M_{\star } a^{\prime }}$, where $dm^{\prime }=2 \pi \Sigma _0 a^{\prime }_0 da^{\prime }$ is the mass of the annulus. Thus, for the variables (Equation (25)) to remain canonical, we must also scale the Hamiltonian itself in a corresponding manner:

$$\begin{eqnarray} \breve{\mathcal {U}}&=&\frac{\mathcal {U}}{2 \pi \int _{a^{\prime }_{\rm {in}}}^{a^{\prime }_{\rm {out}}} \Sigma _0 a^{\prime }_0 \sqrt{G M_{\star } a^{\prime }} da^{\prime }} \nonumber \\ &=& \frac{3 n^{\prime }_{\rm {out}}}{8} \frac{\bar{M}}{M_{\star }} \left(\frac{a^{\prime }_{\rm {out}}}{\bar{a}} \right)^3 \left[ Z^{\prime }-\frac{Z^{\prime 2}}{2} \right]. \end{eqnarray} \tag{ 26 }$$

This expression agrees with the fourth-order Lagrange–Laplace expansion of the disturbing function (Murray & Dermott 2000), where Laplace coefficients are replaced with their leading order hypergeometric series approximations (Batygin & Adams 2013).

The crucial result here is that $\breve{\mathcal {U}}$ does not depend on z', meaning that the disk inclination is exactly preserved in the binary frame:

$$\begin{eqnarray} &&\frac{dZ^{\prime }}{dt}=-\frac{\partial \breve{\mathcal {U}}}{\partial z^{\prime }}=0 \end{eqnarray} \tag{ 27 }$$

while the precession rate depends on both the companion semi-major axis and mass. As shown by Batygin (2012) via a different approach, the constancy of disk–star inclination holds even if an eccentric companion is considered. In this case, however, the precession rate is enhanced by a factor of $(1 + 3 \bar{e}/2)$.

2.3.2. Disk–Central Star Interactions

As already mentioned above, the spin rates of classical T-Tauri stars fall within the characteristic range of 1–10 days (Herbst et al. 2007; Bouvier 2013). This results in substantial rotational deformation of young stars. To an excellent approximation, the dynamical response of a spheroidal star to the gravitational potential of the disk can be modeled by considering the inertially equivalent configuration whereby a ring of mass

$$\begin{equation} \tilde{m} = \left[ \frac{3 M_{\star }^2 \omega ^2 R_{\star }^3 I^4}{4 G k_2 } \right]^{1/3} \end{equation} \tag{ 28 }$$

orbits the star with semi-major axis

$$\begin{equation} \tilde{a} = \left[ \frac{16 \omega ^2 k_2^2 R_{\star }^6}{9 I^2 G M_{\star }} \right]^{1/3}. \end{equation} \tag{ 29 }$$

Within the context of this picture, the standard perturbation techniques of celestial mechanics can be applied (Murray & Dermott 2000; Morbidelli 2002).

In the exploratory study of Batygin & Adams (2013), the gravitational torques were computed using a low mutual inclination approximation to the true potential. This simplification is limiting, especially on the quantitative level, as it forces the topology of the dynamical portrait to be that of the second fundamental model for resonance (Henrard & Lamaitre 1983), for all choices of parameters. In this work, we shall place no restrictions on the mutual inclination and adopt the series of Kaula (1962) as a representation of the potential, which utilizes the semi-major axis ratio $(\tilde{a}/a^{\prime })$ as an expansion parameter. Provided the smallness of the semi-major axis ratio inherent to our problem ($\tilde{a}/a^{\prime } \lesssim \mathcal {O}(10^{-1})$), an octupole-order expansion suffices our needs.

Written in terms of scaled canonical Poincar$\acute{\rm {e}}$ action-angle coordinates (25), the Hamiltonian that governs the dynamics of the stellar spin-axis³ under the gravitational influence of an infinitesimal wire of mass dm' reads

$$\begin{eqnarray} d\mathcal {H} &=& \frac{1}{16}\sqrt{\frac{G M_\star }{a^{\prime 3}}}\frac{dm^{\prime }}{M_\star } \left(\frac{\tilde{a}}{a^{\prime }}\right)^{3/2} \left[ (2 - 6 \tilde{Z} + 3 \tilde{Z}^2)\right. \nonumber \\ &&\nonumber \\ &&\left. \times (2 - 6 Z^{\prime } + 3 Z^{\prime 2}) + 12 \left(\sqrt{\tilde{Z}(2-\tilde{Z})} - \sqrt{\tilde{Z}^3(2-\tilde{Z})} \right)\right. \nonumber\\ &&\left. \times \left(\sqrt{\tilde{Z}(2-\tilde{Z})} - \sqrt{\tilde{Z}^3(2-\tilde{Z})} \right) \cos (\tilde{z} - z^{\prime })\right. \nonumber \\ &&\left. + 3 \tilde{Z} Z^{\prime }(\tilde{Z} - 2) (Z^{\prime } - 2) \cos (2(\tilde{z} - z^{\prime })) \right]. \end{eqnarray} \tag{ 30 }$$

As above, to obtain an expression for the Hamiltonian that incorporates the effect of the entire disk, we imagine the disk to be composed of a series of such aforementioned wires and integrate:

$$\begin{eqnarray} \mathcal {H} &=& \int d\mathcal {H}. \end{eqnarray} \tag{ 31 }$$

Recalling that the mass of each individual wire comprising the disk is $dm^{\prime }=2 \pi \Sigma _0 a^{\prime }_0 da^{\prime }$, we note that the integral (Equation (31)) runs with respect to the disk semi-major axis.

Because of the stiff dependence of Equation (30) on a', the integral (Equation (31)) is not explicitly sensitive to the location of the outer edge of the disk. Rather, it depends upon the disk's total mass, provided that the former is substantial (i.e., 10 s of AU; Anderson et al. 2013). In turn, the disk's mass is predominantly set by the disk's size, $a^{\prime }_{\rm {out}}$:

$$\begin{equation} M_{\rm {disk}} = 2 \pi \int _{a^{\prime }_{\rm {in}}}^{a^{\prime }_{\rm {out}}} a^{\prime } \Sigma _{0} \left(\frac{a^{\prime }_0}{a^{\prime }} \right) da^{\prime } \simeq 2 \pi \Sigma _{0} a^{\prime }_0 a^{\prime }_{\rm {out}}. \end{equation} \tag{ 32 }$$

In addition to the disk's physical properties, we must also prescribe its dynamical behavior to complete the specification of the problem. As already mentioned above, because the Hamiltonian (Equation (26)) is independent of the angles (i.e., is a Birkhoff normal form), the disk inclination with respect to the binary orbital plane (and therefore Z') is conserved, while the disk's nodal precession rate, ν = dz'/dt, is given by

$$\begin{equation} \nu = \frac{\partial \breve{\mathcal {U}}}{\partial Z^{\prime }} = \frac{3 n^{\prime }_{\rm {out}}}{8} \frac{\bar{M}}{M_{\star }} \left(\frac{a^{\prime }_{\rm {out}}}{\bar{a}} \right)^3 [ \,1- Z^{\prime }]. \end{equation} \tag{ 33 }$$

Consequently, $\mathcal {H}$ represents a non-autonomous one degree of freedom system.

Because the time-dependence inherent to the problem at hand is particularly simple (z' = νt), $\mathcal {H}$ can be made autonomous by employing a canonical transformation arising from the following generating function of the second kind (Goldstein 1950):

$$\begin{equation} \mathbb {G}_2 = (\tilde{z}-\nu t) \Phi, \end{equation} \tag{ 34 }$$

where $\phi = (\tilde{z}-\nu t)$ is the new angle and the new momentum is related to the old one through

$$\begin{equation} \tilde{Z} = \frac{\partial \mathbb {G}_2}{\partial \tilde{z}} = \Phi. \end{equation} \tag{ 35 }$$

Accordingly, the Hamiltonian itself is transformed as follows (Lichtenberg & Lieberman 1983):

$$\begin{equation} \mathcal {K} = \mathcal {H} - \frac{\partial \mathbb {G}_2}{\partial t}. \end{equation} \tag{ 36 }$$

Having removed explicit time dependence, we additionally scale the Hamiltonian by ν, which introduces a single dimensionless number; the "resonance proximity parameter," which encompasses the physical properties of the system. An explicit expression for $\tilde{\delta }$ reads

$$\begin{eqnarray} &&\tilde{\delta }=\frac{3}{8}\left(\frac{{n^{\prime }}_{\rm {in}}^2}{\omega \nu } \frac{M_{\rm {disk}}}{M_{\star }} \frac{a^{\prime }_{\rm {in}}}{a^{\prime }_{\rm {out}}} \right). \end{eqnarray} \tag{ 37 }$$

Following the transformations described above, the Hamiltonian takes on the following form:

$$\begin{eqnarray} \mathcal {K} &=& {-} \Phi + \frac{\tilde{\delta }}{12} [ 3 (\Phi -2) \Phi - 3 (2+3(\Phi -2)) \Phi \cos ^2(\beta ^{\prime }) \nonumber \\ &&+ 6 \sin (2 \beta ^{\prime }) (\Phi -1) \sqrt{(2-\Phi) \Phi } \cos (\phi) + 3 \sin ^2(\beta ^{\prime }) \nonumber \\ &&\times (\Phi -2) \Phi \cos (2 \phi) ]. \end{eqnarray} \tag{ 38 }$$

Although the Hamiltonian (Equation (38)) does not exhibit explicit time-dependence, it does possess implicit time-dependence through the evolution of the resonance proximity parameter, $\tilde{\delta }$. As discussed in Batygin & Adams (2013), the time-dependent variation in $\tilde{\delta }$ is primarily brought about as a result of disk mass loss and the physical evolution of the host star (via $\tilde{n}$). While both of these processes can be quite complex, for our purposes, it suffices to note that for any reasonable choice of parameters, $\tilde{\delta } \rightarrow 0$ as the system approaches main sequence.

What are the consequences of the changes in the value of $\tilde{\delta }$? Preliminary progress toward understanding this question can be made by studying the equilibria of the Hamiltonian (Equation (38)). To do so, it is particularly convenient to switch to Cartesian coordinates defined as

$$\begin{eqnarray} &&x= \sqrt{2 \Phi } \cos (\phi) \quad y= \sqrt{2 \Phi } \sin (\phi). \end{eqnarray} \tag{ 39 }$$

The angular dependence of $\mathcal {K}$ shows that all equilibria will have ϕ = 0, also implying y = 0 (Murray & Dermott 2000). The equilibrium values of x are shown as functions of $\tilde{\delta }$ in Figure 1 for three choices of disk–binary orbit inclination.

Zoom In Zoom Out Reset image size

Depending on the value of $\tilde{\delta }$, the Hamiltonian $\mathcal {K}$ possesses between two and four fixed points. Generally, for $\tilde{\delta }<1$ two stable (elliptic) fixed points exist (shown in black and green), while four fixed points (one of which is stable (shown in blue) and the other is unstable, i.e., hyperbolic (shown in red)) are guaranteed for $\tilde{\delta }>2$. There exists a critical bifurcation value $1 \leqslant \tilde{\delta }_{\rm {crit}} \leqslant 2$ where the number of fixed points is three and the stable and unstable fixed points merge into a single unstable equilibrium (Henrard & Lamaitre 1983; Peale 1986). As shown in Figure 1, $\tilde{\delta }_{\rm {crit}}$ depends on the disk–binary orbit inclination: it changes smoothly from $\tilde{\delta }_{\rm {crit}} = 1$ at $\beta ^{\prime } = 0 \,\deg$ to $\tilde{\delta }_{\rm {crit}} = 2$ at $\beta ^{\prime } = 45 \,\deg$, and back to $\tilde{\delta }_{\rm {crit}} = 1$ at $\beta ^{\prime } = 90 \,\deg$.

Physically, $\tilde{\delta }$ represents the ratio of the characteristic precession rates of the star's and disk's angular momentum vectors. It is generally safe to assume that this ratio is well above unity at the epoch when the system emerges from the embedded stage (see, e.g., Williams & Cieza 2011). Moreover, it can be argued that dissipative processes such as accretion will force the stellar spin axis to evolve toward the equilibrium point closest to alignment with the disk's angular momentum vector (see Batygin & Morbidelli 2011 for a related discussion). Provided that the changes in $\tilde{\delta }$ are slow compared to the characteristic precession period of the star, adiabatic theory dictates that the (null) phase-space area (i.e., the adiabatic invariant $\mathcal {J}$) associated with the orbit of the stellar spin axis must be approximately conserved (Henrard 1982): $\mathcal {J}_{\rm {eq}} = 0$. Consequently, as $\tilde{\delta }$ decreases in time the stellar spin axis will reside on the equilibrium solution shown in blue on Figure 1. However, once the evolutionary track of the system reaches $\tilde{\delta } = \tilde{\delta }_{\rm {crit}}$, the associated equilibrium becomes unstable.

To understand the dynamical evolution of the stellar spin axis beyond the aforementioned adiabatic trailing phase, it is useful to consider the geometry of the Hamiltonian. For the three choices of inclination, Figures 1 and 2 show a series of phase-space portraits of $\mathcal {K}$ in Cartesian coordinates. Specifically, the panels of the Figure 2 depict snapshots of the Hamiltonian flow at $\tilde{\delta } = 5$, $\tilde{\delta } = \tilde{\delta }_{\rm {crit}}$, and $\tilde{\delta } = 0$. Here, the equilibrium points of $\mathcal {K}$ are shown as gray dots, while the separatrix (i.e., homoclinic orbit) associated with the secular spin–orbit resonance is depicted as a black curve where it exists (i.e., $\tilde{\delta } \geqslant \tilde{\delta }_{\rm {crit}}$).

Zoom In Zoom Out Reset image size

Qualitatively, the following picture holds. As long as $\tilde{\delta } > \tilde{\delta }_{\rm {crit}}$, the system remains adiabatically frozen on the equilibrium point contained in the inner circulation zone of the separatrix. As $\tilde{\delta } \rightarrow \tilde{\delta }_{\rm {crit}}$, the phase space area associated with the inner circulation zone shrinks and eventually the equilibrium point on which the stellar spin-axis resides is invaded by the separatrix. Because the separatrix is characterized by an infinite period, the adiabatic principle inevitably breaks down and the conservation of $\mathcal {J}$ is momentarily violated. However, immediately after the encounter, the separatrix turns into a regular circulatory orbit and the system returns into the realm of adiabatic theory (Borderies & Goldreich 1984; Henrard 1991; Batygin & Morbidelli 2013). As such, for all subsequent evolution, the value of the adiabatic invariant remains equivalent to that of the separatrix, evaluated at the critical resonance proximity parameter (Peale 1986).

Ultimately, as the disk dissipates, $\tilde{\delta }$ approaches zero and the Hamiltonian (Equation (38)) becomes a trivial one, whose flow is represented by concentric circles on the phase plane. Accordingly, the post-encounter inclination can be calculated from the definition of the adiabatic invariant:

$$\begin{eqnarray} &&\mathcal {J} = 2\pi (1-\cos (\tilde{\beta })) = \left[ \oint \Phi _{\rm {separatrix}} d\phi \right]_{\tilde{\delta }_{\rm {crit}}}, \end{eqnarray} \tag{ 40 }$$

where the separatrix equation at critical $\tilde{\delta }$ can be obtained by substituting the value of $\mathcal {K}$ corresponding to the unstable equilibrium point into Equation (38). A noteworthy property of solution (40) is that it depends exclusively on Z'. In other words, in the adiabatic regime, the excitation of spin–orbit misalignment is independent of all system parameters except the primordial disk–binary plane inclination.

Using the approach delineated above, we have mapped out the relationship between the primordial disk–binary inclination and the final (post-encounter) stellar inclination (also with respect to the binary orbital plane). This function is shown as a purple curve on Figure 3. In the decoupled $\tilde{\delta } = 0$ regime, both the stellar and the disk's inclinations, measured with respect to the binary orbit, are preserved. However, because of differential precession, the mutual disk–star inclination undergoes cyclical variations (Batygin 2012). The maximal and minimal mutual inclinations are obtained when the stellar orbit crosses the y = 0 line with a negative and a positive values of x respectively. The associated range of mutual inclinations is depicted in Figure 3 with red lines. As shown in the Figure, a broad array of spin–orbit angles, ranging from perfectly disk-aligned states to perfectly disk-anti-aligned star states can be produced as a consequence of passage through the secular spin–orbit resonance.

Zoom In Zoom Out Reset image size

As pointed out by Lai et al. (2011; see also Lai 1999) and reviewed in Section 2.2 of this paper, a finite disk–star misalignment can be amplified by magnetic disk–star interactions. Within the context of the picture outlined above, it is tempting to assume that magnetic torques will be of no consequence prior to the encounter with the separatrix, since adiabatic invariance ensures alignment between the disk and the star. However, a more thorough examination shows that the equilibrium on which the star is envisioned to reside at $\tilde{\delta } > \tilde{\delta }_{\rm {crit}}$ deviates away from the exact disk-aligned state by a small amount.⁴ Consequently, the seed misalignment, needed for the magnetic tilting process to become active is ensured to exist from purely gravitational considerations. The amplitude of this equilibrium misalignment is shown as a function of $\tilde{\delta }$ in Figure 4 for the three choices of inclination considered above. Note that even for high values of $\tilde{\delta }$, the misalignment can be consequential (e.g., ${\sim }0.5\hbox{--}5 \,\deg$).

Zoom In Zoom Out Reset image size

The extent to which the purely gravitational picture outlined above will be altered by the incorporation of magnetic fields depends on the assumed parameters inherent to the system. This can be gathered immediately by considering the characteristic magnetic tilting timescale:

$$\begin{eqnarray} &&\tau _{B} = \left[\frac{\zeta }{3} \frac{B_{\star }^2}{\mu _0} \frac{R_{\star }^4}{I M_{\star } \Omega _{\star } {a^{\prime }}_{\rm {in}}^3} \right]^{-1}. \end{eqnarray} \tag{ 41 }$$

Evaluated using the physical parameters quoted in Section 2, a stellar rotation period of ∼5 days (Affer et al. 2013; Bouvier 2013), and a surface field of B_⋆ ≃ 1.5 kG (Johns-Krull 2007; Gregory et al. 2012), this timescale (either assuming a constant surface field or a constant magnetic dipole moment in time), along with the characteristic stellar precession timescale, and a typical disk precession timescale are plotted over a maximal disk lifetime in Figure 5.

Zoom In Zoom Out Reset image size

Owing to gravitational contraction, if an evolutionary track with a constant surface field is assumed, the cumulative effect of the magnetic torque is smaller than that if a constant dipole moment is assumed. The former situation is easily tractable within the context of the picture outlined above. First, let us imagine that system parameters are such that the secular resonance is encountered later than ∼0.5 Myr after the birth of the star (i.e., after the characteristic magnetic tilting timescale becomes considerably longer than the disk lifetime). Under this assumption, the magnetic and gravitational acquisitions of disk–star misalignment occur on separate timescales and can be treated sequentially.

It is worth recalling that the neighborhood of the (nearly) disk-aligned equilibrium of $\mathcal {K}$ (shown on Figure 2) is foliated in elliptical circulatory trajectories. This means that the system can acquire a non-zero value of $\mathcal {J}$ well before resonance crossing. Consequently, as $\tilde{\delta }$ approaches $\tilde{\delta }_{\rm {crit}}$, the trajectory will encounter the separatrix at a moment when the phase-space area occupied by the inner branch of the critical curve matches that of the orbit. In other words, the resonant encounter will take place at $\tilde{\delta } > \tilde{\delta }_{\rm {crit}}$. As a consequence, the resonant excitation of spin–orbit misalignment will occur earlier in the disk's evolution and the acquired misalignment will be somewhat different. However, blunt evaluation shows that barring unreasonable estimates, the quantitative correction is essentially negligible and is of little interest (especially given the substantial uncertainties in other, more essential parameters such as ν).

Considerably more interesting dynamical behavior can be observed if a constant magnetic moment is assumed. As shown in Figure 5, in this case the magnetic tilting timescale is comparable to the disk torquing time. Thus, rather than reasoning through the evolution within the framework of adiabatic theory, we resort to direct numerical integration of equations of motion, accounting for both, magnetic torques (Equation (17)) and gravitational torques arising from Hamiltonian (Equation (38)).

We initialize the system at the gravitational equilibrium point discussed above. The orbit is evolved for 10 Myr, adopting the same choices of disk–binary inclination as before, in addition to an almost orthogonal configuration with $\beta ^{\prime } = 85 \,\deg$. Finally, for a more candid comparison, we ignore the dependence of ν on β' and assume a disk precession period of 2π/ν = 1 Myr for all simulations. The physical evolutions of the star and the disk are assumed to proceed as described in Section 2.1.

The results of the integrations are shown in Figure 6, where the full integrations are plotted in red and solutions that only account for gravitational torques are plotted in gray for comparison. Immediately, a number of interesting features can be observed. First, in all simulations, the magnetically facilitated growth of $\mathcal {J}$ is evident, and the impulsive excitation of spin–orbit misalignment occurs substantially earlier when magnetic tilting is taken into account.

Zoom In Zoom Out Reset image size

For lower inclinations (i.e., $\beta ^{\prime } = 25 \,\deg$ and $\beta ^{\prime } = 50 \,\deg$), it is clear that magnetic tilting complicates the trajectory, although the qualitative behavior of the dynamical evolution is similar to that of the purely gravitational calculation. That is to say that the orbit evolves toward a quasi-circular state (in phase-space) as the disk dissipates. However, unlike the purely gravitational case, here the value of $\mathcal {J}$ continues to grow, even after separatrix crossing. This is largely due to magnetic torques and arises from the fact that a fraction of the orbit resides at $\beta < 90 \,\deg$. The plots of stellar inclination relative to the disk give a physical picture of the above process.

When star–disk inclinations reach values of $\beta _i \ge 90 \,\deg$, the magnetic contribution in our model becomes null and the only effect is that of secular gravitational interactions. This remains true until the oscillatory trajectory brings the star–disk inclination to $\beta _i \le 90 \,\deg$. At this point, the magnetic influence returns and repels the stellar inclination back to values of $\beta _i \ge 90 \,\deg$. This causes the inclination to oscillate in a similar fashion to the purely gravitational case, but with its trajectory is eventually excluded from $\beta _i \le 90 \,\deg$. It is worth noting that a quantitative description of this effect would be significantly altered if there is in fact some magnetic influence (not considered here) for $\beta _i \ge 90 \,\deg$.

The evolution at higher inclinations is qualitatively different, as the orbit evolves toward a fixed point (characterized by a balance between gravitational and magnetic torques) after crossing the separatrix. The $\beta ^{\prime } = 85 \,\deg$ case is particularly striking, as the phase-space plot depicts an initial condition that behaves as a repeller of the trajectory and the binary aligned equilibrium point serves as an attractor. It is interesting to note that similar behavior can be obtained by augmenting the Hamiltonian evolution with dissipative effects of accretion (see Batygin & Adams 2013 for an in-depth discussion).

These results imply that provided an advantageous prescription for the ingredients inherent to the magnetic torquing part of the calculation, the dynamical evolution of the system can be qualitatively altered. However, it is important to note that magnetic effects do not obstruct the acquisition of spin–orbit misalignment within the framework of the disk-torquing model but instead act to accelerate it. In turn, gravitational effects provide the root inclination needed for the magnetic torques to operate. Consequently, it seems reasonable to conclude that the magnetic torquing mechanism proposed by Lai et al. (2011) may play an important but nevertheless secondary role in explaining observed hot Jupiter spin–orbit misalignments.