The Kozai--Lidov Mechanism in Hydrodynamical Disks. II. Effects of binary and disk parameters
Wen Fu{$^{1
2
\cite{Martin14b showed that a substantially misaligned accretion disk around one component of a binary system can undergo global damped Kozai--Lidov (KL) oscillations. During these oscillations, the inclination and eccentricity of the disk are periodically exchanged. However, the robustness of this mechanism and its dependence on the system parameters were unexplored. In this paper, we use three-dimensional hydrodynamical simulations to analyze how various binary and disk parameters affect the KL mechanism in hydrodynamical disks. The simulations include the effect of gas pressure and viscosity, but ignore the effects of disk self-gravity. %consider disks with a mass that is sufficiently small that we can omiignore effects of disk self-gravity. We describe results for different numerical resolutions, binary mass ratios and orbital eccentricities, initial disk sizes, initial disk surface density profiles, disk sound speeds, and disk viscosities. We show that the KL mechanism can operate for a wide range of binary-disk parameters. We discuss the applications of our results to astrophysical disks in various accreting systems.
main.tex
\documentclass{aastex}
%\documentclass[iop]{emulateapj}
\usepackage{amsmath,amssymb,bm}
\begin{document}

\title{The Kozai--Lidov Mechanism in Hydrodynamical Disks. II. Effects of binary and disk parameters}
\author{Wen Fu{$^{1,2}$}, Stephen H. Lubow{$^3$} and Rebecca G. Martin{$^4$}}
\affil{{$^1$}Department of Physics and Astronomy, Rice University, Houston, TX 77005, USA; wf5@rice.edu\\
{$^2$}Los Alamos National Laboratory, Los Alamos, NM 87545, USA\\
{$^3$}Space Telescope Science Institute, Baltimore, MD 21218, USA\\
{$^4$}Department of Physics and Astronomy, University of Nevada, Las Vegas, Las Vegas, NV 89154, USA}


%%%%%%%%%%%%%%%%%%%%%%%%
%%%    ABSTRACT      %%%
%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract}
\cite{Martin14b} showed that a substantially misaligned accretion disk around one component of a binary system can undergo global damped Kozai--Lidov (KL) oscillations. During these  oscillations, the inclination and eccentricity of the disk are periodically exchanged.
However, the robustness of this mechanism and its dependence on the system parameters were unexplored. 
In this paper, we use three-dimensional hydrodynamical simulations to analyze how various
binary and disk parameters affect the KL mechanism in hydrodynamical disks. 
The simulations include the effect of gas pressure and viscosity, but ignore the effects
of disk self-gravity.
%consider disks with a mass that is sufficiently small that we can omiignore effects of disk self-gravity.
We describe results for different numerical resolutions, binary mass ratios and orbital eccentricities, initial disk sizes, 
initial disk surface density profiles, disk sound speeds, and disk viscosities. 
We show that the KL mechanism can operate for a wide range of binary-disk parameters. 
We discuss the applications of our results to astrophysical disks in various accreting systems. 
\end{abstract}

\keywords{accretion, accretion disks -- binaries: general -- hydrodynamics -- planets and satellites: formation}
%%%%%%%%%%%%%%%%%%%%%%%%
%%% INTRODUCTION     %%%
%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

Disks in binary systems may be misaligned with respect
to their binary orbit planes. Alignment is expected
in cases where the disk formation process itself
is coplanar, such as the case of mass-exchange binaries.
In that case, a coplanar gas stream originates from the inner Lagrange
point of the mass-losing star and results in the formation of a coplanar disk.
Coplanarity is expected if the binary forms from disk
fragmentation \citep{Bonnell94, Offner10}. Furthermore, coplanarity is the long term
expected outcome of a misaligned disk, due to the effects of tidal dissipation
in the disk \citep{Papaloizou95, Bate00, Lubow00, King13}.

On the other hand, disk misalignment is expected to persist for some time period in
cases where the binary forms from a merger process,
such as occurs with supermassive black hole (SMBH) binary systems   \citep[e.g.,][]{King06}. 
If a young binary star system accretes material after its
formation process, the material is likely to be misaligned with the binary orbit and so misaligned disks may form around young stars \citep{Bate10}. In addition, a coplanar disk
may become noncoplanar due to an instability, such as
radiation warping \citep{Pringle96, Wijers99, Ogilvie01b}. Be/X-ray binaries  may have transient disks that are misaligned with respect 
to the binary orbit \citep{Martin09,Martin11}. These disks may be produced
by the expulsion of equatorial material 
from a rapidly rotating Be star whose spin axis is misaligned
with the binary rotation axis \citep[e.g.,][]{Porter03}.


There is both direct and indirect observational evidence for 
disk noncoplanarity. For widely separated stars in a binary, greater than 40 AU, 
a misaligned disk may occur because the stellar equatorial inclinations, based on spins, are observationally inferred to be misaligned with respect to the binary orbital planes \citep{Hale94}. 
The young binary system HK Tau provides
direct evidence for noncoplanarity, since disks are observed around
both components with one disk edge-on and the other more face-on
\citep{Stapelfeldt98}. Recent ALMA observations by \cite{Jensen14}
suggest that its disks are mutually misaligned by more than $60^{\circ}$, although
the plane of the binary orbit is not known. Strong mutual disk misalignment ($\sim 72^{\circ}$) has lately also been detected for the two circumstellar disks in V2434 Ori, a binary system in Orion \citep{Williams14}.
Less direct evidence comes from the existence
of extra-solar planets whose orbits are tilted with respect to the spin 
axis of the central star \citep{Albrecht12, Winn14}.  If the planets reside in binary star systems, this evidence suggests
that these planets may have formed in disks that are misaligned with the binary orbital plane \citep[e.g.,][]{Bate10, Batygin12}.

Test particles whose orbits are sufficiently inclined with respect to the plane of a circular orbit binary can undergo Kozai-Lidov (KL) oscillations \citep{Kozai62, Lidov62}. In these oscillations, the particle's orbital inclination and eccentricity evolve in such a manner that inclination is exchanged for eccentricity. Due to this process,  a test particle orbit that is initially circular can achieve a high eccentricity when it evolves to smaller inclination. For example, an initially circular  orbit at an inclination of $60^{\circ}$ achieves an eccentricity of about 0.75. Further extensions of this theory show that even stronger effects can occur for eccentric orbit binaries \citep{Ford00, Lithwick11, Naoz11, Naoz13a, Naoz13b, Teyssandier13, Li14, Liu15}.

The KL effect for ballistic objects has been applied to a wide range of astronomical processes. These include inclinations of asteroids and irregular satellites \citep{Kozai62, Nesvorny03},  high orbital eccentricities of some extra-solar planets and formation of hot Jupiters \citep{Holman97, Mazeh97, Wu03, Takeda05, Xiang13,  Dawson14, Dong14, Petrovich15, Rice15}, misalignment between an exoplanet's orbital axis and stellar spin axis \citep{Storch14, Storch15}, black hole mergers \citep{Blaes02, Miller02},  the formation of close binary stars \citep{Harrington68, Mazeh79, Kiseleva98, Fabrycky07}, tidal disruption events \citep{Chen11}, and the formation of Type Ia supernovae \citep{Kushnir13}.



In a recent {\it Letter} \citep[hereafter Paper I]{Martin14b}, we showed in our Smoothed Particle Hydrodynamics (SPH) simulations 
that a fluid disk with pressure and viscosity that orbits a member of a binary
can also undergo global KL oscillations.
%We described a few SPH simulations that showed that such a disk can undergo global
%KL oscillations.
In this paper, we explore this process in more detail by considering how it is affected
by various disk and binary parameters. Following the approach of Paper I,
we attempt to relate the properties of disk KL oscillations to the properties
of KL particle oscillations.

The outline of the paper is as follows. In Section \ref{sec:part}
we describe the properties of particle orbits in binaries.
The parameters of the particle orbits and the binaries are chosen to be similar to those of the disk simulations.
In Section \ref{sec:kld} we describe the results of SPH simulations
for a range of initial particle numbers, disk viscosities,
disk aspect ratios, initial density profiles (including changes to 
initial inner and outer disk radii),  initial disk inclinations, binary
mass ratios, and binary eccentricities.
Section \ref{sec:sum} contains the discussion and summary.


%%%%%%%%%%%%%%%%%%%%%%%%
%%% KL of a particle %%%
%%%%%%%%%%%%%%%%%%%%%%%%

\section{KL cycle of a test particle \label{sec:part}}
As discussed in Paper I, the properties of a fluid disk undergoing KL oscillations are related
to the properties of test (massless) particles undergoing KL oscillations that orbit at similar radii. 
In both the disk and test particle cases, the inclination and eccentricity are periodically exchanged. However, there are some
fundamental differences.
The period of a test particle oscillation varies with distance from the central object.
In the case of a disk
that is subject to pressure and viscous forces,
the oscillations are global (i.e., the period is the same over all radii)  and subject to damping by dissipation, unlike the test particle case. % But the amplitude range of the initial inclination and eccentricity oscillations 
%are similar in the two cases. 
The global disk oscillation frequency is expected to involve an angular momentum weighted average
of the test particle oscillation frequencies (Equation (4) of Paper I).
Before considering the disk oscillations, we describe in this section the behavior
of test particle orbits. We apply 
binary parameters that are similar to those
we adopt for the disk simulations described in Section \ref{sec:kld}.
 
We consider ballistic particle orbits around one component of a binary system. 
Each particle orbit is initially circular, but substantially inclined 
with respect to the orbital plane of the binary. If the binary orbit is circular, then
over long timescales, the binary companion  effectively acts as a uniform circular ring.
Since this effective potential is time independent,
the energy of the particle about the central mass is conserved. Thus, the semi-major axis $a$ of the particle in a binary is nearly independent of time $t$, that is
\begin{equation}
a(t)  \approx {\rm constant},
\label{eq:adot}
\end{equation}
where the constant is determined by initial conditions.
Because the time-averaged perturbing potential is axisymmetric, 
the component of the  angular momentum of the particle that is perpendicular to the binary orbital plane is approximately conserved over long
timescales. From this conservation principle and Equation (\ref{eq:adot})
we have that
\begin{equation}
\sqrt{1-e^2(t)} \, \cos{\left(i(t) \right)} \approx {\rm constant},
\label{eq:e-i}
\end{equation}
where the constant is determined by initial conditions, $e$  is the orbital eccentricity of the particle, and $i$ is the inclination of the particle orbit with respect to the binary orbital plane.

We denote the masses of the binary component objects as $M_{\rm c}$ and $M_{\rm p}$, where $M_{\rm c}$ ($M_{\rm p}$)  is the mass of the central (perturbing) object to the particle orbit.
Analytic calculations have determined that the KL oscillation period for a particle subject to the effect of a binary in an eccentric orbit is approximately given by
\begin{equation}
\frac{\tau_{\rm KL}}{P_{\rm b}}\approx \frac{M_{\rm c}+M_{\rm p}}{M_{\rm p}}\frac{P_{\rm b}}{P}(1-e^2_{\rm b})^{3/2}\label{eq:klperiod}
\end{equation}
\citep{Holman97, Innanen97, Kiseleva98}, where $P=\sqrt{G M_{\rm c}/a^3}$ is the orbital period of the particle with orbital semimajor axis $a$,  $P_{\rm b}=2\pi/\Omega_{\rm b}$ is the orbital period of the binary,  
$e_{\rm b}$ is the eccentricity of the binary orbit, and
$\Omega_{\rm b}=\sqrt{G(M_{\rm c}+M_{\rm p})/a_{\rm b}^3}$ is the binary orbital frequency for binary semi-major axis $a_{\rm b}$. 

For the KL oscillations to occur, the initial inclination of the test particle orbit $i_0$ must satisfy the condition that $\cos^2{(i_0)} < \cos^2{(i_{\rm cr})}  = 3/5$,
where $i_{\rm cr}$ is the critical angle for KL oscillations.
This condition requires that $39^\circ \la i_0 \la 141^\circ$.
During the KL cycle, the  inclination $i$ of an initially circular prograde particle orbit ($i_0 < 90^\circ$) oscillates between the initial value $i_0$ and the
critical angle $i_{\rm cr} \simeq 39^{\circ}$, while the eccentricity oscillates between $e_{\rm min}=0$ and the maximum value of
\begin{equation}
e_{\rm max}=\sqrt{1-\frac{5}{3}{\rm cos}^2(i_0}),
\label{eq:maxe}
\end{equation}
that is achieved for $i=i_{\rm cr}$,
which follows from Equation (\ref{eq:e-i}) with the right-hand side constant value of 
$\cos{(i_0)}$
\citep[e.g.,][]{Innanen97}.


For an eccentric binary orbit,  the vertical component of a particle's angular
momentum is not conserved. Equation (\ref{eq:e-i}) breaks down over long timescales and the maximum eccentricity of the particle 
can approach unity \citep{Ford00, Lithwick11, Naoz11, Naoz13a, Naoz13b, Liu15}.
The orbit can flip its orientation from prograde to retrograde 
with respect to the binary. This flip occurs as the orbit passes through a radial state ($e=1$). 
In passing through a nearly radial orbit, the particle may collide with the central object.


In our numerical calculations, we compute the eccentricity of a test particle with coordinates $\bm{r}(t)$ by considering its
specific angular momentum $\bm{j}(t)$  relative to the disk central object at position $\bm{r}_{\rm c}(t)$,
\begin{equation}
\bm{j}(t) = ( \bm{r} -\bm{r}_{\rm c})\bm{\times}(\bm{\dot r}-\bm{\dot r}_{\rm c}),
\label{j}
\end{equation}
where the dot denotes differentiation in time.
The relative specific energy of the particle is
\begin{equation}
E(t) = \frac{1}{2}|\bm{\dot r}-\bm{\dot r}_{\rm c}|^2-\frac{G M_{\rm c}}{|\bm{r}-\bm{r_c}|}
\end{equation}
and thus its eccentricity is
\begin{equation}
e(t) =\sqrt{1+\frac{2 E |\bm{j}|^2}{(GM_{\rm c})^2}} .
\label{eq:e}
\end{equation}
The evolving inclination of the particle is
\begin{equation}
i(t) =\arccos \left( \frac{j_z}{|\bm{j}|}\right),
\label{eq:i}
\end{equation}
where $j_z$ is the component of the angular momentum in
Equation~(\ref{j}) that is perpendicular to the orbital plane of the binary.

\begin{figure}
\centering
\includegraphics[width=0.95\textwidth]{testpart.eps}
\caption{Eccentricity (upper row) and inclination (bottom row) evolution of a test particle for different values of the initial orbital inclination $i_0$, initial orbital radius $r_0$, binary mass ratio $M_{\rm p}/M_{\rm c}$, and binary eccentricity $e_{\rm b}$. The binary orbit has a semi-major axis of $a_{\rm b}$. Time is in units of binary orbital period $P_{\rm b}$. These results are obtained from three-body numerical simulations. 
 (Color online) \label{fig:testpart}}
\end{figure}

Figure \ref{fig:testpart} plots the numerical results for a test particle undergoing KL oscillations from simple three-body simulations. 
%Here the binary is initially on circular orbit ($e_{\rm b}=0$). 
The upper and lower rows show the evolution of the particle's orbital eccentricity and inclination, respectively. 
We investigate the effects of 
changes to four system parameters: the initial orbit inclination $i_0$, the initial orbital radius of the particle $r_0$, the binary mass ratio $M_{\rm p}/M_{\rm c}$, and
the binary eccentricity $e_{\rm b}$.  In all cases the initial orbital eccentricity of the particle is taken to be zero, $e_0=0$.
Each column of the figure corresponds to runs where we vary one of these parameters and fix the other three at fiducial (reference standard)
values. 
The fiducial parameters that we adopt for this figure 
are  $i_0 = 60^\circ, r_0 = 0.2 a_{\rm b}, e_{\rm b}=0,$ and  $M_{\rm p}/M_{\rm c}=1$.  
The fiducial value of $r_0$ is chosen as a characteristic radius 
that typically lies within the disks that we simulate  in Section 
\ref{sec:kld}. 



The left hand panels show that when the orbit of the particle is highly inclined (well above critical KL angle 
$i_{\rm cr}$; see the green and blue lines), the KL cycle period 
is relatively insensitive to inclination, as expected by Equation (\ref{eq:klperiod}) that has no dependence on inclination.  The simulated period is sometimes close to that
predicted by Equation (\ref{eq:klperiod}). 
 For $i_0=70^\circ$ (blue line), the eccentricity oscillation has $e_{\rm max}\simeq 0.93$ and $\tau_{\rm KL}\simeq15 P_{\rm b}$, while Equations (\ref{eq:klperiod}) and (\ref{eq:maxe}) 
 predict that $e_{\rm max}=0.90$ and $\tau_{\rm KL}=16 P_{\rm b}$. 
 When the initial orbital tilt is just above the critical KL angle, the KL cycle period given by Equation (\ref{eq:klperiod}) becomes less accurate. 
 For $i_{0}=45^{\circ}$ (red line), the period from the simulation is $\tau_{\rm KL}\simeq 33 P_{\rm b}$, which is more than twice the analytical value. %Therefore, even though Equation (\ref{eq:klperiod}) does not contain  a dependence on initial inclination that is found in simulations, it is in general valid up to a factor of a few. 
 
 The second column shows the effect of changing the initial orbital radius of the particle $r_0$. The KL periods for $r_0=0.15 a_{\rm b}$, $0.2 a_{\rm b}$, and $0.25 a_{\rm b}$ are $\sim 30 P_{\rm b}$,
  $\sim 15 P_{\rm b}$, 
 and $\sim 9 P_{\rm b}$, respectively. These values agree well (within $\sim 25\%$) with  the predictions of Equation (\ref{eq:klperiod}) that scales as $\tau_{\rm KL} \propto r^{-3/2}_0$ and produces KL periods of $24 P_{\rm b}$, $16 P_{\rm b},$ and $11 P_{\rm b}$, respectively.
 Note that with $r_0=0.25 a_{\rm b}$, the KL cycle period becomes shorter and the oscillation amplitude is reduced  at later time ($e_{\rm min}$ does not  reach zero). This is likely because at this orbital radius, the perturbing potential from the companion is more significant.
 As a result, the stronger influence makes the analytic description as a secular (gradual) effect less accurate and causes the oscillations to be irregular. 
 
 The effect of varying the binary mass ratio is shown in the third column.
  Equation (\ref{eq:klperiod}) predicts that $\tau_{\rm KL}/P_{\rm b} \propto (M_{\rm c}/M_{\rm p}+1)^{1/2}(M_{\rm c}/M_{\rm p})^{1/2}$ such that $\tau_{\rm KL}=10P_{\rm b}$, $16P_{\rm b}$, and $27P_{\rm b}$ for binary mass ratios $M_{\rm p}/M_{\rm c}=2$, $1$, and $1/2$, respectively. The results of the numerical simulations give us $\tau_{\rm KL}= 7P_{\rm b}$, $15P_{\rm b}$, and $32P_{\rm b}$. 
  The agreement is then within about 30\%.
  
 In the fourth column, we show the effect of varying the binary eccentricity. %for a fixed value of binary semi-major axis $a_{\rm b}$. CHANGE to $i_0=60^\circ$.
 %The value of initial inclination is taken to be $i_0=50^\circ$, rather than the fiducial value used above of $60^\circ$. The other quantities are set to the fiducial values adopted for the first three columns. This lower value of $i_0$ was used because the effects of the KL oscillations become stronger
%with increasing binary eccentricity. By reducing the initial inclination, we obtain milder orbital changes with increased binary eccentricity, which is important when we compare these particle results with disk simulations in Section \ref{sec:kld}. 
As discussed above, previous studies have shown that for eccentric orbit binaries,
 the maximum eccentricity can grow beyond $e_{\rm max}$ given by Equation~(\ref{eq:maxe}) and the maximum inclination
 can grow above the initial value of $i_0$.
As expected, for a somewhat higher binary eccentricity $e_{\rm b}=0.2$ (green lines), the maximum and minimum particle eccentricity sometimes increases somewhat
in each KL oscillation. %The maximum inclination increases while the minimum inclination is approximately constant.
The oscillation period is shorter than in the circular case, 
as predicted by Equation  (\ref{eq:klperiod}), since the companion object interacts more strongly with the particle at binary
periastron. According to Equation (\ref{eq:klperiod}), for $e_{\rm b}=0.2$, the KL oscillation period is approximately $14P_{\rm b}$, whereas the numerically determined period in Figure \ref{fig:testpart} is about $10P_{\rm b}$. Therefore, the reduction in the KL oscillation period due to binary orbital eccentricity is somewhat more severe than what is predicted analytically . 




%%%%%%%%%%%%%%%%%%%%%%%%
%%% KL of a disk     %%%
%%%%%%%%%%%%%%%%%%%%%%%%

\section{KL cycle of a hydrodynamic disk \label{sec:kld}}




In this section, we describe the KL cycles of three-dimensional hydrodynamic disks. We use the SPH code  \texttt{PHANTOM} \citep{Lodato10, Price10, Nixon12a, Nixon12b,Price12,  Nixon13}.  The disk is initially circular, centered around the binary component of mass $M_{\rm c}$, and subject to perturbations by the other binary component 
of mass $M_{\rm p}$. The disk orbital plane is initially misaligned with respect to the binary orbital plane.  
%Particles that enter a region at smaller radii than the inner boundary of radius $r_{\mathrm{in}}$ are removed from the simulations. 
The disk initially extends from radius $r_{\mathrm{in}}$ to radius
$r_{\mathrm{out}}$. The inner boundary of the simulated region is set to the disk initial inner disk radius 
$r_{\mathrm{in}}$. As particles move to $r \le r_{\mathrm{in}}$, they are removed from the simulation. 
%This inner radius is also called the accretion radius.
We also impose an inner boundary radius around the perturbing companion, since some disk mass can be transferred to that component. The \texttt{PHANTOM} code adopts cubic spline kernel as the smoothing kernel. The softening length in the gravitational force is taken to be the same as the smoothing length. The number of neighbors is roughly constant at $N_{\rm neigh} \approx 58$. 
%In equal mass binary systems, we keep accretion radii for both stars the same. In binary systems with unequal masses, 
%The inner boundary radius around that component is set to $r_{\mathrm{in}} (M_{\rm p}/M_{\rm c})^{1/3}$.
%The initial disk outer radius $r_{\rm{out}}$ is chosen to be the tidal truncation radius for the disk assuming a coplanar binary. 



We define a certain set of fiducial (reference) model parameters. The disk is locally isothermal  with sound speed $c_{\rm s} \propto r^{-3/4}$ and disk aspect ratio $H/r = 0.035$ at the initial inner disk radius $r_{\mathrm{in}}$. 
These parameters allow both $\alpha$ and the smoothing length $\left<h\right> /H$ to be constant throughout the disk radius \citep{Lodato07}. 
We employ an explicit accretion disk viscosity which corresponds to an approximately constant  
$\alpha$ parameter \citep{SS73} over the disk \citep[see][]{Lodato10}.
The viscous stresses include a nonlinear term with a coefficient $\beta_{\rm AV} = 2$ (AV stands for artificial viscosity) that suppresses interparticle penetration, as is standard in SPH codes. 
The fiducial binary mass ratio is unity and the binary orbital eccentricity is zero (circular orbit).
%The fiducial disk model is resolved with shell-averaged smoothing length per scale height $\left<h\right> /H = ?$.
This fiducial setup is similar to  that presented in Paper I. The only difference is that the initial number of particles in the fiducial model is set to $N=3\times10^5$, while Paper I
used $N=1 \times 10^6$. The lower value adopted here reduces computational overhead in carrying out the various simulations, while providing results that are well converged in $N$, as we discuss below.

\begin{deluxetable}{l c r r}
\tabletypesize{\footnotesize}
\tablecolumns{4}
\tablewidth{0pt}
\tablecaption{Parameters of the SPH simulations for Binary Systems with Semi-major Axis of $a_{\rm b}$\label{table:parameters}}
\tablehead{
\colhead{Binary and Disk parameters} & 
\colhead{Symbol} & 
\colhead{Fiducial Value} & 
\colhead{Values} 
}
\startdata
Mass ratio of binary components & $M_{\rm p}/M_{\rm c}$ & 1 & [0.25, 0.5, 1, 2]\\
Binary orbital eccentricity & $e_{\rm b}$ & 0 & [0, 0.2, 0.5]\\
Initial number of particles & $N/10^5$ & 3 & [2, 3, 5, 10]\\
Initial disk mass & $M_{\mathrm{disk}}/(M_{\rm c}+M_{\rm p})$ & 0.001 & 0.001 \\
Initial disk outer radius & $r_{\mathrm{out}}/a_{\rm b}$ & 0.25 & [0.15, 0.25]\\
Initial disk inner radius & $r_{\mathrm{in}}/a_{\rm b}$ & 0.025 & [0.015, 0.02, 0.025]\\
Mass accretion radius & $r_{\mathrm{acc}}/a_{\rm b}$ & 0.025 & [0.015, 0.02, 0.025]\\
Disk viscosity parameter & $\alpha$ & 0.1 & [0.01, 0.1]\\
Disk aspect ratio & $H/r\, (r=r_{\mathrm{in}})$ & 0.035 & [0.02, 0.035, 0.05, 0.065]\\
Initial disk surface density profile $\Sigma \propto r^{-\gamma}$ & $\gamma$ & 1.5 & [0.5, 1.0, 1.5]\\
Initial disk inclination & $i_{\rm 0}$ & $\mathrm{60^{\circ}}$ & [$\mathrm{45^{\circ}}$, $\mathrm{50^{\circ}}$, $\mathrm{55^{\circ}}$,$\mathrm{60^{\circ}}$]\\
\enddata
%\tablenotetext{a}{With different mass ratios, we also adjust the inner and outer disk radii accordingly.}
\label{tab:params}
\end{deluxetable}

Table \ref{tab:params} lists the simulation parameters.  We consider several models whose parameters deviate  from the fiducial model, but not all combinations of parameters were simulated.
Instead, we consider variations of one parameter with all others fixed at the fiducial values, as we did for particle orbits in Figure \ref{fig:testpart}.
For different mass ratios, we also adjust the initial inner and outer disk radii.
We scale the initial disk inner radius from the central binary component by a factor of $(M_{\rm c}/(M_{\rm c}+M_{\rm p}))^{1/3}$. As discussed above, the accretion radius of the central object is taken to be the same as the initial inner disk radius. The inner boundary condition is that any SPH particle that goes inside accretion radius is removed from the simulation. 
The inner boundary radius of the perturbing component is set to $r_{\mathrm{in}} (M_{\rm p}/M_{\rm c})^{1/3}$.
The initial disk outer radius $r_{\mathrm{out}}$ is chosen to be that of a tidally truncated disk in a coplanar binary \citep{Paczynski77}. 
Thus, when we change the mass ratio of the binary, we also change the initial disk outer radius. 
However, in a misaligned binary the tidal torque on the disk is weakened and the disk can expand somewhat beyond this radius \citep{Lubow15}. 
%We always set $R_{\mathrm{acc}}=R_{\mathrm{in}}$ in our runs.
We ignore the effects of disk self-gravity in the simulations presented here. That means the disk does not feel its own gravity, but the stars can feel it.
 We set the ratio of the initial disk mass to the binary total mass to a very low value
of $0.1\%$ such that the binary orbit is hardly affected by disk's gravity.
This mass corresponds to a large initial Toomre Q of $\sim 30$. We will devote a separate paper to studying the effects of higher disk mass and stronger self-gravity.


In post-processing the simulations, we compute the  orbital eccentricity and inclination for each particle using its position and velocity information according to Equations (\ref{eq:e})
and (\ref{eq:i}). We divide the disk into 100 radial bins and calculate the mean properties of the particles within each bin, such as the inclination and eccentricity. The properties of the particles
in these radial bins are not fully independent of each other
because the disk is generally eccentric. Particles residing one bin 
follow streamlines that cover other nearby bins.
Therefore, some correlation of mean properties 
occurs in nearby bins. %A more accurate but complicated approach would be to use eccentric bins that vary in time.
In most of our analysis, we apply fixed radial bins that are fairly widely separated to reduce these correlation effects.
An alternative approach is to apply bins based on particle semi-major axis $a$. We have found that binning in $r$ and $a$ give
similar results for most of our analysis and so we generally report results with bins in $r$
for simplicity.
For the determination of disk warping, however, we adopt bins based on particle semi-major axis, as described later. 


The KL period of a test particle (Equation (\ref{eq:klperiod})) depends on the particle's orbital radius. Following Paper I, we determine an analytic estimate of the KL global period of a disk by taking an angular momentum weighted average involving the local KL period of particle orbits
as follows
\begin{equation}
\left < \tau_{\rm KL} \right > \approx \frac{\int_{r_{\rm in}}^{r_{\rm out}}\Sigma r^3 \sqrt{\frac{GM_{\rm c}}{r^3}} dr}{\int_{r_{\rm in}}^{r_{\rm out}}\tau_{\rm KL}^{-1}\Sigma r^3 \sqrt{\frac{GM_{\rm c}}{r^3}} dr},
\label{eq:klperiod2}
\end{equation}
where $\tau_{\rm KL}$ in the denominator is given by Equation (\ref{eq:klperiod}).
This estimate is equivalent to determining the disk global precession frequency as the integrated torque divided by the total disk angular momentum.
This estimate is crude because we apply the angular momentum of circular
orbits, while the disk is generally eccentric. In addition, the analytic particle orbit period
in Equation  (\ref{eq:klperiod}) is itself somewhat approximate, as we found in Section 
 \ref{sec:part}.








Figure \ref{fig:resolution} shows the evolution of the disk eccentricity and inclination at three different radii in the disk for various initial numbers of SPH particles. In this figure and following similar figures, disk quantities are taken from three individual radial bins spanning $0.1$, $0.2$ and $0.3 a_b$, respectively. Also, we plot results up to a time of 
$t=28 P_{\rm b}$. At this time, the disk has typically lost $90\%$ of its initial particles. Most of the lost particles are accreted onto the central object. A small fraction of the particles are ejected from the disk, typically less than 10\%. Most of ejected particles are accreted onto the companion and others remain in a circumbinary orbit.
From Figure \ref{fig:resolution}, we see that the effect of the numerical resolution on the disk KL mechanism is in general quite small, especially for the first two cycles. 
%At later times, the lowest resolution simulation has fewer particles remaining and thus the fluctuations becomes larger (spikes in red curves of the left column). 
Although the cycle period at later times in the low resolution run is slightly shorter than that in high resolution run, all the curves damp to nearly the same level of eccentricity $e \simeq 0.2$ and inclination $i \simeq 28^\circ$  at  a time $t=28 P_{\rm b}$. Given this insensitivity of the overall disk KL cycle pattern to the numerical resolution, we adopt $N=3 \times 10^5$ as the fiducial value
that we use throughout the rest of the paper. 
%This value is lower than that used in paper I ($N=1\times10^6$), but it allows us to finish the survey of a large parameter space within a reasonable  time-frame. 




Figure \ref{fig:resolution} shows some differences from the test particle simulation results in Figure \ref{fig:testpart}. 
We consider how
Equations (\ref{eq:adot}) and (\ref{eq:e-i}) might apply to a disk.
We may expect $a$ and $e$ to represent a measure of the disk's
global semi-major axis and eccentricity, respectively, at some representative radius in the disk, such as the $r=0.2 a_{\rm b}$  that we applied in 
Figure \ref{fig:testpart}.
Locally, a disk can transport angular momentum through waves and viscous
torques. However, globally the vertical (perpendicular to the binary orbit plane) component of the disk angular momentum is conserved in the presence of an axisymmetric ring, apart from losses due to particles leaving the disk.  On the other hand, the disk can dissipate energy, unlike
the case of test particles. Therefore, the evolution of the disk is expected
to depart from that of test particles, as given by Equations (\ref{eq:adot}) and (\ref{eq:e-i}). In Figure \ref{fig:resolution}, we
see that the KL oscillations in the fluid case
undergo damping in which the maxima of eccentricity decrease in time. In addition, 
the eccentricity minima do not drop to zero as in the particle case. The inclination 
drops below the minimum tilt of $39^\circ$ for test particle oscillations. After the KL
oscillations substantially damp, there is a significant residual disk eccentricity $e \sim 0.2$. %The estimated disk oscillation period for the fiducial model according to Equation~(\ref{eq:klperiod2}) is given below in the discussion of Figure~\ref{fig:viscosity}.





The disk viscous timescale is estimated from Equation
($18'$) of \cite{Lynden-Bell74} 
as $t_v \sim r^2/(12 \nu) \sim 20 P_{\rm b}$
for the fiducial model at $r=0.25 a_{\rm b}$.
We have also verified that the disk density evolution of the fiducial model roughly follows the
results obtained from of a 1D viscous disk model (omits pressure effects)
\cite[e.g.,][]{Lynden-Bell74} for the same initial disk profile and level of viscosity.  Therefore, over the course of the simulations, we
expect the density profile to evolve from the initial one.

The first peak value of the eccentricity seen in the upper row of the plot $e \simeq 0.70$ is approximately consistent with the test particle results 
%at a radius $r_0=0.2a_{\rm b}$
shown in Figure \ref{fig:testpart} for maximum $e \simeq 0.78$ with an initial tilt of $60^\circ$, although the test particles have a somewhat shorter KL period. %ADD INFORMATION HERE: The test particle KL period is  $\sim xx P_{\rm b}$, while
%the disk KL period is  $\sim xx P_{\rm b}$.

%The first KL cycle occurs at an slightly earlier time at $r=0.2$ than $r=0.1 a$. Nevertheless, the KL cycle at $d=0.1a$ manages to quickly catch up with the cycle at $d=0.3a$ and peak at around the same time ($T\simeq 11 P_{\rm b}$). 

\begin{figure}
\centering
\includegraphics[width=0.95\textwidth]{effect_resolution.eps}
\caption{Eccentricity (upper row) and inclination (bottom row) evolution of the disk at three different radii from the central star for different initial numbers of particles. $a_{\rm b}$ is
the binary semimajor axis and $P_{\rm b}$ is the binary orbital period. The other parameter values are given in the third column of Table \ref{tab:params}.
(color online) \label{fig:resolution}}
\end{figure}



\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{effect_viscosity.eps}
\includegraphics[width=0.4\textwidth]{T_kl_viscosity.eps}
%\includegraphics[width=0.45\textwidth]{sigma_viscosity_11.eps}
\caption{The upper figure is similar to Figure \ref{fig:resolution}, but showing the effect of varying the disk viscosity parameter $\alpha$. The lower figure shows the estimated disk KL period as a function of time (according to Equation (\ref{eq:klperiod2})). (color online) \label{fig:viscosity}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.95 \textwidth]{edge_on.eps}
\caption{View of the fiducial disk toward the $x-z$ plane at times of 0 (left) and 13 $P_{\rm b}$ (right). The binary orbit is in the $x-y$ plane (i.e., the perturbing object moves into and out of the page). 
The central mass is denoted by the black dot. The color coding is for the logarithm base 10 of the column density (i.e., density integrated along the line of sight) in
units of $(M_{\rm c}+M_{\rm p})/a^2_{\rm b}$.
 (Color online) \label{fig:edge_on}}
\end{figure}

%\begin{figure}
%\centering
%\includegraphics[width=0.48\textwidth]{warp_diff.eps}
%\includegraphics[width=0.48\textwidth]{tilt_diff.eps}
%\caption{Left panel: The magnitude of the difference between the 
%average disk angular momentum unit vectors for particles in narrow bins centered on $r=0.1a_{\rm b}$ and $r=0.2a_{\rm b}$ 
%plotted as a function of time for runs in Figure \ref{fig:viscosity}.
%Right panel: The difference between disk tilts for particles in narrow bins centered on $r=0.1a_{\rm b}$ and $r=0.2a_{\rm b}$
%plotted as a function of time for runs in Figure \ref{fig:viscosity}.
%\label{fig:warp_tilt}}
%\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.48\textwidth]{warp_diff_a_aver.eps}
\includegraphics[width=0.48\textwidth]{tilt_diff_a_aver.eps}
\caption{Left panel: the magnitude of the difference between the 
average disk angular momentum unit vectors for particles that lie in narrow bins
of semi-major axis values  centered on $a=0.1a_{\rm b}$ and $a=0.2a_{\rm b}$, averaged in time over one binary orbit, 
plotted as a function of time for runs in Figure \ref{fig:viscosity}.  
Right panel: the difference between average disk tilts for particles in narrow bins
of semi-major axis values  centered on $a=0.1a_{\rm b}$ and $a=0.2a_{\rm b}$, averaged in time over one binary orbit,  
plotted as a function of time for runs in Figure \ref{fig:viscosity}.
\label{fig:warp_tilt_a}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{effect_hoR.eps}
\includegraphics[width=0.4\textwidth]{T_kl_hoR.eps}
\caption{Similar to Figure \ref{fig:viscosity}, but showing the effect of varying the initial disk aspect ratio. (color online) \label{fig:hoR}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{effect_gamma.eps}
\includegraphics[width=0.4\textwidth]{T_kl_gamma.eps}
\includegraphics[width=0.45\textwidth]{sigma_gamma.eps}
\caption{Similar to Figure \ref{fig:viscosity}, but showing the effect of varying the initial disk surface density profile. The lower right figure plots the azimuthally averaged surface density profiles at $t=0P_{\rm b}$ and $t=11P_{\rm b}$. 
The surface density is in units of $(M_{\rm c}+M_{\rm p})/a_{\rm b}^2$
and $r$ is in units of $a_{\rm b}$.
(color online)\label{fig:gamma}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{effect_rin.eps}
\includegraphics[width=0.4\textwidth]{T_kl_rin.eps}
\caption{Similar to Figure \ref{fig:viscosity}, but showing the effect of varying the initial disk inner radius. (color online)\label{fig:rin}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{effect_rout.eps}
\includegraphics[width=0.4\textwidth]{T_kl_rout.eps}
\caption{Similar to Figure \ref{fig:viscosity}, but showing the effect of varying the initial disk outer radius. (color online) \label{fig:rout}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{effect_inc.eps}
%\includegraphics[width=0.4\textwidth]{T_kl_inc.eps}
\caption{Similar to Figure \ref{fig:viscosity}, but showing the effect of varying the initial disk inclination angle. (color online) \label{fig:inc}}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{effect_mratio.eps}
\includegraphics[width=0.4\textwidth]{T_kl_mratio.eps}
\caption{Similar to Figure \ref{fig:viscosity}, but showing the effect of varying the binary mass ratio. Note that here we show results up to $t=40P_{\rm b}$ in the upper figure. (color online) \label{fig:mratio}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.95\textwidth]{effect_eb.eps}
\caption{Similar to Figure \ref{fig:resolution}, but showing the effect of varying the binary eccentricity. (color online) \label{fig:eb}}
\end{figure}





Figure \ref{fig:viscosity} shows the effects of the disk viscosity on the KL mechanism. For lower viscosity (red curve; $\alpha=0.01$), the time at which the first eccentricity cycle peaks is basically unaffected. 
%HOW DIFFERENT ARE THE DENSITY DISTRIBUTIONS AT T =10Pb? 
However, both the eccentricity and the inclination damp at a lower rate because the disk is less dissipative. For $\alpha=0.01$ at later times, the midpoints of the inclination KL oscillations
settle to a tilt that is  near the critical value of $39^\circ$ for KL particle oscillations. Note that there is a substantial remnant disk eccentricity $e \sim 0.3$
after the KL oscillations significantly damp.
The bottom panel of the figure plots the expected KL oscillation periods based on Equation (\ref{eq:klperiod2}). The estimated period is $\sim 15 P_{\rm b}$ that
is to be compared with the period obtained in the SPH simulations (upper panel)
of $\sim 7 P_{\rm b}$. The agreement is crude (factor of 2), likely for reasons
noted below Equation  (\ref{eq:klperiod2}).
%The disk tilt is expected to damp by two effects.
%The first is due  is due to the standard viscous disk
%tilt decay that would occur in the absence of KL oscillations \citep[see e.g.][]{King13}.
%The second is due to enhanced dissipation caused by KL oscillations. 
%We see from the figure that KL oscillations stop once the mean disk tilt  decays to
%the critical KL angle [{\bf It seems the disk tilt decays to BELOW critical angle?}]. 
%For case of the higher value of $\alpha=0.1$ in Figure
%The standard tilt decay plays some role in the tilt decay below this angle that is more pronounced in the higher
%$\alpha$ case. The lower $\alpha=0.01$ case appears to be dominated by
%by the latter effect, while some influence of the former effect
%is likely also at work in the higher $\alpha=0.1$ case.



The evolution of the disk eccentricity and inclination are well coordinated across different radii of the disk, as is evident from the
similarity of the eccentricity and inclination curves across
different radii in Figures \ref{fig:resolution} and \ref{fig:viscosity}.
Figure \ref{fig:edge_on} shows an edge-on view of the fiducial disk at $t=0$ (left panel) and at $t=13 P_{\rm b}$
(right panel) when
the line of nodes has rotated by about $180^\circ$.
This figure shows that the disk remains fairly flat, but quite lopsided, after substantial nodal and KL phase oscillation changes.

%The KL cycle periods at these two radii are also similar $\tau_{\rm KL}\simeq 8 P_{\rm b}$. 
We consider how inclination is radially communicated across the a disk.
A KL disk undergoes tilt oscillations, along with nodal precession.
A disk that is in good radial communication would be expected
to remain nearly flat, that is, have a small degree of warping.
The disk warping is often measured by $r \partial_r \bm{\ell}$
where $\bm{\ell}$ is a unit vector that lies perpendicular to
the local disk orbital plane. In the application to our simulation results, $\bm{\ell}(r)$ would be determined by the average value of $\bm{\ell}$ for particles in a narrow radial bin centered on radius $r$. 
However, as noted above, there are artificial 
correlations in quantities such as $\bm{\ell}(r)$. 
The reason is that eccentric streamlines carry a single value of $\bm{\ell}$ and extend radially in the disk.
A better measure of warping is $a \partial_a \bm{\ell}$ where $\bm{\ell}(a)$
is determined by the average value of $\bm{\ell}$ for particles whose semi-major axes lie in a narrow bin of size $3.5 \times 10^{-3} a_{\rm b}$ centered about semi-major axis  $a$.   The reason for using $a$ is that (noncrossing) streamlines in a fluid disk 
can be uniquely labeled by their $a$ values.
Therefore, $a \partial_a \bm{\ell}$ measures
warping across streamlines.

As an approximate measure of disk warping, 
we compute
\begin{equation}
\delta \hat{j}= |\bf{\hat{j}_2} - \bf{\hat{j}_1}|,
\label{dj}
\end{equation}
where $\bf{\hat{j}_1}$ and $\bf{\hat{j}_2}$ are 
the average of the disk unit angular momentum vectors in bins of particle semi-major axis centered
on $a_1$ and $a_2$, respectively, that are separated by $a_1$ apart.  We choose these values to be $a_1=0.1 a_{\rm b}$ and $a_2=0.2 a_{\rm b}$. 
For the KL effect on test particles in a circular orbit binary,  $a$ is conserved (Equation (\ref{eq:adot})).
By a time $\sim 10 P_{\rm b}$, the value of $\delta \hat{j}$ for test
particles is of order unity, since the outer particle that begins on a circular orbit
with $r_0=a=0.2 a_b$ undergoes a substantial change in tilt, as seen in Figure \ref{fig:testpart}, as well as substantial nodal precession.

Radial communication in a disk is diffusive 
for $\alpha > H/r$ and wave-like for $\alpha < H/r$.
In the diffusive regime, the degree of warping is estimated by equating
the binary
nodal precessional torque to the horizontal (parallel to the disk) component  
of viscous torque  that involves the so-called $\nu_2$ viscosity \citep{Natarajan98}.
The $\nu_2$ effective viscosity owes its existence to the internal flows driven by pressure gradients and limited by (e.g. turbulent) viscosity, as described by \cite{Papaloizou83}. 
We roughly estimate then that the warping due to nodal precession is
\begin{equation}
\delta \hat{j} \simeq \lambda \sin{(i)} \, \alpha \, \left( \frac{\omega_{\rm p}}{\Omega} \right) \, \left(\frac{r}{H} \right)^2,
\label{warp}
\end{equation}
where $\lambda$ is a coefficient of order unity, $\Omega$ is the
circular frequency at the disk radius $r=0.2 a_{\rm b}$, and $\omega_{\rm p}$ is the
nodal precession frequency. 
We have estimated that $\lambda \sim  0.3$ by solving the 1D viscous warp
evolution equation given by Equation (10) of \cite{King13} for the conditions in the fiducial model
that has  $\omega_{\rm p} \sim 0.04 \Omega_{\rm b}$ and find that
$\delta \hat{j} \sim 0.1$, which is roughly similar in value to
the plotted results  in left panel of Figure \ref{fig:warp_tilt_a} for $\alpha =0.1$ (blue curve).





The lower viscosity $\alpha=0.01$ case plotted in left panel Figure \ref{fig:warp_tilt_a}
lies in the wave-like  regime.
The communication is a consequence of radial pressure forces. 
In this simulation
the radial sound crossing timescale
is substantially shorter than the KL oscillation or nodal precessional timescale. 
We have that $r/(c_{\rm s} P_{\rm KL}) \sim 0.1.$ 
As described in Paper 1,  because this ratio is small
the disk undergoes global rigid tilt oscillations \citep{Larwood97, Lubow01}. In this case, the level of warping 
(plotted as the red curve  on the left panel of Figure \ref{fig:warp_tilt_a})
is similar to or weaker than in the viscous regime (blue curve).


%As noted above, the disk warp is subject to the effects of nodal precession
%and KL tilt oscillations. 
%In comparing left panel of Figure \ref{fig:warp_tilt} 
%with Figure \ref{fig:viscosity}, we see 
%that the degree of warping varies somewhat with phase
%of the KL cycle. At times of $t\sim 18 P_{tm b}$ and  $t\sim 25 P_{tm b}$ When the disk inclination is minimum and eccentricity is maximum,
%the warping is tends to be smaller.   
 
 The warping $\delta \hat{j}$ depends on the variation of both the disk tilt 
  and direction of the line of nodes. It is affected by both nodal precession
 and KL disk tilt oscillations.
 We consider here the variations of only the disk tilt.
The KL  inclination oscillations can lead to variations in disk tilt across streamlines in
the disk.  
These variations can be expected because the KL oscillation timescale
varies with particle orbit period $P$ and therefore with particle orbit semi-major axis $a$. 
According to Equation (\ref{eq:klperiod}), the KL oscillation period should differ by about a factor of $\sim 3$  in a change from
$a=0.1a_{\rm b}$ to $a=0.2 a_{\rm b}$. Radial communication via viscous and pressure effects
as described above play a role in maintaining a constant tilt.
As seen in the right panel of Figure \ref{fig:warp_tilt_a}, the change in
tilt between $a=0.1 a_{\rm b}$ and $a=0.2 a_{\rm b}$ is small.
The difference in tilt angles between them is typically less than $3^\circ$. This
change is similar to the disk aspect ratio $H/r = 0.035 \simeq 2.0^\circ$, while the disk oscillates over a range of tilt angles
$\simeq 20^\circ$.
In comparing the right panel of Figure \ref{fig:warp_tilt_a} 
with Figure \ref{fig:viscosity}, we see 
that tilt difference is phased with 
the KL cycle.  When the disk inclination is minimum and eccentricity is maximum,
the outer inclination $i_2$ (that has $a = 0.2 a_{\rm b}$) is smaller than the inner inclination $i_1$ (that has $a = 0.1 a_{\rm b}$) .  





The effect of varying the disk aspect ratio is shown in Figure \ref{fig:hoR}. Most noticeably, as the disk aspect ratio decreases (from the blue curve to the red curve), the first disk KL cycle gets delayed to a later time. Because the kinematic viscosity $\nu$ varies with $(H/r)^2$, lower values of $H/r$ lead to a slower disk expansion and therefore a weaker perturbation from the companion object.
Lower $\nu$ values also lead to slower damping in both eccentricity and inclination. %Notice that even at $H/r=0.02$ (low disk sound speed, thus less efficient communication), the disk can still maintain a high level of rigidity and does not exhibit significant warping. The disk is in good radial communication even in this case, since the sound crossing time is only a few times $P_{\rm b}$ that is less that
%the KL oscillation time of more than $11 P_{\rm b}$. %We do not show the results for the warmest disk that has $H/r = 0.065$ beyond
%$T=25 P_{\rm b}$ because the number of disk particles is very low, resulting in large fluctuations.
%DO NOT SHOW YELLOW LINE PAST T=25 Pb.
The bottom panel of the figure shows the expected KL oscillation periods based on
Equation (\ref{eq:klperiod2}). As we saw in the case of Figure \ref{fig:viscosity},
the accuracy of this equation is only to about a factor of 2. %However,
%the results may show better relative agreement for different values of $H/r$.
Based on the bottom panel, the ratio of KL periods for $H/r=0.065$ to $H/r=0.02$
is about 1.3, while the ratio in the SPH simulations (upper panel) is about 1.5.


We present in Figure \ref{fig:gamma} three runs using different power-law indices for the initial disk surface density profile. The differences in the results for these three cases is found to be quite small.  The bottom left panel of this figure
shows that the expected KL periods based on Equation
(\ref{eq:klperiod2}) are about the same after a time of 
about $10 P_{\rm b}$, as is in agreement with the SPH simulations in the upper panels. The bottom right panel shows that after initial adjustment, the
density distributions are very similar.


In Figure \ref{fig:rin}, we show the effects of changing the disk inner radius. 
The first disk KL cycle does not seem to be affected by the location of the disk inner boundary. However, the cycle period does get slightly longer as the disk inner radius decreases. 
The bottom panel of this figure
shows that the expected KL periods based on Equation
(\ref{eq:klperiod2}) are in qualitative agreement with this trend. 
This effect is likely be the result of having more material located further from the companion, closer to the central object, resulting
in weaker global perturbations. 

Figure \ref{fig:rout} shows the effects of changing the disk outer radius. 
%With smaller initial disk outer radius, the outer parts of the disk take longer to expand to the point when the perturbation from the secondary can efficiently excite KL cycles. 
With smaller initial disk outer radius, the disk initially responds more slowly since the perturbations from the companion object are weaker. 
Thus, the effect in this case is simply a time delay in the cycles.
The bottom panel of this figure
shows that the expected KL periods based on Equation
(\ref{eq:klperiod2}) are in agreement with this description. 
The case with the smaller initial radius has an initial KL period
that is much longer than for disk of large initial radius. By a time of $\sim 15 P_{\rm b}$,
both cases have the nearly same KL periods, as seen in the SPH simulations results
of the upper panels. 

%CONSISTENT WITH VISCOUS EXPANSION?? NEED TO SEE SIGMA(R).




Figure \ref{fig:inc} shows the effects of changing the initial disk inclination. As $i_0$ varies from $60^{\circ}$ to $45^{\circ}$, the KL cycles occur at a progressively later time with a longer cycle period and smaller oscillation magnitude. This trend seems to be in line with the test particle results
discussed in the previous section, where the KL oscillation period increased with decreasing inclination. Equation (\ref{eq:maxe}) predicts eccentricity maxima of $0.76$, $0.67$, and $0.56$ for initial tilts of $60^\circ$, $55^{\circ}$, and $50^{\circ}$, respectively, while Figure \ref{fig:inc} exhibits maximum $e$ values during the first KL cycle of $0.7$, $0.6$, and $0.4$. In the test particle case, for an initial tilt of $i_0=45^{\circ}$, the
eccentricity can reach values of $0.4$, based on Equation (\ref{eq:maxe}), or 0.3, based on simulations (top-left panel of Figure \ref{fig:testpart}). 
However, for this initial inclination angle, the KL mechanism barely operates in the disk whereas it is still relatively strong in the test particle runs (see the left column of Figure \ref{fig:testpart}). Thus, the critical initial tilt angle for the onset of KL oscillations in a hydrodynamical disk is somewhat higher than for a test particle ($39^\circ$). On the other hand, in the first disk KL cycle, the inclination can drop  below $39^{\circ}$, which does not occur in the particle case.  

Figure \ref{fig:mratio} demonstrates the effects of varying the binary mass ratio. For $M_{\rm p}/M_{\rm c} \leq1$,  a smaller mass ratio leads to a later occurrence of the KL cycle and also a longer cycle period (from the green curve to the red curve), which again is in line with the test particle results. 
The bottom panel of this figure
shows that the expected KL periods based on Equation
(\ref{eq:klperiod2})  at a time of $25 P_{\rm b}$ are about 
$30 P_{\rm b}$, $20 P_{\rm b}$, $15 P_{\rm b}$, and $7 P_{\rm b}$ for mass ratios of $1/4$, $1/2$, $1$, and 2, respectively, whereas the SPH simulation results show 
$15 P_{\rm b}$, $11 P_{\rm b}$, $8 P_{\rm b}$, and $11 P_{\rm b}$. 
For both approaches for mass ratios less than or equal to unity, the KL period decreases with mass ratio, as expected. The period for the corresponding
cases are again different by about a factor or 2. 
In the case where
the perturbing companion is more massive with a mass ratio of 2, the period in the SPH simulation
does not follow this trend and instead increases.

 The effects of varying the binary eccentricity are shown in Figure \ref{fig:eb}.
 The initial disk outer radius is $r_{\rm out} = 0.2 a_{\rm b}$ instead of the fiducial value of $0.25 a_{\rm b}$. 
 %ARE THE INITIAL DISK RADII THE SAME AS BEFORE?
 This smaller initial  disk outer radius was used because the increased binary eccentricity
 causes a stronger response in the disk and a smaller disk truncation radius, as is known to occur for the coplanar case
 \citep{Artymowicz94}.
 Other parameters are set to their fiducial values.
 We do not plot the results for $e_{\rm b} = 0.5$ at $r=0.2 a_{\rm b}$ and $r=0.3 a_{\rm b}$
 because the disk is truncated within this radius due to the smaller periastron
 distance. %The effect of binary eccentricity is to cause the KL oscillation period
 %to be about 10\% shorter for the $e_{\rm b}=0.5$ case relative to the $e_{\rm b}=0$ case,
 %while the KL period for the $e_{\rm b}=0.2$ was about the same as the $e_{\rm b}=0$ case.
 The KL oscillation periods are about the same in the three cases, within about 10\% of each other.
 The maximum eccentricities are also about the same in the three cases, as we found for test particles
 in Section \ref{sec:part}. In the $e_{\rm b}=0.5$ case, the inclination decays to a value close to $i_{\rm cr}=39^\circ$,
 while the lower binary eccentricity cases decay to a smaller tilt angle.
% Equation \ref{eq:klperiod} predicts KL period of $16P_{\rm b}$, $15P_{\rm b}$ and $10P_{\rm b}$ for $e_{\rm b}=0$, $0.2$ and $0.5$, respectively, whereas Figure \ref{fig:eb} gives $9P_{\rm b}$, $8P_{\rm b}$, and $8P_{\rm b}$.  FOR WHAT RADIUS? Equation \ref{eq:maxe} predicts maximum eccentricity of $0.56$, whereas Figure \ref{fig:eb} shows $e_{\rm max} \sim 0.5$ for $e_{\rm b}=0.2$ and $0.5$, $e_{\rm max} \sim 0.4$ for $e_{\rm b}=0$. 
 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%      Summary      %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

\section{ Discussion and Summary \label{sec:sum}}





In this paper, we investigated the conditions under which KL oscillations
can operate on a fluid disk, as was first described in \cite{Martin14b}.
Such a disk undergoes a complicated evolution involving tilt and eccentricity
oscillations, as well as the usual nodal and apsidal precession.
We carried out this analysis
by means of SPH simulations that include the effects of disk pressure and viscosity. We considered a range of disk viscosities (Figure \ref{fig:viscosity}),
sound speeds (Figure \ref{fig:hoR}), initial density profiles (Figures \ref{fig:gamma}, \ref{fig:rin}, \ref{fig:rout}), initial inclinations 
(Figure \ref{fig:inc}), binary mass ratios (Figure \ref{fig:mratio}),
and binary eccentricities (Figure \ref{fig:eb}).
We found that the KL effect generally operates over a range of
of disk and binary properties. 

The general picture we have is that an initially circular circumstellar disk that is sufficiently inclined ($\sim 45^{\circ}$ -- $135^{\circ}$) with respect to the orbit plane of a binary undergoes damped KL eccentricity and tilt oscillations. During these oscillations, the disk can develop substantial eccentricities, $e \ga 0.5$. Once the KL oscillations have mostly damped, the disk achieves a tilt in the prograde disk case  (initial tilt less than $90^{\circ}$)   $i \simeq 39^{\circ}$ or somewhat smaller and an eccentricity of up to a few tenths. Over longer timescales, if the disk survives that long, we expect that the disk eccentricity will eventually damp \citep[e.g.,][]{Ogilvie01a} and the disk will become coplanar with the binary orbit plane \citep[e.g.,][]{King13} due to viscous effects.

The properties of KL disk oscillations are somewhat similar to
the properties of KL test particle oscillations (Figure \ref{fig:testpart}).
In both cases, inclination and eccentricity are interchanged.
The peak eccentricity achieved by a disk in the first KL oscillation
is generally close to the eccentricity achieved by test particles for the same initial inclination.
In addition, for smaller mass perturbers, the KL oscillations are slower for
both test particles and disks
(Figure \ref{fig:mratio}).

However, there are some significant differences between a fluid disk and test particles.
The KL oscillation frequency of test particles varies with particle orbit semi-major axis
as $a^{3/2}$.
On the other hand, a disk undergoes a large scale global oscillation with a small
level of warping for typical parameters (Figures \ref{fig:edge_on} and \ref{fig:warp_tilt_a}). 
For a circular
orbit binary, the inclination of prograde test particle orbits relative to the binary orbit plane  
must be greater than the critical angle $i_{\rm cr}=\arccos{(\sqrt{3/5})} \simeq 39^\circ$
for KL oscillations to occur. In the case of a disk, we find that the minimum
inclination angle appears to be somewhat larger ($\ga 45^{\circ}$)  (see Figure  
\ref{fig:inc}). 
For a circular
orbit binary, the inclination of prograde test particle orbits oscillates
between the initial value $i_0 > i_{\rm cr}$ and $i_{\rm cr}$.
For the case of a disk, however, the oscillation can dip below $i_{\rm cr}$
(e.g., Figure \ref{fig:viscosity}). Therefore, even though a disk
requires a larger tilt to operate than a test particle orbit,
the disk oscillations can reach lower tilt angles.

Another major difference between test particles and disks is that
disks can dissipate energy. As a result, the KL oscillations damp,
but leave a residual disk eccentricity that we find to be typically a few tenths. 
We find that the damping rate
depends on the level of disk viscosity (Figure \ref{fig:viscosity}).
The disk eccentricities may also damp through a decay process involving
a parametric instability \citep{Papaloizou05a, Papaloizou05b, Barker14, Ogilvie14}.
Our simulations may not be capable of resolving this instability that
occurs on scales smaller than the disk thickness.
It is also possible that the initial KL oscillation damping involves shocks,
since the streamlines may attempt to cross at high eccentricity.
We have found indirect evidence for such effects in some simulations
that we carried out with a low value of viscosity parameter $\beta_{\rm AV}$.
As mentioned in Section \ref{sec:kld},
this parameter is used in the SPH code to prevent interparticle penetration that
artificially occurs in the code where shocks are present. For small values of $\beta_{\rm AV}$,
the SPH simulations displayed unphysical behavior, symptomatic
of interparticle penetration.
Over longer timescales in the post-KL oscillation phase, the disk is expected to eventually align
with the binary orbit plane through viscous effects \cite[e.g.,][]{King13}.


For an eccentric orbit binary, the KL test particle oscillations can
reach to higher inclination angles than the initial tilt. Inclination angles can
reach $90^{\circ}$ and even higher, resulting flipping from prograde to retrograde
orbits \citep{Lithwick11, Naoz11, Naoz13a, Naoz13b, Liu15}. In reaching $90^{\circ}$, the orbit becomes radial.
%As in the case of test particles, binary eccentricity (at fixed binary semi-major axis) shortens the KL oscillation period
%(Figure \ref{fig:eb}). 
KL disk oscillations cannot achieve a similar increased
level of inclination (certainly not to eccentricity close to unity) because of dissipation. 

In this paper, we have omitted the effects of disk self-gravity.
Such effects may suppress KL oscillations for sufficiently massive disks
\citep{Batygin12, Martin14b}. We intend to describe the influence of self-gravity
in a future paper.

  
Confirmation of the KL mechanism acting on a disk would come from an observation of an eccentric and misaligned disk in a binary system. The KL oscillations are applicable to binary systems on all scales. For example, the oscillations could play a role in planet formation and evolution around one component of a binary star.   There is observational evidence that some exoplanets may be undergoing KL cycles in wide binaries \citep[e.g.][]{Wu03,Takeda05}. In future work, we plan to investigate the evolution of a planet--disk system that orbits around one component of a binary system in order to understand how these planets can form. 

Be/X-ray binaries may also be influenced by the KL oscillations of a misaligned disk \citep{Martin14a,Martin14b}. Be stars may form disks from material that is ejected from the equator of the rapidly rotating star \citep{Lee91,Hanuschik96,Carciofi11}. Be/X-ray binaries contain a Be star with a rotation axis that is misaligned to that of the orbit of a companion neutron star \citep[see Table 1 in][for some observations]{Martin11}. We found in our simulations that roughly 10\% of the initial disk mass
is captured by the companion object for the fiducial model. As a disk undergoes KL oscillations to maximum eccentricity,
its apastron radius grows  which allows this capture to occur. 
This mass exchange may explain the so-called giant Type II
outbursts as being due the accretion of disk mass by the companion neutron star \citep{Martin14a}.



%The two types of X-ray outburst observed can be understood with a misaligned decretion disk model. The type~I outbursts occur periodically during the periastron passage of the eccentric orbit neutron star. The giant type~II X-ray outbursts can be explained by eccentricity growth of the disk by the KL effect that results in a large mass transfer to the neutron star \citep{Martin14a}.   Because the radial extent of a Be star disk is small (only about 7 times the stellar radius), the eccentricity growth of the disk can cause the complete depletion of the disk as the majority is accreted onto the Be star, while some is accreted onto the neutron star. Thus, in these systems we don't expect the disc to display more than one KL oscillation. However, the disk will be replenished with more material from the Be star, and we expect this single oscillation cycle to repeat. 

Finally, the KL oscillations of a disk may also be relevant to SMBH binaries. It is possible that the circumbinary disk around a SMBH binary will be misaligned \citep[e.g.][]{Nixon11a,Nixon11b}. Thus, as gas accretes inwards from the circumbinary disk, misaligned disks will form around each SMBH. These disks will be unstable to KL oscillations if the misalignment angle is sufficiently large. This effect could have implications for the subsequent disk evolution and star formation in these systems.
 

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Acknowledgements %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

W.F. and S.H.L. acknowledge support from NASA grant NNX11AK61G.  Computing resource supporting this work was provided by the institutional computing program at Los Alamos National Laboratory. We thank Daniel Price for providing the \texttt{PHANTOM} code for SPH simulations and \texttt{SPLASH} code \citep{Price07} for data analysis and rendering of figures. 



%%%%%%%%%%%%%%%%%%%%%%%%%
%%%     References    %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{thebibliography}{}

\bibitem[Artymowicz \& Lubow (1994)]{Artymowicz94} Artymowicz, P., \& Lubow,
S. H. 1994, \apj, 421, 651 

\bibitem[Albrecht et al. (2012)]{Albrecht12}
Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2012, \apj, 757, 18


\bibitem[Barker \& Ogilvie (2014)]{Barker14}
Barker, A. J., \& Ogilvie, G. I. 2014, \mnras, 445, 2637

\bibitem[Bate et al. (2000)]{Bate00}
Bate, M. R., Bonnell, I. A., Clarke, C. J., et al. 2000, \mnras, 317, 773

\bibitem[Bate et al. (2010)]{Bate10}
Bate, M. R., Lodato, G., \& Pringle, J. E. 2010, \mnras, 401, 1505

\bibitem[Batygin (2012)]{Batygin12}
Batygin, K. 2012, \nat, 491, 418

\bibitem[Blaes et al. (2002)]{Blaes02}
Blaes, O., Lee, M. H., \& Socrates, A. 2002, \apj, 578, 775

\bibitem[Bonnell \& Bate(1994)]{Bonnell94} 
Bonnell, I.~A., \& Bate, M.~R.\ 1994, \mnras, 269, L45  

\bibitem[Carciofi (2011)]{Carciofi11}
Carciofi, A. C. 2011, in IAU Symp. 272, Active OB Stars: Structure, Evolution,
Mass Loss, and Critical Limits, ed. C. Neiner, G. Wade, G. Meynet, \& G.
Peters (Cambridge: Cambridge Univ. Press), 325

\bibitem[Chen et al. (2011)]{Chen11}
Chen, X., Sesana, A., Madau, P., \& Liu, F. K. 2011, \apj, 729, 13

\bibitem[Dawson \& Chiang (2014)]{Dawson14}
Dawson, R. I. \& Chiang, E. 2014, Science, 346, 212

\bibitem[Dong et al. (2014)]{Dong14}
Dong, S., Katz, B., \& Socrates, A. 2014, ApJL, 781, L5

\bibitem[Fabrycky \& Tremaine (2007)]{Fabrycky07}
Fabrycky, D., \& Tremaine, S. 2007, \apj, 669, 1298

\bibitem[Ford et al. (2000)]{Ford00}
Ford, E. B., Kozinsky, B., \& Rasio, F. A. 2000, \apj, 535, 385

\bibitem[Hale (1994)]{Hale94}
Hale, A. 1994, \aj, 107, 306

\bibitem[Hanuschik (1996)]{Hanuschik96}
Hanuschik, R. W. 1996, A\&A, 308, 170

\bibitem[Harrington (1968)]{Harrington68}
Harrington, R. S. 1968, \aj, 73, 190

\bibitem[Holman et al. (1997)]{Holman97}
Holman, M. J., Touma, J., \& Tremaine, S. 1997, \nat, 386, 254

\bibitem[Innanen et al. (1997)]{Innanen97}
Innanen, K. A., Zheng, J. Q., Mikkola, S., \& Valtonen, M. J. 1997, \aj, 113, 1915

\bibitem[Jensen \& Akeson (2014)]{Jensen14}
Jensen, E. L. N. \& Akeson, R. 2014, \nat, 511, 567

%\bibitem[Jensen et al. (2004)]{Jensen04}
%Jensen, E. L. N., Mathieu, R. D., Donar, A. X., \& Dullighan, A. 2004, \apj, 600, 789

\bibitem[King 
\& Pringle(2006)]{King06} King, A.~R., \& Pringle, J.~E.\ 2006, \mnras, 373, L90 

\bibitem[King et al. (2013)]{King13}
King, A. R., Livio, M., Lubow, S. H., \& Pringle, J. E. 2013, \mnras, 431, 2655

\bibitem[Kiseleva et al. (1998)]{Kiseleva98}
Kiseleva, L. G., Eggleton, P. P., \& Mikkola, S. 1998, \mnras, 300, 292

\bibitem[Kozai (1962)]{Kozai62}
Kozai, Y. 1962, \aj, 67, 591

\bibitem[Kushnir et al. (2013)]{Kushnir13}
Kushnir, D., Katz, B., Dong, S., Livne, E., \& Fern\'{a}ndez, R. 2013, ApJL, 778, L37

%\bibitem[Larwood et al. (1996)]{Larwood96}
%Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., \& Terquem, C. 1996, \mnras, 282, 597

\bibitem[Larwood \& Papaloizou (1997)]{Larwood97}
Larwood, J. D., \& Papaloizou, J. C. B. 1997, \mnras, 285, 288

\bibitem[Lee, Osaki \& Saio (1991)]{Lee91}
Lee, U., Osaki, Y., \& Saio, H. 1991, \mnras, 250, 432

\bibitem[Li et al. (2014)]{Li14}
Li, G., Naoz, S., Kocsis, B., \& Loeb, A. 2014, \apj, 785, 116

\bibitem[Lidov (1962)]{Lidov62}
Lidov, M. L. 1962, P\&SS, 9, 719

\bibitem[Lithwick \& Naoz (2011)]{Lithwick11}
Lithwick, Y., \& Naoz, S. 2011, \apj, 742, 94

\bibitem[Liu et al. (2015)]{Liu15}
Liu, B., Mu\~{n}oz, D., \& Lai, D. 2015, \mnras, 447, 747

\bibitem[Lodato \& Pringle (2007)]{Lodato07}
Lodato, G., \& Pringle, J. E. 2007, MNRAS, 381, 1287

\bibitem[Lodato \& Price (2010)]{Lodato10}
Lodato, G., \& Price, D. J. 2010, \mnras, 405, 1212

\bibitem[Lubow, Martin \& Nixon (2015)]{Lubow15}
Lubow, S. H., Martin, R. G., \& Nixon, C. J. 2015, ApJ, 800, 96

\bibitem[Lubow \& Ogilvie (2000)]{Lubow00}
Lubow, S. H., \& Ogilvie, G. I. 2000, \apj, 538, 326

\bibitem[Lubow \& Ogilvie (2001)]{Lubow01}
Lubow, S. H., \& Ogilvie, G. I. 2001, \apj, 560, 997

\bibitem[Lynden-Bell \& Pringle (1974)]{Lynden-Bell74}
Lynden-Bell, D. \& Pringle, J. E., MNRAS, 168, 603

\bibitem[Martin et al. (2009)]{Martin09}
Martin R. G., Tout, C. A., \& Pringle J. E. 2009, MNRAS, 397, 1563

\bibitem[Martin et al. (2011)]{Martin11}
Martin R. G.,  Pringle J. E., Tout C. A., \& Lubow S. H. 2011, MNRAS, 416, 2827

\bibitem[Martin et al. (2014a)]{Martin14a}
Martin, R. G., Nixon, C., Armitage, P. J., Lubow, S. H., \& Price, D. J. 2014a, ApJL, 790, L34

\bibitem[Martin et al. (2014b)]{Martin14b}
Martin, R. G., Nixon, C., Lubow, S. H., et al. 2014b, \apj, 792, L33

\bibitem[Mazeh et al. (1997)]{Mazeh97}
Mazeh, T., Krymolowski, Y., \& Rosenfeld, G. 1997, ApJL, 477, L103

\bibitem[Mazeh \& Shaham (1979)]{Mazeh79}
Mazeh, T., \& Shaham, J. 1979, A\&A, 77, 145

\bibitem[Miller \& Hamilton (2002)]{Miller02}
Miller, M. C., \& Hamilton, D. P. 2002, \apj, 576, 894

\bibitem[Naoz et al. (2011)]{Naoz11}
Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., \& Teyssandier, J. 2011, \nat, 473, 187

\bibitem[Naoz et al. (2013a)]{Naoz13a}
Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., \& Teyssandier, J. 2013a, \mnras, 431, 2155

\bibitem[Naoz et al. (2013b)]{Naoz13b}
Naoz, S., Kocsis, B., Leob A., \& Yunes, N. 2013b, \apj,773, 187

\bibitem[Natarajan \& Pringle (1998)]{Natarajan98}
Natarajan, P., \& Pringle, J. E. 1998, ApJL, 506, L97

\bibitem[Nesvorn\'{y} et al. (2003)]{Nesvorny03}
Nesvorn\'{y}, D., Alvarellos, J. L. A., Dones, L., \& Levison, H. F. 2003, \aj, 126, 398

\bibitem[Nixon et al. (2011a)]{Nixon11a}
Nixon, C. J., Cossins, P. J., King, A. R.,  Pringle, J. E. 2011a, MNRAS,
412, 1591

\bibitem[Nixon et al. (2011b)]{Nixon11b}
Nixon, C. J., King, A. R.,  Pringle, J. E. 2011b, MNRAS, 417, L66

\bibitem[Nixon (2012)]{Nixon12a}
Nixon, C. J. 2012, MNRAS, 423, 2597

\bibitem[Nixon \& King (2012)]{Nixon12b}
Nixon, C. J., \& King, A. R. 2012, MNRAS, 421, 1201

\bibitem[Nixon et al.(2013)]{Nixon13}
Nixon, C., King, A., \& Price, D. 2013, MNRAS, 434, 1946

\bibitem[Offner et al. (2010)]{Offner10}
Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz, M. R., \& Klein, R. I. 2010, \apj, 725, 1485

%\bibitem[Ogilvie (1999)]{Ogilvie99}
%Ogilvie, G. I. 1999, \mnras, 304, 557

\bibitem[Ogilvie (2001)]{Ogilvie01a}
Ogilvie, G. I. 2001, \mnras, 325, 231

\bibitem[Ogilvie \& Barker (2014)]{Ogilvie14}
Ogilvie, G. I., \& Barker, A. J. 2014, \mnras, 445, 2621

\bibitem[Ogilvie \& Dubus (2001)]{Ogilvie01b}
Ogilvie, G. I., \& Dubus, G. 2001, \mnras, 320, 485

\bibitem[Paczy\'{n}ski (1977)]{Paczynski77}
Paczy\'{n}ski, B. 1977, ApJ, 216, 822

\bibitem[Papaloizou (2005a)]{Papaloizou05a}
Papaloizou, J. C. B. 2005a, A\&A, 432, 743

\bibitem[Papaloizou (2005b)]{Papaloizou05b}
Papaloizou, J. C. B. 2005b, A\&A, 432, 757

\bibitem[Papaloizou \& Pringle (1983)]{Papaloizou83}
Papaloizou, J. C. B., \& Pringle, J. E. 1983, \mnras, 202, 1181

\bibitem[Papaloizou \& Terquem (1995)]{Papaloizou95}
Papaloizou, J. C. B., \& Terquem, C. 1995, \mnras, 274, 987

\bibitem[Petrovich (2015)]{Petrovich15}
Petrovich, C. 2015, \apj, 799, 27

\bibitem[Porter \& Rivinus (2003)]{Porter03}
Porter, J M., \& Rivinius, T. 2003, PASP, 115, 1153

\bibitem[Price (2007)]{Price07}
Price, D. J. 2007, \pasa, 24, 159

\bibitem[Price \& Federrath (2010)]{Price10}
Price, D. J., \& Federrath, C. 2010, \mnras, 406, 1659

\bibitem[Price (2012)]{Price12}
Price, D. J. 2012, JCoPh, 231, 759

%\bibitem[Pringle (1981)]{Pringle81}
%Pringle, J. E. 1981, ARA\&A, 19, 137

%\bibitem[Pringle (1992)]{Pringle92}
%Pringle, J. E. 1992, MNRAS, 258, 811

\bibitem[Pringle (1996)]{Pringle96}
Pringle, J. E. 1996, \mnras, 281, 357

\bibitem[Rice (2015)]{Rice15}
Rice, K. 2015, \mnras, 448, 1729

\bibitem[Storch et al. (2014)]{Storch14}
Storch, N., Kassandra, R. A., \& Lai, D. 2014, Science, 345, 1317

\bibitem[Storch et al. (2015)]{Storch15}
Storch, N., \& Lai, D. 2015, \mnras, 448, 1821

\bibitem[Shakura \& Sunyaev (1973)]{SS73}
Shakura, N. I., \& Sunyaev, R. A., 1973, A\&A, 24, 337

\bibitem[Stapelfeldt et al. (1998)]{Stapelfeldt98}
Stapelfeldt, K. R., Krist, J. E., M\'{e}nard, F., et al. 1998, ApJL, 502, L65

\bibitem[Takeda \& Rasio (2005)]{Takeda05}
Takeda, G., \& Rasio, F. A. 2005, \apj, 627, 1001

\bibitem[Teyssandier et al. (2013)]{Teyssandier13}
Teyssandier, J., Naoz, S., Lizarraga, I., \& Rasio, F. A. 2013, \apj, 779, 166

\bibitem[Wijers \& Pringle (1999)]{Wijers99}
Wijers, R. A. M. J., \& Pringle, J. E. 1999, \mnras, 308, 207

\bibitem[Williams et al. (2014)]{Williams14}
Williams, J. P., Mann, R. K., Di Francesco, J., et al. 2014, \apj, 796, 120

\bibitem[Winn \& Fabrycky (2014)]{Winn14}
Winn, J. N., \& Fabrycky, D. C. 2014, arXiv:1410.4199

\bibitem[Wu \& Murray (2003)]{Wu03}
Wu, Y., \& Murray, N. 2003, \apj, 589, 605

\bibitem[Xiang-Gruess \& Papaloizou (2013)]{Xiang13}
Xiang-Gruess, M., \& Papaloizou, J. C. B. 2013, \mnras, 431, 1320


\end{thebibliography}


\end{document}

1504.00393.extracted_bib
[Artymowicz \& Lubow (1994)]{Artymowicz94 Artymowicz, P., \& Lubow, S. H. 1994, \apj, 421, 651
1504.00393.extracted_bib
[Albrecht et al. (2012)]{Albrecht12 Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2012, \apj, 757, 18
1504.00393.extracted_bib
[Barker \& Ogilvie (2014)]{Barker14 Barker, A. J., \& Ogilvie, G. I. 2014, \mnras, 445, 2637
1504.00393.extracted_bib
[Bate et al. (2000)]{Bate00 Bate, M. R., Bonnell, I. A., Clarke, C. J., et al. 2000, \mnras, 317, 773
1504.00393.extracted_bib
[Bate et al. (2010)]{Bate10 Bate, M. R., Lodato, G., \& Pringle, J. E. 2010, \mnras, 401, 1505
1504.00393.extracted_bib
[Batygin (2012)]{Batygin12 Batygin, K. 2012, \nat, 491, 418
1504.00393.extracted_bib
[Blaes et al. (2002)]{Blaes02 Blaes, O., Lee, M. H., \& Socrates, A. 2002, \apj, 578, 775
1504.00393.extracted_bib
[Bonnell \& Bate(1994)]{Bonnell94 Bonnell, I.~A., \& Bate, M.~R.\ 1994, \mnras, 269, L45
1504.00393.extracted_bib
[Carciofi (2011)]{Carciofi11 Carciofi, A. C. 2011, in IAU Symp. 272, Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits, ed. C. Neiner, G. Wade, G. Meynet, \& G. Peters (Cambridge: Cambridge Univ. Press), 325
1504.00393.extracted_bib
[Chen et al. (2011)]{Chen11 Chen, X., Sesana, A., Madau, P., \& Liu, F. K. 2011, \apj, 729, 13
1504.00393.extracted_bib
[Dawson \& Chiang (2014)]{Dawson14 Dawson, R. I. \& Chiang, E. 2014, Science, 346, 212
1504.00393.extracted_bib
[Dong et al. (2014)]{Dong14 Dong, S., Katz, B., \& Socrates, A. 2014, ApJL, 781, L5
1504.00393.extracted_bib
[Fabrycky \& Tremaine (2007)]{Fabrycky07 Fabrycky, D., \& Tremaine, S. 2007, \apj, 669, 1298
1504.00393.extracted_bib
[Ford et al. (2000)]{Ford00 Ford, E. B., Kozinsky, B., \& Rasio, F. A. 2000, \apj, 535, 385
1504.00393.extracted_bib
[Hale (1994)]{Hale94 Hale, A. 1994, \aj, 107, 306
1504.00393.extracted_bib
[Hanuschik (1996)]{Hanuschik96 Hanuschik, R. W. 1996, A\&A, 308, 170
1504.00393.extracted_bib
[Harrington (1968)]{Harrington68 Harrington, R. S. 1968, \aj, 73, 190
1504.00393.extracted_bib
[Holman et al. (1997)]{Holman97 Holman, M. J., Touma, J., \& Tremaine, S. 1997, \nat, 386, 254
1504.00393.extracted_bib
[Innanen et al. (1997)]{Innanen97 Innanen, K. A., Zheng, J. Q., Mikkola, S., \& Valtonen, M. J. 1997, \aj, 113, 1915
1504.00393.extracted_bib
[Jensen \& Akeson (2014)]{Jensen14 Jensen, E. L. N. \& Akeson, R. 2014, \nat, 511, 567 %
1504.00393.extracted_bib
[Jensen et al. (2004)]{Jensen04 %Jensen, E. L. N., Mathieu, R. D., Donar, A. X., \& Dullighan, A. 2004, \apj, 600, 789
1504.00393.extracted_bib
[King \& Pringle(2006)]{King06 King, A.~R., \& Pringle, J.~E.\ 2006, \mnras, 373, L90
1504.00393.extracted_bib
[King et al. (2013)]{King13 King, A. R., Livio, M., Lubow, S. H., \& Pringle, J. E. 2013, \mnras, 431, 2655
1504.00393.extracted_bib
[Kiseleva et al. (1998)]{Kiseleva98 Kiseleva, L. G., Eggleton, P. P., \& Mikkola, S. 1998, \mnras, 300, 292
1504.00393.extracted_bib
[Kozai (1962)]{Kozai62 Kozai, Y. 1962, \aj, 67, 591
1504.00393.extracted_bib
[Kushnir et al. (2013)]{Kushnir13 Kushnir, D., Katz, B., Dong, S., Livne, E., \& Fern\'{andez, R. 2013, ApJL, 778, L37 %
1504.00393.extracted_bib
[Larwood et al. (1996)]{Larwood96 %Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., \& Terquem, C. 1996, \mnras, 282, 597
1504.00393.extracted_bib
[Larwood \& Papaloizou (1997)]{Larwood97 Larwood, J. D., \& Papaloizou, J. C. B. 1997, \mnras, 285, 288
1504.00393.extracted_bib
[Lee, Osaki \& Saio (1991)]{Lee91 Lee, U., Osaki, Y., \& Saio, H. 1991, \mnras, 250, 432
1504.00393.extracted_bib
[Li et al. (2014)]{Li14 Li, G., Naoz, S., Kocsis, B., \& Loeb, A. 2014, \apj, 785, 116
1504.00393.extracted_bib
[Lidov (1962)]{Lidov62 Lidov, M. L. 1962, P\&SS, 9, 719
1504.00393.extracted_bib
[Lithwick \& Naoz (2011)]{Lithwick11 Lithwick, Y., \& Naoz, S. 2011, \apj, 742, 94
1504.00393.extracted_bib
[Liu et al. (2015)]{Liu15 Liu, B., Mu\~{noz, D., \& Lai, D. 2015, \mnras, 447, 747
1504.00393.extracted_bib
[Lodato \& Pringle (2007)]{Lodato07 Lodato, G., \& Pringle, J. E. 2007, MNRAS, 381, 1287
1504.00393.extracted_bib
[Lodato \& Price (2010)]{Lodato10 Lodato, G., \& Price, D. J. 2010, \mnras, 405, 1212
1504.00393.extracted_bib
[Lubow, Martin \& Nixon (2015)]{Lubow15 Lubow, S. H., Martin, R. G., \& Nixon, C. J. 2015, ApJ, 800, 96
1504.00393.extracted_bib
[Lubow \& Ogilvie (2000)]{Lubow00 Lubow, S. H., \& Ogilvie, G. I. 2000, \apj, 538, 326
1504.00393.extracted_bib
[Lubow \& Ogilvie (2001)]{Lubow01 Lubow, S. H., \& Ogilvie, G. I. 2001, \apj, 560, 997
1504.00393.extracted_bib
[Lynden-Bell \& Pringle (1974)]{Lynden-Bell74 Lynden-Bell, D. \& Pringle, J. E., MNRAS, 168, 603
1504.00393.extracted_bib
[Martin et al. (2009)]{Martin09 Martin R. G., Tout, C. A., \& Pringle J. E. 2009, MNRAS, 397, 1563
1504.00393.extracted_bib
[Martin et al. (2011)]{Martin11 Martin R. G., Pringle J. E., Tout C. A., \& Lubow S. H. 2011, MNRAS, 416, 2827
1504.00393.extracted_bib
[Martin et al. (2014a)]{Martin14a Martin, R. G., Nixon, C., Armitage, P. J., Lubow, S. H., \& Price, D. J. 2014a, ApJL, 790, L34
1504.00393.extracted_bib
[Martin et al. (2014b)]{Martin14b Martin, R. G., Nixon, C., Lubow, S. H., et al. 2014b, \apj, 792, L33
1504.00393.extracted_bib
[Mazeh et al. (1997)]{Mazeh97 Mazeh, T., Krymolowski, Y., \& Rosenfeld, G. 1997, ApJL, 477, L103
1504.00393.extracted_bib
[Mazeh \& Shaham (1979)]{Mazeh79 Mazeh, T., \& Shaham, J. 1979, A\&A, 77, 145
1504.00393.extracted_bib
[Miller \& Hamilton (2002)]{Miller02 Miller, M. C., \& Hamilton, D. P. 2002, \apj, 576, 894
1504.00393.extracted_bib
[Naoz et al. (2011)]{Naoz11 Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., \& Teyssandier, J. 2011, \nat, 473, 187
1504.00393.extracted_bib
[Naoz et al. (2013a)]{Naoz13a Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., \& Teyssandier, J. 2013a, \mnras, 431, 2155
1504.00393.extracted_bib
[Naoz et al. (2013b)]{Naoz13b Naoz, S., Kocsis, B., Leob A., \& Yunes, N. 2013b, \apj,773, 187
1504.00393.extracted_bib
[Natarajan \& Pringle (1998)]{Natarajan98 Natarajan, P., \& Pringle, J. E. 1998, ApJL, 506, L97
1504.00393.extracted_bib
[Nesvorn\'{y et al. (2003)]{Nesvorny03 Nesvorn\'{y, D., Alvarellos, J. L. A., Dones, L., \& Levison, H. F. 2003, \aj, 126, 398
1504.00393.extracted_bib
[Nixon et al. (2011a)]{Nixon11a Nixon, C. J., Cossins, P. J., King, A. R., Pringle, J. E. 2011a, MNRAS, 412, 1591
1504.00393.extracted_bib
[Nixon et al. (2011b)]{Nixon11b Nixon, C. J., King, A. R., Pringle, J. E. 2011b, MNRAS, 417, L66
1504.00393.extracted_bib
[Nixon (2012)]{Nixon12a Nixon, C. J. 2012, MNRAS, 423, 2597
1504.00393.extracted_bib
[Nixon \& King (2012)]{Nixon12b Nixon, C. J., \& King, A. R. 2012, MNRAS, 421, 1201
1504.00393.extracted_bib
[Nixon et al.(2013)]{Nixon13 Nixon, C., King, A., \& Price, D. 2013, MNRAS, 434, 1946
1504.00393.extracted_bib
[Offner et al. (2010)]{Offner10 Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz, M. R., \& Klein, R. I. 2010, \apj, 725, 1485 %
1504.00393.extracted_bib
[Ogilvie (1999)]{Ogilvie99 %Ogilvie, G. I. 1999, \mnras, 304, 557
1504.00393.extracted_bib
[Ogilvie (2001)]{Ogilvie01a Ogilvie, G. I. 2001, \mnras, 325, 231
1504.00393.extracted_bib
[Ogilvie \& Barker (2014)]{Ogilvie14 Ogilvie, G. I., \& Barker, A. J. 2014, \mnras, 445, 2621
1504.00393.extracted_bib
[Ogilvie \& Dubus (2001)]{Ogilvie01b Ogilvie, G. I., \& Dubus, G. 2001, \mnras, 320, 485
1504.00393.extracted_bib
[Paczy\'{nski (1977)]{Paczynski77 Paczy\'{nski, B. 1977, ApJ, 216, 822
1504.00393.extracted_bib
[Papaloizou (2005a)]{Papaloizou05a Papaloizou, J. C. B. 2005a, A\&A, 432, 743
1504.00393.extracted_bib
[Papaloizou (2005b)]{Papaloizou05b Papaloizou, J. C. B. 2005b, A\&A, 432, 757
1504.00393.extracted_bib
[Papaloizou \& Pringle (1983)]{Papaloizou83 Papaloizou, J. C. B., \& Pringle, J. E. 1983, \mnras, 202, 1181
1504.00393.extracted_bib
[Papaloizou \& Terquem (1995)]{Papaloizou95 Papaloizou, J. C. B., \& Terquem, C. 1995, \mnras, 274, 987
1504.00393.extracted_bib
[Petrovich (2015)]{Petrovich15 Petrovich, C. 2015, \apj, 799, 27
1504.00393.extracted_bib
[Porter \& Rivinus (2003)]{Porter03 Porter, J M., \& Rivinius, T. 2003, PASP, 115, 1153
1504.00393.extracted_bib
[Price (2007)]{Price07 Price, D. J. 2007, \pasa, 24, 159
1504.00393.extracted_bib
[Price \& Federrath (2010)]{Price10 Price, D. J., \& Federrath, C. 2010, \mnras, 406, 1659
1504.00393.extracted_bib
[Price (2012)]{Price12 Price, D. J. 2012, JCoPh, 231, 759 %
1504.00393.extracted_bib
[Pringle (1981)]{Pringle81 %Pringle, J. E. 1981, ARA\&A, 19, 137 %
1504.00393.extracted_bib
[Pringle (1992)]{Pringle92 %Pringle, J. E. 1992, MNRAS, 258, 811
1504.00393.extracted_bib
[Pringle (1996)]{Pringle96 Pringle, J. E. 1996, \mnras, 281, 357
1504.00393.extracted_bib
[Rice (2015)]{Rice15 Rice, K. 2015, \mnras, 448, 1729
1504.00393.extracted_bib
[Storch et al. (2014)]{Storch14 Storch, N., Kassandra, R. A., \& Lai, D. 2014, Science, 345, 1317
1504.00393.extracted_bib
[Storch et al. (2015)]{Storch15 Storch, N., \& Lai, D. 2015, \mnras, 448, 1821
1504.00393.extracted_bib
[Shakura \& Sunyaev (1973)]{SS73 Shakura, N. I., \& Sunyaev, R. A., 1973, A\&A, 24, 337
1504.00393.extracted_bib
[Stapelfeldt et al. (1998)]{Stapelfeldt98 Stapelfeldt, K. R., Krist, J. E., M\'{enard, F., et al. 1998, ApJL, 502, L65
1504.00393.extracted_bib
[Takeda \& Rasio (2005)]{Takeda05 Takeda, G., \& Rasio, F. A. 2005, \apj, 627, 1001
1504.00393.extracted_bib
[Teyssandier et al. (2013)]{Teyssandier13 Teyssandier, J., Naoz, S., Lizarraga, I., \& Rasio, F. A. 2013, \apj, 779, 166
1504.00393.extracted_bib
[Wijers \& Pringle (1999)]{Wijers99 Wijers, R. A. M. J., \& Pringle, J. E. 1999, \mnras, 308, 207
1504.00393.extracted_bib
[Williams et al. (2014)]{Williams14 Williams, J. P., Mann, R. K., Di Francesco, J., et al. 2014, \apj, 796, 120
1504.00393.extracted_bib
[Winn \& Fabrycky (2014)]{Winn14 Winn, J. N., \& Fabrycky, D. C. 2014, arXiv:1410.4199
1504.00393.extracted_bib
[Wu \& Murray (2003)]{Wu03 Wu, Y., \& Murray, N. 2003, \apj, 589, 605
1504.00393.extracted_bib
[Xiang-Gruess \& Papaloizou (2013)]{Xiang13 Xiang-Gruess, M., \& Papaloizou, J. C. B. 2013, \mnras, 431, 1320