Tidal disruptions of white dwarfs from ultraclose encounters
with intermediate mass spinning black holes
Abstract
We present numerical relativity results of tidal disruptions of white dwarfs from ultraclose encounters with a spinning, intermediate mass black hole. These encounters require a full general relativistic treatment of gravity. We show that the disruption process and prompt accretion of the debris strongly depend on the magnitude and orientation of the black hole spin. However, the latetime accretion onto the black hole follows the same decay, , estimated from Newtonian gravity disruption studies. We compute the spectrum of the disk formed from the fallback material using a slim disk model. The disk spectrum peaks in the soft Xrays and sustains Eddington luminosity for yrs after the disruption. For arbitrary black hole spin orientations, the disrupted material is scattered away from the orbital plane by relativistic frame dragging, which often leads to obscuration of the inner fallback disk by the outflowing debris. The disruption events also yield bursts of gravitational radiation with characteristic frequencies of Hz and strain amplitudes of for galactic intermediate mass black holes. The optimistic rate of considered ultraclose disruptions is consistent with no sources found in ROSAT allsky survey. The future missions like WideField Xray Telescope (WFXT) could observe dozens of events.
Subject headings:
accretion – black hole physics – gravitational waves – hydrodynamics – relativity – radiation mechanisms: general – Xrays: bursts1. Introduction
Tidal disruptions of stars by black holes (BHs) are fascinating and violent cosmic events, releasing copious amounts of energy in electromagnetic radiation and, in some cases, also accompanied by potentially detectable gravitational wave emission. The scattered debris from a stellar disruption will yield different radiation patterns depending on the orientation of the orbit. This variety of radiation patterns will provide clues about both the BH and the internal structure of the disrupted star. As the debris from a tidal disruption showers back into the BH to form an accretion disk, for instance, it will radiate close to or above the Eddington luminosity for some duration, with the spectrum peaking at UV/Xray frequencies depending on the BH mass (Rees, 1988; Evans & Kochanek, 1989; Cannizzo et al., 1990; Bogdanović et al., 2004; Strubbe & Quataert, 2009). In addition, the unbound debris will quickly scatter over a large volume, and may be illuminated by the accretion disk, leading to an optical irradiation spectrum dominated by lines (Bogdanović et al., 2004; Sesana et al., 2008; Strubbe & Quataert, 2009; Clausen & Eracleous, 2011).
When a main sequence star is disrupted, the superEddington outflow is by itself hot enough to be a significant source of optical emission resembling supernova radiation (Kasen & RamirezRuiz, 2009; Strubbe & Quataert, 2009). In these circumstances, instead of outflowing debris, a steady optically thick envelope may reprocess the inner disk radiation (Loeb & Ulmer, 1997; Ulmer et al., 1998). As a consequence, neither the inner disk nor the irradiated debris emission is visible. Another tidal disruption effect arises from the tidal compression perpendicular to the orbital plane. The compression leads to the formation of a strong shock that triggers a powerful short outburst (Kobayashi et al., 2004; Guillochon et al., 2009; Rosswog et al., 2009). In some instances, e.g. ultraclose disruptions, the tidal compression could be strong enough to produce a thermonuclear ignition of the star. In such cases, the detonation of the star would be observed as an underluminous supernova (Rosswog et al., 2009).
The variety of proposed radiative signatures often disagree with each other, calling for precise dynamical modeling of tidal disruptions. Early hydrodynamic simulations by Evans & Kochanek (1989) confirmed a flat distribution of debris mass over its energy range, yielding the predicted law of fallback accretion rate (Rees, 1988). Other numerical investigations addressed the tidal compression and related shock formation; some of these studies used smoothed particle hydrodynamics (SPH) techniques (Kobayashi et al., 2004; Rosswog et al., 2009) and others hydrodynamic gridbased codes (Guillochon et al., 2009). The studies by Rosswog et al. (2009) included nuclear reactions with simplified networks. On the other hand, Bogdanović et al. (2004) looked at the formation of spherical envelopes Loeb & Ulmer (1997) using SPH simulations, but the study did not find any conclusive evidence of a spherical structure. At present tidal disruption simulations that capture all the relevant physics are quite challenging. The simulations require handling, in addition to hydrodynamics, effects from general relativity, nuclear reactions, and magnetic fields. Thus, many important questions are still left unanswered. In particular, there are very few simulations that account for general relativistic effects (Laguna et al., 1993a, b; Kobayashi et al., 2004; Bogdanović et al., 2004). The inclusion of magnetic fields may lead to jet formation.
Very likely, the BHs involved in a tidal disruption will have spins misaligned with the orbital angular momentum of the incoming star. Thus the disk formed from the tidal debris will have angular momentum that is also misaligned with the BH spin direction (Rees, 1988). For thin disks, the BardeenPetterson effect will align the BH spin with the angular momentum of the inner region of the disk (Bardeen & Petterson, 1975), but for the thick disks expected from tidal disruptions, the alignment may not happen (Papaloizou & Pringle, 1983; Fragile et al., 2007; Stone & Loeb, 2011). This misalignment could have important consequences on the emission channels from the disruption, e.g. jet formation. The spin of the BH will also produce distinct effects in an ultraclose encounter, when the pericentric radius is comparable to the gravitational radius of the BH. In this strong gravity regime, the details of the disruption will depend not only on the stellar radius, stellar equation of state, and orbital energy, but also on the magnitude and orientation of the spin of the BH. In particular, the properties of the accretion will depend on the BH spin since its magnitude determines the location of the innermost stable circular orbit (ISCO) radius , ranging from for a counterrotating maximally spinning BH to for a corotating maximally spinning BH. The value of spin determines whether the star following the same orbit gets disrupted (Kesden, 2011). Furthermore, a misaligned rotating BH will drag around nearby debris (Bardeen et al., 1972; Shapiro & Teukolsky, 1986), and will effectively push this material away from the orbital disruption plane. The net effect will be the formation of a shelllike disruption debris engulfing the BH, as opposed to the traditional Sshape debris observed in Newtonian gravity simulations (Evans & Kochanek, 1989).
The present work aims at exploring the exciting regime of ultraclose encounters (i.e. few ) of a white dwarf (WD) with an intermediate mass black hole (IMBH). Our main goal is to investigate observational signatures due to strong gravity effects that may shed light on the presence of IMBHs. We consider a carbonoxygen WD with mass and radius (Nauenberg, 1972)
(1)  
with the Chandrasekhar mass and km the Earth mean radius. We fix the mass of the IMBH to . With these choices, . Therefore, the WD and IMBH can be numerically modeled with comparable grid resolutions. The other important scale in the problem is the tidal disruption radius . For our setup (Rees, 1988),
(2) 
gave us ample room to carry out deep penetration encounters.
Currently, there is a large body of evidence for the existence of both solar mass as well as supermassive BHs with masses in the range of – . IMBHs are, however, the missing link between stellar mass and supermassive BHs. Today, only tentative evidence exists (Irwin et al., 2010; Farrell et al., 2009; Davis et al., 2011). On the other hand, WDs are thought to be abundant in spiral galaxies (Evans et al., 1987; Reid, 2005) and globular clusters (Gerssen et al., 2002). Thus the identification of distinct signals in both gravitational waves (GWs) and the electromagnetic spectrum (Gould, 2011) from a disruption event could potentially guide observations that provide evidence for the existence of IMBH as well as insights into the structure of the WD involved.
Our study uses the full machinery of numerical relativity to solve the Einstein equations of general relativity and hydrodynamics. We do not include nuclear reactions. For the ultraclose encounters of our interest, with few , a general relativistic description of gravity is needed. However, given the mass ratio of the bodies involved and encounters not driven by GW emission, one does not need to account for dynamical general relativistic gravity. The present study could have been carried out by doing hydrodynamics on the fixed spacetime background provided by the IMBH. However, this would have implied developing a general relativistic hydro code on a fixed background, a code similar to the SPH code developed by one of us to investigate the tidal disruption of main sequence stars by supermassive BHs (Laguna et al., 1993a, b). We decided instead to take advantage of our numerical relativity code Maya. The code has demonstrated excellent performance in handling fluid flows in the vicinity of BHs in our studies of binary BH mergers in astrophysical environments (Bode et al., 2008, 2009, 2011). The other advantage of using the Maya code is the ability to obtain the GW signal directly from the simulation, without having to take recourse to the possibly inaccurate quadrupole formula or the application of post Newtonian (PN) approximations to spacetimes containing a spinning BH.
The paper is organized as follows. In § 2, we describe how selfconsistent general relativistic initial data is constructed. Results from the dynamics of the disruption events are discussed in § 3. We estimate electromagnetic transient radiation in § 4 using a slim disk model to compute the spectrum during the fallback phase. We find that the sources shine at Eddington luminosity for about yrs, and then fade approximately as . We also find that there is a 50/50 chance that the inner disk will be obscured by outflowing debris for a fully misaligned BH spin. In § 5, we calculate the gravitational wave signal produced during the encounters. In § 6 we estimate the event rates and discuss the observational prospects. We discuss the progress and the limitations of our study in § 7.
2. Initial data and code tests
When using a fully general relativistic code, such as our Maya code, that is based on a 3+1 formulation of the Einstein equations, the initial data for the dynamical spacetime are comprised of the spatial metric and the extrinsic curvature of the initial spacetime hypersurface. All tensor indices are spatial indices and, in this section only, we use units for which . The metric characterizes the gravitational potentials of the IMBH and WD, and the tensor provides the “velocity” of the metric, or the embedding of the spacelike hypersurface, in the spacetime.^{1}^{1}1For a review on numerical relativity, we recommend the textbooks by Alcubierre (2008) and Baumgarte & Shapiro (2010)
Although and are dominated by the IMBH, it is crucial to account for the contributions from the WD in order to correctly include the selfgravity of the WD. The components of and cannot all be freely specifiable. They need to satisfy the following elliptic equations
(3)  
(4) 
where
(5)  
(6)  
(7)  
(8) 
Above, is a conformal factor, the flat spatial metric and its associated Laplace operator. In writing Eqs. (3) and (4), we have made the customary assumption that the initial physical metric is conformally flat and is tracefree. Eqs. (3) and (4) are, respectively, the conformal versions of the socalled Hamiltonian and momentum constraints of general relativity (York, 1979). The sources and are, respectively, the total energy and momentum densities.
Our approach to solve the Hamiltonian and momentum constraints is to specify the “free” data in Eqs. (3) and (4) from the solutions to an isolated, spinning BH, and a boosted WD (Löffler et al., 2006).
Let us consider first the solution for a boosted WD. The sources and in Eqs. (3) and (4) are obtained from the stressenergy tensor for a perfect fluid (see chapter 5 in Baumgarte & Shapiro 2010). They read
(9)  
(10) 
where is the Lorentz factor and the boost velocity of the WD. In (9) and (10), is the restmass density, the pressure, and the specific internal energy density. The values for those quantities are obtained from TolmanOppenheimerVolkoff (TOV) solutions (Oppenheimer & Volkoff, 1939; Tolman, 1939) with a polytropic equation of state . We use appropriate for a nondegenerate gas, which is sufficient to capture the dynamics during the inspiral and disruption phases. During the evolution, we switch to a gammalaw equation of state . We denote the solutions to Eqs. (3) and (4) in the absence of the BH by and .
Now we consider the solution to the constraints for a spinning BH. In the absence of the WD, and . In this case, the momentum constraint (4) can be solved analytically (Bowen & York, 1980). The solution reads
(11) 
with and the threedimensional LeviCivita symbol. Above, is the spin vector of the BH.
With the solutions and at hand, it is clear that the linear superposition is the solution to the momentum constraint (4) in the presence of both a WD and a BH. Next, we solve the Hamiltonian constraint (4) with and given by (9). We use the popular ansatz
(12) 
Thus, equation (3) becomes an equation for . In (12), the first term is the flat space solution to (3) ; the second is the BH solution; and the third term includes both the gravitational potential of the WD as well as the interaction effects with the BH.
To test our initial data method, and in particular to assess the importance of the WDBH interaction effects, we compute initial data for a WD initially at rest; that is, and thus . In addition, we assume a nonspinning BH, i.e. . The Hamiltonian constraint (3), with the help of (12), takes the form
(13) 
The solution to (13) is shown in Figure 1. For color versions of the figures in this paper see the electronic edition of the Journal. Also in Figure 1 is the solution to equation (13) without the BH; that is, with in (12). The difference between and is entirely due to the gravitational interaction of the WD with the BH.
We also tested the ability of our method to produce stable isolated WD models and to conserve rest mass. We found that, over six dynamical times, the central density of the WD did not change more than of its initial value and rest mass was conserved at the level. This is in spite of the adaptive mesh refinement infrastructure we use (Carpet) not being designed to accurately preserve rest mass conservation when interpolating between different refinement levels.
3. Simulation Results
3.1. Tidal Disruption Simulation Parameters
Our study consists of a series of six simulations. As mentioned before, the IMBH has a mass , and the WD has a mass and radius (Hamada & Salpeter, 1961). The WD mass and radius correspond to a central density , , and polytropic constant . All the simulations start with the WD in a parabolic orbit, in the Newtonian gravity sense, a distance along the negative axis away from the IMBH. The orbital angular momentum of the WD is aligned with the axis, as depicted in Figure 2.
Table 1 lists the configurations used in the simulations where is the penetration factor for an orbit with pericentric distance , and is the magnitude of the dimensionless spin parameter of the IMBH. The angles and determine the direction of the BH spin. They are the standard angles in a righthanded spherical coordinate system as depicted in Figure 2. We use the following convention to label the cases we studied. If a case is labeled ”BxSy”, the WD had penetration factor . In addition, the value of y denotes a nonspinning BH (y=0) or the BH’s spin orientation: (up, down, inplane, arbitrary) for y = (u,d,i,a), respectively.
Run  

B6S0  6  0.0  0°  0° 
B6Su  6  0.6  0°  0° 
B6Sd  6  0.6  180°  0° 
B6Si  6  0.6  90°  0° 
B6Sa  6  0.6  63°  90° 
B8Sa  8  0.6  63°  90° 
In all the simulations, we employ eight levels of mesh refinement with radii () centered at the BH location. In addition to these refinement levels, we surround the WD during the predisruption phase with five additional nested boxes of radii (). The resolution on the finest refinement is , the resolution on the level covering the WD is . This mesh refinement setup becomes insufficient when the WD begins to be disrupted. At this point, we turn off the boxes tracking the WD and construct a larger mesh with resolution that includes both the IMBH and the highly distorted WD. Once the expanding debris has cleared the inner region and the density has dropped, we turn off this level again to speed up the simulation.
The mesh refinements affect our ability to conserve rest mass throughout the computational domain, in particular if the meshes are moving, created, or removed. Figure 3 displays the relative gain/loss of rest mass in the computational domain, , where is the total relativistic rest mass and the initial WD rest mass. The computed value of takes into account the mass loss through the BH horizon and the outer boundary of the computational domain. We find tolerable mass conservation, with violations of less than over the course of the simulations. Runs B6Su and B6S0 display the worst mass conservation. As we shall see, these are the cases in which most of the material of the star () escapes the BH, forming an expanding cloud that reaches the coarsest mesh refinements. Notice also that most of the errors in mass conservation accumulate after the disruption. Run B6Sd is the only run that loses mass. This is the case in which the IMBH accretes almost all of the star during the initial passage of the WD.
3.2. Tidal compression
The first stage in a tidal disruption encounter consists of the deformations that the incoming star experiences due to the tidal forces. In the orbital plane of the star, there are two competing tidal deformations. One of them is tidal stretching along the radial direction joining the IMBH with the WD. The other is tidal compression perpendicular to this radial direction. In the direction perpendicular to the orbital plane, there is also tidal compression. This compression, or squeezing, is stronger than that in the orbital plane because in the latter the radial direction changes rapidly as the star reaches periapsis, thus modifying the “contact point” in the star of the inplane compression. We concentrate here on the compression in the direction perpendicular to the orbital plane. This compression could, in principle, generate enough heat to detonate the star (Luminet, 1985; Luminet & Marck, 1985; Luminet & Pichon, 1989) by igniting nuclear reactions at the core of the WD (Rosswog et al., 2009). To investigate whether a nuclear ignition is likely to happen from tidal compression, we plot in Figure 4 the instantaneous maximum temperature of the star as a function of the density at that location. To a good approximation, this density is also the instantaneous maximum density. The instantaneous maximum temperature is computed as the massdensityweighted temperature average of the portion of the star whose numerical cells that have the largest internal energy and make up for of its mass.
We calculate temperatures from the thermal specific energy as (see Lee et al. (2005))
(14) 
with
(15) 
where the second term in equation (15) is the “cold” specific internal energy in the absence of shocks. is the number of electrons per nucleon making up the gas, is the number of ions per nucleon, is the number density of nucleons in the gas, is the proton (nucleon) mass, is the Boltzmann constant, and here is the radiation constant. For simplicity we assume , for an equal mixture of oxygen and carbon atoms. For , which depends on the ionization state of the plasma, we use for and for . That is, we ignore any partial ionization states of the constituent atoms at low temperatures.
In all cases investigated, the maximum temperature and density displayed in the tracks of Figure 4 reach high enough values to trigger nuclear burning in the star (cf. Figure 2 of Dursi & Timmes (2006)). Interestingly, the temperaturedensity tracks in Figure 4 are qualitatively similar to those by Rosswog et al. (2009), who included nuclear networks.
Run  ()  

B6S0  
B6Su  
B6Sd  
B6Si  
B6Sa  
B8Sa 
In Table 2, we list the maximum temperature and compression . Here and are the maximum temperature and density attained over the course of the evolution. We also list the “actual” penetration factor . Recall that the orbital parameters, and in particular , used to set the WD in a parabolic orbit were obtained using Newtonian gravity. Because of relativistic effects, the value of will differ from the Newtonian estimate , and thus . In our simulations, we define as the distance of closest approach of the point within the WD where , or equivalently , is found. Notice that for the B6Su and B6Sd cases, almost doubles and triples when going from a spin aligned with the orbital angular momentum to one that is antialigned. This is because in the antialigned case (B6Sd), the WD penetrates closer, , than in the aligned case (B6Sd), . Figure 5 shows the maximum temperature and maximum density as a function of the actual penetration factor from Table 2. Because of the small range of values covered by , it is not possible to verify the scaling and suggested by Carter & Luminet (1982)
3.3. Outflowing debris and spin
All cases showed that the mass and orientation of the debris after disruption is strongly affected by the BH spin. This is because the spin determines the location of ISCO radius, and it is also responsible for framedragging effects. In particular, framedragging pulls the material out of the orbital plane, thereby changing the orientation of the resulting disk.
We first discuss three disruptions: disruption by a nonspinning IMBH (case B6S0), by an IMBH with spin aligned (case B6Su), and by an IMBH with spin antialigned (case B6Sd) relative to the orbital angular momentum (see Table 2). For the spinning BH cases, the magnitude of the dimensionless spin parameter was kept at . With this spin magnitude, for the case B6Su, and for the case B6Sd. In the nonspinning case, (McClintock et al., 2011). All three cases had penetration factors . From equation (2), , thus a penetration factor of 6 translates to a pericentric distance . We then expect that the WD in the retrograde case B6Sd will pass within the ISCO; in the prograde case B6Su the star will mostly stay outside of the ISCO; and in the case B6S0 the WD will graze by the ISCO. Because in the B6Sd case, almost all of the star is accreted by the BH soon after the disruption. For the cases B6S0 and B6Su, on the other hand, most of the stellar debris escapes direct capture, expanding away from the hole. The debris cloud has a shape resembling a thick circular arch. This archlike cloud eventually closes up and forms a disk around the BH.
In all cases, the leading edge of the tidallydisrupted star wraps around the BH and crashes into the trailing material creating a hot region. For the B6Su, B6Sd and B6S0 cases, this hot region is mostly contained within the equatorial plane and moves outwards from the central region in the form of a spiraling disturbance in the debris disk, eventually appearing on the debris surface. Figure 6 shows snapshots in the orbital plane of the density (top row) and temperature (bottom row) for runs B6Si and B6Su. The snapshots in the left column are for the case B6Si seconds after disruption. The middle and right columns correspond to case B6Su seconds and seconds after disruption, respectively. In the aligned case B6Su, the hot region forms a spiral which moves outwards from the point where the leading edge of the WD intersected with the tail material. At late times this feature has been washed out. In run B6Si the hot region is much smaller; frame dragging by the BH has pulled material out of the plane into a shell surrounding the BH.
For the cases in which the spin of the BH is not aligned with the orbital angular momentum of the binary (i.e. B6Si, B6Sa and B8Sa cases), frame dragging is effectively able to wrap the debris material around the BH. This effect – which is not captured by Newtonian calculations – greatly changes the geometry of the debris, blanketing the central region with a shelllike cloud of gas. We will address the observational consequences of the material surrounding the BH in the next section. The effect of the misalignment is evident from Figure 6 in the remarkable differences of the spiral patterns observed between the left and middle snapshots at seconds after disruption. The left panels are for the B6Si case in which the BH spin direction is perpendicular to the orbital angular momentum. On the other hand, the middle panels are for the B6Su case in which the BH spin is aligned with orbital angular momentum. The observed differences are a consequence of how the debris responds to the orientation of the BH spin.
3.4. Prompt and late accretion
We observe an initial prompt phase of accretion when the star swings by the BH. This prompt accretion phase is followed by an intermediate phase in which material initially deflected by the BH slows down and begins to accrete. At late times, long after the end of our simulations, the characteristic fallback behavior takes over.
Figure 7 shows the flow of matter through the IMBH’s horizon during the first six seconds after the disruption takes place. The instant labeled is the time at which the star enters the tidal radius. The horizontal dashed line labeled “noise” denotes the level of accretion due to the atmosphere used to model the vacuum regions. We observe that the magnitude of the prompt accretion rate decreases the more prograde an orbit is. Both B6S0 and B6Su disrupt at comparable distances from the central BH, but the peak accretion rate differs by almost two orders of magnitude.
As already pointed out in Section 3.3, the BH spin strongly affects what fraction of the star is accreted onto the BH soon ( s) after the star disrupts. Table 2 shows this accreted fraction for each case. In the B6Sd case, the star is completely accreted. On the other hand, the debris in the B6Su and B6Sa cases escapes almost in its entirety during the flyby. Notice also from Table 2 that the prompt accretion is also correlated with the maximum density the WD is able to reach as a consequence of the tidal compression.
The latetime behavior consists of material raining back onto the BH with a characteristic rate of accretion (Rees, 1988; Phinney, 1989). This rate is maintained as long as , with being the specific kinetic energy, is approximately constant (Lodato et al., 2009). We follow Rosswog et al. (2009) and compute the fallback time for each fluid element at a time long after disruption, based on data available at the end of the simulation. At this time, hydrodynamical interactions between the fluid elements are small, and each fluid element moves on an almost geodesic orbit. Figure 8 shows the results of these calculations. Clearly we recover the expected behavior for return times longer than . Notice that earliest times in Figure 8 are comparable to the late time accretion rates found during the simulation in Figure 7.
4. Electromagnetic Signatures
Ultraclose (i.e. few ) tidal disruptions are violent events. The conversion of even a small fraction of gravitational energy into light would result in a powerful electromagnetic signature. In this section we will use the hydrodynamics results from our simulations to make order of magnitude estimates of the tidal disruption flare’s luminosity, its characteristic photon energy, and the fallback material’s electromagnetic signature as it eventually forms a slim accretion disk near the BH.
When the WD descends into the potential well of the BH, with the typical dynamical timescale of s, it acquires energy in excess of . Even if a tiny fraction of this energy gets converted into radiation on this dynamical timescale, it would lead to a huge luminosity which would be noticeable at cosmological distances. Furthermore, since the WD undergoes a substantial acceleration in the deep gravitational potential of the BH, the streams of matter from the disruption will approach speeds , comparable to the speed of light. The streams of tidal debris later collide and transform a portion of their kinetic energy into thermal energy. Being heated by this process in excess of K, the dense hot plasma quickly produces photons (Krolik & Piran, 2011). These photons, however, are effectively trapped by the debris, capping the luminosity around the Eddington luminosity (Eddington, 1926)
(16) 
where is the proton mass and is the Thomson crosssection. Notice that the equilibrium between the gravitational and spherical radiation fields is reached at a luminosity due to composition effects. The system can by far exceed in the presence of a collimated outflow, especially in the first several minutes after the disruption. We elaborate on this tidal disruption phenomenon in a subsequent paper (Shcherbakov et al. 2012, in prep). In the present manuscript we limit ourselves to electromagnetic signatures at late times of months to years.
When the accretion rate falls below the Eddington rate , the accretion becomes radiatively efficient, radiating mostly in Xrays with luminosity . Assuming an efficiency of , this luminosity translates into an accretion rate
(17) 
The time elapsed until the accretion rate drops to can be estimated based on the fallback times reported in § 3. We have demonstrated that the famous fallback law (Rees, 1988) is also present in our ultraclose encounters. From Figure 8, we observe that the accretion rate in all but the B6Sd case drops below at times of yrs. Thus the tidal events we are considering can, in principle, shine at Eddington luminosities for many years. The actual luminosity could be substantially subEddington, though, if outflowing debris produced while the BH swallows matter at superEddington rates obscures the accretion flow. In this section we will first estimate the disk emission within a slim disk model and then consider the added effects of obscuration by the outflowing debris.
4.1. Slim Disk Model
For our study we employ a slim disk model, developed and summarized in Abramowicz et al. (1988); Lasota (1994); Kawaguchi (2003); Sadowski (2009), of a simple onezone verticallyintegrated dynamical model. We compute modified blackbody bremsstrahlung emission from the disk. For our radiation estimates, we assume that the angular momentum of the marginally bound debris accreting at late times, and thus the angular momentum of the disk, are aligned with the orbital angular momentum of the incoming WD. The disk model will then, in principle, be misaligned with the spin of the BH. We presume that, despite the misalignment, the slim disk extends down to ISCO radius computed for the aligned case. In our situation it is challenging to find a more reasonable set of assumptions for a geometrically thick slim disk. The ultimate answers require performing full numerical relativity simulations of the fallback regime, a task requiring modeling dynamical timescales times longer than those addressed in the present study.
We take our heavy white dwarf to consist of approximately carbon and oxygen as suggested by stellar evolution computations by Mazzitelli & Dantona (1986). The same heavy white dwarf composition was assumed by both Rosswog et al. (2009) and Krolik & Piran (2011). Thus there is nuclei and electrons per one nucleon with mass . In terms of nucleon density, the ion and electron densities are
(18) 
respectively. The corresponding nonrelativistic heat capacity of neutral and fullyionized gas per single nucleon of temperature are
(19) 
respectively. The energy of full ionization is for the assumed composition, where is the charge of the nucleus. A CNO WD is thus fully ionized at . The molecular mass of the fullyionized plasma is , and its density is .
A fully selfconsistent model of a slim disk requires an estimate of the disk scaleheight . For a disk in the superEddington regime, (Abramowicz et al., 1988; Kawaguchi, 2003; Strubbe & Quataert, 2009). In our simulations, circularization of debris happens within a distance from the BH. We therefore consider a truncated disk located between the ISCO radius and the maximum radius . The viscosity parameter is taken to be . We solve the equations for the central temperature as detailed below, taking into account the effects of viscous heating, emission, and the finite heat capacity of matter. For simplicity, we follow Shakura & Sunyaev (1973) and set the viscous stress at the ISCO to zero. The energy release per unit disk area is then
(20) 
The accretion rate is expressed as
(21) 
where is the onezone density and is the inflow velocity.
We adopt a beta disk model where the viscous stress is given by with being the gas pressure (Sakimoto & Coroniti, 1981). A beta disk model is one of three options. In the standard Shakura & Sunyaev (1973) model, , whereas in the “mean field” models (see e.g. Done & Davis 2008 for a review). There is no definitive answer yet about which model better describes stress in a radiationdominated accretion disk. The beta disk model assumption leads to a radial velocity of plasma
(22) 
where is the isothermal sound velocity. The expression (22) for radial velocity, though not fully selfconsistent, helps to avoid the infinite increase of the density near the ISCO inherent to early slim disk models (Abramowicz et al., 1988).
The accretion disk emits predominantly bremsstrahlung radiation. Since the scattering crosssection in the disk is much higher than the absorption crosssection, the emission spectrum is described by a modified black body spectrum (Rybicki & Lightman, 1979; Czerny & Elvis, 1987; Kawaguchi, 2003). The photon production rate via bremsstrahlung is (Rybicki & Lightman, 1979; Katz et al., 2010)
(23) 
where is the fine structure constant and is the Thomson crosssection. With the high densities, the production rate is high so the radiation quickly dominates internal energy and pressure.
Let us estimate the equilibrium temperature of this radiationdominated flow. The typical nucleon density in the simulation varies from near the BH, to a minimum of The typical virial temperature is K. As photons are produced, the equilibrium implies
(24) 
leading to photon temperatures K. The timescale for reaching equilibrium can be estimated as
(25) 
which is significantly shorter than the dynamical timescale s.
We work in the diffusion approximation of high optical depth to find the surface temperature. The corresponding scattering opacity, , is much larger than the bremsstrahlung absorption opacity . The electron scattering opacity is
(26) 
whereas the absorption opacity is
(27)  
where (Shapiro & Teukolsky, 1986). The total emitted power through one side of the disk with surface temperature is then
(28) 
where, in a onezone model, the central density is taken as a proxy for the surface density. Here denotes mean opacities. We use here the Rosseland mean opacity (Shapiro & Teukolsky, 1986), .
The energy flux delivered to the surface is
(29) 
To sustain the radiation flow through the surface, this energy flux is equal to the radiated flux. Thus from (28) and (29) we obtain that . In a slim disk, cooling does not balance heating, . The residual energy is stored within the accreting material. The energy balance can then be written as
(30) 
where the residual energy per unit surface area is given as
(31) 
where . We solve the equations (21) and (22) for the central density . The central temperature is obtained from equation 30 with the help of (20), (29) and (31). Given the central density and temperature, we calculate the surface temperature from (28) and (29) as well as the corresponding spectrum.
Figure 9 shows the dependence of surface temperature on radius for the accretion rates inferred from the output of the B6Su simulation. Later times correspond to lower accretion rates. The temperature is, in general, higher at higher accretion rates, i.e. earlier times. This is because the energy release per unit time goes up with accretion rate while the radiating surface area stays constant. At very low , our slim disk accretion model resembles the standard thin disk model (Shakura & Sunyaev, 1973) so at the ISCO. In contrast, at the ISCO saturates at the highest accretion rates. With little energy release, , the energy density is proportional to the accretion rate . In turn, the scattering opacity is also proportional to the accretion rate . The energy density of radiation is proportional to . In addition there is a factor of in both the numerator and the denominator of the expression (29) for the surface flux. In the end, in the radiatively inefficient regime.
The corresponding disk spectra are shown in Figure 10. From top to bottom, as a function of energy are shown at months, then and years, respectively. Like , the hard tail of the spectrum also saturates at the highest accretion rates. On the other hand, the lowenergy tail does not saturate at . This is because the radiation time is still shorter then or comparable to the advection timescale at large distances from the BH, where the lowenergy tail is emitted. As the accretion rate drops with time, the luminosity (the area under the curve) also drops.
The integrated luminosity is shown in Figure 11 for the case B6Su. At late times, i.e. years, the luminosity is proportional to the accretion rate and decreases as . Specifically, we estimate a luminosity for our disk model. Notice that the estimated luminosity in a fully general relativistic model is for a thin disk without radiative transfer effects (Bardeen et al., 1972). At early times years the luminosity saturates at about the Eddington value . The luminosity at times months is not reliable, since the slim disk model may not be an adequate representation of accretion at rates . For instance, when the outflow from the inner disk is present, the luminosity may further increase (Owocki & Gayley, 1997; Shaviv, 2001; Begelman, 2001; Dotan & Shaviv, 2011). Such a computation is, however, beyond the scope of the paper.
We consider now the radiative properties of disks from the other cases with different penetration parameters , spins magnitudes , and spin orientations. Figure 12 shows the disk spectra at yr after the disruption for all the cases studied. In all cases, we assumed a slim disk with inner boundary at and outer boundary at .
The spectra in Figure 12 are different because of two factors. The first contributing factor is the accretion rate , which we parameterize by the fallback time. More debris falls back onto the BH in the B6Sa case than in the other cases (see Figure 8). Therefore the spectrum for the B6Sa case is systematically brighter. Similarly, the lower accretion rate for the B6Si case makes its spectrum systematically weaker than the others. We found that one can achieve coincidence of spectra among all cases with a time shift. That is, at times when the accretion rates are identical for different simulations, the spectra are also identical. The spectrum for the B6Sa case at about months is identical to the spectrum for the B6Su case at months and to the spectrum from the B6Si case at about months. Also, more matter is swallowed by the BH initially in the B8Sa case when compared with the B6Sa case, but at later times, less matter that will be eventually accreted is scattered. Thus the fall back time is lower for the B8Sa case, resulting in a less luminous spectrum. The spectra also vary because of changes in the ISCO radius. For instance, the ISCO is further from the BH in the B6S0 case with spin . Since little energy release is taking place inside , its spectrum is significantly softer.
4.2. Obscuration
The edgeon view of a slim disk is selfobscured. The spectrum of a selfobscured disk is not easy to determine. We can approximate the spectrum by assuming that the thick outer edge at emits modified black body radiation the way the top of the slim disk does. The corresponding edgeon spectrum is included in Figure 12 as the black dotted line. It is softer, since only the material far from the BH contributes. The surface area of the thick outer edge is large, though, resulting in a raised lowenergy tail.
Besides exploring the properties of a slim disk, we investigated whether the emission from the disk is visible through the outflowing debris surrounding the disk. Figure 13 shows the directions in which the debris is scattered in the “sky” from the BH point of view. We considered only the material outside a sphere of radius , corresponding to the boundary between the bound and the unbound matter at time s after the disruption. The unbound matter travels essentially along straight lines at s so Figure 13 adequately represents the angular distribution of matter at late times. White color in Figure 13 indicates there is no matter traveling in that direction whereas all darker colors indicate the presence of gas.
The optical depth is very high at directions that are not denoted by white. The inner disk would be completely obscured along those directions. The left panel in Figure 13 shows the distribution of the scattered debris for the B6Si case. Notice that although the debris is scattered all over the sky, there are some preferential directions manifest as a sinusoid in the plane. On the other hand, in the case B6Su shown on the right panel, the debris is scattered preferentially in the disruption plane. Thus the orientation of BH spin determines the direction of debris outflow for ultraclose encounters.
We consider now the situation in which we view a WDBH system with an unknown orientation. Recall that we have assumed that the disk at late times is fully aligned with the initial disruption plane. As discussed before, the disk can, in principle, be viewed both edgeon or faceon, although the edgeon disk is selfobscured and has a peculiar spectrum (see Figure 12). The probability the disk is obscured by the outflowing debris is proportional to the solid angle subtended by the debris. This solid angle corresponds to the normalized surface area on the plane (see Figure 13). There may in principle be four different types of obscuration: i) edgeon selfobscured disk also obscured by outflowing debris, ii) edgeon selfobscured disk unobscured by outflowing debris, iii) faceon disk obscured by outflowing debris, and iv) faceon completely unobscured disk.
When the height of the photosphere equals the disk height (inclination ), the probability to encounter a selfobscured edgeon disk is . The obscuration probability by outflowing debris is depicted in Figure 14 for both edgeon and faceon disks. The panels show the probability as a function of debris location at s after the disruption. At that time, the boundary between the bound and the unbound matter is located at . Thus, this obscuration probability is also the obscuration probability for days to years later. The obscuration probability of an edgeon disk can be inferred from the left panel. Since most debris is scattered in the disruption plane, this probability is relatively high. The values reach , practically independent of spin value, spin orientation or penetration factor. The obscuration probability of a faceon disk can be inferred from the right panel. Only the fully misaligned case B6Si scatters significant amount of debris perpendicular to the disruption plane, which leads to relatively high for this simulation. Partially and fully aligned simulations all scatter matter within from the disruption plane, and thus cannot obscure the faceon disk, so that . The total obscuration probability is
(32) 
The values are for the B6Si case and for the B6Su, B6Sa, B8Sa and B6S0 cases. The misaligned simulation with yields a marginally higher compared to a similar simulation with for which . As expected, deeper penetration results in wider scattering of debris. However, surprisingly, partially aligned and fully aligned simulations with and yield almost the same obscuration probabilities.
We have previously stated that if there is debris along the line of sight, then the inner disk is obscured. However, the optical depth of debris decays as the debris expands and eventually becomes transparent. Let us estimate the time when the optical depth of Xray absorption equals unity . The actual time will of course depend on the line of sight. The debris consists mostly of neutral carbon and oxygen with high absorption crosssection within keV photon energy range (Henke et al., 1993). The typical outflow velocity of the debris is a substantial fraction of the speed of light . The transparency time is then
(33) 
where is the fraction of outbound material. From our simulations, we found that () for the B6Si (B6Su) case. With our obscuration probability estimate ( ) for the B6Si (B6Su) case, we find that yr ( yr). The transparency time given by equation (33) is an underestimate, since marginally unbound matter moves at speeds lower than the average .
5. Gravitational wave emission
In addition to producing the electromagnetic signatures discussed in previous sections, tidal disruptions events are potential sources of GW radiation. A disruption event generates a burst of GWs lasting (Kobayashi et al., 2004; Rosswog et al., 2009). The GW burst produced will have a characteristic frequency and amplitude , with the distance to the BH. For the case of a WD disrupted by an IMBH (Rosswog et al., 2009),
(34)  
(35)  
These estimates place the GWs from these events at the low frequency edge of the advanced LIGO design sensitivity range. At this edge, the equivalent strain noise is (LIGO Scientific Collaboration, 2009).
In Figure 15, we compare the GW strain extracted from the spacetime (solid red) to that obtained using the quadrupole formula Rosswog et al. (2009) applied to the WD’s center of mass trajectory (dashed blue) for B6Sa. While the qualitative behavior of the two polarizations, (top) and (middle), are similar, the fully nonlinear GW strain differs in both amplitude and phase.
In Figure 16, we show the GW strain as a function of time for disruption events at a distance of 20 kpc. Depicted are the two polarizations of GWs: (solid red line) and (dashed green line) for cases B6Sa (top) and B8Sa (bottom). Time denotes approximately the time when the disruption takes place. Evident is the burstlike nature of the GWs. Notice that both cases exhibit comparable burst duration and strain amplitude. This is expected since the only difference between B6Sa and B8Sa is their value of . Although we have stated that case B6Sa has and B8Sa a , from Table 2, we have that the actual penetration factor for both cases is . Interestingly, the other cases also showed comparable burst durations, s, suggesting a weak dependence with the orientation of the spin and penetration factor.
Figure 17 shows the power spectrum for all the cases investigated where we assume that the events took place at a distance of 20 kpc. All cases indicate a characteristic frequency of Hz. This is about an order of magnitude smaller that the estimated frequency using (34) with , the average real penetration factor (see Table 2). This should not be too surprising since rough estimates using (34) and (35) do not take into account, among other things, spin effects. Similarly, Figure 17 shows characteristic strain amplitudes of , an order of magnitude higher than the estimate using (35) with .
At low frequencies, Hz, we note that the B6Sa and B6Su cases in particular have higher strain amplitude than the rest. Interestingly, these are the same cases which, as seen in Fig. 7, exhibit significantly lower prompt accretion in the first second after disruption and higher accretion in the s after disruption. The third case with higher prompt accretion in the s interval, B8Sa, however, exhibits the lowest lowfrequency tail. This suggests that relativistic effects create an entangled parameter space involving spin and penetration factor during the disruption phase. Given that the actual penetration factors in all the cases were comparable, , further study is required to explore this parameter space.
Synergistic observations of electromagnetic emission from the prompt accretion discussed in Section 4, together with GW detections, may provide ways to measuring the spin of the BH since in both instances the emission seems to depend, although not strongly, on the spin orientation. We are currently expanding the parameter space covered by our simulations to investigate this dependence in more detail. Multimessenger observations will be, of course, only possible if the prompt electromagnetic emission is not obscured by the debris surrounding the BH. If such obscuration prevents observing prompt accretion, the detection of a GW burst would, in that case, still serve as an identifying precursor to the observations of the late fallback electromagnetic emission.
6. Rates and Observational Prospects
6.1. Tidal Disruption Rates
The disruptions of stars by BHs are rare events. More studied are tidal disruptions of stars by SMBHs in galactic centers. Their rate per galaxy was predicted theoretically by (Magorrian & Tremaine, 1999; Wang & Merritt, 2004) and found consistent with observations by ROSAT (Donley et al., 2002), XMMNewton (Esquej et al., 2007), and GALEX (Gezari, 2010). The volume density of SMBHs is about (Tundo et al., 2007)
(36) 
which leads to tidal disruption rate per unit volume
(37) 
The rate of disruption of stars in globular cluster (GC) is much less certain. Baumgardt et al. (2004) predicts that, for a GC with BH, the optimistic disruption rate of stars is per GC. According to McLaughlin (1999), the mean GC mass is , which is larger than clusters considered in Baumgardt et al. (2004). For a largest cluster simulated by Baumgardt et al. (2004) with stars they find of disrupted stars to be WDs. Therefore the optimistic rate of WDIMBH tidal disruptions is
(38) 
per cluster. The number of GCs changes dramatically from one galaxy to another (Harris, 1991). However, GCs have about constant formation efficiency (McLaughlin, 1999), which is the ratio of GC mass to the initial mass of gas in a galaxy. We assume half of baryonic matter in galaxies, take the Hubble constant and current baryonic content (Komatsu et al., 2011) and find the number density of GCs as
(39) 
This density overpredicts the number of GCs in the Milky Way, but is consistent with the overpopulated with GCs elliptical galaxies such as M87. It is not known, what percentage of GCs contains an IMBH. For an optimistic estimates, we take one IMBH per cluster. Then, the disruption rate per unit volume is
(40) 
This rate is even higher than , but the WDIMBH disruptions are much less luminous, so they are harder to find. Also, the rates of tidal disruptions in GCs are quite uncertain (Baumgardt et al., 2004), so even disproving the rate given by equation (40) does not necessarily disprove the hypothesis that every GC has an IMBH.
Out of all tidal disruptions, we focused on ultraclose ones with the ratio of tidal radius to pericenter distance about . These constitute a fraction of all disruptions. For a realistic triaxial potential of a GC, chaotic feeding would prevail (Merritt & Poon, 2004). Chaotic feeding leads to tidal disruption probability proportional to for objects coming within pericenter distance , as a clear consequence of gravitational focusing. Therefore, tidal disruptions with constitute of all events. The correspondent optimistic rate is
(41) 
6.2. Observations of Xray Flares
Considered ultraclose tidal disruptions manifest with a yearlong Xray flare at an Eddington luminosity off the center of the galaxy. Tidal disruptions by SMBHs rather produce an optical flare (Strubbe & Quataert, 2009) in the center of the galaxy. Let us explore if the optimistic high rates of tidal disruptions are consistent with past observations and if the future observations have a chance of catching such an event.
ROSAT survey probed all sky down to flux within keV band (Voges et al., 2000). The peak of Xray tidal disruption spectrum at yr (see Figure 12) is at about keV, so that more than half of photons are captured. Taking luminosity of B6Su simulation at yr we obtain the maximum distance of
(42) 
Therefore, the expected number of observed WDIMBH disruptions is , which is consistent with no events. With Chandra, photons can be obtained over ks from the source Mpc away and over Ms from the source Gpc away, according to the Portable, Interactive MultiMission Simulator (PIMMS; Mukai 1993). However, owing to a small field of view of the satellite, the careful choice of field is required to observe at least one event.
In turn, the future mission like Wide Field Xray Telescope (WFXT) (Conconi et al., 2010) can detect a substantial number of events. It is much more sensitive than ROSAT and observes in the same waveband. The flux limits of , , and for deep, medium and wide surveys with correspondent solid angles , , and yield maximum event distances Gpc, Mpc, and Mpc. This translates into , , and WDIMBH disruptions or times fewer ultraclose disruptions.
Tidal disruptions can also produce powerful prompt superEddington Xray flare with the duration of minutes in addition to the Eddingtonlimited Xray flux on the timescale of yr. The lightcurve, spectrum and detection prospects will be discussed in the next paper Shcherbakov et al., (2012, in prep.).
7. Discussion and Conclusions
Our work addresses disruptions of WDs by IMBHs. This is the first fully generalrelativistic study that includes both strong gravity (the metric of the spinning BH) and dynamical gravity effects (GWs). We focused on ultraclose disruptions, where the periapsis radius is comparable to the BH gravitational radius. Our study presents for the first time results on the influence of the BH spin on the disruption.
Our study made some simplifying assumptions. Instead of employing a degenerate equation of state for the WD, we used an ideal gas law. This is, however, expected to have a very minor effect on dynamics. If the entropy is approximately constant, the evolution of a degeneracy pressure supported star is identical to that with an ideal gas pressure. The choice of pressure prescription does not matter until shocks and turbulence take over. Similarly, the absence of explicit radiation pressure does not influence the initial stages of the disruption. Cold outflowing debris may not have significant radiation pressure support, and only the fallback disk will be radiation pressure dominated. Radiation emitted from the inner disk may halt the infall of the marginally bound debris. However, the disk is selfobscured near the disruption plane where most of debris is scattered, and the photon mean free path is quite small. The photons cannot thus effectively transfer energy from the inner flow to the outer flow to unbind the debris. Should radiation pressure be included, the same internal energy contribute a radiation pressure smaller than the gas pressure, and thus the total pressure in the medium would be lower.
Our study does not include nuclear reactions, believed to be important for deep WD disruptions (Rosswog et al., 2009). Temperaturedensity tracks of the matter in our simulations as depicted in Figure 4 are comparable to those found in Figure 14 of Rosswog et al. (2009). The tracks suggest that nuclear burning of carbon and oxygen should occur in the cases we studied. The additional energy gained in C and O burning into Fe is MeV per nucleon (Dursi & Timmes, 2006). This translates to a gain of energy of . On the other hand, the descent into the gravitational potential of the BH yields energies of for a periapsis radius of . The dispersion of debris energies after ultraclose disruptions is also large. Thus, the extra energy from the burning will not likely have a dramatic effect on the debris. The optical lightcurve may indeed change after an initially cold white dwarf is heated by nuclear reactions, but the influence on the fallback disk is expected to be small.
A realistic fallback disk simulation should incorporate magnetic fields. Magnetic fields give rise to large viscosity and allow the disk to selfconsistently accrete. The assumed parameter will be a proxy for viscosity until magnetized simulations can be run for a sufficiently long time. The presence of magnetic fields may readily lead to an outflow (Blandford & Payne, 1982), which may boost luminosity far above the Eddington limit or radiation pressure itself may drive an outflow under certain circumstances (Shaviv, 2001). The duration of our simulations ( s) was too limited to directly model the fallback. The temporal dependence of accretion rate depicted in Figure 8 is a ballistic guess, based on a confirmed flat distribution of mass over energy. Any outflow or halted inflow at late times may result in deviations of the fallback law from .
We make a distinction in the paper between the faceon and edgeon disks. In reality, only minor differences between these orientations may exist. The disk with extreme accretion rate may drive a wind or have a corona, which might completely obscure the very inner portions of the disk. In such a case, the spectrum would be similar to an edgeon spectrum regardless of orientation. The precise quantitative predictions of spectra will only be possible with detailed numerical simulations of the fallback accretion.
Our study should be viewed as a step towards realistic tidal disruptions of WDs by IMBH. Some of the conclusions may hold up despite our limited understanding of matter fallback and superEddington accretion flows. Our main conclusions are:

For a nonspinning BH and a BH with a spin aligned or misaligned with the orbital angular momentum, the debris after disruption forms a thick accretion disk.

For misaligned spins, framedragging effects scatter the debris around the BH and will often obscure the inner region from observation.

There is a qualitatively different behavior before and after a timescale of yr, when the accretion rate approaches Eddington accretion rate. The accretion flow luminosity stays around before yr and starts dropping as afterwards.

The spectrum peaks at soft Xrays. The spectra are similar to thin disk spectra at low accretion rates ,

Selfobscuration and obscuration by debris are present in the system. Selfobscuration leads to softer spectrum while obscuration by debris may make a fallback disk invisible during most of the active accretion period.

The GW signal depends weakly on the orientation of the spin. The GW burst will be challenging to be detected for extragalactic sources.
8. Acknowledgements
The results reported in this paper were obtained with the the Maya code of the numerical relativity group at Georgia Tech. The Maya code solves the Einstein equations of general relativity using the BaumgarteShapiroNakamuraOoharaKojima (BSSNOK) formulation, and BHs are modeled following the moving puncture method (Campanelli et al., 2006; Baker et al., 2006). The source code is partially generated using Kranc (Husa et al., 2006). The Maya code works under Cactus, with adaptive mesh refinement provided by Carpet (Schnetter et al., 2004). Hydrodynamics in the Maya code is based on the public version of the Whisky code developed by the European Union Network on Sources of Gravitational Radiation (Baiotti et al., 2005). We construct initial data using the 2Punctures code of Ansorg et al. (2004), and use AHFinderDirect (Thornburg, 2004) to locate the apparent horizon. Both codes are part of the Einstein Toolkit (EinsteinToolkit, 2011). Details of the Maya code can be found in Vaishnav et al. (2007); Hinder et al. (2008); Herrmann et al. (2007b, a); Bode et al. (2008, 2009, 2011).
The authors are thankful to Ramesh Narayan for fruitful discussions. RVS is supported by NASA Hubble Fellowship grant HSTHF51298.01. RH gratefully acknowledges support by the Natural Sciences and Engineering Council of Canada. Work supported by NSF grants 0653443, 0855892, 0914553, 0941417, 0903973, 0955825. Computations at Teragrid TGMCA08X009 and the Georgia Tech FoRCE cluster.
References
 Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, Ap. J., 332, 646
 Alcubierre (2008) Alcubierre, M. 2008, Introduction to 3+1 numerical relativity (Oxford: Oxford University Press)
 Ansorg et al. (2004) Ansorg, M., Brügmann, B., & Tichy, W. 2004, Phys. Rev. D, 70, 064011
 Baiotti et al. (2005) Baiotti, L., Hawke, I., Montero, P. J., Löffler, F., Stergioulas, N., Rezzolla, L., Font, J. A., & Seidel, E. 2005, Phys. Rev., D71, 024035
 Baker et al. (2006) Baker, J. G., Centrella, J., Choi, D.I., Koppitz, M., & van Meter, J. 2006, Phys. Rev. Lett., 96, 111102
 Bardeen & Petterson (1975) Bardeen, J. M., & Petterson, J. A. 1975, Astrophys.J., 195, L65
 Bardeen et al. (1972) Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, Ap. J., 178, 347
 Baumgardt et al. (2004) Baumgardt, H., Makino, J., & Ebisuzaki, T. 2004, Ap. J., 613, 1143
 Baumgarte & Shapiro (2010) Baumgarte, T. W., & Shapiro, S. L. 2010, Numerical Relativity: Solving Einstein’s Equations on the Computer (Cambridge University Press)
 Begelman (2001) Begelman, M. C. 2001, Ap. J., 643, 1065
 Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
 Bode et al. (2011) Bode, T., Bogdanovic, T., Haas, R., Healy, J., Laguna, P., & Shoemaker, D. 2011, ArXiv eprints
 Bode et al. (2008) Bode, T., Shoemaker, D., Herrmann, F., & Hinder, I. 2008, Phys. Rev. D, 77, 044027
 Bode et al. (2009) Bode, T., et al. 2009, Phys. Rev., D80, 024008
 Bogdanović et al. (2004) Bogdanović, T., Eracleous, M., Mahadevan, S., Sigurdsson, S., & Laguna, P. 2004, Ap. J., 610, 707
 Bowen & York (1980) Bowen, J. M., & York, Jr., J. W. 1980, Phys. Rev. D, 21, 2047
 Campanelli et al. (2006) Campanelli, M., Lousto, C. O., Marronetti, P., & Zlochower, Y. 2006, Phys. Rev. Lett., 96, 111101
 Cannizzo et al. (1990) Cannizzo, J. K., Lee, H. M., & Goodman, J. 1990, Ap. J., 351, 38
 Carter & Luminet (1982) Carter, B., & Luminet, J. P. 1982, Nature, 296, 211
 Clausen & Eracleous (2011) Clausen, D., & Eracleous, M. 2011, Ap. J., 726, 34
 Conconi et al. (2010) Conconi, P., Campana, S., Tagliaferri, G., Pareschi, G., Citterio, O., Cotroneo, V., Proserpio, L., & Civitani, M. 2010, MNRAS, 405, 877
 Czerny & Elvis (1987) Czerny, B., & Elvis, M. 1987, Ap. J., 321, 305
 Davis et al. (2011) Davis, S. W., Narayan, R., Zhu, Y., Barret, D., Farrell, S. A., Godet, O., Servillat, M., & Webb, N. A. 2011, ApJ, 734, 111
 Done & Davis (2008) Done, C., & Davis, S. W. 2008, Ap. J., 683, 389
 Donley et al. (2002) Donley, J. L., Brandt, W. N., Eracleous, M., & Boller, T. 2002, AJ, 124, 1308
 Dotan & Shaviv (2011) Dotan, C., & Shaviv, N. J. 2011, MNRAS, 413, 1623
 Dursi & Timmes (2006) Dursi, L. J., & Timmes, F. X. 2006, Astrophys. J., 641, 1071
 Eddington (1926) Eddington, A. S. 1926, The Internal Constitution of the Stars (Cambridge: Cambridge University Press)
 EinsteinToolkit (2011) EinsteinToolkit. 2011, Einstein Toolkit: Open software for relativistic astrophysics
 Esquej et al. (2007) Esquej, P., Saxton, R. D., Freyberg, M. J., Read, A. M., Altieri, B., SanchezPortal, M., & Hasinger, G. 2007, A&A, 462, L49
 Evans et al. (1987) Evans, C. R., Iben, I., & Smarr, L. 1987, Astrophys.J., 323, 129
 Evans & Kochanek (1989) Evans, C. R., & Kochanek, C. S. 1989, Ap. J., 346, L13
 Farrell et al. (2009) Farrell, S., Webb, N., Barret, D., Godet, O., & Rodrigues, J. 2009, Nature, 460, 73
 Fragile et al. (2007) Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, Ap. J., 668, 417
 Gerssen et al. (2002) Gerssen, J., van der Marel, R. P., Gebhardt, K., Guhathakurta, P., Peterson, R. C., & Pryor, C. 2002, AJ, 124, 3270
 Gezari (2010) Gezari, S. 2010, in IAU Symposium, Vol. 267, IAU Symposium, 319–324
 Gould (2011) Gould, A. 2011, ApJ, 729, L23+
 Guillochon et al. (2009) Guillochon, J., RamirezRuiz, E., Rosswog, S., & Kasen, D. 2009, Ap. J., 705, 844
 Hamada & Salpeter (1961) Hamada, T., & Salpeter, E. E. 1961, ApJ, 134, 683
 Harris (1991) Harris, W. E. 1991, ARA&A, 29, 543
 Henke et al. (1993) Henke, B. L., Gullikson, E. M., & Davis, J. C. 1993, Atomic Data and Nuclear Data Tables, 54, 181
 Herrmann et al. (2007a) Herrmann, F., Hinder, I., Shoemaker, D., & Laguna, P. 2007a, Classical and Quantum Gravity, 24, 33
 Herrmann et al. (2007b) Herrmann, F., Hinder, I., Shoemaker, D., Laguna, P., & Matzner, R. A. 2007b, Ap. J., 661, 430
 Hinder et al. (2008) Hinder, I., Vaishnav, B., Herrmann, F., Shoemaker, D. M., & Laguna, P. 2008, Phys. Rev. D, 77, 081502
 Husa et al. (2006) Husa, S., Hinder, I., & Lechner, C. 2006, Computer Physics Communications, 174, 983
 Irwin et al. (2010) Irwin, J. A., Brink, T. G., Bregman, J. N., & Roberts, T. P. 2010, ApJ, 712, L1
 Kasen & RamirezRuiz (2009) Kasen, D., & RamirezRuiz, E. 2009, Ap. J., 714, 155
 Katz et al. (2010) Katz, B., Budnik, R., & Waxman, E. 2010, Astrophys. J., 716, 781
 Kawaguchi (2003) Kawaguchi, T. 2003, Ap. J., 593, 69
 Kesden (2011) Kesden, M. 2011, ArXiv eprints
 Kobayashi et al. (2004) Kobayashi, S., Laguna, P., Phinney, E. S., & Meszaros, P. 2004, Astrophys. J., 615, 855
 Komatsu et al. (2011) Komatsu, E., et al. 2011, ApJS, 192, 18
 Krolik & Piran (2011) Krolik, J. H., & Piran, T. 2011, astroph/1106.0923
 Laguna et al. (1993a) Laguna, P., Miller, W. A., & Zurek, W. H. 1993a, Ap. J., 404, 678
 Laguna et al. (1993b) Laguna, P., Miller, W. A., Zurek, W. H., & Davies, M. B. 1993b, Ap. J. Lett., 410, L83
 Lasota (1994) Lasota, J. P. 1994, in NATO ASI Series C, Vol. 417, Theory of Accretion Disks, ed. W. J. D. et al., Garching, 341
 Lee et al. (2005) Lee, W. H., RamirezRuiz, E., & Page, D. 2005, Astrophys.J., 632, 421
 LIGO Scientific Collaboration (2009) LIGO Scientific Collaboration. 2009
 Lodato et al. (2009) Lodato, G., King, A. R., & Pringle, J. E. 2009, Monthly Notices of the Royal Astronomical Society, 392, 332
 Loeb & Ulmer (1997) Loeb, A., & Ulmer, A. 1997, Ap. J., 489, 573
 Löffler et al. (2006) Löffler, F., Rezzolla, L., & Ansorg, M. 2006, Phys. Rev. D, 74, 104018
 Luminet (1985) Luminet, J.P. 1985, Annales de Physique, 10, 101
 Luminet & Marck (1985) Luminet, J.P., & Marck, J.A. 1985, MNRAS, 212, 57
 Luminet & Pichon (1989) Luminet, J.P., & Pichon, B. 1989, A&A, 209, 85
 Magorrian & Tremaine (1999) Magorrian, J., & Tremaine, S. 1999, MNRAS, 309, 447
 Mazzitelli & Dantona (1986) Mazzitelli, I., & Dantona, F. 1986, Ap. J., 311, 762
 McClintock et al. (2011) McClintock, J. E., et al. 2011, Classical and Quantum Gravity, 28, 114009
 McLaughlin (1999) McLaughlin, D. E. 1999, AJ, 117, 2398
 Merritt & Poon (2004) Merritt, D., & Poon, M. Y. 2004, ApJ, 606, 788
 Mukai (1993) Mukai, K. 1993, Legacy, vol. 3, p.2131, 3, 21
 Nauenberg (1972) Nauenberg, M. 1972, ApJ, 175, 417
 Oppenheimer & Volkoff (1939) Oppenheimer, J., & Volkoff, G. 1939, Phys.Rev., 55, 374
 Owocki & Gayley (1997) Owocki, S. P., & Gayley, K. G. 1997, ASP Conf. Proc., 120, 121
 Papaloizou & Pringle (1983) Papaloizou, J. C. B., & Pringle, J. E. 1983, MNRAS, 202, 1181
 Phinney (1989) Phinney, E. S. 1989, in IAU Symposium, Vol. 136, The Center of the Galaxy, ed. M. Morris, 543–+
 Rees (1988) Rees, M. J. 1988, Nature, 333, 523
 Reid (2005) Reid, I. N. 2005, ARA&A, 43, 247
 Rosswog et al. (2009) Rosswog, S., RamirezRuiz, E., & Hix, R. 2009, Astrophys. J., 695, 404
 Rybicki & Lightman (1979) Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley)
 Sadowski (2009) Sadowski, A. 2009, Ap. J. Supp., 183, 171
 Sakimoto & Coroniti (1981) Sakimoto, P. J., & Coroniti, F. V. 1981, Ap. J., 247, 19
 Schnetter et al. (2004) Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Class. Quant. Grav., 21, 1465
 Sesana et al. (2008) Sesana, A., Vecchio, A., Eracleous, M., & Sigurdsson, S. 2008, MNRAS, 391, 718
 Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astron. and Astrophys., 24, 337
 Shapiro & Teukolsky (1986) Shapiro, S. L., & Teukolsky, S. A. 1986, Black Holes, White Dwarfs and Neutron Stars: The Physics of Compact Objects (Wiley: New York)
 Shaviv (2001) Shaviv, N. J. 2001, MNRAS, 326, 126
 Stone & Loeb (2011) Stone, N., & Loeb, A. 2011, ArXiv eprints
 Strubbe & Quataert (2009) Strubbe, L. E., & Quataert, E. 2009, MNRAS, 400, 2070
 Thornburg (2004) Thornburg, J. 2004, Class. Quant. Grav., 21, 743
 Tolman (1939) Tolman, R. C. 1939, Phys.Rev., 55, 364
 Tundo et al. (2007) Tundo, E., Bernardi, M., Hyde, J. B., Sheth, R. K., & Pizzella, A. 2007, ApJ, 663, 53
 Ulmer et al. (1998) Ulmer, A., Paczynski, B., & Goodman, J. 1998, Astron. and Astroph., 333, 379
 Vaishnav et al. (2007) Vaishnav, B., Hinder, I., Herrmann, F., & Shoemaker, D. 2007, Phys. Rev., D76, 084020
 Voges et al. (2000) Voges, W., et al. 2000, IAU Circ., 7432, 1
 Wang & Merritt (2004) Wang, J., & Merritt, D. 2004, Ap. J., 600, 149
 York (1979) York, J. W. 1979, in Sources of Gravitational Radiation, ed. L. L. Smarr (Cambridge, UK: Cambridge University Press), 83–126