The Kozai–Lidov Mechanism in Hydrodynamical Disks. III. Effects of Disk Mass and Self-Gravity

The Kozai-Lidov Mechanism in Hydrodynamical Disks. III. Effects of Disk Mass and Self-Gravity


Previously we showed that a substantially misaligned viscous accretion disk with pressure that orbits around one component of a binary system can undergo global damped Kozai-Lidov (KL) oscillations. These oscillations produce periodic exchanges of the disk eccentricity with inclination. The disk KL mechanism is quite robust and operates over a wide range of binary and disk parameters. However, the effects of self-gravity, which are expected to suppress the KL oscillations for sufficiently massive disks, were ignored. Here, we analyze the effects of disk self-gravity by means of hydrodynamic simulations and compare the results with the expectations of analytic theory. The disk mass required for suppression in the simulations is a few percent of the mass of the central star and this roughly agrees with an analytical estimate. The conditions for suppression of the KL oscillations in the simulations are close to requiring that the disk be gravitationally unstable. We discuss some implications of our results for the dynamics of protoplanetary disks and the related planet formation.

accretion, accretion disks – binaries: general – hydrodynamics – planets and satellites: formation

1 Introduction

Estimates suggest that 40% to 50% of observed exoplanets are in binary systems (Horch et al., 2014). As the birthplace of the planets, protoplanetary disks are fundamental to explaining planet formation and thus it is important to understand their evolution in a binary system. Recent ALMA observations revealed large mutual inclinations (greater than ) between the circumstellar disks around binary system components (e.g., Jensen & Akeson, 2014; Williams et al., 2014). Although the binary orbital planes in these systems have not yet been identified, it is possible that at least one of the disks in each system is significantly inclined ( ) with respect to the binary orbit.

The Kozai-Lidov (KL) mechanism is a well known process that occurs when a ballistic object orbits around one component of a binary system and the object’s orbital plane is substantially misaligned () with respect to the binary orbital plane (Kozai, 1962; Lidov, 1962). In this process, the orbital eccentricity and inclination of the object undergo periodical exchange. An object on an initially circular orbit periodically gains and loses eccentricity. There have been a large number of works on this mechanism since its first discovery in the 1960s (e.g., Holman et al., 1997; Innanen et al., 1997; Kiseleva et al., 1998; Ford et al., 2000; Lithwick & Naoz, 2011; Naoz et al., 2013a; Liu et al., 2015; Antognini, 2015) and the KL mechanism has found applications in various astronomical processes, most notably here on the formation of hot Jupiter planets (e.g., Wu & Murray, 2003; Takeda & Rasio, 2005; Fabrycky & Tremaine, 2007; Naoz et al., 2013b; Petrovich, 2015; Rice, 2015).

While the KL cycle of orbiting objects has been extensively studied, only very recently was it found to operate on a fluid disk with pressure and viscosity by means of hydrodynamic simulations, (Martin et al., 2014, hereafter Paper I). Thus, an initially circular disk can become highly eccentric if its initial tilt exceeds about with respect to the binary orbital plane (Fu et al., 2015, hereafter Paper II). Due to the efficient radial communication via either disk pressure or viscosity, the disk KL cycle for a typical protostellar disk operates in a global, coherent fashion, and the disk remains quite flat (unwarped) throughout the process. The eccentricity growth is fairly uniform in radius across the disk. Unlike the ballistic object case, the disk oscillations are damped due to viscous dissipation, likely involving shocks in the disk. When the oscillations fully damp, the disk is tilted at the critical KL angle () or somewhat less and the disk is approximately circular. After this stage, the circular disk evolves towards alignment with the binary orbital plane through tidal dissipation effects associated with turbulent viscosity (e.g., King et al., 2013). Paper II extended the work of Paper I by surveying a large range of binary and disk parameter space, including the initial disk inclination, temperature, viscosity, size, surface density profile, and the binary mass ratio and eccentricity. The disk KL cycle is a fairly robust process that can occur under many different binary-disk conditions. However, Paper II did not consider the effects of self-gravity.

At early times in the protoplanetary disk evolution, the disk mass is large, several percent of the mass of the central star, and the disk may experience dynamical effects of self-gravity. Depending on the properties of the disk, it is possible that self–gravity may be sufficiently strong to prevent KL oscillations from operating. Instead, the disk will globally precess as a circular object (as it does for lower inclination disks) (Batygin et al., 2011; Batygin, 2012; Lai, 2014). This suppression is expected when the disk apsidal precession rate due to self–gravity dominates over the apsidal precession rate due to the gravitational effects of the binary that cause the KL oscillations. The suppression may not continue as the disk evolves. As material drains on to the central star, the self-gravity weakens and the KL cycles can begin if the disk remains sufficiently misaligned.

Differential precession due to the binary acts to disrupt a disk. For a disk to precess coherently, internal torques are required to operate globally and communicate faster than the differential precession timescale (Larwood et al., 1996). The disk model of Batygin (2012) maintained disk coherence by means of disk self-gravity and omitted the effects of pressure and viscosity. Xiang-Gruess & Papaloizou (2014) pointed out that pressure alone allows a disk to precess coherently and thus self–gravity is not necessarily required. In this paper, for the first time we present results on the long term evolution of misaligned self-gravitating disks that have pressure and viscosity.

2 Numerical Simulations

In this section, we describe three-dimensional simulations of an inclined fluid disk around one component of an equal mass binary system, where the binary orbit is circular. We use the smoothed particle hydrodynamic (SPH) code PHANTOM (Lodato & Price, 2010; Price & Federrath, 2010; Nixon, 2012; Nixon & King, 2012; Price, 2012; Nixon et al., 2013). The simulation setup is similar to that in Paper II. An initially circular disk orbits around a central binary component of mass under the influence of the perturber binary component of mass . The disk and stars interact gravitationally and respond to these interactions. In particular, the orbit of the binary is affected by its gravitational interactions with the disk, although the changes to the binary orbit are small. The orbital plane of the disk is initially inclined to the binary orbital plane by . We adopt a locally isothermal equation of state and an explicit accretion disk viscosity. We include a nonlinear term with a coefficient (AV stands for artificial viscosity) in order to suppress interparticle penetration, as is standard in SPH codes. The disk sound speed is and the initial surface density distribution is , such that both (Shakura & Sunyaev, 1973) and the smoothing length are constant over the disk radius, (Lodato & Pringle, 2007). We start the simulations with SPH particles in the disk. The PHANTOM code adopts a cubic spline kernel as the smoothing kernel. The number of neighbors is roughly constant at . The disk initially extends from radius to radius . The inner boundary of the simulated region is set to the initial disk inner disk radius . As particles move to , they are removed from the simulation while their mass and momentum are added to the central star. We also impose an inner boundary radius around the perturbing companion, since some disk mass can be transferred to that component. In terms of the binary orbital size , the initial disk inner and outer radii in our simulations are located at and , respectively. The initial outer disk radius is chosen to be the tidal truncation radius of a coplanar disk (Paczynski, 1977). However, because the tidal torques on a misaligned disk are weaker (Lubow, Martin & Nixon, 2015), the disk initially expands somewhat beyond this radius.

We describe simulations with two sets of parameters for the disk aspect ratio and viscosity. The first is and and the second is and ( is evaluated at disk inner edge ). The disk is resolved with shell-averaged smoothing length per scale height and , respectively. For each set of and , we consider three different initial disk masses , and (in units of the total binary mass ). In Papers I and II, we ignored disk self-gravity given the low disk mass used there (). In this paper, we take into account the effect of disk self-gravity and compare simulations with and without self-gravity for the same disk mass. The algorithm for the SPH implementation of self-gravity in PHANTOM is described in Price & Monaghan (2007), which discusses numerical tests based on the radial oscillations and the static structure of a polytrope. The self-gravity algorithm of PHANTOM has also been used to study the formation of giant molecular clouds (Dobbs, 2008) and star clusters (Price & Bate, 2008).

Binary and Disk parameters Symbol Value
Mass of each binary component 0.5
Binary orbital eccentricity 0
Initial number of particles
Initial disk mass [0.001, 0.01, 0.03]
Initial disk outer radius 0.25
Initial disk inner radius 0.025
Mass accretion radius 0.025
Disk viscosity parameter [0.01, 0.1]
Disk aspect ratio [0.06, 0.1]
Initial disk surface density profile 1.5
Initial disk inclination
Table 1: Parameters of the SPH simulations for equal mass binary systems with total mass of and separation of
Figure 1: Evolution of the eccentricity (upper row), inclination (middle row), and total mass (bottom row) of the circumprimary disk. The left column is for simulations without disk self-gravity, whereas the right column is for simulations that take into account disk self-gravity. The eccentricity and inclination are measured at a radius of (the initial disk radial range is ), where is the binary semi–major axis and is the binary orbital period. Within each plot, each line represents a simulation with a different initial disk mass. The units of mass are that of the total binary mass, . All lines are averaged over one binary orbital period. The initial number of SPH particles is . The disk aspect ratio at the inner edge is . The disk viscosity parameter is and the initial disk inclination is .

Figure 1 shows the evolution of the disk eccentricity (upper row), inclination (middle row), and mass (bottom row) for three different initial disk masses. The eccentricity and inclination are measured at a representative radius , since the disk remains flat and the eccentricity is roughly constant with radius. The left and right columns show simulations without and with disk self-gravity, respectively. In all cases, disk gravity acts on the stars and slightly affects the binary orbit. In post-processing the simulation outputs, we compute the orbital eccentricity and inclination for each particle using its position and velocity information, then divide the disk into radial bins and calculate the mean properties of the particles within each bin. Results for the eccentricity and inclination are taken from the radial bin centered on radius . This is similar to the method used in Papers I and II, except now we average the properties over one binary orbit to smooth out the fluctuations on the binary orbital timescale. Note that the initial disk eccentricity is non-zero () in the upper row, even though the disk is initially set up to be circular. This non-zero eccentricity value is a consequence of the way we calculate the particle eccentricity and inclination (Equations [7] and [8] in Paper II) that treats particles as ballistic, i.e., assuming particles only feel the gravitational force of the central object, whereas particles actually are also influenced by the disk pressure force. Therefore, the non-zero initial disk eccentricities we see here are purely an artifact of our ignoring the disk pressure in the orbital elements calculation. The effect is more prominent here than in the simulations shown in Papers I and II because here the disk has higher temperature, (compared with ). However, it is still small enough not to affect our general interpretation, since we are focused on the major changes of the disk properties during KL oscillations.

As seen in the left column of Figure 1, where disk self-gravity is not included, increasing the disk mass does not significantly affect the disk orbital elements during the KL cycle. The period and amplitude of the oscillations are similar, although there is some variation in the time of the first peak. The difference in the onset of the initial KL cycle is due to the fact that disk mass slightly affects the companion’s orbit. The timing of the onset is sensitive to the initial conditions. Even with , the peak eccentricity value of is almost the same as that for the simulations with and with a peak . As seen in the right column of the figure, for the lower mass disks, self-gravity generally has little effect on disk eccentricity and inclination. However, the effects of the KL cycle are significantly reduced at the highest disk mass of , which has an initial disk minimum Toomre Q . The peak eccentricity in that case is compared to with and . The bottom row shows that by the end of the simulation at time , the disk in general has lost more than half of its initial mass. The majority of the lost mass is accreted on to the central object. A small fraction of the mass is ejected from the disk, mostly ending up around the companion binary star. Comparing simulations with the same initial disk mass, the run with self-gravity has more mass remaining at the end of the simulation than the one without self-gravity. This effect may be a consequence of the disk self-gravity helping to retain particles, and also reducing the particle loss at inner disk boundary, since the KL cycle is weakened.

Figure 2: View of the disk towards the plane for the two simulations with , and , with and without self–gravity. The binary orbit is in the plane (i.e., the perturbing object moves into and out of the page). The left panel shows the initial disk setup that is the same for both simulations. The middle (right) panel shows the disk without (with) self-gravity after it has undergone nodal precession by . The central mass is denoted by the black dot. The color coding is for the logarithm base of the column density (i.e., density integrated along the line of sight) in units of . (Color online)

Figure 2 shows edge-on views of the disk in the two simulations with , one with and the other without disk self-gravity. The left panel shows the initial conditions that were applied to both simulations. The middle and right panels are edge-on views after the disks have undergone nodal precession by . The middle panel is from the run without disk self-gravity (the blue curves in the left column of Figure 1) at a time of , while the right panel is from the run that includes disk self-gravity (the blue curves in the right column of Figure 1) at a time of . These runs have different disk nodal precession rates and as a result the middle and right panels are shown at slightly different epochs. Without self-gravity, the middle panel shows a substantially lopsided disk due to the KL-driven eccentricity growth. In the right panel, with disk self-gravity, the disk is still eccentric, but not as eccentric as in the middle panel. Although at an earlier time () the disk has a slightly higher eccentricity, the average disk eccentricity is clearly lower when self-gravity is included (see comparison of the blue curves in the top two panels of Figure 1). The main point here is that in both panels, the disk remains fairly flat (no significant warping) due to the efficient radial communication (provided mainly by pressure in these cases; see discussions in Paper II).

Figure 3 is similar to Figure 1 except for a different disk aspect ratio, , and viscosity, . When disk self-gravity is not included, we see again that the disk mass has only a minor effect on the disk KL cycle (left column). The oscillation timescales and amplitudes are similar. As we discussed in Figure 1, the reason why the disk mass plays a role in the absence of self-gravity is that the disk mass slightly affects the orbit of the companion, which in turn affects the onset of the disk KL cycle. With self-gravity included, we see a greater suppression of the disk KL cycle by self-gravity than that in Figure 1. With , the KL mechanism barely operates. For the disk mass of , there is a significant reduction. The lowest mass disk remains relatively unaffected by self-gravity. Note that for the highest disk mass, the initial disk has a minimum Toomre Q , which puts the disk on the verge of becoming gravitationally unstable to non-axisymmetric disturbances. In this paper, we have restricted our attention to simulations in which we do not see any sign of gravitational instability, such as fragmentation or clump formation.

Figure 3: Similar to Figure 1 except for different disk aspect ratio and viscosity parameter .

3 Discussion and Conclusions

Self-gravity introduces a source of disk apsidal precession, in addition to that due to the binary companion that causes the KL oscillations. If the self-gravity induced apsidal disk precession rate is faster than the binary induced precession rate, then the disk KL cycle can be suppressed (e.g., Holman et al., 1997). We estimate the magnitude of the local self-gravity induced disk apsidal precession rate from Equation (1) of Tremaine (1998), which is consistent with the rate given in Batygin et al. (2011). We determine the binary companion induced disk local precession rate by approximating the companion’s potential as arising from a uniform ring of mass and radius , where is binary separation. Since the disk maintains its coherence, its global precession rate is crudely estimated as the angular-momentum-weighted average of the local rates (see Equation (9) of Paper II). For the equal mass binary system we consider, we take , (due to disk expansion from the initial ). The global disk precession rate caused by a binary companion is , where is the binary orbital frequency. This precession rate is similar in value to the nodal precession rate implied by Figure 2 that shows that the disks undergo half the nodal precession period at a time of about . The global disk precession rate due to self-gravity is . Requiring provides an estimate of the minimum disk mass to suppress KL oscillations as . According to Figure 7 of Batygin et al. (2011), the critical disk mass for suppressing KL oscillations is generally somewhat larger than this value, depending on the disk inclination. Based on that work, for the initial disk inclination we have used, the estimated disk mass for KL suppression is about a factor of 5 times larger, that is .

In the case of , the blue curve in the right column of Figure 1 clearly shows that the disk KL cycle is weakened for the case of an initial disk mass of . In spite of the reduced oscillation amplitude, the disk still undergoes KL cycles. There is an eccentricity peak and inclination valley at a time when the disk has a mass of about which we adopt as a crude estimate for the minimum mass for KL suppression.

In the case of , the blue curve in the right column of Figure 3 shows that the KL cycle is suppressed, for the case of an initial disk mass of . In the case of the lower initial disk mass of , the KL oscillation appears to be present, but weakened. A reasonable estimate of the critical disk mass suppression of KL oscillations is similar to our estimate for the other disk model, . Therefore, the simulations show that the critical disk mass for suppressing KL oscillations is in the range of ( 4% of the mass of the central star) that is about a factor of 5 larger than the value based on equating the self-gravity induced precession rate with the binary induced precession rate. It is in agreement with the value implied by Batygin et al. (2011) that is discussed above.

Papers I and II showed that global disk KL oscillations damp. After damping, the disk inclination is at or below critical KL angle. If a planet forms after this stage, then its initial orbital tilt is not large enough to trigger the KL mechanism. This reduced tilt poses a challenge to the KL model of planets which assumes initially highly inclined planet orbits. The tilt evolution of such planet–disk configurations requires further study.

Suppression of the KL cycle by means of self-gravity requires substantial disk masses. At the critical disk mass we found, the disk is quite close to being gravitationally unstable. In other words, we find a relatively narrow window of disk masses where the disk KL oscillation is significantly subdued, while the disk is still stable to gravitational instability.

The highest initial disk mass we consider in this paper, , shows more warping with self–gravity included than in cases without self–gravity. Initially, the inclination of the inner parts of the disk increases, while that of the outer parts decreases. In Paper II, we found that the critical (minimum) tilt angle for KL oscillations to operate in a disk is about , somewhat higher than the critical angle for free particles (). In this paper, we have limited our analysis to an initial disk tilt of that is just above the minimum value for KL oscillations. Such initial conditions produce much more moderate KL effects on a disk than occurs at higher tilt angles. We also reported in Paper II that strong density enhancements are found during KL oscillations of a non-self-gravitating disk with a larger initial tilt. These density enhancements appear to involve shocks. With larger initial tilts than we have considered here, a self-gravitating disk that undergoes KL oscillations may undergo disk fragmentation. Fragmentation of a disk is a possible mechanism for planet formation and may be aided by binary perturbations (Boss, 1997, 2006). The disk KL cycle may then provide an alternative means by which a binary companion can promote disk fragmentation/clumping. Such effects should be explored in the future.

W.F. and S.H.L. acknowledge support from NASA grant NNX11AK61G. Computing resources supporting this work were provided by the institutional computing program at Los Alamos National Laboratory. We thank Daniel Price for providing the PHANTOM code for SPH simulations and SPLASH code (Price, 2007) for data analysis and rendering of figures.


  1. Antognini, J. M. O. 2015, arXiv:1504.05957
  2. Batygin, K., Morbidelli, A., & Tsiganis, K. 2011, A&A, 533, A7
  3. Batygin, K. 2012, Nature, 491, 418
  4. Boss, A. P. 1997, Science, 276, 1836
  5. Boss, A. P. 2006, ApJ, 641, 1148
  6. Dobbs, C. L. 2008, MNRAS, 391, 844
  7. Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298
  8. Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385
  9. Fu, W., Lubow, S. H., & Martin, R. G. 2015, ApJ, 807, 75
  10. Horch, E. P., Howell, S. B., Everett M. E., Ciardi D. R., 2014, ApJ, 795, 60
  11. Holman, M. J., Touma, J., & Tremaine, S. 1997, Nature, 386, 254
  12. Innanen, K. A., Zheng, J. Q., Mikkola, S., & Valtonen, M. J. 1997, AJ, 113, 1915
  13. Jensen, E. L. N. & Akeson, R. 2014, Nature, 511, 567
  14. King, A. R., Livio, M., Lubow, S. H., & Pringle, J. E. 2013, MNRAS, 431, 2655
  15. Kiseleva, L. G., Eggleton, P. P., & Mikkola, S. 1998, MNRAS, 300, 292
  16. Kozai, Y. 1962, AJ, 67, 591
  17. Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., & Terquem, C. 1996, MNRAS, 282, 597
  18. Lai, D. 2014, MNRAS, 440, 3532
  19. Lidov, M. L. 1962, P&SS, 9, 719
  20. Lithwick, Y., & Naoz, S. 2011, ApJ, 742, 94
  21. Liu, B., Muñoz, D., & Lai, D. 2015, MNRAS, 447, 747
  22. Lodato, G., & Pringle, J. E. 2007, MNRAS, 381, 1287
  23. Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212
  24. Lubow S. H., Martin R. G., & Nixon C. J., 2015, ApJ, 800, 96
  25. Martin, R. G., Nixon, C., Lubow, S. H., et al. 2014, ApJL, 792, L33
  26. Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013a, MNRAS, 431, 2155
  27. Naoz, S., Kocsis, B., Leob A., & Yunes, N. 2013b, ApJ,773, 187
  28. Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946
  29. Nixon, C. J. 2012, MNRAS, 423, 2597
  30. Nixon, C. J., & King, A. R. 2012, MNRAS, 421, 1201
  31. Paczynski B., 1977, ApJ, 216, 822
  32. Petrovich, C. 2015, ApJ, 799, 27
  33. Price, D. J. 2007, PASA, 24, 159
  34. Price, D. J., & Bate, M. R. 2008, MNRAS, 385, 1820
  35. Price, D. J., & Monaghan, J. J. 2007, MNRAS, 374, 1374
  36. Price, D. J., & Federrath, C. 2010, MNRAS, 406, 1659
  37. Price, D. J. 2012, JCoPh, 231, 759
  38. Rice, K. 2015, MNRAS, 448, 1729
  39. Shakura, N. I., & Sunyaev, R. A., 1973, A&A, 24, 337
  40. Takeda, G., & Rasio, F. A. 2005, ApJ, 627, 1001
  41. Tremaine, S. 1998, ApJ, 116, 2015
  42. Williams, J. P., Mann, R. K., Di Francesco, J., et al. 2014, ApJ, 796, 120
  43. Wu, Y., & Murray, N. 2003, ApJ, 589, 605
  44. Xiang-Gruess, M. & Papaloizou, J. C. B. 2014, MNRAS, 440, 1179
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description