Free Access
Issue
A&A
Volume 523, November-December 2010
Article Number A8
Number of page(s) 16
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201014969
Published online 10 November 2010

© ESO, 2010

1. Introduction

Despite the fundamental role they play in gravitational-wave astronomy, no undisputed observational evidence of the existence of supermassive binary black holes (SMBBHs) systems has been found yet. However, circumstantial evidence does exist for a number of potential candidates. This is the case, for instance, of the radio galaxy 0402+379, which shows a projected separation between the two black holes of 7.3pc and a total mass of  ~1.5 × 108   M (Rodriguez et al. 2006). Similarly, the ultraluminous infrared galaxy NGC 6240 shows two optical nuclei and is thought to be in an early merger phase (Komossa et al. 2003). Finally, potential candidates have been suggested in few other cases where the two galaxies are more widely separated, like in the pair IC 694/NGC 3690 hosting two active nuclei as revealed by the presence of two distinct Kα lines in their X-ray spectra (Ballo et al. 2004).

In addition, there is a large family of even more uncertain SMBBH candidates, that are spatially unresolved and whose ultimate nature is, of course, a matter of strong debate. They include the class of X-shaped radio galaxies, in which the observed changes in the orientation of the black hole spin axis could be due to an ongoing merger with a second black hole (Gopal-Krishna et al. 2003); or the class of double-double radio galaxies presenting a pair of double-lobed radio structures that could be the remnants of a SMBBH merger event (Liu et al. 2003); or, finally, the class of sources showing periodicities in the light curves, like in BL Lac Object OJ287 (Komossa 2006). Quite recently, also the quasar SDSSJ0927 (at a redshift z ≈ 0.7) has been identified as a promising SMBBH candidate, with a mass for the primary black hole of M ≈ 2 × 109   M and a semi-major axis of 0.34pc (Dotti et al. 2009).

A strong motivation for studying the merger of supermassive binary black hole systems comes from the fact that their gravitational signal will be detected by the planned Laser Interferometric Space Antenna (LISA), whose optimum sensitivity is placed in the range (10-4 ÷ 0.1)  Hz. Considerable attention has therefore been recently attracted by the possibility of detecting also the electromagnetic (EM) counterpart of these events through the emission coming from the circumbinary accretion disc that is expected to form when the binary is still widely separated. Such a detection will not only act as a confirmation of the gravitational wave (GW) detection, but it will also provide a new tool for testing a number of fundamental astrophysical issues (Haiman et al. 2009). More specifically, it will offer the possibility of testing models of galaxy mergers and accretion discs perturbation, probing basic aspects of gravitational physics, and allowing for the measurement of the Eddington ratio and for the determination of cosmological parameter once the redshift is known (Phinney 2009). In spite of the presence of the disturbing effects due to weak-lensing errors, Kocsis et al. (2006) have computed the average number of quasars in the three-dimensional LISA error volume and have shown that for mergers with masses in the range  ~4 × (105 ÷ 107)   M at redshift z ~ 1, the error box may contain 1 quasar with a luminosity LB ~ (1010 ÷ 1011)   L (see also Kocsis et al. 2008; Kocsis & Loeb 2008).

As a product of this increased interest, a number of studies have been recently carried out to investigate the properties of these EM counterparts either during the stages that precede the merger, or in those following it. As an example, recent work has considered the interaction between the binary and a dense gas cloud (Armitage & Natarajan 2002; van Meter et al. 2010; Bode et al. 2009; Farris et al. 2010; Lodato et al. 2009; Chang et al. 2010) even though astrophysical considerations seem to suggest that during the very final stages of the merger the SMBBH will inspiral in a rather tenuous intergalactic medium. At the same time, other scenarios not involving matter have also been considered. In these cases the SMBBH is considered to be inspiralling in vacuum but in the presence of an external magnetic field which is anchored to the circumbinary disc (Palenzuela et al. 2009; Mösta et al. 2010). The extensive analysis of Mösta et al. (2010), in particular, has shown that, even though the electromagnetic radiation in the lowest  = 2 and m = 2 multipole accurately reflects the gravitational one, the energy emitted in EM waves is 13 orders of magnitude smaller than the one emitted in GW, thus making the direct detection of the two different signals very unlikely. The situation changes in the post-merger phase. In this case, in fact, the EM counterpart is supposed to be mainly due to the radiation from the circumbinary accretion disc, and, because of that, it will contain an imprint of any strong dynamical change produced on the disc by the merger event. There are indeed two major such dynamical effects. The first one is the abrupt reduction of the rest-mass of the binary, emitted away in GWs, which is a function of the binary mass ratio, amounting up to  ≃10% for equal-mass spinning systems (Reisswig et al. 2009). The second one is the recoil velocity of the merged system, resulting in a kick velocity of the resulting black hole with respect to the hosting galaxy (see Bekenstein 1973; Redmount & Rees 1989; for a first discussion of the process and Rezzolla 2009, for a recent review). Leaving aside possible problems due to the actual value of the kicked velocity, which in some cases could be even larger than the escape velocity, it is clear that both events mentioned above can significantly affect the dynamics of the circumbinary disc, mainly as they contribute to the formation and propagation of shocks, thus enhancing the possibility of a strong EM counterpart.

After the first smooth-particle-hydrodynamics approach to the dynamical evolution of circumbinary discs performed by Artymowicz & Lubow (1994), several additional numerical investigations have been proposed in the very recent past. MacFadyen & Milosavljević (2008), for example, performed two dimensional hydrodynamical simulations and studied in detail the evolution of the binary separation and of the disc eccentricity. By perturbing Keplerian orbits of collisionless test particles, on the other hand, Lippai et al. (2008) found a clear spiral shock pattern in the plane of the disc as a response to the kick. By performing pseudo-Newtonian numerical simulations of Keplerian discs O’Neill et al. (2009) have recently questioned the contribution of the shocks to the expected bremsstrahlung emissivity, while Megevand et al. (2009) showed that the intensity of bremsstrahlung luminosity is not much affected by the magnitude of the kick velocity, provided this is less than the smallest orbital speed of the fluid. Although they represent the first fully general relativistic calculations of this process, the simulations of Megevand et al. (2009) used unrealistically small discs which were also placed extremely close to the recoiling black hole. As a considerable improvement over all the previous investigations, Corrales et al. (2010) have carried out a systematic study of the effects of the mass-loss and recoil over a number of α-discs in Newtonian gravity and two-dimensions. While confirming the existence of spiral shocks, they also provided a first realistic estimate of the resulting enhanced luminosity, which can be as large as few  × 1043   erg/s when the disc is assumed to be extremely efficient in radiating any local increase of the temperature. Very interesting results have also been obtained by Rossi et al. (2010), who estimated the maximum disc-to-hole mass ratio that would be stable against fragmentation due to self-gravity to be Md/M ~ 6 × 10-4 for a supermassive black hole with mass M = 106   M. In addition, by performing three-dimensional but Newtonian SPH simulations of geometrically thin discs, they found that the emitted luminosity corresponding to such small disc-to-hole mass ratios is unlikely to make the EM counterpart visible via wide-area sky surveys.

In this paper we present the results of two-dimensional relativistic numerical simulations of extended circumbinary discs in the post-merger phase of the merger, when the disc reacts to the mass loss of the central black hole and to the received kick velocity. By accurately capturing the dynamics of the perturbed disc in the relativistic regime, we investigate the dependence of the accretion rate on the black-hole spin and on the kick velocity. At the same time, we introduce a new technique to locate the shocks that are potentially produced by the recoil and can therefore assess under what conditions a spiral pattern can develop, producing a variability in the accretion rate and, hence, in the luminosity. Our “shock detector” is based on the analysis of the initial states of the Riemann problem solved at each cell interface and can therefore determine the location of the shock with extreme precision, thus revealing that the previously proposed criteria for the occurrence of the shock are often inaccurate.

To compare with the general-relativistic calculations performed by Megevand et al. (2009), our initial models consider small-size discs with an inner radius at r ~ 40 M and an outer one at r ~ 120   M. In addition, however, we also study the dynamics of large-size discs with an inner radius at r ~ 400   M and an outer one at r ~ 4700   M. These configurations have almost Keplerian distributions of angular momentum and are therefore closer to what is believed to be a realistic configuration for a circumbinary disc. Furthermore, because the mass in the discs is always much smaller than the mass of the black hole (i.e. with a mass ratio  ~10-3), we solve the equations of relativistic hydrodynamics in the fixed spacetime of the final black hole.

At first sight it may appear that the use of general-relativistic hydrodynamics is unnecessary when simulating astrophysical systems such as the ones considered here and especially for the case of large-size discs. Such a view, however, does not take into account that much of the dynamics in these EM counterparts takes place near the black-hole horizon, where general relativistic effects are not only large but essential for a correct physical description. Moreover, we do not have any firm theoretical basis to exclude small-size discs and for which the relativistic corrections are non-negligible. Finally, even in a scenario in which gravity could be approximated by the Newton law, we cannot exclude the importance of special relativistic effects.

As all of the above mentioned investigations, also our approach suffers from the absence of a a fully consistent treatment of the radiation transfer, thus allowing only for tentative conclusions about the energetics involved in circumbinary accretion discs. However, we extend to the relativistic framework the strategy reported in Corrales et al. (2010) of performing an isothermal evolution as a tool to extract luminosity curves more realistic than those obtained from thermal bremsstrahlung, although it exaggerates some features of the dynamics, such as the formation of shocks.

The paper is organized as follows. In Sect. 2 we provide the essential information about the numerical code adopted in the simulations. Section 3 describes the physical properties of the initial models, while Sect. 4 highlights the most relevant diagnostic quantities used in the rest of the paper. Section 5 is devoted to the presentations of the results, and, finally, Sect. 6 contains a summary of our work. We assume a signature  {−, +, +, +}  for the space-time metric and we will use Greek letters (running from 0 to 3) for four-dimensional space-time tensor components, while Latin letters (running from 1 to 3) will be employed for three-dimensional spatial tensor components. Moreover, we set c = G = 1 and we extend the geometric units by setting mp/kB = 1, where mp is the mass of the proton, while kB is the Boltzmann constant. In this way the temperature is a dimensionless quantity.

2. Numerical methods

In the stationary spacetime of a Schwarzschild or Kerr black hole we consider the time evolution of a perfect fluid described by the usual energy momentum tensor (1)where uμ is the four velocity of the fluid, gμν is the space-time metric tensor, ρ is the rest-mass density, h = 1 + ϵ + p/ρ the specific enthalpy (including rest-mass energy contribution), ϵ the specific internal energy, p the thermal pressure, related to ρ and ϵ through the usual ideal-gas equation of state (EOS) (2)where γ is the (constant) adiabatic ratio of the gas. We solve the corresponding equations of general relativistic non-dissipative hydrodynamics through the ECHO code (Del Zanna et al. 2007). Because the dynamics of the EM emission takes place on a timescale which is of the order of the orbital one and because the latter is much shorter than the viscous timescale1, the use of inviscid hydrodynamics is indeed a very good approximation.

ECHO adopts a 3 + 1 split of spacetime in which the space-time metric is decomposed according to (3)where α is the lapse function, βi is the shift vector, and γij is the spatial metric tensor. The general-relativistic hydrodynamical equations are written in the following conservative form (4)which is appropriate for numerical integration via standard high-resolution shock-capturing (HRSC) methods developed for the Euler equations. The conservative variables and the corresponding fluxes in the i direction are respectively given by (5)whereas the sources, in any stationary background metric, can be written as (6)where only purely spatial quantities are present. We note that is the determinant of the spatial metric. The relation between the evolved conservative variables (D,Sj,U) and the primitive variables is given by where Γ = (1−v2)−1/2 is the Lorentz factor of the bulk flow with respect to the Eulerian observer associated to the 3 + 1 splitting of the spacetime, and (10)is the fully spatial component of the energy-momentum tensor. In our setup for two dimensional disc simulations we assume the Kerr spacetime metric in Boyer-Lindquist coordinates (i.e. only βφ ≠ 0), with the limiting case of Schwarzschild metric for vanishing black-hole spins, and lay our coordinates in the equatorial plane of the disc (i.e. θ = π/2). The radial numerical grid is discretised by choosing Nr points from rmin to rmax, non-uniformly distributed according to the following scheme where a1 = (rmax − rmin)/a0, a2 = arctana0, while are the coordinate points of the uniform grid from rmin to rmax. In practice, the free parameter a0 controls the extent to which the gridpoints of the original uniform grid are concentrated towards rmin, and we have chosen a0 = 5 in most of our simulations. The actual value of Nr depends on the size of the disc, and it varies between Nr = 600 and Nr = 1200. Outflow boundary conditions are adopted both at rmin and rmax. The azimuthal grid extends from 0 to 2π, with periodic boundary conditions, and Nφ = 200. All runs are performed with a Courant-Friedrichs-Lewy coefficient CFL = 1/2.

The set of hydrodynamics equations is discretised in time with the method of lines and the evolution is performed with a second-order modified Euler scheme. A fifth-order finite-difference algorithm based on an upwind monotonicity preserving filter is employed for spatial reconstruction of primitive variables, whereas a two-wave HLL Riemann solver is used to ensure the shock-capturing properties  (see Del Zanna et al. 2007, for further details). The timestep is generically chosen to be sufficiently small so that the second-order truncation error in time is comparable with the fifth-order one in space.

As a final remark we note that as customary in HRSC methods, we introduce a tenuous and static “atmosphere” in the regions of the fluid outside the initial model for the disc and follow the prescription detailed in Baiotti et al. (2005) for its evolution. In practice we set to zero the velocity field and reset to a pre-defined floor value the rest-mass density of any cell whose density falls below the chosen threshold value. Such a threshold is set to be 8 orders of magnitude below the maximum rest-mass density and we have checked that essentially identical results are obtained when changing this value of one or more orders of magnitude.

Table 1

Main properties of the initial models.

3. Initial models

As initial models we adopt stationary and axisymmetric configurations that are consistent solutions of the relativistic Euler equations and describe a fluid in sub-Keplerian rotation around a Kerr black hole of prescribed mass and spin (Abramowicz et al. 1978). The resulting discs are geometrically thick and they can either have a constant or, more generally, a non-constant radial distribution of the specific angular momentum . In our simulations we have considered both “small-size” models, for which we adopt a constant distribution of , and “large-size’’ models, with a distribution of that, on the equatorial plane, obeys a power law (13)where 𝒮 is chosen to be positive, thus providing a disc rotation that is prograde with respect to the black hole rotation. A detailed description of the equilibrium models for non-constant specific angular momentum discs can be found Daigne & Font (2004). In particular, we have chosen a value of 𝒮 such that the resulting thick discs possess two well defined Keplerian points, namely the “cusp” (which is where matter can accrete onto the black hole) and the “centre” (which is where the pressure has zero gradient); cf. Table 1.

When the exponent q in (13) is chosen close to 1/2, the rotation law tends to the Keplerian one, and the disc flattens towards the equatorial plane. In these circumstances the vertical structure of the disc can be essentially neglected and two-dimensional simulations are therefore indicative of the full three-dimensional dynamics. It is worth mentioning that discs with a rotation law given by (13) have been the subject of a long-standing debate about whether they are subject to the so called “runaway instability” (Abramowicz et al. 1983), which would lead to an exponentially rapid accretion onto the black hole (Font & Daigne 2002b,a; Zanotti et al. 2003; Daigne & Font 2004; Zanotti et al. 2005; Montero et al. 2007). Because the onset and development of this instability depends on the response of the torus to the increased mass of the black hole, simulating this instability accurately requires also the evolution of the Einstein equations. Recent calculations of this type have been performed by Montero et al. (2010) and reveal that the tori are indeed stable irrespective of the angular momentum distribution, thus excluding any role of the runaway instability in the dynamics of the discs simulated here. However, as we will further comment in Sect. 5.1.1, other non-axisymmetric instabilities are possible and have been indeed found.

Table 1 reports the main properties of the models chosen, where the naming convention used allows to easily distinguish the small-size models (S*) from the large-size ones (L*), and where the number * in the name refers to the spin of the black hole, thus ranging between 0.00 and 0.99. As already commented in the Introduction, while the small-size models are particularly suitable for investigating any effect of the black hole spin, the large-size models are those that are (astro)physically more realistic. The inner radius of these large-size models is typically of a few hundreds of gravitational radii and represents the size of the cavity produced by the torque of the SMBBH as estimated from the expression deduced from Table 1 of  Milosavljeć & Phinney (2005) (14)where α-1 ≡ α/0.1, η-1 ≡ η/0.1 ( and being the effective -parameter of thin accretion discs and the radiative efficiency, respectively), Edd is the mass accretion rate in Eddington units and q is the mass ratio between the two coalescing black holes.

By construction, the recoil velocities that can be studied in our setup are those contained in the equatorial (i.e. (r,φ)) plane and because it is much more advantageous to study the dynamics of the disc in a reference frame comoving with the black hole, we impose a net velocity field in addition to the equilibrium orbital one. In practice, at time t = 0 we perform a Lorentz boost of the fluid velocity along the radial direction with φ = 0, thus mimicking a recoil velocity of the black hole in radial direction but with φ = π. We treat the imparted recoil velocity Vk essentially as a free parameter ranging from Vk = 100 km s-1 to Vk = 3000 km s-1, where the latter values are not realistic and serve here only to appreciate the disc dynamics under extreme conditions. We recall, in fact, that the recoil velocities in the orbital plane are expected to be  ≲450 km s-1 (Koppitz et al. 2007; Herrmann et al. 2007; Pollney et al. 2007).

In addition to the recoil, in some initial models we also consider the effects of the mass lost to gravitational waves and which we account by first computing the initial model in the gravitational potential of the full black hole mass, and then evolving it in the gravitational potential of the reduced mass. As a reference value we consider a decrease in the mass of  ~3%, thus corresponding to that obtained from the typical merger of equal-mass spinning black holes with spins anti-aligned with the orbital angular momentum  (see Reisswig et al. 2009, Fig. 11). Higher values of the mass loss do not introduce qualitative changes in the overall dynamics.

A final comment is devoted to the EOS of the initial model, that we chose to be that of a polytrope p = κργ, with γ = 4/3 or γ = 5/3. We recall that a peculiar property of these equilibrium models is that the ratio p/ρ, and therefore the temperature T and the sound speed cs, do not depend on the polytropic constant κ2, which, on the other hand, determines the mass of the disc (see Rezzolla et al. 2003b, Appendix B). As a result, the size of the torus is fixed, also the temperature is uniquely determined and cannot be rescaled further. The last column in Table 1 reports such a temperature at the centre of the torus, rc, as estimated from the ideal-gas EOS via the expression (15)where, we recall, kB is the Boltzmann constant and mp the rest-mass of the proton. In geometric units and with mp/kB = 1 the transformation of the temperature from the dimensionless values to Kelvin degrees is given by (16)

4. Methodology of the analysis

In what follows we discuss in detail the physical quantities computed during the evolution either as representatives of the global evolution or of the local one.

4.1. Global quantities

In addition to the local Eulerian fluid variables, during the evolution we also monitor a few global quantities that are very helpful for interpreting the main properties of the dynamics. These are: the rest-mass, the internal energy and the accretion rate at the innermost radial point of the grid, each of which is computed as Note that when computing the volume integral we consider the discs to have half thickness H which is assumed to be constant in radius, i.e. with H ~ cs/Ω as in the standard thin disc approximation, with cs the sound speed and Ω the orbital velocity.

In addition to (17)–(19) we also compute a few more diagnostic quantities that, on the contrary, rely on simplified assumptions reflecting the fact that the implementation does not account for processes such as radiation transfer and viscous dissipation. In particular, we compute the bremsstrahlung emissivity of the electron-proton collision as (Rybicki & Lightman 1986) where ne and ni ≃ ne are the number densities of electrons and ions (protons), respectively, while T is the equilibrium temperature of both electrons and protons. The bremsstrahlung luminosity is then obtained after performing the volume integral (22)where the large multiplicative factor comes from the fact that both T and ρ in (22) are expressed in geometrized units.

4.2. A relativistic “shock detector”

An obvious expectation, which has been confirmed by all of the numerical simulations to date, is that the as the recoiling black hole will move in the plane of the accretion disc it will introduce spiral shocks which will move outwards on a timescale which is comparable with the orbital one. Because determining the accurate position of the shocks is important to correlate the latter to the EM emission, a number of suggestions have been made in the literature, which have a varying degree of precision. In particular, Lippai et al. (2008); O’Neill et al. (2009); Megevand et al. (2009), all just looked at density and/or pressure gradients to infer the propagation of a spiral caustic and, therefore, of a possible shock (we note that in the collisionless particles treatment of Lippai et al. (2008), the existence of a shock is purely indicative as no shocks can be produced in this approximation). On the other hand, Rossi et al. (2010) used the introduction of an artificial viscosity, which is itself related to local density increases, to identify the location of shocks. Finally, Corrales et al. (2010) used a shock detector present in the FLASH code, which marks a given region as a shocked one if ·v < 0 and if the pressure difference between the monitored zone and at least one of its neighbors exceeds the difference expected from the Rankine-Hugoniot jump condition for a shock of a pre-specified minimum Mach number. While more robust than those considered by the other authors, also this prescription is a delicate one as we will discuss in Sect. 5.1.1.

All of the methods mentioned above contain rather empirical criteria and can fail to detect shocks unless they are very strong. To improve the determination of the location of the shock, even when the latter are arbitrarily weak, we have devised a relativistic “shock detector” which exploits an idea discussed in all its details in Rezzolla & Zanotti (2002) and Rezzolla et al. (2003c), and which consists essentially in the possibility of predicting the outcome of the wave pattern in a Riemann problem. (We note that a similar detector can be prescribed also for non-relativistic flows; the interested reader can find a detailed discussion in Sect. 100 of Landau & Lifshitz 2004.)

To illustrate the logic of our shock detector let us suppose that along a given direction, say the x-direction, two adjacent fluid elements 1 and 2 manifest a jump in the hydrodynamical quantities, such as pressure, density and velocity, thus reproducing the typical conditions of a local Riemann problem. In the absence of magnetic fields, the time evolution of a Riemann problem consists in the propagation along opposite directions of two nonlinear waves, either rarefactions or shocks, separated by a third wave, the contact discontinuity. As a result, a shock front will be produced if the wave pattern generated by the Riemann problem contains at least one shock wave, while the other wave can be a rarefaction wave. As shown by Rezzolla & Zanotti (2002), there is a simple criterion for predicting the occurrence of a wave pattern containing a shock wave and this amounts to the requirement that the relative velocity between the two states 1 and 2 (i.e. between two adjacent fluid cells) is larger than a threshold value (23)where is a function of the thermodynamic states of 1 and 2, while v12 ≡ (v1 − v2)/(1 − v1v2) is the special relativistic expression for the relative velocity. When there are nonzero velocities in the direction tangential to the discontinuity front, the analytic form of is given by (Rezzolla et al. 2003c) (24)where while cs is the sound speed. If, on the contrary, the relative velocity v12 is smaller than , then no shock wave can be produced and the wave pattern of the corresponding Riemann problem consists of two rarefaction waves propagating in opposite directions.

With this idea in mind we have constructed a sensitive shock detector to locate the regions of the disc which are producing a spiral shock. In practice we first select the direction along which we want to monitor the propagation of shock waves. Secondly, since (23) and (24) have been derived in a flat spacetime, we project the velocity field vj in a local tetrad in Boyer-Lindquist coordinates so as to obtain the new components vĵ Thirdly, we calculate and vŷ from and through a simple rotation. Finally, we compute the integral (24) in terms of the hatted Cartesian components and compare the result with v12.

Note that the integral (24) effectively provides the minimum value for the occurrence of a wave pattern containing a single shock wave. In the limit of , in fact, the pressure jump across the shock wave becomes vanishingly small and a single rarefaction wave joining p1 and p2 propagates in the direction opposite to that of the vanishing shock wave. Therefore, when computing (24) we are actually integrating inside the rarefaction wave, that is notoriously a self-similar solution and hence isentropic. This means that in evaluating (24) we can use the isentropic expression for the sound speed (27)where the density ρ is given in terms of p from p = p1(ρ/ρ1)γ.

thumbnail Fig. 1

Rest-mass density distributions (left columns) and shock structure (right columns) and at three different times (i.e. t = 6.07,47.90 and 99.46h) for model S.00 and a recoiling velocity Vk = 300 km s-1. Note that the last panel refers to almost 25 orbital revolutions. The rest-mass density is plotted on logarithmic scale and in cgs units, while the shock structure is obtained by plotting the quantity Sd (see beginning of Sect. 5.1.1 for a definition); shock waves can form in regions where Sd > 0. Note that is very hard to locate a shock by simply looking at the density distribution.

Open with DEXTER

thumbnail Fig. 2

Sound speed normalized to the kick velocity cs/Vk (left panel) and relativistic Mach number ℳ (right panel) at t = 99.46h for model S.00 when subject to a recoil of Vk = 300 km s-1. Note that Vk ≥ cs is not a good criterion for the localization of the shock (cf. left panel) and that no obvious correlation is present between the supersonicity of the flow and the appearance of the shock (cf. right panel).

Open with DEXTER

The procedure described above is completely general and can be proposed as an efficient shock detector for numerical relativistic hydrodynamics. However, two subtleties should also to be taken into account. The first subtlety is that, for more complicated spacetimes or coordinates systems, the flat-spacetime projection (25)–(26) should be replaced by the more general form (28)with given by (Pons et al. 1998) (29)The second subtlety concerns the fact that, because the shock detector validates the inequality (23), it can be arbitrarily sensitive. Although this certainly represents an advantage, one often wishes to disregard the whole class of weak shocks, for which the contribution to the entropy jump is of higher order and Δs ∝ (Δp)3 (Thorne 1973). In these cases the weakest shocks can be filtered out by making the condition (23) somewhat more restrictive and require therefore that a shock is detected if (30)where (31)with being the velocity of the shock wave propagating towards state 2 (see Rezzolla et al. 2003c, for the explicit expression). Because , any value of χ between 0 and 1 will effectively raise the threshold for the detection of the shocks, filtering out the weakest ones; the shocks encountered in the simulations reported here were all rather weak and we have therefore always used χ = 0. The whole procedure is repeated for as many directions as the dimensions of the problem. Finally, a prescription of the relativistic shock detector as adapted for Newtonian fluids is presented in Appendix A.

5. Results

5.1. Small-size models

Although the small-size models are not astrophysically very realistic as they presume the existence of small tori in equilibrium near the recoiling black hole, they serve to set a comparison with the other general-relativistic calculations of Megevand et al. (2009), where similar tori were considered. In addition, by being so close to the black hole, they are helpful in capturing those features of the dynamics that are most influenced by the regions of strong gravity. However, because of their limited extensions and high densities/temperatures (as an example, the model S.00 has ρc = 3.38 × 10-3   g/cm3 and Tc = 7.9 × 108 K) they will not be used to draw any conclusion on the emitted luminosity, which will be instead discussed in more detail when analyzing the large-size models in Sect. 5.2.

5.1.1. Shock Dynamics

The different panels in the left column of Fig. 1 show the rest-mass density at three different times (i.e. t = 6.07,47.90 and 99.46h) for model S.00 and a recoiling velocity Vk = 300 km s-1. Although the imparted velocity is rather small (but close to the maximum possible in the orbital plane), the disc undergoes large variations in size and density, with a shock front that expands from the inner parts of the disc in an initially axisymmetric manner. This is essentially due to the reduction in the black-hole mass and which moves all of the equilibrium orbits to larger radii. As the influence area of the black hole becomes larger and the orbital velocities become comparable with that of the recoil, the disc develops shocks with the characteristic spiral structure discussed also in previous works (Lippai et al. 2008; Corrales et al. 2010; Rossi et al. 2010) and that transports angular momentum outwards. This is shown by the panels in the right column of Fig. 1, which report the location of the shocks as obtained with the procedure illustrated in Sect. 4.2. More specifically, they show the quantity , whereby any positive value of Sd marks a shocked region. Note that the region very close to the black-hole horizon, namely at radii smaller than  ~10   M, is always a highly shocked one. Furthermore, as the evolution proceeds and the disc expands first in response to the mass loss and subsequently to the shocks, very little (if any) of the computational domain is filled by atmosphere, thus removing de-facto any role it can play in the dynamics of the disc.

Figure 2 provides additional information about the dynamics of the disc by showing the local sound speed normalized to the kick velocity, and the relativistic Mach number, the latter computed as ℳ ≡ Γvscs, where is the Lorentz factor of the sound speed (Konigl 1980). Our results shows that the criterion suggested by  Corrales et al. (2010) for the occurrence of shocks, namely Vk ≥ cs, can be rather misleading in the relativistic context. Indeed, as shown by the snapshots at time t = 99.46 h in Fig. 1, and which refers to almost 25 orbital revolutions, a clear spiral shock forms even if the sound speed is more than 30 times larger than the kick velocity. Indeed, the left panel of Fig. 2 seems to suggest that if any, cs/Vk ≫ 1 is possibly a reasonable necessary condition for the approximate location of the shock. In addition, the correlation between the occurrence of a shock and the local sound speed is very weak and this is apparent in the right panel of Fig. 2, where it is clear that the flow is highly supersonic in the inner regions of the disc and mildly subsonic in the outer regions. Yet, the spiral-shock structure extends continuously across the whole disc. We also note that although the precise morphology of the spiral shocks will depend on the spin of the black hole, this dependence is only very weak and all the considerations made above for model S.00 hold true qualitatively also for spinning black holes.

It is finally worth remarking that the shocks formed here are very mild and not relativistic. Even for Vk = 3000 km s-1, the shock velocity maintains a typical value Vs ~ 0.15 and the velocity jump at the shocks is also rather limited, producing a vv ~ 30. This means, for instance, that such shocks are unable of accelerating electrons through the classical mechanism of Bell (1978).

5.1.2. Mass loss and quasi-periodic dynamics

The dynamics of the disc can change considerably if the black hole is assumed to be recoiling with negligible velocity in the orbital plane and only mass loss is taken into account. By considering mass losses in the range 1−10%,  O’Neill et al. (2009) showed that shocks can form even in the absence of a recoil velocity, provided that the mass loss is larger than the half thickness of the disc. The perturbation induced by the mass loss is spherically symmetric and it causes the disc to expand as each fluid element will want to move to the larger radii corresponding to the equilibrium orbit for the given initial angular momentum. Together with this expansion, however, restoring forces will also induce the disc to contract in the effective-potential well of the black hole with the characteristic frequency of the lowest order p-mode, and which is not too different from the epicyclic frequency at the disc centre (Rezzolla et al. 2003b). We recall that the restoring force responsible for the appearance of such p modes is a combination of pressure gradients, centrifugal and gravitational forces, with the last two playing the dominant roles for the discs considered here (Kato 2001; Rezzolla et al. 2003b).

thumbnail Fig. 3

Time evolution of the internal energy when normalized to the its initial value (top panel) and the corresponding power spectrum (bottom panel) in a model with Vk = 0 and with a mass loss of 1% the initial mass of the black hole.

Open with DEXTER

The oscillating behavior induced by the sudden change of the potential well and the subsequent development of the instability is shown in Fig. 3, where the top panel reports the time evolution of a typical global quantity, i.e. the internal energy, when normalized to its initial value. Interestingly, the remarkable periodicity that characterizes the dynamics is the same as found by Zanotti et al. (2003) when studying global modes of oscillation of thick discs around black holes. The corresponding power spectrum is shown in the bottom panel of Fig. 3, obtained through a FFT of the time series for t ≲ 3d reveals the presence of a fundamental mode of oscillations at f ~ 4.17 × 10-5Hz and of two overtones. The first overtone is at o1 ~ 6.28 × 10-5Hz, while the second one, very close to twice the fundamental frequency, o2 ~ 8.37 × 10-5Hz ~ 2f, is produced by nonlinear coupling of the fundamental mode with itself (see Zanotti et al. 2005). Collectively, these modes of oscillations provide a series of modes in the same ratio of the integer numbers 2:3:4 observed in the QPOs of low-mass X-ray binaries containing a black hole (see Remillard & McClintock 2006, for a recent review) and for which a simple model based on basic p-mode oscillations of a small accretion torus orbiting close to the black hole was recently proposed (Rezzolla et al. 2003a), and which remains one of the most convincing explanations of the observed phenomenology (Schnittman & Rezzolla 2006).

It is difficult not to note that the harmonic behaviour shown in the top panel of Fig. 3 is lost at t ≃ 3d and that the internal energy increases monotonically after that. This is due to the onset of a non-axisymmetric instability which produces spiral arms that rapidly spread to cover the whole disc. An instability of this type was not pointed out by Megevand et al. (2009) although they have used similar models and we believe that this is probably because their simulations were interrupted after  ~11h (or  ~6 orbital periods as measured at the point of maximum rest mass density), which is too early for the development of the instability. On the other hand, such type of instabilities in non-Keplerian discs have been discussed by a number of authors, starting from the pioneering work by Papaloizou & Pringle (1984). A detailed comparison between the linear perturbative analysis of these instabilities and two-dimensional numerical simulations in a Schwarzschild spacetime was already proposed more than twenty years ago by Blaes & Hawley (1988), who found the development of the same spiral structures, which transport both mass and angular momentum outwards even in the absence of mass loss3. While we cannot concentrate here on a detailed discussion of these instabilities, it is sufficient to remark that a spiral-shock pattern and all of the associated phenomenology, can be generated even when the recoil velocity is zero.

thumbnail Fig. 4

Mass accretion rate measured at r = rmin for different values of the black hole spin-parameter in the small-size models, i.e. S.00–S.99. Shown in the inset as a function of the black-hole spin is the stationary accretion rate reached after 5 d.

Open with DEXTER

Overall, the dynamics observed for model S.00 suggests that transient oscillating phenomena may exist in the post-merger phase of SMBBH. In this case, the occurrence of QPOs in the accretion and thus in the luminosity of potential EM counterparts, followed then by the development of non-axisymmetric instabilities would be a unique and convincing signature that a SMBBH merger with small recoil velocities has taken place.

5.1.3. Accretion rates

As originally pointed out by Kozlowski et al. (1978), in non-Keplerian discs no viscosity is needed in order to support accretion in the vicinity of the cusp. Figure 4 reports the baryon mass accretion rate measured at r = rmin for different values of the black hole spin-parameter in the small-size models, i.e. S.00–S.99, and provides further support to the interpretation of the shock dynamics discussed above. In particular, it is easy to realize that because of the sudden mass loss and hence of the reduced gravitational attraction, the disc reacts to the excess of angular momentum by expanding. This effect only lasts for a couple of orbital periods after the merger and leads to the large decrease in mass accretion rate shown in Fig. 4 for t ≲ 0.6 d.

thumbnail Fig. 5

Mass accretion rate measured at r = rmin for the small-size model S.00 and different values of the recoil velocity. Shown in the inset as a function of the recoil velocity is the relative baryon mass accreted after 3d.

Open with DEXTER

Subsequently, as the effect of the kick velocity extends to regions of the flow with smaller orbital velocities and becomes dominant, non-axisymmetric density structures form and the perturbed disc starts filling the low density central cavity while increasing the accretion rate. This is reflected by the the large increase in the mass accretion rate shown in Fig. 4 for t ≳ 0.6 d, which is essentially independent of the black-hole spin.

Since our treatment does not include the compensating effect of the radiation drag exerted by the photons on the accreting matter, the accretion rate increases undisturbed reaching values that are  ~6 orders of magnitude above the Eddington limit. After about 6 orbital periods (≈1d)  saturates and after  ~5 d all of the models have lost  ~32% of their mass. Figure 4 also shows that the spin of the black hole has little influence on the dynamics of the disc, although not as little as inferred from Fig. 13 of Megevand et al. (2009). In particular, we have found that after  ~3d the accretion rate slightly decreases when increasing the spin-parameter (see the inset of Fig. 4), so that  |a = 0 ~ 2.5a = 0.99 after  ~5 d of evolution. Interestingly, this dependence of the accretion rate on the spin of the black hole is rather generic and has been found also in other simulations (not reported here) where no perturbation on the flow is introduced.

Of course, the effect of an increasingly large kick velocity on the disc dynamics is much more pronounced in this case, as it emerges from Fig. 5, which reports the accretion rate for different values of Vk. As already found by Megevand et al. (2009), larger recoil velocities tend to anticipate the occurrence of the burst in the accretion rate and, simultaneously, produce larger absolute values of at the burst. After the initial burst, however, and no later than  ~2 days, the accretion rates become nearly independent of the recoil velocity. The differences in are illustrated in the inset of Fig. 5 which shows the relative baryon mass accreted after 3d and which shows a variation of  ~50% at most. It is also interesting to note that there is no clear imprint of the spiral shock dynamics onto the accretion rate. This is due to the fact that the spiral shocks essentially redistribute angular momentum and thus modify the disc structure and dynamics far from the black hole.

thumbnail Fig. 6

Left column: rest-mass density after 60 and 180d for a recoil velocity Vk = 500 km s-1 applied to model L.00. The scale is logarithmic and expressed in cgs units. Right column: shock structure as presented in Fig. 1 for the same panels in the left column. Once again we remark that is very hard to locate a shock by simply looking at the density distribution, especially when they are very weak (here we have used χ = 0 in Eq. (30)). Note also that the temperature distribution is inversionally proportional to that of the density (not shown here).

Open with DEXTER

thumbnail Fig. 7

The same as in Fig. 6 but for a recoil velocity Vk = 3000 km s-1. Note that the spiral-shock structure is never present and that the inner cavity is rapidly filled by accreting gas. In addition, an oblique shock comprising a low-density region is formed in the inner parts of the flow.

Open with DEXTER

5.2. Large-size model

We next switch to discussing the dynamics of the large-size model L.00 for different values of the recoil velocity. We recall that there are at least two different reasons to consider this second class of models. The first one is that they reflect the expectations of a circumbinary disc: they are quasi-Keplerian, extended and at a large distance from the binary. The second one is that by having a lower density, i.e. ρc = 1.38 × 10-10   g/cm3, and hence a lower temperature, i.e. Tc = 8.7 × 106   K, they lead to much more reasonable values of the recoil-enhanced luminosity. In what follows we briefly review the overall dynamics and then concentrate on how to compute a realistic estimate for the luminosity.

5.2.1. Shock dynamics

The overall dynamics of the large-size models is qualitatively similar to that of the smaller counterparts but with three important differences. The first one is a much weaker dependence of the dynamics on the recoil velocity. This is essentially due to the fact that the inner edge of the discs is so far from the recoiling black hole that only extremely large (and unrealistic) values of the recoil velocity induce a significant modification of the orbital velocity. This is shown in Fig. 6, which reports in the left column the rest-mass density after 60 and 180d for a recoil velocity Vk = 500 km s-1 applied to model L.00. The right column shows instead the corresponding shock structure with the same notation used in Fig. 1 and highlights that a clear spiral structure is lost already after 180d, which now corresponds to almost 18 orbital periods. Note that the temperature distribution is inversionally proportional to that of the density (not shown in Fig. 6) and that the central region is filled only with the atmosphere fluid and thus appears as white in the left column. However, several small shocks are produced in this cavity and these are clearly revealed by the shock detector images on the right column which reports the central region as dark; this is just an artifact that does not have a dynamical impact. The behaviour in Fig. 6 should also be contrasted with the spiral-shock structure of model S.00, which instead persisted intact after almost 25 orbital periods (cf. bottom right panel of Fig. 1). This behaviour and the rapid disappearance of the spiral shock structure is further pronounced as the recoil velocity is increased to Vk = 1000 km s-1 (not shown here).

The second difference is that for sufficiently large recoils, namely for Vk ≥ 2000 km s-1, the initial cavity between the black hole horizon and the inner edge of the disc can be filled rapidly by infalling material. This is evident when looking at Fig. 7, which reports the rest-mass density and shock structure after 60 and 180d for a recoil velocity Vk = 3000 km s-1 applied to model L.00. When contrasted with Fig. 6, which has the same spatial extent, it is easy to notice that already after  ~60d, or  ~5 orbital revolutions, the central cavity is filled with high-density material, some of which extends right onto the black hole. Similarly to what already seen for smaller recoils, also the right column of Fig. 7 shows that in this case the spiral structure is not present, although spatially extended shocks are formed both in the inner regions and in the outer parts of the disc. Note that the velocity jump at the shocks is again rather small, being at most Δv ~ 2.5 × 10-4, hence insufficient to accelerate particles through the various acceleration mechanisms involving shock waves.

Interestingly, because the black hole is moving at very large velocities in an ambient fluid which also has a non-negligible angular momentum, a low-density region which resembles a horn-shaped “cavity” is produced in the downstream part of the flow when Vk = 3000 km s-1 at t = 180 d after the merger. This cavity is shown in Fig. 8 and its orientation is directly related to the direction of the recoil velocity. A simple change of sign in the recoil velocity, in fact, would rotate the cavity of 180 degrees around the black hole (not shown in Fig. 8). The formation of such a cavity is noticeable also in the simulations of Rossi et al. (2010), where however it is not discussed. The cavity leads to quasi-periodic variation of the accretion rate as clumps of matter in the downstream of the flow enter the cavity and streams onto the black hole (cf. the oscillations of in Fig. 9 after t ~ 200d for VK = 3000 km s-1). The generation of this flow pattern has an interest of its own, being a non-trivial variant of the Bondi-Hoyle accretion flow onto a moving black hole and will be investigated with greater detail in a distinct work.

The third and final difference is really a combination of the phenomenology discussed above as it manifests itself in terms of the accretion rate. This is shown in Fig. 9, which reports the mass accretion rate at r = rmin for different values of the kick velocity in the large-size model L.00. As already discussed with Fig. 5 for the small-size models, also in these larger discs the accretion rate shows a rapid increase when the black hole “meets” the disc, reaching values  ≃ 10   M/yr, that are well above the Eddington one because of the absence of the radiation-drag contribution. Differently from the smaller discs, however, the lag between the merger and the increase of the accretion rate is not linear but rather increases nonlinearly as the recoil velocity is decreased. This is clearly shown in Fig. 9, where the accretion rate jumps to high values after  ~35,d for Vk = 3000 km s-1, after  ~100,d for Vk = 2000 km s-1 and after almost one year, i.e.  ~330,d for Vk = 1000 km s-1. The reason for this nonlinear response of the disc has to be found in the fact that for smaller recoil velocities the time required by the black hole to cross the initial cavity τcross will be much larger than the typical orbital revolution time. As an example, for Vk = 1000 km s-1 the timescale ratio is τcross/τc ~ 33, so that the disc will have sufficient time to readjust itself to the new gravitational field and thus redistribute its orbital angular momentum. In practice the inner edge of the disc will move to larger radii as a result of the mass loss and of the varied potential, thus delaying nonlinearly the contact with the black hole and thus the steep increase in the accretion rate.

thumbnail Fig. 8

Magnification of the central region of the rest-mass density in model L.00 after 180d and for a recoil velocity Vk = 3000 km s-1. Note that the colormap is slightly different from the one in Fig. 7 to highlight the presence of a cavity.

Open with DEXTER

As a final remark we note that all of the phenomenology discussed here for the model L.00 has been found also for a large-size model around a rapidly spinning black hole, namely L.90. Because the overall differences are minute and of the order of few percent at most, they will not be discussed here. The reasons behind these similarities are rather obvious: the large-size discs have inner radii that are too far from the black hole to be sensitive to the spin-induced corrections which decay much more rapidly and as 1/r3.

thumbnail Fig. 9

Mass accretion rate at r = rmin for different values of the kick velocity in the large-size model L.00.

Open with DEXTER

5.2.2. EM luminosities

Because of the relatively high temperature of the gas and of the generation of a shock pattern, thermal bremsstrahlung is thought to be an efficient emission mechanism through which circumbinary discs may become visible in the electromagnetic spectrum (Megevand et al. 2009; Corrales et al. 2010; Bode et al. 2009; Anderson et al. 2010). However, thermal-bremsstrahlung emission from circumbinary is affected by a serious problem which has been so far underestimated or not sufficiently emphasized. This has to do with the fact that bremsstrahlung cooling time is too short (Corrales et al. 2010) or, stated differently, that the internal energy budget of the emitting gas is not large enough to allow for the bremsstrahlung emission to last but for a few seconds. This can be easily estimated as tcool = Eint/LBR, with Eint and LBR obtained from (18) and (22), respectively. For the large model L.00 we have Eint ~ 3.4 × 1050   erg and LBR ~ 2.8 × 1049   erg/s at time t = 0 and this estimate remains of the same order of magnitude during the evolution. This yields to tcool ≃ 12 s. The situation is even worse if we consider the transition to the relativistic regime. In this case, in fact, not only the bremsstrahlung emissivity is increased by a factor4  ∝  [1 + 4.4T/(1010   K)]  (Rybicki & Lightman 1986), but also the collisions between particles of the same species start contributing significantly to the bremsstrahlung emission (Svensson 1982) through radiation in moments other than the electric dipole (which is strictly zero for particles of the same species, Krolik 1999).

Of course, there are also other factors that can work in favour of a bremsstrahlung emission and which we have not taken into account. A first one is that we have neglected the thermal bremsstrahlung absorption, which is likely to enhance significantly the bremsstrahlung cooling time by acting as a source of additional internal energy. Moreover, it is also possible that the spiral shock originating from the very central region dissipates considerably as it propagates outwards, hence confining the bulk of the bremsstrahlung luminosity from within a very small portion of the disc (we recall that the bremsstrahlung luminosity is proportional to the volume integral over the emitting source). Overall, while we cannot rule out thermal bremsstrahlung as an emission mechanism from the circumbinary disc, it is also evident to us that the luminosity estimates made so far without a proper treatment of radiation transfer are excessively optimistic.

A second possible estimate of the luminosity is given by the accretion-powered luminosity, Lacc = ηc2, where η is the radiative efficiency. However, lacking any treatment of the radiation pressure effects, such an estimate would only provide misleading conclusions and therefore we will note use it hereafter.

A third and possibly more accurate estimate of the luminosity can be made by assuming that all the changes in the temperature that are due to a local compression will be dissipated as radiation. This idea, proposed in Newtonian physics by Corrales et al. (2010), can be summarized and implemented in a general relativistic context as detailed below.

Consider the evolution of the disc with an equation of state p(T) = ρkbT/mp and a specific internal energy given by , where the last equality has been obtained for γ = 5/3. In general, there is no necessity to evolve the energy equation in an isothermal evolution since the energy can be computed directly from the temperature and the latter is constant by construction. However, the internal energy can be nevertheless evolved in time with the only aim of computing the difference ρ [ϵ − ϵ(T)] , which is then assumed to be radiated instantaneously. The relativistic equation for the evolution of the total internal energy density e ≡ ρ(1 + ϵ) is (Anile 1990) (32)where Θ ≡ ∇μuμ is the expansion of the fluid. The continuity equation ∇μ(ρuμ) = 0 can then be used to rewrite Eq. (32) as (33)One aspect to note is that Eq. (33) is not written in a conservative form because of the derivatives on the right hand side acting on the flow variables. While this is not ideal within our formulation of the hydrodynamics equations, the modifications are minimal. Indeed, since the Lorentz factor does not change significantly during the evolution, it is reasonable to neglect the term  ∝ , while the spatial derivatives of the term can be treated with standard finite-difference methods without a significant loss of accuracy across the discontinuities.

In practice we have assumed an initial temperature for the disc which is uniform in space and set it to be , where Tc is the maximum temperature at the center of the disc (cf. Table 1). An estimate of the luminosity is then trivially computed by performing at each timestep a volume integration of the difference ρ [ϵ − ϵ(T)]  and by dividing it by the simulation timestep. Finally, we reset the specific internal energy to ϵ = ϵ(T0), so as to guarantee that the evolution is effectively isothermal.

thumbnail Fig. 10

Top panel: luminosity computed in the isothermal evolution Lisot for different values of the kick velocity for model L.00 and for a polytropic index γ = 4/3. Note the presence of a peak at about  ~20d after the merger, of a binary with total mass M ≃ 106   M and of a persistent luminosity for several days at values which are a factor of a few smaller. Bottom panel: comparison of Lisot for model L.00 when computed with a polytropic index γ = 4/3 (thick lines) or γ = 5/3 (thin lines). The comparison is made for two reference recoil velocities and shows that the results are very similar, although a stiffer EOS leads to slightly larger luminosities.

Open with DEXTER

The top panel of Fig. 10 shows the luminosity computed in this way for model L.00 with a polytropic index γ = 4/3 and for different values of the recoil velocity. Following Corrales et al. (2010), the values reported have been computed by neglecting the negative contributions to the luminosity that are produced in regions experiencing rarefactions. Because when accounted for these negative contributions typically yield values that are one order of magnitude smaller, the values in Fig. 10 should be taken as upper limits to the emitted luminosity (see Corrales et al. 2010). Clearly, the evolution of the emitted energy has a peak that is larger for stronger recoils and that appears at t ~ 33,18, and 7 d after the merger for Vk = 300,1000 and 3000 km s-1, respectively. While the peaks are the consequence of the strong shocks that are produced in the inner parts of the disc as the latter approaches the black hole, the asymptotic values of the isothermal luminosity are instead produced by the local compressions in the disc. As such, the peaks in the luminosity are not related to the encounter of the black hole with the disc and therefore they are not correlated with the increase in the mass accretion rate, which in general takes place at later times (cf. right panel of Fig. 11). The bottom panel of Fig. 10, on the other hand, shows a comparison of Lisot for model L.00 when computed with a polytropic index γ = 4/3 (thick lines) or γ = 5/3 (thin lines). The comparison is made for two reference recoil velocities of Vk = 300 km s-1 and Vk = 1000 km s-1 and it shows that the results are very similar, although a stiffer EOS leads to slightly larger luminosities.

Overall, our estimates of the luminosity computed within the isothermal evolution approximation confirm those by Corrales et al. (2010) even though the temperature of our large size model is one order of magnitude larger than that reported in Fig. 7 of Corrales et al. (2010). While we believe that these estimates of the luminosity are the most reasonable ones that can be obtained with a code that is intrinsically unable to account for radiation losses, we should also stress that the isothermal evolution by itself provides a less reliable description of the overall dynamics. This is evident, for example, from Fig. 11, whose left panel offers a comparison of the rest-mass density profiles along φ = 0 as obtained with the standard evolution of the energy equation (red solid line) and with the isothermal evolution (blue dashed line) for the large-size model L.00 and a recoil of Vk = 1000 km s-1. Note that the isothermal evolution tends to increase the density gradients, especially during the initial phases, which is also when the peak in the luminosity appears in the light curves. Also shown in the right panel of Fig. 11 is the mass accretion rate at r = rmin for the isothermal evolution and different values of the kick velocity. This panel shows that since the discs are more sensitives to compressions, their response to the recoil is different and, in particular, takes place earlier than in the full-evolution case (cf. with Fig. 9). Furthermore, with the exception of the very large kick, the mass accretion rate does not show a very large jump, but rather a first small jump followed by smooth increase to the final asymptotic behaviour which is reached at times comparable to those of the full evolution. This is clearly the case for Vk = 1000 km s-1 and is due to the fact that accretion starts earlier for this matter whose energy has been decreased by the radiative losses.

thumbnail Fig. 11

Left panel: comparison of the rest-mass density along the φ = 0 direction as obtained with the standard evolution of the energy equation (red solid line) and with the isothermal evolution (blue dashed line), at two different times as shown in the two panels. The data refers to model L.00 with a recoil of 1000 km s-1. Right panel: mass accretion rate at r = rmin for the isothermal evolution and different values of the kick velocity in the large-size model L.00. This panel should be compared with Fig. 9 and which refers to an evolution of the energy equation.

Open with DEXTER

In conclusion, and in spite of the caveats made above, we believe that luminosities as large as few L ≃ 1043   erg/s should be expected at about  ~30 d after the merger of a binary with total mass M ≃ 106   M and that these luminosities should persist for several days to values which are a factor of a few smaller.

6. Conclusions

We have presented the results of two-dimensional general relativistic numerical simulations of small and extended circumbinary discs in the post-merger phase of the merger, when the disc reacts to the mass loss of the central black hole and to its recoil velocity. Our analysis benefitted from being able to capture accurately the dynamics of the perturbed disc in the relativistic regime, thus allowing us to investigate the dependence of the accretion rate on the black-hole spin and on the kick velocity. Furthermore, by considering discs that are quasi-Keplerian, extended and at a large distance from the binary, we were able to consider realistic scenarios even in the general relativistic framework.

Another important aspect of our work is the use of a novel and accurate technique to construct a “shock detector” and hence to determine where, within the flow, the shocks produced by the recoil are located. This, in turn, has allowed us to assess under what conditions a spiral-shock pattern can develop, produce a variability in the accretion rate and, hence, in the luminosity. Our relativistic shock detector (for which we also present a Newtonian equivalent in Appendix A) is based on the analysis of the initial states of the Riemann problem solved at each cell interface and can therefore determine the location of the shock with the same resolution as that of the spatial grid, revealing that the previously proposed criteria for the occurrence of the shock are often inaccurate.

Overall, we can confirm within a general relativistic regime many of the results found previously in Newtonian or pseudo-Newtonian gravity. More specifically, we find that for discs that are sufficiently small and close to the black-hole, a regular spiral-shock develops as a result of the recoiling black hole. The strength, shape and persistence of the shocks, however, depend sensitively on both the size of the tori and on the intensity of the recoil. As a result, while the spiral shock is stable over many orbital periods in the case of small discs subject to small recoils, it never develops or is rapidly destroyed in discs that are large and subject to large recoil velocities. It is worth noting that the typical velocity jumps at the shocks are Δv ≲ 2.5 × 10-4, even with a kick velocity Vk = 3000   km s-1. It is therefore possible that such shocks may be damped by dissipative viscous processes or radiative losses.

Besides the interesting shock properties described above, we also found that the disc dynamics is only very weakly dependent on the black-hole spin. The latter influences only the small-size tori by modifying the accretion rate and by leading to smaller accreted masses per unit time for more rapidly spinning black holes. Finally, we found that even in the limit of a vanishing kick velocity and as long as a mass loss is present, the disc goes through a phase of regular oscillations characterized by the excitation of the p modes of the disc; this is then followed by the appearance of a spiral-shock pattern generated by non-axisymmetric instabilities affecting the disc. This opens the door to the possibility that quasi-periodic oscillations are observed in the post-merger phase of SMBBH.

Computing the EM counterpart to the merger event represents of course a fundamental aspect of our investigation. In contrast with other works, however, we have questioned the estimates of the bremsstrahlung luminosity when computed without properly taking into account the radiation transfer. The energetic reservoir available for the EM emission is, in fact, too small to yield but unrealistically short cooling times. While we cannot rule out thermal bremsstrahlung as the main EM emission mechanism, we believe that the bremsstrahlung-luminosity estimates made so far without a proper treatment of radiation transfer are excessively optimistic. A somewhat realistic estimates of the emitted luminosity can be obtained by assuming that the internal-energy enhancements due to local compressions in the perturbed disc are immediately radiated away, i.e. by considering an “isothermal” evolution like the one recently investigated in Newtonian physics by Corrales et al. (2010). In this case, and despite the fact that the isothermal evolution tends enhance the compressibility of the fluid, we find that the energy emitted can reach a peak value above L ≃ 1043  erg/s at about  ~30 d after the merger of a binary with total mass M ≃ 106   M and persist for several days at values which are a factor of a few smaller. If confirmed by more sophisticated calculations such a signal could represent an interesting EM counterpart of the merger of binary black-hole system.

As a final remark we note that while a rather robust picture is emerging from the collective work done so far on the post-merger dynamics of the circumbinary disc around a SMBBH, much remains to be done to compute realistically the resulting EM emission. Important improvements to the treatment considered here must include the presence of a magnetic field, a proper treatment of the radiation transfer and the extension to three-dimensional calculations. All of these will be the focus of our future work on this subject.


1

We recall that in geometrically thin accretion discs the local viscous timescale is given by , where is the standard alpha parameter ad H is the half-thickness of the disc.

2

The argument consists in proving that the function Ψ ≡ κ(n + 1)Θ, where γ = 1 + 1/n and ρ = Θn, does not depend on κ. From this it follows that p/ρ = κΘ does not depend on κ either.

3

We have verified that the spiral arms do indeed develop also in the absence of a mass loss or of a recoil and that also in this case the instability takes place after  ~4d.

4

It should be remarked, however, that when the electron become relativistic, i.e. for T ≥ 5.9 × 109 K, other emission mechanisms, such as inverse Compton or synchrotron (if a magnetic field is also present), are generally more efficient.

Acknowledgments

We are grateful to Marek Abramowicz and Constanze Roedig for useful discussion. The computations were performed on the Damiana cluster at the AEI and on the IBM/SP6 of CINECA (Italy) through the “INAF-CINECA” agreement 2008-2010. This work was supported in part by the DFG grant SFB/Transregio 7.

References

  1. Abramowicz, M., Jaroszynski, M., & Sikora, M. 1978, A&A, 63, 221 [NASA ADS] [Google Scholar]
  2. Abramowicz, M. A., Calvani, M., & Nobili, L. 1983, Nature, 302, 597 [NASA ADS] [CrossRef] [Google Scholar]
  3. Anderson, M., Lehner, L., Megevand, M., & Neilsen, D. 2010, Phys. Rev. D, 81, 044004 [NASA ADS] [CrossRef] [Google Scholar]
  4. Anile, A. M. 1990, Relativistic Fluids and Magneto-fluids (Cambridge University Press) [Google Scholar]
  5. Armitage, P. J., & Natarajan, P. 2002, ApJ, 567, L9 [NASA ADS] [CrossRef] [Google Scholar]
  6. Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [NASA ADS] [CrossRef] [Google Scholar]
  7. Baiotti, L., Hawke, I., Montero, P. J., et al. 2005, Phys. Rev. D, 71, 024035 [NASA ADS] [CrossRef] [Google Scholar]
  8. Ballo, L., Braito, V., Della Ceca, R., et al. 2004, ApJ, 600, 634 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bekenstein, J. D. 1973, ApJ, 183, 657 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bell, A. R. 1978, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
  11. Blaes, O. M., & Hawley, J. F. 1988, ApJ, 326, 277 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bode, T., Haas, R., Bogdanovic, T., Laguna, P., & Shoemaker, D. 2009 [Google Scholar]
  13. Chang, P., Strubbe, L. E., Menou, K., & Quataert, E. 2010, MNRAS, 407, 2007 [NASA ADS] [CrossRef] [Google Scholar]
  14. Corrales, L. R., Haiman, Z., & MacFadyen, A. 2010, 404, 947 [Google Scholar]
  15. Daigne, F., & Font, J. A. 2004, MNRAS, 349, 841 [NASA ADS] [CrossRef] [Google Scholar]
  16. Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, A&A, 473, 11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Dotti, M., Montuori, C., Decarli, R., et al. 2009, MNRAS, 398, L73 [NASA ADS] [Google Scholar]
  18. Farris, B. D., Liu, Y. T., & Shapiro, S. L. 2010, Phys. Rev. D, 81, 084008 [NASA ADS] [CrossRef] [Google Scholar]
  19. Font, J. A., & Daigne, F. 2002a, ApJ, 581, L23 [NASA ADS] [CrossRef] [Google Scholar]
  20. Font, J. A., & Daigne, F. 2002b, MNRAS, 334, 383 [NASA ADS] [CrossRef] [Google Scholar]
  21. Gopal-Krishna,Biermann, P. L., & Wiita, P. J. 2003, ApJ, 594, L103 [NASA ADS] [CrossRef] [Google Scholar]
  22. Haiman, Z., Kocsis, B., Menou, K., Lippai, Z., & Frei, Z. 2009, Classical and Quantum Gravity, 26, 094032 [NASA ADS] [CrossRef] [Google Scholar]
  23. Herrmann, F., Hinder, I., Shoemaker, D., Laguna, P., & Matzner, R. A. 2007, ApJ, 661, 430 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kato, S. 2001, Publications of the Astronomical Society of Japan, 53, 1 [NASA ADS] [CrossRef] [Google Scholar]
  25. Kocsis, B., & Loeb, A. 2008, Phys. Rev. Lett., 101, 041101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  26. Kocsis, B., Frei, Z., Haiman, Z., & Menou, K. 2006, ApJ, 637, 27 [NASA ADS] [CrossRef] [Google Scholar]
  27. Kocsis, B., Haiman, Z., & Menou, K. 2008, ApJ, 684, 870 [NASA ADS] [CrossRef] [Google Scholar]
  28. Komossa, S. 2006, Mem. Soc. Astron. Ital., 77, 733 [NASA ADS] [Google Scholar]
  29. Komossa, S., Burwitz, V., Hasinger, G., et al. 2003, ApJ, 582, L15 [NASA ADS] [CrossRef] [Google Scholar]
  30. Konigl, A. 1980, Phys. Fluids, 23, 1083 [NASA ADS] [CrossRef] [Google Scholar]
  31. Koppitz, M., Pollney, D., Reisswig, C., et al. 2007, Phys. Rev. Lett., 99, 041102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  32. Kozlowski, M., Jaroszynski, M., & Abramowicz, M. A. 1978, A&A, 63, 209 [NASA ADS] [Google Scholar]
  33. Krolik, J. H. 1999, Active galactic nuclei: from the central black hole to the galactic environment (Princeton University Press) [Google Scholar]
  34. Landau, L. D., & Lifshitz, E. M. 2004, Fluid Mechanics, Course of Theoretical Physics (Oxford: Elsevier Butterworth-Heinemann), 6 [Google Scholar]
  35. Lippai, Z., Frei, Z., & Haiman, Z. 2008, ApJ, 676, L5 [NASA ADS] [CrossRef] [Google Scholar]
  36. Liu, F., Wu, X., & Cao, S. 2003, MNRAS, 340, 411 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lodato, G., Nayakshin, S., King, A. R., & Pringle, J. E. 2009, MNRAS, 398, 1392 [NASA ADS] [CrossRef] [Google Scholar]
  38. MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83 [NASA ADS] [CrossRef] [Google Scholar]
  39. Megevand, M., Anderson, M., Frank, J., et al. 2009, Phys. Rev. D, 80, 024012 [NASA ADS] [CrossRef] [Google Scholar]
  40. Milosavljeć, M., & Phinney, E. S. 2005, ApJ, 622, L93 [NASA ADS] [CrossRef] [Google Scholar]
  41. Montero, P. J., Font, J. A., & Shibata, M. 2010, Phys. Rev. Lett. D, 104, 191101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  42. Montero, P. J., Zanotti, O., Font, J. A., & Rezzolla, L. 2007, MNRAS, 378, 1101 [NASA ADS] [CrossRef] [Google Scholar]
  43. Mösta, P., Palenzuela, C., Rezzolla, L., et al. 2010, Phys. Rev. D, 81, 064017 [NASA ADS] [CrossRef] [Google Scholar]
  44. O’Neill, S. M., Miller, M. C., Bogdanović, T., Reynolds, C. S., & Schnittman, J. D. 2009, ApJ, 700, 859 [NASA ADS] [CrossRef] [Google Scholar]
  45. Palenzuela, C., Anderson, M., Lehner, L., Liebling, S. L., & Neilsen, D. 2009, Phys. Rev. Lett., 103, 081101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  46. Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 [NASA ADS] [CrossRef] [Google Scholar]
  47. Phinney, E. S. 2009, Astron., 2010, 235 [NASA ADS] [Google Scholar]
  48. Pollney, D. et al. 2007, Phys. Rev. D, 76, 124002 [NASA ADS] [CrossRef] [Google Scholar]
  49. Pons, J. A., Font, J. A., Ibanez, J. M., Marti, J. M., & Miralles, J. A. 1998, A&A, 339, 638 [NASA ADS] [Google Scholar]
  50. Redmount, I. H., & Rees, M. J. 1989, Comments on Astrophysics, 14, 165 [NASA ADS] [Google Scholar]
  51. Reisswig, C., Husa, S., Rezzolla, L., et al. 2009, Phys. Rev. D, 80, 124026 [NASA ADS] [CrossRef] [Google Scholar]
  52. Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49 [NASA ADS] [CrossRef] [Google Scholar]
  53. Rezzolla, L. 2009, Class. Quant. Grav., 26, 094023 [NASA ADS] [CrossRef] [Google Scholar]
  54. Rezzolla, L., & Zanotti, O. 2002, Phys. Rev. Lett., 89, 114501 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  55. Rezzolla, L., Yoshida, S., Maccarone, T. J., & Zanotti, O. 2003a, MNRAS, 344, L37 [NASA ADS] [CrossRef] [Google Scholar]
  56. Rezzolla, L., Yoshida, S., & Zanotti, O. 2003b, MNRAS, 344, 978 [NASA ADS] [CrossRef] [Google Scholar]
  57. Rezzolla, L., Zanotti, O., & Pons, J. A. 2003c, J. Fluid Mech., 479, 199 [NASA ADS] [CrossRef] [Google Scholar]
  58. Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49 [NASA ADS] [CrossRef] [Google Scholar]
  59. Rossi, E. M., Lodato, G., Armitage, P. J., Pringle, J. E., & King, A. R. 2010, MNRAS, 401, 2021 [NASA ADS] [CrossRef] [Google Scholar]
  60. Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics (Wiley-VCH) [Google Scholar]
  61. Schnittman, J. D., & Rezzolla, L. 2006, ApJ, 637, L113 [NASA ADS] [CrossRef] [Google Scholar]
  62. Svensson, R. 1982, ApJ, 258, 335 [NASA ADS] [CrossRef] [Google Scholar]
  63. Thorne, K. S. 1973, ApJ, 179, 897 [NASA ADS] [CrossRef] [Google Scholar]
  64. van Meter, J. R., Wise, J. H., Miller, M. C., et al. 2010, ApJ, 711, L89 [NASA ADS] [CrossRef] [Google Scholar]
  65. Zanotti, O., Rezzolla, L., & Font, J. A. 2003, Mon. Not. Roy. Soc., 341, 832 [NASA ADS] [CrossRef] [Google Scholar]
  66. Zanotti, O., Font, J. A., Rezzolla, L., & Montero, P. J. 2005, MNRAS, 356, 1371 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Newtonian version of the shock detector

In this Appendix we provide the basic expressions that allow to build the Newtonian version of the relativistic shock detector described in Sect. 4.2. The logic, of course, is exactly the same and it consists of the following steps.

  • 1.

    First choose the direction, say the x-direction, along which to monitor the generation of shock waves.

  • 2.

    Given any two adjacent cells, label with 1 the cell with higher pressure and with 2 the other one.

  • 3.

    Compute the relative velocity v12 = v1 − v2 along the x-direction and compare it with the value (see Landau & Lifshitz 2004, Sect. 100) (A.1)where γ is the adiabatic index of the gas and cs(p1) is the sound speed in the state 1. It is interesting to note that in the Newtonian case the tangential velocities do not affect the actual value of the threshold , which, therefore, maintains the same expression as for one-dimensional flows. Indeed, the fact that in the relativistic regime the threshold given by (24) does depend on the tangential velocities is at the origin of the appearance of new relativistic effects discussed in Rezzolla & Zanotti (2002).

  • 4.

    A shock is therefore detected if (A.2)

The procedure is repeated for as many directions as the dimensions of the problem. As commented in the main text, we note here too that it may be necessary at times to filter-out the smallest shocks and therefore make the condition for shock detection more stringent. In analogy with Eq. (30), the following condition can then be implemented (A.3)where (see Landau & Lifshitz 2004, Sect. 100) (A.4)and where the parameter χ ∈  [0,1]  needs to be tuned to the desired level of accuracy.

All Tables

Table 1

Main properties of the initial models.

All Figures

thumbnail Fig. 1

Rest-mass density distributions (left columns) and shock structure (right columns) and at three different times (i.e. t = 6.07,47.90 and 99.46h) for model S.00 and a recoiling velocity Vk = 300 km s-1. Note that the last panel refers to almost 25 orbital revolutions. The rest-mass density is plotted on logarithmic scale and in cgs units, while the shock structure is obtained by plotting the quantity Sd (see beginning of Sect. 5.1.1 for a definition); shock waves can form in regions where Sd > 0. Note that is very hard to locate a shock by simply looking at the density distribution.

Open with DEXTER
In the text
thumbnail Fig. 2

Sound speed normalized to the kick velocity cs/Vk (left panel) and relativistic Mach number ℳ (right panel) at t = 99.46h for model S.00 when subject to a recoil of Vk = 300 km s-1. Note that Vk ≥ cs is not a good criterion for the localization of the shock (cf. left panel) and that no obvious correlation is present between the supersonicity of the flow and the appearance of the shock (cf. right panel).

Open with DEXTER
In the text
thumbnail Fig. 3

Time evolution of the internal energy when normalized to the its initial value (top panel) and the corresponding power spectrum (bottom panel) in a model with Vk = 0 and with a mass loss of 1% the initial mass of the black hole.

Open with DEXTER
In the text
thumbnail Fig. 4

Mass accretion rate measured at r = rmin for different values of the black hole spin-parameter in the small-size models, i.e. S.00–S.99. Shown in the inset as a function of the black-hole spin is the stationary accretion rate reached after 5 d.

Open with DEXTER
In the text
thumbnail Fig. 5

Mass accretion rate measured at r = rmin for the small-size model S.00 and different values of the recoil velocity. Shown in the inset as a function of the recoil velocity is the relative baryon mass accreted after 3d.

Open with DEXTER
In the text
thumbnail Fig. 6

Left column: rest-mass density after 60 and 180d for a recoil velocity Vk = 500 km s-1 applied to model L.00. The scale is logarithmic and expressed in cgs units. Right column: shock structure as presented in Fig. 1 for the same panels in the left column. Once again we remark that is very hard to locate a shock by simply looking at the density distribution, especially when they are very weak (here we have used χ = 0 in Eq. (30)). Note also that the temperature distribution is inversionally proportional to that of the density (not shown here).

Open with DEXTER
In the text
thumbnail Fig. 7

The same as in Fig. 6 but for a recoil velocity Vk = 3000 km s-1. Note that the spiral-shock structure is never present and that the inner cavity is rapidly filled by accreting gas. In addition, an oblique shock comprising a low-density region is formed in the inner parts of the flow.

Open with DEXTER
In the text
thumbnail Fig. 8

Magnification of the central region of the rest-mass density in model L.00 after 180d and for a recoil velocity Vk = 3000 km s-1. Note that the colormap is slightly different from the one in Fig. 7 to highlight the presence of a cavity.

Open with DEXTER
In the text
thumbnail Fig. 9

Mass accretion rate at r = rmin for different values of the kick velocity in the large-size model L.00.

Open with DEXTER
In the text
thumbnail Fig. 10

Top panel: luminosity computed in the isothermal evolution Lisot for different values of the kick velocity for model L.00 and for a polytropic index γ = 4/3. Note the presence of a peak at about  ~20d after the merger, of a binary with total mass M ≃ 106   M and of a persistent luminosity for several days at values which are a factor of a few smaller. Bottom panel: comparison of Lisot for model L.00 when computed with a polytropic index γ = 4/3 (thick lines) or γ = 5/3 (thin lines). The comparison is made for two reference recoil velocities and shows that the results are very similar, although a stiffer EOS leads to slightly larger luminosities.

Open with DEXTER
In the text
thumbnail Fig. 11

Left panel: comparison of the rest-mass density along the φ = 0 direction as obtained with the standard evolution of the energy equation (red solid line) and with the isothermal evolution (blue dashed line), at two different times as shown in the two panels. The data refers to model L.00 with a recoil of 1000 km s-1. Right panel: mass accretion rate at r = rmin for the isothermal evolution and different values of the kick velocity in the large-size model L.00. This panel should be compared with Fig. 9 and which refers to an evolution of the energy equation.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.