A unified treatment of tidal disruption by Schwarzschild black holes

A unified treatment of tidal disruption by Schwarzschild black holes

Juan Servin Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA    Michael Kesden Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA
July 29, 2019
Abstract

Stars on orbits with pericenters sufficiently close to the supermassive black hole at the center of their host galaxy can be ripped apart by tidal stresses. Some of the resulting stellar debris becomes more tightly bound to the hole and can potentially produce an observable flare called a tidal-disruption event (TDE). We provide a self-consistent, unified treatment of TDEs by non-spinning (Schwarzschild) black holes, investigating several effects of general relativity including changes to the boundary in phase space that defines the loss-cone orbits on which stars are tidally disrupted or captured. TDE rates decrease rapidly at large black-hole masses due to direct stellar capture, but this effect is slightly countered by the widening of the loss cone due to the stronger tidal fields in general relativity. We provide a new mapping procedure that translates between Newtonian gravity and general relativity, allowing us to better compare predictions in both gravitational theories. Partial tidal disruptions in relativity will strip more material from the star and produce more tightly bound debris than in Newtonian gravity for a stellar orbit with the same angular momentum. However, for deep encounters leading to full disruption in both theories, the stronger tidal forces in relativity imply that the star is disrupted further from the black hole and that the debris is therefore less tightly bound, leading to a smaller peak fallback accretion rate. We also examine the capture of tidal debris by the horizon and the relativistic pericenter precession of tidal debris, finding that black holes of solar masses and above generate tidal debris precessing by or more per orbit.

I Introduction

Supermassive black holes (SBHs) with masses in the range from are found at the centers of most large galaxies Kormendy and Richstone (1995); Magorrian et al. (1998). Observations of quasars Schmidt (1963) and active galactic nuclei (AGN) have aided our understanding of these massive celestial objects, but these constitute only a small fraction Kauffmann et al. (2003) of all SBHs in the local universe, making it difficult to take a fully representative census of galactic nuclei. Studying quiescent SBHs, visible via their interactions with nearby stellar objects through the phenomena of tidal disruption events (TDEs) Rees (1988), would allow one to fill in this population gap, particularly for lower mass SBHs. Quiescent SBHs also have dynamical effects on neighboring stars and gas, but these can only be observed locally where the influence radius can be resolved.

TDEs occur when tidal forces due to a SBH’s gravitational field overcome an orbiting star’s self gravity, causing roughly half of the resulting debris to become more tightly bound to the SBH Rees (1988). This bound debris quickly forms a disk that is accreted by the SBH, releasing large amounts of energy in what is known as a tidal flare. Flares can be found in the optical van Velzen et al. (2011); Gezari et al. (2012); S. B. Cenko, et al. (2012a); I. Arcavi, et al. (2014); R. Chornock, et al. (2014); T. W.-S. Holoien, et al. (2014); J. Vinkó, et al. (2015), ultraviolet D. Stern, et al. (2004); S. Gezari, et al. (2006, 2008, 2009), soft X-ray Bade et al. (1996); Grupe et al. (1999); Donley et al. (2002a); P. Esquej, et al. (2008); Maksym et al. (2010); R. D. Saxton, et al. (2012), and hard X-ray J. S. Bloom, et al. (2011); D. N.Burrows, et al. (2011); S. B. Cenko, et al. (2012b) bands. A recent survey Donley et al. (2002b) discovered many candidate TDEs, showing rates consistent with those previously predicted by early models, resulting in an increased interest in TDEs in recent years.

Analysis of TDEs has been conducted under both Newtonian gravity Rees (1988); Magorrian and Tremaine (1999); Wang and Merritt (2004); Merritt (2013a); Cohn and Kulsrud (1978); Evans and Kochanek (1989); Lodato et al. (2009) and general relativity Luminet and Marck (1985); Laguna et al. (1993); Diener et al. (1997); Ivanov et al. (2003); Ivanov and Chernyakova (2006); Kesden (2012a, b); Sa̧dowski et al. (2016). In this paper, we provide a unified, self-consistent treatment of several aspects of the tidal-disruption process that are sensitive to general relativity. We explore how relativity changes the boundaries of the loss cone in phase space that determines TDE rates and the distribution of orbital elements for tidally disrupted stars. We then examine how relativity affects the distributions of the specific energy and angular momentum of the tidal debris and combine these results with recent hydrodynamical simulations to predict the rate at which this debris falls back onto the SBH to fuel a tidal flare Guillochon and Ramirez-Ruiz (2013).

We restrict ourselves in this paper to the Schwarzschild spacetime Schwarzschild (1916) corresponding to non-spinning SBHs. The spherical symmetry of the Schwarzschild metric implies that only the magnitudes of the orbital energy and angular momentum are needed to define the boundaries of the loss cone. However, in the Kerr spacetime Kerr (1963) of spinning SBHs, the orientation of the tidally disrupted star’s orbit with respect to the SBH’s equatorial plane affects the magnitude of the tidal stresses Ivanov and Chernyakova (2006). This creates a higher-dimensional parameter space, the analysis of which is beyond the scope of this paper. We will, however, provide some comments on SBH spin when relevant throughout the paper and lay the foundation for future work that will focus on incorporating spin dependence into our procedure.

The rest of this paper is structured as follows. Sec. II provides a brief overview of how loss-cone theory is used to calculate TDE rates in Newtonian gravity. Sec. III details how relativity changes the boundaries of the loss cone, allowing us to compare the TDE rates predicted in Newtonian gravity and general relativity. We then focus on trying to compare individual TDEs in the two gravitational theories, which introduces the somewhat subtle issue of which parameter to hold constant in these comparisons. We introduce a new procedure to map between Newtonian orbits and relativistic geodesics in Sec. IV that allows us to estimate how relativity would modify simulations performed in Newtonian gravity. In Sec. V, we use this mapping and loss-cone theory to predict distributions of physical quantities like the peak fallback accretion rate and amount of apsidal precession that affect observed TDE properties. A summary of our key results, their implications for TDE observations, and our plans to generalize to the Kerr spacetime are provided in Sec. VI.

Ii Newtonian loss-cone theory

We begin by providing a brief summary of Newtonian loss-cone theory Cohn and Kulsrud (1978); Magorrian and Tremaine (1999); Wang and Merritt (2004); Merritt (2013a). We define the tidal radius as

 rt≡(M∙M⋆)1/3R⋆ (1)

where is the mass of the SBH, is the mass of the tidally disrupted star, and is the stellar radius. This definition provides an order-of-magnitude estimate of the distance from the SBH at which tidal disruption occurs. Stars on parabolic orbits with pericenters equal to will have orbital angular momentum

 Lt≡(2GM∙rt)1/2 . (2)

Any star on an orbit with angular momentum will become tidally disrupted once it reaches a distance from the SBH. We can refine these estimates to account for stellar structure by introducing the parameter which measures the duration of the star’s pericenter passage relative to the hydrodynamical time of the star Press and Teukolsky (1977); Evans and Kochanek (1989). This parameter is a function of the polytropic index of the star’s equation of state Chandrasekhar (1939). We may then define

 rd ≡η2/3rt , (3a) Ld ≡η1/3Lt , (3b)

the radial distance and maximum orbital angular momentum at which a star described by parameter is tidally disrupted. Throughout this paper, we use the value consistent with the threshold for full tidal disruption in Guillochon and Ramirez-Ruiz Guillochon and Ramirez-Ruiz (2013) for a Solar-type star with .

The TDE rate thus depends on the rate at which stars are driven onto such orbits, which are said to lie in the loss cone in phase space because the velocity vector lies in a cone about the radial direction. We assume that the stellar distributions in galactic centers are described by isothermal spheres with density profile

 ρ(r)=σ22πGr2 (4)

where is the velocity dispersion given by the relation Gebhardt et al. (2000); Ferrarese and Merritt (2000). This is a reasonable approximation for galaxies with cusps in the surface brightness at their galactic centers Wang and Merritt (2004). We assume that loss-cone orbits are refilled by two-body non-resonant relaxation Cohn and Kulsrud (1978); Magorrian and Tremaine (1999); Wang and Merritt (2004); Merritt (2013a), in which case the differential flux of stars into the loss cone per unit dimensionless specific binding energy , following the approach and notation of Merritt Merritt (2013b), is given by

 F(ε∗)=4π2q(ε∗)L2d[∫yd0f(ε∗,y)dy] . (5)

Here is the orbital period divided by the time a star takes to diffuse across the loss cone, is a dimensionless angular-momentum variable, , and is the stellar phase-space density as a function of our dimensionless variables. For the stellar density profile of Eq. (4),

 q(ε∗)=32π23√2lnΛ(M⋆M∙)h(ε∗)ψ∗(rd)−ε∗(rdrh)−2 (6)

where is the Coulomb logarithm Chandrasekhar (1943a, b, c), is the negative dimensionless potential energy, is the influence radius Peebles (1972), and (Eq. (6.78) of Merritt Merritt (2013b)) is a function derived from the isotropic distribution function consistent with our isothermal stellar density profile . We set , appropriate for stars whose velocity distribution is consistent with the virial theorem Spitzer and Hart (1971); Wang and Merritt (2004).

The ratio helps provide intuition about how two-body relaxation feeds stars into SBHs. Portions of phase space with are said to belong to the empty loss-cone (ELC) regime because diffusion is not efficient enough to repopulate loss-cone orbits on their dynamical timescale. Conversely, corresponds to the full loss-cone (FLC) regime where diffusion is effective enough to keep orbits in these parts of phase space filled. SBHs with have for and thus belong primarily to the ELC regime, implying lower TDE rates despite their larger tidal radii Merritt (2013a). Approximately of TDEs occur in the FLC regime for the galaxy sample considered in Stone and Metzger Stone and Metzger (2016), where result in full rather than partial disruptions.

The integrand in Eq. (5) for the isotropic distribution function consistent with the isothermal density profile of Eq. (4) is given by

 f(y)=f(yd)[1−2√yd∞∑m=1e−γ2m/4γmJ0(γm√y)J1(γm√yd)] (7)

where and are Bessel functions of the first kind and is defined such that is the -th zero of Merritt (2013a). The flux of stars into the loss cone is then,

 F(ε∗)=4π2L2df(yd)ξ(q) (8)

where is defined to be times the integral in Eq. (5) and is well approximated by . The total TDE rate is obtained by integrating over all binding energies,

 ˙N=∫F(ε∗)dε∗, (9)

and depends implicitly on the SBH mass .

Iii Relativistic tidal forces

In Newtonian gravity, the tidal force acting on a fluid element of a star is simply the gravitational force in a frame freely falling with the star’s center of mass. In general relativity however, gravity affects particle motion through the curvature of spacetime. In this section, we review how tides arise in general relativity so that we can determine the value of the angular momentum that sets the boundary of the loss cone in this gravitational theory. We restrict our analysis in this paper to the Schwarzschild metric Schwarzschild (1916) describing non-spinning SBHs, which in Boyer-Lindquist coordinates Boyer and Lindquist (1967) is

 ds2 =−(1−2M∙r)dt2+r2r2−2M∙dr2 +r2(dθ2+sin2θdϕ2). (10)

Massive test particles experiencing no non-gravitational forces travel on timelike geodesics of this metric, the generalization of straight lines in flat space. We can choose to parameterize these geodesics by the proper time . Particles traveling on these geodesics have conserved specific energy and specific angular momentum since the Schwarzschild metric is both time independent and spherically symmetric. Note that the specific energy

 E≡−gμνUμ(∂∂t)ν=(1−2M∙r)dtdτ, (11)

where

 Uμ=dxμdτ=(dtdτ,drdτ,dθdτ,dϕdτ) (12)

is the 4-velocity, asymptotes to unity for particles at rest far from the SBH since it contains the rest-mass energy. We approximate for the orbits of tidally disrupted stars because the velocity dispersion of their host galaxies is much less than the speed of light. This 4-velocity satisfies the geodesic equation

 dUμdτ+ΓμναUνUα=0 (13)

where are the Christoffel symbols for the Schwarzschild metric. Note the resemblance to Newton’s second law for force-free motion for vanishing Christoffel symbols as can be chosen for flat space.

In general relativity, tidal forces arise because of the tendency of parallel geodesics to deviate from each other in the presence of spacetime curvature. If two neighboring particles with 4-velocity are separated by an infinitesimal spacelike deviation 4-vector , this deviation will evolve with proper time according to the geodesic deviation equation

 d2Xβdτ2 =Uμ∇μ(Uα∇αXβ)=−Rβ μανUμXαUν =−Cβ αXα (14)

where are covariant derivatives, is the derivative with respect to the proper time, is the Riemann curvature tensor, and

 Cβ α≡Rβ μανUμUν (15)

is the tidal tensor. In Newtonian gravity, a fluid element is tidally stripped from the surface of a star when the tidal force exerted by the SBH exceeds the force exerted on this fluid element by that star’s self gravity. If the star is non-degenerate like the Solar-type stars considered in this paper, we can approximate its self gravity as a Newtonian force even if general relativity is needed to describe the geodesics of the SBH. In this approximation, the fluid element is tidally stripped when the acceleration due to this self gravity is exceeded by the proper acceleration given by the geodesic deviation equation (14) Luminet and Marck (1983).

We can solve this equation more easily by transforming from the global Boyer-Lindquist coordinates given by Eq. (10) to local Fermi normal coordinates () valid in a neighborhood of spacetime about the center of mass of the star at an arbitrarily chosen proper time Manasse and Misner (1963). These coordinates define an orthonormal tetrad: the star’s 4-velocity provides the timelike 4-vector , and three spacelike 4-vectors are chosen that are parallel transported along the geodesic. Spatial indices that run from 1 to 3 are denoted by Latin indices, unlike Greek indices that run from 0 to 3. All of the tensors appearing in the geodesic deviation equation can be projected into this basis:

 Uμ =λμ (0) (16a) Xα =X(i)λα (i) (16b) Rβ μαν =R(γ)(δ)(κ)(ξ)λβ (γ)λ (δ)μλ (κ)αλ (ξ)ν (16c) Cβ α =C(i)(j)λβ (i)λ (j)α=R(i)(0)(j)(0)λβ (i)λ (j)α. (16d)

The Boyer-Lindquist indices are raised and lowered by the Schwarzschild metric (10) while the Fermi normal coordinate indices in brackets are raised and lowered by the Lorentz metric . Note that the symmetry of the Riemann tensor implies that the tidal tensor is symmetric and only has spatial components in Fermi normal coordinates. Inserting Eq. (16) into (14) yields the geodesic deviation equation in Fermi normal coordinates

 d2X(i)dτ2=−C(i)(j)X(j). (17)

Because is a real, symmetric tensor, it has 3 real eigenvalues which in Boyer-Lindquist coordinates are , , and . The negative sign in Eq. (17) implies that the eigenvectors associated with the positive eigenvalues correspond to directions along which the star is compressed, while the eigenvector associated with the negative eigenvalue corresponds to the direction in which the star is stretched. Equating the magnitude of this negative eigenvalue times the stellar radius with the acceleration due to the star’s self gravity yields

 (1+3L22r2)2M∙R⋆r3=2M⋆(ηR⋆)2. (18)

Combining this result with the relation

 L2=2M∙r2r−2M∙, (19)

between the orbital angular momentum and Boyer-Lindquist coordinate at pericenter provides two equations that can be solved for the the relativistic tidal radius and angular-momentum threshold for tidal disruption. Note that in the limit these reduce to the Newtonian expressions of Eqs. (2) and (3).

Fig. 1 shows how the angular-momentum threshold depends on SBH mass in both Newtonian gravity and general relativity. All orbits with angular momentum below this threshold are considered to be inside the loss cone. The Newtonian tidal radius according to Eq. (1), while the gravitational radius scales linearly with SBH mass. This implies that relativistic effects become negligible for small SBH masses and the two curves in Fig. 1 converge towards the left edge of the plot. The fact that the term in brackets on the left-hand side of Eq. (18) is greater than unity implies that tidal forces are stronger in general relativity than in Newtonian gravity for orbits with the same angular momentum. This explains why the relativistic curve in Fig. 1 is above the Newtonian curve; stronger tides allow the SBH to disrupt stars at larger angular momentum in relativity.

Relativity also introduces the phenomena of stellar capture; stars on orbits with angular momentum less than will plunge directly into the event horizon leaving no debris to emit photons and generate an observable TDE. This angular-momentum threshold for direct capture, shown by the horizontal red dashed line in Fig. 1, grows linearly with SBH mass and eventually exceeds the threshold for tidal disruption at . This prevents observable TDEs by SBHs with masses above this value as shown by the intersection of the solid green and dashed red curves in Fig. 1.

Fig. 2 shows the rates for observable TDEs and direct capture predicted by Eq. (5) when two-body relaxation is responsible for refilling the loss cone. Relativity modifies the Newtonian prediction in two ways: (1) it changes the thresholds and in Eqs. (6) and (7) to the maximum of the two thresholds set by tidal disruption and direct capture, and (2) the limits of the integral in Eq. (5) must be modified. If (), the integral must be decomposed into two parts: an integral from 0 to giving the capture rate and an integral from to giving the observable TDE rate. For (), there are no observable TDEs and a single integral from 0 to gives the capture rate.

Iv Comparing tidal disruption in Newtonian gravity and general relativity

Having introduced tidal disruption in Newtonian gravity and general relativity in Secs. II and III, we now seek to compare the predictions of the two theories. This point is more subtle than it might initially appear, as there is no unique way to map Keplerian orbits in Newtonian gravity to geodesics of the Schwarzschild metric in general relativity. In the limit that the pericenter , it is natural to map a geodesic to its identical Keplerian counterpart, but this limiting behavior is insufficient to fully specify the mapping. The appropriate mapping to use depends on the nature of the problem one is trying to solve. We consider in this section three distinct mappings that all possess the desired behavior in the Newtonian limit; these three mappings identify:

• orbits with equal pericenter coordinates ,

• orbits with equal angular momenta ,

• orbits on which a star experiences equal tidal forces at pericenter.

Mapping (1) is perhaps the most obvious choice and was used for the corrections to the orbital constants provided in Kesden Kesden (2012b). However, according to the Schwarzschild metric in Boyer-Lindquist coordinates given in Eq. (10), the physical significance of this mapping is that it identifies parabolic orbits such that the circular orbits in the two gravitational theories with the same pericenter would have the same circumference . It is unclear why this choice of mapping would be particularly useful for analyzing the tidal disruption of stars on non-circular orbits.

Mapping (2) identifies orbits in the two theories with the same values of the gauge-invariant orbital angular momentum, defined for the Schwarzschild metric as

 L≡gμνUμ(∂∂ϕ)ν=r2dϕdτ (20)

for equatorial orbits (). This mapping seems like a useful choice, as TDE properties do depend on the orbital angular momentum as seen in the recent simulations of Guillochon and Ramirez-Ruiz Guillochon and Ramirez-Ruiz (2013). However, tides are stronger on Schwarzschild geodesics than on Keplerian orbits with the same value of as established in Sec. III, so it seems unlikely that TDEs on orbits in the two theories identified in this manner would have the same properties.

We conjecture that TDEs resulting from stars initially on orbits identified by mapping (3) will have similar properties because the stars were subjected to similar tidal forces. This conjecture is supported by the ”freezing” model of Lodato, King, and Pringle Lodato et al. (2009) which proposed that the orbital energy distribution of tidal debris is frozen in following an instantaneous tidal disruption of an unperturbed star. This model was used with reasonable success to describe the light curve of the TDE PS1-10jh Gezari et al. (2012). We will use this freezing assumption later in Sec. V to determine the appropriate relativistic correction to the energy distribution, but we will apply the correction at the disruption radius rather than pericenter. This is consistent with the analytic model of Stone, Sari, and Loeb Stone et al. (2013) and some hydrodynamical simulations of TDEs Hayasaki et al. (2013); Guillochon and Ramirez-Ruiz (2013) which showed that the spread in debris energy was largely independent of the orbital angular momentum for full disrupted stars.

To implement mapping (3), we must find the angular momentum of the orbit in Newtonian gravity that has the same peak tidal force as that experienced by a star on a Schwarzschild geodesic with angular momentum and pericenter . We accomplish this by setting the magnitude of the negative eigenvalue of the tidal tensor equal to its Newtonian limit for an orbit with angular momentum :

 (1+3L22r2)1r3=(2M∙L2N)3. (21)

This equation, combined with Eq. (19) relating the angular momentum to the pericenter in Boyer-Lindquist coordinates, can be solved to determine for mapping (3).

Instead of using the angular momentum to identify orbits, the TDE literature often uses the penetration factor

 β≡L2tL2, (22)

where is defined in Eq. (2). This choice is convenient because stars on orbits with are tidally disrupted, where is of order unity. Our mapping implies an equivalent mapping , where

 βN≡L2tL2N=(r+M∙r−2M∙)1/3rtr (23)

and is an implicit function of and hence through Eq. (19). We show this mapping for several values of the SBH mass in Fig. 3. These mappings are undefined for since such geodesics plunge directly into the event horizon and therefore do not have pericenters at which one can calculate the tidal force. The mappings are above the diagonal and have positive curvature because tides are stronger in general relativity, requiring Newtonian orbits to have deeper penetration factors to match the tides on Schwarzschild geodesics at pericenter. The mapping for SBH mass terminates at , precisely the threshold for the full disruption of Solar-type stars Guillochon and Ramirez-Ruiz (2013) seen in Fig. 1 for .

To further illustrate how relativity affects the minimum penetration factor for tidal disruption, we show its dependence on SBH mass in Fig. 4 for both full and partial tidal disruptions. Newtonian hydrodynamical simulations Guillochon and Ramirez-Ruiz (2013) indicate that Solar-type stars are fully disrupted for , while partial disruptions occur for . As increases and approaches , the curves fall increasingly below these Newtonian thresholds because of the stronger tides in general relativity. The values of at which these curves intersect the curve are the maximum masses and capable of full and partial tidal disruption. These limits are consistent with recent work suggesting that observable TDEs by SBHs with masses above result exclusively from giant stars MacLeod et al. (2012). Determinations of that failed to account for the stronger relativistic tides, effectively using mapping (2) shown by the horizontal lines in Fig. 4, would underestimate by a factor of for both full and partial disruptions. Future work will explore how SBH spin can push to even higher values for stars on prograde orbits.

V Distributions of relativistic TDE properties

In this section, we use concepts developed earlier in the paper to predict distributions of several quantities that affect the observed properties of TDEs. Guillochon and Ramirez-Ruiz Guillochon and Ramirez-Ruiz (2013) performed Newtonian hydrodynamical simulations to determine the peak accretion rate at which tidal debris falls back on the SBH following a TDE and the delay between the tidal disruption and the time at which this peak accretion occurs. They then provided analytic fits to these quantities as functions of the Newtonian penetration factor . Under the approximation that these quantities depend only on the peak tidal forces along each orbit, we use the inverse of the mapping shown in Fig. 3 to derive fits to and as functions of the relativistic penetration factor in the Schwarzschild spacetime. We then use the loss-cone theory discussed in Sec. II, modified by the relativistic corrections to the boundaries of the loss cone shown in Fig. 1, to predict the distributions of these quantities for different SBH masses.

In Subsec. V.1, we determine the differential TDE rate per unit relativistic penetration factor . In Subsec. V.2, we use this rate and our corrected fits and to derive initial estimates for the differential TDE rates and per unit peak fallback accretion rate and time delay respectively. These initial estimates do not account for an additional relativistic effect discussed in Kesden Kesden (2012b), that the gradient of the potential well that determines the width of debris energy distribution differs in general relativity and Newtonian gravity. We review this result and then derive revised estimates for and in Subsec. V.3. The stellar debris produced following tidal disruption will also have a distribution of orbital angular momentum, but this often receives less attention in Newtonian gravity where the orbital periods that determine the fallback accretion rate depend on energy but not angular momentum. However, the distribution of orbital angular momentum can be very important when considering highly relativistic TDEs for which much of the debris can lose enough specific angular momentum to fall below the threshold for capture even if the initial star was above this threshold. We discuss this regime in Subsec. V.4 and show that this effect sharply limits the relativistic correction to the peak fallback accretion rate. Finally, in Subsec. V.5 we use our relativistically corrected loss-cone theory to determine the differential TDE rate per unit relativistic shift in argument of pericenter, an important effect responsible for tidal stream crossings and the prompt circularization of debris Hayasaki et al. (2013); Dai et al. (2013); Guillochon et al. (2014); Guillochon and Ramirez-Ruiz (2015); Shiokawa et al. (2015).

v.1 Distribution of penetration factor β

We begin by calculating the differential TDE rate in Newtonian gravity. The angular-momentum variable can be expressed in terms of the penetration factor as

 y=L2tqL2dβ . (24)

This relation allows us to differentiate Eq. (5) with respect to , which when combined with the distribution function in Eq. (7) yields

 dFdβ =(2πLtβ)2f(yd)× [1−2√q∞∑m=1e−γ2m/4γmJ0(γmLt/Ld√qβ)J1(γm/√q)]. (25)

Integrating this expression with respect to the dimensionless specific binding energy gives us the Newtonian differential TDE rate shown in the left panel of Fig. 5 above. The total TDE rate is the area under this curve for which decreases with increasing SBH mass consistent with the TDE rates seen in Fig. 2. The blue curve in the left panel of Fig. 5 extends smoothly beyond because much of the phase space beyond this value is still in the full loss-cone regime for which Stone and Metzger (2016). However, as increases towards , the differential TDE rate becomes exponentially suppressed for reflecting that much of the phase space beyond this point now lies in the empty loss-cone regime Stone and Metzger (2016).

Relativity introduces two changes into the calculation of the differential TDE rate : (1) the stronger tides increase the numerical value of from the Newtonian result given by the black curve in Fig. 1 to the relativistic result given by the green curve in this figure, and (2) capture by the event horizon cause the differential rate to fall discontinuously to zero for . As most stars diffusing into the loss cone have apocenters near the influence radius Frank and Rees (1976), the Newtonian treatment of stellar diffusion encapsulated in the ratio remains an accurate approximation. Both of the changes listed above leave signatures in the relativistic differential TDE rate seen in the right panel of Fig. 5. The bends in these curves marking the threshold for full disruption in the empty loss-cone regime migrate to lower values of consistent with the SBH mass dependence of the mapping shown in Fig. 3. The curves also end abruptly at as can be seen in the red curve corresponding to SBH mass for which .

v.2 Distributions of peak fallback accretion rate ˙Mpeak and time delay tpeak

Guillochon and Ramirez-Ruiz Guillochon and Ramirez-Ruiz (2013) performed a series of Newtonian hydrodynamical simulations at different penetration factors and used them to derive analytic fits to the peak fallback accretion rate and time delay for Solar-type stars with polytropic index :

 ˙Mpeak =A4/3(M∙106M⊙)−1/2, (26a) tpeak =B4/3(M∙106M⊙)1/2, (26b)

where,

 A4/3 =exp(27.3−27.5βN+3.87β2N1−3.26βN−1.39β2N), (27a) B4/3 =−0.387+0.573√βN−0.312βN1−1.27√βN−0.9βN, (27b)

and . Given these fits, the mapping derived in Sec. IV, and the differential TDE rate obtained in the previous subsection, we can calculate differential rates

 d˙NdX=d˙Ndβ(dβNdβdXdβN)−1 (28)

for . We show these differential TDE rates in Fig. 6 below. A complication arises from the fact that the dependence on penetration factor in Eq. (27) is not monotonic; () increases (decreases) with until the threshold for full disruption is reached, then decreases (increases) for TDEs that penetrate more deeply into the SBH’s potential well. This dependence contradicts the naive predictions of freezing models which assume that the energy distribution of the tidal debris is frozen in at the disruption radius , but does not fatally compromise our analyis. We address this complication by plotting the TDEs from partial () and full () disruptions separately in the left and right panels of Fig. 6. The differential rates diverge at the extremum shown by the vertical black dashed lines in these plots, but the total rates (areas under the curves) remain finite. We see that the stronger tides of relativity suppress the rates of both full and partial TDEs, particularly for SBH masses for which portions of phase space with lie primarily in the empty loss-cone regime. This emptiness even reduces the rate of partial TDEs because of the sharper gradients driving stronger diffusion across the boundary as seen by the dip near the threshold in the top left panel of Fig. 6.

v.3 Relativistic correction to the debris energy distribution

The estimated differential TDE rates and in the previous subsection accounted for relativistic corrections to the boundaries of the loss cone, but in the interest of a systematic exploration of relativistic effects we did not simulataneously include relativistic corrections to the energy distribution of the tidal debris. We turn our attention to this effect in this subsection.

A key approximation often used in the analysis of TDEs is that the tidal debris travels on ballistic trajectories following disruption, and that orbit circularization followed by viscous accretion occurs promptly after the bound debris falls back to pericenter Rees (1988); Evans and Kochanek (1989); Kochanek (1994). The fallback accretion rate onto the SBH, and thus its peak value and the time between disruption and when this peak occurs, is determined by the distribution of orbital periods of the tidal debris, which is in turn set by the energy distribution through Kepler’s third law

 τ=2π(a3M∙)1/2=2πM∙(2E)−3/2. (29)

This Newtonian relation remains an excellent approximation even for relativistic TDEs, because most of the debris is on highly eccentric orbits prior to circularization and spends most of its time near apocenter where Eq. (29) holds. However, if the tidal disruption itself occurs near pericenter where relativistic effects are strongest, Newtonian predictions for the specific binding energy entering into Eq. (29) may not be accurate. TDE simulations in general relativity Hayasaki et al. (2013); Cheng and Bogdanović (2014); Bonnerot et al. (2016); Shiokawa et al. (2015); Cheng and Evans (2013) naturally yield proper relativistic energy distributions, but our goal in this section is to study how relativistic corrections might alter the predictions of Newtonian simulations.

We addressed this issue in Kesden Kesden (2012b), where we found that for an undistorted star of radius , the width of the potential-energy distribution before disruption and thus the width of the debris energy distribution after disruption is given in general relativity by the expression

 σE,GR=|λα (i)∇αE|R⋆=|gβγλβ (0)λα (i)Γγαt|R⋆, (30)

where is the metric tensor, are the Christoffel symbols, and and are the orthonormal tetrad of basis 4-vectors determined by our choice of Fermi normal coordinates in Sec III. For corresponding to partial disruptions, tidal debris is assumed to be liberated at pericenter where the tidal forces are strongest. Evaluating Eq. (30) at pericenter for such orbits yields

 σE,GR,PD =M∙R⋆r2(1−2M∙r)−1/2 =Et(rtr)2(1−2M∙r)−1/2. (31)

The pericenter in Boyer-Lindquist coordinates can be expressed in terms of the angular momentum using Eq. (19),

 Et≡M∙R⋆r2t=(M∙M⋆)1/3M⋆R⋆=(M∙M⋆)1/3E⋆ (32)

is an order-of-magnitude estimate for the width of this energy distribution, and is a similar estimate for the specific self-binding energy of the star. The mass hierarchy implies that , supporting the assumption of freezing models that the relative velocities of fluid elements at the time of disruption can be neglected when determining the debris energy distribution.

In trying to compare this relativistic result to Newtonian predictions, we encounter the same issue of choosing a mapping between the two gravitational theories that we addressed in Sec. IV. Mapping (1), which was used in Kesden Kesden (2012b), compared Eq. (31) to Newtonian orbits with the same pericenter coordinate and found

 σE,N,PD(1)=M∙R⋆r2=Et(rtr)2. (33)

Dividing Eq. (31) by Eq. (33) yields a peak correction of for , the minimum pericenter for an orbit that avoids direct capture by the horizon. As argued in Sec. IV however, this mapping is a somewhat unnatural way to identify orbits in the two theories.

Mapping (2) compared Schwarzschild geodesics to Newtonian orbits with the same angular momentum (and hence penetration factor ). For this choice, the width of the Newtonian energy distribution is

 σE,N,PD(2)=M∙R⋆r2L=Et(2M∙L2)2=β2Et , (34)

where is the pericenter of a Newtonian orbit with angular momentum . Freezing models posit that the energy distribution is frozen in for , implying that for full disruptions in mapping (2). Although there is no nice analytic expression for the widths for full disruption in general relativity, they can be calculated numerically using Eq. (30). We show the widths of these distributions in mapping (2) for both partial and full disruptions in the left panel of Fig. 7. We see that for partial disruptions , the energy distribution is broader in relativity than Newtonian gravity. The ratio reaches a maximum value of at , even larger than the peak correction in mapping (1). However, the weaker tides in Newtonian gravity imply that stars can reach larger penetration factors before being fully disrupted. SBHs with masses can have deeply penetrating encounters where the star still manages to avoid direct capture. For such penetration factors, the relativistic prediction can fall below that in Newtonian gravity as can be seen near the right edge of the left panel of Fig. 7. It is interesting to note that the relativistic predictions retain some mild dependence for full disruptions , unlike in Newtonian freezing models. This is because the tidal tensor and energy width are velocity-dependent in relativity, even if evaluated at the fixed tidal acceleration set by the star’s self gravity.

Mapping (3) solves the problem of stars being fully disrupted at different values of in the two theories by identifying the Schwarzschild geodesic with penetration factor with the Newtonian orbit with penetration factor on which a star experiences the same peak tidal force at pericenter as shown in Fig. 3. In this case, the width of the energy distribution for partial disruptions is

 σE,N,PD(3)=M∙R⋆r2F=β2N(β)Et=M∙R⋆r2(r+M∙r−2M∙)2/3 , (35)

where is the pericenter of a Newtonian orbit with penetration factor and as previously is the pericenter in Boyer-Lindquist coordinates of a Schwarzschild geodesic with penetration factor . This width again freezes out at for full disruptions as in mapping (2), but these full disruptions now occur for as in general relativity. We show and as functions of penetration factor in the right panel of Fig. 7. We see that full disruption occurs at in both theories with this mapping, but that for a given penetration factor , the relativistic predictions are now below the Newtonian predictions . The deeper penetration needed for tidal disruption in Newtonian gravity leads to steeper potential gradients and thus broader tidal debris energy distributions. The ratio reaches a minimum value of at . The Newtonian orbit mapped to the Schwarzschild geodesic with has , but this is not a problem in principle since Newtonian point masses have no horizons to capture stars.

Which of the three expressions given by Eqs. (33), (34), (35) is the ”right” one to use when comparing TDEs in Newtonian gravity with those in general relativity? The answer to this question depends on which of the three mappings discussed in Sec. IV you are using to relate orbits in the two theories. If you want to compare orbits with the same angular momentum , Eq. (34) provides the width of the tidal energy distribution for a Newtonian orbit with the same penetration factor , leading to more tightly bound debris in relativity for partial disruptions as seen in the left panel of Fig. 7. However, this choice implies that different fractions of the stellar mass are stripped away by tides as seen in Fig. 8. More material is lost by the partially disrupted star in relativity, further enhancing the fallback accretion rate in addition to the more tightly bound material falling back more quickly. For , full disruptions occur in general relativity but not Newtonian gravity, while for stars are fully disrupted in both theories but the tidal debris is more tightly bound in Newtonian gravity.

If you are instead trying use Newtonian TDE simulations like those in Guillochon and Ramirez-Ruiz Guillochon and Ramirez-Ruiz (2013) to predict what this process might be like in relativity, you should use mapping (3) to relate orbits experiencing the same peak tidal forces and thus the same amount of mass loss by partially disrupted stars. Eq. (35) should be used to determine the relativistic correction as shown in the right panel of Fig. 7. We use this correction to predict the peak fallback accretion rate and time delay

 ˙Mpeak(β) =(σE,GRσE,N(3))3/2×˙Mpeak[βN(β)], (36a) tpeak(β) =(σE,GRσE,N(3))−3/2×tpeak[βN(β)], (36b)

in general relativity, where the exponent of 3/2 follows from the energy dependence in Kepler’s third law (29). We show the corrected differential TDE rates and in Fig. 9. We see that the peak fallback accretion rates have been reduced below their Newtonian values by for a SBH, consistent with Eq. (36) and the ratios of the curves shown in the right panel of Fig. 7. The time delays show a corresponding increase. Fig. 9 conveys the qualitative message of this subsection: the stronger tides in general relativity allow SBHs to tidally disrupt stars at lower penetration factors leading to less tightly bound debris and lower fallback accretion rates.

v.4 Capture of tidal debris by the event horizon

In this subsection, we focus on an issue that was only briefly addressed in Kesden Kesden (2012b): tidal debris can be captured by the event horizon of an SBH even if the disrupted star is not initially on a capture orbit. In general relativity, the tidal debris will have a distribution of specific angular momentum with a width given by

 σL,GR=|λα (i)∇αL|R⋆=|gβγλβ (0)λα (i)Γγαϕ|R⋆ (37)

analogous to Eq. (30) giving the width of the specific energy distribution. This width is set at pericenter for partial disruptions, allowing us to evaluate Eq. (37) as

 σL,GR,PD=R⋆(2M∙r)1/2=Lt(M⋆M∙)1/3(rtr)1/2. (38)

The mass hierarchy between the star and SBH implies that and thus that it is usually a good approximation to assume that the specific angular momentum of the tidal debris is equal to that of the initial star. This contrasts with the debris energy distribution for which implying that the tidal debris is much more tightly bound to the SBH than the initial star. However, for highly relativistic TDEs, the small amount of specific angular momentum lost in the disruption process may be enough for some of the tidal debris to be captured by the horizon. We examine this possibility by defining the dimensionless parameter

 fL≡L+ΔL−Lcap|ΔL|=L−4M∙|ΔL|−1, (39)

where is the specific angular momentum of the initial star, is the change in specific angular momentum of a fluid element in the disruption process, and is the specific angular momentum threshold for direct capture by the horizon. An element of the tidal debris with will be captured by the SBH, and is the minimum value of this parameter for a star that is not originally on a capture orbit ().

We need an estimate for to evaluate our new parameter . In the freezing model, an element of the tidal debris located at in Fermi normal coordinates at tidal disruption has its specific energy and angular momentum changed by an amount

 ΔE ≡X(i)rE(i)=X(i)λα (i)∇αE, (40a) ΔL ≡X(i)rL(i)=X(i)λα (i)∇αL, (40b)

in the tidal-disruption process. The element of the star that will become the most tightly bound element of the tidal debris will have of magnitude anti-aligned with , leading to the specific binding energy given by Eq. (30). The element of the star that loses the most angular momentum in the tidal-disruption process will have of magnitude anti-aligned with and have its angular momentum reduced by an amount given by Eq. (37). If and are parallel to each other (as at pericenter in Newtonian gravity, where both point in the radial direction), the same element of the tidal debris will be both the mostly tightly bound (the first to fall back onto the SBH) and have lost the most angular momentum (and thus be at greatest risk for direct capture). We will assume this to be true, as it is exactly true in general relativity at pericenter and a good approximation for those fluid elements with most at risk of direct capture. Under this assumption, TDEs for which for in Eq. (39) will have suppressed fallback accretion rates and greater delays between disruption and the beginning of emission because direct capture will have removed the most tightly bound debris. Furthermore, if and are parallel, there is a one-to-one relationship

 ΔL=|rL(i)||rE(i)|E (41)

between the specific binding energy of an element of tidal debris and the specific angular momentum it loses during tidal disruption. If we use Kepler’s third law (29) to relate the time delay to the specific binding energy of debris accreted at that time, we can use Eq. (41) to estimate and thus for such debris. When this estimate of becomes negative, all of the debris accreted before the fallback accretion rate reaches its peak will be captured by the horizon.

We show for both peak and most tightly bound fluid elements as a function of SBH mass in Fig. 10, for stars with penetration factors and corresponding to the thresholds for full disruption and nonzero partial disruption. For the most tightly bound elements with , Eq. (39) yields

 fL=(M∙M⋆)1/3(rrt)1/2[1√β−(M∙M⋆)1/3√8E⋆]−1. (42)

For , the term in the square brackets goes to and , however as increases this term decreases and ultimately reaches zero () for when the entire star is directly captured. This increase then subsequent decrease in as a function of is seen in the left panel of Fig. 10 indicating that tidal debris is safe from direct capture for all but the most massive SBHs. Only when is within of does the direct capture of tidal debris become significant as seen in the right panel of Fig. 10. Although such events constitute only a small fraction of the total TDE rate, they are the TDEs subject to the most extreme relativistic effects such as the pericenter precession that will be considered in the next subsection.

v.5 Relativistic pericenter precession

Early work on TDEs assumed that tidal debris would promptly circularize after falling back to pericenter Rees (1988), but Newtonian hydrodynamical simulations demonstrated that energy dissipation at pericenter was inefficient for Guillochon et al. (2014). Early analytic work suggested that relativistic pericenter precession might promote orbit circularization, because this precession would lead to steam crossings at which the inelastic collisions of tidal elements would transform orbital kinetic energy into heat that could be subsequently radiated away. This suggestion was later supported by relativistic TDE simulations in which such precession was indeed shown to generate tidal stream crossings Hayasaki et al. (2013); Dai et al. (2013); Guillochon and Ramirez-Ruiz (2015); Shiokawa et al. (2015).

The angular coordinate specifying the location of pericenter with respect to a reference axis in the orbital plane is known as the argument of pericenter. In Newtonian gravity, the argument of pericenter is a constant of motion, but this is not true for Schwarzschild geodesics of the metric (10) in Boyer-Lindquist coordinates. At lowest post-Newtonian (PN) order, the argument of pericenter changes by an amount

 Δω1PN=6π(M∙L)2=6πβ(M∙M⋆)2/3E⋆ . (43)

This PN approximation breaks down for for which pericenter precession must be integrated numerically along Schwarzschild geodesics,

 ΔωS=2∫∞r˙ϕ