Kozai-Lidov Disc Instability

TL;DR

This paper analyzes the linear stability of Kozai-Lidov modes in tilted, inviscid accretion discs within binary systems, revealing conditions under which these discs become unstable to oscillations based on their aspect ratio and binary parameters.

Contribution

It introduces a 1D framework to study Kozai-Lidov instability in large-radius discs, extending understanding beyond previous multi-dimensional simulations.

Findings

KL instability occurs within specific disc aspect ratio windows.

Multiple unstable eccentric modes can coexist in a disc.

Instability can arise for any nonzero tilt angle depending on disc properties.

Abstract

Recent results by Martin et al. (2014) showed in 3D SPH simulations that tilted discs in binary systems can be unstable to the development of global, damped Kozai-Lidov (KL) oscillations in which the discs exchange tilt for eccentricity. We investigate the linear stability of KL modes for tilted inviscid discs under the approximations that the disc eccentricity is small and the disc remains flat. By using 1D equations, we are able to probe regimes of large ratios of outer to inner disc edge radii that are realistic for binary systems of hundreds of AU separations and are not easily probed by multi-dimensional simulations. For order unity binary mass ratios, KL instability is possible for a window of disc aspect ratios H/r in the outer parts of a disc that roughly scale as (n_b/n)^2 < H/r < n_b/n, for binary orbital frequency n_b and orbital frequency n at the disc outer edge. We present a framework for understanding the zones of instability based on the determination of branches of marginally unstable modes. In general, multiple growing eccentric KL modes can be present in a disc. Coplanar apsidal-nodal precession resonances delineate instability branches. We determine the range of tilt angles for unstable modes as a function of disc aspect ratio. Unlike the KL instability for free particles that involves a critical (minimum) tilt angle, disc instability is possible for any nonzero tilt angle depending on the disc aspect ratio.

Tables

Table 1. Table 1: Model Parameters

Model	$p$	$s$	$r_{\mathrm{in}}$
A1	1	0.75	0.01
A2	1	0.75	0.0001
B1	1	0.5	0.01
B2	1	0.5	0.0001

Equations

$$\frac{da}{dt}$$

$$\frac{de}{dt}$$

$$\frac{di}{dt}$$

$$\frac{d\omega }{dt}$$

$$\frac{d\Omega }{dt}$$

$$\frac{dE}{dt}=\frac{3G{M}_2}{8n{a}_{\mathrm{b}}^3}((2-3{\sin }^{2}i)\mathrm{i}E+(5{\sin }^{2}i)\mathrm{i}{E}^{\ast })-\left(\frac{d\Omega }{dt}\cos i\right)\mathrm{i}E.$$

$$A$$

$$B$$

$$\lambda =\frac{3\sqrt{2}G{M}_2}{4n{a}_{\mathrm{b}}^3}\sqrt{5-3{\cos }^{2}i}$$

$$\frac{di}{dt}\int hdm$$

$$\frac{d\Omega }{dt}\int hdm$$

$$\begin{array}{r}\mathrm{i}{\partial }_r\left(b_1(r){\partial }_r\left(\frac{E(r,t)}{c_{\mathrm{s}}^2}\right)\right)+\mathrm{i}{b}_2(r)E(r,t)\\ +\mathrm{i}{b}_3(r)E^{\ast }(r,t)=J(r){\partial }_{t}E(r,t).\end{array}$$

$$J(r)=2r^{3}n(r)\Sigma (r),$$

$$b_1(r)$$

$$b_2(r)$$

$$b_3(r)$$

$$A(r)$$

$$B(r)$$

$$\frac{d\Omega }{dt}$$

$${\partial }_{r}E(r_{\mathrm{in}},t)={\partial }_{r}E(r_{\mathrm{out}},t)=0.$$

$$E(r,t)=E_{+}(r)\exp (\lambda t)+E_{-}(r)\exp ({\lambda }^{\ast }t),$$

$$E_{+}^{\prime \prime }+c_1{E}_{+}^{\prime }+c_2{E}_{+}+c_3{E}_{-}^{\ast }=-\lambda c_4{E}_{+}$$

$$E_{-}^{\prime \prime }+c_1{E}_{-}^{\prime }+c_2{E}_{-}+c_3{E}_{+}^{\ast }=-{\lambda }^{\ast }{c}_4{E}_{-},$$

$$c_1(r)$$

$$c_2(r)$$

$$c_3(r)$$

$$c_4(r)$$

$$E_{+}^{\prime }(r_{\mathrm{out}})$$

$$E_{-}^{\prime }(r_{\mathrm{out}})$$

$$E_{+}^{\prime }(r_{\mathrm{in}})$$

$$E_{-}^{\prime }(r_{\mathrm{in}})$$

$$e(r_{\mathrm{out}})$$

$$E_{+}(r_{\mathrm{out}})$$

$$E_{-}(r_{\mathrm{out}})$$

$$\frac{d\omega }{dt}=0,$$

$$\frac{d\varpi }{dt}=\frac{d\Omega }{dt},$$

$$2n^2-{\nu }^2-{\kappa }^2=0,$$

$$\frac{d\varpi }{dt}\simeq -\frac{d\Omega }{dt}\simeq \frac{3G{M}_2}{4n{a}_{\mathrm{b}}^3},$$

$$\frac{d\varpi }{dt}$$

Full text

††pubyear: 2015††pagerange: Kozai–Lidov Disc Instability–Kozai–Lidov Disc Instability

1 Introduction

Discs in binary systems are sometimes misaligned with respect to their binary orbital planes. Misalignment is expected to be more likely in wider binaries with separations greater than $\sim 100$ AU where the effects of tidal dissipation in the disc are weaker and act on longer timescales (Papaloizou & Terquem, 1995; Bate et al., 2000; Lubow & Ogilvie, 2000; King et al., 2013). There are several ideas on how noncoplanarity could come about in a young binary star system. Noncoplanarity could be the result of initial conditions. For example, if a young binary star system accretes turbulent gas, a second generation of accreted gas is likely to be misaligned with the binary orbit and result in misaligned discs around the young stars (Bate et al., 2010). Alternatively, a coplanar disc may evolve to a noncoplanar state due to an instability, such as radiation warping (Pringle, 1996; Ogilvie & Dubus, 2001).

Observational evidence for disc misalignment has been found in several binaries. In binary stars separated by greater than 40 AU, misaligned discs 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 direct evidence comes from images of discs. The young binary system HK Tau with binary separation $\sim 400$ AU provides direct evidence for noncoplanarity, since discs are observed around both components with one disc edge-on and the other more face-on (Stapelfeldt et al., 1998). Recent ALMA observations by Jensen & Akeson (2014) suggest that its discs are mutually misaligned by more than $60^{\circ}$, although the plane of the binary orbit is not known. Strong mutual disc misalignment ($\sim 72^{\circ}$) has lately also been detected for the two circumstellar discs in V2434 Ori, a binary system in Orion with a similar binary separation (Williams et al., 2014).

Test particles that reside in orbits that are inclined to the orbital plane of a circular binary star can undergo the powerful effects of Kozai–Lidov (KL) oscillations (Kozai, 1962; Lidov, 1962). These oscillations cause particles to exchange inclination for eccentricity. The particle initially gains eccentricity while reducing its orbital tilt and later loses eccentricity while gaining orbital tilt back to its original value. The process repeats in a periodic manner. For KL oscillations to occur, the particle orbit must be misaligned by more than about $39.2^{\circ}$ and less than $141.8^{\circ}$ with respect to the binary orbital plane. The amount of eccentricity gain can be quite large. For example, an initially circular orbit at an inclination of $60^{\circ}$ achieves a maximum eccentricity of about 0.75 during a KL cycle. Further extensions of this theory show that even stronger effects can occur for eccentric orbit binaries (Ford et al., 2000; Lithwick & Naoz, 2011; Naoz et al., 2013; Teyssandier et al., 2013; Li et al., 2014; Liu et al., 2015). The KL effect for ballistic objects has been applied to a wide range of astronomical processes. These range from inclinations of asteroids and irregular satellites (Kozai, 1962; Nesvorný et al., 2003) to tidal disruption events (Chen et al., 2011), and the formation of Type Ia supernovae (Kushnir et al., 2013).

In a recent paper, Martin et al. (2014) found in SPH simulations that a fluid disc can undergo large-scale, coherent, damped KL oscillations. In a subsequent paper Fu et al. (2015a) found that these oscillations can occur over a broad range of disc and binary parameters. Disc self-gravity can suppress KL oscillations (Fu et al., 2015b), although the suppression typically requires that the disc conditions are not far from gravitational instability. Such oscillations would have important consequences on the disc tilt and eccentricity evolution. They may even result in the formation of planets, as a consequence of the strong compression of the disc gas due to the strong shocks produced in a highly eccentric disc (Fu et al., 2016).

The geometry of a KL disc is somewhat complicated. The disc undergoes nodal and apsidal precession, along with tilt and eccentricity oscillations. Such a configuration is challenging to simulate with Eulerian grid-based codes because adequate resolution would require that there be many grid cells, most of which are empty at any instant in time. KL discs are more easily studied with a Lagrangian code, such as SPH. However, the SPH simulations have certain limitations. First they are subject to effects of artificial viscosity. These effects can be partially controlled by including many particles, typically more than $1\times 10^{6}$. However, simulations with such a large number of particles require long running times. During the course of the simulations, viscosity causes the shape of the density distribution to evolve and the particle count to decrease, thereby reducing the resolution. Ideally, one would like to understand how the evolution would occur under conditions of low viscosity, in order to isolate the effects of viscosity from pressure. Another issue is the limited dynamic range of the discs studied in the SPH simulations. Typically, the ratio of the outer to inner disc radii is about a factor of 10. For a disc in a several hundred AU binary, this ratio is orders of magnitude larger.

One approach to overcoming the limitations of multidimensional simulations is to apply 1D reduced equations based on asymptotic methods that assume that the disc is thin: $H/r\ll 1$. Such equations have been developed in Ogilvie (2001) that describes the nonlinear evolution of eccentric discs and Ogilvie (1999, 2006) that describe the nonlinear evolution of warped discs. Unfortunately, there are not currently corresponding equations that describe the nonlinear evolution of both warped and eccentric discs as can occur in the case of KL discs.

Instead, we analyze in this paper the onset of the KL oscillations when the eccentricity is small and nonlinear effects can be ignored. The early states of a KL particle oscillation can be understood as an exponential growth of eccentricity resulting from a linear instability of a tilted circular orbit (Tremaine & Yavetz, 2014). It is the generalization of this instability to a continuous disc that we describe in this paper. The analysis is carried out as a linear stability problem of an inviscid disc. We are interested in determining the conditions under which KL oscillations can occur. The initial stages of KL oscillations are dominated by eccentricity changes over very small changes in disc tilt. Consequently, only the eccentricity evolution needs to be followed. In addition, we assume the disc remains flat (unwarped). The discs found in the SPH simulations are typically quite flat, although some warping by amounts of order the disc aspect ratio are found during the course of KL oscillations (see Fig. 5 of Fu et al., 2015a). The flatness can be maintained by bending waves communicated by pressure in the disc. The condition for flatness is that the sound crossing time is shorter than the warping timescale, which is of order the precession timescale at the outer edge of the disc (Papaloizou & Terquem, 1995; Larwood & Papaloizou, 1997; Lubow & Ogilvie, 2000) .

To carry out the analysis, we apply a linear 1D eccentricity evolution equation of Teyssandier & Ogilvie (2016) and augment that equation with terms that describe the effects of KL oscillations. In a recent preprint, Zanazzi & Lai (2016) took a very similar approach. Our results are in agreement with theirs where they overlap. Our analysis is different in that we map out the regimes of instability and analyze the consequences of having a small inner disc edge radius.

The outline of the paper is as follows. In Section 2, we derive the equations for the eccentricity evolution of a KL disc. Section 3 describes the parameters of the disc models that we analyze. Section 4 describes some results for the time evolution of eccentricity and the methods used in this paper. Section 5 discusses the role of resonances in determining KL instability. Section 6 describes the states of marginal stability for the disc models. Section 7 discusses how these states of marginal stability are related to zones of instability. Section 8 describes some properties of the eccentric modes, including the behaviour at small radii, the sensitivity of the results to the inner boundary location, and an estimate for the scaling of the level of nonlinearity with eccentricity. Section 9 contains a summary.

2 Derivation of Linearized KL Equation

2.1 Particle Evolution Equation for Small Eccentricity

We derive the linear evolution equations for a nearly circular fluid disc by first considering the standard secular KL particle evolution equations in the quadrupole approximation (e.g., Kiseleva et al., 1998). Consider a circular orbit binary with orbital radius $a_{\rm b}$. The particle orbits about a binary member with mass $M_{1}$ and is gravitationally perturbed by a companion member with mass $M_{2}$ on a circular orbit. To lowest order in particle eccentricity $e$, the secular evolution equations are given by

[TABLE]

where $a,e,i,\omega,\Omega$ and $n$ describe the semi-major axis, eccentricity, inclination, argument of periapsis, longitude of ascending node, and orbital frequency of a particle, respectively. In obtaining these evolution equations, we dropped terms of order $e^{2}$ or higher on the RHSs.

We apply a complex eccentricity defined such that $E=e\exp{({\rm i}\,\omega)}$. Consider a Cartesian coordinate system $(x,y)$ that lies in the plane of the orbit such that the positive $x-$axis lies along the instantaneous line of ascending nodes. Complex eccentricity $E$ is related to the eccentricity vector $\mn@boldsymbol{e}$ that points from the origin to periastron by $E=e_{x}+{\rm i}\,e_{y}$, so that $e_{x}(r,t)=Re(E(r,t))$ and $e_{y}(r,t)=Im(E(r,t)).$ In terms of a complex eccentricity, Equations (2) and (4) can be expressed as

[TABLE]

Note that the linear system $dE/dt={\rm i}(AE+BE^{*})$ with real coefficients $A$ and $B$ implies that $d^{2}E/dt^{2}=(B^{2}-A^{2})E$. Complex eccentricity $E$ is unstable if $B^{2}>A^{2}$. An initially nearly circular orbit has

[TABLE]

implying that an instability occurs with growth rate

[TABLE]

when $\cos^{2}{i}<3/5$, in agreement with Tremaine & Yavetz (2014).

2.2 Disc Evolution Equation for Small Eccentricity

We apply the above equations to a system of particles, or to a continuous eccentric disc, in which additional, internal forces of constraint in the normal direction maintain a rigid tilt. The particles then represent disc fluid elements. By the assumption that the disc remains flat, fluid elements then share a common $i$ and $\Omega$, although they have independently varying $a$, $e$, and $\omega$. By Newton’s Third Law, the sum of the internal forces over all particles (or rings) vanishes and we have

[TABLE]

where $h=na^{2}$ is the specific angular momentum to lowest order in $e\ll 1$ and $dm$ is a mass element of the disc. The integrals are taken over the entire disc. We can regard the disc inclination as constant in space and time during this small-$e$ phase.

Equation (A22) of Teyssandier & Ogilvie (2016) provides a 1D linear secular eccentricity evolution equation of a 3D locally isothermal fluid disc with local sound speed $c_{\rm s}(r)$ subject to gravitational interactions with a planet. It includes the effects of pressure and gravity perpendicular to the disc plane that vary along the eccentric orbits. The theory includes the effects of vertical oscillations of the eccentric disc resulting from a lack of vertical hydrostatic equilibrium that in turn leads to a prograde contribution to the precession of the disc (Ogilvie, 2008). We drop the planet-disc potential terms $\Phi_{\rm pd}$ and add the KL terms given by Equation (6), together with $d\Omega/dt$ given by Equation (11).

We consider a disc that extends from inner radius $r_{\rm in}$ to outer radius $r_{\rm out}$ and has surface density $\Sigma(r)$. The evolution equation is given by

[TABLE]

Quantity $J(r)$ is the disc angular momentum per unit radius divided by $\pi$ and is given by

[TABLE]

where $n(r)$ is the Keplerian orbital frequency about mass $M_{1}$ given by $n(r)=\sqrt{GM_{1}/r^{3}}.$ Functions $b_{i}(r)$ are given by

[TABLE]

where $c_{\rm s}(r)$ is the local disc sound speed, $\Sigma(r)$ is the disc surface density distribution, and $A(r)$ and $B(r)$ are gravitational terms due to the binary for small $e$ that are associated with Kozai–Lidov oscillations in an initially slightly eccentric disc and are given by

[TABLE]

with

[TABLE]

where the integrals in Equation (19) extend over $r$ from $r_{\rm in}$ to $r_{\rm out}$. $A(r)$ is the binary-induced apsidal precession rate relative to the line of ascending nodes, $\partial_{t}\,\omega(r,t)$ as would occur for a ballistic particle.

We adopt the inner and outer boundary conditions that

[TABLE]

For a Keplerian disc, the divergence of the velocity is proportional to $\partial_{r}E$. Consequently, this boundary condition is equivalent to requiring that the Lagrangian pressure perturbations near the disc edges vanish.

We consider normal modal solutions to Equation (12). The KL effect introduces a nonanalytic term involving $E^{*}$. We follow the method of solution given by equations (37)–(40) of Lubow & Ogilvie (2000) that involves a similar nonanalytic term but describes tilt, rather than eccentricity. We consider modes of the form

[TABLE]

where $\lambda$ is the complex growth rate. The modal equations are given by

[TABLE]

and

[TABLE]

where

[TABLE]

where a prime denotes differentiation by $r$, $T=c_{\rm s}^{2}$ is proportional to the temperature, and $A$ and $B$ are given by Equations (18) and (18). We apply the following boundary conditions:

[TABLE]

The latter condition imposes a normalization constraint on $e$ that is arbitrary, since the equations are linear. The values of $e(r)$ are then measured relative to $e(r_{\rm out})$.

3 Description of Models

We explore the properties of KL oscillations in a set of models all of which involve an equal mass binary and have outer disc edge $r_{\rm out}$ set to $0.3a_{\rm b}$. This disc outer radius is slightly larger than the value of about $0.25a_{\rm b}$ for a disc in an equal-mass binary assuming coplanar orbits (Paczyński, 1977). A misaligned disc feels a weaker binary torque and thus the outer truncation radius is larger (e.g. Lubow, Martin & Nixon, 2015; Nixon & Lubow, 2015; Miranda & Lai, 2015) and so we adopt a larger value. The other parameters are listed in Table 1, where $p$ and $s$ are defined by $\Sigma(r)\propto r^{-p}$ and $T(r)=c^{2}_{\rm s}\propto r^{-s}$. Models A1 and A2 involve flared discs that have a disc aspect ratios $H/r\propto r^{1/8}$. This case applies to a standard ‘active’ accretion disc (Lynden-Bell & Pringle, 1974). An active protostellar disc is typically dominated by accretional heating within about $1$ AU from the central star, where it would follow such a temperature profile. Models B1 and B2 involve discs with somewhat greater flaring, having disc aspect ratios $H/r\propto r^{1/4}$. Such a level of flaring (or more) is expected in so-called ‘passive’ protostellar discs, where the heating is dominated by the contributions from the central star at all radii, or in active discs on scales greater than $\sim 1$ AU (Chiang & Goldreich, 1997). We explore a range of disc aspect ratios at the disc outer edge $h_{\rm out}=H/r(r_{\rm out})$, typically between 0.02 and 0.15.

Observations of protostellar discs suggest that the surface density parameter power law exponent $p$ is in the range of 0.5 to about 1, although there is considerable uncertainty (Willams & Cieza, 2011). We adopt a value of 1. Since we are interested in fairly wide binaries, we apply two values for the disc inner radius that are small compared to the binary separation. The orbital radii $\sim 10R_{\odot}$ for the so-called hot Jupiter planets suggests a crude estimate for the central hole size in the disc where migrating planets are trapped (Lin, Bodenheimer, & Richardson, 1996). For a binary with separation $\sim 500$ AU and period $\sim 10^{4}$ yr, and a disc with a $\sim 10$ solar radii inner hole, the inner disc radius $r_{\rm in}\sim 10^{-4}a_{\rm b}$. Consequently, it is important to analyze the eccentricity behavior at such small relative radii. We consider two values of inner radii with other parameters fixed in order to probe the sensitivity of the results to $r_{\rm in}$.

4 Time Evolution and Methods of Solution

We determine the eccentricity evolution by applying similar methods used in Lubow (2010) with Mathematica. We compute the complex eccentricity evolution given by Equation (12) in space and time for model A1 with the initial condition that $E(r,0)=1$ for $r_{\rm in}\leq r\leq r_{\rm out}$ using the method of lines. The initial eccentricity satisfies boundary conditions (20). The results for this model with an initial tilt $i=45^{\circ}$ and disc aspect ratios at the disc outer edge of $h_{\rm out}=0.03$ and 0.05 are plotted in Figure 1. After less than 20 binary orbital periods $P_{\rm b}$, $E(r,t)/|E(r_{\rm out},t)|$ settles closely to an eigenfunction in which $\partial_{t}E/E$ is nearly constant in $r$ and is equal to the eigenvalue. As we will see later, there is more than one growing eigenmode in each case. The plotted distributions at $t=50P_{\rm b}$ reflect the fastest growing eigenmode that emerges as dominant for each case. The eccentricity distributions are quite different for the two cases. The final eccentricity distribution has a lower eccentricity at the inner edge than the outer edge for the $h_{\rm out}=0.03$ case, while just the opposite occurs for the $h_{\rm out}=0.05$ case.

We numerically determine eigenfunction $E(r)$ and eigenvalue $\lambda$ in Equations (21)–(23) by a shooting method. Equations (22) and (23) are second order in space and are integrated inward in radius from $r=r_{\rm out}$ to the inner boundary $r=r_{\rm in}$ by using the starting conditions for the radial derivatives given by Equation (28) and (29) at the outer radius $r_{\rm out}$. The normalization condition (32) is replaced by

[TABLE]

for the real weighting factor $w$ and phase $\omega_{\rm out}$ at the disc outer edge. The complex inner boundary conditions (30) and (31) are equivalent to four real conditions. These are satisfied by adjusting the four real parameters $Re(\lambda)$, $Im(\lambda)$, $\omega_{\rm out}$, and $w$. Once an eigenmode is determined for a model, parameter changes from that eigenmode, such as in the disc inclination $i$, can then be made incrementally and iteratively with rapid convergence. Consequently, it is more efficient to solve the eigenvalue problem, rather than the initial value problem.

Although we do not explicitly consider cases with $i>90^{\circ}$, the growth rates are invariant when $i$ is replaced by $180^{\circ}-i$.

5 Role of Resonances

The KL instability can be understood as a consequence of a resonance in which

[TABLE]

or equivalently the matching of apsidal and nodal precession rates

[TABLE]

where $\varpi$ is the longitude of the periapsis. A test particle is subject to perturbations that are solely due to the gravitational forces of the companion. From Laplace’s equation, it follows that near the binary orbital plane

[TABLE]

where $\nu$ is the vertical oscillation frequency and $\kappa$ is the epicyclic oscillation frequency of a test particle. Using this identity and that $\nu\simeq\kappa\simeq n$, it follows that

[TABLE]

for $r\ll a_{\rm b}$, where $d\varpi/dt=n-\kappa$ and $d\Omega/dt=n-\nu.$ Consequently, resonance condition (36) is impossible to satisfy for a nearly coplanar test particle orbit where $i\simeq 0$.

In the case of a test particle, we have from Equation (4) that

[TABLE]

Notice that for small $i$, these equations recover Equation (38). With increasing inclination $i\leq 90^{\circ}$, the apsidal precession rate given by Equation (39) can become negative, while the nodal precession rate in Equation (40) remains negative. Resonance is first possible when $\omega=90^{\circ}$ ($\omega$ is constant in time by Equation (35)), in order to make the apsidal rate as negative as possible. Applying this $\omega$ value and Equations (39) and (40) to resonance condition (36), we then recover the KL instability condition for a test particle given at the end of Section 2.1 that $\cos^{2}{i}<3/5$. We see that KL instability is tied to satisfying the resonance condition.

The effects of pressure on a fluid disc modify the apsidal precession rate and can act to make it negative, even for a coplanar disc. Consequently, as we will see, the coplanar KL resonance condition can be satisfied for certain disc parameters. There are several (formally infinitely many) values of disc aspect ratios for which the coplanar KL condition is satisfied. These coplanar resonances play a key role in delineating the various branches of instability, as we will see in the next section.

In this paper we are concerned with KL instability. Consequently we describe instability that is associated with the condition that $d\omega/dt=0$ or

[TABLE]

in Equation (21). Other forms of eccentric disc instability that do not satisfy this condition are possible (e.g., Lubow, 1991).

6 Marginal Stability

We determine the conditions for marginal KL stability of the disc models described in Section 3. Marginally stable KL modes have $Re(\lambda)=0^{+}$ and $Im(\lambda)=0$ by Equation (41). These modes are stable, but are very close to instability. We determine some properties of the marginally stable modes by using Equation (12) with the RHS set to zero, since $\partial_{t}E=0$. Expressing $E=e_{\rm x}+ie_{\rm y}$ for real $e_{\rm x}$ and $e_{\rm y}$ and taking the real and imaginary parts of that equation, we find that there are modes that satisfy

[TABLE]

with normalization condition

[TABLE]

and other modes that satisfy

[TABLE]

with normalization condition

[TABLE]

The boundary condition given by Equation (20) implies that

[TABLE]

for $j={\rm x},{\rm y}$.

The marginal modes are therefore untwisted and have a constant phase of either $\omega=(0^{\circ},180^{\circ})$ or $\omega=(90^{\circ},270^{\circ})$. We note that the same phase results occur for 2D discs and adiabatic discs using appropriate eccentricity evolution equations given by Teyssandier & Ogilvie (2016). Equations (42)–(44) or Equations (45)–(47), together with Equation (48) are regarded as an eigenvalue problem for the critical tilt angle for the marginally stable mode with $i=i_{\rm crit}$.

We consider all parameters fixed, except the disc aspect ratio at the disc outer edge $h_{\rm out}=H/r(r_{\rm out})$ and the inclination $i$. We determine the critical angle $i_{\rm crit}$ as a function of $h_{\rm out}$ for all marginal modes. To determine the critical angles for modes with $\omega=(0^{\circ},180^{\circ})$, we start with a given value of $h_{\rm out}$ and integrate Equation (42) inward with outer boundary condition (48) and normalization condition (44). We determine the values of $i=i_{\rm crit}$ such that the inner boundary condition (48) is satisfied. That is, by integrating inwards for several assumed $i$ values, we determine the values of $de_{\rm x}/dr(r_{\rm in})$. We find the $i$ values near where $de_{\rm x}/dr(r_{\rm in})$ changes sign and refine the search for $i_{\rm crit}$ with a local root finder. With this technique, we believe that we determine all marginal modes at a given $h_{\rm out}$. The same process is applied to modes with $\omega=(90^{\circ},270^{\circ})$.

In general there exists more than one $i_{\rm crit}$ for the same $h_{\rm out}$ corresponding to different modes with a different number of radial nodes in the eigenfunctions $e(r)$. We follow these marginal modes to higher and lower values of $h_{\rm out}$ and obtain $i_{\rm crit}(h_{\rm out})$.

Figure 2 plots the inclinations $i_{\rm crit}$ for the five marginal modes with the highest $h_{\rm out}$ values for models A2 and B2 over $0.02\la h_{\rm out}\la 0.15$. The branch numbers are not directly related across models. We later determine this relationship. For both models, branches 1, 2, and 3 have phase $\omega_{\rm out}=90^{\circ}$, while the nearly vertical branches 4 and 5 have phase $\omega_{\rm out}=0$ and cross the other branches. The values of $i_{\rm crit}(h_{\rm out})$ along branch 4 are $\sim 3\,i_{\rm crit}(h_{\rm out})$ along branch 2, until branch 4 reaches $i=90^{\circ}$. The same relationship exists between branches 3 and 5. Notice that all branches except for branch 1 extend to coplanarity, $i=0$. For larger $h_{\rm out}$, the marginal stability of branch 1 extends to $90^{\circ}$. More branches occur at smaller $h_{\rm out}$ than are plotted here. There is an another marginally stable branch for model B2 with $i=90^{\circ}$ at $h_{\rm out}\simeq 0.025$ that resembles branch 1 that we do not include in the plot. Notice that the number of branches with phase $\omega_{\rm out}=90^{\circ}$ (branches 1, 2, and 3) at the same $h_{\rm out}$ value increases with decreasing $h_{\rm out}$. The marginal mode branches for models A2 and B2 plotted in Figure 2 have a very similar topology, but are shifted somewhat relative to each other.

Figure 3 plots the eigenfunctions for branches 1, 2, and 3 of models A2 and B2. The vertical scale is logarithmic and the dips in the plots correspond to nodes. The number of nodes in the region of the disc with $r\ga 0.02a_{\rm b}$ is seen to increase with the branch number. The modes along a given branch of a model have the same number of nodes. The behaviour at small $r$ is hard to discern in this plot, but will be described in Section 8. The eigenfunctions for branches 4 and 5 are very similar to those for branches 2 and 3, respectively, but are shifted in phase.

In Figure 4 we plot the $h_{\rm out}$ values for some of the coplanar KL resonances as a function of the temperature exponent $s$ for a fixed density exponent $p=1$. Some points are labeled by model and branch number. For example, $h_{\rm out}=H/r(r_{\rm out})\simeq 0.07$ for the point on the plot that is labeled for branch 2 of model A2 corresponds to the $h_{\rm out}$ value at zero inclination in the top panel of Figure 2 for branch 2. The effect of increasing $s$ is to shift the resonances towards lower $h_{\rm out}$. The direction of this shift is in the sense of reducing changes to the pressure gradient, since an increase in $s$ raises the pressure gradient at fixed sound speed, while the decrease in $h_{\rm out}$ lowers the sound speed. We see from the plot that branch 3 of model A2 lies on the same curve as branch 2 of model B2. Consequently, they correspond to the same mode structure, i.e., have same number of nodes in the eigenfunction. We then expect that the branch numbers for the same mode structure differ by unity across the two models. Therefore, branch 2 of model A2 corresponds to branch 1 of model B2. The mode structure of branch 1 of model A2 is not present in model B2. The mode structure for branch 3 of model B2 exists at smaller $h_{\rm out}$ than is plotted for model A2 in the top panel of Figure 2. In Figure 5 we plot $h_{\rm out}$ for some of the coplanar KL resonances as a function of the density exponent $p$ for a fixed temperature exponent $s=0.75$. The effect of increasing the exponent is again to shift the resonances towards lower $h_{\rm out}$.

7 Zones of Instability

We apply the marginally stable state information to determine the parameter ranges for KL instability (requiring $Im(\lambda)=0$). We expect that modes will be growing in time for inclinations $i$ above the marginal stability curves in Figure 2.

We consider continuous sequences of modes (eigenfunctions and eigenvalues) for some model that start on marginally stable branch 1, 2, or 3 and extend to larger $i$ at fixed $h_{\rm out}$ (vertically upward from $i_{\rm crit}(h_{\rm out})$ on Figure 2). We consider these modes as functions of parameters, that we denote as $M(b,h_{\rm out},i)$, where $b$ is the starting branch number 1, 2, or 3 of the mode sequence. From brevity, we sometimes refer to such a mode as being on an extension of a particular branch. For example, a mode that is obtained from a continuous sequence of modes at fixed $h_{\rm out}$ that start on branch 1 is referred to as lying on an extension of branch 1. For a given $h_{\rm out}$ and inclination $i$, any branches $b$ with $i>i_{\rm crit}$ may provide unstable modes that lie in their extensions. In general, there can be multiple unstable modes that are extensions of different branches.

Figure 6 plots the growth rates and phases as a function of inclination for modes that lie on extensions of branch 2 for model A2, that is $M(2,h_{\rm out},i)$. These cases have $h_{\rm out}$ values of 0.06607, 0.06616, 0.06812, and 0.06941. The peak growth rates increase with starting angle and decrease with $h_{\rm out}$. For the bottom three rows in the figure, the plots begin with phase $\omega_{\rm out}=90^{\circ}$ and end with $\omega_{\rm out}=0^{\circ}$. Recall that modes on branches 2 and 4 have constant $\omega$ with $\omega_{\rm out}=90^{\circ}$ and $\omega_{\rm out}=0^{\circ}$, respectively. Instability terminates with increasing $i$ on branch 4 with $\omega_{\rm out}=0^{\circ}$ at approximately 3 times the starting inclination angle on branch 2 for the lower three cases. The absolute values of the derivatives at the endpoints of the curves in the left and right panels increase with starting angle and reach a near infinite absolute value at $i\simeq 43^{\circ}$ on the second panel from the top. For a slightly smaller aspect ratio, $h_{\rm out}=0.06607$ there is an abrupt change of behaviour in which the terminating angle is no longer limited by branch 4 and is plotted to $90^{\circ}$.

The modes along a marginally stable branch have the same mode structure, i.e., the same number of nodes. However, modes that extend from a marginal mode branch to higher inclinations can undergo major structural changes. Figure 7 plots the eigenfunctions for the model in the top row of Figure 6 for three different tilt angles. These modes lie on extensions of branch 2 of model A2 that correspond to $M(2,0.06607,12.5^{\circ})$, $M(2,0.06607,30^{\circ})$, and $M(2,0.06607,60^{\circ})$ in the notation defined above. The mode at $i=12.5^{\circ}$ is marginally stable and lies on branch 2. Its eigenfunction resembles the branch 2 eigenfunction for a different $h_{\rm out}$ that is shown in the middle left panel of Figure 3. However, at larger tilt angles, the outer node disappears and the eigenfunction more closely resembles a branch 1 marginal mode of model A2, such as the one shown in the upper left of Figure 3. Notice also that the knee in the growth rate curve shown in Figure 6 occurs at about $i\simeq 45^{\circ}$ that is roughly at the critical angle for branch 1, $i_{\rm crit}(0.06607)\simeq 48^{\circ}$. Above this knee, the growth rate climbs substantially and the eigenfunctions have a simpler structure that are similar to a branch 1 mode.

The abrupt change in behaviour across the top two rows of Figure 6 in going from $h_{\rm out}=0.06607$ to 0.06616 is likely related to the crossing of branch 1 with branch 4 in Figure 2. The models plotted in the lower three rows of Figure 6 have a positive growth rate over a much more limited range of angles than the top row. The reason is that these modes are terminated at angles where they intersect with branch 4 at fixed $h_{\rm out}$ in proceeding vertically upward from branch 2 to branch 4 in Figure 2.

Figure 8 plots the growth rates for some sequences of modes involving model A2. For $h_{\rm out}=0.05,0.06$ and $0.065$, the sequences begin on branch 2, while for $h_{\rm out}=0.07,0.08,0.09$ and 0.1 the sequences begin on branch 1. (In the case $h_{\rm out}=0.065$ there is a short second sequence of unstable modes between branches 4 and 1, near $i=50^{\circ}$ that we do not plot in this figure.) In most cases, the growth rate increases monotonically with tilt angle above the critical tilt for marginal stability. The case of $h_{\rm out}=0.07$, which is near the $h_{\rm out}$ value for the crossing of branches 1 and 4, is somewhat different. The plot in that case shows that unstable modes first extend below branch 1 before rising to higher inclinations with higher growth rates. The sequence joins to another sequence of unstable modes with complex growth rates (indicated by the dashed line) that extends to lower inclination. This sequence terminates near $i=40^{\circ}$ at a marginally stable mode that has a nonzero precession frequency. A similar behaviour occurs for $h_{\rm out}=0.08$.

A closeup of the region near branch 4 is shown in Figure 9. Branch 4 is shown as the dotted line that overlaps with the solid line on the lower right. The solid line for branch 4 is an upper boundary of the region of instability above branch 2 from $h_{\rm out}\simeq 0.0665$ to 0.0696. The region at low inclinations below branch 2 is stable. Branch 1 is plotted by a dashed line on the left and the solid line on the right of the dashed short vertical line. In addition, the region between the solid lines for branches 1 and 4 also contains unstable modes. Over the range of the solid line for branch 1 in this figure, sequences of unstable modes extend below branch 1, as discussed above for $h_{\rm out}=0.07$. Such sequences include a complex extension that is terminated by a marginal mode that we have not plotted because it has a non-zero precession frequency. We do not pursue a detailed analysis of the stability of modes in this region in this paper. The triangular region in the figure is an unstable zone below branch 1 that also includes a complex extension that is terminated by a marginal mode that has $Im(\lambda)\neq 0$. It also includes the some of unstable modes from the mode sequence plotted top row in Figure 6.

Branch 4 extends from $h_{\rm out}=0.0614$ at $i=90^{\circ}$ to $h_{\rm out}=0.0696$ at $i=0^{\circ}$. The top row of Figure 6 is for $h_{\rm out}=0.06607$ that lies within the span of $h_{\rm out}$ values of branch 4. Yet, the range of angles over which instability occurs is not terminated by branch 4. The abrupt change in behaviour across the top two rows of Figure 6 in going from $h_{\rm out}=0.06607$ to 0.06615 is likely related to the crossing of branch 1 with branch 4 in Figures 2 and 9.

As branch 4 reaches close to branch 1 from lower $i$, its influence on modes that are an extension of branch 2, that is $M(2,h_{\rm out},i)$, abruptly stops and shifts to modes that are on an extension of branch 1, $M(1,h_{\rm out},i)$. This effect on modes that are an extension of branch 1 is seen in Figure 10. Notice that the smallest value $h_{\rm out}=0.055$ that is plotted (top row in this figure) is smaller than of the minimum $h_{\rm out}$ value spanned by branch 4, while the largest value (top row in this figure) has $h_{\rm out}=0.10$ that lies beyond the largest $h_{\rm out}$ value spanned by branch 4. In both cases, the positive growth rates extend from the critical value of $i_{\rm crit}$ at marginal stability plotted in Figure 2 to $90^{\circ}$. However, the intermediate $h_{\rm out}=0.065$ lies within the span of branch 4 as seen in Figure 9. Its positive growth rates are terminated by branch 4 where its phase reaches $180^{\circ}$, in a somewhat similar manner as occurs for modes that are a continuation of branch 2, as seen in the bottom three panels of Figure 6.

Branches 1 and 4 cross at $h_{\rm out}=0.0654$ corresponding to $i=48.4^{\circ}$. The second row of Figure 6 shows the growth terminating by a crossing with branch 4 at $i=42.9^{\circ}$. Consequently, the boundary for the interaction between modes $M(2,h_{\rm out},i)$ and branch 4 occurs at $h_{\rm out}$ values that are slightly larger than where branches 1 and 4 cross, where the dashed short vertical line lies.

Away from mode crossings, we find that the dominant (fastest growing) mode is the one that extends from the branch having the smallest critical angle $i_{\rm crit}$. For example, for $h_{\rm out}\simeq 0.035$ to 0.07, modes $M(2,h_{\rm out},i)$ dominate over modes $M(1,h_{\rm out},i)$ for the same $h_{\rm out}$ and $i$, except possibly for the regions close to branch 4. This can be seen for example in Figure 11. The plotted growth rates correspond to the same $h_{\rm out}=0.055$ value as for the upper left panel of Figure 10. Positive growth rates for the mode $M(1,0.055,i)$ occur for $i\ga 46^{\circ}$. But at such angles the growth rates plotted in Figure 11 for the modes $M(2,0.055,i)$ are larger.

In Figure 2 we see that the highest $h_{\rm out}$ for a marginal mode occurs for branch 1 that reaches an inclination of $90^{\circ}$ for $h_{\rm out}\sim 0.15$. As a result, the range of unstable angles vanishes for such disc aspect ratios and above. Additional apsidal precession that is not due to the gravitational forces of the binary companion can suppress KL oscillations if that additional apsidal precession rate exceeds that due to the binary (e.g., Fabrycky & Tremaine, 2007). As noted in Martin et al. (2014), the pressure-induced apsidal precession rate is of order $\sim h_{\rm out}^{2}n$. For the pressure apsidal precession rate to exceed the binary precession rate $\sim q\,n_{\rm b}^{2}/n$ requires roughly $h_{\rm out}\ga\sqrt{q}\,n_{\rm b}/n$, where $q$ is the binary mass ratio $M_{2}/M_{1}$ that is assumed to be in the range $q\la 1$, as is similar to the findings of Zanazzi & Lai (2016). In the current case with $q=1$, we have that $n_{\rm b}/n\sim 0.2$. So this relationship is roughly satisfied. As was also noted in Martin et al. (2014), another condition for KL oscillations is that the inverse radial sound crossing time be shorter than the disc precession rate in order that the disc remains flat (Larwood & Papaloizou, 1997). If this condition is not satisfied, the disc will severely warp and may break up. This condition can crudely be written as $n\,h_{\rm out}\ga qn_{\rm b}^{2}/n$. For the models in this paper this condition crudely requires that $h_{\rm out}\ga 0.05$. However, since we have dropped all factors of order unity, the flatness condition might be satisfied at somewhat smaller values as well, as was found in SPH studies (Martin et al., 2014; Fu et al., 2015a). In any case, we then have crudely that KL oscillations can occur because there is a window of $h_{\rm out}$ values for which

[TABLE]

This relationship can also be expressed as a constraint on binary mass ratios as

[TABLE]

where $q_{\rm crit}=(h_{\rm out}n/n_{\rm b})^{2}.$

For a disc that relies on only self-gravity to maintain flatness, the level of self-gravity required for flatness produces an apsidal precession rate that exceeds the precession rate due to the binary companion and therefore suppresses KL oscillations (Batygin et al., 2011). As a result, there is no equivalent KL window analogous to Equation (49) for such a disc.

8 Properties of KL Eccentric Modes

As discussed in the Section 3, the disc inner radius is likely to be about $\sim 10^{-4}a_{\rm b}$ for wider binaries where the KL oscillation timescales are of order $10^{4}$ yr or larger. We see from Figure 3 that the eccentricity at the inner boundary can be much larger than the eccentricity in the outer parts of a disc. We examine the behaviour of the eccentric modes at small radii and the sensitivity of the results to the inner radius.

For $r\ll a_{\rm b}$, the time dependence of $E$ can be ignored and the spatial dependence of $E$ can be determined. Equation (12) reduces to

[TABLE]

where $p$ and $s$ are the pressure and temperature exponents discussed in Section 3. We apply boundary condition

[TABLE]

Equation (51) has power law solutions in which the exponent is generally a complex number. For $\sigma$ real, terms of the form $r^{{\rm i}\,\sigma}$ can be expressed as $\exp{({\rm i}\,\sigma\log{r})}$. For the range of parameters in this paper, the solution can then be expressed in the form

[TABLE]

where $k$, $\phi$, and $z$ are real constants which depend on $p$ and $s$ and are determined analytically. There is also a normalization constant $\beta$ that we adjust, as described below. We find that

[TABLE]

For typical surface density profiles $p<2$ and so $z$ is positive. The eccentricity then varies at small $r$ as $e(r)\sim 1/r^{z}$, ignoring oscillatory terms, and is therefore divergent.

For model A, the analytic solution of the form given in Equation (53) has $k=\sqrt{31}/4$, $\phi=\arctan{(2/\sqrt{31}})$, and $z=1/2$. The only free parameter is the normalization constant $\beta$ that is chosen so that the analytic and numerically determined solutions match at $r=r_{\rm in}$. Figure 12 plots a typical case. The analytic result is barely distinguishable from the numerical solution for $r\la 0.03a_{\rm b}$.

In an adiabatic disc, there is a globally conserved quantity called the angular momentum deficit, which expresses the difference between the angular momentum of elliptical and circular orbits with the same energy (Goodchild & Ogilvie, 2006; Teyssandier & Ogilvie, 2016). In a locally isothermal disc, this quantity is not conserved due to the exchange of heat with the background radiation. However, there is a related globally conserved quantity that we refer to as the modified angular momentum deficit (MAMD) that is given in equation (C23) of Teyssandier & Ogilvie (2016) as

[TABLE]

where the integral extends over the disc. To measure the influence of the large amplitudes near the inner boundary, we determine the fractional contributions to ${\cal A}$ that arise in its vicinity.

From Equation (53), the MAMD per unit radius (the integrand in Equation (55)) can be shown to vary as $\sim r^{s-1/2}$ (ignoring oscillatory terms). Consequently, the integrated contributions of the inner region to the MAMD increase roughly as $\sim r^{s+1/2}$. This integral is then well behaved and nonsingular for small $r$ for positive values of $s$, even though $|E(r)|=e(r)$ may be singular as $r=r_{\rm in}$ goes to zero.

Figure 13 plots $e(r_{\rm in})/e(r_{\rm out})$ as a function of inclination $i$ for modes that extend from branch 2 of model A2 with $h_{\rm out}=0.055$. At smaller inclinations, where the growth rate is positive but smaller (see Figure 11), the eccentricity at the inner boundary is substantially bigger than at the outer boundary. There is a rapid drop in this ratio near $i\simeq 42^{\circ}$. In Figure 14 we plot the cumulative radial distribution starting at $r=r_{\rm in}$ of the MAMD. As expected from Figure 13, there is a rapid change in the distributions between $i=40^{\circ}$ and $45^{\circ}$. This distribution is normalized to unity at $r=r_{\rm out}$. However, the cumulative MAMD distributions show only small contributions coming from the region near the inner boundary. For both models the MAMD distributions are fairly similar at similar inclinations. The MAMD near the inner boundary for $i=13^{\circ}$ in the bottom panel reaches a few percent contribution out to $r\simeq 0.005a_{\rm b}$ and the contributions are smaller at larger inclinations.

We consider the sensitivity of the results to the inner boundary location by comparing model A1 to A2 and models B1 to B2. The models being compared have inner radii that differ by a factor of 100. Figure 15 plots the branches of marginal modes for models A1 and B1. In comparing this figure with Figure 2 we see that there are only small differences between models A1 and A2.

On the other hand, there are major differences between models B1 and B2. Model B1 has some unusual properties. Branch 1 (the branch that does not reach to $i=0$) ranges over lower $h_{\rm out}$ values in model B1 than in model B2. A major difference is that modes that lie on the extension of branch 1 in model B1 never dominate. The reason is that the span of $h_{\rm out}$ values covered by branch 1 is also covered by other branches with lower critical angles for marginal stability. For example, at $h_{\rm out}=0.05$ branches 1 and 2 are marginally stable. As shown in Figure 16, the modes that extend from branch 2 dominate at all inclinations in this case. The termination of positive growth with increasing $h_{\rm out}$ in model B1 is then different from the other models we have presented. Rather than growth being terminated at higher $h_{\rm out}$ by branch 1 at $i_{\rm crit}=90^{\circ}$, as in the other models, growth in model B1 appears to be terminated by branch 2 at $i_{\rm crit}=0^{\circ}$ that occurs for $h_{\rm out}\simeq 0.15$.

For a locally isothermal disc with a very small inner radius, we estimate the scaling of the eccentricity required for nonlinearity to set in. A relevant measure of nonlinearity is $|rdE/dr-E|$ (Ogilvie, 2001). Using Equation (54), we have that this measure of nonlinearity in the inner disc scales like $e(r)\sim e(r_{\rm out})r^{p/2-1}$. So for fixed $i$, the radius within which the mode becomes nonlinear scales as $r_{\rm NL}\sim e(r_{\rm out})^{1/(1-p/2)}$. The fraction of MAMD that is subject to nonlinear modification scales as $\sim r_{\rm NL}^{s+1/2}\sim e(r_{\rm out})^{(2s+1)/(2-p)}$. For models A and B this power is 2.5 and 2 respectively. Nonlinearity of the eccentric mode can lead to enhanced dissipation. When the eccentricity gradient is sufficiently large that neighbouring orbits approach mutual intersection (Ogilvie 2001), the local rate of viscous damping of the mode increases nonlinearly. At slightly lower amplitudes, inertial waves in the disc become strongly destablilized by the eccentric motion (Papaloizou, 2005a, b; Barker & Ogilvie, 2014), and this is also expected to damp the eccentricity in a way that depends nonlinearly on its amplitude. The occurrence of either of these mechanisms in the inner disc provides a possible way to saturate the growth of a global mode through nonlinear damping, regardless of the (very small) value of the inner radius.

9 Summary

We have investigated the linear stability of tilted discs in binary star systems to Kozai–Lidov (KL) oscillations by extending the model of Teyssandier & Ogilvie (2016) that includes 3D effects in a 1D calculation. The results support the existence of disc KL oscillations found in SPH simulations (Martin et al., 2014; Fu et al., 2015a, b). We find that sufficiently tilted discs can undergo KL oscillations in binaries with order unity mass ratios provided that roughly $(n_{\rm b}/n)^{2}\la h_{\rm out}\la n_{\rm b}/n$, for binary orbital frequency $n_{\rm b}$, disc aspect ratio $h_{\rm out}$, and orbital frequency $n$ at the disc outer edge. In agreement with Zanazzi & Lai (2016), we find KL disc instability is not possible for binaries with small mass perturbers in which the binary mass ratio $M_{2}/M_{1}\la(h_{\rm out}\,n/n_{\rm b})^{2}$, but that instability is possible in equal mass binaries for disc inclinations that lie below the critical value required for test particles of $39.2^{\circ}$.

As the disc evolves, its eccentricity follows an eigenmode (see Figure 1). The states of marginal stability lie on mode branches in a $h_{\rm out}-i_{\rm crit}$ diagram, as shown in Figure 2. Generally, discs are unstable if their tilts lie above the critical inclination for any marginal stability branch. There are cases where this is not true, as seen in the lower panel of Figure 15. Modes can be unstable in discs that lie very close to the disc midplane, unlike the case of test particles. However, the growth rates and range of unstable angles is small for nearly coplanar discs (see lower panels of Figure 6). The coplanar ($i=0$) apsidal-nodal disc resonances delineate the marginally stable branches. The disc aspect ratios $h_{\rm out}$ for these resonances are a function of the disc temperature and density gradients (see Figures 4 and 5). In general, more than one unstable mode can be present in a disc. The number of these unstable modes increases in cooler discs. Mode crossings limit the range of unstable tilt angles for modes that lie on extensions of branches whose critical inclinations $\la 12^{\circ}$, as seen in Figure 6. The dominant mode is the one that lies on an extension of a marginally stable branch having with the smallest critical angle (e.g., Figure 16).

Discs at smaller tilt angles tend to have larger eccentricities near the disc inner edge relative to the disc outer edge (Figure 13). The eccentricities are formally divergent as the inner radius goes to zero. However, their effect on the modified angular momentum deficit near the inner edge is small (Figure 14). At the end of Section 8, we describe the scaling of eccentricity for nonlinearity to set in.

We found some sensitivity of the results to the location of the inner boundary in models B1 and B2. The marginal stability curves, especially for branch 1 for model B1 in Figure 15, are shifted relative to model B2 in Figure 2.

We investigated some of the interactions between modes at mode crossings (see Figure 9). But more remains to be explored, particularly involving oscillatory unstable modes. We have idealized the disc as flat, but bending modes may be present that could provide additional mechanisms for instability.

Acknowledgements

We thank Rebecca Martin for useful discussions and informing us about the preprint by Zanazzi & Lai (2016). SHL acknowledges support from NASA grant NNX11AK61G.