Local numerical simulations of warped discs

TL;DR

This study uses numerical simulations to analyze the stability and nonlinear behavior of warped astrophysical discs, revealing how warp amplitudes and viscosity influence angular momentum flux and disc evolution timescales.

Contribution

It provides new insights into the nonlinear saturated states of warped discs and how they affect warp evolution, especially at large amplitudes and low viscosities.

Findings

Linear growth rates are due to parametric resonance of inertial waves.

Nonlinear saturated states can reduce laminar flows, affecting warp evolution.

Reversal of angular momentum flux indicates anti-diffusive mass distribution.

Abstract

We study the hydrodynamical stability of the laminar flows associated with warped astrophysical discs using numerical simulations of warped shearing boxes. We recover linear growth rates reported previously due to a parametric resonance of inertial waves, and show that the nonlinear saturated state can significantly reduce the laminar flows, meaning that the warp would evolve on much longer time scales than would be concluded from the internal torques due to these laminar flows. Towards larger warp amplitudes, we find first of all a reversal of angular momentum flux, indicating that the mass distribution would evolve in an anti-diffusive manner, and second that the linear growth rates disappear, possibly because of the very strong shear in the laminar flows in this regime. For discs with small enough viscosity, a nonlinear state can still be found when linear growth rates are absent by introducing a large enough perturbation, either by starting from a nonlinear state obtained at smaller warp amplitude, or by starting from a state with no laminar flows.

Equations

$${\partial }_t\rho +{\partial }_x(\rho v_x)+{\partial }_y[\rho (v_y+q\tau v_x)]+$$

$${\partial }_z[\rho (v_z+|\psi |\cos \tau v_x)]=0,$$

$${\partial }_t{v}_x+v_x{\partial }_x{v}_x+(v_y+q\tau v_x){\partial }_y{v}_x+$$

$$(v_z+|\psi |\cos \tau v_x){\partial }_z{v}_x+$$

$$({\partial }_{x}p+q\tau {\partial }_{y}p+|\psi |\cos \tau {\partial }_{z}p)/\rho =2\Omega v_y,$$

$${\partial }_t{v}_y+v_x{\partial }_x{v}_y+(v_y+q\tau v_x){\partial }_y{v}_y+$$

$$(v_z+|\psi |\cos \tau v_x){\partial }_z{v}_y+$$

$${\partial }_{y}p/\rho =(q-2)\Omega v_x,$$

$${\partial }_t{v}_z+v_x{\partial }_x{v}_z+(v_y+q\tau v_x){\partial }_y{v}_z+$$

$$(v_z+|\psi |\cos \tau v_x){\partial }_z{v}_z+$$

$${\partial }_{z}p/\rho =-|\psi |\Omega \sin \tau v_x-{\Omega }^{2}z.$$

$$v_x$$

$$v_y$$

$$v_z$$

$$p=c_{\mathrm{s}}^2\rho ,$$

$$\rho =\frac{\Sigma }{\sqrt{2\pi }}\exp \left(-\frac{{\Omega }^2{z}^2}{2c_{\mathrm{s}}^2}\right),$$

$$v_x(z,t)$$

$$v_y(z,t)$$

$$v_z(z,t)$$

$$h(z,t)$$

$$Q_1$$

$$Q_2$$

$$Q_3$$

$${\partial }_t\rho +{\partial }_x(\rho v_x)+{\partial }_y[\rho (v_y+q\tau v_x)]+$$

$${\partial }_z[\rho (v_z+|\psi |\cos \tau v_x)]=0,$$

$${\partial }_t(\rho v_x)+{\partial }_x(\rho v_x^2+p)+{\partial }_y[\rho v_x(v_y+q\tau v_x)+q\tau p]+$$

$${\partial }_z[\rho v_x(v_z+|\psi |\cos \tau v_x)+|\psi |\cos \tau p]=2\Omega \rho v_y,$$

$${\partial }_t(\rho v_y)+{\partial }_x(\rho v_x{v}_y)+{\partial }_y[\rho v_y(v_y+q\tau v_x)+p]+$$

$${\partial }_z[\rho v_y(v_z+|\psi |\cos \tau v_x)]=(q-2)\Omega \rho v_x,$$

$${\partial }_t(\rho v_z)+{\partial }_x(\rho v_x{v}_z)+{\partial }_y[\rho v_z(v_y+q\tau v_x)]+$$

$${\partial }_z[\rho v_z(v_z+|\psi |\cos \tau v_x)+p]=-|\psi |\Omega \sin \tau \rho v_x-{\Omega }^{2}z\rho .$$

$${\partial }_t\mathbf{W}+{\partial }_x{\mathbf{F}}_x+{\partial }_y{\mathbf{F}}_y+{\partial }_z{\mathbf{F}}_z=\mathbf{S},$$

$${\partial }_t\rho +{\partial }_x(\rho v_x)=0,$$

$${\partial }_t(\rho v_x)+{\partial }_x(\rho v_x^2+p)=2\Omega \rho v_y,$$

$${\partial }_t(\rho v_y)+{\partial }_x(\rho v_x{v}_y)=(q-2)\Omega \rho v_x,$$

$${\partial }_t(\rho v_z)+{\partial }_x(\rho v_x{v}_z)=-|\psi |\Omega \sin \tau \rho v_x.$$

$${\partial }_t\rho +{\partial }_y[\rho (v_y+q\tau v_x)]=0,$$

$${\partial }_t(\rho v_x)+{\partial }_y[\rho v_x(v_y+q\tau v_x)+q\tau p]=0,$$

$${\partial }_t(\rho v_y)+{\partial }_y[\rho v_y(v_y+q\tau v_x)+p]=0,$$

$${\partial }_t(\rho v_z)+{\partial }_y[\rho v_z(v_y+q\tau v_x)]=0.$$

$${\partial }_t\rho +{\partial }_z[\rho (v_z+|\psi |\cos \tau v_x)]=0,$$

$${\partial }_t(\rho v_x)+{\partial }_z[\rho v_x(v_z+|\psi |\cos \tau v_x)+|\psi |\cos \tau p]=0,$$

$${\partial }_t(\rho v_y)+{\partial }_z[\rho v_y(v_z+|\psi |\cos \tau v_x)]=0,$$

$${\partial }_t(\rho v_z)+{\partial }_z[\rho v_z(v_z+|\psi |\cos \tau v_x)+p]=-{\Omega }^{2}z\rho .$$

$${\lambda }_1$$

$${\lambda }_2$$

$${\lambda }_3$$

$${\lambda }_4$$

Full text

††pubyear: 2018††pagerange: Local numerical simulations of warped discs–References

1 Introduction

Discs are ubiquitous in astrophysics and can be found for example around young stars, supermassive black holes, in cataclysmic variables and around planets, most famously Saturn. While in the simplest picture a disc has a single orbital plane, the orbital plane can vary with radius if there is some form of misalignment present in the system. In that case, the disc is said to be warped.

The misalignment necessary to induce a warp in a disc can for example come from accreted material that has a different angular momentum vector (Bate et al., 2010), a distant companion on an inclined orbit (Papaloizou & Terquem, 1995) or a misaligned magnetic field (Lai, 1999). Warps can also be excited through radiation (Pringle, 1996), winds (Schandl & Meyer, 1994) or tidal effects (Lubow, 1992). The classic case of misalignment is that of an accretion disc for which the orbital axis does not coincide with the spin axis of the central black hole (Bardeen & Petterson, 1975). Observationally, we know of warped discs in X-ray binaries (e.g. Katz, 1973; Kotze & Charles, 2012), and around black holes in galaxy centres (e.g. Miyoshi et al., 1995; Greenhill, 2005). Understanding the evolution of such warped discs due to internal and external torques is a rich and intricate problem.

The first consistent linear theory of the evolution of a warped viscous Keplerian disc was provided by Papaloizou & Pringle (1983), who improved upon earlier work (e.g. Bardeen & Petterson, 1975; Petterson, 1977, 1978; Hatchett et al., 1981) by including the correct form of internal torques necessary to conserve angular momentum. Parametrising the viscosity through the $\alpha$-parameter (Shakura & Sunyaev, 1973), they found that the warp diffuses on a time scale that is a factor $\alpha^{2}$ shorter than the viscous time scale in the regime $H/r<\alpha\ll 1$ (the diffusive regime), where $H$ is the disc angular semithickness. A transition to wavelike propagation of the warp occurs when $\alpha<H/r$ (Papaloizou & Lin, 1995). A fully nonlinear theory for the diffusive regime in Keplerian discs (and for bending waves in non-Keplerian discs) was developed in Ogilvie (1999). This analysis also provided a means to calculate the internal torques governing the evolution of the warp from the amplitude of the warp and other parameters.

The internal fluid motions in a warped disc, driven by a horizontal pressure gradient stemming from the change of the orbital plane with radius, are usually treated as laminar flows. It has been noted that these oscillatory shear flows could in fact be hydrodynamically unstable (Papaloizou & Terquem, 1995; Gammie et al., 2000). The resulting state may be fully turbulent but will at least involve significant wave activity altering the laminar flows and therefore the internal torques (Ogilvie & Latter, 2013b). This means one has to rely on numerical simulations to determine these torques and therefore the evolution of warped discs.

Global numerical hydrodynamical simulations of warped accretion discs are computationally demanding for several reasons. As is usual for astrophysical discs, the range of time scales and length scales is substantial. For warped discs in particular, the calculations are in addition necessarily three-dimensional. Moreover, the vertical extent of the computational domain has to be significantly larger compared to unwarped discs to allow the orbital plane to vary with radius. This leads to regions in the domain of extremely low density, which poses a challenge to most grid-based codes. Most simulations of warped discs have been done using smoothed particle hydrodynamics (SPH) (e.g. Nelson & Papaloizou, 1999; Lodato & Pringle, 2007; Lodato & Price, 2010; Nixon et al., 2012; Xiang-Gruess & Papaloizou, 2013; Nealon et al., 2016). Grid-based global calculations were done by e.g. Fragile et al. (2007); Fragner & Nelson (2010); Krolik & Hawley (2015).

In none of the global simulations mentioned above was a hydrodynamic instability seen, perhaps due to a lack of resolution. A local model of a warp, harnessing the power of the shearing box formalism (e.g. Goldreich & Lynden-Bell, 1965) to deal with warped astrophysical discs, was developed in Ogilvie & Latter (2013a). The warped shearing box has been used recently by Paris & Ogilvie (2018) to examine the effect of a magnetic field on the internal dynamics of a warped disc and on the propagation of a warp. An analogous model has also been developed by Ogilvie & Barker (2014) for eccentric discs. Barker & Ogilvie (2014) used this model to analyse the parametric instability of inertial waves in an eccentric disc. More recently, Wienkers & Ogilvie (2018) have followed the nonlinear evolution of the instability in 2D numerical simulations in an eccentric shearing box. Breaking of the waves away from the midplane was found to limit the growth of unstable modes, and a competitive dynamics ensued between the inertial waves and zonal flows that developed in the box. While the resulting Reynolds stresses were not very effective in transporting angular momentum, they would have the effect of damping the eccentricity.

In the local model that is the warped shearing box it is easier to obtain high resolution, while the potentially unstable laminar flows are still present. The hydrodynamic laminar flows were computed in Ogilvie & Latter (2013a), while their stability was investigated in Ogilvie & Latter (2013b). They found widespread hydrodynamic instability where a pair of inertial waves can grow exponentially by coupling to the laminar warp motion (see also Gammie et al., 2000). As it is the nonlinear outcome of this instability that in the end will determine the internal torques, numerical simulations are needed to calculate these. This is the subject of this paper.

We start in section 2 with the governing equations and the definitions of the components of the internal torque that we are aiming to measure. In section 3 we present the numerical method, which is tested in section 4. Results are presented in section 5 and we conclude in section 6.

2 Basic equations

We will be using the local, isothermal model of a warped disc as derived in Ogilvie & Latter (2013a). Starting with the standard shearing box equations, a transformation is made to warped shearing coordinates that follow the warped orbital motion111In Ogilvie & Latter (2013a) the warped shearing coordinates are denoted with primes. In this paper, we work exclusively with warped shearing coordinates and therefore we omit the primes for clarity.. In this model, the warp is represented by a dimensionless amplitude $|\psi|$ and the governing equations read (equations (38)-(40) and (42) in Ogilvie & Latter, 2013a):

[TABLE]

Here, $q$ is the dimensionless rate of orbital shear ($q=1.5$ in a Keplerian disc), $\Omega$ the orbital angular velocity of the box, $\tau=\Omega t$ the orbital phase, $\rho$ is the fluid density, ${\bf v}=(v_{x},v_{y},v_{z})^{T}$ is the velocity relative to the warped orbital motion. It is related to the standard velocity in the rotating frame of the shearing box ${\bf u}$ by:

[TABLE]

Note that throughout this paper the components of the velocity vector (and, later, of the viscous stress tensor) are referred to the standard orthonormal Cartesian basis.

We use an isothermal equation of state relating pressure and density:

[TABLE]

with $c_{\mathrm{s}}$ the isothermal sound speed. For brevity, we have omitted viscous forces. The three dimensionless numbers characterising a particular physical setup are the warp amplitude $|\psi|$, the shear parameter $q$ and the (shear) viscosity $\alpha$ (see equation (58)). Note that the aspect ratio $H/r=c_{\mathrm{s}}/(r\Omega)$ is not one of the dimensionless parameters of the local model.

In the absence of a warp ($|\psi|=0$), $v_{x}=v_{y}=v_{z}=0$ is a solution, with a density profile

[TABLE]

where $\Sigma$ is the surface density. A warp introduces a laminar flow in the disc, which in the local model only depends on $z$ and $\tau$ (Ogilvie & Latter, 2013a). A nonlinear separation of variables:

[TABLE]

where $h=c_{\mathrm{s}}^{2}\log\rho+\mathrm{constant}$ is the pseudo-enthalpy, leads to a set of ordinary differential equations for the dimensionless amplitudes $u$, $v$, $w$, $f$ and $g$ (Ogilvie & Latter, 2013a). These amplitudes are periodic in $\tau$ with period $2\pi$ and accurate solutions can be obtained by standard numerical methods for solving ordinary differential equations.

The evolution of both warped and unwarped discs is governed by transport of angular momentum, and therefore by an internal torque $\bm{\mathcal{G}}$. In a warped disc, there are three relevant components of $\bm{\mathcal{G}}$, with dimensionless coefficients (Ogilvie, 1999; Ogilvie & Latter, 2013a)

[TABLE]

where $T_{ij}$ are the components of the viscous stress tensor and the subscript $h,\tau$ indicates that a horizontal and time average over one period is taken. The vertical component of the torque, $Q_{1}$, involves the radial transport of vertical angular momentum and is similar to the usual torque in an unwarped accretion disc. If negative, it causes the mass distribution to evolve diffusively. The horizontal components of the torque, $Q_{2}$ and $Q_{3}$, involve the radial transport of horizontal angular momentum. If $Q_{2}>0$, it causes a diffusion of the warp, while $Q_{3}$ causes a dispersive wavelike propagation of the warp.

While the torque components can be calculated for the laminar flow induced by the warp (Ogilvie & Latter, 2013a), it was found that these flows may in fact be hydrodynamically unstable (Ogilvie & Latter, 2013b). Therefore, in order to measure the torque components in the nonlinear state of the instability one has to resort to numerical simulations. This is the purpose of this paper.

3 Numerical method

In order to use a finite volume volume numerical scheme we write the system (1)-(4) in conservative form:

[TABLE]

Note that, in the absence of a warp ($|\psi|=0$), this system reduces to the standard equations for a shearing box in shearing coordinates, and that it can be written as

[TABLE]

where ${\bf W}$ is the state vector, ${\bf F}_{xyz}$ are the flux vectors in three directions and ${\bf S}$ is the source vector. Then, we use the technique of operator splitting to treat each spatial dimension separately.

3.1 Dimensional splitting

3.1.1 $x$-direction

For the $x$-direction, we keep only the flux term with ${\bf F}_{x}$ and the associated source terms (see Eulderink & Mellema, 1995):

[TABLE]

Source terms are integrated using a Crank-Nicholson scheme allowing for exact integration of epicyclic motion (Stone & Gardiner, 2010). These equations are similar to the ordinary shearing box equations but for the extra source term $\propto|\psi|$, and we refer to Paardekooper (2012) for details of the implementation.

3.1.2 $y$-direction

In this paper, we only consider “axisymmetric" shearing boxes, i.e. no quantity depends on $y$. This situation corresponds to studying warped discs in which the azimuthal variation is on a global scale comparable to the circumference of the disc, rather than on the local length scale $H=c_{\mathrm{s}}/\Omega$. For completeness, the equations for the $y$-direction read:

[TABLE]

Note that for $q\tau\gg 1$ the system will be completely dominated by the shearing of the coordinate frame. Therefore, it is necessary to occasionally remap the solution to $\tau=0$, by shifting the solution in the $y$-direction by an amount $qx\tau_{\mathrm{remap}}$. For a carefully chosen $\tau_{\mathrm{remap}}$, this can be done without interpolation. Basically, if the grid is symmetric with respect to $x=0$, and the row with smallest $|x|$ has to be shifted by 1 grid cell (which can be done without interpolation), then the next row will have to be shifted by two cells (again without interpolation), etc.

3.1.3 $z$-direction

For the $z$-direction, we keep only the flux term with ${\bf F}_{z}$ and the corresponding source term:

[TABLE]

Fluxes now contain terms $\propto|\psi|\cos\tau$, which find their way into the eigenvalues of the Jacobian matrix $\mathsf{A}=d{\bf F}_{z}/d{\bf W}$ through functions $\xi\equiv|\psi|\cos\tau$ and $\eta\equiv\sqrt{1+\xi^{2}}$:

[TABLE]

and the corresponding eigenvectors:

[TABLE]

A flux difference at the interface between cells $i$ and $i-1$ is projected onto the eigenvectors of a suitably averaged Jacobian $\hat{\mathsf{A}}$:

[TABLE]

with projection coefficients:

[TABLE]

where ${\bf F}_{i}-{\bf F}_{i-1}=(\Delta_{\rho},\Delta_{x},\Delta_{y},\Delta_{z})^{T}$ and $\hat{{\bf v}}$ denotes the Roe-averaged velocity. The first order interface flux is then given by

[TABLE]

Second order corrections are obtained using a flux limiter (Paardekooper, 2012).

The dependence on $\tau$ is periodic, so there is no need to remap if simulating an axisymmetric disc (i.e. no $y$-dependence). Otherwise it is advantageous to pick $\tau_{\mathrm{remap}}$ a multiple of $2\pi$. In practice, this can be achieved by slightly modifying the extent of the domain in $x$ and $y$.

3.2 Viscosity

Viscous source terms are integrated separately from the hyperbolic part of the governing equations using the technique of operator splitting. In this step, we solve for the change in velocities due to viscous forces only, which in warped shearing coordinates take the form:

[TABLE]

where $T_{ij}$ denote the components of the viscous stress tensor, which, neglecting bulk viscosity, read:

[TABLE]

where the velocity divergence

[TABLE]

The constant shear viscosity $\nu$ is parameterised in the usual way as

[TABLE]

The viscous terms are integrated using a simple explicit Euler step, which was found to be accurate enough in the regime of interest.

4 Initial conditions and tests

In units of $H=c_{\mathrm{s}}/\Omega$, we take the computational domain to be within $-6<z/H<6$. The horizontal extent of the computational domain is such that it can fit an integer number of inertial waves of interest (see Section 4.2 below), but typically $-4<x/H<4$. This domain is covered by a uniform mesh of at least $(N_{x},N_{z})=(128,192)$, so that a scale height $H$ is resolved by at least $16$ cells. In the figures below, we will refer to this base resolution as “$N=128$", and “$N=256$" refers to $(N_{x},N_{z})=(256,384)$, etc. The boundaries are periodic in $x$, while for the vertical boundaries we use non-reflecting conditions (see Paardekooper & Mellema, 2006).

4.1 Laminar flows

The laminar flows discussed in section 2 provide a first basic test for the implementation of the Roe solver for the warped shearing box. For small $|\psi|$ in a Keplerian $q=1.5$ disc, one expects $u,v\propto|\psi|/\alpha$ and to be sinusoidal, with $w$ and $g-1$ a factor $|\psi|$ smaller (Ogilvie & Latter, 2013a). Such a solution is presented in the top panel of figure 1. Note that even for such a small warp amplitude $|\psi|=0.01$ the horizontal flow speeds exceed the sound speed sufficiently high above the mid plane of the disc (since for example $v_{x}=u(\tau)\Omega z$), which is due to the fact that the Keplerian case suffers from a resonance due to the coincidence of the orbital and epicyclic frequency. This resonance is moderated by viscosity, but since typically $\alpha\ll 1$ the resulting horizontal velocities can be very high. For larger warp amplitudes, the laminar flows become less sinusoidal and more extreme, as illustrated in the bottom panel of figure 1.

The laminar solutions obtained with the hydrodynamic solver compare well with those obtained using a standard ordinary differential equation solver. More challenging is to follow the laminar flow accurately for many (hundreds to thousands) periods, which is necessary when the growth rate of any instability feeding off the laminar flow is very small. Any evolution in the laminar flow can be conveniently observed in the torque component $Q_{2}$, which measures roughly speaking the diffusion of the warp (Ogilvie & Latter, 2013a). For the purely laminar flow one expects (Ogilvie & Latter, 2013a):

[TABLE]

In a one dimensional calculation, where all quantities depend only on $z$ and $t$, the laminar flows should be stable and therefore $Q_{2}$ should be constant. In reality, numerical diffusion can act to reduce the amplitude of the laminar motions. In figure 2 we show the evolution of $Q_{2}$ for four different resolutions and $|\psi|=0.01$, $\alpha=0.001$ and $q=1.5$. At low resolution, the amplitude of the laminar motions is reduced by more than $5\%$ over $2500$ periods. This reduction is mainly due to difficulties at high altitude $(z/H>3)$, where the densities get very low and any slight error in maintaining hydrostatic equilibrium may result in an unphysical state. The use of flux limiters can prevent this from happening at the expense of more numerical diffusion. This has the effect that the very upper regions at low resolution no longer participate in the laminar motion, resulting in a reduction of $Q_{2}$. A similar reduction in $Q_{2}$ would be obtained by choosing a smaller vertical domain. In fact, the reason for choosing the upper boundary at $z/H=6$ is based on getting the correct values of the torque components.

The numerical diffusion of the laminar flow was found to scale with the amplitudes of the flow. For example, for $|\psi|=\alpha=0.01$, the amplitudes are reduced by an order of magnitude relative to the case considered above, leading to the torque $Q_{2}$ component to be preserved to within $0.2\%$ even at our lowest resolution $N=128$.

4.2 Inertial waves

In the absence of a warp, the linearised governing equations allow for solutions with either an inertial or an acoustic character. In particular, it is possible to decompose perturbations to the basic state into Fourier-Hermite modes:

[TABLE]

where $\mathrm{He}_{n}$ is the $n$th Hermite polynomial and $L_{x}$ is the horizontal size of the domain. For a single Fourier-Hermite mode $(m,n)$ the solution $\propto\exp(i\omega\tau)$ where $\omega$ is given, in the inviscid case, by

[TABLE]

where $k_{x}=2\pi mc_{\mathrm{s}}/(\Omega L_{x})$ is a dimensionless wave number. The low frequency branch of (64) is the branch of inertial waves, and in the absence of a warp the various modes are not coupled. The same is true in the presence of viscosity as long as the bulk viscosity $\alpha_{b}=2\alpha/3$. It has been found (Ogilvie & Latter, 2013b) that two inertial waves can couple to the laminar flow induced by the warp, resulting in exponential growth. It is therefore important that inertial waves are captured accurately by the numerical scheme.

In order to test the ability of the method to sustain inertial waves, we set up an unwarped Keplerian box ($|\psi|=0$, $q=1.5$), and set up a linear inertial wave (amplitude $10^{-5}c_{\mathrm{s}}$) with $k_{x}=1.6$ and $n=1$. The horizontal domain is chosen such that two of these waves fit in the domain: $L_{x}=4\pi c_{\mathrm{s}}/(\Omega k_{x})\approx 7.854\,c_{\mathrm{s}}/\Omega$. The results for $\alpha=0$ are shown for various resolutions in figure 3. The temporal frequency of the inertial wave is $\approx 0.48$, so that we follow the wave for $\approx 40$ periods. During this time, at lowest resolution the wave amplitude is reduced by $50\%$. The major source of error appears to be the dimensionally split nature of the code, in combination with a flow pattern that does not follow the numerical grid (see e.g. Paardekooper, 2017). High resolution is necessary in order to follow the growth of inertial waves, in particular for small growth rates. When viscosity is included, the waves are damped on a time scale $c_{\mathrm{s}}^{2}\Omega\tau_{d}\sim k_{x}^{2}/\alpha$. In this case, numerical convergence is easier to obtain, and for $\alpha=0.01$ results with $N=128$ are indistinguishable from those obtained with $N=256$.

5 Results

5.1 Non-Keplerian disc ($q=1.6$)

While the main focus of this paper is on Keplerian discs with $q=1.5$, we start off by briefly showing results for a disc with $q=1.6$. The non-Keplerian case is slightly easier to set up, as it requires no viscosity and the amplitudes of the laminar motions are much lower for a given warp amplitude compared to the Keplerian case. Initial conditions are given by the laminar warp motion, on top of which we introduce two $k_{x}=1.6$ inertial waves with $n=1$ and $n=2$ of velocity amplitude $10^{-5}c_{\mathrm{s}}$. The evolution of the $n=1$ Fourier-Hermite component of the vertical velocity is shown in figure 4 for $|\psi|=0.01$ (top panel) and $|\psi|=0.04$ (bottom panel) for three different resolutions, together with the expected growth rate obtained from linear theory. For $|\psi|=0.01$, the growth rate is relatively small which means that quite high resolution is needed to match the expected growth rate. This is because the instability is battling against numerical diffusion of the inertial waves. For $N=256$ and $N=512$ the linear growth rate is recovered nicely. The case with $|\psi|=0.04$ has a much faster linear growth rate, which reduces the resolution requirements to an extent that the linear growth rate can now be recovered for $N=128$.

The torque components evolve erratically with time, as illustrated in figure 5. As expected from the laminar flows, initially only the $Q_{3}$ component is non-zero. In the end, all three warp amplitudes considered settle on an average value of $Q_{3}$ that is slightly reduced compared to the initial value, but with a large spread. Interestingly, also the $Q_{2}$ component, governing roughly speaking the diffusion of the warp, reaches similar values on average for all three warp amplitudes. Only for $Q_{1}$ do we observe a clear trend with $|\psi|$, with $Q_{1}$ increasing with increasing $|\psi|$. While there are large short time scale oscillations visible, on average $Q_{1}$ settles to a positive value. A positive value of $Q_{1}$ indicates a reversal of the angular momentum flux from a usual accretion disc, which means one could expect the mass distribution to behave in an antidiffusive manner, leading to a breakup of the disc into disjoint rings. However, the effect is very weak in this case.

The absence of viscosity in principle allows for structure on very small scales. In figure 6 we show snapshots of $\sqrt{\rho}|{\bf\delta v}|$, where ${\bf\delta v}$ is the velocity perturbation on top of the laminar motion, which is a measure of wave energy. The wave energy increases with warp amplitude, and there is a tendency for more small-scale structure at larger $|\psi|$. For $|\psi|=0.01$, the final state is dominated by the $k_{x}=1.6$ mode, independent of whether we initialize a single $k_{x}=1.6$ mode or start from white noise. While smaller scale inertial waves (larger wave numbers) have higher growth rates, these tend to saturate at low amplitudes, leaving enough of the laminar motions for the smaller wave number modes to feed off and grow. As a result, the nonlinear state is well captured (as indicated by the torque components) even at modest resolution of $N=128$, despite the fact that this resolution was unable to accurately reproduce the linear growth rates (see figure 4).

As the amplitudes of the laminar flow remain relatively modest for $q=1.6$, larger warp amplitudes can be considered without difficulty. In fact, the main limiting factor, in particular in runs without any explicit viscosity, is the small time steps necessary because of large velocities in the upper layers of the box. Here we note in particular that while no laminar solutions can be found for $|\psi|>0.261$ due to a nonlinear resonance (Ogilvie & Latter, 2013a), with regards to the nonlinear state this value of $|\psi|$ is not special. There exists a nonlinear state for $|\psi|=0.27$, which can be reached either by starting without any laminar flow or by starting from the nonlinear state at $|\psi|=0.26$ and increase the warp amplitude. The final state at $|\psi|=0.27$ is very similar to that at $|\psi|=0.26$ in terms of torque components and velocity structure.

5.2 Keplerian disc ($q=1.5$)

We now turn our attention to the interesting Keplerian case that has $q=1.5$. For $q=1.5$, viscosity is needed in order for laminar flow solutions to exist (Ogilvie & Latter, 2013a). We consider two viscosity coefficients: $\alpha=0.001$ and $\alpha=0.01$, always leaving the bulk viscosity zero.

5.2.1 Growth rates

In figure 7 we show the evolution of the $k_{x}=1.6$, $n=1$ Fourier-Hermite component of the vertical velocity for $|\psi|=\alpha=0.001$ (top panel) and $|\psi|=\alpha=0.01$ (bottom panel), for three different resolutions together with the growth rate obtained from a linear calculation. The linear growth rate is much larger in the case of $\alpha=0.001$, which makes it easier to capture accurately at low resolution, so that even $N=128$ does a reasonable job of reproducing the linear growth rate. The slower growth rate when $\alpha=0.01$ is only captured accurately for $N=256$ and above, even though the saturated state is very similar in all cases.

For $\alpha>0$, growth rates are reduced for large wave numbers, and for a given warp amplitude $|\psi|$ and viscosity $\alpha$, we can calculate the growth rate maximized over $k_{x}$ either by a linear calculation (Ogilvie & Latter, 2013b) or through direct numerical simulation by initializing the velocities with white noise and measuring the growth of the perturbed energy. Note that in the latter case, we obtain the maximum over all $k_{x}$ that fit into the box. In our standard setup, the minimum $k_{x}$ we can consider is $k_{x}=0.8$ (one wave fits in the box), and the allowed wave numbers are therefore $k_{x}=0.8,1.6,2.4,3.2$, etc. While in the picture of the pure parametric instability this may pose a problem, as a particular value of $k_{x}$ is required for the waves to couple (Gammie et al., 2000; Ogilvie & Latter, 2013b) which may not fit into the box, in practice, the instability bands for the values of $\alpha$ considered here are wide enough so that our sampling gives a good estimate of the maximum growth rate.

An issue with solving the linearized equations for the Fourier-Hermite components that was flagged up in Ogilvie & Latter (2013b) is that for strong laminar flows the spectrum fails to converge and no reliable growth rates can be obtained. The difficulty may lie in the Hermite basis being unable to handle warped discs with strong laminar flows, but it may also be that the oscillatory shear of the laminar flow may prevent the formation of axisymmetric eigenmodes altogether. An independent implementation of the numerical method based on Floquet theory discussed in Ogilvie & Latter (2013b) shows the same behaviour: for $\alpha=0.01$ we can find no reliable growth rates for $|\psi|>0.03$, while for $\alpha=0.001$ we can find no reliable growth rates for $|\psi|>0.003$.

The linear results, together with growth rates measured from the hydrodynamic simulations, are displayed in figure 8 for $\alpha=0.01$. For $|\psi|\leq 0.03$, the results obtained with the two methods agree very well. Interestingly, we can numerically find positive growth rates where the linear solver fails to converge. This could indicate that at least part of the problem with the linear solver lies indeed with the use of the Hermite basis. However, the measured growth rates quickly decline towards larger values of $|\psi|$ and in the end for $|\psi|\geq 0.07$ no growing modes are found at all. While it is possible that the hydrodynamic simulations suppress very small growth rates due to numerical damping, it is worth noting that the growth rate at $|\psi|=0.04$ is already reasonably well captured with our lowest resolution $N=128$. In addition, for this relatively large value of $\alpha$ one expects large wave numbers to be strongly damped, so that one would not expect the necessary resolution to be extremely large. The main requirement is for the laminar flows to be accurately represented, for which one needs $N=256$ at least for the larger warp amplitudes. Higher resolution and larger horizontal box sizes did not yield any growth rates for $|\psi|\geq 0.07$, which supports the conclusion that the shear imposed by the laminar flows prevents the formation of axisymmetric eigenmodes.

Lowering the viscosity to $\alpha=0.001$ yields a similar picture, although because of the stronger laminar flows difficulties arise at smaller warp amplitudes. In figure 9 we show again growth rates obtained from both the linearized equations and hydrodynamic simulations. No reliable growth rates were found from the linearized equations for $|\psi|>0.003$, while the hydrodynamic simulations show growing modes up to $|\psi|=0.03$. Beyond $|\psi|=0.03$, no growing modes were found (but see section 5.2.3 where we find a nonlinear instability), again despite increasing the resolution up to $N=2048$ and changing the horizontal size of the box. Note that, unlike in the case with $\alpha=0.01$, growth rates remain high after the linear results break down, before decreasing for $|\psi|>0.01$.

5.2.2 Saturated state

Once the perturbations grow to significant amplitude, the torque components $Q_{1}$, $Q_{2}$ and $Q_{3}$ start to be affected. This is because the perturbations feed off the laminar flow, which can therefore be strongly reduced in the saturated state. The torque components in the saturated state are important indicators of subsequent evolution of the warp. In particular, when $Q_{2}>0$ it can be interpreted as a diffusion coefficient for the warp. From equation (59), it is clear that laminar flows at small viscosity would lead to very fast evolution of the warp because of the factor $\alpha$ in the denominator. However, as these laminar flows were found to be unstable, one should look at the saturated state.

In figure 10 we show the evolution of the three torque components for the high viscosity case $\alpha=0.01$ for three different warp amplitudes. At early times, the laminar flows dominate the torques, and excellent agreement is found with the results of Ogilvie & Latter (2013a). The onset of the nonlinear phase is marked most clearly by a sharp drop in $Q_{2}$, accompanied by a drop in $Q_{1}$, which is barely visible for $|\psi|=0.01$, and an increase in $Q_{3}$. Towards late times, the disc settles into a saturated state, usually with well-defined values for the torque components (the exception being $|\psi|=0.02$, which we will comment on below). The drop in $Q_{1}$ signals that inward angular momentum transport is enhanced, up to a factor of 2 in the case of $|\psi|=0.04$, while for smaller warps the effect is reduced. The drop in $Q_{2}$ means that the diffusion of the warp is reduced by roughly a factor 2 in all cases. The third torque component $Q_{3}$ increases in all cases. While both the reduction in $Q_{2}$ and the increase in $Q_{3}$ are modest, since they are of comparable magnitude in the saturated state one might expect the warp to behave in a more wave-like manner compared to the laminar state, which has $Q_{2}>Q_{3}$. We will see below that this effect is more pronounced in the case where the viscosity is smaller.

The saturated state is usually dominated by modes of a single value of $k_{x}$. This is illustrated in the top panels of figure 11, where we show $\sqrt{\rho}|{\bf\delta v}|$, where ${\bf\delta v}$ is the velocity perturbation on top of any possible (horizontally uniform) laminar flow. This quantity is a measure of wave energy. The smallest warp amplitude $|\psi|=0.01$ is seen to be dominated by $k_{x}=1.6$ (note that because of the absolute value four maxima, as observed in the top left panel of figure 11, mean two wavelengths over the length of the box, which translates to $k_{x}=1.6$), while the two larger warp amplitudes are dominated by modes with $k_{x}=2.4$. In all cases, density perturbations remain below $\sim 10\%$ in the mid plane. Note that the wave energy goes to zero rapidly for $|z|>4c_{s}/\Omega$, far enough away from the vertical boundaries at $|z|=6c_{s}/\Omega$ to be confident that the solution is unaffected by the boundaries.

The torque components for $|\psi|=0.02$ show a second transition around $\Omega t=1000$, while all others settle into a steady saturated state after the onset of nonlinearity with no apparent further changes. In the case of $|\psi|=0.02$, before this second transition, the saturated state is dominated by modes with $k_{x}=2.4$, as shown in figure 11. However, in this state modes with $k_{x}=1.6$ can still couple to the reduced laminar flow and grow, albeit more slowly than they would have in the presence of the full laminar flow. Once these new modes reach nonlinear amplitudes, the second transition happens where the modes of different $k_{x}$ compete for dominance. No such growth of additional modes was observed for $|\psi|=0.01$ or $|\psi|=0.04$.

For the low viscosity case $\alpha=0.001$, we find similar but more extreme results. In figure 12 we show the evolution of the torque components for three different values of $|\psi|$. Since for the laminar flows, $Q_{2}\propto 1/\alpha$, initially $Q_{2}$ is a factor of 10 larger than for $\alpha=0.01$. In the saturated state however, even a modest warp of $|\psi|=0.005$ can significantly reduce the amplitude of the laminar flow, much more so than in the more viscous case of $\alpha=0.01$. Both $|\psi|=0.005$ and $|\psi|=0.01$ show a reduction in $Q_{2}$ of roughly a factor $100$, which should significantly reduce the diffusion of these warps.

The laminar flows associated with the warps shown in figure 12 predict a reversal of angular momentum flux for $|\psi|=0.01$ so that $Q_{1}>0$ initially ($|\psi|=0.005$ being a borderline case). Such a situation would see the mass evolve in an anti-diffusive manner, probably leading to a breakup of the disc into distinct rings, as observed in for example Nixon et al. (2012). In the saturated state, however, all three warp amplitudes have $Q_{1}<0$, where also for $|\psi|=0.001$ the saturated value of $Q_{1}$ is more negative than the initial, laminar, value. The third torque component $Q_{3}$ increases in magnitude roughly by an order of magnitude. This means that while in the laminar state we have $|Q_{2}|\gg|Q_{3}|$, in the saturated state we have $|Q_{2}|<|Q_{3}|$, and one might expect a more wave-like behaviour of the warp.

Similar to the higher viscosity case, the saturated state is dominated by either $k_{x}=1.6$ or $k_{x}=2.4$ (see bottom panels of figure 11). The largest warp amplitude ($|\psi|=0.01$) shows some small-scale structure around $z\approx\pm 2c_{s}/\Omega$. The mid plane perturbations in density remain modest, up to $\sim 5\%$ for $|\psi|=0.01$. Similar to the higher viscosity case, the wave energy is essentially concentrated in $|z|<4c_{s}/\Omega$.

5.2.3 Hysteresis

We are now in a position to survey the parameter space, finding $Q_{i}=Q_{i}(|\psi|,\alpha,q)$, where $Q_{i}$ stands for any of the three torque components, in the saturated state. We will only consider the Keplerian case, $q=1.5$, and consider two values of $\alpha=0.01$ and $\alpha=0.001$ as before.

In figure 13 we show the three torque components as a function of warp amplitude for the case of $\alpha=0.01$. The laminar results show the expected behaviour of $Q_{1}$ increasing quadratically with $|\psi|$, while both $Q_{2}$ and $Q_{3}$ are nearly constant (Ogilvie & Latter, 2013a). For $|\psi|>0.05$ we have that $Q_{1}>0$ and therefore angular momentum flux reversal for the laminar solution. In the saturated state, $Q_{1}$ is always below the laminar value, but all the same turns positive for $|\psi|>0.05$. Moreover, $Q_{1}$ never goes below $-3\alpha/2$, which would be the value for an unwarped laminar disc. Similarly, the value of $Q_{2}$ in the saturated state is always below the laminar value, indicating that the warp will diffuse more slowly than would be the case for purely laminar flow. The minimum value of $Q_{2}=9.2$ is reached for $|\psi|=0.045$.

Several transitions can be identified in the saturated curves in figure 13, most notably at $|\psi|=0.008$, where there is a discontinuity in $dQ_{2}/d|\psi|$, at $|\psi|=0.015$ where there is a jump in $Q_{2}$ and finally at $|\psi|=0.05$ where $Q_{2}$ jumps almost to its laminar value. These transitions are associated with a change in the dominant wave number. For $|\psi|<0.015$ the saturated state is dominated by $k_{x}=1.6$, switching to $k_{x}=2.4$ for $0.015<|\psi|<0.05$. While it looks like at $|\psi|=0.05$ there are no growing modes and the laminar state is recovered, there are growing $k_{x}=1.6$ modes but they saturate at very low amplitude. Beyond $|\psi|=0.05$, the solution is dominated by modes with $k_{x}=0.8$. The growth rate of these modes is relatively slow and they may involve higher-order mode couplings (Ogilvie & Latter, 2013b).

The low level of saturation around $|\psi|=0.05$ does not correlate with a reduced growth rate. As can be seen in figure 8, the growth rate is above $0.01$ for this warp amplitude. The low level of saturation can be understood from the growth rates, however, if we look for specific values of $k_{x}$, as shown in figure 14. Modes with $k_{x}=1.6$ can grow in three distinct regions of $|\psi|$. The amplitude at which the modes saturate depends on a competition between growth from feeding off the laminar flows and viscous or nonlinear damping. A mode with $k_{x}=1.6$ growing in a $|\psi|=0.015$ warp can grow until it becomes nonlinear, at which point the laminar flows are significantly reduced, as indicated by the reduction in $Q_{2}$. The reduced laminar flows correspond to a lower value of $|\psi|$, but from figure 14 we see that to the left of $|\psi|=0.015$ the growth rate is actually larger, making it possible for the mode to grow to larger amplitude. If the laminar flows are reduced too much, however, the growth rate first starts to decline before growth shuts off completely. This puts a constraint on the maximum amplitude these modes can reach. In the case of $k_{x}=1.6$ and $|\psi|=0.015$, the laminar flows can be reduced until they correspond to $|\psi|=0.005$, which is a significant reduction leading to large saturated wave amplitude, as indicated by the low saturated value of $Q_{2}\approx 11$. In contrast, starting at $|\psi|=0.06$, the laminar flows can only be reduced to $|\psi|\approx 0.05$, which is a small relative reduction leading to a low saturated wave amplitude, as indicated by a saturated value of $Q_{2}\approx 23$ that is close to the laminar value. Note that because there is a large single region where modes with $k_{x}=2.4$ can grow, these modes always saturate at large amplitudes.

Such distinct regions of non-zero growth rates as seen in figure 14 are not expected from the simple three-wave coupling model (Gammie et al., 2000; Ogilvie & Latter, 2013b), which is valid for $|\psi|$ small enough so that the shear from the laminar flow is much less than the Keplerian shear. The linear calculations presented in Ogilvie & Latter (2013b) show that for larger warp amplitudes, the instability bands can shift to different values of $k_{x}$(see for example their figure 5, where $k_{x}=1.6$ is inside a band of instability for $|\psi|=0.01$, but not for $|\psi|=0.02$, which is consistent with our results). We find that at even larger warp amplitudes, instability can return to $k_{x}=1.6$, but in the end disappears completely. A shift in frequency of the most unstable modes towards larger warp amplitudes can also be obtained from the incompressible model presented in Gammie et al. (2000), but for this model the instability band widens at the same time so that $k_{x}=1.6$ always remains unstable. It is therefore likely that compressibility plays an important role in this phenomenon.

Lowering the viscosity to $\alpha=0.001$ reveals new interesting behaviour. While we could not find any growing modes for $|\psi|>0.03$ (see figure 9), it is possible to find a nonlinear state for warp amplitudes up to at least $|\psi|=0.1$. This state can be reached for example by starting from the saturated state at lower warp amplitude and slowly increasing $|\psi|$, or by starting without any laminar flow. The system therefore exhibits hysteresis. It takes time to set up the laminar flow either from scratch or from a lower value of $|\psi|$, and as long as modes can grow to sufficient strength before the full laminar flow is reached, the resulting state will be the saturated nonlinear state rather than the laminar flow. It should be noted that hysteresis appears to be limited to the low viscosity case: for $\alpha=0.01$ we always recover the laminar state for $|\psi|\geq 0.07$, no matter the initial conditions.

The results for $\alpha=0.001$ are displayed in figure 15. First thing to note is that since we are considering larger warp amplitudes compared to the more viscous runs with $\alpha=0.01$, the laminar torques start to differ from the quadratic form (for $Q_{1}$) and the constant value (for $Q_{2}$) towards larger $|\psi|$, while $Q_{3}$ remains roughly constant. The torque component $Q_{1}$ saturates at a value of $\sim 0.1$, while $Q_{2}$ starts to decrease for $|\psi|>0.01$. As is the case for the more viscous runs, $Q_{1}>0$ for larger enough $|\psi|$, in this case $|\psi|\geq 0.005$.

The torque components in the saturated state vary much more smoothly compared to the more viscous runs. This is probably due to more unstable modes being available for $\alpha=0.001$. For all warp amplitudes we find $Q_{1}$ in the saturated state to be much smaller than in the laminar state. Angular momentum flux reversal occurs for $|\psi|>0.04$, and $Q_{1}$ keeps increasing for larger warp amplitudes. Except for the lowest warp amplitudes $|\psi|<0.002$, $Q_{2}$ is reduced substantially in the saturated state compared to the laminar state. This is most prominent at $|\psi|\sim 0.01-0.02$, where $Q_{2}$ is reduced by two orders of magnitude. For larger warp amplitudes the effect is less dramatic, mainly due to the lower values of $Q_{2}$ in the laminar state. The third torque component shows a large variation with warp amplitude in the saturated state, unlike the laminar state for which $Q_{3}$ is roughly constant. Note that while $|Q_{2}|\gg|Q_{3}|$ in the laminar state for essentially all values of $|\psi|$, the saturated state has $|Q_{2}|\lesssim|Q_{3}|$, perhaps indicating a more wavelike behaviour of the warp.

For warp amplitudes $|\psi|>0.1$ it becomes more difficult numerically to sustain the nonlinear saturated state. For example, at resolution $N=512$ at $|\psi|=0.12$, a saturated state can be set up but it ultimately decays towards a laminar flow. The same simulation at a higher resolution $N=1024$ shows no signs of decaying. It may be the case that for $\alpha=0.001$ a saturated state exists for all values of $|\psi|$ but the regime $|\psi|>0.2$ will be very expensive to reach computationally.

5.2.4 Inviscid limit

Finally, we briefly consider the inviscid limit. In the limit of vanishing viscosity, laminar solutions cease to exist for a Keplerian rotation profile (Ogilvie & Latter, 2013a). One may first of all wonder whether a saturated state like those found above still exists. We start from purely Keplerian flow with small amplitude white seed noise, and impose a warp $|\psi|$ as before. Results after $\Omega t=200\pi$ are shown in figure 16. For both the small ($|\psi|=0.01$, top row) and large ($|\psi|=0.1$, bottom row) warp amplitude a saturated state very much like the viscous case is set up. Towards higher resolution, smaller scale flow structures can be seen. These did not feature in the viscous runs, compare for example figure 11. However, the measured torque components for $|\psi|=0.01$ ($Q_{1}=10^{-5}$, $Q_{2}=3.0$, $Q_{3}=6.0$) are independent of resolution and agree very well with the viscous case that has $\alpha=0.001$. The small scale features therefore do not affect the internal torque. The agreement with $\alpha=0.001$ suggests that this viscosity is representative of a “low viscosity" regime where the torque components do not sensitively depend on the viscosity.

For the large warp amplitude ($|\psi|=0.1$, bottom row of figure 16), again the measured torque components ($Q_{1}=0.016$, $Q_{2}=8.0$, $Q_{3}=-2.0$) agree well with those obtained for $\alpha=0.001$. However, numerical convergence is more difficult to obtain. While $N=512$ and $N=1024$ give identical results, $N=256$ for example gives $Q_{2}=13.0$. At even lower resolution $N=128$ (not shown in figure 16), the flow settles in an almost laminar state with very low amplitude waves superimposed. This again highlights the strong resolution constraints towards larger warp amplitudes. Even when a completely laminar state does not exist, numerical diffusion will tend to take the flow towards a state where any non-laminar component has very low amplitude.

The large warp amplitude goes hand in hand with strong density perturbations, as illustrated in figure 17. Note that the maximum of the unperturbed density is defined to be unity. While for $|\psi|=0.01$ perturbations due to the nonlinear inertial waves can be seen, for $|\psi|=0.1$ density perturbations are much stronger and feature much smaller scales including some weak shocks. This type of flow persists to at least $|\psi|=0.2$.

6 Discussion and conclusions

We have investigated the hydrodynamic stability of the laminar flows driven by the presence of a warp in the disc. We have confirmed the existence and growth rates of a linear parametric instability (Gammie et al., 2000; Ogilvie & Latter, 2013b), both in inviscid non-Keplerian discs and in Keplerian discs for two values of the viscosity parameter, $\alpha=0.01$ and $\alpha=0.001$. At high viscosity, the saturated state consist mostly of wave activity (at nonlinear amplitudes), while towards lower viscosity, a more turbulent state is found. In the Keplerian case, the internal torques $Q_{1}$ and $Q_{2}$ were reduced, most severely at low viscosity. Towards larger warp amplitudes, we find that $Q_{1}>0$, indicating that the mass distribution would evolve in an anti-diffusive manner, possibly breaking the disc up into distinct rings. The reduction in $Q_{2}$ shows that the warp diffuses on a much longer time scale than one would conclude based on the laminar flows. Since in the saturated state, we no longer have $|Q_{2}|\gg|Q_{3}|$ one might expect the warp to evolve in a more wavelike manner.

This hydrodynamic instability to the best of our knowledge has not been seen in global simulations of warped accretion discs. Armed with the results of section 5, we can shed some more light on why this may be the case, focusing on Keplerian discs. First of all, there is a minimum resolution for the instability to develop. While even at our lowest resolution ($\approx 16$ cells per $H$) we found growing modes, at high viscosity the growth rate was underestimated by the low-resolution simulations (see figure 7). Note that this resolution is still three times higher than used in Fragner & Nelson (2010), who used 5 cells per scale height in the meridional direction, while Lodato & Price (2010) achieved $\sim 3.5$ SPH smoothing lengths per scale height at their highest resolution. For larger warps at low viscosity, the required resolution for capturing the saturated state is more like 128 cells per $H$ at $|\psi|=0.1$, a resolution that is very difficult to achieve in a global calculation.

The amplitude at which the instability saturates depends on the level of viscosity. At $\alpha=0.01$, the amplitude remains relatively small and the effect on the torque components is modest ($Q_{2}$ is reduced by a factor of $\sim 2$ at most), while at $\alpha=0.001$ we find larger amplitudes and a strong effect on the torque components ($Q_{2}$ is reduced by a factor $\sim 100$). If this trend continues towards higher viscosity, the instability may even be difficult to identify for larger values of $\alpha$, which are commonly used in SPH simulations, in combination with a bulk viscosity (e.g. Lodato & Price, 2010), which we have ignored.

In addition to a resolution constraint, there is the interesting result that no growing modes were found towards larger values of $|\psi|$, probably due to the strong shear in the laminar flow. At large viscosity ($\alpha=0.01$), only the laminar state is available to the flow for $|\psi|>0.07$. It may be that many global calculations are in such a regime. At low viscosity ($\alpha=0.001$), growing modes cease to exist for $|\psi|>0.03$. Beyond this warp amplitude, both a laminar and a saturated state are available, and in order to recover the saturated state the initial conditions have to be far enough away from the laminar state. We have found this could be achieved by either starting from a saturated state at lower warp amplitude or by starting from purely Keplerian flow. Note however that the regime $\alpha\lesssim 0.01$ is difficult to achieve with SPH (Lodato & Price, 2010).

We have considered only the simplest form of the warped shearing box: isothermal and no $y$-dependence. These assumptions should be relaxed in future studies. In particular, a stably stratified disc would give rise to inertia-gravity waves rather than pure inertial waves, while a 3D box would allow for non-axisymmetric shearing waves that may be transiently amplified. Both may have an impact on the growth and saturated state of the hydrodynamic instability.

Acknowledgements

SJP is supported by a Royal Society University Research Fellowship. GIO is supported by STFC grant ST/P000673/1.