Inclination Evolution of Protoplanetary Disks Around Eccentric Binaries

TL;DR

This paper presents a theoretical analysis of how inclined protoplanetary disks around eccentric binaries can evolve to a polar alignment, challenging the traditional view of alignment with the binary's orbital plane.

Contribution

The study derives an analytic criterion for polar alignment of disks around eccentric binaries and explores the long-term evolution considering disk warp and dissipation effects.

Findings

Disks with larger aspect ratio than viscosity parameter precess as a quasi-rigid body.

A simple criterion predicts when disks evolve toward polar alignment.

Inclination evolution timescale is shorter than disk lifetime for typical parameters.

Abstract

It is usually thought that viscous torque works to align a circumbinary disk with the binary's orbital plane. However, recent numerical simulations suggest that the disk may evolve to a configuration perpendicular to the binary orbit ("polar alignment") if the binary is eccentric and the initial disk-binary inclination is sufficiently large. We carry out a theoretical study on the long-term evolution of inclined disks around eccentric binaries, calculating the disk warp profile and dissipative torque acting on the disk. For disks with aspect ratio $H/r$ larger than the viscosity parameter $\alpha$, bending wave propagation effectively makes the disk precess as a quasi-rigid body, while viscosity acts on the disk warp and twist to drive secular evolution of the disk-binary inclination. We derive a simple analytic criterion (in terms of the binary eccentricity and initial disk orientation) for the disk to evolve toward polar alignment with the eccentric binary. When the disk has a non-negligible angular momentum compared to the binary, the final "polar alignment" inclination angle is reduced from $90^\circ$. For typical protoplanetary disk parameters, the timescale of the inclination evolution is shorter than the disk lifetime, suggesting that highly-inclined disks and planets may exist orbiting eccentric binaries.

Figures

Tables

Table 1. Table 1: Binary eccentricities e b subscript 𝑒 b e_{\rm b} , with their inclinations I crit subscript 𝐼 crit I_{\rm crit} [Eq. ( 8 )] and I polar subscript 𝐼 polar I_{\rm polar} [Eq. ( 78 )], for the selected eccentric binaries with circumbinary disks. With the exception of the debris disk around 99 Herculis, all binaries have circumbinary disks aligned with the binary orbital plane within a few degrees. Binary eccentricities are from Kennedy et al. ( 2012a ) (99 Herculis), Tomkin & Popper ( 1986 ) ( α 𝛼 \alpha CrB), Pourbaix ( 2000 ) ( β 𝛽 \beta Tri), Czekala et al. ( 2016 ) (DQ Tau), Alencar et al. ( 2003 ) (AK Sco), and Boden et al. ( 2005 ) (HD 98800 B)

Binary	$e_b$	$I_{\mathrm{crit}}$	$I_{\mathrm{polar}}$
99 Herculis	0.77	${20}^{\circ }$	${28}^{\circ }$
$\alpha$ CrB	0.37	${48}^{\circ }$	${58}^{\circ }$
$\beta$ Tri	0.43	${43}^{\circ }$	${53}^{\circ }$
DQ Tau	0.57	${33}^{\circ }$	${42}^{\circ }$
AK Sco	0.47	${40}^{\circ }$	${50}^{\circ }$
HD 98800 B	0.78	${20}^{\circ }$	${27}^{\circ }$

Equations

$${\mathbf{T}}_{\mathrm{b}}=-r^{2}n{\omega }_{\mathrm{b}}\left[(1-e_{\mathrm{b}}^2)(\hat{\mathbf{l}}\cdot {\hat{\mathbf{l}}}_{\mathrm{b}}){\hat{\mathbf{l}}}_{\mathrm{b}}\times \hat{\mathbf{l}}-5(\hat{\mathbf{l}}\cdot {\mathbf{e}}_{\mathrm{b}}){\mathbf{e}}_{\mathrm{b}}\times \hat{\mathbf{l}}\right],$$

$${\omega }_{\mathrm{b}}=\frac{3G{\mu }_{\mathrm{b}}{a}_{\mathrm{b}}^2}{4r^{5}n}$$

$$\frac{\text{d}\hat{\mathbf{l}}}{\text{d}t}=-{\omega }_{\mathrm{b}}\left[(1-e_{\mathrm{b}}^2)(\hat{\mathbf{l}}\cdot {\hat{\mathbf{l}}}_{\mathrm{b}}){\hat{\mathbf{l}}}_{\mathrm{b}}\times \hat{\mathbf{l}}-5(\hat{\mathbf{l}}\cdot {\mathbf{e}}_{\mathrm{b}}){\mathbf{e}}_{\mathrm{b}}\times \hat{\mathbf{l}}\right].$$

$$\Lambda =(1-e_{\mathrm{b}}^2)(\hat{\mathbf{l}}\cdot {\hat{\mathbf{l}}}_{\mathrm{b}}{)}^2-5(\hat{\mathbf{l}}\cdot {\mathbf{e}}_{\mathrm{b}}{)}^2,$$

$${\Phi }_{\mathrm{quad}}=\frac{G{\mu }_{\mathrm{b}}{a}_{\mathrm{b}}^2}{8r^3}(1-6e_{\mathrm{b}}^2-3\Lambda ).$$

$$\hat{\mathbf{l}}$$

$$=(\cos I_e,\sin I_e\sin {\Omega }_e,\sin I_e\cos {\Omega }_e),$$

$$\Lambda =(1-e_{\mathrm{b}}^2){\cos }^{2}I-5e_{\mathrm{b}}^2{\sin }^{2}I{\sin }^2\Omega .$$

$$I_{\mathrm{crit}}={\cos }^{-1}\sqrt{\frac{5e_{\mathrm{b}}^2}{1+4e_{\mathrm{b}}^2}}={\tan }^{-1}\sqrt{\frac{1-e_{\mathrm{b}}^2}{5e_{\mathrm{b}}^2}}.$$

$$\Sigma (r)={\Sigma }_{\text{in}}\left(\frac{r_{\text{in}}}{r}\right).$$

$$L_{\mathrm{b}}={\mu }_{\mathrm{b}}\sqrt{(1-e_{\mathrm{b}}^2)G{M}_{\mathrm{b}}{a}_{\mathrm{b}}},$$

$$M_{\text{d}}$$

$$L_{\text{d}}$$

$$\frac{L_{\text{d}}}{L_{\mathrm{b}}}\simeq 0.067(1-e_{\mathrm{b}}^2{)}^{-1/2}\left(\frac{M_{\text{d}}}{0.01{\mu }_{\mathrm{b}}}\right){\left(\frac{r_{\text{out}}}{100a_{\mathrm{b}}}\right)}^{1/2}.$$

$$\frac{\text{d}{\hat{\mathbf{l}}}_{\mathrm{d}}}{\text{d}t}=⟨\frac{{\mathbf{T}}_{\mathrm{b}}}{r^{2}n}⟩,$$

$$c_{\mathrm{s}}(r)=H(r)n(r)=h\sqrt{\frac{G{M}_{\mathrm{b}}}{r_{\text{in}}}}{\left(\frac{r_{\text{in}}}{r}\right)}^{1/2},$$

$$t_{\mathrm{bw}}{\omega }_{\mathrm{b}}=0.94\left(\frac{0.1}{h}\right)\left(\frac{4{\mu }_{\mathrm{b}}}{M_{\mathrm{b}}}\right){\left(\frac{2a_{\mathrm{b}}}{r_{\text{in}}}\right)}^2{\left(\frac{r_{\text{in}}}{r}\right)}^2$$

$$\frac{\partial\bm{\hat{l}}_{\rm d}}{\partial\ln r}\bigg{|}_{\rm twist}\sim-\frac{4\alpha}{c_{\rm s}^{2}}{\bm{T}}_{\rm b},$$

$$\frac{\partial\bm{\hat{l}}_{\rm d}}{\partial\ln r}\bigg{|}_{\rm warp}\sim-\frac{4}{c_{\rm s}^{2}}\left(\frac{\kappa^{2}-n^{2}}{2n^{2}}\right)\bm{\hat{l}}_{\rm d}{\bm{\times}}{\bm{T}}_{\rm b}.$$

$$\bigg{|}\frac{\text{d}\bm{\hat{l}}_{\rm d}}{\text{d}t}\bigg{|}_{\rm visc}\sim\left\langle\frac{{\bm{T}}_{\rm b}}{r^{2}n}{\bm{\cdot}}\frac{\partial\bm{\hat{l}}}{\partial\ln r}\bigg{|}_{\rm twist}\right\rangle\sim\left\langle\frac{4\alpha}{c_{\rm s}^{2}}(r^{2}n)\omega_{\rm b}^{2}\right\rangle.$$

$${\kappa }^2-n^2=-2{\omega }_{\mathrm{b}}n{f}_{\mathrm{b}},$$

$$f_{\mathrm{b}}=\frac{1}{2}\left\{\left[3(\hat{\mathbf{l}}\cdot {\hat{\mathbf{l}}}_{\mathrm{b}}{)}^2-1\right]\left(1+\frac{3}{2}{e}_{\mathrm{b}}^2\right)-15e_{\mathrm{b}}^2(\hat{\mathbf{l}}\times {\hat{\mathbf{l}}}_{\mathrm{b}}{)}^2\right\}.$$

$$\Sigma r^{2}n\frac{\partial \hat{\mathbf{l}}}{\partial t}=\frac{1}{r}\frac{\partial \mathbf{G}}{\partial r}+\Sigma {\mathbf{T}}_{\mathrm{b}},$$

$$\frac{\partial \mathbf{G}}{\partial t}-{\omega }_{\mathrm{b}}{f}_{\mathrm{b}}\hat{\mathbf{l}}\times \mathbf{G}+\alpha \Omega \mathbf{G}=\frac{\Sigma c_{\mathrm{s}}^2{r}^{3}n}{4}\frac{\partial \hat{\mathbf{l}}}{\partial r},$$

$$\hat{\mathbf{l}}(r,t)$$

$$\mathbf{G}(r,t)$$

$${\mathbf{G}}_0(r_{\text{in}},t)={\mathbf{G}}_0(r_{\text{out}},t)=0,$$

$$\frac{\text{d}{\hat{\mathbf{l}}}_{\mathrm{d}}}{\text{d}t}=-{\bar{\omega }}_{\mathrm{b}}\left[(1-e_{\mathrm{b}}^2)({\hat{\mathbf{l}}}_{\mathrm{d}}\cdot {\hat{\mathbf{l}}}_{\mathrm{b}}){\hat{\mathbf{l}}}_{\mathrm{b}}\times {\hat{\mathbf{l}}}_{\mathrm{d}}-5({\hat{\mathbf{l}}}_{\mathrm{d}}\cdot {\mathbf{e}}_{\mathrm{b}}){\mathbf{e}}_{\mathrm{b}}\times {\hat{\mathbf{l}}}_{\mathrm{d}}\right].$$

$${\bar{\omega }}_{\mathrm{b}}$$

$$=4.97\times 10^{-5}{\left(\frac{2a_{\mathrm{b}}}{r_{\text{in}}}\right)}^2\left(\frac{4{\mu }_{\mathrm{b}}}{M_{\mathrm{b}}}\right)$$

$$\times {\left(\frac{M_{\mathrm{b}}}{2M_{\odot }}\right)}^{1/2}{\left(\frac{r_{\text{out}}}{100\text{AU}}\right)}^{-3/2}\left(\frac{2\pi }{\text{yr}}\right)$$

$${\mathbf{G}}_0(r,t)=g_{\mathrm{b}}\left[(1-e_{\mathrm{b}}^2)({\hat{\mathbf{l}}}_{\mathrm{d}}\cdot {\hat{\mathbf{l}}}_{\mathrm{b}}){\hat{\mathbf{l}}}_{\mathrm{b}}\times {\hat{\mathbf{l}}}_{\mathrm{d}}-5({\hat{\mathbf{l}}}_{\mathrm{d}}\cdot {\mathbf{e}}_{\mathrm{b}}){\mathbf{e}}_{\mathrm{b}}\times {\hat{\mathbf{l}}}_{\mathrm{d}}\right],$$

$$g_{\mathrm{b}}(r)={\int }_{r_{\text{in}}}^r\Sigma {r^{\prime }}^{3}n({\omega }_{\mathrm{b}}-{\bar{\omega }}_{\mathrm{b}})\text{d}{r}^{\prime }.$$

Full text

1 Introduction

To date, 11 transiting circumbinary planets have been detected around 9 binary star systems (Doyle et al., 2011; Kostov et al., 2013, 2014, 2016; Orosz et al., 2012a, b; Schwamb et al., 2013; Welsh et al., 2012, 2015). All planets detected have orbital planes very well aligned with their binary orbital planes, with mutual binary-planet inclinations not exceeding $3^{\circ}$. The circumbinary planet detectability is a very sensitive function of the binary-planet inclination (Martin & Triaud, 2015; Li, Holman, & Tao, 2016). If the mutual inclination is always small ($\lesssim 5^{\circ}$), then the occurance rate of circumbinary planets is comperable to that of planets around single stars, but if modest inclinations ($\gtrsim 5^{\circ}$) are common, the circumbinary planet occurance rate may be much larger (Armstrong et al., 2014). For these reasons, it is important to understand if and how a binary aligns with its circumbinary disk from which these planets form.

Observations show that most circumbinary disks tend to be aligned with their host binary orbital planes. The gas rich circumbinary disks HD 98800 B (Andrews et al., 2010), AK Sco (Czekala et al., 2015), DQ Tau (Czekala et al., 2016), and the debris circumbiniary disks $\alpha$ CrB and $\beta$ Tri (Kennedy et al., 2012b) all have mutual disk-binary inclinations not exceeding $3^{\circ}$. However, there are some notable exceptions. The circumbinary disk around KH 15D is mildly misaligned with the binary orbital plane by $\sim 10^{\circ}$-$20^{\circ}$(Winn et al., 2004; Chiang & Murray-Clay, 2004; Capelo et al., 2012). Shadows (Marino et al., 2015) and gas kinematics (Casassus et al., 2015) of the disks in HD 142527 are consistent with a misalignment of $\sim 70^{\circ}$ between the outer circumbinary disk and binary orbital plane (Lacour et al., 2016). The disks (circumbinary and two circumstellar) in the binary protostar IRS 43 are misaligned with each other and with the binary (Brinch et al., 2016). Most intriguingly, the debris disk around the eccentric $(e_{\rm b}=0.77$) binary 99 Herculis may be highly inclined: By modeling the resolved images from Hershel, Kennedy et al. (2012a) strongly favor a disk orientation where the disk angular momentum vector is inclined to the binary orbital angular momentum vector by $90^{\circ}$ (polar alignment). Kennedy et al. (2012a) also produced a model with a disk-binary inclination of $30^{\circ}$ which fits the observations, but this configuration is unlikely, since differential precession of dust due to the gravitational influence of the binary would rapidly destroy the disk.

Since star/binary formation takes place in turbulent molecular clouds (McKee & Ostriker, 2007), the gas that falls onto the central protostellar core/binary and assembles onto the disk at different times may rotate in different directions (e.g. Bate, Bonnell, & Bromm 2003, see also Bate, Lodato, & Pringle 2010; Fielding et al. 2015). In this scenario, it is reasonable to expect a newly formed binary to be surrounded by a highly misaligned circumbinary disk which forms as a result of continued gas accretion (Foucart & Lai, 2013). The observed orientations of circumbinary disks then depend on the long-term inclination evolution driven by binary-disk interactions.

Foucart & Lai (2013, 2014) studied the warping and the dissipative torque driving the inclination evolution of a circumbinary disk, assuming a circular binary. Foucart & Lai (2013) considered an infinite disk and included the effect of accretion onto the binary, while Foucart & Lai (2014) considered a more realistic disk of finite size and angular momentum, which can precess coherently around the binary. It was shown that under typical protoplanetary conditions, both viscous torque associated with disk warping and accretion torque tend to damp the mutual disk-binary inclination on timescale much shorter than the disk lifetime (a few Myr). By contrast, a circumstellar disc inside a binary can maintain large misalignment with respect to the binary orbital plane over its entire lifetime (Lubow & Ogilvie, 2000; Foucart & Lai, 2014). This is consistent with the observations that most circumbinary disks are nearly coplanar with their host binaries. On the other hand, the observed circumbinary disk misalignment (such as in KH 15D and IRS 43) can provide useful constraints on the uncertain aspects of the disc warp theory, such as non-linear effects (Ogilvie, 2006) and parametric instabilities due to disk warping (Gammie, Goodman, & Ogilvie, 2000; Ogilvie & Latter, 2013)

However, several recent numerical studies using smoothed particle hydrodynamics (SPH) suggest that other outcomes may be possible for disks around eccentric binaries. Aly et al. (2015) showed that disks around binary black holes (which typically lie in the “viscous regime” of disk warps, with the viscosity parameter $\alpha$ larger than the disk aspect ratio $H/r$; Papaloizou & Pringle 1983; Ogilvie 1999; see Sec. 3) around eccentric binaries may be driven into polar alignment. Martin & Lubow (2017) found numerically that a circumbinary protoplanetary disk (typically in the bending wave regime, with $\alpha\lesssim H/r$ Papaloizou & Lin 1995; Lubow & Ogilvie 2000), inclined to an eccentric ($e_{\rm b}=0.5$) binary by $60^{\circ}$ will evolve to a polar configuration. They suggested that this dynamical outcome arises from the combined influence of the gravitational torque on the disk from the binary and viscous torques from disk warping. They also proposed that 99 Herculis (with $e_{\rm b}=0.77$) followed such an evolution to end in the orbital configuration (polar alignment) observed today.

In this paper, we provide a theoretical anaysis to the above numerical results. In particular, we generalize the study of Foucart & Lai (2014) to apply to circumbinary disks with arbitrary disk-binary inclinations and binary eccentricities. We derive the critical condition and calculate the timescale for the disk to evolve toward polar alignment with the binary. In Section 2, we review the secular dynamics of a test particle around an eccentric binary. In Section 3, we calculate the disk warp profile and dissipative disk torques acting on the disk, and derive the requirements for the disk to evolve into polar alignment with the binary. Section 4 considers the situation when the circumbinary disk has a non-negligible angular momentum compared to the inner binary. In Section 5, we examine the back-reaction torque from the disk on the binary and the effect of gas accretion. We discuss our results in Section 6, and summarize in Section 7.

2 Test Particle Dynamics

In preparation for later sections, we review the secular dynamics of a test particle surrounding an eccentric binary (Farago & Laskar, 2010; Li, Zhou, & Zhang, 2014; Naoz et al., 2017). Consider a circular test mass with semimajor axis $r$ and orbital angular momentum unit vector $\bm{\hat{l}}$, surrounding an eccentric binary with orbital angular momentum vector $\bm{\hat{l}}_{\rm b}$, eccentricity vector $\bm{e}_{\rm b}$, semimajor axis $a_{\rm b}$, total mass $M_{\rm b}=M_{1}+M_{2}$ (where $M_{1},M_{2}$ are individual masses), and reduced mass $\mu_{\rm b}=M_{1}M_{2}/M_{\rm b}$. The orbit-averaged torque per unit mass on the test particle is (e.g. Liu, Muñoz, & Lai 2015; Petrovich 2015)

[TABLE]

where $n\simeq\sqrt{GM_{\rm b}/r^{3}}$ is the test particle orbital frequency (mean-motion), and

[TABLE]

characterizes the precession frequency of the test particle around the binary. The torque ${\bm{T}}_{\rm b}$ in Eq. (1) is evaluated to the lowest order in $a_{\rm b}/r$.

The time evolution of the test particle’s orbital angular momentum vector is given by

[TABLE]

Equation (3) can be solved analytically (Landau & Lifshitz, 1969; Farago & Laskar, 2010; Li, Zhou, & Zhang, 2014), but the dynamics may be easily understood by analyzing the energy curves. Equation (3) has an integral of motion

[TABLE]

which is simply related to the quadrupole interaction energy (double-averaged over the two orbits) by (e.g. Tremaine, Touma, & Namouni 2009; Tremaine & Yavetz 2014; Liu, Muñoz, & Lai 2015)

[TABLE]

To plot the energy curves, we set up the Cartesian coordinate system $(x,y,z)$, where $\bm{\hat{l}}_{\rm b}=\bm{\hat{z}}$ and $\bm{e}_{\rm b}=e_{\rm b}\bm{\hat{x}}$. We may write

[TABLE]

where $I$ is the angle between $\bm{\hat{l}}$ and $\bm{\hat{l}}_{\rm b}$, $\Omega$ is the test particle’s longitude of the ascending node (measured in the $xy$ plane from the $x$-axis); similarly $I_{e}$ is the angle between $\bm{\hat{l}}$ and $\bm{e}_{\rm b}$, and $\Omega_{e}$ measures the longitude of the node in the $yz$ plane (see Fig. 1). In terms of $I$ and $\Omega$, we have

[TABLE]

In Figure 2, we plot the test particle trajectories in the $I-\Omega$ (left panel) and $I_{e}-\Omega_{e}$ (right panel) planes for the binary eccentricity $e_{\rm b}=0.3$. The critical separatrix $\Lambda=0$ is displayed in black in both plots. When $\Lambda>0$, $\bm{\hat{l}}$ precesses around $\bm{\hat{l}}_{\rm b}$ with $I\sim\text{constant}$ and $\Omega$ circulating the full range ($0-360^{\circ}$), while $\Omega_{e}$ librates around $0^{\circ}$. When $\Lambda<0$, $\bm{\hat{l}}$ precesses around $\bm{e}_{\rm b}$ with $I_{e}\sim\text{constant}$ and $\Omega_{e}$ circulating the full range ($0-360^{\circ}$), while $\Omega$ librates around $90^{\circ}$ (see Fig. 1).

Thus, the test particle angular momentum axis $\bm{\hat{l}}$ transitions from precession around $\bm{\hat{l}}_{\rm b}$ for $\Lambda>0$ to precession around $\bm{e}_{\rm b}$ for $\Lambda<0$. Because the $\Lambda=0$ separatrix has $\Omega\in[0^{\circ},360^{\circ}]$ (Fig. 2), a necessary condition for $\bm{\hat{l}}$ to precess around $\bm{e}_{\rm b}$ is $I_{\rm crit}<I<180^{\circ}-I_{\rm crit}$, where

[TABLE]

Figure 2 clearly reveals the stable fixed points of the system. In terms of the variables ($\Omega,I$), the stable fixed points (where $\text{d}I/\text{d}t=\text{d}\Omega/\text{d}t=0$) are $I=\pi/2$ and $\Omega=\pi/2,3\pi/2$, corresponding to $\bm{\hat{l}}=\pm\bm{e}_{\rm b}/e_{\rm b}$. In terms of the variables ($\Omega_{e},I_{e}$), the fixed points are $I_{e}=\pi/2$ and $\Omega_{e}=0,\pi$, corresponding to $\bm{\hat{l}}=\pm\bm{\hat{l}}_{\rm b}$. We will see in Section 3 that in the presence of dissipation, the disk may be driven toward one of these fixed points.

3 Circumbinary Disk Dynamics

We now consider a binary (with the same parameters as in Section 2, see also Fig. 3) surrounded by a circular circumbinary disk with inner truncation radius $r_{\text{in}}$, outer truncation radius $r_{\text{out}}$, with unit angular momentum vector $\bm{\hat{l}}=\bm{\hat{l}}(r,t)$, and surface density $\Sigma=\Sigma(r)$. For concreteness, we adopt the surface density profile

[TABLE]

We assume $r_{\text{in}}\ll r_{\text{out}}$ throughout this work. We could assume a more general surface density profile $\Sigma\propto r^{-p}$, with $p$ observationally constrained to lie in the range $0.5-1.5$ (e.g. Weidenschilling 1977; Williams & Cieza 2011; Chiang & Laughlin 2013). A more general $p$ will effect the disk mass [Eq. (11)] and angular momentum [Eq. (12)], as well as the precession [Eq. (28)] and viscous [Eq. (47)] rates, by factors of order unity.

The binary has orbital angular momentum

[TABLE]

while the disk has mass

[TABLE]

and angular momentum (assuming a small disk warp; see below)

[TABLE]

where $n(r)\simeq\sqrt{GM_{\rm b}/r^{3}}$. Comparing $L_{\rm b}$ to $L_{\text{d}}$, we have

[TABLE]

Because $L_{\rm b}\gg L_{\text{d}}$ for typical circumbinary disk parameters, in this section we assume $\bm{\hat{l}}_{\rm b}$ and $\bm{e}_{\rm b}$ are fixed in time, neglecting the back-reaction torque on the binary from the disk. We discuss the system’s dynamics when $L_{\text{d}}$ is non-negligible compared to $L_{\rm b}$ in Section 4, and the effects of the back-reaction torque on the binary from the disk in Section 5.

3.1 Qualitative Discussion

Assuming the disk to be nearly flat, the time evolution of the disk unit angular momentum vector is given by

[TABLE]

where ${\bm{T}}_{\rm b}$ is given in Eq. (1), $\bm{\hat{l}}_{\rm d}(t)$ is a suitably averaged unit angular momentum of the disk [see Eq. (24)], and $\langle\dots\rangle$ implies a proper average over $r$ [see Eq. (27)]. When the disk is flat, the time evolution of $\bm{\hat{l}}_{\rm d}$ is identical to that of a test particle (see discussion at the end of Sec. 3.2).

When $\alpha\lesssim H/r$ ($H$ is the disk scaleheight, $\alpha$ is the viscosity parameter), the main internal torque enforcing disk rigidity and coherent precession comes from bending wave propigation (Papaloizou & Lin, 1995; Lubow & Ogilvie, 2000). As bending waves travel at $1/2$ the sound speed, the wave crossing time is of order $t_{\rm bw}=2r/c_{\rm s}$. When $t_{\rm bw}$ is longer than the characteristic precession time $\omega_{\rm b}^{-1}$ [see Eq. (2)], strong disk warps can be induced. In the extreme nonlinear regime, disk breaking may be possible in circumbinary disks (Larwood & Papaloizou, 1997; Facchini, Lodato, & Price, 2013; Nixon, King, & Price, 2013). To compare $t_{\rm bw}$ with $\omega_{\rm b}^{-1}$, we adopt the disk sound speed profile

[TABLE]

where $h=H/r$. We find

[TABLE]

Thus, we expect that the small warp approximation should be valid everywhere in the disk except the inner-most region. Throughout this paper, we scale our results to $h=0.1$. Real protoplanetary disks can have aspect ratios in the range $h\sim 0.03-0.2$ (e.g. Lynden-Bell & Pringle 1974; Chiang & Goldreich 1997; Williams & Cieza 2011). We normalize $r_{\text{in}}$ to $2a_{\rm b}$, but note that the inner truncation radius of the disk depends non-trivially on the binary’s eccentricity (Miranda, Muñoz, & Lai, 2017).

Although the disk is flat to a good approximation, the interplay between disk twist/warp and viscous dissipation may modify the disk’s dynamics over timescales much longer than $\omega_{\rm b}^{-1}$. When the external torque ${\bm{T}}_{\rm b}$ is applied to the disk in the bending wave regime, the disk’s viscosity causes the disk to develop a small twist of order

[TABLE]

while the precession of bending waves from a non-Keplarian epicyclic frequency $\kappa$ causes the disk to develop a small warp, of order

[TABLE]

The viscous twist [Eq. (17)] interacts with the external torque, effecting the evolution of $\bm{\hat{l}}$ over the viscous timescale. To an order of magnitude, we have

[TABLE]

In the above estimate, we have assumed the relevant misalignment angles (between $\bm{\hat{l}}_{\rm d}$ and $\bm{\hat{l}}_{\rm b}$, or between $\bm{\hat{l}}_{\rm d}$ and $\bm{e}_{\rm b}$) is of order unity.

3.2 Formalism

The torque per unit mass on the disk from the inner binary is given by Eq. (1), with ${\bm{T}}_{\rm b}={\bm{T}}_{\rm b}(r,t)$. In addition, the gravitational potential from the binary induces a non-Keplarian angular frequency (Miranda & Lai, 2015), with

[TABLE]

where

[TABLE]

When the Shakura-Sunaev $\alpha$-viscosity parameter satisfies $\alpha\lesssim H/r$, the disk lies in the bending wave regime (Papaloizou & Lin, 1995; Lubow & Ogilvie, 2000). Any warp induced by an external torque is smoothed by bending waves passing through the disk. Protoplanetary disks typically lie in the bending wave regime. The time evolution of $\bm{\hat{l}}(r,t)$ is governed by the equations (Lubow & Ogilvie 2000; see also Lubow, Ogilvie, & Pringle 2002)

[TABLE]

where ${\bm{G}}$ is the internal torque.

From equation (16), we see that $t_{\rm bw}<\omega_{\rm b}^{-1}$ for standard circumbinary disk parameters, so the disk should be only mildly warped. We may therefore expand

[TABLE]

where $\bm{\hat{l}}_{\rm d}$ is the unit vector along the total angular momentum of the disk, $|{\bm{l}}_{1}|\ll|\bm{\hat{l}}_{\rm d}|=1$ [see Eqs. (31)-(32) below]. As we will see, the internal torque ${\bm{G}}_{0}(r,t)$ maintains the rigid body dynamical evolution of $\bm{\hat{l}}_{\rm d}$, while ${\bm{G}}_{1}(r,t)$ maintains the warp profile ${\bm{l}}_{1}$. Perturbative expansions to study warped disk structure and time evolution have been taken by Lubow & Ogilvie (2000, 2001) and Foucart & Lai (2014). Inserting (24) into Eq. (22), integrating over $r\text{d}r$, and using the zero torque boundary condition

[TABLE]

we find the leading order time evolution of $\bm{\hat{l}}$ is given by

[TABLE]

Here,

[TABLE]

is the characteristic precession frequency of the rigid disk. Equation (27) is equivalent to Equation (3) if one replaces $\bar{\omega}_{\rm b}$ with $\omega_{\rm b}$, and the disk dynamics reduce to those of a test particle with $\bm{\hat{l}}=\bm{\hat{l}}_{\rm d}$ when $c_{\rm s}\to\infty$.

3.3 Disk Warp Profile

With $\bm{\hat{l}}_{\rm d}$ determined with boundary condition (26), we may solve for ${\bm{G}}_{0}(r,t)$:

[TABLE]

where

[TABLE]

With the leading order terms for $\bm{\hat{l}}$ and ${\bm{G}}$, we may solve for ${\bm{l}}_{1}$. We impose the normalization condition

[TABLE]

so that $\bm{\hat{l}}_{\rm d}$ is the unit vector along the total angular momentum of the disk, or

[TABLE]

Inserting Eq. (29) into Eq. (23) and integrating, we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

Here,

[TABLE]

and

[TABLE]

The third term in Eq. (33) arises from the fact that $\bm{\hat{l}}_{\rm d}{\bm{\cdot}}\bm{e}_{\rm b}$ and $\bm{\hat{l}}_{\rm d}{\bm{\cdot}}\bm{\hat{l}}_{\rm b}$ are not constant in time, and is dynamically unimportant. Although it is strait-forward to compute the integrals in Eqs. (36)-(38), this calculation is tedious and unilluminating. Instead, we notice that over most of the region in the integrals, the internal torque radial function $g_{\rm b}(r)$ is of order

[TABLE]

Evaluating the warp functions and using the fact that $r_{\text{in}}\ll r_{\text{out}}$, we obtain the approximate expressions

[TABLE]

In Figure 4, we plot the rescaled warp functions $\tilde{\tau}_{\rm b}=\tau_{\rm b}/[\tau_{\rm b}(r_{\text{in}})-\tau_{\rm b}(r_{\text{out}})]$, $\tilde{V}_{\rm b}=V_{\rm b}/[V_{\rm b}(r_{\text{in}})-V_{\rm b}(r_{\text{out}})]$, and $\tilde{W}_{\rm bb}=W_{\rm bb}/[W_{\rm bb}(r_{\text{in}})-W_{\rm bb}(r_{\text{out}})]$.

3.4 Viscous Torques

The disk twisting due to viscosity $({\bm{l}}_{1})_{\rm twist}$ [Eq. (34)] interacts with ${\bm{T}}_{\rm b}$ [Eq. (1)], effecting the evolution of $\bm{\hat{l}}_{\rm d}$ over viscous timescales. Inserting $\bm{\hat{l}}=\bm{\hat{l}}_{\rm d}+({\bm{l}}_{1})_{\rm twist}$ into Equation (22), integrating over $2\pi r\text{d}r$, and using the boundary condition (26), we obtain

[TABLE]

Here,

[TABLE]

Using the approximate expression of $g_{\rm b}$ [Eq. (42)], one may easily reproduce Eq. (19). The same argument used in the calculation of Eqs. (43)-(45) may be used to calculate the approximate expression of the viscous rate $\gamma_{\rm b}$:

[TABLE]

We choose to normalize $\gamma_{\rm b}$ by $\alpha=0.01$; real protoplanetary disks may have $\alpha$ in the range $10^{-1}-10^{-5}$ (Rafikov, 2017). Since

[TABLE]

the viscous dissipation from disk twisting effects the evolution of $\bm{\hat{l}}_{\rm d}$ according to

[TABLE]

where $\Lambda$ is given by Eq. (4), except we replace $\bm{\hat{l}}$ by $\bm{\hat{l}}_{\rm d}$:

[TABLE]

Equation (50) is the main result of our technical calculation. We see

[TABLE]

Because $-5e_{\rm b}^{2}<\Lambda<(1-e_{\rm b}^{2})$ [Eq. (51)], Equations (52)-(53) show the system has two different end-states depending on the initial value for $\Lambda$:

$\Lambda>0$: The viscous torque (50) pushes $\bm{\hat{l}}_{\rm d}$ towards $\bm{\hat{l}}_{\rm b}$. The final state of $\bm{\hat{l}}_{\rm d}$ is alignment (if $\bm{\hat{l}}_{\rm b}{\bm{\cdot}}\bm{\hat{l}}_{\rm d}>0$ initially) or anti-alignment (if $\bm{\hat{l}}_{\rm b}{\bm{\cdot}}\bm{\hat{l}}_{\rm d}<0$ initially) with $\bm{\hat{l}}_{\rm b}$. 2. 2.

$\Lambda<0$: The viscous torque (50) pushes $\bm{\hat{l}}_{\rm d}$ towards $\bm{e}_{\rm b}$. The final state of $\bm{\hat{l}}_{\rm d}$ is alignment (or anti-alignment) with $\bm{e}_{\rm b}$.

Figure 5 shows several examples of the results for the evolution of disk orientation, obtained by integrating the time evolution of $\bm{\hat{l}}_{\rm d}$, including gravitational [Eq. (27)] and viscous [Eq. (50)] torques. On the left panels, we plot the disk inclination $I$ with time, for the binary eccentricities indicated. We choose the initial $\Omega(0)=90^{\circ}$ for all cases, so that $I<I_{\rm crit}$ ($I>I_{\rm crit}$) corresponds exactly to $\Lambda>0$ ($\Lambda<0$) (see Eqs. (4) and (8)). Thus we expect $I\to 0^{\circ}$ when $I<I_{\rm crit}$, $I\to 90^{\circ}$ when $I_{\rm crit}<I<180^{\circ}-I_{\rm crit}$, and $I\to 180^{\circ}$ when $I>180^{\circ}-I_{\rm crit}$. On the right panels of Figure 5, we plot the disk trajectories on the $I-\Omega$ plane [Eq. (6) with $\bm{\hat{l}}\to\bm{\hat{l}}_{\rm d}$]. Again, we see when $I<I_{\rm crit}$ ($\Lambda>0$), $\bm{\hat{l}}_{\rm d}$ aligns with $\bm{\hat{l}}_{\rm b}$, while when $I>I_{\rm crit}$ ($\Lambda<0$), $\bm{\hat{l}}_{\rm d}$ aligns with $\bm{e}_{\rm b}$, as expected.

4 Secular Dynamics with Massive Inclined Outer Body

Sections 2-3 neglected the circumbinary disk’s angular momentum, a valid assumption as long as $L_{\text{d}}\ll L_{\rm b}$ [Eq. (13)]. When $L_{\text{d}}\gtrsim L_{\rm b}$, the non-zero disk angular momentum will change the locations of the fixed points of the system, and hence may effect its dynamical evolution over viscous timescales.

Consider the setup of Section 2, except we now include the outer body’s mass $m$ and angular momentum ${\bm{L}}=m\sqrt{GM_{\rm b}r}\bm{\hat{l}}$. The evolution equations for $\bm{\hat{l}}$, $\bm{j}_{\rm b}=\sqrt{1-e_{\rm b}^{2}}\bm{\hat{l}}_{\rm b}$, and $\bm{e}_{\rm b}$ are [Liu, Muñoz, & Lai 2015; Eqs. (17)-(19)]

[TABLE]

where

[TABLE]

and $\mu=mM_{\rm b}/(m+M_{\rm b})$. Equation (58) reduces to Eq. (2) when $m\to 0$. The conservations of total quadrupole potential energy [see Eq. (5)] and total angular momentum yield two constants of motion (e.g. Liu, Muñoz, & Lai 2015; Anderson, Lai, & Storch 2017)

[TABLE]

For a given $K$, one may solve Eq. (60) to get $e_{\rm b}^{2}=e_{\rm b}^{2}(I)$. Assuming $0\leq I\leq\pi/2$ and requiring $0\leq e_{\rm b}<1$, we obtain

[TABLE]

Equation (59) then gives $\Psi=\Psi(I,\Omega)$. When $J\sim K^{-1}\ll 1$, Eq. (61) reduces to

[TABLE]

while when $J\gg 1$, Eq. (61) becomes

[TABLE]

The fixed points of the system in the $I-\Omega$ plane are determined by

[TABLE]

The condition $\partial\Psi/\partial\Omega=0$ gives $\Omega=\pi/2$ and $\Omega=3\pi/2$, as before (see Sec. 2). For arbitrary $J$, one must numerically solve $\partial\Psi/\partial I|_{\Omega=\pi/2,3\pi/2}=0$ to calculate the fixed points $I=I_{\rm fp}>0$ ($I=0$ is always a fixed point of the system). However, when $J\ll 1$, one may show analytically that (as found in Sec. 2)

[TABLE]

while when $J\gg 1$,

[TABLE]

where $e_{\rm b}^{2}=e_{\rm b}^{2}(0)$. Notice $I_{\rm fp}$ is the Lidov-Kozai critical inclination when $J\gg 1$ and $e_{\rm b}(0)=0$ (Lidov, 1962; Kozai, 1962).

Figure 6 plots trajectories of the system in the $I-\Omega$ and $e_{\rm b}-\Omega$ planes. When $J\ll 1$, the system’s dynamics reduce to that discussed in Section 2, with $I_{\rm fp}\simeq 90^{\circ}$ (black x’s), $e_{\rm b}\simeq e_{\rm b}(0)$, and trajectories above and below $I=90^{\circ}$ are symmetric. As $J$ increases in magnitude, $I_{\rm fp}$ decreases, $e_{\rm b}$ begins to oscillate, and the inclination symmetry above and below $I=90^{\circ}$ is lost. Although different trajectories may cross in the $I-\Omega$ plane, each still has a unique $\Psi$ value [Eq. (59)], since the binary’s $e_{\rm b}$ value differs from Eq. (61) when $I>\pi/2$. When $J\gg 1$, the system’s dynamics approaches the classic Lidov-Kozai regime (Lidov, 1962; Kozai, 1962). The fixed point $I_{\rm fp}$ of the system approaches Eq. (66), with $e_{\rm b}$ reaching large values when $I(0)>I_{\rm fp}$, and with trajectories symmetric above and below $I=90^{\circ}$.

Figure 7 plots $I_{\rm fp}$ as a function of $J^{-1}$, computed for $\Omega=\pi/2$ with the $e_{\rm b}=e_{\rm b}(0)$ values as indicated. The two limiting cases given by Eqs. (65) and (66) are achieved when $J\ll 1$ and $J\gg 1$, respectively, and $I_{\rm fp}$ generally varies non-monitonically with increasing $J$. Since $e_{\rm b}$ should evolve in time under the influence of viscous dissipation from disk warping, one cannot determine the final value of $I_{\rm fp}$ the system may evolve into starting from initial $I(0)$ and $e_{\rm b}(0)$ values without a detailed calculation similar to Section 3. Nevertheless, Figures 6 and 7 show there exist highly inclined fixed points for any value of $J$. For $J\lesssim 0.1$, the system may evolve into near polar alignment, with $I_{\rm fp}$ somewhat less than $90^{\circ}$.

5 Torque on Binary and Effect of Accretion

In the previous sections, we have studied the evolution of the disk around a binary with fixed $\bm{\hat{l}}_{\rm b}$ and $\bm{e}_{\rm b}$. Here we study the back-reaction torque on the binary from the disk. First consider a circular binary. The viscous back reaction torque on the binary from the disk is [Eq. (46)]

[TABLE]

In addition, accretion onto the binary from the disk adds angular momentum to the binary’s orbit:

[TABLE]

Here, $\dot{M}$ is the mass accretion rate onto the binary, $\lambda\sim 1$ (e.g. Miranda, Muñoz, & Lai 2017), and we have assumed $\bm{\hat{l}}(r_{\text{in}},t)\simeq\bm{\hat{l}}_{\rm d}(t)$ (see below). The torques (67) and (68) are equivalent to those considered in Foucart & Lai (2013), except we give different power-law prescriptions for $\Sigma=\Sigma(r)$ and $H=H(r)$. For disks in steady state, we have

[TABLE]

Using Eqs. (46) (with $e_{\rm b}=0$), (67) and (68), we obtain the net disk-binary alignment timescale for small angle between $\bm{\hat{l}}_{\rm d}$ and $\bm{\hat{l}}_{\rm b}$:

[TABLE]

where

[TABLE]

measures the strength of the accretion torque to the viscous torque on the binary ($\eta/\lambda=f^{-1}$, $\lambda=g$ in the notation of Foucart & Lai 2013).

Since $\bm{\hat{l}}(r_{\text{in}},t)\neq\bm{\hat{l}}(r_{\text{out}},t)$, the disk angular momentum loss through accretion causes $\bm{\hat{l}}_{\rm d}$ to change with time:

[TABLE]

Because the magnitude of the tilt of $\bm{\hat{l}}(r_{\text{in}},t)$ from $\bm{\hat{l}}_{\rm d}$ is of order

[TABLE]

we find

[TABLE]

Detailed calculation shows that the accretion torque (74) is always much less than the viscous torque (19) on the disk. We relegate the calculation and discussion of the accretion torque (74) to the Appendix.

For eccentric binaries, the back-reaction toque from the disk is $\text{d}{\bm{L}}_{\rm b}/\text{d}t=-\text{d}{\bm{L}}_{\rm d}/\text{d}t$ [Eq. (46)]. But this is not sufficient for determining the evolution of $\bm{e}_{\rm b}$ and $\bm{\hat{l}}_{\rm b}$. In addition, how accretion affects the binary eccentricity is also uncertain (e.g. Rafikov 2016; Miranda, Muñoz, & Lai 2017). Nevertheless, as long as $L_{\rm b}\gtrsim L_{\rm d}$, the timescale for the disk-binary inclination evolution should be of order $\gamma_{\rm b}^{-1}$, with an estimate given by Eq. (48).

6 Discussion

6.1 Theoretical Uncertainties

Our theoretical analysis of disks around binaries assumes a linear disk warp. However, we find that at the inner disk region, $|\partial\bm{\hat{l}}/\partial\ln r|$ reaches $\sim 0.1$ for a wide range of binary and disk parameters. Inclusion of weakly non-linear warps in Equations (22)-(23) may introduce new features in the disk warp profile (Ogilvie, 2006). In addition, disk warps of this magnitude may interact resonantly with inertial waves in the disk, leading to a parametric instability which may excite turbulence in the disk (Gammie, Goodman, & Ogilvie, 2000; Ogilvie & Latter, 2013). An investigation of these effects is outside the scope of this paper, but their inclusion is unlikely to change the direction of disk-binary inclination evolution (alignment vs polar alignment).

6.2 Observational Implications

In Section 3.4, we showed that the viscous torque associated with disk twist/warp tends to drive the circumbinary disk axis $\bm{\hat{l}}_{\rm d}$ toward $\pm\bm{\hat{l}}_{\rm b}$ (alignment or anti-alignment) when $\Lambda>0$, and toward $\pm\bm{e}_{\rm b}$ (polar alignment) when $\Lambda<0$. Note that $I_{\rm crit}<I<180^{\circ}-I_{\rm crit}$ is a necissary, but not sufficient condition for polar alignment of the disk [Eq. (8)]. An extreme example is when $\Omega=0^{\circ}$, since $\Lambda\geq 0$ for all inclinations $I$. Because the circumbinary disk probably formed in a turbulent molecular cloud, the disk is unlikely to have a preferred $\Omega$ when it forms. The condition for polar alignment ($\Lambda<0$) requires $\Omega$ to satisfy

[TABLE]

Assuming a uniform distribution of $\Omega$-values from [math] to $2\pi$, the probability of the disk to polar align is (for given $I,e_{\rm b}$)

[TABLE]

where

[TABLE]

We define the inclination $I_{\rm polar}$ through $P_{\rm polar}(I_{\rm polar},e_{\rm b})~{}=~{}0.5$. Solving for $I_{\rm polar}$, we obtain

[TABLE]

In Figure 8, we plot contours of constant $P_{\rm polar}$ in the $I-e_{\rm b}$ space. The $P_{\rm polar}=0$ curve (black) traces out $I_{\rm crit}$ [Eq. (8)], while the $P_{\rm polar}=0.5$ curve (red) traces out $I_{\rm polar}$ [Eq. (78)]. When $I<I_{\rm crit}$, alignment of $\bm{\hat{l}}$ with $\bm{\hat{l}}_{\rm b}$ is inevitable. When $I>I_{\rm polar}$, alignment of $\bm{\hat{l}}$ with $\bm{e}_{\rm b}$ is probable.

Table 1 lists a number of circumbinary systems with highly eccentric binaries. With the exception of 99 Herculis, all the binaries listed have disks coplanar with the binary orbital plane within a few degrees. We also list $I_{\rm crit}$ [Eq. (8)] and $I_{\rm polar}$ [Eq. (78)] for these systems. We do not list the binaries KH 15D (Winn et al., 2004; Chiang & Murray-Clay, 2004; Capelo et al., 2012) and HD 142527B (Marino et al., 2015; Casassus et al., 2015; Lacour et al., 2016) since the orbital elements of these binaries are not well constrained. However, both binaries appear to have significant eccentricities (Chiang & Murray-Clay, 2004; Lacour et al., 2016).

Since planets form in gaseous circumbinary disks, planets may form with orbital planes perpendicular to the binary orbital plane if the binary is sufficiently eccentric. Such planets may be detectable in transit surveys of eclipsing binaries due to nodal precession of the planet’s orbit.

The twist and warp calculated in Section 3.3 is non-negligible. Further observations of (gaseous) circumbinary disks may be able to detect such warps (Juhász & Facchini, 2017), further constraining the orientation and dynamics of circumbinary disk systems.

7 Summary

Using semi-analytic theory, we have studied the warp and long-term evolution of circumbinary disks around eccentric binaries. Our main results and conclusions are listed below.

For protoplanetary disks with dimensionless thickness $H/r$ larger than the viscosity parameter $\alpha$, bending wave propagation effectively couples different regions of the disk, making it precess as a quasi-rigid body. Without viscous dissipation from disk warping, the dynamics of such a disk is similar to that of a test particle around an eccentric binary (Secs. 2 and 3.2). 2. 2.

When the binary is eccentric and the disk is significantly inclined, the disk warp profile exhibits new features not seen in previous works. The disk twist [Eq. (34)] and warp [Eq. (35)] have additional contributions due to additional torques on the disk when the binary is eccentric. 3. 3.

Including the dissipative torque from warping, the disk may evolve to one of two states, depending on the initial sign of $\Lambda$ [Eq. (4)] (Sec. 3.4). When $\Lambda$ is initially positive, the disk angular momentum vector aligns (or anti-aligns) with the binary orbital angular momentum vector. When $\Lambda$ is initially negative, the disk angular momentum vector aligns with the binary eccentricity vector (polar alignment). Note that $\Lambda$ depends on both $I$ (the disk-binary inclination) and $\Omega$ (the longitude of ascending node of the disk). Thus for a given $e_{\rm b}$, the direction of inclination evolution depends not only on the initial $I(0)$, but also on the initial $\Omega(0)$. 4. 4.

When the disk has a non-negligible angular momentum compared to the binary, the system’s fixed points are modified (Sec. 4). The disk may then evolve to a state of near polar alignment, with the inclination somewhat less than $90^{\circ}$. 5. 5.

The timescale of evolution of the disk-binary inclination angle [see Eqs. (52)-(53)] depends on various disk parameters [see Eq. (48)], but is in general less than a few Myrs. This suggests that highly inclined disks and planets may exist around eccentric binaries.

Acknowledgments

JZ thanks Stephen Lubow, Rebbeca Martin, and Gongjie Li for useful conversations. We thank the anonymous referee for their comments which improved the quality of this paper. This work has been supported in part by NASA grants NNX14AG94G and NNX14AP31G, and a Simons Fellowship from the Simons Foundation. JZ is supported by a NASA Earth and Space Sciences Fellowship in Astrophysics.

Appendix: Accretion Torques

If the inner disk is not coplanar with the outer disk, accretion will change the disk angular momentum vector over time. We parameterize this accretion torque according to

[TABLE]

where $\dot{M}$ is the accretion rate onto the binary, and $\lambda\sim 1$ parameterizes the angular momentum loss from the disk to the binary. The time evolution of ${\bm{G}}$, as well as the pericenter precession induced by the non-Keplarian angular frequency, warps the inner edge of the disk in the direction of the binary orbital plane [$({\bm{l}}_{1})_{\rm warp}$; Eq. (35)]. Inserting $\bm{\hat{l}}(r,t)=\bm{\hat{l}}_{\rm d}+({\bm{l}}_{1})_{\rm warp}$ in equation (79), we obtain

[TABLE]

where

[TABLE]

We have assumed the disk to be in a steady state, so

[TABLE]

Equation (80) agrees with the rough magnitude and direction of the accretion torque estimated in Equation (74).

Since

[TABLE]

the radial functions $\tau_{\rm b}(r_{\text{in}}),W_{\rm bb}(r_{\text{in}})<0$, and $f_{\rm b}\sim\Lambda$, Eqs. (83)-(84) drives the disk one of two ways depending on the rough value of $\Lambda$:

$\Lambda\gtrsim 0$: The accretion torque (80) pushes $\bm{\hat{l}}_{\rm d}$ away from $\bm{\hat{l}}_{\rm b}$. 2. 2.

$\Lambda\lesssim 0$: The accretion torque (80) pushes $\bm{\hat{l}}_{\rm d}$ away from $\bm{e}_{\rm b}$.

From Eq. (27), it may be shown that there are no fixed points near the $\Lambda=0$ separatrix. Therefore, the accretion torque drives the disk to a trajectory near the $\Lambda=0$ separatrix.

Figure 9 plots the examples considered in Figure 2 with accretion torques [Eq. (80)]. We take $h$ and $\alpha$ to be significantly higher than our cannonical values of $\alpha=0.01$ and $h=0.1$ so that accretion torques effect the dynamical evolution of the circumbinary disk [Eq. (81)]. In the left panels of Figure 9, we plot the disk inclination with time, for the binary eccentricities indicated. The trajectories which start at $I(0)=20^{\circ},40^{\circ},$ and $80^{\circ}$ all evolve toward the prograde seperatrix, which nutates around $I\sim 50^{\circ}$ when $e_{\rm b}=0.3$, and $I\sim 40^{\circ}$ when $e_{\rm b}=0.6$. The trajectories which start at $I(0)=60^{\circ}$ both evolve to the retrograde seperatrix, which nutates around $I\sim 130^{\circ}$ when $e_{\rm b}=0.3$, and $I\sim 140^{\circ}$ when $e_{\rm b}=0.6$. On the right panels, we plot the disk trajectories on the $I-\Omega$ plane, for the binary eccentricities indicated. All disk trajectories evolve toward the $\Lambda\approx 0$ seperatrix.

The relative strength of the viscous to the accretion torques from disk warping is given by the ratio

[TABLE]

As long as $|\gamma_{\rm b}|\gg|\gamma_{\rm a}W_{\rm bb}(r_{\text{in}})|$, the viscous torque dominates, and $\bm{\hat{l}}_{\rm d}$ aligns with either $\bm{\hat{l}}_{\rm b}$ or $\bm{e}_{\rm b}$, depending on the sign of $\Lambda$ (Sec. 3.4). When $|\gamma_{\rm b}|\lesssim|\gamma_{\rm a}W_{\rm bb}(r_{\text{in}})|$, the accretion torques may dominate, and $\bm{\hat{l}}_{\rm d}$ may be driven to the seperatrix $\Lambda\approx 0$.

Figure 10 is identical to Figure 5, except we include viscous [Eq. (50)] and accretion [Eq. (80)] torques with $\alpha=0.01$, $h=0.1$, and $\lambda=1$. Because $|\gamma_{\rm b}|\gg|W_{\rm bb}(r_{\text{in}})\gamma_{\rm a}|$, the viscous torque dominates the disk’s dynamics. As a result, Figure 10 is almost indistinguishable from Figure 5. Only for unrealistically hot protoplanetary disks with $h\gtrsim 0.5$ may accretion torques significantly effect the disk evolution over viscous timescales.