Table Of Content

Mon.Not.R.Astron.Soc.000,1–15(2014) Printed9January2015 (MNLATEXstylefilev2.2) On the Bardeen-Petterson effect in black hole accretion discs Rebecca Nealon1⋆, Daniel J. Price1† & Chris J. Nixon2‡§ 5 1MonashCentreforAstrophysicsandSchoolofMathematicalSciences,MonashUniversity,Vic3800,Australia 1 2JILA,UniversityofColorado&NIST,BoulderCO80309-0440,USA 0 2 n a J ABSTRACT 7 Weinvestigatetheeffectofblackholespinonwarpedormisalignedaccretiondiscs—inpar- ticulari)whetherornottheinnerdiscedgealignswiththeblackholespinandii)whetherthe ] disc can maintaina smoothtransitionbetweenan alignedinnerdisc anda misalignedouter E disc,knownastheBardeen-Pettersoneffect.Weemployhighresolution3Dsmoothedparti- H clehydrodynamicssimulationsofα-discssubjecttoLense-Thirringprecession,focussingon h. thebendingwaveregimewherethediscviscosityissmallerthantheaspectratioα . H/R. p Wefirstaddressthecontroversyintheliteratureregardingpossiblesteady-stateoscillationsof - thetiltclosetotheblackhole.Wesuccessfullyrecoversuchoscillationsin3Datbothsmall o andmoderateinclinations(.15◦),providedbothLense-ThirringandEinsteinprecessionare r t present,sufficientresolutionisemployed,andprovidedthediscisnotsothicksoastosimply s accretemisaligned.Second,wefindthatdiscsinclinedbymorethanafewdegreesingeneral a [ steepenandbreakratherthanmaintainasmoothtransition,againincontrasttopreviousfind- ings,butonlyoncethediscscaleheightisadequatelyresolved.Finally,wefindthatwhenthe 1 discplaneismisalignedtotheblackholespinbyalargeangle,thedisc‘tears’intodiscrete v ringswhichprecesseffectivelyindependentlyandcauserapidaccretion,consistentwithpre- 7 viousfindingsinthediffusiveregime(α & H/R).Thusmisalignmentbetweenthediscand 8 thespinaxisoftheblackholeprovidesarobustmechanismforgrowingblackholesquickly, 6 1 regardlessofwhetherthediscisthickorthin. 0 Keywords: accretion,accretiondiscs—blackholephysics—hydrodynamics—galaxies: . 1 jets—(galaxies:)quasars:supermassiveblackholes 0 5 1 : v 1 INTRODUCTION Kumar&Pringle (1985) and Pringle (1992) were performed in i X the diffusive regime. The first studies performed in the bending Bardeen&Petterson (1975) first computed the evolution of a waveregime(Ivanov&Illarionov1997;Lubow,Ogilvie&Pringle r warped accretion disc subjected to Lense-Thirring precession a 2002, hereafter LOP02) suggested a conflict with the Bardeen- (Lense&Thirring1918)causedbyframe-dragging fromthespin Pettersonpicture—findingthattheblackholespincoulddrivethe ofacentralblackhole.Althoughtheiroriginalequationswerelater disctiltintoasteadystatethatisoscillatoryandnon-zero.Thisim- shown to be incorrect (Papaloizou&Pringle 1983), their qualita- pliesthattheinneredgeofthediscmaybemisalignedwithrespect tivefindingsofanalignedinnerdiscsmoothlyconnectedtoamis- totheblackholespin.Assimplepicturesofgasaccretionfavour aligned outer disc — the ‘Bardeen-Petterson effect’ — has been alignmentoftheouteraccretiondiscwiththegalaxydisc,thishas confirmed both from subsequent 1D calculations with corrected beensuggestedasanexplanationfortheobservedrandomorienta- equations (Kumar&Pringle1985; Pringle1992) and inthree di- tionofjetswithrespecttotheirhostgalaxies(Kinneyetal.2000; mensionalsmoothedparticlehydrodynamics(SPH)simulationsby LOP02).However,unlikethediffusiveregimeforwhichtherenow Nelson&Papaloizou(2000). existsafullnon-lineartheorydescribingwarpsofarbitraryampli- Papaloizou&Pringle (1983) showed that there were two tudeand α(Ogilvie1999,2000)only linear theory existsfor the regimes for warp propagation in discs depending on the ra- bending wave regime (Papaloizou&Lin 1995; Lubow&Ogilvie tio of the effective viscosity α to the aspect ratio H/R. For 2000;LOP02;thoughseeOgilvie2006);hencethesestudiesapply α & H/R, warps can be described with a diffusion equation, onlytosmallamplitudewarpsandcouldnotaccountfornon-linear whereas for α . H/R they propagate as bending waves at half effects. the sound speed (Papaloizou&Lin 1995). Previous studies by Thefirst3D numerical simulationsof theBardeen-Petterson effect were performed by Nelson&Papaloizou (2000), using a ⋆ [email protected] post-Newtonian description of the central potential. Their simu- † [email protected] ‡ [email protected] lations in both the wavelike and diffusive regimes largely con- § EinsteinFellow firmedtheBardeen-Pettersoneffect,namelyanalignedinneredge, 2 Nealon, Price&Nixon asmoothtransitiontotheouterdiscandimportantly,noevidenceof materialmaybebrokenifthedisctears,ashasbeenobservedinthe oscillationsinthetilt.ThereasonforthediscrepancywithLOP02 diffusiveregime(Nixon&King2012;Nixonetal.2012b).Finally, is unclear, with Nelson&Papaloizou (2000) claiming that ‘non- itmaynotbepossibleforthedisctofindasteadystateifthedisc linear effects lead to the damping of these short wavelength fea- isrelativelythickandtheviscoustimeisshort. tures’ whilst LOP02 show that these features are not small com- paredtothelocaldiscscaleheight.Furthercomplicatingmatters, 2.1 Istheinnerdiscaligned? recent simulations that included a full general relativistic (GR) treatmentmayhavefoundtheseoscillationsandanon-zerotiltat LOP02considered warps in geometrically thin, almost Keplerian theinneredge(Fragileetal.2007;Zhuravlevetal.2014). discs described by a surface density Σ(R) and angular velocity AfurtherchallengetotheBardeen-Pettersondescriptionhas Ω(R). The scale-height of the disc is given by H(R) cs/Ω, arisen from recent simulations of warped discs in the diffusive where cs(R) is the sound speed in the disc. Their descr≡iption is regime (α & H/R). 3D simulations of warps in isolated discs one dimensional inthesense that thetotal angular momentum in have demonstrated closeagreement withthe non-linear theory of the disc L is a function only of the cylindrical radial coordinate, Ogilvie(1999)(seeLodato&Price2010).However,insimulations R. The disc is then discretised into a series of rings, each de- ofwarpsdrivenbyLense-Thirringprecessionitwasfoundthatin- scribed by the orientation of its tilt and twist angle. The tilt an- steadofmaintainingasmoothtransition(asdescribedbythelinear gle β is measured from the z-axis, and if this angle varies with Bardeen-Pettersoneffect),discsbreakwhentheanglebetweenthe radius the disc is considered to be warped. The twist angle γ is discandtheblackholespinismorethanafewdegrees(Nixonetal. measuredfromanaxisthatisperpendiculartothez-axis,andsim- 2012b).Ringsof gaswerefound tobe tornoff thedisc, precess- ilarly,ifthetwistanglevarieswithradiusthediscistwisted.These ingeffectivelyindependentlybeforebeingaccreted.Theneteffect twoanglescanberelatedtotheunitangularmomentumvectorby was a much higher accretion rate on to the black hole. Thus the l=L/L=(cosγsinβ,sinγsinβ,cosβ)(Pringle1996). Bardeen-Petterson idea of asmooth transition does not appear to We assume an α disc viscosity where ν = αcsH holdinthediffusiveregimeforanarbitrarychoiceofparameters.It (Shakura&Sunyaev 1973). For accretion discs with α . H/R, isnotclearwhetherornotmisaligneddiscsinthewavelikeregime the warp propagates as a dispersive wave (Papaloizou&Pringle willalsobreak. 1983; Papaloizou&Lin 1995). Assuming that the disc is nearly HerewereexaminetheBardeen-Pettersoneffectinblackhole Keplerianandnotself-gravitating,theequationsofmotiondescrib- accretiondiscs.Buildingonthesuccessofearlierstudiesweusethe ingthewavepropagationare(Lubow&Ogilvie2000;Lubowetal. SPHcodePHANTOMtomodelwarpeddiscsin3D(Lodato&Price 2002) 2010; Nixon 2012; Nixon,King&Price 2012a; Nixonetal. 2012b;Facchini,Lodato&Price2013;Nixon,King&Price2013; ∂l 1 ∂G Martinetal. 2014a, Martinetal. 2014b). Here we focus on the ΣR2Ω = +T, (1) ∂t R∂R wavelike propagation regime (α . H/R).Westart by considering the possible reasons for thediscrepancy between LOP02 and Nelson&Papaloizou(2000),aswellastwopossiblecomplications ∂G Ω2−κ2 l G+αΩG=ΣR3Ωc2s ∂l . (2) totheBardeen-Pettersonpictureofanalignedinnerdiscsmoothly ∂t − 2Ω × 4 ∂R (cid:18) (cid:19) joinedtoamisalignedouterdiscinSection2.Wepresentthenu- Here κ is the epicyclic frequency, G represents the internal mericalmethodandtestsofwavelikewarppropagationwithSPH horizontal torqueinthediscandTistheexternal torqueper unit inSection3.InSection4wepresentasuiteof3Dsimulationsde- area. LOP02 chose a complex representation where the warp is signedtoconfirmunderwhatcircumstancestheBardeen-Petterson givenbyW =l +il andtheinternaltorqueasG=G +iG . x y x y descriptionholds.InSection5wediscusstheimplicationsofthese ThisallowsEquations1and2toberewrittenas andinSection6weoutlineourconclusions. ∂W Ω2 Ω2 1 ∂G An obvious caveat is that we do not consider magnetic ΣR2Ω i − z W = , (3) fields, even though it is widely accepted that magnetorota- ∂t − 2Ω R∂R (cid:20) (cid:18) (cid:19) (cid:21) tional instability (MRI) is the controlling mechanism for viscosity in the disc (Balbus&Hawley 1991). However, we can ∂G Ω2 κ2 PR3Ω∂W i − G+αΩG= . (4) still capture the dominant behaviour in the disc as magnetic ∂t − 2Ω 4 ∂R (cid:18) (cid:19) fields have been shown to have little effect on the geometri- These equations describe the propagation of a warp in the linear cal evolution (Sorathia,Krolik&Hawley 2013). We also follow regime,andweresolvednumericallybyLOP02tofindthesteady Nelson&Papaloizou(2000)inusingusingapost-Newtonian ap- state shape of the disc around a Kerr black hole. In thiscase the proachinsteadofgeneralrelativity(GR).Aswewillsee,oneofthe apsidalandnodalprecessionfrequenciesinthedisc(scaledbyΩ) findingsofthispaperisthatthisapproximationmustbeconsidered can be approximated to first order from theKerr metricas (Kato carefullyinordertocapturethecombinationofrelativisticeffects 1990) thatleadtotiltoscillations. ηLOP= κ22−Ω2Ω2 =−23RRs, (5) 2 DOESTHEBARDEEN-PETTERSONEFFECTHOLD INTHEWAVELIKEREGIME? ζLOP= Ω2z2−Ω2Ω2 =−√a2 RRs 3/2, (6) WeconsiderthreepossiblewaysthattheBardeen-Pettersoneffect (cid:18) (cid:19) maybeviolatedinwavelikediscs.Firstly,radialoscillationsinthe where Rs = 2GM/c2 and a is the black hole spin. These fre- tiltofthediscmaypreventthediscfromaligningattheinneredge. quencies are used in the solution by inserting them directly into Secondly, the smooth transition between aligned and misaligned Equations3and4.AnexampleofthesolutionbyLOP02isshown On theBardeen-Pettersoneffect 3 11 1 TiltTilt 00..55 Tilt 0.5 00 0 1100 2200 3300 4400 10 20 30 40 RR//RR R/R iinn in Figure1.ThesteadystatetiltprofilefoundbyLOP02,asintheirFigure3andshowntothesametime(ateightequallyspacedtimes).Theoriginalapsidal andnodalfrequenciesareusedintheleftpanel(Equations5and6)withthesamesign.Intherightpanelwehavereversedthesignofthenodalfrequencyso theprecessionfrequencieshavedifferentsigns,andtherearenooscillationspresentinthissteadystate(similartoLOP02,Figure6).Herethetiltisshownas afractionofthemaximumtilt,theinitialconditionisshownwiththedashedlineandRin=4Rg. intheleftofFigure1,withthesameparametersusedintheirwork. sureandviscosity.Theviscoustorquethatactsbetweensuccessive, Thesteadystatesolutionisformedfromtheinteractionoftheingo- discreteringsinthediscisgivenby(Lynden-Bell&Pringle1974) ingandtheoutgoingbendingwaves,wheretheoutgoingwavesare G=3πνΣ(GMR)1/2. (7) createdbythereflectionoftheingoingwavesattheinnerboundary (LOP02).Heretheoscillatorybehaviourofthesteadystatenearthe Lense-Thirring precession causes the rings that make up the disc inneredgeisclear,asisthenon-zerotiltattheinneredge. toprecess. Per unit area on the disc, thistorque is given by (e.g. It is known that the relative signs of the apsidal and nodal Nixonetal.2012b) frequencies determines whether the solution is oscillatory or not (Ivanov&Illarionov1997).ThefrequenciesusedbyLOP02have T= GMΣR2Ω sinθ Rg 3, (8) thesamesign,leadingtoradialoscillationsinthesteadystatetilt 2a | | R (cid:18) (cid:19) profile.Weconfirmthattheoscillatoryprofileisdependentonlyon thesignsofthefrequenciesbychangingthesignofζLOP(equiva- whereRg = GM/c2,aistheblackholespinandθ istheangle lenttomodellingaretrogradeblackhole,seeLOP02,Figure6)in between the plane of the disc and the direction of the black hole theright hand panel of Figure1. Thetwosolutions evolveinthe spin. If the external torque applied to the disc is greater than the samemannerwiththeexceptionoftheoscillationsneartheinner internal torque maintaining the disc, the rings will precess inde- edge. pendently faster thanthe discisableto communicate the preces- While one is free to set the precession frequencies directly sion (Nixonetal. 2012b). This will result in the disc being sep- when solving Equations 3 and 4, in 3D the nodal precession can arated and breaking, perhaps into differentially precessing rings. be induced directly (e.g. using thepost-Newtonian description of Assumingthatthedischasnoinitialwarpandthatinternal com- Lense-Thirringprecessionfromaspinningblackhole)andtheap- munication isdominated by viscosity, acomparison of theabove sidalprecession(i.e.Einsteinprecession)arisesindirectlyfromthe torques predicts a maximum radius that it is possible for this to central potential. Hence, it is possible for the choice of poten- occur(Nixonetal.2013) tial to preclude oscillations from the steady state solution in 3D simulations. It is then not surprising that simulations that do not 4a H −1 2/3 Rbreak . sinθ Rg. (9) take the apsidal precession into account as above also do not re- 3α| |(cid:18)R(cid:19) ! porttiltoscillations(Sorathiaetal.2013).However,simulationsby Thisapproximaterelationshipplacesanupperboundonthebreak- Nelson&Papaloizou (2000) did make use of a potential that re- ing radius of the disc at a given angle. For the typical parame- sulted in apsidal and nodal precession frequencies with the same signbutdidnotfindoscillations.Hereweusehighresolutionsim- ters used in this paper, we have H/R = 0.05, Rout = 40Rin, α = 0.01 and a = 0.9. Attheouter edge of thedisc, theabove ulationsalongwithapotentialthatleadstoprecessionfrequencies relationthenreducesto ofthesamesigntoinvestigatethisdiscrepancy. Rbreak .45Rin(sinθ)2/3Rg. (10) This predicts that tearing may occur in the disc for inclina- 2.2 Whendoesthediscbreak? tionsofmorethan6◦.Atthisinclinationorgreateronewouldex- ThederivationofEquations1and2assumesthattheinclinationof pectthediscstobreakratherthanalign.However,inthebending thediscislinear.Frompreviousresultsinthenon-linearregimewe waveregimethatweconsiderhere,theinternalcommunicationis wouldanticipatethatthediscmaybreakwhentheexternaltorque dominated bypressure. Inthiscasewecanestimatetheradiusat applied to the disc is stronger than the internal torque. Here the which the disc will break by comparing the sound crossing and internaldisccommunicationisgovernedbyacombinationofpres- theprecessiontimescalesinthedisc.FollowingNixonetal.(2013) 4 Nealon, Price&Nixon and assuming that the disc is inviscid (and hence not taking into as a first order (in v/c) correction in the momentum equa- accountanywavedamping)wefindthat tion. Taking this into account, the momentum equation becomes (Nelson&Papaloizou2000) R 2/3 Rbreak,t .(cid:18)4a|sinθ|H(cid:19) Rg. (11) ddvt =−ρ1∇P +v×h−∇Φ+Svisc. (14) Wethereforeonlyexpectthedisctobreakclosertotheblackhole thanEquation10. HereSviscistheviscousforceperunitmassandv×histhegravo- magneticforceperunitmass,withhgivenby 2S 6(S r)r 2.3 Canthediscaccretemisaligned? h= · , (15) R3 − R5 A further assumption made indeveloping Equations 1and 2was andS=aG2M2k/c3,pointinginsamedirectionastheblackhole that the viscous timescale in the disc is much larger than any spin. other timescale, equivalent to assuming that the disc is replen- Assuming a thin disc with linear disturbances, the quantity ishedfromradiioutsidethecomputationaldomain,orα H/R v hcanbeexpressedas ≪ × (Lubowetal. 2002). This implies that the surface density profile does not change during the evolution of these equations, which v h= 2S v×ˆk + 6S(zr×v). (16) isvalid until the warp reaches the outer boundary. Wecan quan- × (cid:16)R3 (cid:17) R5 tify this approximation during the evolution of the equations by This expression is mainly useful for setting up the initial orbital considering the ratio of the wave and viscous timescales (as in velocitiesforthediscparticles(Section3.2). Lodato&Pringle2006;Facchinietal.2013) twave = 2Rν =2αH. (12) 3.1.1 ImplementationinLeapfrog tν csR2 R A complication to the implementation in the code is that we use For α < H/R 1, the above relationimplies that the viscous ≪ a leapfrog integrator in the ‘Velocity-Verlet’ form, where the po- time is much greater than the sound crossing time that the warp sitionsandvelocitiesoftheparticlesareupdated fromtimetn to communicates, sowecanneglectmassaccretion. Indeed, LOP02 tn+1accordingto neglectedtheevolutionofthesurfacedensityprofilecompletelyin their solution, equivalent to assuming no mass is accreted at all. vn+12 =vn+ 1∆tan, (17) HoweverfortheirdiscH/R = 0.1030,andsoitisnotclearthat 2 thisassumptionholds.Additionally,theviscoustimecanbewritten xn+1=xn+∆tvn+21, (18) as 1 H −2 vn+1=vn+12 + 12∆tan+1. (19) t = . (13) ν αΩ R However,theaccelerationcausedbyLense-Thirringprecessionde- (cid:18) (cid:19) Inthisform,itisclearthatincreasingtheaspectratioofadiscre- pendsonvelocity.Thus(19)becomesimplicit.Thiscanbeeasily sultsinasignificantdecreaseintheviscoustime.Atagivenradius solvedbywritingthecorrectorstepintheform R,whentheviscoustimescaleiscomparableto(orsmallerthan) vn+1 =vn+12 + 1∆tan+1+ 1∆t vn+1 hn+1 , (20) theprecessiontimescaleatthatsameradius,accretiondominates. 2 pos 2 × Inthickdiscs(or tori)itmaythen bepossiblefor thematerialin where an+1 contains the position-d(cid:0)ependent terms(cid:1). This forms a theouterdisctobeaccretedbeforeithasachancetoalign.Thetilt pos setofthreelinearequationsforeachcomponent ofvn+1,thatwe profileinthiscasewillnotreachasteadystatebutinsteadbedeter- solveanalyticallybyinvertingtheresulting3x3matrix. minedbytheinwardfluxofangularmomentum.Thusinrelatively thickdiscsnotiltoscillationswouldbeexpected. 3.1.2 Potentials We make use of two gravitational potentials in this work. The 3 NUMERICALMETHOD firstwaspreviouslyintroducedbyNelson&Papaloizou(2000),re- We use 3D simulations to investigate whether the Bardeen- ferredtoastheEinsteinpotential(seetheirEquation8).Inourno- Petterson effect holds in accretion discs focussing on tilt oscilla- tationitisgivenas tαio.ns,Hla/rgRe.iWncelinuasetioPnHsAaNndTOmMis,aalingneeffidcaicecnrteatinodnlionwthmeelimmoirtywShPerHe ΦE(R)=−GRM 1+ 3RRg . (21) code(Lodato&Price2010;Price&Federrath2010;Price2012). (cid:18) (cid:19) Thispotentialwasintroducedbecauseitpreventsthegravitational Thiscodehasbeenwidelyusedtomodelwarpeddiscsinthediffu- forcetendingtoinfinityastheradiusdecreases.However,italsore- siveregimeandtoinvestigaterelatedproblemswithwarpedwave- sultsinthecorrectapsidalprecessionfrequencyatlargedistances likediscs(seeSection1).Thephysicaldiscviscosityinthediscis fromtheblackholeandhasthesamesignasthenodalfrequency representedinthecodeusingtheartificialviscosity(αAV)method (Nelson&Papaloizou2000).ThisisincontrasttothestandardKe- outlinedinLodato&Price(2010). plerianpotential GM 3.1 ModellingLense-ThirringPrecession Φ(R)= . (22) − R The Lense-Thirring precession exerted by the black hole is Thestandardpotential(22)wasusedinallofthenon-linearsimu- modelled using a post-Newtonian approximation, represented lations,exceptforFigures4–7where(21)wasused. On theBardeen-Pettersoneffect 5 3.1.3 PrecessionFrequencies Simulation θ(◦) a α αAV We calculate the apsidal and nodal precession frequencies in our PS1 30 0.1 0.01 0.395 disc using the standard (Newtonian) definitions for the epicyclic PS2 30 0.1 0.03 1.186 andverticalfrequencies, PS3 30 0.3 0.01 0.395 d PS4 30 0.3 0.03 1.186 κ2 =4Ω2+R (Ω2), (23) PS5 30 0.5 0.01 0.395 dR PS6 30 0.5 0.03 1.186 PS7 30 0.7 0.01 0.395 ∂2Φ(R) PS8 30 0.7 0.03 1.186 Ω2 . (24) z ≡ ∂z2 PS9 30 0.9 0.01 0.395 (cid:12)z=0 PS10 30 0.9 0.03 0.186 (cid:12) Following Nels(cid:12)on&Papaloizou (2000), we compute these using (cid:12) A1 0 0.9 0.01 0.395 aneffectivepotentialthattakesintoaccountthesecondandthird A2 15 0.9 0.01 0.395 termsofEquation14, A3 30 0.9 0.01 0.395 4S√GM 6Sz2√GM A4 45 0.9 0.01 0.395 Φeff(R)=Φ(R)+ 5R5/2 − R9/2 . (25) A5 60 0.9 0.01 0.395 A6 90 0.9 0.01 0.395 ThispotentialaccountsforboththenormalKeplerianpotentialand A7 120 0.9 0.01 0.395 thecorrectionduetothev hterm,whereΦ(R)couldberepre- A8 150 0.9 0.01 0.395 × sentedbyeitherEquation21or22.FirstlyconsideringtheEinstein potential, using Equations 23 and 24 the post-Newtonian apsidal andnodalprecessionfrequencies(scaledbyΩ)aregivenby Table 1. Simulation parameters, including the spin (a) and ηE = 2−Ω12 6GRM4Rg − 3SR√9G/2M , (26) SUhnalkesusrao&theSruwniyseaenvo(te1d9,73th)evaiscccoresittiyon(αdi)scasndalsaortihfiacdialHv/iRscos=ity0(.α05A,Va)n. (cid:20) (cid:21) outerradiusof40Rinandmadeuseof107particles. 2S√GM ζE = −R9/2Ω2 , (27) alignedtotheblackholespin,withtheparticlesarrangedusinga MonteCarloplacementmethod.Eachparticlewasthenrotatedby whereΩ2isgivenby theinclinationangleandassignedavelocityaccordingtothefol- Ω2 = GM 1+ 6Rg 2S√GM. (28) lowingexpressionderivedfromEquation16, E R3 (cid:18) R (cid:19)− R9/2 v = vk4 a2+ R3 a cos(θ), (31) Asthesignsoftheapsidalandnodal precessionfrequencieshere φ c3 "s Rg3 − # are the same throughout the disc, this potential will allow an where v is the Keplerian orbital velocity. The discs were there- k oscillatory profile to develop. Considering now a standard post- foreinitiallytiltedtotheblackholespin,butnotwarped.There- Newtonianpotential,givenbyEquation22,wefindtheapsidaland sultspresentedbelowhavetimeshowninorbitsattheinneredge nodalprecessionfrequenciestobe(againscaledbyΩ) and show the tilt as a function of radius only. This was found 3S fromthesimulationsusingthemethodoutlinedinSection3.2.6of ηPN = 2√GMR3/2 4S, (29) Lodato&Price(2010)whereweusedN = 300sphericalshells. − ForallofthesimulationstheinnerradiuswassetasRin = 4Rg, inordertocomparetotheLOP021Dcode.Atsmallradiiwenote 4S ζPN = − . (30) thattheabsenceofGRlimitsthevalidityofourresults,andindeed 2√GMR3/2 4S − the need to carefully account for relativisticeffects isone of our Herewenoteanimportantdifferencewithrespecttothesolution findings. usedbyLOP02.Asthesignsoftheprecessionfrequencieshereare opposite,thesteadystatetiltprofilewillnothaveoscillationsifthe potentialinEquation22isused. 3.3 Testofwavelikewarppropagation To confirm that we can correctly describe the propagation of warps in the wavelike regime, we use the test described by 3.2 InitialConditionsandScope Fragner&Nelson (2010). They simulated a wavelike accretion Unless otherwise stated, the discs presented all made use of 107 disc with a point mass potential and compared to a 1D calcula- particles, and simulations with106 and 105 were also conducted tion,findingagreementatthe 10%level.Wechoosetocompare ∼ to check convergence. We note that each time the resolution is tothis1DsolutioninsteadofthesolutionfromEquations1and2 changedbetweensimulationswithotherwisethesameparameters, because the linear solution from Fragner&Nelson (2010) allows theartificialviscosityisalteredaccordingtothescalingdescribed thesurfacedensitytoevolve,asoccursinour3Dsimulations. inLodato&Price(2010)sothat thediscshavethesameαinde- Weconduct asimulationusingthesameparametersascited pendentofwhichresolutionisused.Thelocallyisothermalsound intheirFigure1,withH/R = 0.03andα = 0.001andaninitial speedinthediscwassettocs(R) = cs,in(R/Rin)−q andthesur- disturbanceof5◦.Inthiscasewedonotdrivetheevolutionwith facedensityprofileΣ(R) = Σin(R/Rin)−p,wherep = 3/2and Lense-Thirringprecession,sothatwecanisolatethebehaviourdue q = 3/4 to give a constant α viscosity in the disc and uniform towarppropagationonly.OurresultsareshowninFigure2using resolution(Lodato&Pringle2007).Eachdiscwasinitiallysetup 106 and 107 particles. As the disc evolves, the disturbance splits 6 Nealon, Price&Nixon 444 4 3 grees)grees)grees) 333 grees) 3 bits) dedede de or Tilt (Tilt (Tilt ( 222 Tilt ( 2 escale ( 2 m 111 1 Ti g o l 000 0 000 555 111000 0 5 10 1 2 4 6 8 10 R/R 4 4 in Figure3.Precessiontimescalemeasuredfromaninviscidandpressureless es) 3 es) 3 3Ddiscasafunctionofradius(blackcircles), comparedtotheexpected gre gre Lense-Thirringprecession(redline).Rin=4Rg. e e d d Tilt ( 2 Tilt ( 2 4 RESULTS 1 1 4.1 TiltOscillations 0 0 0 5 10 0 5 10 Wefirst investigated whether or not thetilt oscillations predicted Radius Radius byLubowetal.(2002)arephysicalusing3Dsimulations.Thedisc Figure 2. Evolution of bending waves in a disc, not subject to Lense- is initiated with a constant misalignment of 3◦, within the linear Thirringprecession.Thegreen(106)andred(107)linesshowtheresults regime required by Equations 1 and 2. We chose parameters for fromour3Dsimulation.Theblacklineshowstheresultsfromthe1Dcode oursimulationsimilartothatofLubowetal.(2002),withtheex- of Fragner&Nelson (Figure 1 2010), using the same initial parameters. ceptionof thesurfacedensityprofile,theblackhole spinandthe TheagreementbetweenthesetwosolutionsconfirmsthatSPHcanbeused discthickness.Additionally,wemadeuseoftheEinsteinpotential todescribetheevolutionofwarppropagationinthewavelikeregime. outlinedinEquation21,inordertogiveprecessionfrequenciesof thesamesign(Equations26and27).Wesetp=1.5andq=0.75 sothatthediscisuniformlyresolved,asdiscussedinSection3.2. Thedisadvantageisthatthisresultsinloweramplitudeoscillations inthe1Dcode.Tocombatthisweencouragelargeramplitudeos- intotwowavestravellingathalfthesoundspeed(aspredictedby cillationsbyincreasingthespintoa=0.9anddecreasingthedisc Papaloizou&Lin1995);oneinwardandtheotheroutward.Bythe thicknesstoH/R = 0.05. Theevolutionofthetiltasafunction endofthesimulationthesehavefullyseparatedandarebeginning of radius is shown in Figure 4 from a simulation employing 107 tointeractwiththeboundaries. particles. TheSPHsolutionshowsthesamebehaviourasthe1Dsolu- Oneofthemaindifferencesbetweenthissimulationandthose tion,and increasingtheresolutionreduces thediscrepancy. How- conductedbyNelson&Papaloizou(2000)istheangleofinclina- ever, at late times the inner edge of the disc there is increasing tion.Whileour3◦initialtiltwaswellinthelinearregimerequired disagreement, most likely due to differences in the inner bound- bytheanalyticdescriptioninEquations1and2,theminimumincli- ary condition. This test confirms that PHANTOM can be used to nationusedbyNelson&Papaloizou(2000)was10◦.Weexplore describethepropagationofwarpsinthewavelikeregime. theeffectofnon-linearinclinationsinthispotentialbymisaligning thesamediscat15◦.Figure5showsacrosssectionofdensityin theinnerdiscfromthiscalculation.Thetiltprofileafter 600or- ∼ bitsisqualitativelysimilartoFigure4,showingthesameevolution. The quantitative tilt evolution is resolution-dependent, but never- 3.4 TestofLense-Thirringprecession thelessanon-zerotiltandoscillationswerefoundatbothmedium We also perform a simple test of the Lense-Thirring precession. (106particles)andhigh(107particles)resolution. Wesimulateadiscconsistingoftestparticleswithnoviscosityand Figure6showsaresolutionstudyofthetiltandsurfaceden- zerosoundspeed(i.e.α= cs = 0)subjecttoLense-Thirringpre- sityprofilesusing105,106and107particles.Themainartefactof cession.TheinitialvelocitiesaresetusingEquation31andthedisc lowresolution isthataccretion occursfaster andasaresult there isinclinedat30◦.Wethencalculatetheprecessioninthediscasa islessmassattheinneredgeinthelowerresolutioncalculations. functionoftheradiususingtheprocedureoutlinedinAppendixA. Comparisonof thesurfacedensityprofilesindicatesthatΣ(R)is Figure3showsthecomparisonbetweentheprecessionmea- notfullyconvergedneartheinneredge,whichhasadramaticeffect suredfromourdiscandthepredictedprecessioninthedisc,given onthetiltprofiles.However,inalldiscsthereisanon-zerotiltat bytp = R3/(2a).Wefindagreementtowithinmeasurementun- theinneredgeandinthe106 and107 discsradialoscillationsare certaintiesthroughoutthedisc. observed. The wavelength of these oscillations is consistent with On theBardeen-Pettersoneffect 7 3 3 2 2 es) es) e e gr gr e e d d Tilt ( Tilt ( 1 1 0 0 10 20 30 40 2 4 6 8 10 R/R R/R in in Figure4.Timeevolutionoftheanglebetweenthediscplaneandtheblackholespinasafunctionofradiusina3DdiscsubjecttoLense-Thirringprecession, usingsimilarparameterstoLubowetal.(2002)andthesametimesasinFigure1.Theshapeofthisprofiledependssensitivelyonthesurfacedensityatthe inneredgeandhenceonresolution(seeFigure6).Therightpanelshowsazoom-inoftheinnerdisc. Figure5.Cross-sectionviewofthesteady-statetiltoscillationformedinadiscinitiallyinclinedat15◦totheblackholespin(spinaxisisverticalwithrespect tothepage,i.e.alongthezaxis).Thecolourscaleshowsdensity.DiscparametersarethesameasinFigure4,butwithalargerinitialinclination. the criteria given by Lubowetal. (2002). Even using the preces- viscosity. To check this we conduct a low-resolution simulation sion frequencies and surface density profiles in the 1D code that equivalent to simulation E1 of Nelson&Papaloizou (2000) and are appropriate to the 3D simulations (Equations 26 and 27) still βAV =0.Figure7showstheresults(blacksolidline),comparedto doesnotprovideaclosematchwiththe3Dresults.Itisnotclearif anequivalentsimulationwithβAV = 2(reddashedline)andalso thisdiscrepancyisduetonon-linearfluideffects,e.g.asdiscussed comparedtoahigherresolutionsimulations.Athighresolutionwe byNelson&Papaloizou(2000),orsimplyrequireshigherresolu- findtiltoscillationsregardlessofthevalueofβAV(greensolidand tioncalculationstoobtainnumericallyconvergedresults.However bluedottedlines)butwefindthatusingβAV =0canindeederase it is clear that with the appropriate potential and system param- thetiltoscillationsatlowresolution.ThelowerpanelofFigure7 eters, the disc can display radial tilt oscillations as predicted by showsthatthisisnotsimplyduetotheeffectonΣ(R),sincetwo Ivanov&Illarionov(1997)andLOP02. of the calculations show very similar surface density profiles but ratherdifferentevolutionsoftheinnerdisctilt. Despite the resolution-dependence of our results, we were still able to observe tilt oscillations at resolutions used by Nelson&Papaloizou(2000)solongasEinsteinprecessionwasac- 4.2 Whendoesthediscbreak? countedfor.Wefurtherinvestigatedwhether thismight bedueto 4.2.1 Bardeen-PettersonAlignment thedifferencesintheartificialviscosityparametersused,asweset theVonNeumann-RichtmyerviscositycoefficientβAV =2.0(see AsecondpossibleviolationoftheBardeen-Pettersonpicturemay Price2012)forallofoursimulationswhilstNelson&Papaloizou bethatthediscbreaksinsteadofmaintainingasmoothtransition (2000) used βAV = 0. A nonzero βAV viscosity is required to betweenanaligned inner discandamisalignedouter disc. Inor- prevent particle penetration (Monaghan 1989) and the absence der toinvestigatethis,wesimulatearange ofdiscsat30◦ whilst of bulk viscosity is known to be problematic in disc simulations varyingαandaaccordingtothelistPS1-10inTable3.1.3.Here, (Lodato&Price 2010). Thus with βAV = 0, the simulations of forsimplicity,wemakeuseofastandardpotentialgivenbyEqua- Nelson&Papaloizou(2000)mightnothavecapturedthewavein- tion22,andhencedonotexpectanyoscillatorybehaviour. teractions that create the tilt oscillations and the absence of bulk Figure8showsa3Drenderingofdensityinonesuchsimu- 8 Nealon, Price&Nixon 3 10 2 Tilt (degrees) 1 Tilt (degrees) 5 0 0 -8 y Densit ensity -(cid:2) log Surface -1-90 log Surface D -8 ssiimmiillaarr ttoo NNPP22000000 EE11,, bbeettaa==02..00 <h>/H ~ 0.2, our simulation, beta=0.0 <h>/H ~ 0.2, our simulation, beta=2.0 -9 10 20 30 40 R/Rin 0(cid:0)(cid:1) 1R/Rin 1(cid:0)(cid:1) 2 Figure6.Resolutionstudyshowingtheinclination(tiltangle)andsurface Figure 7. The effect of bulk viscosity on the final tilt of the accretion densityasafunctionofradiusatthefinaltimeofthediscshowninFigure4, disc.Thesolid,black line usesthesameparameters as simulation E1of using105(longdashedgreen),106(shortdashedred)and107(solidblack) Nelson&Papaloizou(2000),with5.2×104particlesandhasnobulkvis- particles.Increasingtheresolutionbetterresolvesthesurfacedensityprofile cosity.Bulkviscosityisusedwithboththered(shortdashes,5.2×105 attheinneredge,whichstronglyaffectsthefinaltiltprofilefound. particles)andgreen(longdashes,5.2×106particles)lines.Athigherres- olutionandwithbulkviscositytiltoscillationsareresolved,buttheinner- lation with a = 0.9 and α = 0.03 at three different resolutions mostpartsofthediscremainunconverged.HereRg = 0.04andthedisc hasbeenevolveduntilthewarphasreachedtheouteredge. (PS10).Exceptforthepotentialused,thisdischassimilarparam- eterstosimulationE3ofNelson&Papaloizou(2000)whichmade useof52,000particlesacrossalargerradialextentthanoursim- clearlybetweenthelowviscosity,highspincases.Ata=0.5,the ulations, representingalowerresolution thananyofthoseshown tiltsteepenedandthedisctorebeforetheendofthesimulation.For inFigure 8. At our lowest resolution (left panel of Figure 8), we thediscwitha=0.7,thetiltbegansteepeningneartheendofthe alsoobservetheinnerdiscaligningandsmoothlytransitioningto simulationbutwasnotabletoseparate,whilstata = 0.9thedisc anouter,misaligneddisc(seetheirFigure12).However,athigher has not yet begun steepening. We have confirmed that this is the resolutionsthisbehaviourisnolongerobserved—thediscinstead casebyextendingthehighresolutionsimulationsofthea = 0.9 breaksintotwodistinctsections, withtheinnerdiscalignedwith case,andindeedobservedsteepeningtooccuratlatertimes. theblackholespinandtheouterdiscremainingmisaligned. Aswiththeprevioussimulations,Figure10demonstratesthat Figure9showsthatresolvingdiscbreakingismainlyaques- the simulations are not fully converged, especially when consid- tion of resolving the disc scale height. For the lowest resolution ering the low spin cases (a < 0.5). For these discs, the discrep- simulationstheresolutionlengthisgreaterthanthescaleheightof ancyintheinnertiltisagainduetothemassaccretedattheinner thediscandhencediscbreaking(onalengthscalesmallerthanH) edgeofthedisc.Atlowresolutionstheinnerpartofthediscisac- cannot be resolved and a smooth transition is observed. By con- cretedfaster,resultinginlessmassneartheinneredge.Thesame trast,thetwohigher resolutionsareabletoresolvediscbreaking. Lense-Thirringtorquethenactsonlessmass,andisthusnotable Figure10showsthesameresolutionstudyperformedinFigure8 toalignthedisctothesameextent.Atincreasingspinsthiseffect forallofourdiscsat30◦,wherethegreenlineshowssimulations isobservedless,asthehigherspinprovidesalargertorqueandso thatmadeuseof105particles,redshows106andblackshows107 eventhelowestresolutiondiscsareabletoalign.Inthediscswith particles.Thediscsareshownafter1500orbitsattheinneredge, a < 0.5,increasingtheresolutionleadstoamoredistincttearin allowingthewarptopropagateallthewaytotheouterradius.In- the disc suggesting that our results are consistent. Hence we can creasingthespinoftheblackholeincreasestherateatwhichthe beconfidentthatthesediscsdotear,andpresentanupperlimiton innermostpartofthediscaligns. theradiusatwhichthisoccurs.Asthetearingoccursoutsideofthe Acrossalloftheparameterschosenhere,increasingthereso- radiuswhereoscillationswerefoundinSection4.1,usingsimilar lutionchangesthebehaviourfromthesmoothtiltprofileobserved parameters,thisbehaviourshouldnotbeaffectedbyourchoiceof by Nelson&Papaloizou (2000) to a steepening of the tilt profile potential. andultimatelyadiscthatisbrokenintodistinctsections.Thehigher resolution resultsshow analigned inner edge, amisaligned outer edgeandasharptiltprofileconnectingthese,representingabreak 4.2.2 DiscTearing inthedisc.Thediscssimulatedwithlowerspinsappeartosteepen andtearfasterthanthosewithhigherspinasthebreakoccursfur- Toinvestigatethedependenceofdisctearingonthemisalignment therout(andhencealongerprecessiontime).Thisisobservedmost between the disc and the black hole spin in the wavelike regime On theBardeen-Pettersoneffect 9 !""#"""$%&'()*+,- !$.)++)/0$%&'()*+,- !"$.)++)/0$%&'()*+,- Figure8.Structureofthediscwitha=0.9andα=0.03atincreasingnumericalresolution(lefttoright).AtlowresolutiontheBardeen-PettersonEffectis observed,similartotheresultsofNelson&Papaloizou(2000),butathighresolutionthediscdistinctlytearsintotwoseparatesections.Thecolourindicates density,withwhitebeinghighest. we simulated a suite of discs at different inclinations. We again make use of the traditional post-Newtonian approximation given by Equation 25. We held α = 0.01 and a = 0.9 constant and varied theinclination of thedisc between 0◦ (aligned) and 150◦, notedinTable3.1.3withA1-8.Figure11showsthesesimulations aftermorethan1500orbitsmeasuredattheinneredge.Eachdisc 1 wasinitiallytiltedbutnotwarped.Asthesimulationprogressed,a H warpevolvedinresponse totheLense-Thirringtorqueandinthe h>/ < higherinclinationcasesresultedinthediscbreaking. At 15◦ (top right of Figure 11) the disc was observed to smoothlyaligntothespinoftheblackhole.Attheendofthesimu- lation,thetiltofthediscwasconsistentwiththeBardeen-Petterson effectandissimilartoresultsseeninprevioussimulationsat10◦by Nelson&Papaloizou (2000). Extending the lower resolution ver- sionofthissimulation(with106particles)fortwiceaslongshows 0 10 20 30 thatthedisccontinuestoalignwiththeblackholespin,implying R/R in thatthesteadystateforthisdiscisfullalignment.Incliningthedisc at30◦alsodidnotyetresultindisctearing,howeverthisisbecause Figure9.Resolutionlengthasafractionofthediscscaleheight(hhi/H) forthethreeresolutions(105ingreenlongdashes,106inshortreddashes forthisparticularchoiceofviscosityandspinthissimulationhas and107insolidblack)showninFigure8.Asdiscbreakingoccursonlength notbeenrunlongenough(seeSection4.2.1) scalessmallerthanthescaleheightofthedisc,thelowestresolutionsimu- For discs at higher inclinations (& 45◦; second, third and lationsherecannotresolvebreakingbehaviour(astheresolutionlengthis fourthrowsofFigure11),theinnersectionofthediscwasfoundto greaterthanthescaleheightthroughoutthedisc).Thetwohigherresolution alignwithin50orbitsandasmoothtransitionwasformedbetween simulationsareabletoresolvethisbehaviour. thisandtheouterregionofthedisc.Thistransitionthensteepened untilthediscbrokeintotwosectionsthatwereconnectedbypre- Disc tearing has also been observed in the wavelike cessingringsofmaterial.Multipleringsofmaterialweretornoff disc regime for circumbinary discs inclined at high angles fromtheouter,misaligneddiscandeachwasobservedtoprecess (Facchinietal. 2013). In a simulation of a circumbinary disc in- effectivelyindependently.Towardstheendofthesimulations,upto tworingswereprecessingatthesametime(forexample,120◦disc clinedat 60◦,their discseparatesintotwosections and theinner oneprecessedeffectivelyindependentlyoftheouterdisc.Astheir ofFigure11)andwerepresentforupto 400orbits.Eventually ∼ disc is thicker than ours (H/R = 0.1) and has a higher viscos- eachoftheseringssettledwithandincreasedtheinner,alignedre- gionofthedisc.Thediscinclinedat90◦(righthandpanelinthird ity(α = 0.05), wewouldanticipatethatastrongexternaltorque wouldberequiredtotearthedisc,andweobservetheirdiscdoes rowofFigure11)alsodevelopedprecessingringsofmaterialthat tearatasmallerradiusthananyofours. wereaccreted.However,forthisinclinationnoinneraligneddisc wasobservedandtheringsofmaterialwereaccreteddirectlyonto theblackhole. 4.2.3 Locationoftearingradius Figure12showstheinstantaneousmassaccretionratebythe discsatdifferentangles.Itcanbeseenthatincliningthedisctothe The disc is expected to tear when the Lense-Thirring torque is spinoftheblackholeincreasestherateofaccretionbymorethan largerthan theinternal communication inthedisc. If theinternal anorderofmagnitudewhencomparedtoanaligneddisc,similarto communication in the disc is governed by viscosity, the torques previousfindings(Nixonetal.2012b).Thediscsthatformaninner given in Section 2.2 can be used to estimate the upper breaking aligneddiscandprecessingringshaveevenhigheraccretionrates, radius given in Equation 9. However in our simulations the disc astheinnerdisciscontinuallyfedbytheringsastheyalign.When internaldynamicsaredominatedbypressureratherthanviscosity, takenincontextwiththeresultsinthediffusiveregime(Nixonetal. henceEquation11maybemoreappropriate. 2012b), thisimpliesthat regardless of whether the disc isthinor As the Lense-Thirring torque has a radial dependence, it is thick,massaccretionisfasterwhenthediscisinclined. largestintheinnermostpartsofthediscanditisreasonablethat 10 Nealon, Price& Nixon 333000 222000 Tilt 111000 a=0.1, alpha=0.01 a=0.1, alpha=0.03 333000 222000 111000 a=0.3, alpha=0.01 a=0.3, alpha=0.03 30 20 10 30 20 10 a=0.7, alpha=0.01 a=0.7, alpha=0.03 30 20 10 a=0.9, alpha=0.01 a=0.9, alpha=0.03 10 20 30 40 10 20 30 40 R/R R/R in in Figure10.Disctearingfordifferentspinandviscositycombinations,withthesameinitialtiltof30◦and105particlesshowningreen(widedashes),106 particlesshowninred(dashes)and107particlesshowninblack(solid).AtthelowestresolutionweobservetheBardeen-Pettersoneffectcompletewitha smoothtransitionformostdiscs.Increasingtheresolutionresultsindisctearing(Bardeen-Pettersonalignment)independentofourchoiceofviscosityorspin. AswedonotincludetheaffectofEinsteinprecession,wedonotobserveradialtiltoscillationsinthesediscs. thesediscswillbreakataradiussmallerthanpredictedbyEqua- 4.2.4 Widthoftherings tion9.Figure13showsacomparisonoftheestimatedbreakradius forthesimulationsinclinedat30◦comparedtothepredictionfrom The rings that are torn off during the simulations appear to be muchwiderthanthosefoundinthediffusiveregime(Nixonetal. Equation 9(upper line; assuming that α = 0.02, the averagefor 2012b),someupto∆R/H(R) 25(where∆Rrepresents the our simulations) and from Equation11 (lower line). Wefindthat ∼ ringwidth). It ispossible for rings toform when the disc isable thediscdoesbreakatradiilowerthanourpredictionfromthevis- to break and differential precession is present, such as when the coustorquesalone,andthatthebreakingradiusisintermediatebe- discissubjectedtoLense-Thirringprecession.Wethereforeexpect tweenthepredictionsfromEquations9and11,indicatingthatthe thewidthoftheringtobedeterminedbyarelativecomparisonbe- torquesinourdiscsliebetweenthesetwoextremes.Theincreasing tweenthesound crossingandprecessional timescalesinthedisc. uncertaintiesatlowspincorrespondtothedecreasingconvergence Wecanapproximatethisbyletting∆Rbethedistancethatawave ofoursimulationsduetomassaccretionattheinneredge,seenin cantravelinaprecessiontimesuchthat Figure10. The discrepancy between the predicted and the observed 2 dR tp, (32) btiroenak9i,ntgherabdreiuaskianpgpreaadrisustofoorctchuer6a0t◦adllisicncislifnoautinodnsto. UbesiRng Equa- Z cs ∝ break ∼ acrossthering.If weassume thattheinner edge oftheringisat 41RinwhichisgreaterthanRout.Howeverthisdiscisobservedto break(atR . 18Rin),inlinewiththeresultsofFigure13.Ifwe Rin =R−∆R/2,theouteredgeatRout =R+∆R/2anduse now consider the 15◦ disc, it is predicted to break at R theexpressionforthesoundspeed,weget break ∼ 18Rin but from the simulation we do not observe tearing. This couldoccuriftheactualtearingradiusislessthanRin,consistent R3 a 1 Rin q+1 Rout , (33) withthepreviousresults. ∝ (q+1)" −(cid:18)Rout(cid:19) #cs(Rout)