The Kozai-Lidov Mechanism in Hydrodynamical Disks

We use three dimensional hydrodynamical simulations to show that a highly misaligned accretion disk around one component of a binary system can exhibit global Kozai-Lidov cycles, where the inclination and eccentricity of the disk are interchanged periodically. This has important implications for accreting systems on all scales, for example, the formation of planets and satellites in circumstellar and circumplanetary disks, outbursts in X-ray binary systems and accretion on to supermassive black holes.


INTRODUCTION
Disks that orbit objects in a variety of types of astrophysical binary systems can sometimes be misaligned with respect to their binary orbital planes. When material that is misaligned with the binary orbit is accreted, a misaligned accretion disk can form around the primary and/or the secondary masses. Supermassive black hole (SMBH) binaries are likely to accrete material in a chaotic fashion (King & Pringle 2006, 2007 and therefore highly inclined disks around each black hole are expected . Also, warped disks have been observed in active galactic nucleus (AGN) maser disks (e.g., Caproni et al. 2006;Martin 2008).
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 and become the sites of planet formation (Bate et al. 2010). 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 (Hale 1994). More speculatively, highly misaligned disks could form around young giant planets if the planet's inclination to the protoplanetary disk is sufficiently high.
Evolved binary star systems such as low mass X-ray binaries (Nixon & Salvesen 2014), microquasars (e.g., Maccarone 2002;Martin et al. 2008aMartin et al. , 2008b, and Be/X-ray binaries (Martin et al. 2011) are thought to have disks that are misaligned with the binary orbit. When a star in a binary star system undergoes an asymmetric supernova explosion, the explosion can leave the spin of the unexploded star in a highly misaligned state with respect to the binary orbit. In this configuration, a misaligned disk may form from material ejected by the unexploded star (Martin et al. 2009(Martin et al. , 2010. Kozai-Lidov (KL) oscillations occur in highly misaligned test particle orbits around one component of a binary, where the 6 Sagan Fellow. 7 Einstein Fellow.
particle's inclination is periodically exchanged for eccentricity (Kozai 1962;Lidov 1962). During this process, the component of the angular momentum that is perpendicular to the binary orbital plane is conserved, which is expressed as where i p is the inclination of the particle orbital plane relative to the binary orbit plane and e p is the eccentricity of the test particle. A test particle that is initially on a circular and highly misaligned orbit undergoes oscillations of its orbital plane involving closer alignment (higher values of | cos i p |) and therefore higher values of its eccentricity e p . For these oscillations to occur, the initial inclination of the test particle orbit i p0 must satisfy the condition that cos 2 i p0 < cos 2 i cr = 3/5. This condition requires that 39 • i p0 141 • . Furthermore, inclination values i p during the oscillations are bounded by the condition that cos 2 i p0 cos 2 i p cos 2 i cr . Following from Equation (1), the maximum eccentricity that an initially circular particle orbit can achieve is given by (e.g., Innanen et al. 1997). The KL mechanism for orbiting objects has been extensively studied in the literature. For example, it can occur for asteroids (Kozai 1962), artificial satellites (Lidov 1962), triple star systems (Eggleton & Kiseleva-Eggleton 2001;Fabrycky & Tremaine 2007), planet formation with inclined stellar companions (Wu & Murray 2003;Takeda & Rasio 2005), inclined planetary companions (Nagasawa et al. 2008), merging SMBHs (Blaes et al. 2002), stellar compact objects (Thompson 2011), and blue straggler stars (Perets & Fabrycky 2009). Batygin et al. (2011) and Batygin (2012) investigated the evolution of selfgravitating, but pressureless and inviscid, misaligned disks in binary systems. However, no KL oscillations were found. Recently, Terquem & Ajmia (2010) showed that an external gas Figure 1. Eccentricity and inclination evolution of a test particle around the primary of an equal mass binary. The particle is initially at a radius d = 0.2 a, in a circular orbit, from the primary mass at an inclination i p0 = 60 • . disk can induce Kozai oscillations in the orbit of a misaligned companion. However, to our knowledge, the effect of the KL mechanism acting on a hydrodynamical disk has not yet been investigated. In this work we explore this process with threedimensional hydrodynamical simulations of misaligned disks in binary systems.

TEST PARTICLE ORBITS
We first consider ballistic particle orbits around the primary of a circular orbit binary system. The primary has mass M 1 , the secondary has mass M 2 , and they orbit at a separation a. The total mass of the binary is M = M 1 + M 2 .
In Figure 1 we show the eccentricity and inclination evolution of a particle that has an initial distance of d = 0.2 a from the primary and inclination i p0 = 60 • to the binary orbital plane. The maximum eccentricity reached is similar to that predicted by Equation (2), e max ≈ 0.76. The analytic period for KL cycles is (e.g., Kiseleva et al. 1998), where the orbital period of the binary is P b = 2π/Ω b , the eccentricity of the binary is e b (in our circular binary case, e b = 0), and the orbital period of the particle about the primary is P p = 2π/Ω p . However, we find in our test particle simulations that the timescale τ KL depends on the initial inclination of the test particle. This dependence is not included in this formula. We determined that Equation (3) is valid up to a factor of a few, due to this inclination dependence. With M 2 /M 1 = 1 and d = 0.2 a, we find in Equation (3) that τ KL = 15.8 P b , similar to that displayed in the particle orbits in Figure 1. We find that the oscillatory behavior occurs only for inclinations i 40 • , in line with the KL mechanism.
In the next section we investigate the response of a hydrodynamical (pressure and viscous internal forces) disk that satisfies the criteria for KL oscillations to occur on a test particle. As far as we know, this is the first time that this has been investigated.

HYDRODYNAMICAL DISK SIMULATIONS
In this section we consider the evolution of a highly misaligned fluid disk around one component of a circular equal mass binary. We use the smoothed particle hydrodynamics (SPH; e.g., Price 2012) code phantom (Price & Federrath 2010;Lodato & Price 2010). Misaligned accretion disks in binary systems have been modeled previously with phantom (e.g., Nixon 2012;Nixon et al. 2013;Martin et al. 2014). The binary and disk parameters are summarized in Table 1. The equal mass binary, with total mass M = M 1 + M 2 , has a circular orbit in the x-y plane with separation a. We choose the accretion radius for particle removal from the simulation about each object to be 0.025 a. Figure 2 shows the initially flat and circular but tilted disk. The disk has a mass of 10 −3 M with 10 6 particles and is inclined by 60 • to the binary orbital plane. The value of the disk mass has no dynamical significance in the calculation, since the disk self-gravity is ignored, and has too low mass to affect the binary orbit for the timescale simulated. The initial surface density of the disk has a power-law distribution Σ ∝ R −3/2 between R in = 0.025 a and R out = 0.25 a. The outer radius of the disk is chosen to be the tidal truncation radius for the disk assuming a coplanar binary (Paczynski 1977). However, misaligned disks feel a weaker binary torque and thus the outer truncation radius can be much larger than this value (C. Nixon et al., in preparation). We take a locally isothermal disk with sound speed c s ∝ R −3/4 and H/R = 0.02 at R = R out . This is chosen so that both α and h /H are constant over the disk (Lodato & Pringle 2007). The Shakura & Sunyaev (1973) α parameter varies in the small range 0.1-0.12 over the disk (we implement the disk viscosity in the usual manner by adapting the SPH artificial viscosity according to the procedure described in Lodato & Price (2010), using α AV = 1.91 and β AV = 2.0). The disk is resolved with shell-averaged smoothing length per scale height h /H ≈ 0.52.
In Figure 3 we show the time evolution of the eccentricity and inclination of the disk at two radii from the primary, d = 0.1 a and d = 0.2 a. The figure clearly shows damped KL oscillations of the disk. As the eccentricity increases, the inclination decreases and vice versa. In Figure 4 we show the disk at the maximum eccentricity at a time of t = 11 P b . Since there is dissipation within a disk (that is not present in a particle orbit), the eccentricity does not reach the maximum value of 0.76 predicted by Equation (2) (see also Figure 1). In addition, the magnitude of the oscillations decays in time. According to Equation (3) the local oscillation timescales τ KL should differ   by a factor of about 2.8 at the two radii considered. However, we see by comparing the two plots that the evolution of the disk at these two radii is very similar in magnitude and timescale. Thus, the disk is undergoing global KL oscillations.
In our simulation the disk has a small radial extent (R out /R in ≈ 10), whereas in an astrophysical situation the inner disk radius is often much smaller. For the simulations in this work, the inner boundary is a circular mass sink. If the disk extends inward to the surface of a star, then a hard wall condition is appropriate. In that case, the inner boundary condition is that the disk eccentricity vanishes (see Lubow 2010). Another possible inner boundary condition is that the radial derivative of the eccentricity vector vanishes. We discuss this further in Section 4.
The realignment of this highly misaligned disk (which is not accreting any new material) proceeds in two stages. In the first stage, the KL mechanism drives a rapid decay down to i = i cr ≈ 40 • . The rapidity of the decay depends on the parameters we have adopted and should be further explored in the future. In the second stage, the long term evolution of the disk proceeds due to secular processes. In this stage, the alignment torque is due to viscous forces that interact with a disk warp. The warp is due to binary gravitational torques that act to precess the disk differentially, rather than the KL torques of the first stage. The second stage alignment occurs on a much longer timescale than we have simulated (e.g., . We have considered a range of disk parameters and we find qualitatively the same disk behavior. For example, when we reduce α by an order of magnitude to 0.01 (more relevant to protoplanetary disks), the KL oscillations are longer-lived due to weaker dissipation. But shock damping may also play a role. We have also considered a counterrotating disk with initial inclination i = 120 • . In this case, the eccentricity growth is similar to that in Figure 3, but the inclination of the disk evolves toward counteralignment rather than alignment. There are a wide range of disk and binary parameters that should be explored in future work. In the next section we make some order-of-magnitude estimates for binary and disk parameters for which the KL mechanism is important.

DISK TIMESCALES
In this section we consider critical disk and binary parameters for which the KL mechanism operates. For the KL mechanism to operate, the apsidal precession due to the internal disk forces cannot dominate over the binary gravitational induced apsidal rate. Otherwise, the KL effects cancel over time.

Pressure
The precession rate associated with pressure can be quite small (Lubow 1992;Goodchild & Ogilvie 2006). The pressureinduced precession rate |Ω − κ| ∼ (H/R) 2 Ω, where H/R and Ω are evaluated at some point within the disk. Furthermore, if the sound crossing timescale radially across the disk is shorter than the period of KL oscillations, the disk can communicate globally and undergo a coherent large-scale response (see, e.g., Larwood et al. 1996, for the case of differential precession). This condition can be written as 1/τ KL (H/R)Ω, while the apsidal rate is smaller ∼(H/R) 2 Ω. For the parameters in the simulation in Figure 3, the estimated apsidal rate is smaller than the disk KL precession rate, and the coherence constraint is satisfied. Consequently, the required conditions for KL global oscillations to operate in our simulations appear to be satisfied.

Self-gravity
Apsidal disk precession can also occur due to the effects of self-gravity in the gas, but we have not included this effect in our simulations. We determined the local apsidal precession rate due to self-gravity in a disk with surface density profile Σ ∝ 1/R q for q = 1 that extends from R = 0 to R = 0.35 a using Equations (6) and (7) of Lubow & Ogilvie (2001) with a smoothing parameter H whose square is added inside the square root of Equation (5) in that paper. We identify H with the disk thickness. The (global) disk precession rate is obtained by taking the angular-momentum-weighted average of the local rates (Larwood & Papaloizou 1997;Lubow & Ogilvie 2001). Based on this result, we crudely estimate the condition for suppressing KL by self-gravity as M d b(H/R)M 1 M 2 /M, where we estimate b ∼ 1/3. b is fairly insensitive to q for 0 < q < 1.5. Disks where self-gravity is weak are relevant to a wide variety of astrophysical scenarios such as protoplanetary, Be star, X-ray binary, circumplanetary, and AGN disks. However, if the disk self-gravity is sufficiently strong, as sometimes occurs in the early evolution of disks around young stars, the KL oscillations may be suppressed (Batygin 2012).

Global Oscillation Timescale
Since the disk responds globally, we apply the theory of rigid disks (Larwood & Papaloizou 1997;Lubow & Ogilvie 2001) and estimate that the global disk response period of inclination oscillations is Given that Equation (3) for τ KL has an inclination dependence and is therefore only valid up to a factor of a few, this equation is also only valid up to a factor of a few. If we assume the surface density is a power law in radius, Σ ∝ R −p , then we find for R out R in . With M 1 = M 2 = 0.5 M and p = 1.5, we find τ KL = 17.1(0.35a/R out ) 3/2 P b . We have normalized the outer radius by 0.35 a since the disk expands to approximately this value during the oscillations. We find that this estimate for τ KL is commensurate with the timescale observed initially in the simulation in Figure 3 of τ KL ≈ 16 P b .

DISCUSSION
The applications of the KL disk oscillations to various astrophysical systems depends on how these oscillations behave under different disk conditions more generally than we have considered here. In particular, the nature of the longevity of the oscillations needs to be understood. If the results found here hold generally, then the oscillations are damped after a few dozen binary orbital periods and the disk attains a somewhat eccentric state with inclination at the critical value for the KL oscillations, subject to a longer timescale decay by viscous dissipation. But the timescale for the KL decay may not be the same for warmer disks, less viscous disks, or more extreme mass ratio binaries.
In the case of Be/X-ray binaries, the disks are transient and the KL oscillations may play an important role. In Martin et al. (2014), we found significant eccentricity growth in a highly misaligned circumprimary disk of a Be star in an eccentric binary with a neutron star companion. At the time of writing we had only investigated eccentricity growth due to the eccentric companion in the coplanar case. However, the results of this Letter indicate that the eccentricity growth we found for the Be star disk is due to the KL effect. The KL effect explains why the eccentricity growth in the disk was only present for large inclination angles, and why its strength did not decay with increasing tilt angle. The eccentricity growth is a key ingredient of our Type II X-ray outburst model. Thus, we would expect Type II outbursts only in Be/X-ray binary systems that have a Be star spin (or disk plane) misalignment in the range 40 • -140 • .
In the case of a Be star, the ratio of the outer radius of the disk to star is not very large, around 7. In that case, an inner boundary condition, such as requiring that the eccentricity vanish, will play some role in the determining the structure of the disk eccentricity and may have some influence on the evolution of the global structure of the disk.
SMBH binary disks may also be susceptible to KL oscillations. It is often expected that SMBH binaries accrete gas from circumbinary disks and it is possible that these circumbinary disks form highly misaligned with the binary orbit (Nixon et al. 2011a(Nixon et al. , 2011b. It is therefore also possible that misaligned circumprimary and circumsecondary disks can form (see, e.g., Figures 6 and 7 of Nixon et al. 2013) and be subject to KL cycles. In this case, the strong increase in density at the pericenter of the disk orbit could result in strong dissipation or enhanced star formation.
Since binary stars are common, the disk KL mechanism presented here could play a role in the process of planet formation around the components of a binary. In particular, we would like to understand how planets that are apparently undergoing KL oscillations in wide binaries could form in protostellar disks (e.g., Wu & Murray 2003;Takeda & Rasio 2005). According to our calculations, the disk tilt damps after KL oscillations to the critical angle. A planet forming in such a disk after damping would not undergo KL oscillations. They must have formed in a disk that avoided evolution that damped the inclination to the critical inclination angle i cr at the time of their formation. Based on the results of this Letter, a circular binary with separation a = 10 3 AU, with each star having 1 M , would induce KL oscillations in a ∼350 AU disk with oscillation period ∼3 × 10 5 yr. The KL decline in tilt could be avoided, for example, if the disk self-gravity is sufficiently strong (Batygin 2012). This issue requires further study.

CONCLUSIONS
We have found that the KL mechanism, which exchanges inclination for eccentricity in a highly misaligned particle orbit around a component of a binary system, can also operate in a fluid disk. This result has implications for a range of astrophysical systems. Much work remains to be done to understand its general behavior. Simulations should be performed for different system parameters. It would also be desirable to develop a linear model for the disk evolution involving a disk whose initial inclination is just above i cr by extending existing linear models for eccentricity and inclination evolution. Support for R.G.M. was provided under contract with the California Institute of Technology (Caltech) funded by NASA through the Sagan Fellowship Program. Support for C.J.N. was provided by NASA through the Einstein Fellowship Program, grant PF2-130098. S.H.L. acknowledges support from NASA grant NNX11AK61G. P.J.A. acknowledges support from NASA's ATP program under awards NNX11AE12G and NNX14AB42G. D.J.P. is supported by Future Fellowship FT130100034 from the Australian Research Council. We acknowledge the use of SPLASH (Price 2007) for the rendering of the figures. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794), the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research. The Janus supercomputer is operated by the University of Colorado Boulder.