Disc tearing: numerical investigation of warped disc instability
aa r X i v : . [ a s t r o - ph . H E ] J a n Draft version January 18, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Disc tearing: numerical investigation of warped disc instability
A. Raj, C. J. Nixon, and S. Do˘gan School of Physics and Astronomy, University of Leicester, Leicester, LE1 7RH, UK Department of Astronomy and Space Sciences, University of Ege, Bornova, 35100, ˙I zmir, Turkey ABSTRACTWe present numerical simulations of misaligned discs around a spinning black hole covering a rangeof parameters. Previous simulations have shown that discs that are strongly warped by a forcedprecession—in this case the Lense-Thirring effect from the spinning black hole—can break apart intodiscrete discs or rings that can behave quasi-independently for short timescales. With the simulationswe present here, we confirm that thin and highly inclined discs are more susceptible to disc tearing thanthicker or low inclination discs, and we show that lower values of the disc viscosity parameter lead toinstability at lower warp amplitudes. This is consistent with detailed stability analysis of the warpeddisc equations. We find that the growth rates of the instability seen in the numerical simulationsare similar across a broad range of parameters, and are of the same order as the predicted growthrates. However, we did not find the expected trend of growth rates with viscosity parameter. Thismay indicate that the growth rates are affected by numerical resolution, or that the wavelength of thefastest growing mode is a function of local disc parameters. Finally, we also find that disc tearing canoccur for discs with a viscosity parameter that is higher than predicted by a local stability analysis ofthe warped disc equations. In this case, the instability manifests differently producing large changesin the disc tilt locally in the disc, rather than the large changes in disc twist that typically occur inlower viscosity discs.
Keywords:
Accretion, accretion discs — Hydrodynamics — Instabilities — Black hole physics INTRODUCTIONAccretion discs are generally considered to orbit theircentral accretor in a smooth sequence of planar and cir-cular orbits. However, there is a substantial body of ob-servational evidence that shows that real discs are morecomplicated. For example, warped discs are used tomodel the X-ray binary Her X-1 (e.g. Scott et al. 2000),water masers reveal discs warps around supermassiveblack holes (e.g. Miyoshi et al. 1995; Greenhill et al.1995), and the recent spatially resolved observationsof protoplanetary discs have directly revealed complexdisc structures including disc warps (e.g. Casassus et al.2015). Very recently it has been suggested that a mis-alignment between the gas and black hole angular mo-mentum may be required to resolve discrepancies in themass measurements for the supermassive black hole inM87 that was imaged by the Event Horizon Telescope(Jeter & Broderick 2020). In different systems warps [email protected] may form in different ways. For example, warps canform during the formation of an accretion disc (e.g.in chaotic star forming regions; Bate et al. 2010; Bate2018), or the disc may become warped if the orbits areforced to precess in a radially differential manner by thespin of the central object (e.g. a spinning black hole;Bardeen & Petterson 1975) or the gravitational fieldof a binary system (e.g. Papaloizou & Terquem 1995;Larwood & Papaloizou 1997).To a first approximation warps can propagate througha disc in one of two ways dependent on the relativemagnitudes of the Shakura & Sunyaev (1973) dimen-sionless viscosity parameter, α , and the disc angularsemi-thickness, H/R (Papaloizou & Pringle 1983). Forlow viscosity discs, where α < H/R , warps propagatevia bending waves. For high viscosity discs, where α > H/R , warps propagate following a diffusion equa-tion. In this work we are principally interested in thediffusive case, although we do present some simulationswith α . H/R . For a review of the basics of warpeddisc dynamics see Nixon & King (2016).
A. Raj et. al.
Recently it has been possible to explore detailed dy-namics of warped discs with numerical hydrodynamicalsimulations. There is now a broad range of simulationscovering a variety of astrophysical setups, including; (1)isolated warped discs (e.g. Nelson & Papaloizou 1999;Lodato & Price 2010), (2) warped discs in binary sys-tems (e.g. Larwood et al. 1996; Larwood & Papaloizou1997; Fragner & Nelson 2010), and (3) warped discsaround spinning black holes (e.g. Nelson & Papaloizou2000; Fragile et al. 2007). A recurring feature of thewarped discs in simulations is that they can ‘break’ or‘tear’ (e.g. Nixon et al. 2012a), and this typically occurswhen the simulation parameters allow for large enoughwarp amplitudes (either present at the start of the simu-lation, or formed over time) at sufficiently low viscosity.The warp amplitude, ψ , is given by ψ ( R, t ) = R (cid:12)(cid:12)(cid:12)(cid:12) ∂ l ∂R (cid:12)(cid:12)(cid:12)(cid:12) , (1)where R is the orbital radius, and l ( R, t ) is the unit ‘tilt’vector pointing in the direction of the disc angular mo-mentum vector. The unit tilt vector is usually writtenas l = (cos γ sin β, sin γ sin β, cos β ), where β ( R, t ) is thedisc ‘tilt’ and γ ( R, t ) is the disc ‘twist’. From this wecan see that discs can be warped in distinct ways. Forexample, a disc may achieve a warp ( ψ >
0) if β varieswith radius or if γ varies with radius (and β = 0). Weillustrate these two cases in Fig 1. As we shall see below,a large warp amplitude can be caused by a sharp varia-tion in β or a sharp variation in γ . In general a warpeddisc has both β and γ varying with radius. Locally, onscales ∼ H , the disc doesn’t know whether misalignedneighbouring rings are tilted or twisted with respect toeach other, and thus the stability of the disc (as deter-mined by a local stability analysis) is determined by ψ and not β or γ . However, as we shall see below the globalbehaviour of the disc can result in the local instabilitymanifesting as either a sharp change in the disc twist(Nixon et al. 2012a) or a sharp change in the disc tilt(see Section 3.5).The nomenclature for unstable warped discs, i.e.‘break’ and ‘tear’, is ill-defined in the literature. By adisc ‘break’ we mean (see also, for example, Doˇgan et al.2018) that the disc has a sharp transition in either thedisc tilt angle ( β ) or the disc twist angle ( γ ), and thatthis is usually accompanied by a sharp variation in the Note that here, for clarity, we refer to the warp amplitude by thesymbol ψ . This is different to some previous works which refer tothe warp amplitude as | ψ | due to the use of complex variables inOgilvie (1999) in which the warp amplitude was the magnitude, | ψ | , of a complex variable ψ (see equations 44-45 of Ogilvie 1999).Here we take ψ to be the warp amplitude. disc ‘surface’ density at the same radius. Breaks mayarise in a warped disc that achieves a sufficient warp am-plitude and can occur with or without a forced preces-sion of the disc orbits (see, for example, Nixon & King2012). The necessary conditions for disc ‘breaking’ arediscussed by Doˇgan et al. (2018) and Do˘gan & Nixon(2020). By disc ‘tearing’ we mean that an otherwisestable disc has been forced to precess (differentially) suf-ficiently rapidly that it achieves a large enough warpamplitude that renders the disc unstable to disc ‘break-ing’, and that the resulting parts of the broken disc areable to subsequently precess quasi-independently beforeinteracting in a highly dynamic and variable fashion.With this nomenclature, we have that disc ‘breaking’ isthe underlying instability (cf. Doˇgan et al. 2018), anddisc ‘tearing’ refers to the action of breaking a disc witha forced precession and the subsequent dynamical be-haviour of the unstable disc (cf. Nixon et al. 2012a).Larwood et al. (1996) and Larwood & Papaloizou(1997) present numerical simulations using smoothedparticle hydrodynamics (SPH) of misaligned discs in bi-nary systems, including both the circumbinary disc caseand the case where the companion is external to thedisc. In each case they show a simulation, typically thethinnest simulation presented, which exhibits disc tear-ing; the gravitational torque on the disc causes the discorbits to precess and a break is formed between the in-ner and outer regions of the disc. They note that inthese models the sound speed is too low to allow effi-cient communication between the different regions of thedisc. Similar results are presented by Fragner & Nelson(2010) who modelled the hydrodynamics using a gridcode. These investigations exhibit disc tearing (break-ing of the disc caused by a forced precession), but didnot follow the evolution for long enough to explore thesubsequent nonlinear, chaotic behaviour. The dynami-cal behaviour of disc tearing was explored by a series ofpapers in different contexts, namely; discs around spin-ning black holes (Nixon et al. 2012a), circumbinary discs(Nixon et al. 2013) and discs with an external binarycompanion (Do˘gan et al. 2015). These simulations re-veal the quasi-independent precession of rings of matter Note that for a planar disc orbiting in, say, the x - y plane, the‘surface’ density is defined as the integral of the volume densitywith vertical height through the disc and is typically averaged inazimuth to give a single value for each radius. In a warped disc,the ‘surface’ density retains this meaning, but now the ‘vertical’extent of the disc is taken normal to the local orbital plane and isa local quantity measuring the amount of mass orbiting per unitarea at each radius. This is distinct from the column integralwhich would measure the amount of mass along any given sightline as seen by a distant observer. isc tearing tilted twisted tilted R ψ twisted R ψ Figure 1.
Example disc structures depicting the way in which a disc can be warped. On the left hand side we show a threedimensional rendering of a warp imposed on a disc where only the tilt angle, β , varies with radius. On the right hand side weshow a similar structure, but this time only the disc twist angle, γ , is varied with radius. In each case the innner disc radius isat R = 4, and the outer disc radius at R = 120. For the tilted case, we take γ = 0 ◦ at all radii, while β = 0 ◦ for R <
40 and β = 45 ◦ for R >
80. Between these two radii we vary β smoothly from 0 ◦ to 45 ◦ with a cosine bell. This results in the warpamplitude with radius, ψ ( R ), shown in the plot below the image; in this case we have ψ = R | ∂β/∂R | . For the twisted case, wetake β = 45 ◦ at all radii, with γ = 0 ◦ for R <
40 and γ = 45 ◦ for R >
80. Between these two radii we vary γ smoothly from 0 ◦ to 45 ◦ with a cosine bell. This results in the warp amplitude with radius, ψ ( R ), shown in the plot below this image, and in thiscase we have ψ = R | sin β ∂γ/∂R | . We have shown these discs in an orientation in which the outer plane is the same in eachcase. (some of which are radially narrow with ∆ R ∼ H , andothers more extended with ∆ R ≫ H ). These rings arefound to precess until the internal angle between twoneighbouring rings becomes large, at which point theycan interact violently, with shocks between rings pro-ducing gas on low angular momentum orbits which canfall to smaller radii. This material can fall directly onto the central accretor or circularise at a new smallerradius. These processes can dramatically change thedisc evolution and its observational character (see, e.g.,Nixon & Salvesen 2014, for discussion). Over the lastfew years disc tearing has been explored in the lowviscosity case (Nealon et al. 2015), has been found ingeneral relativistic magnetohydrodynamic simulations(Liska et al. 2020), and has been used to model observa-tions of the circumstellar disc in GW Ori (Kraus et al.2020). Recently, Doˇgan et al. (2018) connected the discbreaking and tearing behaviour observed in numericalsimulations with the instability of warped discs thatwas derived by Ogilvie (2000) through a local stabilityanalysis of the warped disc equations. Ogilvie (1999)showed that in some circumstances as the warp ampli-tude is increased the local “viscosity” that resists thewarping (e.g. ν ; Pringle 1992) decreases. Thus themore warped the disc gets, the less it is able to resistwarping further. This is found to occur for viscositiesthat are not much larger than α = 0 .
1, and at a crit-ical warp amplitude, ψ c , that depends on α ; typicallylower α yields lower critical warp amplitudes. WhileDoˇgan et al. (2018) focussed on discs with Keplerian ro-tation, Do˘gan & Nixon (2020) explored the instability inthe case of non-Keplerian rotation. A. Raj et. al.
Here we perform numerical hydrodynamical simula-tions to explore the instability of warped discs and thedisc tearing behaviour. To provide the forced preces-sion we simulate discs around spinning black holes, andwe vary the disc parameters including the viscosity pa-rameter α , the disc angular semi-thickness H/R and theinclination of the disc θ with respect to the rotation axisof the black hole. In Section 2 we discuss the numericalmethodology and expectations for the simulations. InSection 3, we report the results of the numerical simu-lations. In Section 4 we present our conclusions. NUMERICAL CONSIDERATIONSIn this section we discuss the methodology that is em-ployed in our numerical simulations (reported in Sec-tion 3 below), and how numerical effects may impactthe results of our simulations. We use smoothed par-ticle hydrodynamics (SPH) to model the accretion disc(see e.g. Price 2012), and in particular we use the SPHcode phantom (Price et al. 2018). To generate a time-dependent warp in the disc, with which we can analysethe evolution, we perform simulations of discs that aremisaligned to a central spinning black hole. This in-cludes both apsidal (Einstein) precession of eccentricdisc orbits, and nodal (Lense-Thirring) precession ofmisaligned disc orbits.Following Nelson & Papaloizou (2000) we employ apost-Newtonian approximation to model the gravityof the black hole. We use the Einstein potential forthe apsidal precession with a gravito-magnetic forceterm to induce the required Lense-Thirring precession(see Nixon et al. 2012a; Price et al. 2018, for details) .While the Einstein potential provides a good descrip-tion of the expected apsidal precession of the disc or-bits, it does not produce an innermost stable circularorbit (ISCO) to define the inner edge of the disc. In-stead we place an accretion radius at the location ofthe ISCO to truncate the disc there. The ISCO isdefined by the location where circular orbits becomeunstable, which is equivalent to the square of the di-mensionless epicyclic frequency, ˜ κ = κ / Ω , becom-ing negative. In terms of the orbital shear parameter, q = − d lnΩ / d ln R = (4 − ˜ κ ) /
2, this is where q = 2.For the post-Newtonian approximation we employ, q increases as the radius decreases, but q < a = 0 . R ISCO = 4
GM/c , the shear parameter at We refer the reader to Liptai & Price (2019) for a recent imple-mentation of the SPH equations of motion in a General Relativityframework. the ISCO is q ≈ .
75. We plot the shear parameterand the apsidal precession rate for the Keplerian, Ein-stein and Kerr cases in Fig. 2. This shows that theEinstein potential accurately models the apsidal preces-sion of disc orbits, but the local shear rate in the discis only accurate to within approximately 10 per cent.The numerical method implemented in the phantom
SPH code to model Lense-Thirring precession has beenshown to accurately reproduce the required precessionrate (e.g. Fig. 17 of Price et al. 2018). As shown byDo˘gan & Nixon (2020) the growth rates of the instabil-ity depend on the value of the shear rate (see also Fig. 3below).In numerical simulations it is important to distinguishbetween physical and numerical viscosity. By physicalviscosity we refer to the Shakura & Sunyaev (1973) vis-cosity that we employ here to model the angular momen-tum transport in discs that is generally thought to arisedue to disc turbulence created by the magnetorotationalinstability (Balbus & Hawley 1991). In different astro-physical systems the efficiency of angular momentumtransport can differ, and this is usually encapsulated inthe Shakura & Sunyaev (1973) α -parameter. In fullyionised discs α is typically found to be large (King et al.2007; Martin et al. 2019a), and may be lower in discsthat are only partially ionised (see Martin et al. 2019afor a discussion). Here we employ several values of α toexplore the possible range of disc evolution. This physi-cal viscosity is implemented as a direct (Navier-Stokes)viscosity term in the simulations (see Section 3.2.4 ofLodato & Price 2010, and Flebbe et al. 1994).Numerical solution of the equations of hydrodynam-ics also leads to numerical viscosity (as is the case forany numerical method; see, for example, the discussionin the Appendix of Sorathia et al. 2013). In the SPHmethod, numerical viscosity is added explicitly to theequations of motion as it is required to smooth discon-tinuities in the velocity field. We use the standard SPHartificial viscosity terms, with a linear term (with co-efficient α SPH ) and a quadratic term (with coefficient β SPH ). We take β SPH = 2 and allow α SPH to vary fol-lowing the switch proposed by Cullen & Dehnen (2010)with minimum value of 0.01 and maximum value of unity(see Price et al. 2018, for details). These numericalterms contribute a small level of viscosity, which whencomputed in terms of a Shakura & Sunyaev parameter istypically of the order of α AV ∼ .
01, but note that thisis resolution dependent (see equation 2 below). Ideally,this numerical viscosity should be small compared to thephysical viscosity. Unfortunately this is not always thecase, particularly when simulating low physical viscosityor thin discs. In these cases the numerical viscosity may isc tearing log R/R g q log R/R g l og Ω a p -(cid:0)(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6) KeplerianEinsteinEinstein + LTKerr
Figure 2.
Orbital shear rate (left panel) and apsidal precession rate (right panel) for the Keplerian potential (black solid line),Einstein potential (red dashed line), Einstein potential with Lense-Thirring gravito-magnetic force (green long-dashed line) andfrom the Kerr metric (blue dotted line). The Keplerian shear rate is constant, with q = 1 .
5, and the apsidal precession rate iszero in this case (not depicted). The curves have been evaluated for a black hole spin parameter of 0.5585. The equation forthe Kerr shear parameter is taken from Gammie (2004), while the apsidal precession rate is from Kato (1990). The Einsteinpotential provides an accurate description of the apsidal precession rate. However, the shear parameter is only approximatelycorrect with errors of the order of 10 per cent. play a role in determining the dynamics of the insta-bility, especially when the critical warp amplitude andgrowth rates of the instability depend sensitively on thevalue of the disc viscosity (see, for example, Figure 3 ofDoˇgan et al. 2018).It is also worth noting that in SPH simulations, the lo-cal resolution follows the local density of particles as thesmoothing length (roughly the resolution lengthscale) h ∝ ρ − / , where ρ is the density. This is a highlyuseful feature of the SPH method as it allows, for ex-ample, converging flows to be modelled with resolutionthat increases in regions of increasing density (e.g. grav-itational fragmentation of a fluid). However, this fea-ture also provides a challenge in modelling flows wherethe density drops significantly below the mean density.Such a drop in density occurs in warped discs in regionswhere the warp amplitude is high (Nixon & King 2012),and for an unstable disc which breaks into discrete ringsthe ‘surface density’ is significantly reduced between therings. Therefore it is possible that as the disc becomesunstable, the numerical viscosity could increase and be-come comparable to, or larger than, the physical vis-cosity in the simulation. However, as the instability isstabilised by larger viscosity (Doˇgan et al. 2018), we ex-pect this to have a stabilising effect on the simulations– and thus if instability is observed, we expect that the disc is physically unstable . Such instability in low res-olution simulations may not capture the correct growthrates or extent of the instability, but at high enough res-olution (such that the physical viscosity is significantlylarger than the numerical viscosity) it should be possibleto capture the relevant dynamics.In numerical simulations there will be a bulk viscos-ity present, which was not included in, for example,Doˇgan et al. (2018) which employed α b = 0. The bulkviscosity has an effect on the viscosity coefficients andthus on the stability and growth rates. Lodato & Price(2010) report that the magnitude of the bulk viscositypresent in an SPH simulation is α AVb ≈ α AV /
3, where α AV is the shear viscosity arising from the numericalviscosity . The magnitude of the shear viscosity arisingfrom the numerical viscosity in the continuum limit is We note that this may not be true of simulations performed inthe wavelike regime, where α < H/R . In this case, the analyticalwork is not as well-established as in the diffusive case, and thestability criteria is not known. Nealon, Price, & Nixon (2015)report numerical simulations of discs with α < H/R and find thatthe disc can break. However, it is not clear that the numericalviscosity was sufficiently small that these simulations were in thewavelike regime at the time instability was found. We shall returnto this question in the future. We note that an alternative method for modelling the physi-cal viscosity in SPH simulations is to scale the numerical vis-cosity term to the appropriate value for a Shakura & Sunyaev(1973) viscosity (e.g. Murray 1996; Lodato & Price 2010). In thismethod, the bulk viscosity will always be significant compared tothe shear viscosity.
A. Raj et. al. (Meru & Bate 2012) α AV = 31525 α SPH h h i H + 970 π β SPH (cid:18) h h i H (cid:19) , (2)where h h i /H is the shell-averaged smoothing length perdisc scale-height. For a typical value of h h i /H = 0 . h α SPH i ≈ α AV ≈ .
02, and thus α AVb ≈ . α = 0 .
03, 0 . . α b = 0 and α b = 0 .
03. We also show the re-sults for three different values of the disc shear, givenby q = 1 .
5, 1 . .
7. These plots show that thebasic picture is unchanged; the discs are generally un-stable at large warp amplitudes and small α , while gen-erally stable at small warp amplitudes and large α (seeOgilvie 2000; Doˇgan et al. 2018; Do˘gan & Nixon 2020,for details). However, the details are changed with theinclusion of non-zero bulk viscosity and non-Keplerianrotation. For example, the exact value of the criticalwarp amplitude can vary, and the growth rates at largewarp amplitude can vary by factors of several.From the above it is clear that numerical modellingof the disc breaking instability is likely to yield sub-tle effects. In particular ‘convergence’ of the simula-tions with increasing spatial resolution may not be easyto achieve. As we have noted above, the critical warpamplitudes and growth rates can depend sensitively onthe value of the viscosity and thus may depend on theresolution. With modern day resolution the numericalviscosity in SPH simulations is not negligible, and fordiscs that are sufficiently thin ( H/R ≪ .
1) to real-istically model black hole accretion, we find typicallythe numerical viscosity is several 10s of per cent of thephysical viscosity. For understanding the dynamics ofthese discs this is not a grave issue, as the net effect isthat the simulations are somewhat more viscous thananticipated, i.e. if the numbers given above are repre-sentative, then for a physical α = 0 .
03 we have a to-tal simulated viscosity of α ≈ .
05, while for α = 0 . α ≈ . . However, for interpreting the be-haviour of simulations—for example, measuring the in-stability growth rates and comparing them with theo-retical predictions—we must be aware that the numer-ical terms provide an additional, resolution-dependent It has been shown by Lodato & Price (2010) that the standardSPH numerical viscosity behaves in an ‘isotropic’ manner, that isto say it behaves like a Navier-Stokes viscosity, and thus one canadd the viscosities in this way. This may not be true for othernumerical methods. effect. Therefore we anticipate that, for simulations inwhich the disc is ‘resolved’ (that is h h i /H <
1) as theresolution of the simulations is increased we should re-cover the same basic dynamics (such as whether the discis stable or unstable and whether the disc aligns to theblack hole spin or not), but the precise details (such asthe exact location of disc breaks and the exact growthrate in the unstable region) may vary .With these caveats in mind we return to the focus ofthe paper, namely to compare numerical simulations ofwarped discs with the instability of warped discs foundby analytical calculations to occur at low viscosity andlarge warp amplitude (Ogilvie 2000; Doˇgan et al. 2018).We note here, that the analytical work on this insta-bility takes the form of a local linear stability analysiswhich makes use of the assumption that the perturba-tions grow (or decay) on timescales short compared withthe evolution of the unperturbed state. In a dynamicsimulation in which the background state is continu-ously evolving, it is likely that the instability is only re-alised in situations where the growth rate is sufficientlyhigh that the instability has time to act before the discevolves further. In this sense, the solutions presentedby Doˇgan et al. (2018) represent necessary, but not suf-ficient, criteria for instability, with the sufficiency sup-plied by the condition that the background state evolvemore slowly than the growth of unstable modes. Tak-ing the above into account our aims for the simulationspresented in this paper are to show, in agreement withthe analytical predictions of the warped disc instability,that1. discs with warp amplitude above the critical warpamplitude, ψ > ψ c , for extended periods of timelead to instability,2. discs with ψ < ψ c are stable, and3. that the growth rate of the warp amplitude in theunstable regions follow the general trends of thepredicted growth rates, i.e. that the growth ratesare generally higher for smaller viscosity and de-pend on the warp amplitude. NUMERICAL SIMULATIONS It is also worth noting that the disc has ρ → isc tearing R e { s } (cid:7) (cid:8) b = 0.0 (cid:9) = 0.03 q = 1.5, (cid:10) b = 0.03q = 1.6, (cid:11) b = 0.0q = 1.6, (cid:12) b = 0.03q = 1.7, (cid:13) b = 0.0q = 1.7, (cid:14) b = 0.03 (cid:15) (cid:16) = 0.1 (cid:17) (cid:18) = 0.3 Figure 3.
The dimensionless growth rates, ℜ{ s } , plotted against the warp amplitude, ψ , for different values of the viscosityparameter α , the bulk viscosity α b and the shear parameter q . Recall that s = − i ( ω/ Ω)(Ω /c s k ) , and that perturbations grow( ℜ{ s } >
0) or decay ( ℜ{ s } <
0) as exp(
ℜ{− iω } t ). The values are calculated following the method outlined in Doˇgan et al.(2018), and using a code kindly provided by G. I. Ogilvie to calculate the torque coefficients (Ogilvie & Latter 2013). In theleft, middle and right hand panels we plot the values for α = 0 .
03, 0 . . q = 1 . q = 1 . q = 1 .
7. The solid curves represent zero bulk viscosity,while the dashed curves have α b = 0 .
03. The blue dotted line represents zero growth rate, marking the line of instability; forpositive (and real) growth rates the disc becomes unstable. For lower values of the disc viscosity, we see that the disc is unstablefor a broad range of warp amplitude, and the growth rates increase with increasing q (although, it is worth noting this trenddoes not increase to very large q ≈
2; Do˘gan & Nixon 2020). For large disc viscosity, α = 0 .
3, we see that the disc is expectedto be stable. A small level of bulk viscosity does not change the overall picture, but can affect the exact value of the criticalwarp amplitudes and growth rates.
We present three dimensional hydrodynamical sim-ulations of accretion discs with an imposed externaltorque. The fluid dynamics is modelled with thesmoothed particle hydrodynamics technique (SPH; e.g.Price 2012) using the publicly available code phan-tom (Price et al. 2018). This code was first used tomodel warped discs by Lodato & Price (2010), who sim-ulated isolated warped discs with no external torque andfound excellent agreement between the disc evolutionand the evolution predicted by the analytical theory ofOgilvie (1999). Motivated by the broken discs foundin one dimensional calculations of the disc structure byNixon & King (2012), phantom was subsequently usedto model the behaviour of broken discs by Nixon et al.(2012b). Discs with an external torque were modelledin the Lense-Thirring case (Nixon et al. 2012a), the cir-cumbinary case (Nixon 2012; Nixon et al. 2013, see alsoFacchini et al. 2013), and the circumstellar case with abinary companion (Do˘gan et al. 2015). In each case thedynamics proceeds in a similar manner; modest incli-nations lead to a mild warp in the disc and alignmentwith the perturbing torque, while large inclinations gen-erally lead to strong warps and often the disc breaksinto discrete planes which can subsequently precess in-dependently – this process was named as disc tearing byNixon et al. (2012a). Here, we return to the case of discs that are misalignedto the rotation of a spinning black hole, and thus precessdue to the Lense-Thirring effect. As discussed above,we model the central potential with the Einstein poten-tial (e.g. Nelson & Papaloizou 2000) to take account ofthe apsidal precession of disc orbits, and we include theLense-Thirring term through a gravito-magnetic forceterm (e.g. Nelson & Papaloizou 2000). For details of theimplementation see Nealon et al. (2015) and Price et al.(2018).We take the black hole spin to be a = (4 − √ ≈ . R g , where R g = GM/c . The inner edge of the disc, R in , is set to the ISCO radius, and the outer edge ofthe disc is initially at 120 R g . The initial disc surfacedensity follows a power-law that accounts for the zero-torque inner boundary condition as matter accretes onto the black hole, given byΣ( R ) ∝ ( R/R in ) − p (cid:16) − p R in /R (cid:17) , (3)where the normalisation is set by the total disc mass .The disc sound speed follows the locally isothermal ap- As we do not take account of the gas self-gravity in these simu-lations, or allow the back reaction on the black hole spin vector(from accretion or the Lense-Thirring torque), the exact value ofthe disc mass plays no role in the calculation.
A. Raj et. al. proximation with c s ( R ) ∝ ( R/R in ) − q s , (4)with the normalisation fixed by the value of the discangular semi-thickness, H/R , which is specified at R in .For the simulations presented here, we take the power-law indices to be p = 1 and q s = 0 .
5, with the lattergiving a constant disc angular semi-thickness. As dis-cussed above there are two components to the viscos-ity in the simulated disc, a numerical component and aphysical component. For the numerical viscosity we usethe switch proposed by Cullen & Dehnen (2010) withthe linear numerical viscosity coefficient α SPH allowed tovary between 0.01 and unity, and we employ a quadraticnumerical viscosity with β SPH = 2. For the physical vis-cosity, we impose a Navier-Stokes (isotropic) viscosity ofmagnitude given by the Shakura-Sunyaev shear viscosity α (see Section 3.2.4 of Lodato & Price 2010, for details).Unless stated otherwise, the simulations presented be-low employed N p = 10 particles. The simulations areevolved to a time of 10 GM/c , which corresponds toa factor of a few less than the Lense-Thirring (nodal)precession timescale at the outer disc radius (to ensurethe outer regions are undisturbed on the timescale ofthe simulation) and is long enough to ensure the innerdisc regions have evolved significantly (the run time isequivalent to a Lense-Thirring precession timescale at aradius of approximately 50 R g ).In the following sections we present quantities fromthe simulations that require averaging over the SPH par-ticle data. For example, the warp amplitude at eachradius requires taking the radial derivative of the unittilt (orbital angular momentum) vector for the disc. Tocalculate the angular momentum vector as a functionof radius in the disc, we must take an average over aselection of the disc particles. For the plots presentedbelow, we split the disc into N shells spaced logarithmi-cally in radius. For each shell, with spherical radius r s ,we average over the particles that have | r − r s | < ξH ,where H is the disc scale-height at r s . We find thatchoosing N = 150 and ξ = 2 serves to minimise Poissonnoise from the discreteness of the particles, while notoversmoothing features in the disc.We investigate the warped disc dynamics by vary-ing the following parameters; we simulate discs withtwo thicknesses ( H/R = 0 . , . θ = 10 ◦ , ◦ & 60 ◦ ), and three values of the disc vis-cosity ( α = 0 . , . . Varying the disc thickness
The disc thickness is expected to affect the evolutionof warped discs. For a disc subject to an external torque,the disc thickness affects the magnitude of the effec-tive viscosities in the disc (which are proportional to(
H/R ) ) and thus affects the location at which the ex-ternal torque becomes dynamically important. For ex-ample, the disc tearing radius derived by Nixon et al.(2012a) depends on the disc thickness, with a largerextent of the disc expected to be vulnerable to insta-bility if the disc is thinner. Another effect that arisesfrom the disc thickness is that it alters the lengthscaleon which the disc can be warped. It is not possible tobend a disc on a lengthscale shorter than ∼ H , andthus when the same torque is applied thick discs cannotachieve as large a warp amplitude as thin discs. Thisimplies that thick discs are less vulnerable to instabilityas they are less likely to achieve a warp amplitude that isgreater than the critical warp amplitude (Doˇgan et al.2018). We demonstrate this by plotting in Fig. 4 thewarp amplitude with radius for the simulations with α = 0 . H/R = 0 .
02 and
H/R = 0 .
05 (with the left hand panel at θ = 10 ◦ , middlepanel at 30 ◦ and the right hand panel at θ = 60 ◦ ). Inthese cases the discs are all coherently warped, exceptfor the thinnest and most inclined case ( H/R = 0 . θ = 60 ◦ ) which exhibits two breaks at R/R g ≈ R/R g ≈
20. We note that the plots are made after atime of 5 × GM/c , which is halfway through the sim-ulation. This corresponds to a Lense-Thirring precessiontimescale at R ≈ R g . The main results depicted bythese plots is that (1) thicker discs have a peak in thewarp amplitude that occurs at a smaller radius thanthinner discs (see also, for example, Kumar & Pringle isc tearing , and that (2) in gen-eral (but not always ) the peak warp amplitude ishigher for thinner discs and this occurs across a broadsection of the disc. For the right hand panel of Fig. 4we can see that for the same viscosity parameter anddisc inclination, the disc thickness can play a key rolein determining how vulnerable the disc is to instability;this can occur both due to the increased warp amplitude,and also due to the increased local viscous timescale, as-sociated with thinner discs. We also note that for theseparameters, the thinner disc ( H/R = 0 .
02) is aligned tothe black hole spin plane at the ISCO, but that this isnot true for the thicker case (
H/R = 0 .
05) where, whilethere is a downturn in the inclination, typically the discis still significantly misaligned at the ISCO.3.2.
Varying the disc inclination
The disc inclination has an effect on the disc evolu-tion as, for larger inclinations, the warp amplitude istypically higher. This is most easily seen for low incli-nations where the evolution is essentially linear, in thatthe behaviour can be rescaled by the inclination angle θ (cf. Lubow et al. 2002). However, once the inclina-tion angle becomes large enough (typically greater thana few times H/R , which for our parameters is & ◦ )then the evolution becomes noticeably non-linear in theinclination angle. This is evident in Fig. 5 where the discbehaviour shows a strong dependence on inclination an-gles between θ = 10 ◦ , 30 ◦ and 60 ◦ . This figure showsthe warp amplitudes, with the same format as Fig. 4,but this time for the simulations with α = 0 .
03; theleft hand panel shows
H/R = 0 .
02 while the right handpanel shows
H/R = 0 .
05 and in both cases the black-solid(red-dashed)[green-long-dashed] line corresponds to θ = 10 ◦ (30 ◦ )[60 ◦ ]. For the lowest inclination of 10 ◦ thewarp amplitude remains low, and the disc attains a co-herent warped shape. This is also true for the 30 ◦ casewhen H/R = 0 .
05. For the other cases, the disc be- It is worth noting that the fact that the warp amplitude canpeak near the inner edge of the disc is of interest, as it allowsfor the possibility that the disc can produce time-dependent, re-peating, behaviour there. If the warp amplitude is maximal nearthe inner edge, and the value of the warp amplitude is above thecritical warp amplitude there (cf. Do˘gan & Nixon 2020), then theinnermost ring of the disc may tear off, precess and be rapidlyaccreted, and this process may repeat indefinitely (or at least un-til the disc conditions change, causing the warp amplitude thereto evolve). We discuss one example of this from our simulationsin more detail in Raj & Nixon (2021). For example, the simulations corresponding to the left hand panelof Fig. 4 but with α = 0 .
03 (not depicted) show a higher peakwarp amplitude in the thicker case. For these cases, the peak inthe thicker case occurs at a small radius, at which the thinnerdisc has already aligned to the black hole spin. comes sufficiently warped that the warp amplitude risesabove the critical warp amplitude (Doˇgan et al. 2018)and the disc tears into discrete rings (the rapid changein radius of the disc orbital angular momentum vectoris indicated in the warp amplitude by the sharp peakspresent in Figs 4 & 5). These results reinforce the con-clusion from the previous section (3.1) as the thinnerdiscs are subject to more vigorous tearing, and also con-firm the expectation that larger inclination results inhigher warp amplitude and a higher susceptibility to in-stability. We note that to first order the Lense-Thirringprecession that we simulate here has a frequency thatis independent of the inclination angle (i.e. the torqueis proportional to sin θ ), while the precession frequencyinduced by, for example, the gravitational field of a mis-aligned binary system (see, e.g., Bate et al. 2000) hasan angular dependence such that the frequency goes tozero at θ = 90 ◦ (the torque in this case is proportionalto sin 2 θ ). This means that for discs in binary systemsthe increase in warp amplitude with increasing inclina-tion angle peaks at 45 ◦ (and also 135 ◦ ), rather than at90 ◦ in the Lense-Thirring case.3.3. Varying the disc viscosity
The disc viscosity parameter α plays a crucial rolein the dynamics of warped discs. For α < H/R thedisc response to warping is to propagate bending waves,while for α > H/R the disc responds by diffusing thewarp (Papaloizou & Pringle 1983). In the diffusive case,which we focus on here (although note that some ofour simulations are potentially wavelike with α slightlysmaller than H/R ), the exact value of α plays two dis-tinct roles in the dynamics.The first is that the torque coefficients for modestwarp amplitudes vary with α . In the limit of smallwarps and low α , the usual ‘planar’ viscosity is pro-portional to α , while the ‘vertical’ viscosity is propor-tional to 1 /α (Papaloizou & Pringle 1983) . Thus forsmall values of α , there is a significant discrepancy be-tween the timescales associated with the disc resistinga warp by attempting to straighten itself (governed bythe vertical viscosity) and the timescale on which massflows inwards through the disc (governed principally bythe planar viscosity). However, when α is large andapproaches unity, these timescales become comparable.This means that large α makes it difficult to break a Here we have assumed that the underlying mechanism generatingthe viscosity is essentially isotropic such that any shear is dampedat the same rate (given by α ), such as might be expected whenthe viscosity arises due to small scale turbulence. For a discus-sion of the viscosity coefficients in the case that the underlyingmechanism generates an anisotropic viscosity, see Nixon (2015). A. Raj et. al.
R(cid:19)(cid:20) g ψ α = 0.1, H/R = 0.02, θ = 10 o t = 5x10 GM/c α = 0.1, H/R = 0.05, θ = 10 o R/R g ψ
50 10000.511.5 α = 0.1, H/R = 0.02, θ = 30 o t = 5x10 GM/c α = 0.1, H/R = 0.05, θ = 30 o R/R g ψ
50 10005101520 α = 0.1, H/R = 0.02, θ ’ ()o t = 5x10 GM/c α = 0.1, H/R = 0.05, θ * +,. Figure 4.
The warp amplitude, ψ = R | ∂ l /∂R | , plotted against radius, R/R g , at a time of 5 × GM/c for the simulations with α = 0 .
1. Each panel compares the warp amplitude in the simulation with
H/R = 0 .
02 (black solid line) with the correspondingsimulation with
H/R = 0 .
05 (red dashed line), with the left(middle)[right] hand panel showing the simulation with the discinitially inclined by 10 ◦ (30 ◦ )[60 ◦ ]. In each panel the range of warp amplitude that is plotted is varied to show the full extentof the solution. In each case the peak of the warp amplitude moves to smaller radius as the disc thickness is increased, but thepeak warp amplitude is typically higher when the disc is thinner. For the simulation with H/R = 0 .
02 and θ = 60 ◦ , there aretwo clear breaks in the disc at R/R g ≈ R/R g ≈
20, and these are typically accompanied by deficits of warp amplitude inneighbouring regions. This figure also shows that the critical warp amplitude at which the disc breaks is larger than unity for α = 0 . θ = 30 ◦ . For θ = 60 ◦ , the disc typically becomes unstable when ψ &
2, consistent with ψ c ≈ .
25 in this case.
R/R g ψ
50 100024 α = 0.03, H/R = 0.02, θ = 10 o t = 5x10 GM/c α = 0.03, H/R = 0.02, θ = 30 o α = 0.03, H/R = 0.02, θ / 234 R/R g ψ
50 100024 α = 0.03, H/R = 0.05, θ = 10 o t = 5x10 GM/c α = 0.03, H/R = 0.05, θ = 30 o α = 0.03, H/R = 0.05, θ Figure 5.
The warp amplitude, ψ = R | ∂ l /∂R | , plotted against radius, R/R g , at a time of 5 × GM/c for the simulationswith α = 0 .
03. Each panel compares the warp amplitude in the simulation with an inclination of θ = 10 ◦ with the correspondingsimulation with θ = 30 ◦ and θ = 60 ◦ , with the left(right) hand panel showing the simulations with disc thickness of H/R =0 . . α = 0 .
03 the critical warp amplitude is ψ c ≈ . disc; this is because in this case the timescale on whichthe disc can internally rearrange itself, and potentiallylocally flatten the disc into two (or more) distinct pla-nar regions, is similar to the timescale on which matterflows from one region of the disc to another; this lattereffect can keep the disc as a coherent whole. Whereasfor smaller α the radial communication in the disc canbe severely restricted at large warp amplitudes, allowingdiscrete parts of the disc to behave independently (seeNixon & King 2012, for discussion).The second role that the disc viscosity plays in warpeddisc dynamics is more subtle. Ogilvie (1999) showedthat the torque coefficients vary as a function of thewarp amplitude in the disc, and that the way in which they vary also depends on α . Ogilvie (2000) showedthat the variation in torque coefficients with warp ampli-tude causes an instability, and Doˇgan et al. (2018) andDo˘gan & Nixon (2020) explored this instability and con-nected it with the disc breaking/tearing that had beenfound in numerical simulations. These analyses showthat the instability is generally more vigorous (lowercritical warp amplitudes and higher growth rates) forlower values of α .In Fig. 6 we show the effect of varying the value of α inthe simulations. In the left hand panel we show the sim-ulations with H/R = 0 .
02 and θ = 30 ◦ for α = 0 .
03, 0 . .
3. In this case we can clearly see the non-linear ef-fect of changing the disc viscosity on the disc warp. For isc tearing α the disc exhibits a smooth warpedshape with a peak at R/R g ≈
30. As α is decreasedthe general shape across most of the disc remains thesame, and the peak warp amplitude occurs in a similarradial location. However, for lower α the peak warp am-plitude grows. This is a modest increase when α variesfrom 0 . .
1, but reducing further to α = 0 .
03 rendersthe disc unstable and the peak warp amplitude increasessharply because of this. In the right hand panel we showthe simulations with
H/R = 0 .
05 and θ = 30 ◦ for thesame range of α . In this case none of the simulationsreach the critical warp amplitude and they all remainstable. The peak warp amplitude moves inwards slightlywhen α is decreased from 0 . . α is more likely tolead to instability. It is worth reminding the reader herethat the values of viscosity here are the values for thephysical viscosity term in the simulations, and that thetotal viscosity present also includes the numerical vis-cosity that we discussed in detail in Section 2.It is also worth noting that we detect the presence ofpropagating waves in the simulations with α = 0 .
03 and
H/R = 0 .
05 (i.e. with α . H/R ; see the low ampli-tude waves in the warp amplitude at 30 . R/R g . R/R g &
90 in theright hand panel of Fig. 6), but that these effects arenoticeably reduced by the time α & H/R (e.g. compar-ing the results with α = 0 .
03 and
H/R = 0 .
02 to thosewith α = 0 . H/R = 0 . α is increasedfrom below to above H/R may be quite sharply focussedaround α = H/R (cf. Martin et al. 2019b). This is aninteresting question that warrants a focussed investiga-tion in the future. 3.4.
Disc tearing
In this sub-section we take a more detailed look atthe properties of the simulated discs in which instabilityoccurs. The simulations which exhibit clear disc tear-ing, and thus which we analyse in this section are thefollowing:A: α = 0 . H/R = 0 . θ = 30 ◦ ,B: α = 0 . H/R = 0 . θ = 60 ◦ ,C: α = 0 . H/R = 0 . θ = 60 ◦ , andD: α = 0 . H/R = 0 . θ = 60 ◦ .For these simulations we note that B is the same as A,but at higher inclination. C is the same as B, but with a thicker disc. And D is the same as B, but with a higherviscosity. One of the other simulations, with α = 0 . H/R = 0 . θ = 30 ◦ , shows an early break, whichsubsequently decays after propagating outwards throughthe disc, so we do not include this in the analysis. Wealso find instability in the case ofE: α = 0 . H/R = 0 . θ = 60 ◦ .As this case is anomalous from the point of view ofthe local stability analysis presented by Doˇgan et al.(2018) we discuss it separately in Section 3.5 below.The remaining simulations all warp in response to Lense-Thirring precession, but do not exhibit any instability.We discuss in turn each of the simulations that showclear examples of disc tearing behaviour (listed A-Eabove). In each case we identify two regions of the discwhich exhibit a disc break caused by disc tearing. Welook for clean examples where the disc is coherent be-forehand and results in a single clear break (at leastlocally to that patch of the disc). For example, we havenot chosen regions where a peak in warp amplitude ul-timately splits into several distinct breaks, or where thebreak occurs in a region where there is already a com-plex shape to the warp amplitude. In such regions weexpect the effects of different growing modes would com-plicate matters and make our analysis less clear. For thetwo breaks examined, in each case we measure the ra-dial range ( R , R ) within which the peak of the warpamplitude resides while the peak warp amplitude growsfrom ψ at time t to ψ at time t . We then calculatethe time evolution of the peak warp amplitude in thisregion ψ p ( t ). In unstable regions we expect the warpamplitude to grow as ψ ( t ) ∼ ψ exp h ℜ{ s } Ω ( Hk ) ( t − t ) i , (5)where Ω is the orbital frequency, H ≡ c s / Ω, k is thewavenumber, ℜ{ s } is the dimensionless growth rate ofthe instability. In each unstable case, we find that theevolution of the peak warp amplitude with time is ap-proximately exponential while the peak warp amplitudeis in the range 2 . − Note that when searching for the peak warp amplitude we findthat the 150 logarithmically spaced radial points used in the pre-vious sections is insufficient. For these calculations we employ alinear spaced grid with 1200 points. A. Raj et. al.
R/R g ψ
50 100024 α = 0.03, H/R = 0.02, θ = 30 o t = 5x10 GM/c α = 0.10, H/R = 0.02, θ = 30 o α = 0.30, H/R = 0.02, θ = 30 o R/R g ψ
50 10000.20.4 α = 0.03, H/R = 0.05, θ = 30 o t = 5x10 GM/c α = 0.10, H/R = 0.05, θ = 30 o α = 0.30, H/R = 0.05, θ = 30 o Figure 6.
The warp amplitude, ψ = R | ∂ l /∂R | , plotted against radius, R/R g , at a time of 5 × GM/c for the simulationswith θ = 30 ◦ . Each panel compares the warp amplitude in the simulations with varying values of the Shakura-Sunyaev viscosityparameter of α = 0 .
03, 0 . .
3. The left hand panel shows the results for the thinner case of
H/R = 0 .
02 and the righthand case shows the thicker case of
H/R = 0 .
05. The left hand panel shows that decreasing the viscosity can destabilise thedisc, with the peak warp amplitude increasing until it becomes larger than the critical value. The warp amplitude at R ≈ R g in the higher viscosity cases is ≈ . −
1. For α = 0 .
03, the critical warp amplitude ψ c ≈ .
5, and thus the lower viscosity casebecomes unstable when ψ & ψ c as expected. The right hand panel shows that for stable discs, as the viscosity is decreased thesolutions can transition from diffusive to wavelike in character (as shown by Papaloizou & Pringle 1983). estimate the growth rate, we rearrange (5) to give ℜ{ s } ( Hk ) = ln( ψ /ψ )Ω∆ t , (6)where ∆ t = t − t is the time taken to grow from ψ =2 . ψ = 4. We find that in each case the peak of thewarp amplitude does not move significantly while theamplitude grows between these two bounds, allowingus to use a single value for Ω, which for the Einsteinpotential we employ is given byΩ = GMR (cid:18) R g R (cid:19) . (7)Ideally, we would also like to measure the wavenum-ber, k , from the simulations so that we could comparethe value of the dimensionless growth rate directly tothe predicted value. However, there does not appearto be an accurate way to measure this directly fromthe simulation data. Future simulations which insert agrowing mode with a single, well-defined value for k aredesirable. In the dynamic simulations we present hereit is likely that there is significant power at a range ofwavenumbers. In this case we can expect growth to beslower initially, until sufficient power is in the fastestgrowing modes. For the mode to fit in the disc (recallthat we cannot bend the disc on scales . H ) we musthave Hk . Hk less than, but not muchless than, unity. However, to be consistent we report thevalues of χ ≡ ℜ{ s } ( Hk ) which can be measured fromthe simulation data.For simulation A with α = 0 . H/R = 0 .
02 and θ = 30 ◦ we find that the warp amplitude in the inner disc regions rapidly grows to greater than the criticalvalue for R . R g (which for this viscosity is ψ c ≈ . ≈ R g . Outsideof this region ψ remains below the critical value at alltimes, and no breaks are observed in this region. Forthis simulation we analyse two clear breaks, the firstoccurring between R ≈ − R g with a peak in thewarp amplitude at R p ≈ R g while the warp amplitudeis growing from ψ = 2 . ψ = 4. The second breakforms an initial peak in the warp amplitude at R ≈ R g and by the time the peak warp amplitude is ≈ R p ≈ R g . For thefirst break, we measure ∆ t ≈ GM/c )and with the peak at R p ≈ R g we have Ω ≈ . χ ≈ . t ≈ R p ≈ R g , Ω ≈ × − and therefore χ ≈ . . − . ℜ{ s } = 0 . − .
5, thiswould imply k − ≈ (3 − H which seems reasonable,but we also note that it is likely that the total level ofviscosity in the simulation is slightly higher than α =0 .
03 (see discussion in Section 2) and thus we wouldexpect the growth rate to be reduced somewhat (i.e. if α ≈ .
05, then the predicted dimensionless growth rateimplies k − ≈ (2 − H ).For simulation B with α = 0 . H/R = 0 .
02 and θ = 60 ◦ we analyse two breaks in the disc, with the firstoccurring at R p ≈ R g and the second which begins at isc tearing R ≈ R g and settles into a strong peak at R p ≈ R g .For the first break, we measure ∆ t ≈ R p ≈ R g we have Ω ≈ .
01. This yields χ ≈ .
03. For the second peak we have ∆ t ≈ R p ≈ R g , Ω ≈ × − and therefore χ ≈ .
02. Asbefore, when scaled to the local dynamical timescale thegrowth rate in the unstable region is similar between thetwo breaks in the same disc, and are also similar to theprevious simulation which has the same α value.For simulation C with α = 0 . H/R = 0 .
05 and θ = 60 ◦ we also analyse two breaks in the disc, withthe first occurring between R ≈ − R g with thepeak typically at R p ≈ R g and the second occurringbetween R ≈ − R g with the peak typically at R p ≈ R g . For the first break, we measure ∆ t ≈ R p ≈ R g we have Ω ≈ × − . Thisyields χ ≈ .
02. For the second peak we note that thegrowth in warp amplitude slows appreciably for ψ & . ψ = 2 . ψ = 3 .
5, so in this instance we explore this range. Herewe have ∆ t ≈ R p ≈ R g , Ω ≈ . × − andtherefore χ ≈ . H/R = 0 .
05, andin this case the radial extent of the peak of the warpamplitude is broader than in the thinner disc simulationsby a factor of approximately 2.5 which demonstratesthat the break in the disc occurs across a radial lengthscale of a few H independent of the value of H/R .For simulation D with α = 0 . H/R = 0 .
02 and θ = 60 ◦ we again analyse two breaks in the disc, withthe first occurring at R p ≈ R g and the second oc-curring at R p ≈ R g . For the first break, we measure∆ t ≈ R p ≈ R g we haveΩ ≈ .
01. This yields χ ≈ .
03. For the second peak wehave ∆ t ≈ R p ≈ R g , Ω ≈ × − and there-fore χ ≈ .
02. Therefore we find that the growth ratesin this case, with α = 0 .
1, are similar to the growth ratesfound for α = 0 .
03. The predictions from the stabilityanalysis suggest that the growth should be slower in thiscase. This may indicate that for the lower viscosity sim-ulations ( α = 0 .
03) the numerical viscosity is impedingthe growth rate somewhat. Alternatively it might bethat initially growth occurs for small k (on large spa-tial scales) and thus the growth is slow until enoughpower is transferred into large k modes (smaller lengthscales) at which point the instability is sufficiently rapidthat the timescale is essentially independent of disc pa-rameters once the disc becomes unstable. We note that χ ≈ .
02 implies that the instability growth timescaleis ∼ / Ω (i.e., several orbits). It seems unreasonableto expect the instability to grow faster than this. We speculate that this second reason is what is occurringin the simulations as we find that analysis of two of thebreaks in the simulation with α = 0 . H/R = 0 .
02 and θ = 60 ◦ when performed at a resolution correspondingto N p = 10 both also yield χ ≈ .
02, hinting that thegrowth rate (which encompasses the wavenumber of thefastest growing mode) may not depend strongly on res-olution.Finally, we also note that simulation C, with α = 0 . H/R = 0 .
05 and θ = 60 ◦ , performed at a resolutioncorresponding to N p = 10 showed some interesting dy-namics that was not present in any of the other simula-tions. In this simulation the peak of the warp amplitudeis concentrated at small radii. There are two additionalpeaks of warp amplitude, which occur at larger radii.In the N p = 10 simulation, these peaks at larger radiiare resolved into breaks which then alter the inner discevolution, while in the N p = 10 case the outer discremained coherent. This allowed the inner regions, onlonger timescales, to develop a break near the inner edgewhich is quickly accreted. This process repeats with theinnermost ring being broken off, it subsequently pre-cesses, and it is then accreted. It is possible that forthe right combination of parameters, and variation ofthose parameters with radius, that this evolution couldmanifest in a physical disc. For example, if the disc an-gular semi-thickness were a slowly increasing functionof radius, then the disc may only be unstable near theinner disc edge (and this could remain so as resolutionis increased). We discuss this manifestation of the disctearing dynamics, and the implications this may havefor the observational properties of these discs, in moredetail in Raj & Nixon (2021).3.5. Instability at high viscosity
We noted above that simulation E, with α = 0 . H/R = 0 .
02 and θ = 60 ◦ , showed instability. Thiswas not expected as the local stability analysis of thewarped disc equations suggests that the disc is stablefor any value of the warp amplitude at this high a vis-cosity. However, instability in this case had previouslybeen reported by Nixon & King (2012) for discs subjectto the Lense-Thirring precession using a 1D ring code toevolve the warped disc evolution equations. In this caseNixon & King (2012) showed that the disc can breakinto two distinct planes for α ≈ . − .
3. They presentevolution that is distinct from the evolution presentedby most hydrodynamical simulations in that the breaksmanifest as a sharp variation in the disc tilt, while thebreaks in most hydrodynamical simulations manifest (atleast initially) as a sharp variation in the disc twist. Herewe find that for α = 0 .
3, the break in the disc forms with4
A. Raj et. al. a large variation in the disc tilt, consistent with the re-sults in Nixon & King (2012).For comparison with the growth rates reported above,we again analyse two breaks in the disc for this simula-tion with α = 0 .
3. The first occurring at R p ≈ R g andthe second occurring at R p ≈ R g . For the first break,we measure ∆ t ≈
950 and with the peak at R p ≈ R g we have Ω ≈ .
02. This yields χ ≈ .
02. For the secondpeak we have ∆ t ≈ R p ≈ R g , Ω ≈ .
015 andtherefore χ ≈ . α = 0 .
3, are similar to the growthrates found for the lower α cases.There are several possibilities for why the simulateddisc with α = 0 . α for some values of the shear rate q and bulk viscosity α b that was not revealed in the analy-ses presented in Doˇgan et al. (2018) and Do˘gan & Nixon(2020). It is possible that the disc tearing presented inthese simulations is a numerical artefact, but again weconsider this unlikely given the broad range of codesand methods that have now reported such behaviour instrongly warped discs. We speculate that the answeris related to the applicability of the warped disc equa-tions (and in particular the accuracy of calculating thetorque coefficients) and/or the local stability analysis atthe extremes of the parameter space. Here we have largeviscosity and large warp amplitudes, and the disc struc-ture is strongly affected by the external torque. It seemspossible that one or more of the assumptions of the the-ory of warped discs could be violated, or that non-localeffects become important. For example, as discussed byNixon & King (2012), the torques responsible for hold-ing the disc together generally weaken with increasingwarp amplitude. This is in part due to the drop in sur-face density in the region of large warp amplitude. Oncethe surface density has dropped sufficiently far, it maybe that there is simply not a strong enough torque leftto hold the disc together in that region. If the stabilityof the disc is dependent on such a process, in which thelocal mass is slowly reduced over time, this may not becaptured by a local stability analysis which assumes astable background state. CONCLUSIONSWe have presented and discussed the results of numer-ical simulations of discs that are warped by the Lense-Thirring effect from a spinning black hole. In generalwe find that the results of our simulations (Section 3)are consistent with the predictions of the local stabil-ity analysis of the warped disc equations (Doˇgan et al. 2018). We confirm that disc tearing occurs preferen-tially in thin, low viscosity, and highly inclined discs.Thin and highly inclined discs are more likely to yieldlarge warp amplitudes, and low viscosity discs are gen-erally expected to become unstable at lower warp am-plitudes. We show that instability arises in the simu-lated discs at higher values of the warp amplitude forlarger viscosity parameters, and the critical warp am-plitudes are commensurate with the predicted values.We also find that the growth rates of the instabilityare consistent with the predicted growth rates, assum-ing that the growing modes have wavenumbers, k , with kH ∼ . − .
5. We also find some unexpected results.In particular we found instability for a simulation with α = 0 . α value should be stable. We note that wedo find that the results of the simulations vary some-what with increasing resolution; typically we find thatfor stable discs the results are converged for the bulkof the disc with small differences near the inner edgewhere the density, and thus resolution, is lower, but wedo find that higher resolution simulations exhibit moredisc breaks that are sharper (i.e. have a larger peak warpamplitude; see Appendix A). We attribute this to theincrease in spatial resolution which allows smaller phys-ical features to be resolved in the simulations, and tothe decreasing influence of numerical viscosity at higherresolution. Returning to the three aims that we listedat the end of Section 2, namely1. discs with warp amplitude above the critical warpamplitude, ψ > ψ c , for extended periods of timelead to instability,2. discs with ψ < ψ c are stable, and3. that the growth rate of the warp amplitude in theunstable regions follow the general trends of thepredicted growth rates, i.e. that the growth ratesare generally higher for smaller viscosity and de-pend on the warp amplitude.Our simulations support (1) and (2), showing that theinstability seen in the simulations is consistent with thepredicted critical warp amplitudes. However, for (3), isc tearing splash (Price 2007) for the figures.REFERENCES Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Bardeen, J. M., & Petterson, J. A. 1975, ApJL, 195, L65Bate, M. R. 2018, MNRAS, 475, 5618Bate, M. R., Bonnell, I. A., Clarke, C. J., et al. 2000,MNRAS, 317, 773Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS,401, 1505Casassus, S., Marino, S., P´erez, S., et al. 2015, ApJ, 811, 92Cullen, L., & Dehnen, W. 2010, MNRAS, 408, 669Do˘gan, S., Nixon, C., King, A., & Price, D. J. 2015,MNRAS, 449, 1251Do˘gan, S., & Nixon, C. J. 2020, MNRAS, 495, 1148Doˇgan, S., Nixon, C. J., King, A. R., & Pringle, J. E. 2018,MNRAS, 476, 1519Facchini, S., Lodato, G., & Price, D. J. 2013, MNRAS, 433,2142Flebbe, O., Muenzel, S., Herold, H., Riffert, H., & Ruder,H. 1994, ApJ, 431, 754Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson,J. D. 2007, ApJ, 668, 417Fragner, M. M., & Nelson, R. P. 2010, A&A, 511, A77Gammie, C. F. 2004, ApJ, 614, 309Greenhill, L. J., Jiang, D. R., Moran, J. M., et al. 1995,ApJ, 440, 619Jeter, B., & Broderick, A. E. 2020, arXiv e-prints,arXiv:2010.11303Kato, S. 1990, PASJ, 42, 99 King, A. R., Pringle, J. E., & Livio, M. 2007, MNRAS, 376,1740Kraus, S., Kreplin, A., Young, A. K., et al. 2020, Science,369, 1233Kumar, S., & Pringle, J. E. 1985, MNRAS, 213, 435Larwood, J. D., Nelson, R. P., Papaloizou, J. C. B., &Terquem, C. 1996, MNRAS, 282, 597Larwood, J. D., & Papaloizou, J. C. B. 1997, MNRAS, 285,288Liptai, D., & Price, D. J. 2019, MNRAS, 485, 819Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2020,MNRAS, arXiv:1904.08428Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2002,MNRAS, 337, 706Martin, R. G., Nixon, C. J., Pringle, J. E., & Livio, M.2019a, NewA, 70, 7Martin, R. G., Lubow, S. H., Pringle, J. E., et al. 2019b,ApJ, 875, 5Meru, F., & Bate, M. R. 2012, MNRAS, 427, 2022Miyoshi, M., Moran, J., Herrnstein, J., et al. 1995, Nature,373, 127Murray, J. R. 1996, MNRAS, 279, 402Natarajan, P., & Pringle, J. E. 1998, ApJL, 506, L97Nealon, R., Price, D. J., & Nixon, C. J. 2015, MNRAS, 448,1526Nelson, R. P., & Papaloizou, J. C. B. 1999, MNRAS, 309,929 A. Raj et. al. —. 2000, MNRAS, 315, 570Nixon, C. 2015, MNRAS, 450, 2459Nixon, C., & King, A. 2016, Lecture Notes in Physics, 905,45Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946Nixon, C., King, A., Price, D., & Frank, J. 2012a, ApJL,757, L24Nixon, C., & Salvesen, G. 2014, MNRAS, 437, 3994Nixon, C. J. 2012, MNRAS, 423, 2597Nixon, C. J., & King, A. R. 2012, MNRAS, 421, 1201Nixon, C. J., King, A. R., & Price, D. J. 2012b, MNRAS,422, 2547Ogilvie, G. I. 1999, MNRAS, 304, 557—. 2000, MNRAS, 317, 607Ogilvie, G. I., & Latter, H. N. 2013, MNRAS, 433, 2403Papaloizou, J. C. B., & Pringle, J. E. 1983, MNRAS, 202,1181 Papaloizou, J. C. B., & Terquem, C. 1995, MNRAS, 274,987Price, D. J. 2007, PASA, 24, 159—. 2012, Journal of Computational Physics, 231, 759Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA,35, e031Pringle, J. E. 1992, MNRAS, 258, 811Raj, A., & Nixon, C. J. 2021, ApJ in pressScott, D. M., Leahy, D. A., & Wilson, R. B. 2000, ApJ,539, 392Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Sorathia, K. A., Krolik, J. H., & Hawley, J. F. 2013, ApJ,768, 133 isc tearing R/R g ψ
50 10000.10.20.30.40.5 Np = 10 t : ;<=>? GM/c @A = 10 BC = 10 R/R g ψ
50 10000.511.5 DE = 10 F G HIJKL
GM/c MO = 10 PQ = 10 R/R g ψ
50 1000510 ST = 10 U V WXYZ[
GM/c \] = 10 ^_ = 10 Figure 7.
The warp amplitude, ψ = R | ∂ l /∂R | , plotted against radius, R/R g , at a time of 5 × GM/c for the simulations atdifferent inclinations (left to right) with α = 0 . H/R = 0 .
02. Each panel compares the warp amplitude in the simulationsvarying the resolution as given in the legend. The left hand panel shows the results for the inclination of 10 ◦ , the middle panelshows the inclination of 30 ◦ and the right hand panel shows the inclination of 60 ◦ . At low inclinations the disc is stable, andthe results are similar at high particle number. The same is true at inclination of 30 ◦ . At 60 ◦ the disc is unstable for theseparameters, and the solutions all exhibit this but the details are different at different resolutions as discussed in the text. APPENDIX A. RESOLUTIONWe have discussed the effects of resolution in detail in Section 2, and throughout. There are two principle effectsthat numerical resolution has on the results of warped disc simulations. First, higher spatial resolution allows for theformation of sharper features in the disc. For example, an underresolved region of high warp amplitude may result inthe warp amplitude being underestimated and thus the disc remaining (artificially) stable when (physically) the discshould be unstable. Secondly, the numerical viscosity in SPH simulations is a function of the local resolution with thelinear numerical viscosity proportional to the resolution lengthscale and the quadratic numerical viscosity proportionalto the square of the resolution lengthscale (cf. equation 2). As the resolution is increased the contribution to thetotal viscosity from the numerical viscosity decreases, and thus for a fixed physical viscosity the total viscosity in thesimulation decreases. As discussed in Section 2 we expect the magnitude of the numerical viscosity to be of the orderof, but smaller than, the lowest physical viscosity we simulated ( α = 0 . α = 0 .
3. Bothof these effects combine at the inner edge of the disc in SPH simulations where the surface density goes to zero atthe inner boundary. Here, spatial gradients in the disc (for example of the surface density) are typically largest andmay be underresolved by the available spatial resolution, and the decreasing spatial resolution associated with thedecreasing density of the disc results in an increased viscosity. As resolution is increased the solution tends towardsthe expected solution (e.g. equation 3 for the surface density of a planar disc). Finally, we also noted that becausethe instability, including both the critical warp amplitude and the growth rates (see Doˇgan et al. 2018, for details),depends sensitively on the viscosity parameter α , we cannot expect numerical simulations to show converged resultsfor the precise details of unstable regions, for example the exact location at which the instability manifests. To showan example of the differences in the simulation results at different resolutions we plot the warp amplitudes for thesimulations at different inclination angles with α = 0 . H/R = 0 ..