Circularization of Tidally Disrupted Stars around Spinning Supermassive Black Holes

16.7 MB·English

by  Kimitake Hayasaki

Mon.Not.R.Astron.Soc.000,1–??(2014) Printed9September2016 (MNLATEXstylefilev2.2)
Circularization of Tidally Disrupted Stars around Spinning
Supermassive Black Holes
Kimitake Hayasaki1,2,3(cid:63), Nicholas Stone3,4 and Abraham Loeb3
1Department of Astronomy and Space Science, Chungbuk National University, Cheongju 361-763, Korea
6 2Korea Astronomy and Space Science Institute, Daedeokdaero 776, Yuseong, Daejeon 305-348, Korea
1 3Harvard-Smithsonian Center for Astrophysics, 60 GardenStreet, Cambridge, MA02138, USA
0 4Department of Astronomy, Columbia University, 550 W. 120th Street, New York, NY 10027, USA
2

p
e
S
ABSTRACT
8
Westudythecircularizationoftidallydisruptedstarsonboundorbitsaroundspinning

supermassiveblackholesbyperformingthree-dimensionalsmoothedparticlehydrody-
]
E namicsimulationswithPost-Newtoniancorrections.Oursimulationsrevealthatdebris
H circularizationdependssensitivelyontheefficiencyofradiativecooling.Therearetwo
stagesindebriscircularizationifradiativecoolingisinefficient:first,thestellardebris
.
h streams self-intersect due to relativistic apsidal precession; shocks at the intersection
p points thermalize orbital energy and the debris forms a geometrically thick, ring-like
-
structurearoundtheblackhole.Theringrapidlyspreadsviaviscousdiffusion,leading
o
r to the formation of a geometrically thick accretion disc. In contrast, if radiative cool-
t ing is efficient, the stellar debris circularizes due to self-intersection shocks and forms
s
a a geometrically thin ring-like structure. In this case, the dissipated energy can be
[ emittedduringdebriscircularizationasaprecursortothesubsequenttidaldisruption

flare. The circularization timescale is remarkably long in the radiatively efficient cool-
2
ing case, and is also sensitive to black hole spin. Specifically, Lense-Thirring torques
v
7 cause dynamically important nodal precession, which significantly delays debris cir-
0 cularization. On the other hand, nodal precession is too slow to produce observable
2 signatures in the radiatively inefficient case. Since the stellar debris is optically thick
5 and its photon diffusion time is likely longer than the timescale of shock heating, our
0 inefficient cooling scenario is more generally applicable in eccentric tidal disruption
1. events(TDEs).However,inparabolicTDEsforMBH (cid:38)2×106M(cid:12),thespin-sensitive
0 behavior associated with efficient cooling may be realized.
5
Key words: accretion, accretion discs – black hole physics – galactic: nuclei – hy-
1
: drodynamics
v
i
X
r
a
1 INTRODUCTION surroundingtheseSMBHsisnotaccompaniedbysignificant
emission (Genzel et al. 2003).
Most galaxies are thought to harbor supermassive black
holes (SMBHs) with masses from 105 to 109M at their Tidal disruption events (TDEs) provide a distinctive
(cid:12)
centers (Kormendy & Richstone 1995; Kormendy & Ho opportunitytoprobedormantSMBHsatthecentresofsuch
2013). This is inferred from observing proper motions of inactivegalaxies.MostTDEstakeplacewhenastaratlarge
stars bound to the SMBHs (Scho¨del et al 2002), measur- separation (∼ 1pc) is perturbed onto a parabolic orbit ap-
ing stellar velocity dispersions around SMBHs (Magorrian proaching close enough to the SMBH to be ripped apart
et al. 1998) or detecting radiation emitted from gas accret- by the tidal forces, at the radius rt (cid:39) (MBH/m∗)1/3r∗ =
ing onto the SMBHs (Miyoshi et al. 1995). In the last of 24(MBH/106M(cid:12))−2/3(m∗/M(cid:12))−1/3(r∗/R(cid:12))rS. Here we
these signatures, continuous accretion from a gas reservoir denote SMBH mass with MBH, stellar mass with m∗ and
in active galactic nuclei (AGNs) produces intense radiation radius with r∗, and the Schwarzschild radius with rS =
and powerful outflows and jets. On the other hand, gas ac- 2GMBH/c2, where G and c are Newton’s constant and the
cretionproceedsquiescently,atasignificantlylowerrate,in speedoflight,respectively.Thesubsequentaccretionofstel-
the centers of inactive galaxies. The gas poor environment lar debris falling back to the SMBH produces a character-
istic flare with a luminosity large enough to exceed the Ed-
dingtonluminosityforatimescaleofweekstomonth(Rees
(cid:63) E-mail:[email protected] 1988;Evans&Kochanek1989).RecentobservationsofSwift
(cid:13)c 2014RAS
2 K. Hayasaki, N. Stone and A. Loeb
J164449.3+573451 showed that relativistic jets are associ- clump accretion might be able to produce black holes that
ated with some fraction of TDEs (Bloom et al. 2011; Bur- rotate much more slowly. In contrast to generally aligned
rows et al. 2011; Zauderer et al. 2011; Levan et al. 2011). prograde accretion in AGN, spinning SMBHs undergoing
Candidates for TDEs have also been observed at X-ray, ul- TDEscanrotateeitherretrogradeorprogradewithrespect
traviolet, and optical wavebands (Komossa & Bade 1999; to inflowing gas, with a full range of possible inclinations
Gezari et al. 2012; Arcavi et al. 2014; Holoien et al. 2014; for the transient accretion disk. TDEs therefore act as nat-
Vinko´etal.2015),withinferredeventratesof10−5peryear urallaboratoriesfortestingtheoriesaboutaccretionandjet
per galaxy (Donley et al. 2002; van Velzen & Farrar 2014), launching physics over a full range of prograde and retro-
althoughtheobservedlightcurvesandspectra(Gezarietal. grade inclination angles.
2012) do not always match the simplest theoretical expec- Oursimulationshavefocusedprimarilyontidaldisrup-
tations(Loeb&Ulmer1997;Lodato&Rossi2011;Strubbe tionofstarswitheccentric,ratherthanparabolic(Evans&
& Quataert 2011). Kochanek1989;Ayaletal.2000;Ramirez&Rosswog2009;
Black hole spin is one of two fundamental quantities Guillochonetal.2014),centreofmasstrajectories.Although
characterizing astrophysical black holes, which inevitably thestandardtwo-bodyscatteringmechanismforgenerating
acquirespinangularmomentumasaresultofstandardmass TDEs(Magorrian&Tremaine1999;Merritt&Wang2004)
accretion or chaotic accretion (King & Pringle 2006). Mea- predicts effectively parabolic trajectories (1−e (cid:46) 10−6),
∗
suringblackholespinshasprovenmuchmoredifficultthan other mechanisms can feed stars to SMBHs at lower eccen-
black hole mass estimation, because the dynamical effects tricities. Among these non-standard sources of TDEs, the
ofspinoccurmuchclosertotheeventhorizon.Sinceanac- most promising are binary SMBHs, a recoil accompanying
cretion disk can get closer to the black hole when the black a SMBH merger, and the tidal separation of binary stars.
hole is spinning (Bardeen et al. 1972), a detailed spectral Recent numerical simulations have shown that observable
analysisofdiskX-rayemissioncandeterminetheblackhole properties of these “eccentric” TDEs significantly deviate
spin (Tanaka et al. 1995). Such indirect spin measurements from those of standard TDEs; in particular, the rate of
have been so far made for about 30 SMBHs (Miller 2007), mass return is substantially increased by being cut off at
and recently indicated that the SMBH at the centre of the a finitetime, rather than continuing indefinitely asa power
nearbygalaxyNGC1365hasatleast84%ofthemaximum law decay (Hayasaki et al. 2013). Because of their natu-
theoretically allowed value (Risaliti et al. 2013). rally limited dynamic range, simulations of eccentric tidal
SMBH spins are difficult to measure, but are of signifi- disruption were the first to capture relativistic circulariza-
cantastrophysicalimportance.Spinamplitudeanddirection tionofdebrisaroundSMBHs,whichisextremelycomputa-
significantlyaffecttheefficiencyforconvertingrestmassen- tionally challenging for the canonical parabolic case. These
ergyintoradiation.Whilethemass-to-energyconversionef- simulations found that circularization is driven by general
ficiency reaches ≈42% for an extreme Kerr black hole in a relativistic pericentre shift, which causes shocks to form at
prograde rotation, it is only ≈ 4% for the retrograde case stream self-intersections. The orbital energy dissipated at
(Kato et al. 2008), suggesting a wide range of bolometric these self-intersections subsequently circularizes the debris
disk luminosities depending on the relative inclination of into a more compact accretion disk (Hayasaki et al. 2013).
disk and SMBH spin. A SMBH-disk system can also work ThisisincontrasttopastNewtoniansimulationsofcircular-
as an engine to convert the black hole’s rotational energy ization for parabolic orbits around intermediate mass black
into outflows and jets (Blandford & Znajek 1977; Koide et holes (∼ 103M ), where purely hydrodynamic effects cir-
(cid:12)
al.2002).Theoutflowefficiencydependsonspinmagnitude cularizetidallystretcheddebris(Ramirez&Rosswog2009);
and direction via a large-scale magnetic flux threading the althoughtheseareexpectedtobeineffectiveforSMBH-like
black hole and the disk (Tchekhovskoy et al. 2012). mass ratios (Guillochon et al. 2014).
It is still theoretically uncertain whether such jets will The primary motivation for this work is to investi-
alignwiththeblackholespinaxis,withtheangularmomen- gate the effect of SMBH spin on debris circularization.
tum vector of the accretion disk, or with some other aspect We test the hypothesis that nodal precession due to the
of the magnetic field geometry (Stone & Loeb 2012a). A Lense-Thirring effect can delay the onset of stream self-
misalignedaccretiondiskwillundergodifferentialprecession intersections and strongly retard formation of a luminous
duetoLense-Thirringframedragging.Whileageometrically accretion disk (Cannizzo et al. 1990; Kochanek 1994). If
thindiskwarpsbytheBardeen-Pettersoneffect(Bardeen& proven true, TDE circularization delays could be used as
Petterson1975),ageometricallythickdiskcanprecessasa probes of SMBH spin. Such delays could also decrease the
rigid-body rotator, as has been seen in general relativistic average luminosity of many TDEs; if spin-induced circular-
magneto-hydrodynamic (GRMHD) simulations (Fragile et izationdelayiscommonbutnotuniversal,itwouldproduce
al.2007).VeryrecentGRMHDsimulationshavealsoshown a bimodality in TDE optical emission that could explain
thatahighlymagnetizedgeometricallythickdiskcanwarp a discrepancy between theoretically predicted and observa-
due to electromagnetic torques (McKinney et al. 2013). tionally inferred TDE rates (Stone & Metzger 2014).
The present spin of a SMBH records the history of gas In this paper, we study the circularization of a tidally
accretionandmergerswithotherblackholes,andstatistical disrupted star on an eccentric orbit around a spinning
samplesofSMBHspinsencodevaluableinformationonthe SMBH. In section 2, we describe our numerical approach,
growth history of SMBHs in the universe (Volonteri et al. focusing on the Post-Newtonian corrections we make use
2005; Berti & Volonteri 2008). The SMBHs in most AGN of. In section 3, we examine the results of our numerical
are thought to have accreted sufficient gas in their active simulations in two limiting regimes: one is the radiatively
phase to be rotating near the extreme Kerr limit (Doele- efficient cooling case, in which the photon diffusion time is
man et al. 2012), although events of randomly oriented gas much shorter than the energy dissipation timescale. In the
(cid:13)c 2014RAS,MNRAS000,1–??
Circularization of Tidally Disrupted Stars around Spinning Supermassive Black Holes 3
opposite scenario, where the radiative cooling is inefficient,
the debris circularization proceeds in a qualitatively differ-
ent way. In section 4, we examine the effect of black hole
spin on debris circularization, and the nodal precession of
the newly formed accretion disk by the Lense-Thirring ef-
fect.Finally,section5isdevotedtosummaryanddiscussion
of our scenario.
2 METHODS
We start by describing our numerical methods, with a spe-
cialfocusonhowtotreatrelativisticeffectsinthenumerical
code, and summarize the setup of our physical and numeri-
calmodels.First,wedescribeourproceduresfornumerically
modeling the tidal disruption of stars on bound orbits. The
simulations presented below were performed with a three-
dimensional(3D)SmoothedParticleHydrodynamics(SPH)
code, which is a particle method that divides the fluid into
a set of particles, and is flexible in setting various initial
configurations.Thecodeisbasedonaversionoriginallyde-
veloped by Benz (1990); Benz et al. (1990); and Bate et
al. (1995). Figure 1. Initial configuration of our simulations. The dashed
The SPH equations are composed of a mass conserva- whitecircleanditscentralsmallwhitedotshowthetidaldisrup-
tionequation,amomentumequationwiththeSPHstandard tionradiusrt andtheblackholeattheorigin,respectively.The
artificialviscosity,andanenergyequation.Theirdetailswill runtimetinunitsofP∗ andthenumberofSPHparticlesNSPH
areannotatedatthetop-leftcornerandthebottom-rightcorner,
bedescribedlater,insection2.1.Theseequations,withthe
standard cubic-spline kernel, are integrated using a second- respectively.Boththex-axisandthey-axisarenormalizedbyrt.
The star is initially located at (0.3rt,0) for Models 1-3, and is
orderRunge-Kutta-Fehlbergintegratorwithindividualtime
zoomed into the small square inside the main panel. There, the
stepsforeachparticleandavariablesmoothinglength(Bate
whitesmallarrowindicatesthevelocityvectorofthestar.
et al. 1995), resulting in enormous computational savings
when a large range of dynamical timescales are involved.
The variable smoothing length scheme we used gives the black hole is spinning with spin parameters χ=0.9 for
appropriatespatialresolutioninourcode,butweignorethe Model 5 and χ = −0.9 for Model 6. Both models have an
term proportional to the gradient of the smoothing length. inclination angle i = 0◦ between the spin angular momen-
This term is introduced for calculating the gradient of fluid tum and the axis perpendicular to the orbital plane of the
propertieswhenthesmoothinglengthisvariedinspaceand stellar debris. Models 7 and 8 have the same parameters
time, and is important for ensuring energy conservation if as Model 6 but for i = 90◦ for Model 7 and i = 45◦ for
the gradient of any physical quantities varies over a shorter Model 8, respectively. Model 9 has the same parameters as
scalethanthesmoothinglength(seeBate1995forareview). Model 8 but for χ = 0.9. Model 10 has the same simula-
In our simulations, the specific energy is well conserved for tionparametersasModel1butfor(e,β)=(0.8,5),andhas
all the models. This shows that the term plays no crucial beenperformedtocomparewithModel2aofHayasakietal.
role in our simulations. (2013)(seesection3.2).Eachofthesetensetsofsimulation
We have performed 3D SPH simulations self- parametershasbeenruntwice,withtwodifferentequations
consistently modeling a star from before its entry into the of state (adiabatic and polytropic, discussed in more detail
tidalsphereuptolatetimes,whenthestellardebrishascir- in section 3).
cularized into a disk. We model general relativistic effects,
includingleadingorderSMBHspincorrections,byincorpo-
2.1 Treatment of relativistic effects in SPH
ratingPost-Newtonian(PN)forcesupto2PNintotheSPH
code. We have run ten pairs of simulations of tidal disrup- TheformalismofPost-Newtonian(PN)hydrodynamicswas
tion events with different parameters. The common param- constructed by Blanchet et al. (1990) for the approximate
eters through all of simulations are following: m∗ = 1M(cid:12), treatment of relativistic effects in a non-covariant frame-
r∗ = 1R(cid:12), MBH = 106M(cid:12), γ = 5/3, and a unit of run work. Their formalism is applicable to a moderately rela-
(cid:112)
time P =2πΩ−1 =2π r3/Gm (cid:39)2.8hr. The total num- tivisticself-gravitatingfluid(withgravitationalradiationre-
∗ ∗ ∗ ∗
berofSPHparticlesusedineachsimulationis100K,where action,ifdesired),solongasthePNparameterGM /Rc2
BH
K=1000. We also adopt the standard value of the artificial (foratypicalspatialscaleR)neverexceeds≈10%.Itisnot,
viscosity parameters: α = 1 and β = 2 through all however, simple to implement this formalism into existing
SPH SPH
the simulations. Newtonian SPH codes.
Table 1 summarizes each model. Models 1-4 show the For a typical TDE with an orbital speed v, the PN pa-
eccentric TDEs around non-spinning SMBHs with (e,β) = rameter is estimated to be O(v2/c2) = 10−2 at the tidal
(0.9,1), (0.8,1), (0.7,1), and (0.7,2). Models 5 and 6 have disruptionradius.Themagnitudeoftheself-gravitatingpo-
the same simulation parameters as Model 4, except that tential and thermal energy of the star can be similarly
(cid:13)c 2014RAS,MNRAS000,1–??
4 K. Hayasaki, N. Stone and A. Loeb
tionradiusoftheblackholeastheradiusofthemarginally
stable orbit for the non-spinning black hole (Bardeen et al.
1972): r (cid:39)0.12r , in all the models.
ms t
The initial position and velocity of the star is given
by that of a test particle orbiting around the black hole.
In the test-particle limit, the specific energy and angular
momentum with PN corrections are given from equations
(A7) and (A9) by
E (r =0,v =0)
(cid:15) = i BH BH , (1)
tp m
i
J (r =0,v =0)
j = i BH BH , (2)
tp m
i
where the index i refers to a given SPH particle, and r
BH
and v are the position and velocity vector of the black
BH
holeparticle,respectively.Thisenergyandangularmomen-
tumshouldapproximatelyequaltheirrespectiveNewtonian
analoguesatadistancefarawayfromtheblackhole.Given
an initial position and desired pericenter distance, we nu-
merically solve for an initial velocity vector using the PN
constants of motion. The initial velocity and position vec-
tor in our simulation models are summarized in Table 2.
Figure 2. Orbits of two test particles and one SPH particle The number of initial SPH particles NSPH are 100K for all
(from our full disruption simulation with the adiabatic equation the models. Figure 1 shows an initial configuration of our
of state) around a spinning SMBH. The parameters of Model 5 simulations for Models 1-3.
areadoptedforthethreeparticles.Eachaxisisnormalizedbythe Figure2showsorbitsoftwotestparticlesandoneSPH
tidal disruption radius. The initial positions of the particles are
particle (from our full disruption simulation with the adia-
locatedatthecenterofthesmallyellowcircle.Thecentralwhite
baticequationofstate)aroundaspinningSMBH.Thesolid
point shows the black hole with (χ,i) = (0.9,0◦). The solid red
red and dashed white lines show the motion of a SPH par-
anddashedwhitelinesshowthemotionofaSPHparticleanda
ticle and a test particle in the gravitational potential with
test particle in the gravitational potential with Post-Newtonian
corrections (up to 2PN). The dotted line denotes an orbit of a Post-Newtoniancorrections(upto2PN).Thedottedlinede-
testparticlemovingintheKerrmetric. notes an orbit of a test particle moving in the Kerr metric.
The orbit of the test particle in the 2PN potential deviates
slightlyfromthatoftheKerrmetric.Itisinitiallyidentical
parameterized to be O((v2/c2)(m∗/MBH)2/3) = 10−6 and withtheorbitoftheSPHparticleforthefirstfiveorbits,but
O(c2s/c2) = 10−5 where cs is the sound speed (for a stellar the two diverge afterwards because of hydrodynamic forces
temperature ∼ 107K), respectively. These order of magni- on the SPH particle.
tude estimates show that even the lowest PN order terms
for stellar self-gravitation and thermal energy can be self-
consistently neglected, even if up to 2PN precision in the 2.3 Errors of energy and angular momentum
blackhole’sgravityisdesired.Becauseweonlyneedtomod- conservation
ify the SMBH potential, and can continue to treat hydro-
Inordertocheckconvergenceofenergyandangularmomen-
dynamics and gas self-gravity in a Newtonian fashion, it
tumconservationwithPNcorrections,wenumericallysolve
becomesmuchsimplertoimplementthePNformalisminto
thetwo-bodyproblemwithPNcorrectionsinatestparticle
our SPH code.
limit by using a fourth-order Runge-Kutta method. While
In order to treat approximately relativistic effects such
theenergyandangularmomentumarewellconservedinthe
as pericentre shift and spin-induced precession, we have
caseofacircularbinary,theyoscillatewithtimeinaneccen-
incorporated acceleration terms corrected by the Post-
tric binary case (but remain conserved in a time-averaged
Newtonian approximations into the momentum equation of
sense). The oscillation amplitude grows with increasing or-
SPH particles. The detailed formulae can be seen in Ap-
bital eccentricity of the test particle and increasing ratio of
pendix A.
thetidaldisruptionradiustopericentredistanceofthetest
particle.
In Table 2, we compare the error levels of energy and
2.2 Initial conditions
angular momentum conservation between the test-particle
We have performed two-stage simulations: a star is first simulations and the SPH simulations. Each error level is
modeled as a polytropic gas sphere in hydrostatic equilib- measured by
rium.Thetidaldisruptionprocessisthensimulatedbyset-
(cid:15)¯−(cid:15) ¯j−j
ting the star in motion through the gravitational field of a δ(cid:15)= ∗, δj = ∗, (3)
(cid:15) j
black hole. In our simulations, the black hole is represented ∗ ∗
by a sink particle with the appropriate gravitational mass where (cid:15)¯and ¯j are the time-averaged and number-averaged
M . All gas particles that fall within a specified accretion PN values of specific energy and angular momentum dur-
BH
radius are accreted by the sink particle. We set the accre- ing the first ten orbits, while (cid:15) and j are the Newtonian
∗ ∗
(cid:13)c 2014RAS,MNRAS000,1–??
Circularization of Tidally Disrupted Stars around Spinning Supermassive Black Holes 5
Table 1. Tabulatedsimulationparameters.Thefirstcolumnshowseachsimulatedmodelnumber.Thesecondtoseventhcolumnsare
thepenetrationfactorβ=rp/rT,theinitialorbitaleccentricitye∗,theinitialsemi-majoraxisa∗,theradialdistancebetweentheblack
hole andthe initial positionof the star,the specific orbital binding energy ofthe star (cid:15)∗ =−(1/2)β(1−e∗)(cid:15)t where (cid:15)t =GMBH/rt (cid:39)
1.9×1019[erg/g], and the black hole spin parameter χ (with values between 0 and 1), respectively. The eighth column indicates the
anglebetweentheblackholespinaxisandtheaxisperpendiculartotheorbitalplaneofthestellardebris.Theninthandtenthcolumns
aretheperiodsofthemosttightlyandlooselyboundorbits,respectively(seeequations(15)and(16)).Theeleventhcolumnshowseach
(cid:112)
termination time normalized by P∗ =2π rt3/GMBH (cid:39)2.8hr. The last two columns describe the number of SPH particles at the end
ofsimulationsforradiativelyefficient(N )andinefficientcooling(N )cases.
eff ineff
Model β e∗ a∗[rt] r0[rt] (cid:15)∗[(cid:15)t] χ i tmtb[P∗] tmlb[P∗] tend[P∗] Neff Nineff
1 1 0.9 10 3 −0.05 0 0◦ 11 44 100 99142 96854(t =80)
end
2 1 0.8 5 3 −0.1 0 0◦ 4 13 100 99540 91140
3 1 0.7 10/3 3 −0.15 0 0◦ 2.2 6.7 100 99980 85675
4 2 0.7 5/3 2.5 −0.3 0 0◦ 0.76 2.3 40 99830 74520
5 2 0.7 5/3 2.5 −0.3 0.9 0◦ 0.76 2.3 40 99824 81610
6 2 0.7 5/3 2.5 −0.3 −0.9 0◦ 0.76 2.3 40 99687 71148
7 2 0.7 5/3 2.5 −0.3 −0.9 90◦ 0.76 2.3 40 99805 72233
8 2 0.7 5/3 2.5 −0.3 −0.9 45◦ 0.76 2.3 40 99869 69329
9 2 0.7 5/3 2.5 −0.3 0.9 45◦ 0.76 2.3 40 99907 82574
10 5 0.8 1 1.8 −0.5 0 0 0.35 1.03 10 99632 −
Table2. Initialconditionsanderrorsforoursimulations.Thefirstcolumnshowseachsimulatedmodel.Thesecondandthirdcolumns
(cid:112)
denotetheinitialpositionandvelocityvectorforeachmodel.Thenormalizationofthevelocityisgivenbyvt= GMBH/rt.Thefourth
andfifthcolumnshowtheenergyconservationerrorandangularmomentumconservationerrorofatestparticle,respectively.Thelast
twocolumnsdescribetheenergyconservationerrorandangularmomentumconservationerrorofaSPHparticle,respectively.
Model r0[rt] v0[vt] δ(cid:15)tp[%] δjtp[%] δ(cid:15)SPH[%] δjSPH[%]
1 (0.0,−3.0,0.0) (0.447,0.589,0.0) 1.0 0.08 0.054 0.041
2 (0.0,−3.0,0.0) (0.436,0.512,0.0) 1.0 0.07 0.041 0.026
3 (0.0,−3.0,0.0) (0.506,0.474,0.0) 0.27 0.02 0.16 0.05
4 (0.0,−2.5,0.0) (0.359,0.25,0.0) 4.0 0.45 0.25 0.02
5 (0.0,−2.5,0.0) (0.359,0.251,0.0) 3.5 0.35 0.12 0.31
6 (0.0,−2.5,0.0) (0.359,0.248,0.0) 4.8 2.8 0.38 0.23
7 (0.0,−2.5,0.0) (0.359,0.249,0.0) 4.2 0.45 0.25 0.4
8 (0.0,−2.5,0.0) (0.359,0.251,0.0) 4.5 2.0 0.19 0.19
9 (0.0,−2.5,0.0) (0.359,0.249,0.0) 3.5 0.18 0.31 0.6
10 (0.556,−1.71,0) (0.377,0.1225,0) 10.6 1.74 5.1 22.6
specific energy and angular momentum of the initially ap- where a , e , and β are the semi-major axis and orbital
∗ ∗
proachingstar,respectively.ExceptforModel10,theenergy eccentricityoftheinitiallyapproachingstar,andtheratioof
and angular momentum is well conserved at an error level tidaldisruptionradiusr topericenterdistancer =a (1−
t p ∗
of 2%. e ), respectively. The specific binding energy of the stellar
∗
debris measured at the circularization radius can then be
written as
1 β
3 STELLAR DEBRIS CIRCULARIZATION (cid:15) =− (cid:15) , (5)
c 21+e t
∗
Recent numerical simulations have shown that the peri-
where (cid:15) = GM /r is a characteristic specific energy of
centre shift of the stellar debris plays an essential role in t BH t
the tidal disruption radius. On the other hand, the specific
quickly forming an accretion disk around a non-spinning
orbital energy of the initially approaching star is:
SMBH, because it leads to debris orbit self-intersections,
whichdissipateenergyinshocksandcauserapidcirculariza- 1
(cid:15) =− β(1−e )(cid:15) . (6)
tion (Hayasaki et al. 2013). This work also showed that the ∗ 2 ∗ t
angularmomentumofthestellardebrisisconservedduring Notethatthe1PN,1.5PN,and2PNordertermsarepropor-
circularization.Thisangularmomentumconservationallow tionalto(GM/r )/c2 ∼2.1%,(GM/r )3/2/c3 ∼0.31%,and
t t
ustoestimatethecircularizationradiusofthestellardebris, (GM/r )2/c4 ∼ 0.045%, respectively, at the tidal disrup-
t
which is given by tion radius. This shows that equations (4)-(6) are typically
1+e corrected at the ∼2.5% level by our PN approaches.
r =a (1−e2)= ∗r , (4)
c ∗ ∗ β t The difference between m∗(cid:15)∗ and m∗(cid:15)c gives the maxi-
(cid:13)c 2014RAS,MNRAS000,1–??
6 K. Hayasaki, N. Stone and A. Loeb
Figure3. EvolutionofthespecificangularmomentumandenergyforModels1-9intheradiativelyefficientcoolingcases.Theseconstants
are averaged out per SPH particle. In panels (a1)-(a3), the black, blue, and red solid lines denote the specific angular momentum
√
normalized by jt = GMBHr(cid:112)t. The corresponding dashed lines show the specific angular momentum of a test particle moving in
a Newtonian potential, j∗ = a∗(1−e2∗). In panels (b1)-(b3), the black, blue, and red solid lines represent the specific binding
energy normalized by (cid:15)t = GMBH/rt. The corresponding dashed lines show the Newtonian specific binding energy of a test particle,
(cid:15)∗=−(1/2)β(1−e∗)(cid:15)t.ThedottedlineshowstheNewtonianspecificbindingenergymeasuredatthetidaldisruptionradius.Therun
(cid:112)
timetisinunitsofP∗=2π rt3/GM (cid:39)2.8hr.
mumamountofbindingenergypotentiallydissipatedduring (cid:18) M (cid:19)−2/3(cid:16)r (cid:17)−1(cid:16)∆r(cid:17)−1
× BH , (9)
debris circularization: 106M r r
(cid:12) t t
m βe2
δ(cid:15) = m |(cid:15) −(cid:15) |= ∗ ∗ (cid:15) (cid:39)1.9×1052[erg] whereκ =0.4[cm2g−1]istheopacityforelectronscatter-
max ∗ ∗ c 2 (1+e ) t es
∗ ing, and
βe2 (cid:18)m (cid:19)4/3(cid:18) r (cid:19)−1(cid:18) M (cid:19)2/3
× ∗ ∗ bh .(7) (cid:18) (cid:19)−2(cid:18) (cid:19)5/3
(1+e) M(cid:12) R(cid:12) 106M(cid:12) Σ ≡ m∗ (cid:39) 6.5×106[gcm−2] r∗ m∗
0 2πr∆r R M
(cid:12) (cid:12)
Itiscrucialtoconsiderwherethedissipatedenergygoes
(cid:18) M (cid:19)−2/3(cid:16)r (cid:17)−1(cid:16)∆r(cid:17)−1
duringdebriscircularization.Thephotondiffusiontimescale × BH (10)
of the stellar debris is given by (Mihalas & Mihalas 1984) 106M(cid:12) rt rt
H (cid:16) κ (cid:17)(cid:16) Σ (cid:17)(cid:16) H (cid:17)(cid:18) r (cid:19)−2 is the fiducial surface density, where r and ∆r are the ra-
t = τ (cid:39)6.1×108[s] ∗ dialsize andwidthof thedebrisring,respectively.We note
diff c κ Σ ∆r R
es 0 (cid:12) that the stellar debris is clearly optically thick. If the pho-
(cid:18)m (cid:19)5/3(cid:18) M (cid:19)−2/3(cid:16)r (cid:17)−1 tondiffusiontimescaleislongerthantheenergydissipation
× ∗ BH , (8)
M 106M r timescale(i.e.thedebriscircularizationtimescale),thenra-
(cid:12) (cid:12) t
diative cooling is inefficient and dynamically unimportant.
where H, κ, Σ, and τ are the scale height, opacity, surface Otherwise,theradiativecoolingcanbeefficient;thedynam-
density,andopticaldepthofthestellardebris,respectively. icaleffectofthiswillbetoreducethethicknessofthedebris
The optical depth is approximately estimated to be streams.
Itisclearfromequation(8)thatintheeccentricTDEs
τ = κρH ∼κΣ(cid:39)2.6×106(cid:16) κ (cid:17)(cid:16) Σ (cid:17)(cid:18) r∗ (cid:19)−2(cid:18)m∗ (cid:19)5/3we simulate, tdiff will always be very long compared to or-
κ Σ R M bital timescales, and cooling will generally be radiatively
es 0 (cid:12) (cid:12)
(cid:13)c 2014RAS,MNRAS000,1–??
Circularization of Tidally Disrupted Stars around Spinning Supermassive Black Holes 7
Figure4. EvolutionofthespecificangularmomentumandenergyforModels1-9intheradiativelyinefficientcoolingcases.Thefigure
formatsarethesameasFigure3.
inefficient. However, this is not necessarily true for the parabolic TDEs, we obtain the following condition for ra-
parabolic (e ≈ 1) TDEs which dominate the event rate. diatively efficient cooling:
∗
If we approximate the frozen-in specific energy spread of
e∗ ≈1 tidal debris as (Stone et al. 2013) (cid:16) κ (cid:17)6/7(cid:18) r (cid:19)−15/7(cid:18)m (cid:19)16/7
M (cid:38)1.6×106M ∗ ∗ .(14)
∆(cid:15) = GMBHr∗, (11) BH (cid:12) κes R(cid:12) M(cid:12)
t r r
t t
Under our simplifying assumptions, radiative cooling will
then we can calculate the fallback time with Kepler’s third
be, in general, likely to play some dynamical role for de-
law for the most tightly bound debris to be
bris streams of parabolic TDEs. The radiative efficiency of
1 (cid:16)M (cid:17)1/2 (cid:18) M (cid:19)1/2 these streams could be reduced if they cool to the point
t = √ BH P =3.5×106[s] BH
fb 2 2 m∗ ∗ 106M(cid:12) where bound-free absorption dominates electron scattering
as a source of opacity (Kochanek 1994), but it could also
(cid:18) (cid:19)−1(cid:18) (cid:19)3/2
× m∗ r∗ . (12) be increased if magnetically driven turbulent advection of
M R photons enhances cooling rates (Jiang et al. 2014).
(cid:12) (cid:12)
Giventhemanyuncertaintiesinthisdiscussion,andthe
We can estimate the regime where cooling is important by
possibleapplicabilityofbothregimestoparabolicTDEs,we
requiring t < t ; this is conservative because circular-
diff fb present two extreme cases: one involves a set of radiatively
ization likely takes several fallback times to complete. If we
efficientcoolingsimulations,wheretheentropyremainscon-
assume roughly cylindrical debris streams, with H = ∆r,
stant through the simulation, in the polytropic equation of
thenfactorsof∆r cancelandweareleftwithasimplecon-
state with γ =5/3. The other involves the radiatively inef-
ditiononthemaximumextentofthedebrisstream,r,which
ficient cooling simulations, where the entropy is locally in-
we hereafter identify as debris apocenter r . Specifically,
a creasedwithadiabaticequationofstatebutthetotalenergy
r (cid:16) κ (cid:17)(cid:18) M (cid:19)−5/6(cid:18) r (cid:19)−5/2(cid:18)m (cid:19)7/3 is conserved. The detailed parameters of ten (9+1) simula-
a (cid:38)1.7×102 BH ∗ ∗ .(13)tion models are shown in Table 1. Note that Model 10 is
r κ 106M R M
t es (cid:12) (cid:12) (cid:12) done for the purpose of comparison to the fiducial simu-
This illustrates why eccentric TDEs should generally be lation model of Hayasaki et al. (2013), where the pseudo-
in the radiatively inefficient limit, but if we substitute in Newtonian potential was adopted. Further details are de-
the apocenter of the most tightly bound debris streams for scribed in section 3.2.
(cid:13)c 2014RAS,MNRAS000,1–??
8 K. Hayasaki, N. Stone and A. Loeb
3.1 Radiatively efficient cooling cases
First,wedescribetheresultsofourradiativelyefficientcool-
ing simulations, which serve as one extreme of possible ra-
diative cooling. Figure 3 shows the evolution of the spe-
cific angular momentum and specific energy in Models 1-
9, which are averaged per SPH particle. Panels (b1)-(b3)
show the evolution of the specific binding energy in Mod-
els 1-9. The thermal energy is estimated to be the order of
10−5(cid:15) forallthemodels,andthereforeitisnegligiblysmall
t
comparedwiththemagnitudeoftheorbitalbindingenergy.
Panels (a1)-(a3) show the evolution of specific angular mo-
mentumperSPHparticle.Thesmall-amplitudeoscillations
seenthereareduetothePNtermsaddedintheSPHequa-
tionofmotion.Equivalenterrorsarealsoseeninenergyand
angular momentum conservation for a test particle on an
eccentric orbit moving under a gravitational potential with
PNcorrections.WecomparedtheseerrorsbetweentheSPH
simulations and the test particle integrations. The detailed
resultscanbeseenincolumns4-7ofTable2inSection2.3.
NotethattheangularmomentumoftheSPHparticlesinall
themodelsisconservedatalessthan2%errorlevelthrough
the simulations.
Figure5. Asequenceofsnapshotsofthetidaldisruptionprocess
Figures 5-7 show a sequence of snapshots of the sur-
inModel1(a∗ =10rt,e∗ =0.9,β =1,χ=0.0,andi=0◦)in
face density of stellar debris, which is projected on the x-y
theradiativelyefficientcase.Panel(a)to(d)areshowninchrono-
plane in a logarithmic scale, covering two orders of mag- logical order. All thepanels show asurface density projectedon
nitude, for Models 1-3. Each figure progresses from panel thex-yplanefor0(cid:54)t(cid:54)100.Thecolorbarshowsthemagnitude
(a) to panel (d) in chronological order. The central small ofthedensityinalogarithmicscale,whereΣ0=6.5×106gcm−2
point, dashed circle and white small arrows show the black is the fiducial surface density (see equation 10). The black hole
hole,tidaldisruptionradius,andvelocityfieldofthestellar issetattheorigin.TheruntimetisinunitsofP∗ andisanno-
debris, respectively. The run time is noted at the top-right tated at the top-right corner. The number of SPH particles are
indicatedatthebottom-rightcorner.Thewhitesmallarrowsand
corner in units of P , while the number of SPH particles is
∗ the dashed circle indicate the velocity field of the stellar debris
indicated at the bottom-right corner.
andthetidaldisruptionradius,respectively.
The stellar debris moves around the black hole for sev-
eralorbits.Overtime,thedebrisstretchesduetothespread
in its constituent orbital energies, and the debris head in- difference between the orbital periods of most tightly and
teracts with the tail near apocentre, leading to significant loosely bound gas, which shrinks as the orbital eccentricity
energydissipationinshocks.Thebindingenergyofthestel- decreases.Theseperiodsaregivenby(Hayasakietal.2013)
lar debris is substantially reduced by a sequence of orbit
crossings, causing the debris to circularize. From Figures 5 1 (cid:20) 1 (cid:21)3/2
t = √ P , (15)
and6,thestellardebrisclearlycircularizesinModels1and mtb 2 2 β(1−e∗) ∗
2.However,circularizationhasproceededmuchlessrapidly
in Model 3, as we see from Figure 7.
(cid:20) (cid:21)−3/2
1 β(1−e )
Panel (b1) shows the evolution of specific binding en- t = √ ∗ −q1/3 P , (16)
ergyinModels1-3.Sincethespecificbindingenergyhasnot mlb 2 2 2 ∗ ∗
reducedfrom(cid:15) to(cid:15) evenattheendoftheseruns,thecircu-
∗ c whereq =m /M .Panel(b1)showsthattheenergydis-
∗ ∗ BH
larizationprocesshasnotyetcompleted.Adoptingasimple
sipation rate increases with orbital eccentricity for the case
extrapolation from (cid:15) to (cid:15) , the circularization timescales
∗ c of β = 1. This implies that stellar debris should efficiently
can be estimated to be ∼ 120P in Model 1, ∼ 180P in
∗ ∗ circularize in most standard, parabolic TDEs (β = 1 and
Model 2, and ∼ 2500P in Model 3, respectively. These
∗ e =1) around non-spinning SMBHs.
∗
extrapolated timescales indicate strongly varying per-orbit
efficiencies of shock dissipation: if this efficiency were con-
stant, then t ∝ a3/2. In fact, the circularization timescale
c 3.2 Radiatively inefficient cooling cases
behaves in an inverse manner. The prominently long cir-
cularization timescale of Model 3 shows this counterintu- Next, we describe the results of our radiatively inefficient
itive behavior. This declining dissipation efficiency at fixed cooling simulations for Models 1-3. Figures 8-10 show a se-
β and decreasing e is likely because the relative velocity quenceofsnapshotsofdebrissurfacedensityforModels1-3,
∗
between the debris head and the debris tail at their self- withthesamefigureformatsasinFigure.5,butforthera-
intersection decreases as we go from Model 1 to Model 3. diatively inefficient cooling simulations.
This is a general feature of eccentric TDEs at fixed β: low InFigure8,thestellardebrismovesawayfromtheblack
eccentricity produces self-intersections closer to apocentre, holefollowingtidaldisruption,asshowninpanel(a),andis
with lower relative velocities. This is also confirmed by the then stretched during pericentre return in panel (b). Near
(cid:13)c 2014RAS,MNRAS000,1–??
Circularization of Tidally Disrupted Stars around Spinning Supermassive Black Holes 9
Figure6. Asequenceofsnapshotsofthetidaldisruptionprocess Figure8. Asequenceofsnapshotsofthetidaldisruptionprocess
in Model 2 (a∗ =5rt, e∗ =0.8, β =1, χ=0.0, and i=0◦) in inModel1(a∗ =10rt,e∗ =0.9,β =1,χ=0.0,andi=0◦)in
theradiativelyefficientcase.Thefigureformatsarethesameas the radiatively inefficient case. The figure formats are the same
Figure5. asFigure5,butfor0(cid:54)t(cid:54)80.
Figure7. Asequenceofsnapshotsofthetidaldisruptionprocess
Figure9. Asequenceofsnapshotsofthetidaldisruptionprocess
itnheMroaddeialt3iv(ealy∗=effi1c0iernt/t3c,aes∗e.=T0h.e7,fiβgu=re1,foχrm=a0ts.0a,raentdhie=sa0m◦e)ains in Model 2 (a∗ =5rt, e∗ =0.8, β =1, χ=0.0, and i=0◦) in
the radiatively inefficient limit. The figure formats are the same
Figure5.
asFigure5,butfor0(cid:54)t(cid:54)80.
apocentre, the leading “head” of the debris significantly in-
tersectswiththetrailing“tail,”ascanbeseeninpanel(c). 2-3.Whilenoaccretiondiskformsbytheendofsimulation
Afterseveraltensoforbits,thedebrisexpandssignificantly, in Model 1, an accretion disk clearly forms and viscously
as the thermal energy increases from shock heating. This evolvesinModels2and3.Thisviscousevolutioncanbeun-
can be seen in panel (d). derstood as follows: the viscous timescale for an α-viscosity
ThereisacleardifferencebetweenModel1andModels disk is given by
(cid:13)c 2014RAS,MNRAS000,1–??
10 K. Hayasaki, N. Stone and A. Loeb
and H/r ∼ 1 for Model 1. Substituting into equation (8),
t
thephotondiffusiontimescaleis∼107s,clearlylongerthan
theshookheatingtimescale.Itisthereforeclearthateccen-
tric TDEs operate in the radiatively inefficient regime, al-
thoughaswehavearguedabove,theparaboliccaseismore
ambiguous.
3.3 Comparison to pseudo-Newtonian potential
simulation
In this section, we compare our SPH simulation with
2PN corrections to our SPH simulation with a pseudo-
Newtonianpotential.ThePNsimulationofModel10,whose
parameters can be seen in Table 1, was performed for the
purpose of comparing with the pseudo-Newtonian simula-
tionmodel.InModel10,weusedthesameinitialcondition
as that of Model 2a in our previous paper (Hayasaki et al.
2013). The initial position and velocity are set by
r = (r cosφ,r sinφ,0),
0 0 0
v = (r˙(r )cosφ −r φ˙(φ )sinφ ,r˙(r )sinφ
0 0 0 0 0 0 0 0
+ r φ˙(φ )cosφ ,0),
0 0 0
Figure 10. Asequenceofsnapshotsofthetidaldisruptionpro-
wherer =a (1−e )andφ =−0.2πareadopted.Here,the
ice=ss0i◦n)Minotdheelr3ad(iaa∗tiv=ely10in/3effirtc,ieen∗t=lim0i.t7.,Tβhe=fig1u,rχef=orm0.a0t,saanrde radial v0elocit∗y r˙ an∗d angul0ar velocity r˙ and φ˙ are given by
thesameasFigure5,butfor0(cid:54)t(cid:54)80. energyconservationandangularmomentumconservationas
(cid:114)
l2
r˙ = 2((cid:15) −U(r))− pseudo,
r2 5 (cid:16)0.1(cid:17)(cid:16)r (cid:17)3/2(cid:16)H(cid:17)−3/2 pseudo r2
t ∼ ∼ P . (17)
vis ν π αSS rt r ∗ φ˙ = lpseudo,
r2
The viscous timescale can be comparable to the debris or-
bitalperiodbecauseoftheenhancedpressureandresultant where (cid:15)pseudo and lpseudo are the specific energy and the
geometrically thick structure (H/r ∼ 1) in the radiatively specific angular momentum for bound orbits, the pseudo-
inefficientregime,meaningthatrapidviscousredistribution Newtonian potential (Wegg 2012), respectively. They are
of angular momentum and energy is possible. However, our written by
simulations become significantly less reliable at late times, (cid:20) (cid:21)
GM 1−c r
after formation of the accretion disk, due to the lack of ra- U(r) = − BH c + 1 +c S ,
r 1 1−c (r /2r) 32r
diation pressure and magneto-hydrodynamics in our code. 2 S
Wepresenttheselate-timeresultsforcompleteness,butem- (cid:15) = (rp/ra)2U(rp)−U(ra), (18)
phasize that our simulations are primarily designed only to psuedo (r /r )2−1
p a
simulate the circularization process itself. (cid:112) (cid:112)
l = 2r2((cid:15)−U(r ))= 2r2((cid:15)−U(r )),
Figure4showsevolutionofthespecificangularmomen- pseudo p p a a
tum and specific energy in Models 1-9, which are averaged wherer =a (1−e )andr =a (1+e )arethepericenter
p ∗ ∗ a ∗ ∗
perSPHparticle.Panels(d1)-(d3)showtheevolutionofthe distance and the apocenter distance, respectively, and we
√ √
specific energy in Models 1-9. The specific energy increases adopt that c =(−4/3)(2+ 6), c =(4 6−9), and c =
√ 1 2 3
with time, in stark contrast with the radiatively efficient (−4/3)(2 6−3).Theinitialpositionandvelocityvectorfor
case.PreferentialaccretionofhighlyboundSPHparticlesin- Model 10 are seen in Table 2.
creases,overtime,themeanspecificenergyofnon-accreted Figure 12 shows the evolution of specific binding en-
SPHparticles.Panels(c1)-(c3)showtheevolutionofspecific ergy and angular momentum. In panel (a), the blue solid,
angular momentum averaged over all remaining SPH parti- red solid, and black dashed lines show the specific angu-
cles. This also increases with time because of preferential larmomentaofModel10,ourpseudo-Newtoniansimulation
accretion of the lowest angular momentum SPH particles. model, and l respectively. From the panel, we note
pseudo
As a check on our radiative efficiency assumptions, we that the angular momentum is conserved for the pseudo-
compare the photon diffusion timescale with the timescale Newtonian case, while it is shifted at the ∼20% level from
forenergydissipationbyshockheating.Weestimatethislat- the Newtonian specific angular momentum and at the 3%
tertimescaleusingthedifferencebetweentheorbitalperiods levelevenfromthepseudo-Newtoniancase.Inpanel(b),the
ofthemosttightlyandlooselyboundgas,whicharegivenby blue solid, red solid, and black dashed lines show the spe-
equations (15) and (16), respectively. These timescales are cific binding energies of Model 10, our pseudo-Newtonian
shown for each model in Table 1. Model 1 has the longest simulation model, and (cid:15) respectively. From the panel,
pseudo
energy input timescale among all models, on the order of thespecificbindingenergyfirstagreeswellwiththatoftest
105s. From Figure 11, the surface density of the stellar de- particle, but substantially reduces due to debris circular-
bris and its scale height are estimated to be Σ/Σ ∼10−0.8 ization, and eventually saturates at ∼ 7P in the pseudo-
0 ∗
(cid:13)c 2014RAS,MNRAS000,1–??