# Instability of warped discs

Instability of warped discs Abstract Accretion discs are generally warped. If a warp in a disc is too large, the disc can ‘break’ apart into two or more distinct planes, with only tenuous connections between them. Further, if an initially planar disc is subject to a strong differential precession, then it can be torn apart into discrete annuli that precess effectively independently. In previous investigations, torque-balance formulae have been used to predict where and when the disc breaks into distinct parts. In this work, focusing on discs with Keplerian rotation and where the shearing motions driving the radial communication of the warp are damped locally by turbulence (the ‘diffusive’ regime), we investigate the stability of warped discs to determine the precise criterion for an isolated warped disc to break. We find and solve the dispersion relation, which, in general, yields three roots. We provide a comprehensive analysis of this viscous-warp instability and the emergent growth rates and their dependence on disc parameters. The physics of the instability can be understood as a combination of (1) a term that would generally encapsulate the classical Lightman–Eardley instability in planar discs (given by ∂(νΣ)/∂Σ < 0) but is here modified by the warp to include ∂(ν1|ψ|)/∂|ψ| < 0, and (2) a similar condition acting on the diffusion of the warp amplitude given in simplified form by ∂(ν2|ψ|)/∂|ψ| < 0. We discuss our findings in the context of discs with an imposed precession, and comment on the implications for different astrophysical systems. accretion, accretion discs, black hole physics, hydrodynamics, instabilities 1 INTRODUCTION Accretion discs occur throughout astrophysics, from the birthsites of stars and planets to the feeding of supermassive black holes (SMBHs) in galaxy centres. Through several different effects, each important in different contexts, these discs can warp. Discs can be formed in a warped configuration through, for example, chaotic infall on to protostellar discs (Bate, Lodato & Pringle 2010) or SMBHs (Lucas et al. 2013). Alternatively, discs can warp if they orbit in a non-spherical gravitational potential that induces radial differential precession, for example, the Lense–Thirring effect from a spinning black hole (Lense & Thirring 1918) or tides from a companion star (Papaloizou & Terquem 1995). Further, there are warping instabilities that act on planar discs to drive a disc warp from an initially axisymmetric configuration with only a small perturbation, for example, radiation warping (Pringle 1996, 1997), winds (Schandl & Meyer 1994), or resonant tidal effects (Lubow 1992; Lubow & Ogilvie 2000). Once a disc is warped, the local plane of the disc changes with radius. This creates a displacement of the high-pressure mid-plane region from one ring of the disc to its neighbour which oscillates with azimuth. At the descending (azimuthal angle ϕ = 0) and ascending (ϕ = π) nodes, where the neighbouring rings cross, there is no displacement, whereas it is maximal at ϕ = π/2 and (with opposite direction) at 3π/2 (see e.g. fig. 1 of Nixon & King 2016). Thus, a fluid element orbiting in the disc feels a radial pressure gradient that oscillates on the local orbital frequency. The fluid element then exhibits epicyclic motion. In a near-Keplerian disc, the epicyclic frequency is close to the orbital frequency, creating a resonance between the motion and the forcing. The epicyclic motion communicates the disc warp radially, and if undamped, it will lead to the launch of a bending wave. However, in the presence of turbulence, the wave will be damped, and if damped strongly enough such that the motions are dissipated locally, the behaviour takes on a diffusive character. Detailed discussions of the physics of these discs can be found in, for example, Papaloizou & Pringle (1983), Papaloizou & Lin (1995), Ogilvie (1999, 2000), Lubow & Ogilvie (2000), and Lodato & Pringle (2007), with some of the key points summarized in the review by Nixon & King (2016). To make progress, it is often useful to parametrize the turbulence in the disc (Shakura & Sunyaev 1973). In planar, circular, thin discs, this is usually done by taking the R–ϕ component of the vertically integrated stress tensor proportional to the vertically integrated pressure, with the constant of proportionality equal to a dimensionless parameter α. Equivalently, the kinematic viscosity ν = αcsH. For a warped disc, it is also necessary to determine the contributions to the radial communication of angular momentum from the R–z and ϕ–z components of the stress tensor. This was first done using the linearized fluid equations by Papaloizou & Pringle (1983), and subsequently for the non-linear case by Ogilvie (1999, 2000) and Ogilvie & Latter (2013). In these investigations, the damping of shearing motions is assumed to be isotropic, i.e. the same α causing dissipation from azimuthal shear acts to dissipate vertical shear. As discussed in the above paragraph, damping of these motions inhibits radial communication of the misaligned component of angular momentum, leading to an inverse dependence between the α parameter and the strength of the torque. A simple, but informative, analytical description of this relation is provided in section 4.1 of Lodato & Pringle (2007). The isotropy of the α parameter used in almost all warped disc modelling has been the subject of much speculation within several recent MHD simulation studies that include MRI turbulence directly. These investigations reported a perceived departure from the behaviour of an α model (Sorathia, Krolik & Hawley 2013a,b; Zhuravlev et al. 2014; Morales Teixeira et al. 2014; Krolik & Hawley 2015). Unfortunately, none of these works compared their results directly with an α model. But there has been some recent progress on this issue. The α formalism provides a simple method for encapsulating the dominant effects of turbulence within the disc. For the general case of weak magnetic fields, the turbulence is found to be a local fluid property, with little power on scales larger than the disc thickness ( ≲ 10 per cent; Simon, Beckwith & Armitage 2012). As the disc cannot be warped on scales smaller than the disc thickness, it seems reasonable to assume that the turbulence will not affect any particularly directed shear differently – suggesting the assumption of isotropy is reasonable. However, this cannot be taken as a given. A direct comparison between different codes is difficult due to the range of numerical techniques and different physics employed.1 In spite of this, Nealon et al. (2016) endeavoured to match the parameters of the MHD simulation presented in Krolik & Hawley (2015) and perform this calculation with an α viscosity playing the role of the disc turbulence to discern any differences found. However, no differences were found. Therefore, Nealon et al. (2016) directly confirm that, at least for the parameters simulated, the α model provides a suitable description of MRI turbulence for modelling warped discs. It is clear though that the dynamics of α discs will diverge from MHD simulations if the magnetic field strength approaches dynamical importance and can provide local and large-scale stresses that are comparable to the resonant Reynolds stress created by the warp itself. This is highly likely in, for example, simulations of MAD discs (e.g. McKinney, Tchekhovskoy & Blandford 2013), or in discs around highly magnetized stars. In this study, we adopt the α formalism, in which recent numerical investigations show that accretion discs may break or tear into distinct planes when the viscosity is not strong enough to communicate the precession induced in the disc. Disc tearing has been shown to occur in discs inclined to the spin of a central black hole (Nixon et al. 2012; Nealon, Price & Nixon 2015), in circumbinary discs around misaligned central binary systems (Facchini, Lodato & Price 2013; Nixon, King & Price 2013; Aly et al. 2015), and in circumprimary discs misaligned with respect to the binary orbital plane (Doğan et al. 2015). In these initial papers, the criterion for disc tearing has been derived simplistically by comparing the viscous torque with the precession torque induced in the disc, or in the case of wavelike warp propagation, the wave travel time was compared to the precession time-scale. Initially, the precession torque was compared to the torque arising from azimuthal shear (Nixon et al. 2012). For small α and small disc inclination angles, the simulations of Doğan et al. (2015) showed that the inclusion of the effective viscosity arising from vertical shear was required. However, as this term depends on the disc structure through the warp amplitude |ψ|, which cannot be determined a priori, it is not very useful as a predictive criterion. Further to these simulations of disc tearing, it is also possible for warped discs to break without any external forcing (Lodato & Price 2010). For these isolated warped discs, the evolution must be driven by the dependence of the effective viscosities on the disc structure. This instability occurs in a simpler environment than disc tearing, and likely underpins that behaviour. In this work, we seek to understand the physics causing an isolated warped disc to break. We do this by performing a stability analysis of the warped disc equations, finding a dispersion relation and the instability criterion. As we shall see below, instability appears in the form of anti-diffusion, which leads to a discontinuity in the disc angular momentum. This can be achieved by a discontinuity in the surface density or in the disc inclination. In the case of a flat disc, this is the familiar (Lightman & Eardley 1974) viscous disc instability, which underlies the thermal-viscous instability commonly thought to be the basis of dwarf nova outbursts (Meyer & Meyer-Hofmeister 1982), and when combined with central X-ray irradiation, thought to be the basis of soft X-ray transient outbursts (King & Ritter 1998). The thermal-viscous stability of a warped disc was briefly analysed by Ogilvie (2000), who concluded that a disc that is Lightman–Eardley stable (like those we consider here) is usually viscously stable when warped except for α ≲ 0.05 and the warp exceeds a critical amplitude. For a warped disc, anti-diffusion can lead to a separating of the disc into two or several discrete planes, which are connected by a series of sharply changing, low surface-density orbits. The layout of the paper is as follows. In Section 2, we describe the governing equations and their assumptions. In Section 3, we obtain the dispersion relation and derive the instability condition to determine a general criterion for discs to break. Further, we provide numerical evaluation of the growth rates of the instability. In Section 4, we provide a simple physical picture of the instability. In Section 5, we discuss our findings in the context of disc tearing where the disc is subject to an external torque. Finally, in Section 6, we provide our conclusions. 2 WARPED DISC EVOLUTION EQUATIONS In this work, we are considering an isolated, Keplerian, warped disc with α > H/R. In the literature, there exist various notations for the underlying equations with some subtleties that we outline here. 2.1 Definitions from previous works Using the notation of Pringle (1992), the evolution equation for the disc angular momentum is   \begin{eqnarray} \frac{\mathrm{\partial} \boldsymbol {L}}{\mathrm{\partial} t} &=& \frac{1}{R} \frac{\mathrm{\partial} }{\mathrm{\partial} R} \left\lbrace \frac{\left(\mathrm{\partial} / \mathrm{\partial} R \right) \left[\nu _{1}\Sigma R^{3}\left(-\Omega ^{^{\prime }} \right) \right] }{\Sigma \left( \mathrm{\partial} / \mathrm{\partial} R \right) \left(R^{2} \Omega \right)} \boldsymbol {L}\right\rbrace +\ \ {}\frac{1}{R}\frac{\mathrm{\partial} }{\mathrm{\partial} R}\left[\frac{1}{2} \nu _{2}R \left| \boldsymbol {L} \right|\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} R} \right]+\ \ {}\frac{1}{R}\frac{\mathrm{\partial} }{\mathrm{\partial} R} \left\lbrace \left[\frac{\frac{1}{2}\nu _{2}R^{3}\Omega \left|\mathrm{\partial} \boldsymbol {l} / \mathrm{\partial} R \right| ^{2}}{\left( \mathrm{\partial} / \mathrm{\partial} R \right) \left( R^{2} \Omega \right)} + \nu _{1}\left(\frac{R \Omega ^{^{\prime }}}{\Omega } \right) \right] \boldsymbol {L}\right\rbrace \nonumber \\ && +\ \ {}\frac{1}{R}\frac{\mathrm{\partial} }{\mathrm{\partial} R} \left(\nu _{3} R \left|\boldsymbol {L} \right| \boldsymbol {l} \times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} R} \right), \end{eqnarray} (1)where ν1 and ν2 are the effective viscosities, ν3 is the coefficient of a dispersive wavelike torque, Ω(R) is the orbital angular velocity of each annulus of the disc, Σ(R, t) is the disc surface density, $$\boldsymbol {l}\left(R,t\right)$$ is the unit angular momentum vector pointing perpendicular to the local orbital plane, and Ω΄ = ∂Ω/∂R. Also $$\boldsymbol {L} = \Sigma R^2 \Omega \boldsymbol {l}$$ is the angular momentum per unit area. This equation (1) contains four distinct terms responsible for the radial communication of orbital angular momentum (ν1 term), the radial communication of the misaligned component of angular momentum (ν2 term), an advective term resulting from the first two torques, and a precessional torque between misaligned rings (ν3 term). This equation was derived by Pringle (1992) from conservation of mass and angular momentum, but excluding the ν3 term which is not required by conservation. The ν3 term arises as the imaginary part of the complex diffusion coefficient in the fluid analysis of Papaloizou & Pringle (1983). The linearized analysis of Papaloizou & Pringle (1983) provides the leading order corrections to the coefficients αi of the effective viscosities (ν1 and ν2) and the dispersive torque (ν3) in the limit of small α and small warp amplitude. Developing this theory further, Ogilvie (1999), and subsequently Ogilvie (2000) and Ogilvie & Latter (2013), provide a non-linear derivation of this equation valid for general α and warp amplitude. On the whole, these works confirm the above equation, but include the important effect that the local fluid properties have on the torque coefficients. Ogilvie (1999) presents (1) as2  \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\Sigma \bar{\upsilon }_r r^3\Omega \boldsymbol {l}\right) = \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_1\mathcal {I}r^2\Omega ^2\boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\mathcal {I}r^3\Omega ^2\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_3\mathcal {I}r^3\Omega ^2\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right), \end{eqnarray} (2)where $$\bar{\upsilon }_r$$ is the mean radial velocity, $$\mathcal {I}$$ is the azimuthally averaged second vertical moment of the density, and Qi are the dimensionless torque coefficients. The coefficient Q1 represents the usual viscous torque tending to spin the ring up or down. Coefficient Q2 represents the torque diffusing the disc tilt. This torque tends to align the ring with its neighbours and flatten the disc. Coefficient Q3 represents the precession torque causing the rings misaligned with its neighbours to precess. The Q3 term is responsible for the dispersive wave-like propagation of the warp in the non-Keplerian inviscid case, but as we shall see below is relatively unimportant for the usual Keplerian diffusive case. The equivalence between (1) and (2) can be seen through the equation for the radial velocity   \begin{eqnarray} \bar{\upsilon }_r = \frac{\left(\mathrm{\partial} /\mathrm{\partial} r\right)\left(\nu _1\Sigma r^3\Omega ^\prime \right)-\frac{1}{2}\nu _2\Sigma r^3 \Omega \left|\mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r\right|^2}{r\Sigma \left(\mathrm{\partial} /\mathrm{\partial} r\right)\left(r^2\Omega \right)}, \end{eqnarray} (3)and the definitions of the effective viscosities   \begin{eqnarray} \nu _1=\frac{(-Q_1)\mathcal {I}\Omega }{q\Sigma }, \end{eqnarray} (4)  \begin{eqnarray} \nu _2=\frac{2Q_2\mathcal {I}\Omega }{\Sigma }, \end{eqnarray} (5)  \begin{eqnarray} \nu _3=\frac{Q_3\mathcal {I}\Omega }{\Sigma }, \end{eqnarray} (6)where   \begin{eqnarray} q=-\frac{\text{d} \ln \Omega }{\text{d} \ln r}. \end{eqnarray} (7)In Ogilvie (1999), the coefficients Qi are calculated for a polytropic disc. However, Ogilvie (1999) does not provide a complete analysis as the behaviour of $$\mathcal {I}$$ with warp amplitude $$\left|\psi \right| = r|\mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r|$$ is not calculated there, and thus, there is no simple equation relating $$\mathcal {I}$$ with Σ. Instead, this equation takes the form $$\mathcal {I} = f(\left|\psi \right|)\Sigma c_{\rm s}^2/\Omega ^2$$. This additional function f(|ψ|) is missing from the calculations in Nixon & King (2012), who assumed $$\mathcal {I} = \Sigma c_{\rm s}^2/\Omega ^2$$. Ogilvie (2000) provides a more complete thermal treatment including viscous heating and radiative cooling. This included an azimuthally dependent vertical treatment of the warp and thus the effect of vertical squeezing of the disc by the warp. In Ogilvie (2000), the evolution equation is presented as (2), complemented by a relation between $$\mathcal {I}$$ and Σ. In this relation, the dependence on the warp is included with a Q5 parameter. Finally, more recently, Ogilvie & Latter (2013) presented a further set of equations entirely consistent with those above, but with a subtle change to the definitions. They give the evolution equations as general conservation equations, with the internal torque defined separately. Putting them together yields   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\Sigma \bar{\upsilon }_r r^3\Omega \boldsymbol {l}\right) = \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_1\Sigma c_{\rm s}^2 r^2\boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\Sigma c_{\rm s}^2 r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_3\Sigma c_{\rm s}^2 r^3\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right). \end{eqnarray} (8)Inspection of (8) and (2) appears to show that $$\mathcal {I} = \Sigma c_{\rm s}^2/\Omega ^2$$, and that the |ψ| dependence of $$\mathcal {I}$$ has been ignored. However, this is not the case. Instead, this dependence has been included through a change in definition of the Qi coefficients between these works. We have been unable to find direct reference to this change of definition. The motivation is clear though, and provided by Ogilvie & Latter (2013), who remark that this form of the torque is natural for the isothermal discs under consideration there (and here). 2.2 Definitions in this work In this work, we adopt the following formalism based on the above discussion. This is in line with Ogilvie & Latter (2013) and is most relevant to the isothermal discs we consider here. We note that although we are using the coefficients from Ogilvie & Latter (2013), which are calculated assuming a locally isothermal equation of state, our results are more widely applicable. For example, Ogilvie (2000), who calculates the coefficients including an explicit energy equation, assumes that the disc is highly optically thick. This means that the azimuthal variations in physical variables, for example, ρ, T are essentially adiabatic (ttherm ≫ tdyn). It follows that in this highly optically thick case, calculation of the coefficients assuming an isothermal equation of state and those employing a polytropic equation of state or including an explicit energy equation differ only in the adiabatic exponent γ. Thus, we expect the physics in each case to differ only marginally. This appears to be confirmed by inspection of, for example, figs 3 and 4 of Ogilvie (1999) and figs 2 and 3 of Ogilvie (2000). We conclude from this that the calculations we present below, which are exact for the isothermal case (and thus also relevant to the optically thin case), are also applicable to the highly optically thick case. If by chance ttherm ∼ tdyn then the azimuthal variations are no longer adiabatic, but we expect the basic physics to remain unchanged.3 For our purposes, the conservation of mass equation is   \begin{eqnarray} \frac{\mathrm{\partial} \Sigma }{\mathrm{\partial} t}+\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(r\bar{\upsilon}_r \Sigma)=0, \end{eqnarray} (9)and the conservation of angular momentum equation is   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) +\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(\Sigma \bar{\upsilon }_r r^3\Omega \boldsymbol {l}) =\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_1\Sigma c_{\rm s}^2 r^2\boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\Sigma c_{\rm s}^2 r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_3\Sigma c_{\rm s}^2 r^3\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right). \end{eqnarray} (10)The effective viscosities are related to the Qi coefficients by   \begin{eqnarray} \nu _1 = (-Q_1)c_{\rm s}^2/q\Omega , \end{eqnarray} (11)  \begin{eqnarray} \nu _2 = 2Q_2 c_{\rm s}^2/\Omega , \end{eqnarray} (12)  \begin{eqnarray} \nu _3 = Q_3 c_{\rm s}^2/\Omega , \end{eqnarray} (13)and the Qi coefficients are evaluated using a code provided by Ogilvie following the calculations in Ogilvie & Latter (2013). By also defining the specific angular momentum   \begin{eqnarray} h = r^2\Omega , \end{eqnarray} (14)with h΄ = dh/dr, and substituting in the radial velocity, we have   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) =\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[Q_1\Sigma c_{\rm s}^2 r^2\boldsymbol {l} + Q_2\Sigma c_{\rm s}^2 r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r} + Q_3\Sigma c_{\rm s}^2 r^3\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r} - \left(\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[Q_1\Sigma c_{\rm s}^2 r^2\right] - Q_2\Sigma c_{\rm s}^2r\left|\psi \right|^2\right)\frac{h}{h^\prime }\boldsymbol {l}\right]. \end{eqnarray} (15)This equation, supplemented by the Qi coefficients from Ogilvie & Latter (2013), is what we will work with for the remainder of this paper. 3 STABILITY ANALYSIS OF THE WARPED DISC EQUATIONS 3.1 Dispersion relation We consider the local stability of an isolated warped disc, with α > H/R. From (15), we define the internal torque components gi(r, Σ, α, αb, q, |ψ|) that pre-multiply the dimensionless basis vectors $$(\boldsymbol {l},r\mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r,r\boldsymbol {l}\times \mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r)$$ as (Ogilvie 2000)   \begin{eqnarray} g_i=Q_i(\alpha ,\alpha _b,q,|\psi |)\Sigma c_{\rm s}^2r^2. \end{eqnarray} (16)In this work, we assume that the shear viscosity parameter α is a constant, and that the bulk viscosity parameter αb = 0. We focus on Keplerian discs with q = 3/2, and thus Qi = Qi(α, |ψ|). Using (16), we can rewrite (15) as   \begin{eqnarray} h\frac{\mathrm{\partial} }{\mathrm{\partial} t}(\Sigma \boldsymbol {l})=\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[{\it g}_1\boldsymbol {l}+{\it g}_2r\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_3r\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}-\left( \frac{\mathrm{\partial} {\it g}_1}{\mathrm{\partial} r}-\frac{{\it g}_2|\psi |^2}{r}\right)\frac{h}{h^{\prime }}\boldsymbol {l}\right]. \end{eqnarray} (17)Following Ogilvie (2000), we consider the stability of any solution of (17) with respect to linear perturbations ($$\delta \Sigma , \delta \boldsymbol {l}$$). The perturbed form of the angular momentum equation is   \begin{eqnarray} \begin{array}{ll}h\displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} t}(\delta \Sigma \boldsymbol {l}+\Sigma \delta \boldsymbol {l})=&{\displaystyle \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\Bigg [\delta {\it g}_1\boldsymbol {l}+{\it g}_1\delta \boldsymbol {l}+\delta {\it g}_2r\displaystyle \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_2r\displaystyle \frac{\mathrm{\partial} \delta \boldsymbol {l}}{\mathrm{\partial} r}+\delta {\it g}_3r\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_3r\delta \boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_3r\boldsymbol {l}\times \frac{\mathrm{\partial} \delta \boldsymbol {l}}{\mathrm{\partial} r}}\\ &{-\left( \displaystyle \frac{\mathrm{\partial} \delta {\it g}_1}{\mathrm{\partial} r}-\displaystyle \frac{\delta ({\it g}_2|\psi |^2)}{r}\right)\displaystyle \frac{h}{h^{\prime }}\boldsymbol {l}-\left( \displaystyle \frac{\mathrm{\partial} {\it g}_1}{\mathrm{\partial} r}-\displaystyle \frac{{\it g}_2|\psi |^2}{r}\right)\displaystyle \frac{h}{h^{\prime }}\delta \boldsymbol {l}\Bigg ]}. \end{array} \end{eqnarray} (18)Here, the perturbed quantities are defined by   \begin{eqnarray} \delta {\it g}_i=\frac{\mathrm{\partial} {\it g}_i}{\mathrm{\partial} \Sigma }\delta \Sigma +\frac{\mathrm{\partial} {\it g}_i}{\mathrm{\partial} \left|\psi \right|}\delta \left|\psi \right|, \end{eqnarray} (19)  \begin{eqnarray} \delta |\psi |=\frac{r^2}{|\psi |}\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}{\cdot }\frac{\mathrm{\partial} \delta \boldsymbol {l}}{\mathrm{\partial} r}, \end{eqnarray} (20)and $$\boldsymbol {l}\cdot \delta \boldsymbol {l}=0$$. We now search for solutions to the perturbations of the form   \begin{eqnarray} \exp \left(-\text{i}\int \omega \text{d}t + i\int k \text{d}r \right), \end{eqnarray} (21)where ω is the wave frequency and k is the wavenumber. This requires |kr| ≫ 1, and while formally k → ∞ is acceptable, the equations become invalid as k approaches 1/H (Ogilvie 2000). Similarly, |ω| must be large enough such that the perturbations grow or decay on time-scales short compared to the background solution. We adopt an arbitrary scaling $$|\delta \boldsymbol {l}|=O(k^{-1})$$, then δΣ = O(1), δ|ψ| = O(1), and δgi = O(1), whereas ω = O(k2) as expected for a diffusion equation (Ogilvie 2000). If we define an orthonormal basis ($$\boldsymbol {l}, {\boldsymbol {m}}, {\boldsymbol {n}}$$) by $$\boldsymbol {l}\times {\boldsymbol m}={\boldsymbol n}$$ and $$\delta \boldsymbol {l}=\delta m{\boldsymbol {m}} + \delta n{\boldsymbol {n}}$$ (where $$r(\mathrm{\partial} \boldsymbol {l} / \mathrm{\partial} r)=|\psi |{\boldsymbol {m}}$$), we obtain the $$\boldsymbol {l}$$, $${\boldsymbol {m}}$$, and $${\boldsymbol {n}}$$ components of (18), at a leading order in k, as follows:   \begin{eqnarray} -\text{i}\omega hr\delta \Sigma =k^2\frac{h}{h^{\prime }}\delta {\it g}_1, \end{eqnarray} (22)  \begin{eqnarray} -\text{i}\omega hr\Sigma \delta m=\text{i}k|\psi |\delta {\it g}_2-k^2{\it g}_2r\delta m+k^2{\it g}_3r\delta n-\text{i}k|\psi |\frac{h}{rh^{\prime }}\delta {\it g}_1, \end{eqnarray} (23)  \begin{eqnarray} -\text{i}\omega hr\Sigma \delta n=\text{i}k|\psi |\delta {\it g}_3-k^2{\it g}_2r\delta n-k^2{\it g}_3r\delta m, \end{eqnarray} (24)where δ|ψ| = ikrδm. We can express δgi in terms of δΣ and δ|ψ| by differentiating (16), then write δ|ψ| in terms of δm. As a result, we get three linear equations with three unknowns (δΣ, δm, δn), which yield a coefficients determinant:   \begin{eqnarray} {\rm det} {\left[\begin{array}{c{@}{\quad }c{@}{\quad }c}s-a Q_1 & -aQ_1^{\prime } & 0 \\ (Q_2-aQ_1) |\psi | & s+Q_2+Q_2^{\prime }|\psi |-aQ_1^{\prime }|\psi | & - Q_3 \\ Q_3|\psi | & Q_3+Q_3^{\prime }|\psi | & s+ Q_2 \end{array}\right]} =0. \end{eqnarray} (25)Here, the prime on Qi represents differentiation with respect to |ψ|, a = h/rh΄ = dln r/dln h = 1/(2 − q) = 2 (for a Keplerian disc with q = 3/2), s is defined by   \begin{eqnarray} s=-\frac{{\rm i}\omega }{\Omega }\Bigg (\frac{\Omega }{c_{\rm s} k}\Bigg )^2, \end{eqnarray} (26)and we shall call ℜ[s] the dimensionless growth rate. s is related to the physical growth rate (ℜ[ − iω]) through   \begin{eqnarray} -{\rm i}\omega =s{\Omega }\Bigg (\frac{c_{\rm s} k}{\Omega }\Bigg )^2 \end{eqnarray} (27)and the perturbations grow (ℜ[s] > 0) or decay (ℜ[s] < 0) as exp [ℜ(−iω)t]. This shows that the physical growth rate is faster on smaller length-scales as expected for a diffusion equation. As the disc pressure scalelength is defined by Hp ≡ cs/Ω, the dimensionless growth rate can be rewritten as s = −i(ω/Ω)(Hpk)−2. Therefore, s ∼ 1 indicates strong instability on dynamical time-scales and length-scales ∼Hp. We again note that validity of the equations requires k ≲ 1/H. The determinant (25) gives a third-order dispersion relation:   \begin{eqnarray} &&\!\!\!\!&&{s^3-s^2\left[aQ_1-2Q_2+|\psi |\big (aQ^{\prime }_1-Q^{\prime }_2\big )\right]-s\left[2a{Q}_1{Q}_2-Q_2^2-Q_3^2+|\psi |\big (a{Q}_1{Q}^{\prime }_2-Q_2Q^{\prime }_2-Q_3Q^{\prime }_3\big )\right]}\nonumber\\ &&-a\left[Q_1(Q_2^2+Q_3^2)+|\psi |\big (Q_1Q_2Q^{\prime }_2-Q_1^{\prime }Q_2^2+Q_1Q_3Q^{\prime }_3-Q_1^{\prime }Q_3^2\big ) \right]=0. \end{eqnarray} (28)The disc becomes unstable if any of the roots of (28) has a positive real part, i.e. ℜ(s) > 0, as the perturbations then grow exponentially with time. For completeness, we note, following the discussion of the thermodynamics in Section 2.2, that the equivalent of (28) for an analysis explicitly including the energy equation is given by (cf. Ogilvie 2000)   \begin{eqnarray} &&{s^3-s^2\left[apQ_1-2Q_2+|\psi |\big (aQ^{\prime }_1-Q^{\prime }_2\Big )\right]-s\left[2ap{Q}_1{Q}_2-Q_2^2-Q_3^2+|\psi |\big (ap{Q}_1{Q}^{\prime }_2-Q_2Q^{\prime }_2-Q_3Q^{\prime }_3 -a(p-1)Q^{\prime }_1Q_2\big )\right]}\nonumber\\ &&-ap\left[Q_1(Q_2^2+Q_3^2)+|\psi |\big (Q_1Q_2Q^{\prime }_2-Q_1^{\prime }Q_2^2+Q_1Q_3Q^{\prime }_3-Q_1^{\prime }Q_3^2\big ) \right]=0, \end{eqnarray} (29)where p is the exponent of Σ in νΣ (or equivalently the exponent of Σ in $$\mathcal {I}$$). In a planar disc, the usual thermal-viscous stability criterion (Lightman & Eardley 1974) is given by p > 0. For the isothermal discs we consider here, we have p = 1 (for which 28 and 29 are identical). For Thomson scattering opacity, p = 5/3, and for Kramers opacity, p = 10/7, and as we have noted above, the coefficients do not vary strongly between these cases. Therefore, we expect the calculations we present here (with p = 1 and coefficients derived from an isothermal equation of state) to hold at least qualitatively for a wide range of thermodynamic treatments. 3.2 Instability criterion and the growth rates Before arriving at the full criterion (provided below), we first give the roots of (28) for some simplified cases. We will see below that these cases are generally useful for both understanding the more complex cases, and as close approximations to the full solution. We label the roots as s(Qi), where the subscript i denotes the torque coefficient(s) that we include in the solution. First, if we set all Qi terms in the dispersion relation to zero apart from Q1 (i.e. Q2 = Q3 = 0), we find one root which is given as   \begin{eqnarray} s(Q_1)= a\big (Q_1 + Q^{\prime }_1|\psi |\big ) = -a\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(-Q_1\left|\psi \right|\right). \end{eqnarray} (30)For instability, one requires ℜ[s(Q1)] > 0. For |ψ| = 0, this corresponds to the usual criterion for viscous instability in a planar disc (given by ∂(νΣ)/∂Σ < 0; Lightman & Eardley 1974, although we note that the relevant term is simplified here as we have assumed a simple equation of state, cf. equation 83 of Ogilvie 2000). However, in the warped case, the criterion is modified with an extra term determined by the behaviour of the coefficient Q1 with warp amplitude. Fig. 1 (left-hand panel) shows the variation of the dimensionless growth rates of this mode, ℜ[s(Q1)] with |ψ| for α = 0.01, 0.03, and 0.1. The disc becomes unstable after some critical warp amplitude, |ψ|c. The critical warp amplitude required for instability is very small for discs with low α. The growth rates are a modest fraction ≳0.1 of dynamical. As a result, discs can easily become unstable for this case. Figure 1. View largeDownload slide Left-hand panel: This shows the dimensionless growth rate ℜ[s(Q1)] (equation 30) as a function of dimensionless disc warp amplitude |ψ| caused by the term Q1 in the equations, which governs radial angular momentum transfer, for three different values of the viscosity parameter α. Instability occurs when ℜ[s] > 0. We note that the maximal growth rates ≈ℜ[s]Ω ∼ 0.3Ω are comparable with dynamical and can occur for small values of |ψ|, and remain unstable for large values of |ψ|. Right-hand panel: This shows the dimensionless growth rate ℜ[s(Q2)] caused by the action of the term Q2, which governs the viscous flattening of the warped disc. We note that this term shows no instability for large values of α. For smaller values of α, there is a range of disc warp amplitude |ψ| for which instability occurs. When it does occur, the time-scale for instability can be short with growth rates ≈ℜ[s]Ω ∼ Ω being dynamical. In both panels, the grey line represents zero growth rate. We remind the reader that these illustrative calculations are based on the artificial assumption that two of the three Qi are set to zero. Figure 1. View largeDownload slide Left-hand panel: This shows the dimensionless growth rate ℜ[s(Q1)] (equation 30) as a function of dimensionless disc warp amplitude |ψ| caused by the term Q1 in the equations, which governs radial angular momentum transfer, for three different values of the viscosity parameter α. Instability occurs when ℜ[s] > 0. We note that the maximal growth rates ≈ℜ[s]Ω ∼ 0.3Ω are comparable with dynamical and can occur for small values of |ψ|, and remain unstable for large values of |ψ|. Right-hand panel: This shows the dimensionless growth rate ℜ[s(Q2)] caused by the action of the term Q2, which governs the viscous flattening of the warped disc. We note that this term shows no instability for large values of α. For smaller values of α, there is a range of disc warp amplitude |ψ| for which instability occurs. When it does occur, the time-scale for instability can be short with growth rates ≈ℜ[s]Ω ∼ Ω being dynamical. In both panels, the grey line represents zero growth rate. We remind the reader that these illustrative calculations are based on the artificial assumption that two of the three Qi are set to zero. When we take into account only Q2 (i.e. Q1 = Q3 = 0), then the solution of (28) gives   \begin{eqnarray} s(Q_2)=-\big (Q_2+Q^{\prime }_2|\psi |\big ) = -\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_2\left|\psi \right|\right). \end{eqnarray} (31)The disc becomes unstable when ℜ[s(Q2)] > 0. This condition is similar to (30), but is caused by instability in the warp diffusion coefficient (Q2) and its effect on the warp amplitude. This is distinct from the usual viscous instability as it acts on the warp amplitude rather than the disc surface density. This mode becomes unstable at some critical warp amplitude interval. Fig. 1 (right-hand panel) shows the dimensionless growth rates found for this mode. We see that for small α, the growth rates of this instability are significantly higher than those found in the first simplified case where we included only the Q1 term. However, for large α, there is no instability from s(Q2). We note that there is a second root in this case, which is equal to −Q2, but this root is never unstable within the α formalism. We should also note that in these simplified cases, s is real. Thus, the dimensionless growth rates are equal to s. Below we will encounter complex values of s in which case the dimensionless growth rate is given by ℜ(s). If this analysis is repeated for (29) as opposed to (28), then (30) is modified to $$s = apQ_1 + a\left|\psi \right|Q^{\prime }_1$$, and thus, the Lightman & Eardley (1974) criterion for thermal-viscous instability is modified from apQ1 > 0 to $$apQ_1 > -a\left|\psi \right|Q^{\prime }_1$$. In contrast, (31) remains unchanged. Now we include both Q1 and Q2 terms (i.e. Q3 = 0) to find a more general criterion for instability. In this case, (28) gives three roots, one of which corresponds to −Q2 again, and is thus stable for any |ψ|. The remaining roots of the cubic correspond to   \begin{eqnarray} s_{\pm }(Q_1,Q_2)&=&\frac{1}{2}\bigg [a\big (Q_1 + Q^{\prime }_1|\psi |\big )-\big (Q_2 +Q^{\prime }_2 |\psi |\big )\nonumber\\ &&\pm \,\sqrt{\left[a \big (Q_1 + Q^{\prime }_1 |\psi |\big )-\big (Q_2 + Q^{\prime }_2 |\psi |\big ) \right]^2 +4a\left[\left(Q_1Q_2 + (Q_1Q^{\prime }_2 -Q^{\prime }_1 Q_2)|\psi |\right)\right]}\bigg ]. \end{eqnarray} (32)We see that these modes include both instabilities found in the simplified cases above. Inspection of (32) reveals that ℜ(s+) ≥ ℜ(s−). Therefore, when the disc is unstable, s+ represents the mode that shows more rapid growth. Fig. 2 shows the variation of ℜ[s+(Q1, Q2)] with |ψ|. The disc becomes unstable when ℜ[s+(Q1, Q2)] > 0. We plot ℜ[s(Q1)] and ℜ[s(Q2)] in the same graph for comparison. ℜ[s+(Q1, Q2)] shows similar behaviour to ℜ[s(Q2)] given by (31) for |ψ| ≲ 2.5, and shows similar behaviour to ℜ[s(Q1)] given by (30) for |ψ| ≳ 2.5. The disc becomes unstable at some critical warp amplitude, |ψ|c. For some values of α, the disc becomes stable again at larger warp amplitude before becoming unstable again at a second critical warp amplitude. The critical warp amplitudes and the maximum values of the dimensionless growth rates are almost equivalent for s+(Q1, Q2) and s(Q2). The Q1 term tends to become dominant only at high warp amplitudes. We should note that ℜ[s−(Q1, Q2)] also becomes positive at some critical warp amplitude, but is very close to zero. As a result, both modes, i.e. s+(Q1, Q2) and s−(Q1, Q2), are unstable after the critical warp amplitude, but s+(Q1, Q2) is always the most unstable. Figure 2. View largeDownload slide Here, we combine the two panels of Fig. 1 to illustrate the combined effects of the two terms Q1 and Q2. The curves show the dimensionless growth rates ℜ[s+(Q1, Q2)] (black solid line) with |ψ| (equation 32), compared against ℜ[s(Q1)] (green solid line) and ℜ[s(Q2)] (red dashed line). The grey line represents zero growth rate. For small α = 0.01 (left-hand panel), we see that instability occurs for all warp amplitudes |ψ| > 0.2, except for a very narrow region around |ψ| ≈ 2.96. The growth rate is generally large, except around |ψ| ≈ 3.0. The largest growth rates occur between 0.2 < |ψ| ≲ 2.5, where the Q2 term dominates. For larger |ψ|, the Q1 term dominates and the growth rates are sub-dynamical. For intermediate α = 0.03 (right-hand panel), the same applies except that the growth rates becomes smaller, and there is a wider band of stability around |ψ| = 2.2. For large |ψ|, the growth rates in the two cases (α = 0.01 and 0.03) are similar. Figure 2. View largeDownload slide Here, we combine the two panels of Fig. 1 to illustrate the combined effects of the two terms Q1 and Q2. The curves show the dimensionless growth rates ℜ[s+(Q1, Q2)] (black solid line) with |ψ| (equation 32), compared against ℜ[s(Q1)] (green solid line) and ℜ[s(Q2)] (red dashed line). The grey line represents zero growth rate. For small α = 0.01 (left-hand panel), we see that instability occurs for all warp amplitudes |ψ| > 0.2, except for a very narrow region around |ψ| ≈ 2.96. The growth rate is generally large, except around |ψ| ≈ 3.0. The largest growth rates occur between 0.2 < |ψ| ≲ 2.5, where the Q2 term dominates. For larger |ψ|, the Q1 term dominates and the growth rates are sub-dynamical. For intermediate α = 0.03 (right-hand panel), the same applies except that the growth rates becomes smaller, and there is a wider band of stability around |ψ| = 2.2. For large |ψ|, the growth rates in the two cases (α = 0.01 and 0.03) are similar. Finally, if we include all the Qi terms in the dispersion relation, we obtain the full criterion for instability as ℜ[s1(Q1, Q2, Q3)] > 0 with   \begin{eqnarray} s_{1}(Q_1,Q_2,Q_3)= \frac{1}{6}\big [2\mathcal {C}_1+2^{2/3}\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{1/3}+ 2^{4/3}\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{-1/3}\mathcal {C}_4\big ], \end{eqnarray} (33)where   \begin{eqnarray} \mathcal {C}_1 &=&\ \ a \big (Q_1+Q_1^{\prime } |\psi | \big )-2 \big (Q_2+Q_2^{\prime } |\psi |\big ), \end{eqnarray} (34)  \begin{eqnarray} \mathcal {C}_2 &=&\ \ 3 a \bigg [Q_1^{\prime } |\psi | \bigg (2 Q_2^{\prime 2} |\psi |^2 + 5 Q_2 Q_2^{\prime } |\psi | - 4 Q_2^2 - 3 Q_3 \Big (Q_3^{\prime } |\psi | + 4 Q_3\Big )\bigg ) + Q_1 \bigg (-Q_2^{\prime 2} |\psi |^2 + 2 Q_2 Q_2^{\prime } |\psi |\nonumber\\ &&+\, 2 \big (Q_2^2 + 3 Q_3 \big (Q_3^{\prime } |\psi | + Q_3 \big )\big )\bigg )\bigg ]+2 a^3 \big (Q_1 +Q_1^{\prime } |\psi |\big )^3 + 3 a^2 \big (Q_1^2 -2 Q_1^{\prime 2} |\psi |^2 - Q_1 Q_1^{\prime } |\psi | \big ) \big (2 Q_2 + Q_2^{\prime } |\psi | \big )\nonumber\\ &&-\, \big (2 Q_2+Q_2^{\prime } |\psi |\big ) \big (2 Q_2^{\prime 2} |\psi |^2 - Q_2 Q_2^{\prime } |\psi | - Q_2^2 - 9 Q_3 \big (Q_3^{\prime } |\psi | + Q_3 \big )\big ), \end{eqnarray} (35)  \begin{eqnarray} \mathcal {C}_3&=&\ \ 4 \bigg (3 \big (Q_2^2+Q_2^{\prime } Q_2|\psi |-a Q_1 (2 Q_2+Q_2^{\prime } |\psi |)+Q_3 (Q_3+Q_3^{\prime } |\psi |)\Big )-\big (2 Q_2+Q_2^{\prime } |\psi | -a (Q_1+Q_1^{\prime } |\psi |)\big )^2\bigg )^3\nonumber\\ &&+\,\Big [2 a^3 (Q_1+Q_1^{\prime } |\psi |)^3+(2 Q_2+Q_2^{\prime } |\psi |) \left(Q_2^2+Q_2^{\prime } Q_2|\psi |+9 Q_3^2-2 Q_2^{\prime 2} |\psi | ^2+9 Q_3^{\prime } Q_3 |\psi | \right)\nonumber\\ &&+\,3 a \bigg (Q_1 \left(2 Q_2^2+2 Q_2^{\prime } Q_2|\psi |+6 Q_3^2-Q_2^{\prime 2} |\psi | ^2+6 Q_3^{\prime } Q_3 |\psi | \right)+Q_1^{\prime } |\psi | \big (-4 Q_2^2+5 Q_2^{\prime } Q_2|\psi | +2 Q_2^{\prime 2} |\psi | ^2\nonumber\\ &&-\,12 Q_3^2 -3 Q_3^{\prime } Q_3 |\psi | \big )\bigg )+3 a^2 (2 Q_2+Q_2^{\prime } |\psi |) \left(Q_1^2-Q_1^{\prime } Q_1|\psi |-2 Q_1^{\prime 2} |\psi | ^2\right)\big ]^2,\ \rm and \end{eqnarray} (36)  \begin{eqnarray} \mathcal {C}_4&= &\ \ a^2 \big (Q_1+Q_1^{\prime } |\psi | \big )^2+ Q_2^2+ Q_2^{\prime } Q_2|\psi | +Q_2^{\prime 2} |\psi | ^2 +a \big (2 Q_2+Q_2^{\prime } |\psi | \big ) \big (Q_1-2 Q_1^{\prime } |\psi | \big )-3 \big (Q_3^2 + Q_3^{\prime } Q_3 |\psi |\big ). \end{eqnarray} (37)The other two roots of the dispersion relation are given by   \begin{eqnarray} s_{2}(Q_1,Q_2,Q_3)= \frac{1}{6}\big [2\mathcal {C}_1-2^{-1/3}\big (1-\mathrm{i}\sqrt{3}\big )\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{1/3}-2^{1/3}\big (1+\mathrm{i}\sqrt{3}\big )\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{-1/3}\mathcal {C}_4\big ] \end{eqnarray} (38)and   \begin{eqnarray} s_{3}(Q_1,Q_2,Q_3)= \frac{1}{6}\big [2\mathcal {C}_1-2^{-1/3}\big (1+\mathrm{i}\sqrt{3}\big )\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{1/3}-2^{1/3}\big (1-\mathrm{i}\sqrt{3}\big )\Big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{-1/3}\mathcal {C}_4\big ]. \end{eqnarray} (39) These roots are either all real, or there can be one real root and two complex conjugates. For determining instability, we are interested in the real roots and the real parts of the complex roots. While there is always one real root, the other two roots can swap between real or complex conjugates for different values of |ψ|. One root is always stable. Of the other two roots, one always gives the fastest growth (i.e. when the real part is positive it shows the largest value) and this root is given by (33). For some values of |ψ|, this root is a complex conjugate and thus the other root has the same growth rate. In Fig. 3, we show the dimensionless growth rates of the instability found from the full criterion (33) for different values of α. We see that the disc becomes unstable for sufficiently large |ψ| values. The critical warp amplitudes where the disc becomes unstable are smaller and the growth rates of the instability are higher for low α. In Fig. 4, we plot the maximum values of the dimensionless growth rates of instability (ℜ[s]max) against α. We see that ℜ[s]max values are inversely proportional to α. Finally, the critical warp amplitudes for instability are plotted against α in Fig. 5. This figure shows the stable and unstable regions in the (α, |ψ|) parameter space. Figure 3. View largeDownload slide Here, we show the dimensionless growth rates ℜ[s] as functions of |ψ| taking account of all three terms, Q1, Q2, and Q3, for different values of α. The effect of the Q3 term is negligible for the parameters we have considered (see Figs 6 and 7). Figure 3. View largeDownload slide Here, we show the dimensionless growth rates ℜ[s] as functions of |ψ| taking account of all three terms, Q1, Q2, and Q3, for different values of α. The effect of the Q3 term is negligible for the parameters we have considered (see Figs 6 and 7). Figure 4. View largeDownload slide This plot shows the maximum growth rate of the instability as a function of viscosity parameter α. Note that for different values of α the maximum growth occurs for different values of the dimensionless warp amplitude |ψ|. This plot shows that instability occurs for all values of α in this range. For small values of α, the maximum growth rate is dynamical. For larger values of α, the maximum growth rate is smaller, but still comparable to dynamical (∼0.2Ω). Figure 4. View largeDownload slide This plot shows the maximum growth rate of the instability as a function of viscosity parameter α. Note that for different values of α the maximum growth occurs for different values of the dimensionless warp amplitude |ψ|. This plot shows that instability occurs for all values of α in this range. For small values of α, the maximum growth rate is dynamical. For larger values of α, the maximum growth rate is smaller, but still comparable to dynamical (∼0.2Ω). Figure 5. View largeDownload slide Stable (white) and unstable (green) regions in the (α, |ψ|) parameter space. This plot represents the critical warp amplitudes for instability to occur in discs with various α values. We note that a flat disc is always stable (assuming an isothermal equation of state so that there is no Lightman–Eardley-like instability). For 0.01 ≤ α ≤ 0.2, a nearly flat disc with |ψ| ≲ 0.1 is also always stable. For a given value of α, there is always a minimum value of the warp amplitude, which gives rise to instability. Figure 5. View largeDownload slide Stable (white) and unstable (green) regions in the (α, |ψ|) parameter space. This plot represents the critical warp amplitudes for instability to occur in discs with various α values. We note that a flat disc is always stable (assuming an isothermal equation of state so that there is no Lightman–Eardley-like instability). For 0.01 ≤ α ≤ 0.2, a nearly flat disc with |ψ| ≲ 0.1 is also always stable. For a given value of α, there is always a minimum value of the warp amplitude, which gives rise to instability. If we compare (32) and (33), we see that the inclusion of Q3 produces a big difference in the formula, but does not have a strong impact on the numerical results. In Fig. 6, we compare two criteria, (32), where Q3 = 0, and (33), where Q3 ≠ 0. We see that the inclusion of Q3 has a very small effect on the critical warp amplitudes and the growth rates. We also plot the ratios of the |ψ|c and ℜ[s]max found for the two different cases, s(Q1, Q2) and s(Q1, Q2, Q3). Fig. 7 shows the values of X and Y as a function of α, where X = |ψ|cℜ[s(Q1, Q2, Q3)]/|ψ|cℜ[s(Q1, Q2)] and Y = ℜ[smax(Q1, Q2, Q3)]/ℜ[smax(Q1, Q2)]. The solutions are in excellent agreement. Therefore, (32) can be used to decide whether the disc is stable or unstable, at least for the parameters we have explored. Figure 6. View largeDownload slide Comparison of the full (33) and simplified criteria (32) for instability for α = 0.01 and 0.1. The grey line represents zero growth rate. These plots show that the simplified criteria are an extremely good approximation to both the location of the critical warp amplitude(s) and the dimensionless growth rate as a function of warp amplitude. Figure 6. View largeDownload slide Comparison of the full (33) and simplified criteria (32) for instability for α = 0.01 and 0.1. The grey line represents zero growth rate. These plots show that the simplified criteria are an extremely good approximation to both the location of the critical warp amplitude(s) and the dimensionless growth rate as a function of warp amplitude. Figure 7. View largeDownload slide X = |ψ|cℜ[s(Q1, Q2, Q3)]/|ψ|cℜ[s(Q1, Q2)] and Y = ℜ[smax(Q1, Q2, Q3)]/ℜ[smax(Q1, Q2)] ratios for different α values. The grey line represents unity. Thus, the critical warp amplitudes and the maximum values of the dimensionless growth rates are captured by the simplified root (32) to high accuracy for all of the α values we have explored. Figure 7. View largeDownload slide X = |ψ|cℜ[s(Q1, Q2, Q3)]/|ψ|cℜ[s(Q1, Q2)] and Y = ℜ[smax(Q1, Q2, Q3)]/ℜ[smax(Q1, Q2)] ratios for different α values. The grey line represents unity. Thus, the critical warp amplitudes and the maximum values of the dimensionless growth rates are captured by the simplified root (32) to high accuracy for all of the α values we have explored. In conclusion, the criteria for instability can be expressed from (32) as below:   \begin{eqnarray} {\rm If} \ \ \,\left[a\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_1\left|\psi \right|\right)-\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_2\left|\psi \right|\right)\right]>0,\,\rm the \,\rm disc \,\rm is \,\rm unstable, \end{eqnarray} (40)  \begin{eqnarray} {\rm or \ \rm if} \ \ \ \left[a\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_1\left|\psi \right|\right)-\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_2\left|\psi \right|\right)\right]<0,\ \ \ \ \ {\rm and} \ \ \ \ \ 4a\left[\left(Q_1Q_2 + (Q_1Q^{\prime }_2 -Q^{\prime }_1 Q_2)|\psi |\right)\right]>0,\,\rm the \,\rm disc \,\rm \,is \,\rm also\, \rm unstable. \end{eqnarray} (41) 4 HEURISTIC PICTURE OF THE VISCOUS-WARP INSTABILITY To elucidate the fundamental physics of this viscous-warp disc instability from the complex analysis described above, we present here a simplified analysis. Starting from (15), we simplify the picture aggressively. In a warped disc with α ≪ 1, the smoothing of the disc warp occurs faster than mass is transported radially (Papaloizou & Pringle 1983). So we assume that the important torque is the diffusive torque relating to ν2. Our exploration of the analysis above confirms this is a sensible approach in many, but not all, cases. Thus,   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) \sim \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\Sigma c_{\rm s}^2r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right). \end{eqnarray} (42)Interpreting this as a diffusion equation implies that the flux of misaligned angular momentum is   \begin{eqnarray} \boldsymbol {F}_{\rm mis} = Q_2\Sigma c_{\rm s}^2r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r} \end{eqnarray} (43)with rate of transfer $${\sim } \left|\boldsymbol {F}_{\rm mis}/\Sigma r^2 \Omega \right| = (1/2)\nu _2\left|\psi \right|$$. Connecting this with classical instability of a diffusion equation yields the criterion for instability as   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(\nu _2\left|\psi \right|\right) < 0. \end{eqnarray} (44)This can also be trivially recovered by performing a stability analysis on (42), with $$\boldsymbol {l} \rightarrow \boldsymbol {l}+\delta \boldsymbol {l}$$, δQ2 = (∂Q2/∂|ψ|)δ|ψ| and assuming δΣ = 0. The criterion (44) can be interpreted as stating if ν2(|ψ|) is such that the maximum diffusion rate is not located at maxima in warp amplitude, then local maxima in warp amplitude will grow, and the disc will break. As this will result in rapid transfer of mass out of this region due to the large warp amplitude implying large torques, this will also be realized by a significant drop in local surface density. This resembles the Lightman–Eardley viscous instability but for a warped disc with the warp amplitude playing the role of the surface density. Such behaviour is impossible if the effective viscosities are constant with varying warp amplitude. This is confirmed in our analysis above, which yields no unstable solutions if Qi are independent of |ψ|. This is consistent with the results of Nixon & King (2012).4 Finally, we note that when the correct limits are applied, the full stability analysis presented above reproduces the ad hoc answer derived here (see equation 31). 5 APPLICATION TO DISC TEARING In this work, we have focused on disc ‘breaking’, where an isolated already-warped disc (i.e. with no external torque) evolves to create a discontinuity between two or more planes, with a series of low-density orbits connecting them over a small radial extent. We have revealed above that this occurs due to a viscous-warp instability, which occurs when either (40) or (41) is satisfied. In many physical systems, the disc is subject to an external torque often occurring due to precession (e.g. the Lense–Thirring effect). In this case, the warp initially grows with time. If the disc ‘breaks’ under these conditions creating disconnected rings that keep precessing effectively independently, then this is what we have called disc ‘tearing’ in previous papers (e.g. Nixon et al. 2012; Doğan et al. 2015; Nealon et al. 2015). It is of great interest to connect the analysis in this work to the idea of disc tearing. Unfortunately at first appearances, we cannot perform the same analysis as above for disc tearing, where an external forced precession is imposed, as the background unperturbed state is no longer necessarily an equilibrium solution on short time-scales. It may be that this is unimportant, and an approximate solution can be found that is accurate enough. We will check this in future work. For our purposes here, we repeat that the analysis above is a local analysis that takes a given warped disc shape and determines its stability. For disc tearing, we need to consider how the disc gets into a warped configuration, and at what stage during that process instability occurs. Here, we can make a reasoned guess at what the disc tearing criterion is likely to be. Above we have derived, for various system parameters, the critical warp amplitude |ψ|c at which the disc becomes unstable. It would seem reasonable that a forced disc would tear up if there is no |ψ| < |ψ|c at which the torque attempting to flatten the disc (∝ Q2|ψ|) is stronger than the torque induced by, for example, precession. As discussed in Doğan et al. (2015), without the knowledge of the value of |ψ|c, this calculation cannot be done a priori without performing a full numerical simulation. However, our analysis enables exactly that. We will test this in the future by comparing this reasoning with numerical simulations. We note that given the viscous torque and precession torque (in the Lense-Thirring case) have no dependence on inclination angle, this criterion at first appears at odds with previous simulations that show, for example, a strong inclination angle dependence. There is also no dependence on the disc thickness, but previous simulations show that thinner discs are easier to tear. The answer to this is likely to be that the warp amplitude cannot be larger than a value dependent on the inclination angle and H/R. This can be seen by the following argument. First, the warp amplitude is   \begin{eqnarray} \left|\psi \right| = r\left|\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right|, \end{eqnarray} (45)and the disc unit tilt vector is   \begin{eqnarray} \boldsymbol {l} = \left(\cos \gamma \sin \beta ,\sin \gamma \sin \beta ,\cos \beta \right), \end{eqnarray} (46)where γ is the disc twist and β is the disc tilt (e.g. Nixon 2012). Let us assume γ = 0, and a maximal warp from β to zero (the same result, equation 48, is arrived at if β is assumed a constant and γ varies from 0 to π). Thus,   \begin{eqnarray} \delta \boldsymbol {l} = \left(\sin \beta ,0,\cos \beta -1\right), \end{eqnarray} (47)which has magnitude (2 − 2cos β)1/2. We can also assume that δr cannot become smaller than H, and thus the maximum warp amplitude is given by   \begin{eqnarray} \left|\psi \right|_{\rm max} \sim \left(2-2\cos \beta \right)^{1/2}\frac{R}{H} \sim \beta \frac{R}{H}. \end{eqnarray} (48)Therefore, discs that are already close to alignment, or are thick, are difficult to tear, as the maximum warp amplitude is small and potentially less than the critical warp amplitude required for instability. 6 CONCLUSIONS We have derived the criterion for isolated warped discs to break. This occurs physically due to viscous anti-diffusion of the warp amplitude. The resulting large torques that are induced lead to the disc surface density lowering in the unstable regions. Thus, the instability leads to a succession of orbits with planes that change rapidly with radius, with only a small amount of mass present. This is a local instability so it is not likely to directly affect global disc properties such as the central accretion rate. However, this instability underlies the process of disc tearing which has the capacity to dramatically alter the instantaneous accretion rate and the observable properties of the disc on short time-scales. The instability is governed by a complicated criterion (33), but we have shown that the criterion can be simplified in most cases to a more manageable criteria given by (40) and (41). This simply implies that the diffusion of the warp must be maximal at the location of the maximum warp amplitude otherwise the warp will grow until the disc disconnects on scales of order the disc thickness ∼H. In future works, we will perform this analysis for non-Keplerian discs, connect our results with numerical simulations, and include a more sophisticated thermal treatment. For discs with moderate values of α ∼ 0.01–0.1 and reasonably thin H/R ≲ α, our work supports the idea of discs being viscously unstable to disconnecting into discrete planes. This process underpins the state-transition cycle suggested by Nixon & Salvesen (2014). The natural link between the viscous instability in planar discs, which underlies the thermal-viscous outburst cycle observed in many types of binary systems (many of which are thought to contain warped discs), and the viscous instability of warped discs is appealing. This may be a way to introduce hysteresis through a change in the disc stability with warping. We will explore this in future work. ACKNOWLEDGEMENTS We thank the referee for providing a helpful and detailed report. SD gratefully acknowledges the warm hospitality of the Theoretical Astrophysics Group at University of Leicester during her visit. SD is supported by the Turkish Scientific and Technical Research Council (TÜBİTAK - 117F280). CJN is supported by the Science and Technology Facilities Council (grant number ST/M005917/1). The Theoretical Astrophysics Group at the University of Leicester is supported by an STFC Consolidated Grant. Footnotes 1 To date only the SPH code phantom (Price et al. 2017) has been tested for correctly modelling warped discs against the non-linear theory of Ogilvie (1999), finding excellent agreement in the simulated parameter space. 2 Note that both r here and R in (1) represent the local cylindrical radius of the annulus. 3 In a steady disc, ttherm ∼ α−1tdyn (Pringle 1981), and so, in general, ttherm ≳ tdyn with rough equality only for α ≈ 1. 4 Note that in Nixon & King (2012) when using the Qi coefficients from Ogilvie (1999), we smoothed the viscosity coefficients over the neighbouring few rings to avoid what was perceived there as numerical issues. Actually, this behaviour was the physical instability discussed here. REFERENCES Aly H., Dehnen W., Nixon C., King A., 2015, MNRAS , 449, 65 https://doi.org/10.1093/mnras/stv128 CrossRef Search ADS   Bate M. R., Lodato G., Pringle J. E., 2010, MNRAS , 401, 1505 https://doi.org/10.1111/j.1365-2966.2009.15773.x CrossRef Search ADS   Doğan S., Nixon C., King A., Price D. J., 2015, MNRAS , 449, 1251 https://doi.org/10.1093/mnras/stv347 CrossRef Search ADS   Facchini S., Lodato G., Price D. J., 2013, MNRAS , 433, 2142 https://doi.org/10.1093/mnras/stt877 CrossRef Search ADS   King A. R., Ritter H., 1998, MNRAS , 293, L42 https://doi.org/10.1046/j.1365-8711.1998.01295.x CrossRef Search ADS   Krolik J. H., Hawley J. F., 2015, ApJ , 806, 141 https://doi.org/10.1088/0004-637X/806/1/141 CrossRef Search ADS   Lense J., Thirring H., 1918, Phys. Z. , 19, 156 Lightman A. P., Eardley D. M., 1974, ApJ , 187, L1 https://doi.org/10.1086/181377 CrossRef Search ADS   Lodato G., Price D. J., 2010, MNRAS , 405, 1212 Lodato G., Pringle J. E., 2007, MNRAS , 381, 1287 https://doi.org/10.1111/j.1365-2966.2007.12332.x CrossRef Search ADS   Lubow S. H., 1992, ApJ , 398, 525 https://doi.org/10.1086/171877 CrossRef Search ADS   Lubow S. H., Ogilvie G. I., 2000, ApJ , 538, 326 https://doi.org/10.1086/309101 CrossRef Search ADS   Lucas W. E., Bonnell I. A., Davies M. B., Rice W. K. M., 2013, MNRAS , 433, 353 https://doi.org/10.1093/mnras/stt727 CrossRef Search ADS   McKinney J. C., Tchekhovskoy A., Blandford R. D., 2013, Science , 339, 49 https://doi.org/10.1126/science.1230811 CrossRef Search ADS PubMed  Meyer F., Meyer-Hofmeister E., 1982, A&A , 106, 34 Morales Teixeira D., Fragile P. C., Zhuravlev V. V., Ivanov P. B., 2014, ApJ , 796, 103 https://doi.org/10.1088/0004-637X/796/2/103 CrossRef Search ADS   Nealon R., Price D. J., Nixon C. J., 2015, MNRAS , 448, 1526 https://doi.org/10.1093/mnras/stv014 CrossRef Search ADS   Nealon R., Nixon C., Price D. J., King A., 2016, MNRAS , 455, L62 Nixon C. J., 2012, MNRAS , 423, 2597 https://doi.org/10.1111/j.1365-2966.2012.21072.x CrossRef Search ADS   Nixon C. J., King A. R., 2012, MNRAS , 421, 1201 https://doi.org/10.1111/j.1365-2966.2011.20377.x CrossRef Search ADS   Nixon C., King A., 2016, Lecture Notes in Physics, Vol. 905, Astrophysical Black Holes . Springer International Publishing, Switzerland, p. 45 Nixon C., Salvesen G., 2014, MNRAS , 437, 3994 https://doi.org/10.1093/mnras/stt2215 CrossRef Search ADS   Nixon C., King A., Price D., Frank J., 2012, ApJ , 757, L24 https://doi.org/10.1088/2041-8205/757/2/L24 CrossRef Search ADS   Nixon C., King A., Price D., 2013, MNRAS , 434, 1946 https://doi.org/10.1093/mnras/stt1136 CrossRef Search ADS   Ogilvie G. I., 1999, MNRAS , 304, 557 https://doi.org/10.1046/j.1365-8711.1999.02340.x CrossRef Search ADS   Ogilvie G. I., 2000, MNRAS , 317, 607 https://doi.org/10.1046/j.1365-8711.2000.03654.x CrossRef Search ADS   Ogilvie G. I., Latter H. N., 2013, MNRAS , 433, 2403 https://doi.org/10.1093/mnras/stt916 CrossRef Search ADS   Papaloizou J. C. B., Lin D. N. C., 1995, ApJ , 438, 841 https://doi.org/10.1086/175127 CrossRef Search ADS   Papaloizou J. C. B., Pringle J. E., 1983, MNRAS , 202, 1181 https://doi.org/10.1093/mnras/202.4.1181 CrossRef Search ADS   Papaloizou J. C. B., Terquem C., 1995, MNRAS , 274, 987 Price D. J. et al.  , 2017, PASA , preprint (arXiv:1702.03930) Pringle J. E., 1981, ARA&A , 19, 137 CrossRef Search ADS   Pringle J. E., 1992, MNRAS , 258, 811 https://doi.org/10.1093/mnras/258.4.811 CrossRef Search ADS   Pringle J. E., 1996, MNRAS , 281, 357 https://doi.org/10.1093/mnras/281.1.357 CrossRef Search ADS   Pringle J. E., 1997, MNRAS , 292, 136 https://doi.org/10.1093/mnras/292.1.136 CrossRef Search ADS   Schandl S., Meyer F., 1994, A&A , 289, 149 Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Simon J. B., Beckwith K., Armitage P. J., 2012, MNRAS , 422, 2685 https://doi.org/10.1111/j.1365-2966.2012.20835.x CrossRef Search ADS   Sorathia K. A., Krolik J. H., Hawley J. F., 2013a, ApJ , 768, 133 https://doi.org/10.1088/0004-637X/768/2/133 CrossRef Search ADS   Sorathia K. A., Krolik J. H., Hawley J. F., 2013b, ApJ , 777, 21 https://doi.org/10.1088/0004-637X/777/1/21 CrossRef Search ADS   Zhuravlev V. V., Ivanov P. B., Fragile P. C., Morales Teixeira D., 2014, ApJ , 796, 104 https://doi.org/10.1088/0004-637X/796/2/104 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Instability of warped discs

, Volume 476 (2) – May 1, 2018
13 pages

/lp/ou_press/instability-of-warped-discs-vuLcXHkn7a
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty155
Publisher site
See Article on Publisher Site

### Abstract

Abstract Accretion discs are generally warped. If a warp in a disc is too large, the disc can ‘break’ apart into two or more distinct planes, with only tenuous connections between them. Further, if an initially planar disc is subject to a strong differential precession, then it can be torn apart into discrete annuli that precess effectively independently. In previous investigations, torque-balance formulae have been used to predict where and when the disc breaks into distinct parts. In this work, focusing on discs with Keplerian rotation and where the shearing motions driving the radial communication of the warp are damped locally by turbulence (the ‘diffusive’ regime), we investigate the stability of warped discs to determine the precise criterion for an isolated warped disc to break. We find and solve the dispersion relation, which, in general, yields three roots. We provide a comprehensive analysis of this viscous-warp instability and the emergent growth rates and their dependence on disc parameters. The physics of the instability can be understood as a combination of (1) a term that would generally encapsulate the classical Lightman–Eardley instability in planar discs (given by ∂(νΣ)/∂Σ < 0) but is here modified by the warp to include ∂(ν1|ψ|)/∂|ψ| < 0, and (2) a similar condition acting on the diffusion of the warp amplitude given in simplified form by ∂(ν2|ψ|)/∂|ψ| < 0. We discuss our findings in the context of discs with an imposed precession, and comment on the implications for different astrophysical systems. accretion, accretion discs, black hole physics, hydrodynamics, instabilities 1 INTRODUCTION Accretion discs occur throughout astrophysics, from the birthsites of stars and planets to the feeding of supermassive black holes (SMBHs) in galaxy centres. Through several different effects, each important in different contexts, these discs can warp. Discs can be formed in a warped configuration through, for example, chaotic infall on to protostellar discs (Bate, Lodato & Pringle 2010) or SMBHs (Lucas et al. 2013). Alternatively, discs can warp if they orbit in a non-spherical gravitational potential that induces radial differential precession, for example, the Lense–Thirring effect from a spinning black hole (Lense & Thirring 1918) or tides from a companion star (Papaloizou & Terquem 1995). Further, there are warping instabilities that act on planar discs to drive a disc warp from an initially axisymmetric configuration with only a small perturbation, for example, radiation warping (Pringle 1996, 1997), winds (Schandl & Meyer 1994), or resonant tidal effects (Lubow 1992; Lubow & Ogilvie 2000). Once a disc is warped, the local plane of the disc changes with radius. This creates a displacement of the high-pressure mid-plane region from one ring of the disc to its neighbour which oscillates with azimuth. At the descending (azimuthal angle ϕ = 0) and ascending (ϕ = π) nodes, where the neighbouring rings cross, there is no displacement, whereas it is maximal at ϕ = π/2 and (with opposite direction) at 3π/2 (see e.g. fig. 1 of Nixon & King 2016). Thus, a fluid element orbiting in the disc feels a radial pressure gradient that oscillates on the local orbital frequency. The fluid element then exhibits epicyclic motion. In a near-Keplerian disc, the epicyclic frequency is close to the orbital frequency, creating a resonance between the motion and the forcing. The epicyclic motion communicates the disc warp radially, and if undamped, it will lead to the launch of a bending wave. However, in the presence of turbulence, the wave will be damped, and if damped strongly enough such that the motions are dissipated locally, the behaviour takes on a diffusive character. Detailed discussions of the physics of these discs can be found in, for example, Papaloizou & Pringle (1983), Papaloizou & Lin (1995), Ogilvie (1999, 2000), Lubow & Ogilvie (2000), and Lodato & Pringle (2007), with some of the key points summarized in the review by Nixon & King (2016). To make progress, it is often useful to parametrize the turbulence in the disc (Shakura & Sunyaev 1973). In planar, circular, thin discs, this is usually done by taking the R–ϕ component of the vertically integrated stress tensor proportional to the vertically integrated pressure, with the constant of proportionality equal to a dimensionless parameter α. Equivalently, the kinematic viscosity ν = αcsH. For a warped disc, it is also necessary to determine the contributions to the radial communication of angular momentum from the R–z and ϕ–z components of the stress tensor. This was first done using the linearized fluid equations by Papaloizou & Pringle (1983), and subsequently for the non-linear case by Ogilvie (1999, 2000) and Ogilvie & Latter (2013). In these investigations, the damping of shearing motions is assumed to be isotropic, i.e. the same α causing dissipation from azimuthal shear acts to dissipate vertical shear. As discussed in the above paragraph, damping of these motions inhibits radial communication of the misaligned component of angular momentum, leading to an inverse dependence between the α parameter and the strength of the torque. A simple, but informative, analytical description of this relation is provided in section 4.1 of Lodato & Pringle (2007). The isotropy of the α parameter used in almost all warped disc modelling has been the subject of much speculation within several recent MHD simulation studies that include MRI turbulence directly. These investigations reported a perceived departure from the behaviour of an α model (Sorathia, Krolik & Hawley 2013a,b; Zhuravlev et al. 2014; Morales Teixeira et al. 2014; Krolik & Hawley 2015). Unfortunately, none of these works compared their results directly with an α model. But there has been some recent progress on this issue. The α formalism provides a simple method for encapsulating the dominant effects of turbulence within the disc. For the general case of weak magnetic fields, the turbulence is found to be a local fluid property, with little power on scales larger than the disc thickness ( ≲ 10 per cent; Simon, Beckwith & Armitage 2012). As the disc cannot be warped on scales smaller than the disc thickness, it seems reasonable to assume that the turbulence will not affect any particularly directed shear differently – suggesting the assumption of isotropy is reasonable. However, this cannot be taken as a given. A direct comparison between different codes is difficult due to the range of numerical techniques and different physics employed.1 In spite of this, Nealon et al. (2016) endeavoured to match the parameters of the MHD simulation presented in Krolik & Hawley (2015) and perform this calculation with an α viscosity playing the role of the disc turbulence to discern any differences found. However, no differences were found. Therefore, Nealon et al. (2016) directly confirm that, at least for the parameters simulated, the α model provides a suitable description of MRI turbulence for modelling warped discs. It is clear though that the dynamics of α discs will diverge from MHD simulations if the magnetic field strength approaches dynamical importance and can provide local and large-scale stresses that are comparable to the resonant Reynolds stress created by the warp itself. This is highly likely in, for example, simulations of MAD discs (e.g. McKinney, Tchekhovskoy & Blandford 2013), or in discs around highly magnetized stars. In this study, we adopt the α formalism, in which recent numerical investigations show that accretion discs may break or tear into distinct planes when the viscosity is not strong enough to communicate the precession induced in the disc. Disc tearing has been shown to occur in discs inclined to the spin of a central black hole (Nixon et al. 2012; Nealon, Price & Nixon 2015), in circumbinary discs around misaligned central binary systems (Facchini, Lodato & Price 2013; Nixon, King & Price 2013; Aly et al. 2015), and in circumprimary discs misaligned with respect to the binary orbital plane (Doğan et al. 2015). In these initial papers, the criterion for disc tearing has been derived simplistically by comparing the viscous torque with the precession torque induced in the disc, or in the case of wavelike warp propagation, the wave travel time was compared to the precession time-scale. Initially, the precession torque was compared to the torque arising from azimuthal shear (Nixon et al. 2012). For small α and small disc inclination angles, the simulations of Doğan et al. (2015) showed that the inclusion of the effective viscosity arising from vertical shear was required. However, as this term depends on the disc structure through the warp amplitude |ψ|, which cannot be determined a priori, it is not very useful as a predictive criterion. Further to these simulations of disc tearing, it is also possible for warped discs to break without any external forcing (Lodato & Price 2010). For these isolated warped discs, the evolution must be driven by the dependence of the effective viscosities on the disc structure. This instability occurs in a simpler environment than disc tearing, and likely underpins that behaviour. In this work, we seek to understand the physics causing an isolated warped disc to break. We do this by performing a stability analysis of the warped disc equations, finding a dispersion relation and the instability criterion. As we shall see below, instability appears in the form of anti-diffusion, which leads to a discontinuity in the disc angular momentum. This can be achieved by a discontinuity in the surface density or in the disc inclination. In the case of a flat disc, this is the familiar (Lightman & Eardley 1974) viscous disc instability, which underlies the thermal-viscous instability commonly thought to be the basis of dwarf nova outbursts (Meyer & Meyer-Hofmeister 1982), and when combined with central X-ray irradiation, thought to be the basis of soft X-ray transient outbursts (King & Ritter 1998). The thermal-viscous stability of a warped disc was briefly analysed by Ogilvie (2000), who concluded that a disc that is Lightman–Eardley stable (like those we consider here) is usually viscously stable when warped except for α ≲ 0.05 and the warp exceeds a critical amplitude. For a warped disc, anti-diffusion can lead to a separating of the disc into two or several discrete planes, which are connected by a series of sharply changing, low surface-density orbits. The layout of the paper is as follows. In Section 2, we describe the governing equations and their assumptions. In Section 3, we obtain the dispersion relation and derive the instability condition to determine a general criterion for discs to break. Further, we provide numerical evaluation of the growth rates of the instability. In Section 4, we provide a simple physical picture of the instability. In Section 5, we discuss our findings in the context of disc tearing where the disc is subject to an external torque. Finally, in Section 6, we provide our conclusions. 2 WARPED DISC EVOLUTION EQUATIONS In this work, we are considering an isolated, Keplerian, warped disc with α > H/R. In the literature, there exist various notations for the underlying equations with some subtleties that we outline here. 2.1 Definitions from previous works Using the notation of Pringle (1992), the evolution equation for the disc angular momentum is   \begin{eqnarray} \frac{\mathrm{\partial} \boldsymbol {L}}{\mathrm{\partial} t} &=& \frac{1}{R} \frac{\mathrm{\partial} }{\mathrm{\partial} R} \left\lbrace \frac{\left(\mathrm{\partial} / \mathrm{\partial} R \right) \left[\nu _{1}\Sigma R^{3}\left(-\Omega ^{^{\prime }} \right) \right] }{\Sigma \left( \mathrm{\partial} / \mathrm{\partial} R \right) \left(R^{2} \Omega \right)} \boldsymbol {L}\right\rbrace +\ \ {}\frac{1}{R}\frac{\mathrm{\partial} }{\mathrm{\partial} R}\left[\frac{1}{2} \nu _{2}R \left| \boldsymbol {L} \right|\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} R} \right]+\ \ {}\frac{1}{R}\frac{\mathrm{\partial} }{\mathrm{\partial} R} \left\lbrace \left[\frac{\frac{1}{2}\nu _{2}R^{3}\Omega \left|\mathrm{\partial} \boldsymbol {l} / \mathrm{\partial} R \right| ^{2}}{\left( \mathrm{\partial} / \mathrm{\partial} R \right) \left( R^{2} \Omega \right)} + \nu _{1}\left(\frac{R \Omega ^{^{\prime }}}{\Omega } \right) \right] \boldsymbol {L}\right\rbrace \nonumber \\ && +\ \ {}\frac{1}{R}\frac{\mathrm{\partial} }{\mathrm{\partial} R} \left(\nu _{3} R \left|\boldsymbol {L} \right| \boldsymbol {l} \times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} R} \right), \end{eqnarray} (1)where ν1 and ν2 are the effective viscosities, ν3 is the coefficient of a dispersive wavelike torque, Ω(R) is the orbital angular velocity of each annulus of the disc, Σ(R, t) is the disc surface density, $$\boldsymbol {l}\left(R,t\right)$$ is the unit angular momentum vector pointing perpendicular to the local orbital plane, and Ω΄ = ∂Ω/∂R. Also $$\boldsymbol {L} = \Sigma R^2 \Omega \boldsymbol {l}$$ is the angular momentum per unit area. This equation (1) contains four distinct terms responsible for the radial communication of orbital angular momentum (ν1 term), the radial communication of the misaligned component of angular momentum (ν2 term), an advective term resulting from the first two torques, and a precessional torque between misaligned rings (ν3 term). This equation was derived by Pringle (1992) from conservation of mass and angular momentum, but excluding the ν3 term which is not required by conservation. The ν3 term arises as the imaginary part of the complex diffusion coefficient in the fluid analysis of Papaloizou & Pringle (1983). The linearized analysis of Papaloizou & Pringle (1983) provides the leading order corrections to the coefficients αi of the effective viscosities (ν1 and ν2) and the dispersive torque (ν3) in the limit of small α and small warp amplitude. Developing this theory further, Ogilvie (1999), and subsequently Ogilvie (2000) and Ogilvie & Latter (2013), provide a non-linear derivation of this equation valid for general α and warp amplitude. On the whole, these works confirm the above equation, but include the important effect that the local fluid properties have on the torque coefficients. Ogilvie (1999) presents (1) as2  \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\Sigma \bar{\upsilon }_r r^3\Omega \boldsymbol {l}\right) = \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_1\mathcal {I}r^2\Omega ^2\boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\mathcal {I}r^3\Omega ^2\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_3\mathcal {I}r^3\Omega ^2\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right), \end{eqnarray} (2)where $$\bar{\upsilon }_r$$ is the mean radial velocity, $$\mathcal {I}$$ is the azimuthally averaged second vertical moment of the density, and Qi are the dimensionless torque coefficients. The coefficient Q1 represents the usual viscous torque tending to spin the ring up or down. Coefficient Q2 represents the torque diffusing the disc tilt. This torque tends to align the ring with its neighbours and flatten the disc. Coefficient Q3 represents the precession torque causing the rings misaligned with its neighbours to precess. The Q3 term is responsible for the dispersive wave-like propagation of the warp in the non-Keplerian inviscid case, but as we shall see below is relatively unimportant for the usual Keplerian diffusive case. The equivalence between (1) and (2) can be seen through the equation for the radial velocity   \begin{eqnarray} \bar{\upsilon }_r = \frac{\left(\mathrm{\partial} /\mathrm{\partial} r\right)\left(\nu _1\Sigma r^3\Omega ^\prime \right)-\frac{1}{2}\nu _2\Sigma r^3 \Omega \left|\mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r\right|^2}{r\Sigma \left(\mathrm{\partial} /\mathrm{\partial} r\right)\left(r^2\Omega \right)}, \end{eqnarray} (3)and the definitions of the effective viscosities   \begin{eqnarray} \nu _1=\frac{(-Q_1)\mathcal {I}\Omega }{q\Sigma }, \end{eqnarray} (4)  \begin{eqnarray} \nu _2=\frac{2Q_2\mathcal {I}\Omega }{\Sigma }, \end{eqnarray} (5)  \begin{eqnarray} \nu _3=\frac{Q_3\mathcal {I}\Omega }{\Sigma }, \end{eqnarray} (6)where   \begin{eqnarray} q=-\frac{\text{d} \ln \Omega }{\text{d} \ln r}. \end{eqnarray} (7)In Ogilvie (1999), the coefficients Qi are calculated for a polytropic disc. However, Ogilvie (1999) does not provide a complete analysis as the behaviour of $$\mathcal {I}$$ with warp amplitude $$\left|\psi \right| = r|\mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r|$$ is not calculated there, and thus, there is no simple equation relating $$\mathcal {I}$$ with Σ. Instead, this equation takes the form $$\mathcal {I} = f(\left|\psi \right|)\Sigma c_{\rm s}^2/\Omega ^2$$. This additional function f(|ψ|) is missing from the calculations in Nixon & King (2012), who assumed $$\mathcal {I} = \Sigma c_{\rm s}^2/\Omega ^2$$. Ogilvie (2000) provides a more complete thermal treatment including viscous heating and radiative cooling. This included an azimuthally dependent vertical treatment of the warp and thus the effect of vertical squeezing of the disc by the warp. In Ogilvie (2000), the evolution equation is presented as (2), complemented by a relation between $$\mathcal {I}$$ and Σ. In this relation, the dependence on the warp is included with a Q5 parameter. Finally, more recently, Ogilvie & Latter (2013) presented a further set of equations entirely consistent with those above, but with a subtle change to the definitions. They give the evolution equations as general conservation equations, with the internal torque defined separately. Putting them together yields   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(\Sigma \bar{\upsilon }_r r^3\Omega \boldsymbol {l}\right) = \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_1\Sigma c_{\rm s}^2 r^2\boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\Sigma c_{\rm s}^2 r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_3\Sigma c_{\rm s}^2 r^3\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right). \end{eqnarray} (8)Inspection of (8) and (2) appears to show that $$\mathcal {I} = \Sigma c_{\rm s}^2/\Omega ^2$$, and that the |ψ| dependence of $$\mathcal {I}$$ has been ignored. However, this is not the case. Instead, this dependence has been included through a change in definition of the Qi coefficients between these works. We have been unable to find direct reference to this change of definition. The motivation is clear though, and provided by Ogilvie & Latter (2013), who remark that this form of the torque is natural for the isothermal discs under consideration there (and here). 2.2 Definitions in this work In this work, we adopt the following formalism based on the above discussion. This is in line with Ogilvie & Latter (2013) and is most relevant to the isothermal discs we consider here. We note that although we are using the coefficients from Ogilvie & Latter (2013), which are calculated assuming a locally isothermal equation of state, our results are more widely applicable. For example, Ogilvie (2000), who calculates the coefficients including an explicit energy equation, assumes that the disc is highly optically thick. This means that the azimuthal variations in physical variables, for example, ρ, T are essentially adiabatic (ttherm ≫ tdyn). It follows that in this highly optically thick case, calculation of the coefficients assuming an isothermal equation of state and those employing a polytropic equation of state or including an explicit energy equation differ only in the adiabatic exponent γ. Thus, we expect the physics in each case to differ only marginally. This appears to be confirmed by inspection of, for example, figs 3 and 4 of Ogilvie (1999) and figs 2 and 3 of Ogilvie (2000). We conclude from this that the calculations we present below, which are exact for the isothermal case (and thus also relevant to the optically thin case), are also applicable to the highly optically thick case. If by chance ttherm ∼ tdyn then the azimuthal variations are no longer adiabatic, but we expect the basic physics to remain unchanged.3 For our purposes, the conservation of mass equation is   \begin{eqnarray} \frac{\mathrm{\partial} \Sigma }{\mathrm{\partial} t}+\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(r\bar{\upsilon}_r \Sigma)=0, \end{eqnarray} (9)and the conservation of angular momentum equation is   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) +\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(\Sigma \bar{\upsilon }_r r^3\Omega \boldsymbol {l}) =\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_1\Sigma c_{\rm s}^2 r^2\boldsymbol {l}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\Sigma c_{\rm s}^2 r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right) + \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_3\Sigma c_{\rm s}^2 r^3\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right). \end{eqnarray} (10)The effective viscosities are related to the Qi coefficients by   \begin{eqnarray} \nu _1 = (-Q_1)c_{\rm s}^2/q\Omega , \end{eqnarray} (11)  \begin{eqnarray} \nu _2 = 2Q_2 c_{\rm s}^2/\Omega , \end{eqnarray} (12)  \begin{eqnarray} \nu _3 = Q_3 c_{\rm s}^2/\Omega , \end{eqnarray} (13)and the Qi coefficients are evaluated using a code provided by Ogilvie following the calculations in Ogilvie & Latter (2013). By also defining the specific angular momentum   \begin{eqnarray} h = r^2\Omega , \end{eqnarray} (14)with h΄ = dh/dr, and substituting in the radial velocity, we have   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) =\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[Q_1\Sigma c_{\rm s}^2 r^2\boldsymbol {l} + Q_2\Sigma c_{\rm s}^2 r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r} + Q_3\Sigma c_{\rm s}^2 r^3\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r} - \left(\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[Q_1\Sigma c_{\rm s}^2 r^2\right] - Q_2\Sigma c_{\rm s}^2r\left|\psi \right|^2\right)\frac{h}{h^\prime }\boldsymbol {l}\right]. \end{eqnarray} (15)This equation, supplemented by the Qi coefficients from Ogilvie & Latter (2013), is what we will work with for the remainder of this paper. 3 STABILITY ANALYSIS OF THE WARPED DISC EQUATIONS 3.1 Dispersion relation We consider the local stability of an isolated warped disc, with α > H/R. From (15), we define the internal torque components gi(r, Σ, α, αb, q, |ψ|) that pre-multiply the dimensionless basis vectors $$(\boldsymbol {l},r\mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r,r\boldsymbol {l}\times \mathrm{\partial} \boldsymbol {l}/\mathrm{\partial} r)$$ as (Ogilvie 2000)   \begin{eqnarray} g_i=Q_i(\alpha ,\alpha _b,q,|\psi |)\Sigma c_{\rm s}^2r^2. \end{eqnarray} (16)In this work, we assume that the shear viscosity parameter α is a constant, and that the bulk viscosity parameter αb = 0. We focus on Keplerian discs with q = 3/2, and thus Qi = Qi(α, |ψ|). Using (16), we can rewrite (15) as   \begin{eqnarray} h\frac{\mathrm{\partial} }{\mathrm{\partial} t}(\Sigma \boldsymbol {l})=\frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left[{\it g}_1\boldsymbol {l}+{\it g}_2r\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_3r\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}-\left( \frac{\mathrm{\partial} {\it g}_1}{\mathrm{\partial} r}-\frac{{\it g}_2|\psi |^2}{r}\right)\frac{h}{h^{\prime }}\boldsymbol {l}\right]. \end{eqnarray} (17)Following Ogilvie (2000), we consider the stability of any solution of (17) with respect to linear perturbations ($$\delta \Sigma , \delta \boldsymbol {l}$$). The perturbed form of the angular momentum equation is   \begin{eqnarray} \begin{array}{ll}h\displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} t}(\delta \Sigma \boldsymbol {l}+\Sigma \delta \boldsymbol {l})=&{\displaystyle \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\Bigg [\delta {\it g}_1\boldsymbol {l}+{\it g}_1\delta \boldsymbol {l}+\delta {\it g}_2r\displaystyle \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_2r\displaystyle \frac{\mathrm{\partial} \delta \boldsymbol {l}}{\mathrm{\partial} r}+\delta {\it g}_3r\boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_3r\delta \boldsymbol {l}\times \frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}+{\it g}_3r\boldsymbol {l}\times \frac{\mathrm{\partial} \delta \boldsymbol {l}}{\mathrm{\partial} r}}\\ &{-\left( \displaystyle \frac{\mathrm{\partial} \delta {\it g}_1}{\mathrm{\partial} r}-\displaystyle \frac{\delta ({\it g}_2|\psi |^2)}{r}\right)\displaystyle \frac{h}{h^{\prime }}\boldsymbol {l}-\left( \displaystyle \frac{\mathrm{\partial} {\it g}_1}{\mathrm{\partial} r}-\displaystyle \frac{{\it g}_2|\psi |^2}{r}\right)\displaystyle \frac{h}{h^{\prime }}\delta \boldsymbol {l}\Bigg ]}. \end{array} \end{eqnarray} (18)Here, the perturbed quantities are defined by   \begin{eqnarray} \delta {\it g}_i=\frac{\mathrm{\partial} {\it g}_i}{\mathrm{\partial} \Sigma }\delta \Sigma +\frac{\mathrm{\partial} {\it g}_i}{\mathrm{\partial} \left|\psi \right|}\delta \left|\psi \right|, \end{eqnarray} (19)  \begin{eqnarray} \delta |\psi |=\frac{r^2}{|\psi |}\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}{\cdot }\frac{\mathrm{\partial} \delta \boldsymbol {l}}{\mathrm{\partial} r}, \end{eqnarray} (20)and $$\boldsymbol {l}\cdot \delta \boldsymbol {l}=0$$. We now search for solutions to the perturbations of the form   \begin{eqnarray} \exp \left(-\text{i}\int \omega \text{d}t + i\int k \text{d}r \right), \end{eqnarray} (21)where ω is the wave frequency and k is the wavenumber. This requires |kr| ≫ 1, and while formally k → ∞ is acceptable, the equations become invalid as k approaches 1/H (Ogilvie 2000). Similarly, |ω| must be large enough such that the perturbations grow or decay on time-scales short compared to the background solution. We adopt an arbitrary scaling $$|\delta \boldsymbol {l}|=O(k^{-1})$$, then δΣ = O(1), δ|ψ| = O(1), and δgi = O(1), whereas ω = O(k2) as expected for a diffusion equation (Ogilvie 2000). If we define an orthonormal basis ($$\boldsymbol {l}, {\boldsymbol {m}}, {\boldsymbol {n}}$$) by $$\boldsymbol {l}\times {\boldsymbol m}={\boldsymbol n}$$ and $$\delta \boldsymbol {l}=\delta m{\boldsymbol {m}} + \delta n{\boldsymbol {n}}$$ (where $$r(\mathrm{\partial} \boldsymbol {l} / \mathrm{\partial} r)=|\psi |{\boldsymbol {m}}$$), we obtain the $$\boldsymbol {l}$$, $${\boldsymbol {m}}$$, and $${\boldsymbol {n}}$$ components of (18), at a leading order in k, as follows:   \begin{eqnarray} -\text{i}\omega hr\delta \Sigma =k^2\frac{h}{h^{\prime }}\delta {\it g}_1, \end{eqnarray} (22)  \begin{eqnarray} -\text{i}\omega hr\Sigma \delta m=\text{i}k|\psi |\delta {\it g}_2-k^2{\it g}_2r\delta m+k^2{\it g}_3r\delta n-\text{i}k|\psi |\frac{h}{rh^{\prime }}\delta {\it g}_1, \end{eqnarray} (23)  \begin{eqnarray} -\text{i}\omega hr\Sigma \delta n=\text{i}k|\psi |\delta {\it g}_3-k^2{\it g}_2r\delta n-k^2{\it g}_3r\delta m, \end{eqnarray} (24)where δ|ψ| = ikrδm. We can express δgi in terms of δΣ and δ|ψ| by differentiating (16), then write δ|ψ| in terms of δm. As a result, we get three linear equations with three unknowns (δΣ, δm, δn), which yield a coefficients determinant:   \begin{eqnarray} {\rm det} {\left[\begin{array}{c{@}{\quad }c{@}{\quad }c}s-a Q_1 & -aQ_1^{\prime } & 0 \\ (Q_2-aQ_1) |\psi | & s+Q_2+Q_2^{\prime }|\psi |-aQ_1^{\prime }|\psi | & - Q_3 \\ Q_3|\psi | & Q_3+Q_3^{\prime }|\psi | & s+ Q_2 \end{array}\right]} =0. \end{eqnarray} (25)Here, the prime on Qi represents differentiation with respect to |ψ|, a = h/rh΄ = dln r/dln h = 1/(2 − q) = 2 (for a Keplerian disc with q = 3/2), s is defined by   \begin{eqnarray} s=-\frac{{\rm i}\omega }{\Omega }\Bigg (\frac{\Omega }{c_{\rm s} k}\Bigg )^2, \end{eqnarray} (26)and we shall call ℜ[s] the dimensionless growth rate. s is related to the physical growth rate (ℜ[ − iω]) through   \begin{eqnarray} -{\rm i}\omega =s{\Omega }\Bigg (\frac{c_{\rm s} k}{\Omega }\Bigg )^2 \end{eqnarray} (27)and the perturbations grow (ℜ[s] > 0) or decay (ℜ[s] < 0) as exp [ℜ(−iω)t]. This shows that the physical growth rate is faster on smaller length-scales as expected for a diffusion equation. As the disc pressure scalelength is defined by Hp ≡ cs/Ω, the dimensionless growth rate can be rewritten as s = −i(ω/Ω)(Hpk)−2. Therefore, s ∼ 1 indicates strong instability on dynamical time-scales and length-scales ∼Hp. We again note that validity of the equations requires k ≲ 1/H. The determinant (25) gives a third-order dispersion relation:   \begin{eqnarray} &&\!\!\!\!&&{s^3-s^2\left[aQ_1-2Q_2+|\psi |\big (aQ^{\prime }_1-Q^{\prime }_2\big )\right]-s\left[2a{Q}_1{Q}_2-Q_2^2-Q_3^2+|\psi |\big (a{Q}_1{Q}^{\prime }_2-Q_2Q^{\prime }_2-Q_3Q^{\prime }_3\big )\right]}\nonumber\\ &&-a\left[Q_1(Q_2^2+Q_3^2)+|\psi |\big (Q_1Q_2Q^{\prime }_2-Q_1^{\prime }Q_2^2+Q_1Q_3Q^{\prime }_3-Q_1^{\prime }Q_3^2\big ) \right]=0. \end{eqnarray} (28)The disc becomes unstable if any of the roots of (28) has a positive real part, i.e. ℜ(s) > 0, as the perturbations then grow exponentially with time. For completeness, we note, following the discussion of the thermodynamics in Section 2.2, that the equivalent of (28) for an analysis explicitly including the energy equation is given by (cf. Ogilvie 2000)   \begin{eqnarray} &&{s^3-s^2\left[apQ_1-2Q_2+|\psi |\big (aQ^{\prime }_1-Q^{\prime }_2\Big )\right]-s\left[2ap{Q}_1{Q}_2-Q_2^2-Q_3^2+|\psi |\big (ap{Q}_1{Q}^{\prime }_2-Q_2Q^{\prime }_2-Q_3Q^{\prime }_3 -a(p-1)Q^{\prime }_1Q_2\big )\right]}\nonumber\\ &&-ap\left[Q_1(Q_2^2+Q_3^2)+|\psi |\big (Q_1Q_2Q^{\prime }_2-Q_1^{\prime }Q_2^2+Q_1Q_3Q^{\prime }_3-Q_1^{\prime }Q_3^2\big ) \right]=0, \end{eqnarray} (29)where p is the exponent of Σ in νΣ (or equivalently the exponent of Σ in $$\mathcal {I}$$). In a planar disc, the usual thermal-viscous stability criterion (Lightman & Eardley 1974) is given by p > 0. For the isothermal discs we consider here, we have p = 1 (for which 28 and 29 are identical). For Thomson scattering opacity, p = 5/3, and for Kramers opacity, p = 10/7, and as we have noted above, the coefficients do not vary strongly between these cases. Therefore, we expect the calculations we present here (with p = 1 and coefficients derived from an isothermal equation of state) to hold at least qualitatively for a wide range of thermodynamic treatments. 3.2 Instability criterion and the growth rates Before arriving at the full criterion (provided below), we first give the roots of (28) for some simplified cases. We will see below that these cases are generally useful for both understanding the more complex cases, and as close approximations to the full solution. We label the roots as s(Qi), where the subscript i denotes the torque coefficient(s) that we include in the solution. First, if we set all Qi terms in the dispersion relation to zero apart from Q1 (i.e. Q2 = Q3 = 0), we find one root which is given as   \begin{eqnarray} s(Q_1)= a\big (Q_1 + Q^{\prime }_1|\psi |\big ) = -a\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(-Q_1\left|\psi \right|\right). \end{eqnarray} (30)For instability, one requires ℜ[s(Q1)] > 0. For |ψ| = 0, this corresponds to the usual criterion for viscous instability in a planar disc (given by ∂(νΣ)/∂Σ < 0; Lightman & Eardley 1974, although we note that the relevant term is simplified here as we have assumed a simple equation of state, cf. equation 83 of Ogilvie 2000). However, in the warped case, the criterion is modified with an extra term determined by the behaviour of the coefficient Q1 with warp amplitude. Fig. 1 (left-hand panel) shows the variation of the dimensionless growth rates of this mode, ℜ[s(Q1)] with |ψ| for α = 0.01, 0.03, and 0.1. The disc becomes unstable after some critical warp amplitude, |ψ|c. The critical warp amplitude required for instability is very small for discs with low α. The growth rates are a modest fraction ≳0.1 of dynamical. As a result, discs can easily become unstable for this case. Figure 1. View largeDownload slide Left-hand panel: This shows the dimensionless growth rate ℜ[s(Q1)] (equation 30) as a function of dimensionless disc warp amplitude |ψ| caused by the term Q1 in the equations, which governs radial angular momentum transfer, for three different values of the viscosity parameter α. Instability occurs when ℜ[s] > 0. We note that the maximal growth rates ≈ℜ[s]Ω ∼ 0.3Ω are comparable with dynamical and can occur for small values of |ψ|, and remain unstable for large values of |ψ|. Right-hand panel: This shows the dimensionless growth rate ℜ[s(Q2)] caused by the action of the term Q2, which governs the viscous flattening of the warped disc. We note that this term shows no instability for large values of α. For smaller values of α, there is a range of disc warp amplitude |ψ| for which instability occurs. When it does occur, the time-scale for instability can be short with growth rates ≈ℜ[s]Ω ∼ Ω being dynamical. In both panels, the grey line represents zero growth rate. We remind the reader that these illustrative calculations are based on the artificial assumption that two of the three Qi are set to zero. Figure 1. View largeDownload slide Left-hand panel: This shows the dimensionless growth rate ℜ[s(Q1)] (equation 30) as a function of dimensionless disc warp amplitude |ψ| caused by the term Q1 in the equations, which governs radial angular momentum transfer, for three different values of the viscosity parameter α. Instability occurs when ℜ[s] > 0. We note that the maximal growth rates ≈ℜ[s]Ω ∼ 0.3Ω are comparable with dynamical and can occur for small values of |ψ|, and remain unstable for large values of |ψ|. Right-hand panel: This shows the dimensionless growth rate ℜ[s(Q2)] caused by the action of the term Q2, which governs the viscous flattening of the warped disc. We note that this term shows no instability for large values of α. For smaller values of α, there is a range of disc warp amplitude |ψ| for which instability occurs. When it does occur, the time-scale for instability can be short with growth rates ≈ℜ[s]Ω ∼ Ω being dynamical. In both panels, the grey line represents zero growth rate. We remind the reader that these illustrative calculations are based on the artificial assumption that two of the three Qi are set to zero. When we take into account only Q2 (i.e. Q1 = Q3 = 0), then the solution of (28) gives   \begin{eqnarray} s(Q_2)=-\big (Q_2+Q^{\prime }_2|\psi |\big ) = -\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_2\left|\psi \right|\right). \end{eqnarray} (31)The disc becomes unstable when ℜ[s(Q2)] > 0. This condition is similar to (30), but is caused by instability in the warp diffusion coefficient (Q2) and its effect on the warp amplitude. This is distinct from the usual viscous instability as it acts on the warp amplitude rather than the disc surface density. This mode becomes unstable at some critical warp amplitude interval. Fig. 1 (right-hand panel) shows the dimensionless growth rates found for this mode. We see that for small α, the growth rates of this instability are significantly higher than those found in the first simplified case where we included only the Q1 term. However, for large α, there is no instability from s(Q2). We note that there is a second root in this case, which is equal to −Q2, but this root is never unstable within the α formalism. We should also note that in these simplified cases, s is real. Thus, the dimensionless growth rates are equal to s. Below we will encounter complex values of s in which case the dimensionless growth rate is given by ℜ(s). If this analysis is repeated for (29) as opposed to (28), then (30) is modified to $$s = apQ_1 + a\left|\psi \right|Q^{\prime }_1$$, and thus, the Lightman & Eardley (1974) criterion for thermal-viscous instability is modified from apQ1 > 0 to $$apQ_1 > -a\left|\psi \right|Q^{\prime }_1$$. In contrast, (31) remains unchanged. Now we include both Q1 and Q2 terms (i.e. Q3 = 0) to find a more general criterion for instability. In this case, (28) gives three roots, one of which corresponds to −Q2 again, and is thus stable for any |ψ|. The remaining roots of the cubic correspond to   \begin{eqnarray} s_{\pm }(Q_1,Q_2)&=&\frac{1}{2}\bigg [a\big (Q_1 + Q^{\prime }_1|\psi |\big )-\big (Q_2 +Q^{\prime }_2 |\psi |\big )\nonumber\\ &&\pm \,\sqrt{\left[a \big (Q_1 + Q^{\prime }_1 |\psi |\big )-\big (Q_2 + Q^{\prime }_2 |\psi |\big ) \right]^2 +4a\left[\left(Q_1Q_2 + (Q_1Q^{\prime }_2 -Q^{\prime }_1 Q_2)|\psi |\right)\right]}\bigg ]. \end{eqnarray} (32)We see that these modes include both instabilities found in the simplified cases above. Inspection of (32) reveals that ℜ(s+) ≥ ℜ(s−). Therefore, when the disc is unstable, s+ represents the mode that shows more rapid growth. Fig. 2 shows the variation of ℜ[s+(Q1, Q2)] with |ψ|. The disc becomes unstable when ℜ[s+(Q1, Q2)] > 0. We plot ℜ[s(Q1)] and ℜ[s(Q2)] in the same graph for comparison. ℜ[s+(Q1, Q2)] shows similar behaviour to ℜ[s(Q2)] given by (31) for |ψ| ≲ 2.5, and shows similar behaviour to ℜ[s(Q1)] given by (30) for |ψ| ≳ 2.5. The disc becomes unstable at some critical warp amplitude, |ψ|c. For some values of α, the disc becomes stable again at larger warp amplitude before becoming unstable again at a second critical warp amplitude. The critical warp amplitudes and the maximum values of the dimensionless growth rates are almost equivalent for s+(Q1, Q2) and s(Q2). The Q1 term tends to become dominant only at high warp amplitudes. We should note that ℜ[s−(Q1, Q2)] also becomes positive at some critical warp amplitude, but is very close to zero. As a result, both modes, i.e. s+(Q1, Q2) and s−(Q1, Q2), are unstable after the critical warp amplitude, but s+(Q1, Q2) is always the most unstable. Figure 2. View largeDownload slide Here, we combine the two panels of Fig. 1 to illustrate the combined effects of the two terms Q1 and Q2. The curves show the dimensionless growth rates ℜ[s+(Q1, Q2)] (black solid line) with |ψ| (equation 32), compared against ℜ[s(Q1)] (green solid line) and ℜ[s(Q2)] (red dashed line). The grey line represents zero growth rate. For small α = 0.01 (left-hand panel), we see that instability occurs for all warp amplitudes |ψ| > 0.2, except for a very narrow region around |ψ| ≈ 2.96. The growth rate is generally large, except around |ψ| ≈ 3.0. The largest growth rates occur between 0.2 < |ψ| ≲ 2.5, where the Q2 term dominates. For larger |ψ|, the Q1 term dominates and the growth rates are sub-dynamical. For intermediate α = 0.03 (right-hand panel), the same applies except that the growth rates becomes smaller, and there is a wider band of stability around |ψ| = 2.2. For large |ψ|, the growth rates in the two cases (α = 0.01 and 0.03) are similar. Figure 2. View largeDownload slide Here, we combine the two panels of Fig. 1 to illustrate the combined effects of the two terms Q1 and Q2. The curves show the dimensionless growth rates ℜ[s+(Q1, Q2)] (black solid line) with |ψ| (equation 32), compared against ℜ[s(Q1)] (green solid line) and ℜ[s(Q2)] (red dashed line). The grey line represents zero growth rate. For small α = 0.01 (left-hand panel), we see that instability occurs for all warp amplitudes |ψ| > 0.2, except for a very narrow region around |ψ| ≈ 2.96. The growth rate is generally large, except around |ψ| ≈ 3.0. The largest growth rates occur between 0.2 < |ψ| ≲ 2.5, where the Q2 term dominates. For larger |ψ|, the Q1 term dominates and the growth rates are sub-dynamical. For intermediate α = 0.03 (right-hand panel), the same applies except that the growth rates becomes smaller, and there is a wider band of stability around |ψ| = 2.2. For large |ψ|, the growth rates in the two cases (α = 0.01 and 0.03) are similar. Finally, if we include all the Qi terms in the dispersion relation, we obtain the full criterion for instability as ℜ[s1(Q1, Q2, Q3)] > 0 with   \begin{eqnarray} s_{1}(Q_1,Q_2,Q_3)= \frac{1}{6}\big [2\mathcal {C}_1+2^{2/3}\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{1/3}+ 2^{4/3}\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{-1/3}\mathcal {C}_4\big ], \end{eqnarray} (33)where   \begin{eqnarray} \mathcal {C}_1 &=&\ \ a \big (Q_1+Q_1^{\prime } |\psi | \big )-2 \big (Q_2+Q_2^{\prime } |\psi |\big ), \end{eqnarray} (34)  \begin{eqnarray} \mathcal {C}_2 &=&\ \ 3 a \bigg [Q_1^{\prime } |\psi | \bigg (2 Q_2^{\prime 2} |\psi |^2 + 5 Q_2 Q_2^{\prime } |\psi | - 4 Q_2^2 - 3 Q_3 \Big (Q_3^{\prime } |\psi | + 4 Q_3\Big )\bigg ) + Q_1 \bigg (-Q_2^{\prime 2} |\psi |^2 + 2 Q_2 Q_2^{\prime } |\psi |\nonumber\\ &&+\, 2 \big (Q_2^2 + 3 Q_3 \big (Q_3^{\prime } |\psi | + Q_3 \big )\big )\bigg )\bigg ]+2 a^3 \big (Q_1 +Q_1^{\prime } |\psi |\big )^3 + 3 a^2 \big (Q_1^2 -2 Q_1^{\prime 2} |\psi |^2 - Q_1 Q_1^{\prime } |\psi | \big ) \big (2 Q_2 + Q_2^{\prime } |\psi | \big )\nonumber\\ &&-\, \big (2 Q_2+Q_2^{\prime } |\psi |\big ) \big (2 Q_2^{\prime 2} |\psi |^2 - Q_2 Q_2^{\prime } |\psi | - Q_2^2 - 9 Q_3 \big (Q_3^{\prime } |\psi | + Q_3 \big )\big ), \end{eqnarray} (35)  \begin{eqnarray} \mathcal {C}_3&=&\ \ 4 \bigg (3 \big (Q_2^2+Q_2^{\prime } Q_2|\psi |-a Q_1 (2 Q_2+Q_2^{\prime } |\psi |)+Q_3 (Q_3+Q_3^{\prime } |\psi |)\Big )-\big (2 Q_2+Q_2^{\prime } |\psi | -a (Q_1+Q_1^{\prime } |\psi |)\big )^2\bigg )^3\nonumber\\ &&+\,\Big [2 a^3 (Q_1+Q_1^{\prime } |\psi |)^3+(2 Q_2+Q_2^{\prime } |\psi |) \left(Q_2^2+Q_2^{\prime } Q_2|\psi |+9 Q_3^2-2 Q_2^{\prime 2} |\psi | ^2+9 Q_3^{\prime } Q_3 |\psi | \right)\nonumber\\ &&+\,3 a \bigg (Q_1 \left(2 Q_2^2+2 Q_2^{\prime } Q_2|\psi |+6 Q_3^2-Q_2^{\prime 2} |\psi | ^2+6 Q_3^{\prime } Q_3 |\psi | \right)+Q_1^{\prime } |\psi | \big (-4 Q_2^2+5 Q_2^{\prime } Q_2|\psi | +2 Q_2^{\prime 2} |\psi | ^2\nonumber\\ &&-\,12 Q_3^2 -3 Q_3^{\prime } Q_3 |\psi | \big )\bigg )+3 a^2 (2 Q_2+Q_2^{\prime } |\psi |) \left(Q_1^2-Q_1^{\prime } Q_1|\psi |-2 Q_1^{\prime 2} |\psi | ^2\right)\big ]^2,\ \rm and \end{eqnarray} (36)  \begin{eqnarray} \mathcal {C}_4&= &\ \ a^2 \big (Q_1+Q_1^{\prime } |\psi | \big )^2+ Q_2^2+ Q_2^{\prime } Q_2|\psi | +Q_2^{\prime 2} |\psi | ^2 +a \big (2 Q_2+Q_2^{\prime } |\psi | \big ) \big (Q_1-2 Q_1^{\prime } |\psi | \big )-3 \big (Q_3^2 + Q_3^{\prime } Q_3 |\psi |\big ). \end{eqnarray} (37)The other two roots of the dispersion relation are given by   \begin{eqnarray} s_{2}(Q_1,Q_2,Q_3)= \frac{1}{6}\big [2\mathcal {C}_1-2^{-1/3}\big (1-\mathrm{i}\sqrt{3}\big )\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{1/3}-2^{1/3}\big (1+\mathrm{i}\sqrt{3}\big )\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{-1/3}\mathcal {C}_4\big ] \end{eqnarray} (38)and   \begin{eqnarray} s_{3}(Q_1,Q_2,Q_3)= \frac{1}{6}\big [2\mathcal {C}_1-2^{-1/3}\big (1+\mathrm{i}\sqrt{3}\big )\big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{1/3}-2^{1/3}\big (1-\mathrm{i}\sqrt{3}\big )\Big (\mathcal {C}_2+\mathcal {C}_3^{1/2}\big )^{-1/3}\mathcal {C}_4\big ]. \end{eqnarray} (39) These roots are either all real, or there can be one real root and two complex conjugates. For determining instability, we are interested in the real roots and the real parts of the complex roots. While there is always one real root, the other two roots can swap between real or complex conjugates for different values of |ψ|. One root is always stable. Of the other two roots, one always gives the fastest growth (i.e. when the real part is positive it shows the largest value) and this root is given by (33). For some values of |ψ|, this root is a complex conjugate and thus the other root has the same growth rate. In Fig. 3, we show the dimensionless growth rates of the instability found from the full criterion (33) for different values of α. We see that the disc becomes unstable for sufficiently large |ψ| values. The critical warp amplitudes where the disc becomes unstable are smaller and the growth rates of the instability are higher for low α. In Fig. 4, we plot the maximum values of the dimensionless growth rates of instability (ℜ[s]max) against α. We see that ℜ[s]max values are inversely proportional to α. Finally, the critical warp amplitudes for instability are plotted against α in Fig. 5. This figure shows the stable and unstable regions in the (α, |ψ|) parameter space. Figure 3. View largeDownload slide Here, we show the dimensionless growth rates ℜ[s] as functions of |ψ| taking account of all three terms, Q1, Q2, and Q3, for different values of α. The effect of the Q3 term is negligible for the parameters we have considered (see Figs 6 and 7). Figure 3. View largeDownload slide Here, we show the dimensionless growth rates ℜ[s] as functions of |ψ| taking account of all three terms, Q1, Q2, and Q3, for different values of α. The effect of the Q3 term is negligible for the parameters we have considered (see Figs 6 and 7). Figure 4. View largeDownload slide This plot shows the maximum growth rate of the instability as a function of viscosity parameter α. Note that for different values of α the maximum growth occurs for different values of the dimensionless warp amplitude |ψ|. This plot shows that instability occurs for all values of α in this range. For small values of α, the maximum growth rate is dynamical. For larger values of α, the maximum growth rate is smaller, but still comparable to dynamical (∼0.2Ω). Figure 4. View largeDownload slide This plot shows the maximum growth rate of the instability as a function of viscosity parameter α. Note that for different values of α the maximum growth occurs for different values of the dimensionless warp amplitude |ψ|. This plot shows that instability occurs for all values of α in this range. For small values of α, the maximum growth rate is dynamical. For larger values of α, the maximum growth rate is smaller, but still comparable to dynamical (∼0.2Ω). Figure 5. View largeDownload slide Stable (white) and unstable (green) regions in the (α, |ψ|) parameter space. This plot represents the critical warp amplitudes for instability to occur in discs with various α values. We note that a flat disc is always stable (assuming an isothermal equation of state so that there is no Lightman–Eardley-like instability). For 0.01 ≤ α ≤ 0.2, a nearly flat disc with |ψ| ≲ 0.1 is also always stable. For a given value of α, there is always a minimum value of the warp amplitude, which gives rise to instability. Figure 5. View largeDownload slide Stable (white) and unstable (green) regions in the (α, |ψ|) parameter space. This plot represents the critical warp amplitudes for instability to occur in discs with various α values. We note that a flat disc is always stable (assuming an isothermal equation of state so that there is no Lightman–Eardley-like instability). For 0.01 ≤ α ≤ 0.2, a nearly flat disc with |ψ| ≲ 0.1 is also always stable. For a given value of α, there is always a minimum value of the warp amplitude, which gives rise to instability. If we compare (32) and (33), we see that the inclusion of Q3 produces a big difference in the formula, but does not have a strong impact on the numerical results. In Fig. 6, we compare two criteria, (32), where Q3 = 0, and (33), where Q3 ≠ 0. We see that the inclusion of Q3 has a very small effect on the critical warp amplitudes and the growth rates. We also plot the ratios of the |ψ|c and ℜ[s]max found for the two different cases, s(Q1, Q2) and s(Q1, Q2, Q3). Fig. 7 shows the values of X and Y as a function of α, where X = |ψ|cℜ[s(Q1, Q2, Q3)]/|ψ|cℜ[s(Q1, Q2)] and Y = ℜ[smax(Q1, Q2, Q3)]/ℜ[smax(Q1, Q2)]. The solutions are in excellent agreement. Therefore, (32) can be used to decide whether the disc is stable or unstable, at least for the parameters we have explored. Figure 6. View largeDownload slide Comparison of the full (33) and simplified criteria (32) for instability for α = 0.01 and 0.1. The grey line represents zero growth rate. These plots show that the simplified criteria are an extremely good approximation to both the location of the critical warp amplitude(s) and the dimensionless growth rate as a function of warp amplitude. Figure 6. View largeDownload slide Comparison of the full (33) and simplified criteria (32) for instability for α = 0.01 and 0.1. The grey line represents zero growth rate. These plots show that the simplified criteria are an extremely good approximation to both the location of the critical warp amplitude(s) and the dimensionless growth rate as a function of warp amplitude. Figure 7. View largeDownload slide X = |ψ|cℜ[s(Q1, Q2, Q3)]/|ψ|cℜ[s(Q1, Q2)] and Y = ℜ[smax(Q1, Q2, Q3)]/ℜ[smax(Q1, Q2)] ratios for different α values. The grey line represents unity. Thus, the critical warp amplitudes and the maximum values of the dimensionless growth rates are captured by the simplified root (32) to high accuracy for all of the α values we have explored. Figure 7. View largeDownload slide X = |ψ|cℜ[s(Q1, Q2, Q3)]/|ψ|cℜ[s(Q1, Q2)] and Y = ℜ[smax(Q1, Q2, Q3)]/ℜ[smax(Q1, Q2)] ratios for different α values. The grey line represents unity. Thus, the critical warp amplitudes and the maximum values of the dimensionless growth rates are captured by the simplified root (32) to high accuracy for all of the α values we have explored. In conclusion, the criteria for instability can be expressed from (32) as below:   \begin{eqnarray} {\rm If} \ \ \,\left[a\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_1\left|\psi \right|\right)-\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_2\left|\psi \right|\right)\right]>0,\,\rm the \,\rm disc \,\rm is \,\rm unstable, \end{eqnarray} (40)  \begin{eqnarray} {\rm or \ \rm if} \ \ \ \left[a\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_1\left|\psi \right|\right)-\frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(Q_2\left|\psi \right|\right)\right]<0,\ \ \ \ \ {\rm and} \ \ \ \ \ 4a\left[\left(Q_1Q_2 + (Q_1Q^{\prime }_2 -Q^{\prime }_1 Q_2)|\psi |\right)\right]>0,\,\rm the \,\rm disc \,\rm \,is \,\rm also\, \rm unstable. \end{eqnarray} (41) 4 HEURISTIC PICTURE OF THE VISCOUS-WARP INSTABILITY To elucidate the fundamental physics of this viscous-warp disc instability from the complex analysis described above, we present here a simplified analysis. Starting from (15), we simplify the picture aggressively. In a warped disc with α ≪ 1, the smoothing of the disc warp occurs faster than mass is transported radially (Papaloizou & Pringle 1983). So we assume that the important torque is the diffusive torque relating to ν2. Our exploration of the analysis above confirms this is a sensible approach in many, but not all, cases. Thus,   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} t}\left(\Sigma r^2\Omega \boldsymbol {l}\right) \sim \frac{1}{r}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(Q_2\Sigma c_{\rm s}^2r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right). \end{eqnarray} (42)Interpreting this as a diffusion equation implies that the flux of misaligned angular momentum is   \begin{eqnarray} \boldsymbol {F}_{\rm mis} = Q_2\Sigma c_{\rm s}^2r^3\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r} \end{eqnarray} (43)with rate of transfer $${\sim } \left|\boldsymbol {F}_{\rm mis}/\Sigma r^2 \Omega \right| = (1/2)\nu _2\left|\psi \right|$$. Connecting this with classical instability of a diffusion equation yields the criterion for instability as   \begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} \left|\psi \right|}\left(\nu _2\left|\psi \right|\right) < 0. \end{eqnarray} (44)This can also be trivially recovered by performing a stability analysis on (42), with $$\boldsymbol {l} \rightarrow \boldsymbol {l}+\delta \boldsymbol {l}$$, δQ2 = (∂Q2/∂|ψ|)δ|ψ| and assuming δΣ = 0. The criterion (44) can be interpreted as stating if ν2(|ψ|) is such that the maximum diffusion rate is not located at maxima in warp amplitude, then local maxima in warp amplitude will grow, and the disc will break. As this will result in rapid transfer of mass out of this region due to the large warp amplitude implying large torques, this will also be realized by a significant drop in local surface density. This resembles the Lightman–Eardley viscous instability but for a warped disc with the warp amplitude playing the role of the surface density. Such behaviour is impossible if the effective viscosities are constant with varying warp amplitude. This is confirmed in our analysis above, which yields no unstable solutions if Qi are independent of |ψ|. This is consistent with the results of Nixon & King (2012).4 Finally, we note that when the correct limits are applied, the full stability analysis presented above reproduces the ad hoc answer derived here (see equation 31). 5 APPLICATION TO DISC TEARING In this work, we have focused on disc ‘breaking’, where an isolated already-warped disc (i.e. with no external torque) evolves to create a discontinuity between two or more planes, with a series of low-density orbits connecting them over a small radial extent. We have revealed above that this occurs due to a viscous-warp instability, which occurs when either (40) or (41) is satisfied. In many physical systems, the disc is subject to an external torque often occurring due to precession (e.g. the Lense–Thirring effect). In this case, the warp initially grows with time. If the disc ‘breaks’ under these conditions creating disconnected rings that keep precessing effectively independently, then this is what we have called disc ‘tearing’ in previous papers (e.g. Nixon et al. 2012; Doğan et al. 2015; Nealon et al. 2015). It is of great interest to connect the analysis in this work to the idea of disc tearing. Unfortunately at first appearances, we cannot perform the same analysis as above for disc tearing, where an external forced precession is imposed, as the background unperturbed state is no longer necessarily an equilibrium solution on short time-scales. It may be that this is unimportant, and an approximate solution can be found that is accurate enough. We will check this in future work. For our purposes here, we repeat that the analysis above is a local analysis that takes a given warped disc shape and determines its stability. For disc tearing, we need to consider how the disc gets into a warped configuration, and at what stage during that process instability occurs. Here, we can make a reasoned guess at what the disc tearing criterion is likely to be. Above we have derived, for various system parameters, the critical warp amplitude |ψ|c at which the disc becomes unstable. It would seem reasonable that a forced disc would tear up if there is no |ψ| < |ψ|c at which the torque attempting to flatten the disc (∝ Q2|ψ|) is stronger than the torque induced by, for example, precession. As discussed in Doğan et al. (2015), without the knowledge of the value of |ψ|c, this calculation cannot be done a priori without performing a full numerical simulation. However, our analysis enables exactly that. We will test this in the future by comparing this reasoning with numerical simulations. We note that given the viscous torque and precession torque (in the Lense-Thirring case) have no dependence on inclination angle, this criterion at first appears at odds with previous simulations that show, for example, a strong inclination angle dependence. There is also no dependence on the disc thickness, but previous simulations show that thinner discs are easier to tear. The answer to this is likely to be that the warp amplitude cannot be larger than a value dependent on the inclination angle and H/R. This can be seen by the following argument. First, the warp amplitude is   \begin{eqnarray} \left|\psi \right| = r\left|\frac{\mathrm{\partial} \boldsymbol {l}}{\mathrm{\partial} r}\right|, \end{eqnarray} (45)and the disc unit tilt vector is   \begin{eqnarray} \boldsymbol {l} = \left(\cos \gamma \sin \beta ,\sin \gamma \sin \beta ,\cos \beta \right), \end{eqnarray} (46)where γ is the disc twist and β is the disc tilt (e.g. Nixon 2012). Let us assume γ = 0, and a maximal warp from β to zero (the same result, equation 48, is arrived at if β is assumed a constant and γ varies from 0 to π). Thus,   \begin{eqnarray} \delta \boldsymbol {l} = \left(\sin \beta ,0,\cos \beta -1\right), \end{eqnarray} (47)which has magnitude (2 − 2cos β)1/2. We can also assume that δr cannot become smaller than H, and thus the maximum warp amplitude is given by   \begin{eqnarray} \left|\psi \right|_{\rm max} \sim \left(2-2\cos \beta \right)^{1/2}\frac{R}{H} \sim \beta \frac{R}{H}. \end{eqnarray} (48)Therefore, discs that are already close to alignment, or are thick, are difficult to tear, as the maximum warp amplitude is small and potentially less than the critical warp amplitude required for instability. 6 CONCLUSIONS We have derived the criterion for isolated warped discs to break. This occurs physically due to viscous anti-diffusion of the warp amplitude. The resulting large torques that are induced lead to the disc surface density lowering in the unstable regions. Thus, the instability leads to a succession of orbits with planes that change rapidly with radius, with only a small amount of mass present. This is a local instability so it is not likely to directly affect global disc properties such as the central accretion rate. However, this instability underlies the process of disc tearing which has the capacity to dramatically alter the instantaneous accretion rate and the observable properties of the disc on short time-scales. The instability is governed by a complicated criterion (33), but we have shown that the criterion can be simplified in most cases to a more manageable criteria given by (40) and (41). This simply implies that the diffusion of the warp must be maximal at the location of the maximum warp amplitude otherwise the warp will grow until the disc disconnects on scales of order the disc thickness ∼H. In future works, we will perform this analysis for non-Keplerian discs, connect our results with numerical simulations, and include a more sophisticated thermal treatment. For discs with moderate values of α ∼ 0.01–0.1 and reasonably thin H/R ≲ α, our work supports the idea of discs being viscously unstable to disconnecting into discrete planes. This process underpins the state-transition cycle suggested by Nixon & Salvesen (2014). The natural link between the viscous instability in planar discs, which underlies the thermal-viscous outburst cycle observed in many types of binary systems (many of which are thought to contain warped discs), and the viscous instability of warped discs is appealing. This may be a way to introduce hysteresis through a change in the disc stability with warping. We will explore this in future work. ACKNOWLEDGEMENTS We thank the referee for providing a helpful and detailed report. SD gratefully acknowledges the warm hospitality of the Theoretical Astrophysics Group at University of Leicester during her visit. SD is supported by the Turkish Scientific and Technical Research Council (TÜBİTAK - 117F280). CJN is supported by the Science and Technology Facilities Council (grant number ST/M005917/1). The Theoretical Astrophysics Group at the University of Leicester is supported by an STFC Consolidated Grant. Footnotes 1 To date only the SPH code phantom (Price et al. 2017) has been tested for correctly modelling warped discs against the non-linear theory of Ogilvie (1999), finding excellent agreement in the simulated parameter space. 2 Note that both r here and R in (1) represent the local cylindrical radius of the annulus. 3 In a steady disc, ttherm ∼ α−1tdyn (Pringle 1981), and so, in general, ttherm ≳ tdyn with rough equality only for α ≈ 1. 4 Note that in Nixon & King (2012) when using the Qi coefficients from Ogilvie (1999), we smoothed the viscosity coefficients over the neighbouring few rings to avoid what was perceived there as numerical issues. Actually, this behaviour was the physical instability discussed here. REFERENCES Aly H., Dehnen W., Nixon C., King A., 2015, MNRAS , 449, 65 https://doi.org/10.1093/mnras/stv128 CrossRef Search ADS   Bate M. R., Lodato G., Pringle J. E., 2010, MNRAS , 401, 1505 https://doi.org/10.1111/j.1365-2966.2009.15773.x CrossRef Search ADS   Doğan S., Nixon C., King A., Price D. J., 2015, MNRAS , 449, 1251 https://doi.org/10.1093/mnras/stv347 CrossRef Search ADS   Facchini S., Lodato G., Price D. J., 2013, MNRAS , 433, 2142 https://doi.org/10.1093/mnras/stt877 CrossRef Search ADS   King A. R., Ritter H., 1998, MNRAS , 293, L42 https://doi.org/10.1046/j.1365-8711.1998.01295.x CrossRef Search ADS   Krolik J. H., Hawley J. F., 2015, ApJ , 806, 141 https://doi.org/10.1088/0004-637X/806/1/141 CrossRef Search ADS   Lense J., Thirring H., 1918, Phys. Z. , 19, 156 Lightman A. P., Eardley D. M., 1974, ApJ , 187, L1 https://doi.org/10.1086/181377 CrossRef Search ADS   Lodato G., Price D. J., 2010, MNRAS , 405, 1212 Lodato G., Pringle J. E., 2007, MNRAS , 381, 1287 https://doi.org/10.1111/j.1365-2966.2007.12332.x CrossRef Search ADS   Lubow S. H., 1992, ApJ , 398, 525 https://doi.org/10.1086/171877 CrossRef Search ADS   Lubow S. H., Ogilvie G. I., 2000, ApJ , 538, 326 https://doi.org/10.1086/309101 CrossRef Search ADS   Lucas W. E., Bonnell I. A., Davies M. B., Rice W. K. M., 2013, MNRAS , 433, 353 https://doi.org/10.1093/mnras/stt727 CrossRef Search ADS   McKinney J. C., Tchekhovskoy A., Blandford R. D., 2013, Science , 339, 49 https://doi.org/10.1126/science.1230811 CrossRef Search ADS PubMed  Meyer F., Meyer-Hofmeister E., 1982, A&A , 106, 34 Morales Teixeira D., Fragile P. C., Zhuravlev V. V., Ivanov P. B., 2014, ApJ , 796, 103 https://doi.org/10.1088/0004-637X/796/2/103 CrossRef Search ADS   Nealon R., Price D. J., Nixon C. J., 2015, MNRAS , 448, 1526 https://doi.org/10.1093/mnras/stv014 CrossRef Search ADS   Nealon R., Nixon C., Price D. J., King A., 2016, MNRAS , 455, L62 Nixon C. J., 2012, MNRAS , 423, 2597 https://doi.org/10.1111/j.1365-2966.2012.21072.x CrossRef Search ADS   Nixon C. J., King A. R., 2012, MNRAS , 421, 1201 https://doi.org/10.1111/j.1365-2966.2011.20377.x CrossRef Search ADS   Nixon C., King A., 2016, Lecture Notes in Physics, Vol. 905, Astrophysical Black Holes . Springer International Publishing, Switzerland, p. 45 Nixon C., Salvesen G., 2014, MNRAS , 437, 3994 https://doi.org/10.1093/mnras/stt2215 CrossRef Search ADS   Nixon C., King A., Price D., Frank J., 2012, ApJ , 757, L24 https://doi.org/10.1088/2041-8205/757/2/L24 CrossRef Search ADS   Nixon C., King A., Price D., 2013, MNRAS , 434, 1946 https://doi.org/10.1093/mnras/stt1136 CrossRef Search ADS   Ogilvie G. I., 1999, MNRAS , 304, 557 https://doi.org/10.1046/j.1365-8711.1999.02340.x CrossRef Search ADS   Ogilvie G. I., 2000, MNRAS , 317, 607 https://doi.org/10.1046/j.1365-8711.2000.03654.x CrossRef Search ADS   Ogilvie G. I., Latter H. N., 2013, MNRAS , 433, 2403 https://doi.org/10.1093/mnras/stt916 CrossRef Search ADS   Papaloizou J. C. B., Lin D. N. C., 1995, ApJ , 438, 841 https://doi.org/10.1086/175127 CrossRef Search ADS   Papaloizou J. C. B., Pringle J. E., 1983, MNRAS , 202, 1181 https://doi.org/10.1093/mnras/202.4.1181 CrossRef Search ADS   Papaloizou J. C. B., Terquem C., 1995, MNRAS , 274, 987 Price D. J. et al.  , 2017, PASA , preprint (arXiv:1702.03930) Pringle J. E., 1981, ARA&A , 19, 137 CrossRef Search ADS   Pringle J. E., 1992, MNRAS , 258, 811 https://doi.org/10.1093/mnras/258.4.811 CrossRef Search ADS   Pringle J. E., 1996, MNRAS , 281, 357 https://doi.org/10.1093/mnras/281.1.357 CrossRef Search ADS   Pringle J. E., 1997, MNRAS , 292, 136 https://doi.org/10.1093/mnras/292.1.136 CrossRef Search ADS   Schandl S., Meyer F., 1994, A&A , 289, 149 Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Simon J. B., Beckwith K., Armitage P. J., 2012, MNRAS , 422, 2685 https://doi.org/10.1111/j.1365-2966.2012.20835.x CrossRef Search ADS   Sorathia K. A., Krolik J. H., Hawley J. F., 2013a, ApJ , 768, 133 https://doi.org/10.1088/0004-637X/768/2/133 CrossRef Search ADS   Sorathia K. A., Krolik J. H., Hawley J. F., 2013b, ApJ , 777, 21 https://doi.org/10.1088/0004-637X/777/1/21 CrossRef Search ADS   Zhuravlev V. V., Ivanov P. B., Fragile P. C., Morales Teixeira D., 2014, ApJ , 796, 104 https://doi.org/10.1088/0004-637X/796/2/104 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial