- Split View
-
Views
-
Cite
Cite
Kimitake Hayasaki, Nicholas Stone, Abraham Loeb, Finite, intense accretion bursts from tidal disruption of stars on bound orbits, Monthly Notices of the Royal Astronomical Society, Volume 434, Issue 2, 11 September 2013, Pages 909–924, https://doi.org/10.1093/mnras/stt871
- Share Icon Share
Abstract
We study accretion processes for tidally disrupted stars approaching supermassive black holes on bound orbits, by performing three-dimensional smoothed particle hydrodynamics simulations with a pseudo-Newtonian potential. We find that there is a critical value of the orbital eccentricity below which all the stellar debris remains bound to the black hole. For high but subcritical eccentricities, all the stellar mass is accreted on to the black hole in a finite time, causing a significant deviation from the canonical t−5/3 mass fallback rate. When a star is on a moderately eccentric orbit and its pericentre distance is deeply inside the tidal disruption radius, there can be several orbit crossings of the debris streams due to relativistic precession. This dissipates orbital energy in shocks, allowing for rapid circularization of the debris streams and formation of an accretion disc. The resultant accretion rate greatly exceeds the Eddington rate and differs strongly from the canonical rate of t−5/3. By contrast, there is little dissipation due to orbital crossings for the equivalent simulation with a purely Newtonian potential. This shows that general relativistic precession is crucial for accretion disc formation via circularization of stellar debris from stars on moderately eccentric orbits.
INTRODUCTION
Supermassive black holes (SMBHs) larger than 105 M⊙ are now known to exist ubiquitously in galactic nuclei. SMBHs in nearby galaxies can be studied dynamically, whereas in more distant galaxies only the ∼1 per cent of SMBHs undergoing major accretion episodes can be easily observed. The tidal disruption and subsequent accretion of a star by a SMBH, although intrinsically a rare event, is of observational interest as a way to probe dormant SMBHs in the local Universe, because it can produce a powerful flare in excess of the Eddington luminosity (Carter & Luminet 1983; Rees 1988; Evans & Kochanek 1989; Strubbe & Quataert 2009). These events are also of interest to low-frequency gravitational wave astronomy, as the stellar dynamical processes which funnel stars into the low angular momentum orbits necessary for tidal disruption events (TDEs) are similar to those which produce the extreme mass ratio inspiral of a stellar-mass compact object on to an SMBH (Frank & Rees 1976; Wang & Merritt 2004; Madigan, Hopman & Levin 2011). TDEs have also been considered as part of the means by which a seed black hole grows into an SMBH (Zhao, Haehnelt & Rees 2002; Miralda-Escudé & Kollmeier 2005; Bromley et al. 2012). All of these aspects of stellar tidal disruption have motivated past studies of TDEs.
The traditional picture of a TDE proceeds as follows: a star at large separation (∼1 pc) approaches an SMBH on a nearly parabolic orbit. After the star is tidally disrupted by the black hole, half the stellar debris becomes gravitationally bound to the SMBH, because it loses orbital energy inside the tidal disruption radius. The bound debris falls back, and, after circularizing due to collisional shocks with other gas streams, accretes on to the black hole. Kepler's third law implies that the mass return rate decays with a −5/3 power of time (Rees 1988; Phinney 1989), asymptotically approaching zero. This simple analytic picture has been validated to the first order of approximation by hydrodynamical simulations (Evans & Kochanek 1989), albeit with deviations from this law at early times (Lodato, King & Pringle 2009). Similar power-law behaviour for the flare light curve is often assumed (i.e. |$L \propto \dot{M}$|), although here the theoretical evidence is less clear (Lodato & Rossi 2011).
All-sky surveys in the X-ray and UV have so far observed 13 candidate tidal disruption flares (Komossa & Bade 1999; Grupe, Thomas & Leighly 1999; Greiner et al. 2000; Komossa et al. 2004; Halpern, Gezari & Komossa 2004; Dogiel et al. 2009; Maksym, Ulmer & Eracleous 2010; Bloom et al. 2011; Burrows et al. 2011). The observed light curves are in reasonable agreement with the theoretically predicted mass fallback rate of t−5/3, although some show deviations (Burrows et al. 2011) and the number of samples is sufficiently small to make detailed testing of theoretical models difficult. A smaller number of strong TDE candidates have been found at optical wavelengths (van Velzen et al. 2011; Gezari et al. 2012). Notably, two of the best-sampled TDEs differ strongly from the canonical theoretical picture: one possesses relativistic jets (Bloom et al. 2011; Zauderer et al. 2011), and the other lacks hydrogen lines in its spectra (Gezari et al. 2012).
It has been inferred from observations that the event rate of tidal disruption is ∼10−5 yr−1 per galaxy (Donley et al. 2002). This observed rate is in rough agreement with uncertain theoretical rate estimates. In the standard scenario, stars are supplied to the SMBH by two body scattering, on a stellar relaxational time-scale. Stars with angular momentum less than a critical value are in the phase space ‘loss cone’ and are tidally destroyed on a dynamical time. Theoretical calculations indicate this rate to be ∼10−4–10−6 yr−1 for Milky Way-like galaxies, and that the peak flux into the loss cone comes from parsec scales. This motivates the assumption of nearly parabolic orbits (Magorrian & Tremaine 1999; Wang & Merritt 2004).
However, there are alternate sources of TDEs distinct from the standard two-body scattering model. Many of these feed stars to the SMBH at lower eccentricity. Our aim here is to quantify through hydrodynamical numerical simulations how the observable properties of tidal disruptions of stars on eccentric orbits deviate from the standard ones. In Section 2, we describe numerical procedures for simulating TDEs. We then analyse results of our simulations in Section 3. In Section 4, we consider non-standard sources of TDEs and whether they can supply stars on low eccentricity orbits to an SMBH. Finally, Section 5 is devoted to summary and discussion of our scenario.
NUMERICAL METHOD
In this section, we describe our procedures for numerically modelling the tidal disruption of stars on bound orbits. The simulations presented below were performed with a three-dimensional (3D) smoothed particle hydrodynamics (SPH) code, which is a particle method that divides the fluid into a set of discrete ‘fluid elements’ (i.e. particles), and is flexible in setting various initial configurations. The code is based on a version originally developed by Benz (1990); Benz et al. (1990) and Bate, Bonnel & Price (1995) and has been extensively used by many papers (e.g. Okazaki et al. 2002; Hayasaki, Mineshige & Sudou 2007).
The SPH equations with the standard cubic-spline kernel are integrated using a second-order Runge–Kutta–Fehlberg integrator with individual time steps for each particle (Bate et al. 1995), which results in enormous computational time savings when a large range of dynamical time-scales are involved. We also use variable smoothing lengths scheme to calculate relevant spatial resolution in our code, but ignore the term proportional to the gradient of the smoothing length.1 In the simulations shown in this paper, we adopt standard SPH artificial viscosity parameters αSPH = 1 or βSPH = 2 and use no switch to limit the artificial viscosity.
We have performed two-stage simulations. We model a star as a polytropic gas sphere in hydrostatic equilibrium. The tidal disruption process is modelled by setting the star in motion through the gravitational field of an SMBH.
Formation of polytropes
Initial configuration
Fig. 2 shows a density map of the polytropic gas sphere at |$4.4\Omega _*^{-1}$|. The left-hand panel shows the column density of polytropic sphere with n = 1.5 over a range of five orders of magnitude. The right-hand panel shows the star-black hole system on the x–y plane, where both axes are normalized by rt and the black hole is put on the origin of the system. The initial position of the star is given by |$\boldsymbol {r}_0=(r_{\rm 0}\cos \phi _0, r_{\rm 0}\sin \phi _0, 0)$|, where |$|\boldsymbol {r}_0|=r_0$| is the radial distance from the black hole and ϕ0 shows the angle between x-axis and |$\boldsymbol {r}_0$|. The star within a small square is zoomed out to the whole left-hand panel. Note that the stellar size on the plot has been scaled up by a factor of a few for visual clarity.
In our simulations, the black hole is represented by a sink particle with the appropriate gravitational mass MBH. All gas particles that fall within a specified accretion radius are accreted by the sink particle. We set the accretion radius of the black hole as equal to the Schwarzshild radius rS = 2GMBH/c2, with c being the speed of light.
Treatment of relativistic precession
Fig. 3 tests how accurately the SPH particles composing the star follow test particle motion in the Schwarzschild metric. Specifically, the figure shows the motion of three particles orbiting the black hole with (e, β) = (0.8, 5), where β = rt/rp is the penetration factor which determines how deeply the star plunges into the black hole potential inside the tidal disruption radius. The central point and dashed circle show the black hole and the tidal disruption radius, respectively. The small crosses, dotted line and solid line represent the geodesics of the Schwarzschild metric, orbits of the test particle under the pseudo-Newtonian potential, and those of a SPH particle in the simulation of Model 2a, respectively.
During the first one and a half orbits, the test particle is in good agreement with the corresponding geodesic. The subsequent deviation from the geodesic is due to the relatively low eccentricity and high β, since equation (10) is tailored to parabolic orbits with lower β (Wegg 2012) (see also Fig. 4). The first three and half orbits of the test particle are in rough agreement with corresponding orbits of the SPH particle. Subsequent large deviations originate from orbital circularization via shock-induced energy dissipation (see Section 3.2 in detail), and although our model becomes unreliable at this point, it is clear that stream crossing is leading to rapid energy dissipation.
Fig. 4 shows the dependence of the error rate on the eccentricity of a test particle moving under the pseudo-Newtonian potential given by equation (10). The error rate is defined by (ΔϕGR − ΔϕPN)/ΔϕGR × 100, where ΔϕGR and ΔϕPN are the precession angles corresponding to equation (2) and (4) of Wegg (2012), respectively, with a given specific angular momentum, specific orbital energy and equation (10). The solid line, dashed line and dotted line show error rates of an orbit with rp = 3rS (β ≈ 7.9), rp = 5rS (β ≈ 4.7) and rp = 10rS (β ≈ 2.4), respectively, where rt ≈ 23.6rS is adopted. The error rate increases as the eccentricity is lower and the penetration factor is higher. From the figure, the error rate is estimated to be roughly 15 per cent for e = 0.8 with β = 5, while it is less than 1 per cent for e ≳ 0.95 with arbitrary β.
Numerical models
We have performed eight simulations of TDEs with different parameters. The common parameters through all of simulations are following: m* = 1 M⊙, r* = 1 R⊙, MBH = 106 M⊙, ϕ0 = −0.4π and γ = 5/3. The total number of SPH particles used in each simulation are 105, and the termination time of each simulation is |$4\Omega _{*}^{-1}$|. Table 1 summarizes each model. Model 1a shows the standard TDE under the Newtonian potential, while Model 1b has the same simulation parameters as Model 1a, except that the star moves under the pseudo-Newtonian potential given by equation (10). Model 1c and 1d have the same parameters as Model 1a and Model1b, respectively, but for β = 5. Model 2a has the same parameters as Model 1d except that the star is on an eccentric orbit, with e = 0.8. Model 2b is the same parameters as Model 2a except that the star moves under the Newtonian potential. Model 3a has the same parameters as Model 2a but for e = 0.98 and β = 1, and Model 3b has the same parameters as Model 3a but for β = 5.
Model . | β = rt/rp . | e* . | a* [rt] . | r0 [rt] . | Remarks . |
---|---|---|---|---|---|
1a | 1 | 1.0 | – | 3 | Newtonian |
1b | 1 | 1.0 | – | 3 | Pseudo-Newtonian |
1c | 5 | 1.0 | – | 3 | Newtonian |
1d | 5 | 1.0 | – | 3 | Pseudo-Newtonian |
2a | 5 | 0.8 | 1.0 | 1.8 | Pseudo-Newtonian |
2b | 5 | 0.8 | 1.0 | 1.8 | Newtonian |
3a | 1 | 0.98 | 50.0 | 3 | Pseudo-Newtonian |
3b | 5 | 0.98 | 10.0 | 3 | Pseudo-Newtonian |
Model . | β = rt/rp . | e* . | a* [rt] . | r0 [rt] . | Remarks . |
---|---|---|---|---|---|
1a | 1 | 1.0 | – | 3 | Newtonian |
1b | 1 | 1.0 | – | 3 | Pseudo-Newtonian |
1c | 5 | 1.0 | – | 3 | Newtonian |
1d | 5 | 1.0 | – | 3 | Pseudo-Newtonian |
2a | 5 | 0.8 | 1.0 | 1.8 | Pseudo-Newtonian |
2b | 5 | 0.8 | 1.0 | 1.8 | Newtonian |
3a | 1 | 0.98 | 50.0 | 3 | Pseudo-Newtonian |
3b | 5 | 0.98 | 10.0 | 3 | Pseudo-Newtonian |
Model . | β = rt/rp . | e* . | a* [rt] . | r0 [rt] . | Remarks . |
---|---|---|---|---|---|
1a | 1 | 1.0 | – | 3 | Newtonian |
1b | 1 | 1.0 | – | 3 | Pseudo-Newtonian |
1c | 5 | 1.0 | – | 3 | Newtonian |
1d | 5 | 1.0 | – | 3 | Pseudo-Newtonian |
2a | 5 | 0.8 | 1.0 | 1.8 | Pseudo-Newtonian |
2b | 5 | 0.8 | 1.0 | 1.8 | Newtonian |
3a | 1 | 0.98 | 50.0 | 3 | Pseudo-Newtonian |
3b | 5 | 0.98 | 10.0 | 3 | Pseudo-Newtonian |
Model . | β = rt/rp . | e* . | a* [rt] . | r0 [rt] . | Remarks . |
---|---|---|---|---|---|
1a | 1 | 1.0 | – | 3 | Newtonian |
1b | 1 | 1.0 | – | 3 | Pseudo-Newtonian |
1c | 5 | 1.0 | – | 3 | Newtonian |
1d | 5 | 1.0 | – | 3 | Pseudo-Newtonian |
2a | 5 | 0.8 | 1.0 | 1.8 | Pseudo-Newtonian |
2b | 5 | 0.8 | 1.0 | 1.8 | Newtonian |
3a | 1 | 0.98 | 50.0 | 3 | Pseudo-Newtonian |
3b | 5 | 0.98 | 10.0 | 3 | Pseudo-Newtonian |
TIDAL DISRUPTION OF STARS ON BOUND ORBITS
We first describe the evolution of a tidally disrupted star for standard TDEs (e = 1). Next, the process of tidal disruption for a star on a fairly eccentric orbit (e = 0.98) is presented. Finally, accretion disc formation in relatively low eccentricity TDEs (e = 0.8) is presented.
In Model 1a, the mass distribution broadens with time as the stellar debris is put on near-ballistic orbits by tidal interaction with the black hole. In panel (a), the central peak formed after t ≳ 1.0 is due to mass congregation, from the self-gravity of the stellar debris (in this scenario, the star is barely disrupted). The energy spread corresponds to Δϵ before and after the tidal disruption. The corresponding mass fallback rates are proportional to t−5/3 except for the solid line in panel (b). This is in good agreement with the literatures (Rees 1988; Evans & Kochanek 1989). The very slight deviation from t to the −5/3 power originates from the convexity around Δϵ due to re-congregation of stellar mass under self-gravity (see Lodato et al. 2009).
Model 1c has the same simulation parameters as Model 1a except for β = 5, while Model 1d has the same simulation parameters as Model 1c except for adopting the pseudo-Newtonian potential. For higher value of β than unity, the tidal forces acting on the star become stronger because the pericentre distance is smaller than the tidal disruption radius. From equation (11), the energy spread of the stellar debris is formally proportional to the square of β, and should become 25 times larger for Model 1c and Model 1d than that of β = 1. In these two models, however, the energy spread is only about 1.4 times wider than that of standard model, showing that the energy spread weakly depends on β. Such a weak dependence of the energy spread on β was found in the context of the binary disruption by Sari, Kobayashi & Rossi (2010). This is also in agreement with recent analytic (Stone, Sari & Loeb 2012) and numerical work (Guillochon & Ramirez-Ruiz 2013). Since the tidal disruption is less marginal as β is higher, the re-congregation of the mass due to self-gravity of the stellar debris is prevented. This leads to the mildly sloped mass distribution, and therefore the peak of the mass fallback rate also smooths. There is no remarkable difference between Model 1c and Model 1d except that the energy spread of Model 1d is slightly wider than that of Model 1c following equation (18).
Eccentric tidal disruption
Figs 9 and 10 show the evolution of differential mass distributions and corresponding mass fallback rates in Model 3a and Model 3b, respectively. The figure formats are the same as Figs 5–8. The mass is not distributed around zero but around −Δϵ in Model 3a, and around −5Δϵ in Model 3b. This is because the specific energy of the initial stellar orbit is originally negative (see equation 19). It is clear from the negative shift of the mass distribution's centre that most of the mass in Model 3a and all of the mass in Model 3b are bounded. As shown in panel (a) of Fig. 9, the resultant energy spread is slightly larger than we analytically expected. This suggests that the critical eccentricity is slightly smaller than the value in equation (20).
From equations (21) and (22), tmin = 0.04 yr and tmax → ∞ for Model 3a, and tmin = 0.0035 yr and tmax = 0.014 yr for Model 3b. These values deviate from simulation results, because they were derived assuming a smaller spread in energy than in our simulations. Furthermore, tmin is shorter than that of the parabolic case. This makes the mass fallback rate about one order of magnitude higher. Notably, the mass accretion is completely finite in Model 3b and its rate is enhanced by about two orders of magnitude, because its time-scale is much shorter than that of the standard model.
Accretion disc formation
Here we describe how an accretion disc forms rapidly around the black hole in Model 2a. The orbital angular momentum of a star passing inside the tidal disruption radius should be conserved before and after the tidal disruption, if there is no mechanism to redistribute the angular momentum. The orbits of the stellar debris can then be circularized by dissipation of orbital energy, primarily due to shocks from orbit crossing, which conserve orbital angular momentum. Debris semimajor axes will thus approach the circularization radius of the initial stellar orbit.
The specific energy and specific angular momentum of a test particle with the orbital parameters of Model 2a, moving under the pseudo-Newtonian potential, are given as ϵtp = −0.497 and ltp = 0.714, respectively, by substituting (rp, ra) = (0.2rt, 1.8rt) with (e, β) = (0.8, 5) into equations (8)–(10). Fig. 11 shows the evolution of an averaged specific energy, specific internal energy, and specific angular momentum of SPH particles in Model 2a. In panel (a), the solid line and the dashed line show the averaged specific energy of SPH particles and the specific energy of a test particle, respectively. The small inside panel shows the evolution of the averaged specific internal energy of SPH particles, which is dissipated and radiated away after t = 2.4Ω−1 due to the shocks from self-intersections of the debris orbits induced by relativistic precession. Since the internal energy contributes negligibly to the total energy budget, the specific energy approximately equals the specific binding energy. On the other hand, panel (b) shows the evolution of the averaged specific angular momentum of SPH particles normalized by ltp. Its marginal fluctuation is a numerical error with a magnitude of 0.15 per cent, indicating that the specific angular momentum is almost conserved. Assuming that the circularized disc has an axisymmetric Keplerian rotation, we can evaluate that the circularization radius rcirc is ≈2.5rp.
Fig. 12 shows sequential snapshots of the surface density of stellar debris (projected on the x–y plane in five orders of magnitude, in a logarithmic scale) for Model 2a. The central small point, dashed circle and dotted line show the black hole, tidal disruption radius, and orbits of a test particle moving under the pseudo-Newtonian potential, respectively. The run time is noted at the top-right corner in units of |$\Omega _*^{-1}$|, while the number of SPH particles are indicated at the bottom-right corner. The star is tidally disrupted before it passes through the first pericentre at t = 0.5. Afterward, the stellar debris expands along the orbit of the test particle as shown in panel (a) and (b). It is stretched towards the second pericentre as shown in panel (c). At the second pericentre, the stretched debris weakly intersects with the remnant of the debris in panel (d) and then reaches the second apocentre while continuing to expand, in panel (e). When it reaches the third pericentre in panel (f), two stretched orbits clearly cross over. The orbital energy is significantly dissipated by the shock from orbit crossing between the two stretched debris streams. This is consistent with the decrease of specific energy in panel (a) of Fig. 11. The stellar debris is therefore rapidly circularized as shown in panels (g) and (h). Finally, in panel (i), a disc like structure is formed around the black hole sufficiently inside the tidal disruption radius. Note that the number of SPH particles are slightly reduced from panel (e) to panel (i), since a very small fraction of the total number of SPH particles is accreted on to the black hole. The trajectory of an SPH particle picked out of this simulation for 0 ≤ t ≲ 4 is drawn in the solid line of Fig. 3. The size of the accretion disc is in rough agreement with rcirc ≈ 2.5rp.
Although the pseudo-Newtonian potential does not accurately model the precession of near-circular orbits, we note that by the time the pseudo-Newtonian potential becomes significantly inaccurate, the debris streams have entered a regime of frequent orbit crossings, guaranteeing further rapid circularization. Furthermore, at a qualitative level the pseudo-Newtonian potential probably overestimates the circularization time-scale, as it underestimates the true rate of relativistic precession for eccentric orbits. However, the pseudo-Newtonian potential still plays a crucial role in orbital circularization processes, because the motion of a test particle under the Newtonian potential take place in a closed path and therefore causes no orbit crossing. In order to test this, we have performed the simulation for Model 2b, which has the same simulation parameters as Model 2a except that the Newtonian potential is adopted. Fig. 13 shows sequential snapshots of the surface density of stellar debris for Model 2b. It has the same format as Fig. 12. We note that there is no significant evidence for orbital crossings during the time-scale of the simulation, since SPH particles orbit around the black hole on fixed eccentric orbits. This is also confirmed by the fact that the percentage change of specific orbital energy from t = 0 to t = 4 for Model 2b is less than 0.4 per cent.
Fig. 14 shows the time-dependence of the number of SPH particles inside the tidal disruption radius Nacc and its first derivative (the mass capture rate). After the tidal disruption of the star, its orbit passes through the first apocentre, going completely outside of the tidal disruption radius. The first peak of Nacc in panel (a) comes when the debris streams pass from the first apocentre to the second apocentre via the second pericentre, and the stretched debris re-enters the tidal disruption radius. Part of it exits once more, but the fractional remaining part is still inside the tidal disruption radius. A sequence of these events can be seen in panels (c)–(e) of Fig. 12. The stretched debris returns again to the tidal disruption radius, moving towards the third pericentre. Afterwards, most of debris circularizes and remains inside of the tidal disruption radius. Panel (b) of Fig. 14 shows the rate of mass being captured inside the tidal disruption radius. The first three peaks are formed as stellar debris passes in and out of the tidal disruption radius, while the final peak shows mass circularizing into the accretion disc around the black hole. The mass accretion rate is more than five orders of magnitude higher than the Eddington rate (see equation 17). This shows that the accretion flow is extremely supercritical in the case of moderately eccentric TDEs.
Figs 15 and 16 show the evolution of differential mass distributions and corresponding mass fallback rates in Model 2a and Model 2b. The figure formats are the same as Figs 5–10, except that solid (black) line shows the differential mass distribution at t = 1.35 and corresponding mass fallback rate. The mass is distributed around −50Δϵ in Model 2b as is expected, whereas the differential mass distribution of Model 2a slightly deviates from −50Δϵ to the positive direction because of corrections of the specific energy from the pseudo-Newtonian potential (see equations 8 and 10). The resultant mass-fallback rate of Model 2a corresponds to that of Model 2b except for the shape of peak of mass fallback rate. This is consistent with our estimate of the mass fallback rate, ∼104 M⊙ yr−1, calculated by taking the ratio of the bound mass, 1 M⊙, to the fallback time Δt for Models 2a and 2b (see also equation 24).
It is noted by definition that the second peak of mass capture rate in panel (b) of Fig. 14 corresponds to the mass fallback rate of Model 2a. It completely deviates from the canonical t−5/3 law, since tmax/tmin for Model 2a is estimated to be ∼3 from equations (21) and (22). This is because the dependence of t−5/3 on the mass fallback rate can be seen only when tmax ≫ tmin.
Since the viscous time-scale measured at rcirc is estimated to be |$t_{{\rm vis}}\sim 3.5\times 10^3(0.1/\alpha _{{\rm SS}})\Omega _*^{-1}$| where αSS is the Shakura–Sunyaev viscosity parameter, the accretion flow is clearly in a non-steady state. Here, we assume that the accretion disc is a geometrically thick : rcirc/H ∼ 1, where H is the disc scale height. Although the fate of the circularized debris is unclear because of much shorter simulation time than the viscous time-scale, the super-Eddington accretion flow will likely drive a powerful outflow (Ohsuga et al. 2005) as it becomes radiation-pressure dominated. This may increase the optical luminosity of the flare by orders of magnitude (Strubbe & Quataert 2009). It is even possible that a radiation-pressure supported envelope could be formed (Loeb & Ulmer 1997).
SOURCES OF ECCENTRIC TDEs
Although the canonical source for TDEs is two-body scattering at parsec scales, a variety of alternate mechanisms exist. In this section we briefly review these, and highlight those which could supply stars at eccentricities e < ecrit to their central SMBH.
Resonant relaxation: at small spatial scales close to an SMBH, the approximately Keplerian potential causes stellar angular momentum to relax in a correlated, coherent way (as opposed to the uncorrelated changes from two-body relaxation). It has been suggested that the rapid changes in orbital eccentricity produced by resonant relaxation (RR) may enhance the TDE rate at small scales (Rauch & Ingalls 1998); however, recent work has indicated that the complicated interplay between RR and general relativistic precession gives rise to a ‘Schwarzschild barrier,’ which quenches RR for stars with small semimajor axis (Merritt et al. 2011). RR is therefore unlikely to produce TDEs with e < ecrit.
Nuclear triaxiality: conversely, the triaxial, non-Keplerian potential of a surrounding star cluster results in fewer constants of motion for stellar orbits; in particular, angular momentum is not conserved and some stars can wander arbitrarily close to the SMBH, eventually being disrupted (Merritt & Poon 2004). However, this effect only arises at large radii where the potential from the star cluster can induce a significant triaxial correction to the SMBH potential, and is therefore unlikely to produce low eccentricity TDEs.
A nuclear stellar disc: there is a rotating stellar disc composed of young massive stars at the centre of the Milky Way (Bartko et al. 2010). These stars orbit with moderate eccentricity around the SMBH SgrA*. The non-RR time-scale drops as a star cluster becomes flattened and disc-like, due to greater encounter rates between stars; simple estimates of the two-body eccentricity relaxation time-scale in the Milky Way's stellar disc (∼5000 M⊙) find values trlx ∼ 4 × 108 yr (Kocsis & Tremaine 2011). Assuming a typical stellar mass of 10 M⊙ gives a disc star TDE rate of ∼10−6 yr−1, implying that if analogous discs are common in other galaxies, they could contribute non-trivially to the total TDE rate. However, in order to produce TDEs with e < ecrit, disc stars would need to be scattered in from ∼ mpc scales, implying that more compact discs than the Milky Way's are needed to produce significantly eccentric TDEs.
Binary SMBHs: a binary SMBH system will, for a period of ∼105 yr to ∼106 yr, see a TDE rate enhancement up to 10−1 yr−1 from a combination of chaotic orbital evolution and Kozai cycles (Ivanov, Polnarev & Saha 2005; Chen et al. 2009). Recent work suggests that chaotic orbits in particular are the dominant contribution to the rate enhancement (Chen et al. 2011; Wegg & Nate Bode 2011). Most of the stellar flux originates from within one order of magnitude of the semimajor axis of the binary; since the inspiral of binary SMBHs stalls when the lower mass secondary reaches a radius inside of which is contained its own mass in stars, a low binary mass ratio q will produce a flux of less eccentric TDEs. In particular, fig. 17 of Chen et al. (2011) indicates that for q = 1/81 and a primary black hole mass of 107 M⊙, some stellar flux into the primary's loss cone originates from spatial scales ≤1015 cm, which implies e < ecrit even for β = 1 events. However, the chaotic orbits which dominate the rate enhancement produce TDEs sampling a wide range of β (Merritt & Poon 2004), so the situation is even more favourable. Binary SMBHs on a eccentric orbit appear to produce even more TDEs from small apocentres. The ‘hardening radius’ (the orbital radius within which the stellar cusp of the primary, larger SMBH contains the mass in stars of the secondary, smaller SMBH) effectively sets the scale at which stellar disruption rates are enhanced; to produce eccentric TDEs requires relatively low mass ratios (<1/50) which could in some cases be characterized as SMBH-intermediate mass black hole (IMBH) inspirals. Because binary SMBHs are expected to provide up to 10 per cent of the total TDE rate (Wegg & Nate Bode 2011), it is important to note that a subset of these events may deviate from canonical ‘parabolic’ light curves.
SMBH recoil: the gravitational wave recoil accompanying an SMBH merger will scramble the orbits of surrounding stars and partially refill the kicked merger remnant's loss cone, briefly increasing the TDE rate. If one considers the bulk of the stellar population surrounding the merging binary black holes, the excavation of a ‘binary loss cone’ results in too few stars at small separation for the resulting burst of TDEs to involve any events with e < ecrit (Stone & Loeb 2011). However, at low binary mass ratios, it is possible for mean motion resonances to migrate stars to small semimajor axes during the binary inspiral (Schnittman 2010; Seto & Muto 2010, 2011); any TDEs produced from this stellar subpopulation would have e < ecrit, and would be described by this paper.
- Binary separation: recent theoretical studies on rates of tidal separation of binary stars by SMBHs suggest that a significant fraction of tidal disruption flares may occur from stars approaching the black hole from subparsec scales (Amaro-Seoane et al. 2012; Bromley et al. 2012). When a binary star passes sufficiently close to the black hole to be tidally separated (without immediate disruption of either component), one star becomes tightly bound to the black hole while the other is flung away at a high speed. The subsequent orbital evolution of the bound star will represent a competition between gravitational radiation and stellar relaxation processes. If gravitational radiation dominates, the orbit will circularize and spiral in, likely leading to a phase of steady mass transfer that is unlike the eccentric and violent disruptions which we have simulated in this paper. On the other hand, if stellar relaxational processes dominate, the star can eventually diffuse into the loss cone with a non-zero eccentricity that is still significantly smaller than that expected for TDEs generated by two-body relaxation at parsec scales. For a bound star with rp > rt and semimajor axis a, living in a stellar cusp with density profile ρ(r) ∝ r−b, the eccentricity dividing gravitational-wave-dominated evolution from relaxation-dominated evolution is (Amaro-Seoane et al. 2012)In practice, requiring egw > ecrit for a < rt/(1 − ecrit) is only possible for steep stellar cusps b ≈ 2 and very low mass SMBHs (MBH ≤ 105.5 M⊙). A further complication is the presence of the Schwarzschild barrier, which may prevent stars with long gravitational wave inspiral time-scales from entering the loss cone. Stars bound to an SMBH from binary separation events can become TDEs with lower eccentricity than in the canonical scenario, but equation (25) implies their eccentricity will still tend to be greater than ecrit. On the other hand, equation (25) was derived assuming the Fokker–Planck diffusion limit for two-body relaxation, and anomalous diffusion from strong scattering events (Bar-Or, Kupi & Alexander 2013) may allow violent tidal disruption of separated binaries even from eccentricities e < egw.(25)\begin{eqnarray} e_{\rm gw} & \approx &1-0.016\times (8\times 10^{-4})^{(2b-3)/5}(3-b)^{2/5} \nonumber \\ && \times \displaystyle\left(\frac{M_{\rm BH}}{10^7\,\mathrm{M}_{{\odot }}}\right)^{(8-b)/5} \left( \frac{a}{1000\, {\rm au}} \right)^{(2b-11)/5}. \end{eqnarray}
Of the possibilities we have considered, two stand out as particularly promising: TDEs generated as binary SMBHs harden and stall, and TDEs generated during the coalescence of binary SMBHs with a low mass ratio by stars that were brought inward via mean motion resonances. The first possibility may account for a significant fraction of the total TDE rate; the second possibility will be much less intrinsically common, but could serve as an electromagnetic counterpart to an eLISA/NGO2-band gravitational wave signal. It is also possible that strong scatterings could lead to eccentric TDEs from the bound stars produced by tidal separation of stellar binaries. Other possible enhancements to theoretical TDE rates, such as perturbations from infalling giant molecular clouds, likely occur at too large of a spatial scale to create eccentric TDEs.
SUMMARY AND DISCUSSION
We have carried out numerical simulations of mass fallback and accretion processes around an SMBH. Specifically, we have examined the tidal disruption of a star on a bound orbit, considering relativistic effects with a pseudo-Newtonian potential. Using a polytropic gas sphere as our initial conditions, we have considered both parabolic orbits (e = 1.0) and eccentric ones (e = 0.8 and 0.98), varying the penetration factor β as well.
We have found that a non-steady, non-axisymmetric accretion disc is formed around the black hole in the case of e = 0.8 and β = 5. The formation of an accretion disc occurs as follows: a segment of the stellar debris returning to pericentre, and a different one exiting pericentre, intersect each other due to relativistic precession. The orbital energy is then dissipated by shocks due to orbit crossing. Since the orbital angular momentum of the stellar debris is conserved before and after the tidal disruption, the debris orbits are rapidly circularized during a few orbit crossings. This shows that the initial size of the accretion disc is only determined by the orbital angular momentum of the initial stellar orbit. In our simulations, the circularization radius is estimated as rcirc ≈ 2.5rp, where rp is the pericentre distance. Furthermore, the expected accretion rate is extremely super-Eddington.
As shown in the final peak of Fig. 14, panel (b), the mass accretion rate exhibits variability during the formation of the accretion disc around the black hole. After some viscous evolution of the accretion disc, however, the variation of the mass accretion rate would settle into that of the mass fallback rate as shown in panel (b) of Fig. 15. Observations of TDE light curves roughly following t−5/3 power laws may suggest that the more common parabolic disruption scenario can circularize efficiently as well, perhaps through general relativity (GR) precession as in the eccentric TDE case. This will need to be demonstrated although it takes much computational time.
The striking difference between moderate-eccentricity simulations with and without relativistic precession highlights the importance of general relativistic effects for debris circularization. Specifically, very little energy dissipation was seen in Model 2b (Newtonian) of Table 1, while its pseudo-Newtonian equivalent, Model 2a, saw rapid accretion disc formation. For reasons of computational cost, we did not follow the longer time-scales required for debris circularization in other models, but in future work we hope to investigate if these results generalize to higher eccentricity and parabolic TDEs. Sufficiently rapid SMBH spin could delay disc formation, as Lense–Thirring torques will precess the orbital planes of individual debris streams and limit or prevent orbit-crossings. If a disc is able to eventually form, its luminosity will be periodically modulated by Lense–Thirring precession (e.g. Stone & Loeb 2012). Even considering this complication, the energy dissipation process is crucial for determining the formation and structure of an accretion disc. In these exploratory simulations, we have adopted the polytropic equation of state instead of solving an energy equation. In a subsequent paper, we will study the detailed disc formation and structure by solving more realistic energy equations with and without a radiative cooling term.
Circularization has rarely been seen in past simulations of main-sequence stars being disrupted by massive black holes, due to the difficulty in following orbits over their full dynamic range in parabolic TDEs. While we have avoided this difficulty by simulating eccentric centre of mass trajectories, there are two other examples of circularizing TDE debris in the literature (Ayal, Livio & Piran 2000; Ramirez-Ruiz & Rosswog 2009). In Ramirez-Ruiz & Rosswog (2009), debris circularization is seen in SPH simulations of stars on parabolic orbits being disrupted by 103 M⊙ IMBHs. In this case, circularization is driven by vertical shocks at pericentre, which serves as a ‘nozzle’ in vertically compressing returning gas. This redistributes orbital elements of the gas and act as a source of ‘effective precession’, causing some gas leaving its second pericentre passage to impact gas returning from its first apocentre passage. Whether or not this nozzle-driven circularization scales up to parabolic orbits around SMBHs is unclear. The SPH simulations of Ayal et al. (2000) also appear to capture nozzle-driven circularization of a tidally disrupted solar-type star on a parabolic orbit around a 106 M⊙ black hole. However, the reliability of these results may be limited by the low particle resolution of these simulations (approximately 5000 SPH particles), and the particle-splitting algorithm employed to address resolution problems, since the vertical compression of the star will be sensitive to spurious changes in stream geometry. If vertical compression at pericentre is under-resolved, the net effect will be unphysical velocity perturbations in the orbital plane (Stone et al. 2012; Guillochon & Ramirez-Ruiz 2013), correspondingly unphysical apsidal precession, and unphysical circularization.
Eccentric TDEs are a subpopulation of all TDEs, but a potentially interesting one. The two most promising means of producing them are in the dynamical friction stage of binary SMBHs, and immediately following gravitational wave-driven black hole coalescence. Both of these are especially interesting subsets of TDEs: the former occurs during a time of greatly enhanced TDE rates, when it may be possible for a single galaxy to produce multiple TDEs in ∼10 yr, while the latter would serve as a delayed electromagnetic counterpart to a low-frequency gravitational wave signal. It is also possible that tidal separation of binary stars could produce TDEs of interestingly low-to-moderate eccentricity, if a large population of MBH ≤ 105.5 M⊙ SMBHs reside in steep stellar cusps. Eccentric TDEs are also of interest because of the relatively short delay time between disruption and disc formation, allowing numerical simulations to bypass costly apocentre passages and directly approach open problems in debris circularization. The pseudo-Newtonian potential proposed by Wegg (2012) is, however, not applicable for the very low eccentricities of fully circularized gas because it was derived with the assumption of small binding energy. Therefore, we hope to employ post-Newtonian corrections in a subsequent paper.
The ongoing tidal disruption of a large gas cloud (G2) falling into the centre of the Milky Way Galaxy is another scenario within the ‘eccentric tidal disruption’ regime. When discovered, the G2 gas cloud was estimated to have a mass ∼ 3 M⊕, and to possess a radius ∼1.8 × 1015 cm (Gillessen et al. 2012). More recent observations have resulted in a finer determination of the cloud's orbital parameters; it will soon reach its pericentre distance of 2.8 × 1015 cm, on an orbit with an eccentricity of 0.966 (Gillessen et al. 2013). This value of eccentricity is far less than the critical eccentricity for the cloud, and consequently all the debris will remain bound to SgrA*. Using equation (24), we find a characteristic mass fallback rate of 6.7 × 10−8 M⊙ yr−1. This rate is likely an underestimate of the true mass fallback rate, which may be enhanced dramatically by ram pressure stripping near pericentre (Anninos et al. 2012; Burkert et al. 2012; Schartmann et al. 2012). Assuming immediate accretion of the returning matter would give a bolometric luminosity of ∼4 × 1038 erg/s−1 with an assumption that the mass-to-energy conversion efficiency is 0.1. This would be several orders of magnitude brighter than the observed X-ray luminosity of 1033 − 34 erg/s−1, detectable with the current X-ray satellites such as Chandra, Suzaku and NuSTAR. However, the rate of matter fallback to pericentre will not necessarily equal the rate of accretion on to the black hole. The viscous time-scale at pericentre is tvisc = (αSSΩp)−1(H/rp)−2, where Ωp is the orbital frequency at the pericentre distance rp, and H is the disc scale height at that radius. For the latest calibration of the orbital pericentre, and the approximate values αSS = 0.1, H/rp = 0.5, tvisc ≈ 7 yr, long enough to significantly reduce the flare luminosity. It is also plausible that the accretion of G2 will proceed with radiative efficiency significantly under 0.1 (Anninos et al. 2012; Mościbrodzka et al. 2012), leading to a more modest increase in luminosity.
An additional source of luminosity during the tidal disruption of G2 could arise from tidal compression of the cloud in the direction perpendicular to its orbital plane. Although most hydrodynamical simulations of the disruption have been two dimensional, some three-dimensional work has found significant (factors of ∼100) levels of vertical compression and density enhancement, and correspondingly large increases in Brγ line luminosity (Saitoh et al. 2012). The lack of change in observed Brγ luminosity (Gillessen et al. 2013) indicates that there should be room for improvement on the simple initial conditions of Saitoh et al. (2012) (i.e. a hydrostatic sphere, initialized in 1995). Although a detailed analysis of G2's tidal deformation is beyond the scope of this paper, we note that strong density enhancements from tidal compression are generic to the tidal disruptions of bodies initially in hydrostatic equilibrium (Stone et al. 2012). If G2's origin is in a tidally disrupting protoplanetary disc (Murray-Clay & Loeb 2012), the gravity of the central star could dominate the central regions and alter the tidal compression process; it is also possible that a shell-like configuration for the cloud (Schartmann et al. 2012) would avoid strong density enhancements until the upper and lower portions of the shell collide near pericentre.
As we have stressed in this paper, mass accretion following tidal disruption is mediated by circularization processes. However, the dramatic GR precession seen in our simulations is unlikely to be relevant for G2, as its pericentre is at a distance of 2200 Schwarzschild radii. Since the Galactic Center environment includes complex features such as an atmosphere of diffuse hot gas, and warped stellar discs composed of young stars (Genzel, Eisenhauer & Gillessen 2010), dynamical interactions could dissipate the orbital energy and angular momentum of the G2 cloud and bring it closer to SgrA*. The most important interaction is likely to be ram pressure stripping shortly after pericentre passage, which has been seen to cause rapid circularization in multiple hydrodynamical simulations. Furthermore, observational and theoretical efforts will be necessary to clarify the future accretion of G2 on to SgrA*, which will differ in many important ways (ram pressure stripping, lack of GR precession, uncertain initial conditions) from the tidal disruption of a star on a comparably eccentric orbit.
For eccentric stellar orbits, there is a critical value of the orbital eccentricity below which all stellar debris is bounded to the black hole. It can be seen from our simulations that the critical eccentricity is slightly lower than expected from our analytical predictions. This might be because of the effects of stellar structure on the tidal disruption. There are three important implications for light curves of eccentric TDEs:
The mass fallback rate will not asymptote to |$\dot{M} \propto t^{-5/3}$| at late times but be finite with cut-off time t = tmax.
The lack of unbound debris will eliminate observational signatures associated with emission lines (e.g. Strubbe & Quataert 2009).
A larger amount of mass will return to pericentre in a much shorter time than in the standard parabolic picture, considerably increasing the ratio of |$\dot{M}/M_{\rm Edd}$| for the subsequent flare.
Even for eccentricities e > ecrit, the centre of the differential mass distribution will shift in a negative direction, providing weaker versions of the above effects. These signatures should be searched for when large samples of TDE candidates from next generation optical transient surveys such as the Large Synoptic Survey Telescope3 and next generation all-sky X-ray surveys such as extended Roentgen Survey with an Imaging Telescope Array4 become available.
The authors thank to the anonymous referee for fruitful comments and suggestions, and thank to Ramesh Narayan for helpful discussions about the recent G2 problem. KH is also grateful to Atsuo. T Okazaki, Takahiro Tanaka and Shin Mineshige for helpful discussions and their continuous encouragement. The numerical simulations reported here were performed using the computer facilities of Yukawa Institute of Theoretical Physics, Kyoto University and Harvard-Smithsonian Center for Astrophysics, Harvard University. This work was supported in part by the Grants-in-Aid of the Ministry of Education, Science, Culture, and Sport and Technology [21540304, 22340045, 22540243, 23540271 KH], Education and Research Promotion Foundation in Kyoto University [KH], and NSF grant AST-0907890 and NASA grants NNX08AL43G and NNA09DB30A [NS, AL].
This term is introduced for calculating the gradient of fluid properties when the smoothing length is varied in space and time, and is important for ensuring energy conservation if the gradient of any physical quantities varies over a shorter length scale than the smoothing length (see Bate 1995 for a review). In our simulations, the specific energy is conserved for all models except for Model 2a. This shows that the term plays no crucial role in our simulations.