A Fast-growing Tilt Instability of Detached Circumplanetary Disks

The size of a circumplanetary disk is determined by tidal truncation effects due to the star (Paczynski 1977; D'Angelo et al. 2002; Ayliffe & Bate 2009a). A coplanar circumplanetary disk is initially truncated at a radius of about r_out = 0.4 r_H (Martin & Lubow 2011), where the Hill sphere radius

$$\begin{eqnarray}&&{r}_{{\rm{H}}}={a}_{{\rm{p}}}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{3{M}_{{\rm{s}}}}\right)}^{1/3}.\end{eqnarray} \tag{ 1 }$$

For typical parameters the Hill sphere radius is

$$\begin{eqnarray}&&{r}_{{\rm{H}}}=0.36\left(\displaystyle \frac{{a}_{{\rm{p}}}}{5.2\,\mathrm{au}}\right){\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{10}^{-3}{M}_{{\rm{s}}}}\right)}^{1/3}\mathrm{au}.\end{eqnarray} \tag{ 2 }$$

The Hill sphere and therefore the disk size scale with the orbital semimajor axis of the planet. We note that the tidal truncation radius of a misaligned circumplanetary disk is larger than that of a coplanar disk (Lubow et al. 2015; Miranda & Lai 2015).

The viscosity of the disk is

$$\begin{eqnarray}&&\nu =\alpha {c}_{{\rm{s}}}H=\alpha {\left(\displaystyle \frac{H}{r}\right)}^{2}{r}^{2}{\rm{\Omega }},\end{eqnarray} \tag{ 3 }$$

where α is the Shakura & Sunyaev (1973) viscosity parameter, ${c}_{{\rm{s}}}$ is the disk sound speed, and H/r is the disk aspect ratio. The density of the disk falls off exponentially away from the midplane $\rho \propto \exp [-{(1/2)(z/H)}^{2}]$ (Pringle 1981), where H is the disk scale height. The disk aspect ratio of a circumplanetary disk can be significantly larger than that of a circumstellar disk with typical values in the range 0.1–0.3 (Ayliffe & Bate 2009b; Martin & Lubow 2011). The viscous timescale at the outer edge of the disk is

$$\begin{eqnarray}&&{t}_{\nu }=\displaystyle \frac{{r}_{\mathrm{out}}^{2}}{\nu ({r}_{\mathrm{out}})},\end{eqnarray} \tag{ 4 }$$

which for our typical parameters is

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\nu }}{{P}_{\mathrm{orb}}}=232{\left(\displaystyle \frac{\alpha }{0.01}\right)}^{-1}{\left(\displaystyle \frac{H/r}{0.1}\right)}^{-2}{\left(\displaystyle \frac{{r}_{\mathrm{out}}}{0.4{r}_{{\rm{H}}}}\right)}^{3/2}.\end{eqnarray} \tag{ 5 }$$

This viscous timescale is shown as a function of H/r and the disk outer radius in the right panel of Figure 1 for α = 0.01. The viscous timescale is a measure of the lifetime of the disk if there is no further accretion on to it. We do not include accretion on to the disk in our analysis of the evolution of the tilt of the disk. Accretion that occurs at close to zero inclination will act to lower the tilt of the disk.

Zoom In Zoom Out Reset image size

The angular momentum of each ring of the disk with surface density Σ(r) is defined by the complex variable W = l_x + il_y, where ${\boldsymbol{l}}$ is a unit vector parallel to the disk angular momentum. The internal torque of the disk is represented by the complex variable G. We begin with Equations (37) and (38) in Lubow & Ogilvie (2000) that describe the evolution of a disk in the binary frame that rotates with angular velocity ${{\boldsymbol{\Omega }}}_{{\bf{b}}}={{\boldsymbol{\Omega }}}_{{\bf{b}}}{{\boldsymbol{e}}}_{{\boldsymbol{z}}}$. The equations are based on a linear theory that is valid for small warps and tilts. The equations are

$$\begin{eqnarray}&&{\rm{\Sigma }}{r}^{2}{\rm{\Omega }}\left(\displaystyle \frac{\partial W}{\partial t}+i{{\rm{\Omega }}}_{{\rm{b}}}W\right)=\displaystyle \frac{1}{r}\displaystyle \frac{\partial G}{\partial r}-{Ti}(W+{W}^{* })\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial G}{\partial t}+i{{\rm{\Omega }}}_{{\rm{b}}}G-\displaystyle \frac{T}{{\rm{\Sigma }}{r}^{2}{\rm{\Omega }}}{iG}+\alpha {\rm{\Omega }}G=\displaystyle \frac{{\rm{\Sigma }}{H}^{2}{r}^{3}{{\rm{\Omega }}}^{3}}{4}\displaystyle \frac{\partial W}{\partial r},\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&T=\displaystyle \frac{{{GM}}_{{\rm{s}}}}{4{a}_{{\rm{p}}}^{2}}\left[{b}_{3/2}^{(1)}\left(\displaystyle \frac{r}{{a}_{{\rm{p}}}}\right)\right]{\rm{\Sigma }}r,\end{eqnarray} \tag{ 8 }$$

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

We seek normal modes of the form l_x, l_y, G_x, and ${G}_{y}\propto {e}^{i\omega t}$ where ω is a complex eigenvalue. We can write

$$\begin{eqnarray}&&W={W}_{+}{e}^{i\omega t}+{W}_{-}{e}^{-i{\omega }^{\star }t}\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&G={G}_{+}{e}^{i\omega t}+{G}_{-}{e}^{-i{\omega }^{\star }t},\end{eqnarray} \tag{ 10 }$$

where ω^⋆ is the complex conjugate of the eigenvalue ω. We follow Lubow & Ogilvie (2000) and expand the equations in powers of the tidal potential such that

$$\begin{eqnarray}&&T={T}^{(1)}+{T}^{(2)}+...,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{W}_{+}={{W}_{+}}^{(0)}+{{W}_{+}}^{(1)}+{{W}_{+}}^{(2)}+...,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{W}_{-}={{W}_{-}}^{(1)}+{{W}_{-}}^{(2)}+...,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{G}_{+}={{G}_{+}}^{(1)}+{{G}_{+}}^{(2)}+...,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{G}_{-}={{G}_{-}}^{(1)}+{{G}_{-}}^{(2)}+...,\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&\omega ={\omega }^{(0)}+{\omega }^{(1)}+{\omega }^{(2)}+...\end{eqnarray} \tag{ 16 }$$

To lowest order, the rigid tilt mode has

$$\begin{eqnarray}&&{\omega }^{(0)}=-{{\rm{\Omega }}}_{{\rm{b}}}\,\,\mathrm{and}\,\,{W}_{+}^{(0)}=\mathrm{const}.\end{eqnarray} \tag{ 17 }$$

To first order Equations (6) and (7) become

$$\begin{eqnarray}&&i{\omega }^{(1)}{\rm{\Sigma }}{r}^{2}{\rm{\Omega }}{W}_{+}^{(0)}=\displaystyle \frac{1}{r}\displaystyle \frac{{{dG}}_{+}^{(1)}}{{dr}}-{{iT}}^{(1)}{W}_{+}^{(0)},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&2i{{\rm{\Omega }}}_{{\rm{b}}}{\rm{\Sigma }}{r}^{2}{\rm{\Omega }}{W}_{-}^{(1)}=\displaystyle \frac{1}{r}\displaystyle \frac{{{dG}}_{-}^{(1)}}{{dr}}-{{iT}}^{(1)}{W}_{+}^{(0)\star },\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\alpha {\rm{\Omega }}{G}_{+}^{(1)}=\displaystyle \frac{{\rm{\Sigma }}{H}^{2}{r}^{3}{{\rm{\Omega }}}^{3}}{4}\displaystyle \frac{{{dW}}_{+}^{(1)}}{{dr}},\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&(2i{{\rm{\Omega }}}_{{\rm{b}}}+\alpha {\rm{\Omega }}){G}_{-}^{(1)}=\displaystyle \frac{{\rm{\Sigma }}{H}^{2}{r}^{3}{{\rm{\Omega }}}^{3}}{4}\displaystyle \frac{{{dW}}_{-}^{(1)}}{{dr}}.\end{eqnarray} \tag{ 21 }$$

The precession rate to first order can be found by integrating Equation (18) over the whole disk to find

$$\begin{eqnarray}&&{\omega }^{(1)}=-\displaystyle \frac{{\int }_{{r}_{\mathrm{in}}}^{{r}_{\mathrm{out}}}{T}^{(1)}{W}_{+}^{(0)}r\,{dr}}{{\int }_{{r}_{\mathrm{in}}}^{{r}_{\mathrm{out}}}{\rm{\Sigma }}{r}^{2}{\rm{\Omega }}{W}_{+}^{(0)}r\,{dr}}.\end{eqnarray} \tag{ 22 }$$

We use the boundary conditions that the internal disk torque vanishes at the boundaries

$$\begin{eqnarray}&&{G}_{\pm }^{(1)}({r}_{\mathrm{in}})={G}_{\pm }^{(1)}({r}_{\mathrm{out}})=0.\end{eqnarray} \tag{ 23 }$$

We solve Equations (18) to (21) for a fixed surface density profile ${\rm{\Sigma }}\propto {r}^{-3/2}$ and choose ${W}_{+}^{(0)}=1$ to find ${W}_{\pm }^{(1)}$ and ${G}_{\pm }^{(1)}$. By considering the second-order terms, the change to the inclination of the disk is determined by the tilt growth rate

$$\begin{eqnarray}&&{\mathfrak{I}}({\omega }^{(2)})=\displaystyle \frac{{\int }_{{r}_{{\rm{in}}}}^{{r}_{{\rm{out}}}}\left(\tfrac{4\alpha }{{\rm{\Sigma }}{H}^{2}{r}^{4}{{\rm{\Omega }}}^{2}}\right)\left(| {G}_{+}^{(1)}{| }^{2}-| {G}_{-}^{(1)}{| }^{2}\right)r\,{dr}}{{\int }_{{r}_{{\rm{in}}}}^{{r}_{{\rm{out}}}}{\rm{\Sigma }}{r}^{2}{\rm{\Omega }}| {W}_{+}^{(0)}{| }^{2}r\,{dr}}\end{eqnarray} \tag{ 24 }$$

(Lubow & Ogilvie 2000). The sign of this determines whether the disk tilt increases or decreases. The $| {G}_{+}^{(1)}|$ term is caused by the m = 0 component of the potential and leads to damping while the $| {G}_{-}^{(1)}|$ term is caused by the m = 2 component and causes pure growth.

We calculate the tilt growth timescale as

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\mathrm{growth}}}{{P}_{\mathrm{orb}}}=\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{b}}}}{2\pi | {\mathfrak{I}}({\omega }^{(2)})| }.\end{eqnarray} \tag{ 25 }$$

We consider a planet with a mass of ${M}_{{\rm{p}}}={10}^{-3}\,{M}_{{\rm{s}}}$. The disk extends from ${r}_{\mathrm{in}}=0.03\,{r}_{{\rm{H}}}$ up to r_out, which we vary. Note that because we parameterize the disk size in terms of the Hill sphere radius, the results presented are independent of the planet semimajor axis. The viscosity parameter is α = 0.01. The left panel of Figure 1 shows a contour plot of the tilt growth timescale as a function of the disk outer radius and the disk aspect ratio. The darkest region shows the location of the primary resonance. The resonance is where the frequency of the lowest-order global bending mode in the disk matches the tidal forcing frequency, which is twice the binary orbital frequency (2Ω_b). The trough becomes lower and narrower as α decreases (Lubow & Ogilvie 2000). Growth of the disk tilt occurs when the outer radius of the disk is close to the resonance. Even in the case of zero viscosity, disk tilt growth may occur if the disk is within the resonant band (Lubow & Ogilvie 2000). Note that larger disk aspect ratios correspond to larger viscous torques, for which we expect larger disks.

The white region shows where, for small H/r, the growth rate is positive, meaning that the disk tilt decays. The disk moves toward coplanar alignment with the orbital plane of the planet. However, for larger H/r, the growth rate is negative and the disk inclination increases. For a fixed disk outer radius, there is a narrow range of values for H/r for which the growth rate is very large.

Figure 2 shows a contour plot of the ratio of the growth timescale to the viscous timescale. The dark region shows where a disk is most unstable to the tilt instability. With a smaller α, the width of the unstable region is narrower (but at the same location) while the viscous timescale is longer. In this case, whether the disk can tilt depends sensitively upon the distribution of material within the disk. For larger α, the viscous timescale decreases and so while the disk is more unstable to tilting, there may not be sufficient time for it to occur.

Zoom In Zoom Out Reset image size

We first use the SPH code Phantom (Lodato & Price 2010; Price & Federrath 2010; Price et al. 2018) to model a misaligned circumplanetary disk. Phantom has been used extensively to model misaligned disks in binary systems (e.g., Nixon et al. 2013; Franchini et al. 2019; Smallwood et al. 2019). The simulation has two sink particles, one representing the star with mass M_s and the other representing the planet with mass M_p = 10⁻³ M_s. The accretion radius of the star is 1.4 r_H, while that of the planet is 0.03 r_H. The simulation does not have any dependence on the planet's orbital separation since all lengths are scaled to the Hill radius. The planet is in a circular orbit. The disk is initially tilted to the binary orbital plane by 10°.

The surface density of the disk is initially distributed as a power law Σ ∝ r^−3/2 between r_in = 0.03 r_H up to r_out = 0.4 r_H. We note that the initial truncation radii of the disk do not significantly affect the evolution since the density evolves quickly. We have tested both smaller and larger initial outer truncation radii. The mass of the disk does not affect the evolution since we do not include disk self-gravity, and we take the initial total disk mass to be 10⁻⁶ M_s. We consider three different initial SPH particle numbers, 250,000, 500,000, and 10⁶. The disk is locally isothermal with sound speed c_s ∝ r^−3/4. This is chosen so that α and the smoothing length $\left\langle h\right\rangle /H$ are constant over the disk (Lodato & Pringle 2007). We take the aspect ratio at r_in to be H/r = 0.2. This corresponds to H/r = 0.1 at the initial outer truncation radius. We take the Shakura & Sunyaev (1973) α parameter to be 0.01. This is the lowest value that can be physically resolved in SPH simulations at the lowest resolution we employed (Price et al. 2018). We implement the disk viscosity by adapting the SPH artificial viscosity according to the procedure described in Lodato & Price (2010) with α_AV = 0.36, 0.46, 0.58, in order of increasing resolution, and ${\beta }_{\mathrm{AV}}=2$. The circumplanetary disk is initially resolved with shell-averaged smoothing length per scale height $\left\langle h\right\rangle /H=0.28$, 0.22, and 0.17 in order of increasing resolution.

Figure 3 shows the inclination and nodal precession angle for the circumplanetary disk. The disk nodally precesses at roughly a constant rate as a result of the m = 0 tidal torque component. There are small-scale oscillations on a timescale P_orb/2 as a result of the m = 2 tidal torque component. The inclination of the disk increases in time meaning that the dissipation of the m = 2 component dominates. This is in agreement with the analytic predictions in the previous section that circumplanetary disks are unstable to tilting. We have also run simulations with varying disk aspect ratio and find that growth occurs for H/r ≳ 0.05.

Zoom In Zoom Out Reset image size

We have carried out independent grid-based simulations with Athena++ (Stone et al. 2020), using a planet-centered spherical–polar coordinate system. The simulation domain extends from 0.0286 r_H to r_H in the radial direction, from 0.2 to π-0.2 in the θ direction, and a full 2π in the ϕ direction. One level of mesh refinement is applied at θ = [0.885, 2.256]. The third-order reconstruction scheme has been adopted. The circumplanetary disk's density and temperature profiles are the same as the SPH simulations in Section 3.1. The disk surface density has an exponential tail of $\exp (-{(r/{r}_{\mathrm{cut}})}^{2})$ with r_cut = 0.3 r_H. The detailed disk and grid setup is presented in Zhu (2019). Two simulations with different resolutions have been carried out. In the low-resolution run, we have 72 × 56 × 128 grids at the base level in the r × θ × ϕ domain. With one level of mesh refinement, this is equivalent to 8 grid cells per scale height at the inner boundary and 5 grid cells per scale height at r_cut. Based on Zhu (2019), this is the minimum resolution required to study disk precession. To verify convergence, we have doubled the resolution in every direction and carried out a high-resolution simulation. Due to the computational cost, we have only run the high-resolution simulation for ∼9 P_orb.

Figure 3 shows the evolution of tilt in the grid-based simulations. The high-resolution run is almost identical to the low-resolution run, and displays similar but somewhat faster tilt growth than the corresponding SPH simulations. Figure 4 shows the disk's isodensity contours at 3 and 19 P_orb from the low-resolution simulation. Clear growth of the disk inclination is observed at 19 P_orb. We can also see that the disk is significantly depleted with time due to accretion.

Zoom In Zoom Out Reset image size