Gravitational-wave recoil of merging supermassive black holes (SMBHs) may influence the co-evolution of SMBHs and their host galaxies. We examine this possibility using smoothed particle hydrodynamic/N-body simulations of gaseous galaxy mergers, including SMBH accretion, in which the merged BH receives a recoil kick. This enables us to follow the trajectories and accretion of recoiling BHs in self-consistent, evolving merger remnants. In contrast to recent studies on similar topics, we conduct a large parameter study, generating a suite of over 200 simulations with more than 60 merger models and a range of recoil velocities (vk). With this, we can identify systematic trends in the behaviour of recoiling BHs. Our main results are as follows. (1) BHs kicked at nearly the central escape speed (vesc) may oscillate on large orbits for up to a Hubble time, but in gas-rich mergers, BHs kicked with up to ∼0.7vesc may be confined to the central few kpc of the galaxy, owing to gas drag and steep central potentials. (2) vesc in gas-rich mergers may increase rapidly during final coalescence, in which case recoil trajectories may depend sensitively on the timing of the BH merger relative to the formation of the central potential well. Delays of even a few times 107 yr in the BH merger time may substantially suppress recoiling BH motion for a fixed kick speed. (3) Recoil events generally reduce the lifetimes of bright active galactic nuclei (AGNs), but short-period recoil oscillations may actually extend AGN lifetimes at lower luminosities. (4) Bondi–Hoyle accretion dominates recoiling BH growth for vk/vesc≲ 0.6–0.8; for higher vk, the BH accretes primarily from its ejected gas disc, which we model as a time-dependent accretion disc. (5) Kinematically offset AGNs may be observable either immediately after a high-velocity recoil or during pericentric passages through a gas-rich remnant. In either case, AGN lifetimes may be up to ∼10 Myr (for Δv > 800 km s−1), though the latter generally have higher luminosities. (6) Spatially offset AGNs can occur for vk≳ 0.5–0.7vesc (for ΔR > 1 kpc); these generally have low luminosities and lifetimes of ∼1–100 Myr. (7) Rapidly recoiling BHs may be up to approximately five times less massive than their stationary counterparts. These mass deficits lower the normalization of the MBH–σ* relation and contribute to both intrinsic and overall scatter. (8) Finally, the displacement of AGN feedback after a recoil event enhances central star formation rates in the merger remnant, thereby extending the starburst phase of the merger and creating a denser, more massive stellar cusp.
Within the hierarchical structure formation paradigm, a significant fraction of galaxy growth occurs via successive mergers. Given the ample evidence that most, if not all, local galaxies host central supermassive black holes (SMBHs; e.g. Kormendy & Richstone 1995; Richstone et al. 1998; Ferrarese & Ford 2005), major galaxy mergers (those with mass ratios q≳ 0.25) should inevitably result in the formation of SMBH binaries. The fate of these binaries is somewhat uncertain and likely depends on the conditions within their host galaxies. In highly spheroidal, gas-poor galaxies, these binaries may ‘stall’ at separations of about a parsec for up to a Hubble time (e.g. Begelman, Blandford & Rees 1980; Milosavljević & Merritt 2001; Yu 2002). In gas-rich galaxy mergers, however, gas is driven to the central region of the merging galaxies simultaneously with the formation of the SMBH binary. Numerical simulations indicate that BH merger time-scales may be much shorter than galaxy merger time-scales in this case (∼106–107 yr from the hard binary stage; e.g. Escala et al. 2005; Dotti et al. 2007). This implies that the SMBH binaries most able to accrete gas and produce electromagnetic signatures are likely to be short-lived.
To date, observations seem to confirm this scenario. Several quasar pairs with kpc-scale separations have been found in merging galaxies (Komossa et al. 2003; Bianchi et al. 2008; Comerford et al. 2009b; Green et al. 2010). Recently, spectroscopic and photometric observations have found that between 10−3–10−2 of active galactic nuclei (AGNs) at z≲ 0.7 may in fact contain dual BHs with separations of ∼kpc (Comerford et al. 2009a; Fu et al. 2010; Liu et al. 2010a,b; Smith et al. 2010). However, only one unambiguous case of a SMBH binary has been found, with a separation of ∼7 pc (Rodriguez et al. 2006). Several additional objects have been identified as candidate SMBH binaries but are still unconfirmed (Sillanpaa et al. 1988; Komossa, Zhou & Lu 2008; Bogdanović, Eracleous & Sigurdsson 2009; Boroson & Lauer 2009; Dotti et al. 2009). Thus, it is clear that the large majority of AGNs do not contain binary SMBHs, which suggests that any substantial population of long-lived SMBH binaries must exist in gas-poor environments where they are quiescent.
If most binaries are therefore assumed to merge on reasonably short time-scales, then gravitational-wave (GW) recoil must be a common phenomenon throughout the merger history of galaxies. GW recoil is a natural consequence of gravitational radiation from BH binary mergers (Peres 1962; Bekenstein 1973; Fitchett & Detweiler 1984). If the binary system has any asymmetry – unequal masses, spins or spin orientations – then GWs are radiated asymmetrically, resulting in a net linear momentum flux from the final BH at the time of merger. This causes the BH to recoil in the opposite direction. Whether these recoil kicks were large enough to be astrophysically interesting was uncertain until a few years ago, because accurate calculations of the recoil velocity require BH merger simulations using full general relativity. These simulations have revealed that GW recoils may be quite large. Recoil velocities up to 4000 km s−1 are possible for special configurations, which is far greater than galactic escape speeds (Campanelli et al. 2007a,b). The implication that SMBHs may spend substantial time in motion and off-centre has opened a new line of inquiry into the ramifications for SMBHs and their host galaxies.
A useful starting point for such inquiries would be the distribution of recoil kick velocities, but this is quite difficult to ascertain in practice. The final velocity depends sensitively on not only the mass ratio of the progenitor BHs, but also their spin magnitudes and orientations. The distributions of SMBH binary mass ratios and spins at various redshifts have been estimated using halo merger trees and semi-analytical models of SMBH growth (e.g. Volonteri, Haardt & Madau 2003; Volonteri et al. 2005; Berti & Volonteri 2008; King, Pringle & Hofmann 2008). These distributions depend on a number of model assumptions, however, and the BH spin orientations prior to merger are far more uncertain still. Thus, based on recent results from BH merger simulations using full numerical relativity (NR), several groups have calculated kick probability distributions as a function of the BH mass ratio for either fixed or random values of BH spin, with the assumption that the spins are randomly oriented (Campanelli et al. 2007a; Schnittman & Buonanno 2007; Baker et al. 2008; Lousto et al. 2010a,b; van Meter et al. 2010). Their results are in good agreement with each other and imply that the fraction of high-velocity GW recoils is substantial. This underscores the potential for recoil events to be an important component of galaxy mergers.
It is quite possible, of course, that the spins of SMBH binaries are not randomly oriented, but are preferentially aligned in some way. Bogdanović, Reynolds & Miller (2007) have suggested that torques in a circumbinary gas disc may align the BH spins with the orbital axis of the disc. In this case, the resulting in-plane kicks would have a maximum recoil velocity of <200 km s−1, although spins that became anti-aligned by the same mechanism could result in recoil velocities up to 500 km s−1 (e.g. Campanelli et al. 2007a; González et al. 2007; Baker et al. 2008). Additionally, simulations by Dotti et al. (2010) demonstrate efficient spin alignment of merging BHs in the presence of a highly coherent accretion flow, although it is unclear how efficient this process might be in a gaseous environment that includes, for example, star formation. Kesden, Sperhake & Berti (2010) have recently demonstrated that BH spin alignment may instead occur via relativistic spin precession, regardless of whether a gas disc is present. The aforementioned recoil kick probability distributions are therefore upper limits on the actual distributions.
Numerous possible consequences of GW recoil events have been discussed in the literature. GW recoils may have a large effect at high redshifts, where escape velocities of galaxies are smaller (e.g. Madau & Quataert 2004; Merritt et al. 2004; Volonteri 2007). This is a concern for attempts to understand, for example, the origin of the z = 6 SDSS quasars (e.g. Fan et al. 2001, 2003). Volonteri & Rees (2006) suggest that growth of SMBHs at high z must occur only in highly biased haloes. Using cosmological hydrodynamic simulations, Sijacki, Springel & Haehnelt (2009) investigate BH growth in massive high-z haloes, including the ejection of BHs with recoil velocities above vesc. Their findings are consistent with the observed populations of bright quasars at z = 6, despite the effects of GW recoil.
At lower redshifts, recoiling BHs may produce electromagnetic signatures. The main signatures we will focus on here are recoiling AGNs that are either spatially or kinematically offset from their host galaxies (Madau & Quataert 2004; Loeb 2007; Blecha & Loeb 2008, hereinafter BL08; Komossa & Merritt 2008a; Guedes et al. 2009). Other possible signatures include flares from shocks induced by fallback of gas marginally bound to the ejected BH (Lippai, Frei & Haiman 2008; Schnittman & Krolik 2008; Shields & Bonning 2008), enhanced rates of stellar tidal disruptions (Komossa & Merritt 2008b; Stone & Loeb 2010) and compact stellar clusters around ejected BHs (Merritt, Schnittman & Komossa 2009; O’Leary & Loeb 2009).
Thus far, no confirmed GW recoil events have been observed. An inherent challenge in observing offset quasars is that larger spatial or kinematic offsets are easier to resolve, but less gas will be bound to the recoiling BH at higher recoil velocities, so its AGN lifetime will be shorter. Bonning, Shields & Salviander (2007) conducted a search for kinematic offsets in SDSS quasar spectra and found a null result at a limit of 800 km s−1. Several recoil candidates have been proposed, but their extreme inferred velocities should be exceedingly rare. Indeed, the recoil candidate with a 2600 km s−1 offset, proposed by Komossa et al. (2008), may in fact be a superposition of two galaxies (Heckman et al. 2009; Shields, Bonning & Salviander 2009a) or a binary SMBH system (Bogdanović et al. 2009; Dotti et al. 2009). Another candidate with an even higher (3500 km s−1) offset is most likely a double-peaked emitter (Shields et al. 2009b). Recently, Civano et al. (2010) have suggested that an unusual galaxy discovered in the COSMOS survey by Comerford et al. (2009b), which those authors proposed to be a dual SMBH system, may in fact be a recoiling BH, as new spectra indicate a kinematic offset of 1200 km s−1. This candidate has a less extreme velocity than the others, but further observations are needed to confirm the nature of this object. Additionally, the SMBH in M87 has recently been observed to be spatially offset by ∼7 pc, which could possibly be explained by a past recoil event (Batcheldor et al. 2010).
In addition to producing direct observational signatures, GW recoil may play a role in the co-evolution of SMBHs and their host galaxies. Strong empirical correlations exist between SMBH mass and properties of the host galaxy bulge, including the bulge luminosity, mass and stellar velocity dispersion (e.g. Kormendy & Richstone 1995; Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Merritt & Ferrarese 2001; Tremaine et al. 2002; Marconi & Hunt 2003). These correlations are well reproduced by galaxy formation models in which much of BH and galaxy bulge growth occurs via successive mergers and in which merger-triggered BH fuelling is self-regulated via AGN feedback (e.g. Silk & Rees 1998; Wyithe & Loeb 2003; Di Matteo, Springel & Hernquist 2005; Hopkins et al. 2006a). However, because GW recoil events may occur simultaneously with this rapid BH accretion phase, recoils could significantly disrupt the coordinated growth of BHs and galaxy bulges. In particular, GW recoil may contribute to scatter in the MBH–σ* relation caused by ejected (Volonteri 2007) or bound recoiling (BL08; Sijacki et al. 2010) BHs; we will examine the latter possibility in greater detail. It is also unclear a priori what effects, if any, GW recoil may have on the host galaxies themselves. In purely collisionless galaxies, Boylan-Kolchin, Ma & Quataert (2004) and Gualandris & Merritt (2008) have shown that bound recoiling BHs may scour out a stellar core. In gas-rich galaxy mergers, copious of amounts of cold gas are driven to the central galactic region during coalescence, triggering a luminous starburst such that the galaxy may appear as a ultraluminous infrared galaxy (ULIRG) (Sanders et al. 1988a; Sanders & Mirabel 1996). Feedback from a central AGN may terminate the starburst phase by expelling gas and dust from the central region (e.g. Hopkins et al. 2006a, 2008b; Somerville et al. 2008). We will investigate whether observable starburst properties or gas and stellar dynamics may be affected by the sudden displacement of this central AGN via recoils (see also recent work by Sijacki et al. 2010).
Our present work was initiated as a follow-up of BL08. BL08 explored the trajectories and accretion of BHs on bound orbits in a static potential with stellar bulge and gaseous disc components. They found that recoil kicks at ≲vesc could produce long-lived (≳Gyr) oscillations of the BH; similar results were found by Gualandris & Merritt (2008), Kornreich & Lovelace (2008) and Guedes et al. (2009). BL08 also demonstrated that GW recoil can affect SMBH growth even when the BH is not ejected entirely from the host galaxy. By removing accreting BHs from the central dense region and thereby limiting the BH’s fuel supply, BHs receiving recoil kicks at ≳0.5vesc may be less massive than stationary BHs.
In this study, we use hydrodynamic simulations to self-consistently calculate the dynamics and accretion of recoiling BHs in realistic merger remnant potentials. We simulate galaxy mergers using gadget-3 (Springel 2005), a smoothed particle hydrodynamic (SPH)/N-body code, and apply a recoil kick to the BH at the time of BH merger. Due to the detailed initial conditions and many free parameters involved in such merger simulations, we conduct a large parameter study with dozens of galaxy merger models and a wide range of kick velocities. Each merger model is simulated with at least one kick velocity and also with no recoil kick, for comparison. With this suite of recoil simulations, we are able to observe trends in the behaviour of recoiling BHs in different environments. When paired with our time-dependent, subresolution models for recoiling BH accretion, this approach allows us to, for example, determine velocity-dependent recoiling AGN lifetimes and estimate the effect of GW recoil on scatter in the MBH–σ* relation.
As this paper was in the final stages of preparation, two papers appeared that also involve hydrodynamic simulations of recoiling BHs in galaxy merger remnants (Guedes et al. 2010; Sijacki et al. 2010). However, each study used only a few simulations, probing only a small part of parameter space. Guedes et al. (2010) used results of three merger simulations as initial conditions and ran hydrodynamic simulations of recoiling BHs in these merger remnants for a short time. These initial trajectories were used to calibrate semi-analytic calculations of the recoil trajectory in the remnant potential. A range of kick velocities were tested in each of the three merger models used. Sijacki et al. (2010) simulate recoils in an isolated, stable disc galaxy; they also use this galaxy model as the initial condition for a full merger simulation with GW recoils. Both studies also include prescriptions for accretion on to recoiling BHs. We provide a comparison of our results to the findings of these papers in Section 6.
Our simulation methods, GW recoil treatment and galaxy merger models are outlined in Section 2. Sections 3–5 contain our results, organized as follows. In Section 3, we discuss the general characteristics of our galaxy merger simulations with stationary central BHs (i.e. no recoil kicks). We also discuss the variation in merger dynamics and remnant morphologies between models. In Section 4, the dynamics of recoiling BHs are discussed. We describe the general characteristics of recoil trajectories and trends between models in Section 4.1. In Section 4.2, we analyse the sensitivity of recoil trajectories to the BH merger time and in Section 4.3, we investigate the dependence of recoil trajectories on the direction of the kick. Section 5 contains our results for the accretion and feedback of recoiling BHs. In Section 5.1, we examine accretion and feedback of recoiling BHs, compare to the fuelling of stationary BHs and calculate AGN lifetimes from Bondi–Hoyle accretion. In Section 5.2, we describe our time-dependent analytic model for calculating accretion on to a recoiling BH from a disc of gas ejected with the BH. Using this accretion model, along with accretion from ambient gas, we present in Section 5.3 the velocity-dependent AGN luminosities and active lifetimes for recoiling BHs. In particular, we calculate lifetimes for spatially offset and kinematically offset recoiling AGNs. In Section 5.4, we compare the MBH–σ* relations derived from our set of no-recoil simulations and our set of high-velocity recoil simulations, and discuss implications for the observed MBH–σ* relation. Finally, in Section 5.5, we explore the effects of GW recoil on star formation rates (SFRs) and the host galaxy structure. We summarize and discuss our results in Section 6.
For our simulations of galaxy mergers, we employ the SPH code gadget-3 (Springel 2005), which conserves both energy and entropy (Springel & Hernquist 2002). The code includes radiative cooling as well as a subresolution model for a multiphase interstellar medium (Springel & Hernquist 2003) that accounts for star formation and supernova feedback. In addition, the code models BHs as gravitational ‘sink’ particles that contain a BH seed and a gas reservoir. The reservoir is replenished by stochastic accretion of neighbouring gas particles, but the actual accretion rate on to the BH is calculated smoothly using the Bondi–Hoyle–Lyttleton formula with locally averaged values for the density and sound speed. Because the gas around the BH cannot be resolved at higher densities below the spatial resolution limit, the accretion rate calculated from these values is multiplied by a constant factor. Following other authors, we assume a value of 100 for this factor (e.g. Springel, Di Matteo & Hernquist 2005b; Hopkins et al. 2006b; Johansson, Naab & Burkert 2009). Angular momentum is conserved during accretion of gas particles, but because this is a stochastic process owing to the finite mass resolution, we also introduce an analytic accretion drag force () calculated from the Bondi accretion rate at each time-step. These prescriptions are described in more detail in Springel et al. (2005b).
We note that, as in all hydrodynamic simulations involving BH accretion, the exact nature of the accretion on subresolution scales is unknown, so it is necessary to assume a subresolution model. While the validity of the Bondi–Hoyle–Lyttleton formula (with a multiplicative factor) in our framework is indeed an assumption, few constraints exist on the gas flow from large scales down to the BH. The Bondi–Hoyle–Lyttleton model provides a physically motivated approximation of the gas captured by a BH, including a BH in motion, and its extensive use in the literature helps to place our work in a larger context of galaxy evolution studies. We further assume that the accretion factor remains constant even after a GW recoil event. The latter assumption is justified in that the spatial resolution limit, which motivates the use of this factor, always remains constant. Furthermore, owing to the sink-particle treatment of BHs in the code, the accretion factor does not determine the amount of gas bound to the BH at each time-step; it influences only how this gas is accreted on to the BH. From a practical standpoint, using the same accretion prescription throughout enables a more direct comparison between recoiling and stationary BHs – one of the main goals of this study.
Because we are presently interested in GW recoil, we allow for an arbitrary velocity to be added to the remnant BH at the time of the BH merger, tmrg. To more easily compare recoil events in different merger models, we scale this velocity to the central escape speed at the time of merger, vesc(tmrg). We define at each time-step, that is, the escape speed at the position of the BH, which is close to the position of the potential minimum in the absence of a recoil kick. Campanelli et al. (2007a), Baker et al. (2008) and van Meter et al. (2010) have all fitted distributions to the data from full NR simulations of BH mergers, with results in good agreement. For major mergers (0.25 < q < 1, the range we consider in our simulations), high spins (a1=a2 = 0.9) and randomly distributed spin orientations, van Meter et al. (2010) find that 69.8 per cent of recoil kicks will be above 500 km s−1 and 35.3 per cent will be above 1000 km s−1. If the BH spin magnitudes are also randomly distributed (0 ≤a1, a2≤ 1), then these fractions are 41.6 and 13.5 per cent, respectively. Our fiducial-mass merger simulations (described below) have vesc = 689–1248 km s−1 and we also use low-mass simulations with vesc as low as 338 km s−1. The NR results indicate that there is a substantial probability of BHs in these models receiving kicks up to ∼vesc. In our high-mass simulations, however, vesc = 1737–2555 km s−1, so the probability of recoil events with vk∼vesc is much lower in these models.
We also slightly modify the treatment of BH mergers in the code. In the standard gadget prescription, the BHs merge when they are separated by less than a gravitational softening length (asep < Rsoft) and have a relative velocity vrel < 0.5cs, where cs is the local sound speed. If, as can happen in gas-rich mergers, the central escape velocity is increasing rapidly when the BHs merge, small variations in tmrg can result in very different BH trajectories (see Fig. 8 shown later). In order to reliably compare results from different simulations, we therefore force the BHs to merge at a given time that is predetermined as follows. We define the ‘coalescence time’tcoal as the earliest time after which the BH binary is tight enough that it could plausibly merge. Given the numerical uncertainties near the spatial resolution limit, we generously choose tcoal to be the time when the BH separation falls below asep = 10Rsoft (and vrel < 0.5cs, as before). For reference, Rsoft = 80 pc in our fiducial-mass simulations (see below). Then, restarting the simulation slightly before tcoal, we force the BHs to merge at a predetermined time tmrg≥tcoal. tmrg can then be varied as a free parameter, allowing us to systematically probe the range of possible BH merger times. In practice, this is only necessary for nearly-equal-mass, gas-rich mergers in which vesc varies significantly throughout the merger. In other cases, the results are insensitive to the choice of merger time, so tmrg is simply chosen to be tcoal.
For the majority of our simulations, we use the same total mass, 1.36 × 1012 M⊙, for the primary galaxy and scale the secondary to yield the desired mass ratio. The exception is a small subset of simulations with lower and higher total mass, described below. For all of our ‘fiducial-mass’ simulations, we use a gravitational softening length of Rsoft = 80 pc for baryons and Rsoft, DM = 240 pc for the dark matter (DM). In each galaxy, 4.1 per cent of the total mass is in a baryonic disc component. The primary galaxy has Nhalo = 4.8 × 105 DM halo particles and for fgas = 0.4, it has Ndisc = 3.2 × 105 disc particles, with equal numbers of gas and star particles. In all other models, Nhalo and Ndisc are set such that the same particle mass resolution is maintained for each component (mstar = 4.2 × 105 M⊙, mgas = 2.8 × 105 M⊙ and mhalo = 5.4 × 106 M⊙).
Each galaxy is given a single BH particle, which, as previously stated, consists of a ‘seed’ BH along with a gas reservoir from which an accretion rate is calculated smoothly. We use a small initial seed mass for the BH, 1.43 × 105 M⊙, in accordance with the MBH–Mbulge (Magorrian et al. 1998) relation, because our initial galaxies are pure discs with no bulge component. We find that our results are insensitive to the choice of seed mass. The total mass of this hybrid particle is the ‘dynamical’ mass of the BH; this is the mass that is used for gravitational force calculations. This dynamical BH mass is set to 10−5 of the total galaxy mass initially. This value is chosen such that the seed mass and dynamical mass have similar values by the time of the recoil kick, to avoid a disparity between the mass used for gravitational forces and that used for accretion physics of the recoiling BH. To help ensure that the BH remains in the centre of the galaxy prior to the BH merger and recoil kick, we set the accretion drag force equal to its value for Eddington-limited accretion prior to the BH merger (though the actual accretion rate is still calculated self-consistently based on the Bondi–Hoyle formula). Once the BHs merge, the accretion drag force is calculated using the actual Bondi–Hoyle accretion rate for the remainder of the simulation.
Our set of fiducial-mass simulations covers a wide range of galaxy mass ratios and gas fractions, but in order to probe the effects of GW recoil on black hole–host galaxy relationships, as we do in Section 5.4, we require a larger range in total galaxy mass. We have run a sample of seven high- and low-mass simulations with a primary galaxy mass ranging from 0.05 to 20 times the fiducial primary mass, which we refer to as M0. For the lower mass runs, we maintain the same number of particles as in the fiducial simulations, such that higher mass resolution is achieved. Accordingly, we reduce the softening lengths in these simulations by (M/M0)1/3. For the higher mass runs (M0× 10 and M0× 20), an increase in the particle number by this factor would be prohibitively computationally expensive, so we instead increase the particle number by a factor of 5 in each case, compromising a factor of 2 and 4 in mass resolution, respectively. In the latter case, we maintain a reasonable MBH/Mparticle ratio by increasing the initial dynamical mass by a factor of 10.
In order to evaluate the sensitivity of our results to numerical artefacts, such as the choice of mass resolution, softening length and integration accuracy, we have run a number of additional simulations in which we systematically vary these parameters. In addition to these tests, two of the models included in our results have higher mass resolution (and correspondingly higher spatial resolution) than our fiducial runs: our low-mass merger models have the same number of particles as our fiducial-mass models and thus have 10 and 20 times higher mass resolution. Higher mass and spatial resolution reduce the noisiness of the potential, and varying the resolution changes the detailed structure of the gas, especially in the highest density regions. Higher integration accuracy also reduces numerical noise. Consequently, the exact trajectory, and hence the settling time, of a given recoil event will be affected to some extent by these numerical factors. However, the variation in recoil trajectories for different choices of resolution and integration accuracy are small compared to the differences for varying kick speed and merger remnant properties. We also define the settling radius of the BH to be 4Rsoft to avoid sensitivity of our conclusions to the BH motion in the innermost galaxy region where the softened potential may dominate. Most importantly, for all choices of mass resolution, spatial resolution and integration accuracy parameters tested, the qualitative behaviour of the recoiling BHs and the relative differences between simulations are robust. Therefore, our main results do not depend on these numerical factors.
Table 1 lists the parameters for the galaxy merger models we have constructed. We have tested a total of 62 different galaxy merger models in which we vary the galaxy mass ratio (q), the gas fraction (fgas) and the orbital configuration. 55 of these models use the ‘fiducial’ primary galaxy mass and for easier reference, we assign each of these a name given by q[value]fg[value][orb], where ‘q’ denotes the galaxy mass ratio, ‘fg’ denotes the initial gas fraction and each letter orb is identified with a specific orbital configuration. The remaining seven models have higher or lower total mass and are referred to by q[value]fg[value]M[factor][orb], where ‘M’ denotes the total galaxy mass relative to the equivalent fiducial-mass galaxy.
|Model||Initial model parameters||Resulting merger quantities|
|q||Mtot (1010 M⊙)||Mdisc (1010 M⊙)||fgas||Rperi (kpc)||θ1 (°)||ϕ1 (°)||θ2 (°)||ϕ2 (°)||tmrg (Gyr)||vesc (tmrg) (km s−1)||vesc, max (km s−1)||Per cent difference|
|Model||Initial model parameters||Resulting merger quantities|
|q||Mtot (1010 M⊙)||Mdisc (1010 M⊙)||fgas||Rperi (kpc)||θ1 (°)||ϕ1 (°)||θ2 (°)||ϕ2 (°)||tmrg (Gyr)||vesc (tmrg) (km s−1)||vesc, max (km s−1)||Per cent difference|
In all cases, the orbital configuration is specified by six parameters and the galaxies are initially set on parabolic orbits. The choice of the impact parameter for the first pericentric passage, Rperi, and initial galaxy separation, ai, determine the initial orbital energy and angular momentum. In our fiducial-mass models, we set ai = 143 kpc and Rperi = 7.1 kpc, except for orbit ‘b’, which has Rperi = 14.3 kpc. In our high- and low-mass models, ai and Rperi are scaled such that ai is the same fraction (0.625) of the virial radius, R200, and ai/Rperi = 20. We denote these orbits with the letters w–z to differentiate them from the orbits a–k of our fiducial-mass mergers. Note that the orbital parameters do not remain constant as the merger progresses, owing to energy losses via dynamical friction. The angles (θ1, ϕ1) and (θ2, ϕ2) determine the initial orbital phase and inclination, respectively, of each galaxy. These angles, as well as Rperi for each galaxy model, are given in Table 1. Our fiducial merger orbit (orbit ‘a’) is a ‘generic’ orbit, in that no symmetries exist in the initial orientation angles of the discs. Note that many of these orbits are identical to those used in Cox et al. (2006) and Robertson et al. (2006a).
For each model, we run a merger simulation in which the BHs are not allowed to merge. From this we can determine tcoal, which is used to set tmrg. Then, restarting slightly before tmrg, we simulate both a merger with no recoil kick and a merger with vk/vesc = 0.9. For the models shown in boldface in Table 1, we also simulate recoil kicks with vk/vesc = 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.1 and 1.2. Each simulation is run for 2.9 Gyr after tmrg. Finally, using the q0.5fg0.3a model, we simulate a small sample of recoil kicks with varying kick orientation (θk, ϕk), the results of which are discussed in Section 4.3.
3 MERGERS WITHOUT GRAVITATIONALWAVE RECOIL
We begin by examining galaxy merger simulations with no GW recoil kick given to the BHs. Some basic characteristics of a galaxy merger simulation with q = 0.5 and fgas = 0.3 (our ‘fiducial’ merger, model q0.5fg0.3a) are illustrated in Figs 1 and 2. Fig. 1 shows the gas distribution at the time of BH merger, as well as the paths of the two BHs prior to merger. The gas distribution is very irregular, with distinct tidal streams. The inner region of the merger remnant is also lumpy and irregular. In other words, the initial disc structure of the progenitors has been destroyed by the merger and the remnant is still highly disturbed at the time of BH coalescence. Fig. 2 shows the BH accretion and star formation throughout the simulation as well as the evolution of the BH separation prior to merger. This example illustrates some characteristics that are generic to our galaxy merger simulations. Simultaneous bursts of star formation and BH accretion occur after the first close pericentric passage and at final coalescence (Mihos & Hernquist 1994b, 1996) as the result of gas inflows caused by gravitational torques (Barnes & Hernquist 1991, 1996). As we shall demonstrate, these gas inflows can greatly influence the dynamics of ejected black holes by deepening the central potential and by providing additional drag force. Note that because our simulations include star formation throughout, the gas fraction at the time of BH merger is significantly lower than its initial value; in the example shown, 60 per cent of the initial gas has been consumed by tmrg. By the end of the simulation, 2.9 Gyr after the merger, almost 80 per cent of the gas has been depleted by star formation and the BH accretion rate is very low. In our other merger models, 47–80 per cent of the initial gas is depleted by the time of BH merger and 62–88 per cent is depleted by the end of the simulation. The final black hole mass in this example (2.2 × 107 M⊙) is also typical of our fiducial-mass mergers.
Fig. 3 shows gas and stellar density distributions of this fiducial (q0.5fg0.3a) model and five other merger models. Each plot shows gas and stellar density from three different projections. These examples are chosen to illustrate the generic nature of the morphological features mentioned above. In all cases, the merger remnants are visibly disturbed and lumpy. Tidal tails are ubiquitous, though their size and density vary. Higher fgas mergers have more compact remnants (Dekel & Cox 2006; Robertson et al. 2006b; Hopkins, Cox & Hernquist 2008a; Hopkins & Hernquist 2010a). Lower q mergers are less strongly disrupted; a disc-like structure can be seen in the x–z orientation of model q0.25fg0.3a. A special case is shown in the lower right-hand corner of Fig. 3. Model q0.5fg0.3h has a coplanar orbit with both galaxies rotating prograde to the orbital motion. Due to the highly aligned orbit, the disc structure of the progenitor galaxies is preserved in this merger remnant; tidal features are visible only from the face-on (x–y) projection. Our other ‘highly aligned’ orbits (i, j and k) result in similar remnant morphologies; they are included in our suite of merger models mainly for comparison, as mergers generally do not have such careful alignment.
As mentioned in Section 2, in nearly-equal-mass, gas-rich mergers, the potential well in which the BH sits may deepen rapidly during final coalescence. Once the merger is complete, the remnant central potential will become shallower as the central stellar region begins to relax. The phase of rapid vesc fluctuation coincides with the time of BH coalescence and may therefore affect recoiling BH dynamics. Fig. 4 shows the evolution of the BH escape speed (vesc) for four merger simulations without recoil kicks. The solid lines denote tmrg and vesc(tmrg), while the dotted lines denote the maximum (post-merger) value of vesc and the time at which it occurs. The difference between vesc(tmrg) and vesc, max is small in all except the top left-hand panel, which shows the q1fg0.4a merger. It is clear that the escape speed, and hence the trajectory, of a kicked BH in this model will depend on whether the kick occurs at tcoal or at some slightly later time. However, the sharp increase in vesc occurs only in gas-rich, nearly-equal-mass mergers. The other three panels in Fig. 4 show vesc versus t for a sample of simulations with lower q and fgas; in these models, vesc is much less volatile during and after the BH merger. Moreover, the example shown in Fig. 4 is not the most gas rich of our merger models but is chosen because of the large difference between vesc(tcoal) and vesc, max, 34 per cent. Larger differences between these quantities yield larger uncertainty in our results for BH dynamics when the recoil kick is assigned at tmrg=tcoal. The final four columns in Table 1 give tmrg, vesc(tmrg), vesc, max, and the fractional difference between vesc(tmrg) and vesc, max for each merger model. All mergers with fgas < 0.3 have ≲10 per cent difference between vesc(tmrg) and vesc, max, and mergers with q < 0.5 have <1 per cent difference. Therefore, while our results for nearly-equal-mass, gas-rich mergers are subject to the assumption that the BH merger occurs on a short time-scale, the results for our other models are insensitive to the merger time. We explore further the case of gas-rich, q∼ 1 mergers in Section 4.2.
4 RECOILING BLACK HOLE DYNAMICS
4.1 Characterization of recoil trajectories
A universal feature of our recoil trajectories is that they have low-angular-momentum orbits, which occurs because the centrally concentrated baryonic component of the galaxy dominates the BH trajectories even when they extend far into the halo. Consequently, we refer, throughout this paper, to the trajectories as ‘oscillations’ of the BH about the galactic centre. We see this clearly in Fig. 5, which shows the trajectories of BHs kicked with vk/vesc = 0.4–1.0 for eight different merger models. Also readily apparent in these plots is the variation in BH oscillation amplitude and duration between models. For a fixed value of vk/vesc, BHs travel farther from the galactic centre in mergers with lower q and fgas. In the q1fg0.6a model, BHs kicked with vk = 0.7vesc ( = 878 km s−1) travel <2 kpc from the galactic centre, while in the collisionless q1fg0a model, the same vk/vesc allows the BH to travel 10 times farther. Note that scaling the kick speeds to vesc is important for determining the trends between galaxy models; vesc(tmrg) for all of our models ranges from 338 to 2555 km s−1(689 to 1258 km s−1 for the fiducial-mass models). The trend towards smaller recoil oscillation amplitude for higher fgas occurs because higher gas fractions result in more compact remnants that have higher central densities and these remnants also have larger available supplies of gas. Both factors contribute to steeper central potentials and increased gas drag and dynamical friction. We can also see from Fig. 5 that recoil oscillations are slightly larger for lower q, though the trend is more evident in higher fgas remnants. Higher mass ratio mergers create stronger perturbations that drive gas more efficiently to the galaxies’ central regions, further contributing to the steep potentials and high densities.
We can see these trends more clearly in Fig. 6. The top plot shows the maximum galactocentric distance of each orbit in Fig. 5, normalized to the half-mass effective radius of the galaxy (Rmax/Reff). The bottom plot shows the settling time, tsettle, defined as the time when the apocentric distance of the BH oscillations, relative to the stellar centre of mass, falls below 4Rsoft. This definition is chosen to avoid following BH trajectories below a scale where their oscillations may be dominated by numerical noise or gravitational softening. We find that our conclusions are not sensitive to the exact definition of tsettle. Note that for some high-velocity recoils (for all recoils in the collisionless simulations), we have only lower limits on tsettle, because in these cases the BH did not settle by the end of the simulation.
In the top panel of Fig. 6, the curves for a given fgas lie nearly on top of each other. This supports our assertion that the initial gas fraction is more important than the mass ratio in determining Rmax/Reff for a given vk/vesc. By normalizing these quantities, we essentially isolate the effects of gas drag or potential shape from variation in galaxy size or central potential depth. One can easily see from Fig. 6 that, especially for lower vk/vesc, a gas fraction ≥30 per cent can reduce the amplitude of BH oscillations by up to an order of magnitude with respect to a purely collisionless system.
Note that Rmax/Reff is at least as large for fgas = 0.6 as it is for fgas = 0.3; opposite the trend we expect, if gas drag dominates the BH trajectories. However, galaxies with higher fgas generally have higher SFRs, so systems with initial gas fractions of 0.3–0.6 will have similar amounts of gas by the time of the BH merger. Furthermore, as previously stated, higher gas fractions yield more concentrated merger remnants (smaller Reff). Combined, these effects result in similar Rmax/Re for recoils in mergers with initial fgas = 0.3–0.6.
As discussed in Section 2, the recoil trajectories may be sensitive to numerical parameters, such as resolution or integration accuracy. In particular, because Reff is of the order of a few kpc, recoiling BHs with Rmax/Reff≲ 0.1 spend substantial time at or below the spatial resolution limit. Rmax and tsettle are therefore approximate in these cases. The BH settling times are also susceptible to possible inaccuracies in dynamical friction forces due to our finite mass resolution, which is likely to result in underprediction of the true amount of dynamical friction. However, we can predict that the effect of this uncertainty on our results is likely to be small; in our lower mass galaxy mergers, which have 10 and 20 times better mass resolution, BHs kicked near the escape speed still experience very little dynamical friction and have very long wandering times. Even if dynamical friction is still underestimated in these higher resolution runs, it is clear that, in general, recoils near the escape speed will result in wandering times of at least a few Gyr and that in gaseous mergers, recoiling BHs with vk≲ 0.5vesc remain within the central kpc of the galaxy. Finally, because we use the same mass resolution in all of our simulations (excepting the low-/high-mass runs), the relative trends that we find between different galaxy models are robust.
We expect some variation in Rmax for different merger orbits as well. In the set of 11 different orbital configurations we simulated with q = 0.5 and fgas = 0.3, the central escape speeds at tmrg vary from 805 to 1046 km s−1. Since the total mass in all of these simulations is constant, the variation in vesc reflects differences in the central concentration and hence the steepness of the potential well. We simulate recoil trajectories in each of these with vk/vesc = 0.9 and find that Rmax/Reff = 4–20. Based on this example, we see that the merger orbit alone can influence the remnant central potential well enough to change the amplitude of recoil trajectories by a factor of ∼5.
Fig. 7 shows two examples of recoil trajectories in three different projections as well as R(t). The top plot shows a recoil event with vk/vesc = 0.6, which in this model (q0.5fg0.3a) corresponds to vk≈ 600 km s−1. Note that the orbit is highly non-planar; it looks similar in all three projections. The orbit is also centrophilic, which as we have noted is common to all of our recoil simulations. The average orbital period is short, ∼7 × 105 yr, and tsettle∼ 400 Myr, such that the BH settles to the centre well before the end of the simulation.
Recoil events with velocities near vesc produce long-lived (>Gyr), large-amplitude BH oscillations. Many of these oscillations still have large amplitudes at the end of the simulation, 2.9 Gyr after the BH merger, so we have only lower limits on tsettle. The bottom plot in Fig. 7 illustrates this type of recoil event, again for the q0.5fg0.3a model. The trajectory shown has a kick speed of vk/vesc = 0.9, so vk≈ 900 km s−1. The amplitude of these oscillations is ∼50 kpc, almost two orders of magnitude greater than in the lower vk example. Although the BH only completes five orbits before the end of the simulation, we can see that again there is no evidence for a preferred orbital plane; the orbit is three-dimensional. Because the oscillation amplitude is so large and shows no sign of damping out by the end of the simulation, it is likely that such an orbit will not return to the galactic centre within a Hubble time and certainly not within the mean time between major galaxy mergers. Such BHs essentially can be considered ‘lost’ to the galaxy. If in other cases the BHs do return to the galactic centre after several Gyr and another galaxy merger has taken place in the meantime, then the returning and new BHs could form another binary.
In recoil events with v > vesc, the rapidly recoiling BH can only take a small amount of gas with it; only gas that is orbiting the BH with vorb≳vk stays bound to the BH when it is ejected. The accretion time-scale and luminosity of these ejected discs are discussed in Section 5.2. Once this supply is exhausted, the naked BH is unlikely to become active again, as this would require a serendipitous passage through a dense gas cloud, as well as a sufficiently low relative velocity that substantial accretion could occur. (Recall that Bondi–Hoyle accretion scales as v−3BH.) These escaping BHs are therefore likely to wander indefinitely, undetected, through intergalactic space.
4.2 BH merger time
As illustrated in Fig. 4, vesc can increase significantly during coalescence in equal-mass, gas-rich mergers. To better understand how this may affect BH trajectories, we have run a series of simulations with the BH merger (and recoil kick) occurring at sequentially later times using model q1fg0.4a. The vesc evolution of model q1fg0.4a is replotted in the top panel of Fig. 8, with the fiducial merger time, tmrg=tcoal, marked by the solid black line and the subsequent merger times tested marked by black dashed lines. We have tested simulations with each of these merger times twice, once with vk = 1000 km s−1 and then with vk/vesc(tmrg) = 0.9. By conducting the second experiment, we are normalizing to the evolving depth of the central potential well, which helps us to separate the gravitational and drag effects on the trajectories.
Fig. 8 shows Rmax for the BH orbits in each of these varying-tmrg simulations, with the vk = 1000 km s−1 simulations in the middle panel and the vk/vesc(tmrg) = 0.9 simulations in the bottom panel. As vesc increases from 1100 km s−1 at tcoal to almost 1500 km s−1 at its maximum, the amplitudes of BH trajectories from constant-vk recoil events decrease by more than two orders of magnitude. However, the bottom panel of Fig. 8 shows that when vk is set to a constant fraction of the escape speed, the BHs kicked from deeper potential wells (larger vesc) still have smaller amplitude trajectories, but the variation is much less severe than when vk = constant. The variation that persists even for fixed vk/vesc can be attributed to the enhanced gas drag and dynamical friction, as well as the steeper central potential well, that result from the rapid formation of the central dense cusp. At very late merger times, 100–500 Myr after tcoal, the opposite trend is seen. The amplitude of constant-vk trajectories increases by almost an order of magnitude from the minimum at tcoal+ 70 Myr to the value at tcoal+ 500 Myr. This is partly because the central potential becomes shallower as the merger remnant begins to relax. The same trend occurs at late times in the constant-vk/vesc simulations, because the potential well also becomes less steep, and gas drag becomes less efficient as more of the gas is consumed in star formation. We note that merger times up to tcoal+ 500 Myr are included here only for illustrative purposes; in reality, such delayed mergers are probably unrealistic in gas-rich systems.
If the BHs are able to merge rapidly, before the central potential reaches its maximum depth, another interesting effect can occur. During the period of rapid central potential evolution, the time-scale on which vesc increases is typically much shorter than the time-scale for a recoiling BH orbit to decay via gas drag or dynamical friction. Consequently, on its first several pericentric passages through the central region, the BH encounters an increasingly deep potential well, and its velocity is actually boosted by a small amount, such that the ratio vBH/vesc is roughly constant on subsequent pericentric passages. Because only vBH increases, not vBH/vesc, the BH’s galactocentric distance is not boosted beyond its initial maximum, Rmax, and once the rapid potential evolution ceases, drag forces become dominant and the velocity begins to decay. This effect is by design short-lived and it is strongest when the BH orbital time-scale is short, such that the BH orbit is dominated by the dynamics of the inner galaxy region. These velocity boosts are therefore mostly an interesting illustration of the sensitivity of recoil dynamics to the BH merger time, although we will see some implications of this effect for the recoiling AGN lifetimes calculated in Section 5.3.
The primary implication of our BH merger time analysis is that in nearly-equal-mass, gas-rich mergers, such as the example shown in Fig. 8, a delay of even a few times 107 yr in the BH merger may significantly reduce the amplitude and duration of recoiling BH oscillations. Based on results from numerous hydrodynamic simulations of binary BH inspirals (e.g. Escala et al. 2005; Dotti et al. 2007), we have good reason to believe that BHs in gas-rich systems merge efficiently in general, but the merger time of a given BH binary is impossible to predict with this precision. By assuming tmrg=tcoal in our simulations, we exclude the possibility of delayed mergers, but we also take steps to ensure that this has a minimal effect on our results. First, only a small number of our simulations have large variation in vesc. In these cases, we see by comparing Figs 6 and 8 that the variation in oscillation amplitudes for a given vk/vesc due to evolving central density (less than a factor of 10) is generally much less than the variation for different values of vk/vesc (up to three orders of magnitude). All of the analysis in this study is presented in terms of vk/vesc. Therefore, the dramatic suppression of recoil trajectories with fixed vk, while an important result in itself, does not affect our other conclusions.
4.3 Recoil kick direction
Most of our recoiling BHs have vk≳ 500 km s−1, in which case their kicks will be oriented out of the orbital plane of the binary (Campanelli et al. 2007a,b; Lousto & Zlochower 2009). However, because the BH binary inspiral occurs on subresolution scales in our simulations, the binary orbital angular momentum vector cannot be predicted. Furthermore, the galactic disc structure of the merger progenitors is disrupted by major mergers (see Fig. 3), except when the galaxies have nearly coplanar orbits. In generic mergers, there is no obvious reference direction in which to orient the recoil kick with respect to the progenitors’ initial structure. Clearly, on certain recoil trajectories (i.e. along tidal tails or central clumps), a BH will encounter denser gaseous regions than on others. However, it is unclear a priori whether these density variations are sufficient to cause significant variation in BH recoil trajectories. In most of our simulations, we simply orient the recoil kick along the z-axis of the simulation, (θk, ϕk) = (0, 0), in the coordinate system with respect to which the initial galaxy orbits are assigned. Here, we examine a subset of simulations with varying kick orientation to determine the sensitivity of our results to this choice.
We test a sample of different kick orientations with three different values of vk/vesc (0.6, 0.8 and 0.9) in our q1fg0.3a and q0.5fg0.3a merger models. As a gauge of how the BH trajectories differ in this sample, we calculate the maximum galactocentric distance (Rmax) and average orbital period (〈Porb〉) of the BH in each simulation. Fig. 9 shows the distributions of these quantities, normalized to Rmax and 〈Porb〉 of the (θk, ϕk) = (0, 0) simulations for each corresponding vk/vesc value. The q0.5fg0.3a distributions peak slightly below 1.0 and the q1fg0.3a distributions peak slightly above 1.0, but the combined histogram is distributed fairly evenly around unity. In other words, the (θk, ϕk) = (0, 0) kick orientation is not a ‘special’ direction. The overall variation in Rmax/Rmax(0, 0) for each model is less than a factor of 4. In contrast, Fig. 5 shows that Rmax is almost two orders of magnitude higher for vk/vesc = 0.9 than for 0.6. Clearly, the amplitude and duration of BH oscillations depend much more strongly on the kick magnitude than its direction and we need not concern ourselves too much with our choice of kick orientation.
5 BLACK HOLE ACCRETION AND FEEDBACK
5.1 Bondi accretion and AGN lifetimes
Fig. 10 compares the BH accretion and star formation histories for simulations of the q0.5fg0.3a merger model with vk = 0 and vk/vesc = 0.5–0.8. The top panels in Fig. 10 show that even kicks with vk/vesc = 0.5–0.6 produce qualitatively different AGN light curves (as inferred from ) from those in the zero-kick case, although the final BH masses are similar. In these cases, the peaks in accretion rate at coalescence (t≈ 2 Gyr) are slightly lower due to the sudden removal of the BHs from the highest-density region. However, the BHs remain within 1 kpc of the galactic centre and accretion continues through the end of the simulation at higher rates than in the zero-kick case. Thus, although the total mass accreted is about the same for vk/vesc = 0–0.6, the active lifetime is longer when a recoil kick occurs.
In the vk/vesc = 0.7 simulation, the BH accretion cuts off more sharply at tmrg and is subsequently highly variable, due to the larger amplitude, longer lived BH oscillations. Unlike the BHs with lower velocity recoils, this BH attains a notably lower final mass than its zero-kick counterpart.
For recoil kicks ≳0.8vesc, the BH accretion cuts off sharply at tmrg and the BH subsequently accretes very little gas as it travels on large, long-period orbits that extend well into the halo. These high-velocity recoil events therefore reduce the AGN lifetime. Additionally, the BH mass is essentially frozen at its value at tmrg, which is about half that of the stationary BH mass. If the BH does return to the galactic centre after several Gyr, most of the galaxy’s gas will be consumed in star formation, so the BH will remain undermassive. In the galaxy merger models we have used, the BH mass deficits for high-velocity recoils versus no recoil range from a factor of ∼1 (almost no deficit) to ∼5. These BH mass deficits and their implications are discussed in Section 5.4.
To quantify the amount of enhancement or reduction in AGN activity owing to GW recoil, we calculate AGN lifetimes (tAGN) for both stationary and recoiling BHs in six different galaxy models for which we have varied vk. Fig. 11 shows the results, with coloured lines indicating tAGN for recoiling BHs as a function of vk and black horizontal lines indicating tAGN for a stationary BH in the corresponding merger model. In order to differentiate between bright AGNs and low-luminosity AGNs, we calculate tAGN, assuming three different definitions of an ‘active’ BH. The thick solid lines in the figure are calculated assuming the BH is an AGN when Lbol > 10 per cent LEdd and the thin solid lines assume Lbol > 3 per cent LEdd. Because these quantities depend on the BH mass, we use an absolute luminosity value for our third AGN definition: Lbol > Lmin (denoted by the dashed lines in Fig. 11). We choose Lmin = 3.3 × 109 L⊙, which is 0.1LEdd for a 106-M⊙ BH.
In these models, the low-luminosity AGN lifetime is generally enhanced by moderate-velocity recoil events and the bright AGN lifetime is generally reduced. Although the details depend on the model parameters and kick speed, in all models, tAGN falls off steeply at high vk/vesc, when the BH is ejected far enough from the centre on large enough orbits that Bondi accretion becomes ineffective. This corresponds to the accretion behaviour seen in the last column of Fig. 10 (vk/vesc = 0.8).
For smaller vk, such that the BH is confined to the central few kpc, the dependence of tAGN on the nature of the merger remnant is more complex. We examine first the relatively gas-poor models (those with initial fgas = 0.1). In the q1fg0.1a model, tAGN is enhanced only for vk/vesc≤ 0.5 and falls precipitously at higher vk. In the q0.5fg0.1a model, recoil events always reduce the AGN lifetime. Merger remnants with less gas have smaller central gas reservoirs to fuel the BH during pericentric passages and their shallower central potentials result in larger BH orbits, which have fewer pericentric passages. Consequently, in gas-poor remnants, tAGN will generally be lower for recoiling BHs than for stationary BHs.
In gas-rich merger remnants, tAGN is more readily enhanced by GW recoil. Stationary BHs in gas-rich mergers typically experience a bright AGN phase due to the rapid inflow of cold gas to the central region, which is terminated when the AGN feedback energy heats the surrounding gas enough to quench the BH fuelling. Fig. 11 shows that although GW recoil shortens this bright AGN phase, the low-luminosity tAGN is enhanced for recoil kicks as high as vk/vesc = 0.7. In these cases, the recoiling BHs are on tightly bound orbits and can accrete gas during their numerous passages through the central dense region. Because they remain in motion, their feedback energy is spread over a much greater volume than for stationary BHs. As a result, the gas encountered by the recoiling BHs never heats up enough to completely cut off the BH fuel supply and the BH has a longer lifetime as a low-luminosity AGN. In contrast to the traditional ‘feast or famine’ model of merger-triggered AGN fuelling, these short-period recoiling BHs are more inclined to ‘nibble’.
We can see evidence of this displaced AGN feedback through its effect on star formation. The bottom row of Fig. 10 shows the depletion of gas throughout each simulation, normalized to the initial gas fraction, and the dotted lines show the vk = 0 curve in each of the other panels. It is evident that slightly less gas remains at the end of each higher vk simulation. In the merger with a stationary BH, this gas is accreted or expelled from the galactic centre via feedback, but when the central BH is removed, the gas continues to be consumed by star formation. Although in this example the SFR enhancement is small and difficult to discern from examining the SFR evolution alone, the effect can be much larger in higher q, fgas mergers. We explore in Section 5.5 the implications of enhanced SFRs for our equal-mass, gas-rich mergers.
Although we have found that tAGN may be enhanced by GW recoil for low-luminosity AGNs, we emphasize that the bright AGN phase (L > 10 per cent LEdd) is always shorter for recoiling BHs than for recoiling stationary BHs or at best roughly equal. This is a direct result of the peak accretion episode being disrupted at the time of the recoil kick. The important consequence of this is that the total accreted BH mass is never enhanced by GW recoil in our simulations. For moderate recoil speeds, the BH mass deficit relative to a stationary BH is generally slight, but we find that the effect of recoil is always to decrease the final BH mass. The implications of this finding are discussed in Section 5.4.
Finally, we note a possible observational consequence of longer AGN lifetimes. Fig. 10 clearly shows that when vk = 0, the two main accretion episodes are associated with the periods of rapid star formation. This means that the AGN luminosity may be difficult to distinguish from the starburst luminosity and the AGN may also be dust obscured for at least part of its active phase. If GW recoil allows the AGN to remain active long after the starburst is complete, such an AGN might be more readily observable, because it would no longer be competing with the starburst luminosity and may also be less dust obscured.
5.2 Ejected gas disc
We have seen that moderate-velocity recoils in gas-rich remnants may allow efficient accretion from ambient gas. However, a BH moving rapidly or in a low-gas-density region is unlikely to sweep up much gas from its surroundings and we would expect the probability of observing this type of active BH to be low. Indeed, in our simulations, the Bondi accretion rate of the BH rapidly declines following a large recoil kick (see Fig. 10). The possibility not accounted for in our simulations is that, if the BH is surrounded by an accretion disc at the time of recoil, some amount of this gas may remain bound to the recoiling BH. We therefore implement an analytic model to calculate the BH accretion rate from this ejected gas disc, based on the BH mass and accretion rate at the moment of recoil.
The BH carries with it the part of the accretion disc where vorb≳vk; that is, it will carry a disc with radius approximately Rej∼GMBH/v2k. This is a reasonable approximation if the BH is kicked directly out of the orbital plane of the disc, as is the case when vk≳ 500 km s−1 and the BH orbital plane is aligned with that of the disc. If the accretion rate from this ejected disc is high enough, then such a system could potentially be seen as an offset AGN via either resolved spatial offsets or offset spectral lines. BL08 predicted accretion time-scales of ≲1–10 Myr for the type of recoiling BHs we consider here (MBH∼ 106–108 M⊙; vk∼ 400–1000 km s−1), assuming a constant- model. Here, we improve on these estimates by calculating a time-dependent accretion rate for the ejected disc.
We assume the innermost region of the gas disc is well modelled by a viscous, Keplerian, thin disc (specifically, an α-disc, Shakura & Sunyaev 1973). The evolution of the disc surface density, Σ, can be described by a diffusion equation:Lynden-Bell & Pringle 1974; Pringle 1981). The steady-state solution is recovered when the left-hand side of the equation vanishes, but this is not the case for our ejected disc, which is no longer fed at its outer radius. Exact analytic solutions to equation (1) can be found when ν scales only with radius (Lynden-Bell & Pringle 1974), but in the α-disc model, ν is a function of both R and Σ. Pringle (1991) has derived self-similar solutions to this equation for viscosity scaling as ν∝ΣmRn. These solutions assume a disc with no inner radius, so for our case, we must neglect the finite inner disc radius at the innermost stable circular orbit (ISCO) of the BH (RISCO), as well as the possibility of a circumbinary gap in the disc that persists after the BH merger. This should not affect our results, however, as the ejected disc will be much larger than the ISCO (Rej/RISCO > 103) and the gap is expected to refill on a time-scale of years (Milosavljević & Phinney 2005; Tanaka & Menou 2010). For an α-disc in which viscosity is assumed to scale with gas pressure and the opacity equals the electron-scattering value (κes≈ 0.4 cm2 g−1), the viscosity is given by Pringle (1991) with m = 2/3, n = 1 (see also Cannizzo, Lee & Goodman 1990): equation (3) yields a time-dependent expression for disc mass: equation (4) along with the conditions that Md(t′) =Mej and , where t′=t0+tsim−tmrg, tsim is the time in the simulation, Mej is the ejected disc mass and is the Bondi accretion rate calculated in the simulation at the time of merger. The mass evolution of the disc is therefore given by
We can calculate the initial radius and mass of the ejected disc using the multicomponent, self-gravitating disc model of BL08. The calculation takes as input the kick velocity, BH accretion rate and BH mass at the time of the merger. Note that because we use the accretion rate calculated in the code as input, the normalization of the ejected-disc accretion rate is subject to the same assumptions as is the Bondi–Hoyle accretion rate (see Section 2). We refer the reader to BL08 for the details of the self-gravitating disc model, but we note that for the BH masses and accretion rates considered here, Rej lies in the region of the disc that is self-gravitating but still strongly dominated by the BH potential (i.e. Zone II of the BL08 model). Quantitatively, Rej for our fiducial-mass simulations is between 103.9–105.7 Schwarzschild radii or ∼0.02–0.6 pc. The corresponding mass of the ejected disc, Mej, is 103.5–107.5 M⊙ (104.6–106.3 M⊙ for our fiducial-mass models).
We use these quantities as initial conditions for the ejected disc model outlined above. The time-scale for the accretion rate (and mass) decay, t0, ranges between 105.6–107.6 yr. t0 depends on Mej and , which in turn depend on vk and the galaxy model parameters. Specifically, t0 decreases with increasing vk, fgas and q. After a time t0, the accretion rate has dropped by a factor of 2−19/16≈ 0.44 (note that the time in equation 9 is defined such that t′=t0 at tmrg). By the end of our simulations, about 55–80 per cent of the disc mass has been consumed. This corresponds to an increase in the BH mass of ΔM = 0.3–6 per cent (0.3–3 per cent for vk/vesc = 0.9 simulations only). As expected, these values of ΔM are somewhat lower than the result of BL08, ΔM≈ 10 per cent MBH. BL08 considered kick speeds as low as 100 km s−1; their result also assumed that the entire disc mass was accreted and that the accretion was Eddington-limited at the time of BH merger.
Fig. 10 demonstrates that except in high-vk simulations where the average orbital time-scale is long (vk≳ 0.8vesc), the ejected-disc accretion is typically insignificant compared to the Bondi accretion of the oscillating BH. Of course, the BH may still carry an accretion disc during this phase, but the accretion will be more accurately characterized by the highly variable Bondi rate calculated from the local density and sound speed than by an isolated-disc model. Therefore, our analytic ejected-disc model becomes important only for high-velocity recoils in which the Bondi accretion rate drops sharply after the kick. Below, we discuss the AGN lifetimes resulting from the combined accretion models, as well as the feasibility of observing offset AGNs.
5.3 Offset AGNs
In Fig. 12, we show two examples of the ΔR, Δv phase space occupied by recoiling BHs. The left-hand and right-hand panels show simulations of the q0.5fg0.3a model with vk/vesc = 0.6 and 0.9, respectively. ΔR is the binned spatial offset of the BH from the stellar centre of mass and Δv is the binned kinematic offset from the stellar centre-of-mass velocity. The contour shading shows, for each ΔR, Δv bin, the amount of time spent at that location (Δt, top panels), the total AGN energy output (ΔEbol, middle panels) and the average bolometric luminosity 〈Lbol〉 (bottom panels). Lbol is determined at each time-step from the relation , where the radiative efficiency εrad is 0.1, unless , in which case the system is considered radiatively inefficient, and (Narayan & McClintock 2008).
The vk/vesc = 0.6 recoiling BH, which settles back to the galactic centre in ∼400 Myr, fills a corner of the phase space. The Δt contours peak at low velocities, as these represent apocentric passages where the BH spends the most time. 〈Lbol〉 is above 108 L⊙ throughout the simulation, with peak values >1011 L⊙. Unlike the Δt contours, 〈Lbol〉 peaks at high velocities, because the BH has the highest accretion rates immediately after the kick and during pericentric passages. The total bolometric energy output per ΔR, Δv bin, which is the product of these two distributions, accordingly peaks at both high ΔR, Δv and low ΔR, Δv. Instead of emitting all of its AGN feedback energy from a single central point in the galaxy, the BH distributes this energy throughout the central kpc, as discussed in Section 5.1. Consequently, this BH has a longer AGN lifetime than the vk = 0 BH in the same merger model (see Fig. 11).
In contrast to this example, the vk/vesc = 0.9 recoiling BH occupies a narrow track in the phase space, as it completes only a few orbits and does not experience much drag. The Δt and 〈Lbol〉 contours in this example illustrate an inherent difficulty in observing rapidly recoiling BHs: the brightest luminosities occur just after the kick, when the BH velocities are highest, but this is also where the BH spends the least amount of time. Still, the BH has luminosities ≳108 L⊙ for ΔR up to 10 kpc. Further analysis is required to determine the observability of the recoiling BHs in both these examples, as well as those in our other simulations.
Fig. 13 shows the AGN lifetimes, tAGN, in six merger models for which we have varied vk. AGN lifetimes in the top two panels include only the Bondi accretion rate, while tAGN in the bottom two panels includes the accretion from our ejected-disc model as well. Unlike the lifetimes in Fig. 11, here, we also apply a cut to the AGN lifetime such that only activity that occurs while the BH is recoiling is shown: tAGN < tsettle. We calculate tAGN using each of the three AGN definitions described previously. The solid lines in the figure correspond to AGNs defined as Lbol > (3 per cent, 10 per cent) LEdd and the dashed lines correspond to Lbol > Lmin. The shaded regions of the panels show the range of values between these definitions. As before, we choose Lmin = 3.3 × 109 L⊙ = 1.3 × 1043 erg s−1.
A notable feature in Fig. 13 is that the peak AGN lifetimes for recoiling BHs are quite long, > 100 Myr, even for kick speeds near 1000 km s−1 in some cases. Note also that in the lower fgas mergers, tAGN peaks at the lowest vk/vesc, but in the gas-rich mergers, the peak value is vk/vesc = 0.7–0.8. This reflects the fact that we have defined tAGN < tsettle in these panels. In the gas-rich mergers, gas drag quickly damps out the low-vk trajectories and tAGN is limited by tsettle. For higher vk/vesc, gas drag is less efficient and tAGN is instead limited by the gas supply available to the BH.
In all cases, the AGN lifetimes that include only Bondi accretion (Fig. 13, top panels) fall sharply to <1 Myr for large vk/vesc. This demonstrates that Bondi accretion is inefficient when the BH is kicked far from the dense central region. However, when we also allow for BH accretion from its ejected disc (Fig. 13, bottom panels), the ejected-disc accretion enhances the active lifetime for high-vk/vesc recoils by more than an order of magnitude. In the examples shown, tAGN for low-luminosity AGNs is 3–14 Myr even for kicks as high as 1.2vesc. For even higher kick speeds, tAGN will slowly decay as the trend in the plot indicates, owing to the decreasing ejected disc mass. The relatively slow decay of tAGN with increasing vk occurs partly because, as opposed to a constant- accretion model, here, tAGN depends on the fraction of disc accreted before the BH luminosity falls below the ‘AGN’ limit. In addition, the disc mass in the self-gravitating regime has a fairly weak (r1/2) dependence on radius, so Mej∝ 1/vk. Recall, however, that the probability of recoil events with kick velocities ≫1000 km s−1 is quite small (Campanelli et al. 2007a; Schnittman & Buonanno 2007; Baker et al. 2008; Lousto et al. 2010b; van Meter et al. 2010); this more than the BH’s fuel supply limits the chances of observing a recoil event with vk≫vesc.
These findings illustrate the importance of BH accretion from the ejected disc for the lifetimes of rapidly recoiling AGNs. Also crucial, however, is the time-dependent nature of this accretion. If we were to assume a simple constant- model for the ejected-disc accretion, we would greatly overpredict the lifetime of the bright AGN phase and in some cases we would underpredict the lifetime of the low-luminosity phase. We can calculate the discrepancies between tAGN for a constant- model and the lifetimes shown in Fig. 13, restricting the comparison to vk/vesc≥ 0.9 to capture only the simulations in which ejected-disc accretion is the dominant mode. –50 per cent for all the models shown in Fig. 13 except q0.5fg0.1a. If we were to assume the accretion continues at this rate for a time (tmrg) in these models, we would overpredict tAGN for bright (L > 10 per cent LEdd) AGNs by a factor of ∼4 –40. Conversely, tAGN for L > Lmin would be underpredicted by a factor of up to ∼7. In the q0.5fg0.1a model, (tmrg) . In this case, tAGN for L > 3 per cent LEdd would be overpredicted by a factor of ∼10–40, depending on kick speed, and tAGN for L > Lmin would also be overpredicted by a factor of ∼2–4. Therefore, it is clear that inclusion of a more realistic, time-dependent model for ejected-disc accretion, such as the one presented here, is important for calculating the lifetimes and luminosities of rapidly recoiling AGNs. Similarly, for lower velocity recoils, the BH + ejected-disc system clearly cannot be treated as isolated, because accretion from ambient gas is the dominant mode of AGN fuelling.
To determine whether these recoiling AGNs could be distinguished from their stationary counterparts, we need more information than simply the active lifetimes. Recoiling AGNs could be identified as such if during their active phase they are either kinematically or spatially offset from their host galaxies; we will consider both possibilities.
The gas carried along with the recoiling BH is expected to include at least part of the BH’s broad-line (BL) region, because BL clouds reside deep within the potential well of the BH. The narrow-line (NL) region lies much farther from the BH, so, in general, this gas will not be bound to the recoiling BH and will remain at the redshift of the host galaxy. A kinematically offset AGN could therefore be seen in a spectrum as a BL feature offset from a NL feature, if the recoiling AGN is able to illuminate NL clouds as it leaves the centre of the galaxy. Otherwise, such a system might have a spectrum with no NLs and BLs that are offset from the redshift of the host galaxy’s stellar light. Because the physics of the NL region are not included in our analysis, we consider the velocity offset between the recoiling BH and the stellar centre of mass. We calculate the AGN lifetime as before, based on the BH luminosity at each time-step, but, in this case, we are only interested in the lifetime of offset AGNs; that is, we calculate the time for which the BH is active and exceeds a minimum velocity offset Δvmin. We note that the limit Δvmin is imposed on the physical velocity difference between the BH and the host galaxy. The line-of-sight (LOS) velocity seen by a typical observer will be lower, so the velocity-limited AGN lifetimes shown here are upper limits.
Fig. 14 shows the results of our analysis for Δvmin = 500 km s−1 (top panels) and 800 km s−1 (middle panels). For Δvmin = 500 km s−1, the offset AGN lifetime for low-fgas, low-vk mergers is much smaller than the total AGN lifetime in Fig. 13. vesc is relatively low in these models, so kick speeds are also lower for fixed vk/vesc. In the q0.5fg0.1a merger, for example, vesc = 777 km s−1, so the BH never exceeds the 500 km s−1 limit for vk/vesc≤ 0.6. The offset AGN lifetimes for low vk/vesc in gas-rich models are a factor of 5–10 lower than the total AGN lifetimes, but the peak lifetimes are still ∼10–50 Myr. At high vk/vesc(≥0.8), there is little difference between tAGN with or without the velocity cut. In these cases, the ejected disc dominates the BH accretion and the BH speed is above 500 km s−1 for nearly all of the AGN phase.
When Δvmin is increased to 800 km s−1, the AGN lifetimes drop sharply for most merger models and kick speeds (see the middle panels in Fig. 14). Only in the equal-mass, gas-rich models shown here (q1fg0.6a and q1fg0.3a) do the BHs appear as bright offset AGNs and for only a narrow range of kick speeds. These two models have vesc > 1100, such that BHs kicked with ∼800 km s−1 may have tightly bound orbits. In these cases, the BHs make many rapid pericentric passages through the dense central region, where they encounter an ample gas supply from which to accrete. Interestingly, in the q1fg0.6a and q1fg0.3a models, kicks <800 km s−1 (≤0.6vesc) have non-zero lifetimes in these panels. These are examples of the phenomenon discussed in Section 4.2, in which the BHs are kicked before vesc reaches its maximum, so the BHs receive mild velocity boosts on subsequent passages through the increasingly deep central potential well.
In the other (lower-q, fgas) models, the central region is less dense, such that bound BHs have lower kick speeds and lower accretion rates during pericentric passages. In these models, only ejected-disc accretion is efficient at producing an offset AGN with Δv > 800 km s−1, which occurs only when the BHs are ejected entirely from the host galaxy. These ejected AGNs have lifetimes <10 Myr and low luminosities (tAGN < 1 Myr for L > 10 per cent LEdd).
To be observed as a spatially offset recoiling AGN, a recoiling BH must travel far enough from the galactic centre to be resolved as distinct point source before exhausting its fuel supply. As with the velocity offsets, we define the spatial offsets relative to the stellar centre of mass. The maximum galactocentric distance achieved by a recoiling AGN in any of our fiducial-mass simulations is ∼15 kpc by our lowest-luminosity AGN definition and only ∼1 kpc by the strictest definition. In our high-mass simulations, which have larger ejected discs, offsets ≳50–100 kpc may be achieved while L > Lmin, but the maximum offset while L > 3–10 per cent LEdd is still only a few kpc. Offsets of ∼1 kpc correspond to an angular separation of ∼0.55 arcsec at z = 0.1; these could be resolved with the HST/JWST as long as they are not observed edge-on. Offsets ≳3–4 kpc could be resolved at this redshift with the SDSS and Chandra. To determine the lifetimes of the spatially offset AGN in our simulations, we impose a minimum-distance cut of ΔRmin = 1 kpc from the stellar centre of mass and identify an offset AGN as one that simultaneously has ΔR > ΔRmin and meets our AGN-luminosity criteria. The bottom panels of Fig. 14 show the results of this analysis. We find that above vk/vesc = 0.5–0.7 (physically, the kick speed above which Rmax exceeds 1 kpc) the offset AGN lifetimes can be long, but only for low-luminosity AGNs. In these models, tAGN ranges from 2 to 110 Myr for L > Lmin and from 0.3 to 8 Myr for L > 3 per cent LEdd. In none of our simulations does the offset exceed Rmin for more than 1 Myr while its luminosity is >10 per cent LEdd. The reason that spatially offset AGNs are only seen at low luminosities is as follows. For recoiling BHs with short-period oscillations, Bondi accretion is dominant and most of the bright AGN activity occurs during pericentric passages with small ΔR. For BHs on long-period orbits or escaping BHs, the ejected disc accretion dominates. In this case, the brightest AGN phase occurs as the BH leaves the central dense region and once the BH passes ΔRmin, its luminosity has already begun a monotonic decay.
We can conclude from Fig. 14 that in most of our models, low-luminosity recoiling AGNs may be distinguishable via spatial offsets >1 kpc for kick speeds ≳0.5–0.7vesc, but as kinematically offset AGNs with Δvmin = 800 km s−1, they may be distinguishable only for kick speeds ≳vesc. These offset AGNs generally have luminosities ≲3 per cent LEdd and their luminosity owes mainly to accretion from the disc ejected with the BH. The exceptions are equal-mass, gas-rich mergers (here, q1fg0.6a and q1fg0.3a), in which recoiling BHs on bound trajectories may experience multiple phases as bright (L > 10 per cent LEdd), kinematically-offset AGNs accreting from ambient gas during pericentric passages. This scenario can only occur in merger remnants with large central escape speeds (≳1100–1300 km s−1), such that BHs may receive kicks ≳800 km s−1 and still undergo short-period oscillations. In lower q, lower fgas mergers, the central gas density is lower and the kick speeds for bound trajectories are smaller, so this bright, velocity-offset phase does not occur. If kinematic offsets can be observed with a resolution of 500 km s−1 instead of 800 km s−1, then velocity-offset-AGN lifetimes increase significantly and measurable velocity offsets can be attained for a much wider range of kick speeds.
In the lowest q, fgas merger for which we have varied vk (the q0.5fg0.1a model), the offset AGN lifetimes are always <10 Myr for L > Llim and <2 Myr for L > 3 per cent LEdd. Our q0.5fg0.1a model has 10 per cent gas initially, but by the time of the recoil kick only ∼4 per cent of the baryonic mass is in gas. We can conclude that merger remnants with gas fractions lower than this and with q≲ 0.5 are unlikely to fuel bright offset AGNs according to our criteria, at least for recoiling BHs with MBH∼ 107 M⊙. Larger BHs will have larger gas reservoirs and longer AGN lifetimes, and the reverse will be true of smaller BHs. For example, in the q0.5fg0.1a model, a BH with vk/vesc = 0.9 has a velocity-offset (Δvmin = 800 km s−1) AGN lifetime <1 Myr by all AGN criteria, while in the 10 times more massive q0.5fg0.1M10x model, the velocity-offset tAGN is 27 Myr (for L > Lmin). In general, larger BHs will have longer offset-AGN lifetimes for fixed vk/vesc. However, larger BHs are found in larger galaxies, where a higher kick speed is required to reach an appreciable fraction of the escape speed, and such kicks have a lower probability of occurring. We discuss recoil kick probabilities further in Section 6.1.
As a final caveat, we note that when Bondi accretion is dominant, the offset-AGN lifetimes are subject to the assumption that no significant delay occurs between when material becomes bound to the BH and when it is actually accreted. The brightest AGN phases of this type occur during pericentric passages, when the velocity is high and the BH is subsequently travelling away from the galactic centre. Therefore, substantial delays in gas accretion would likely decrease lifetimes of velocity-offset AGNs and increase lifetimes of spatially-offset AGNs.
5.4 MBH–σ* relation
Several authors have suggested that GW recoil may contribute to scatter in the observed black hole–host galaxy bulge relations due to ejected BHs (Volonteri 2007) or bound recoiling BHs that leave the central galactic region (BL08; Sijacki et al. 2010). Volonteri (2007) used semi-analytic models to estimate the effect of ejected BHs at high redshifts on the BH mass–bulge stellar velocity dispersion relation at z = 0. In this work, we have already demonstrated that, even when they remain bound to the host galaxy, recoiling BHs can undergo substantially different accretion histories from their stationary counterparts and that the resulting AGN light curves may also vary significantly as a function of kick speed (see Figs 10, 13 and 14). Here, we attempt to quantify the effect of individual GW recoil events on the MBH–σ* relation. Because we have simulated all of our merger models both with no recoil kick and with vk = 0.9vesc, we are able to compare the MBH–σ* relations that result from each of these samples. The final BH masses and LOS-averaged stellar velocity dispersions from both sets of simulations are plotted in Fig. 15; the right-hand panel also shows simulations with vk/vesc = 0.4–0.8. We include only models with the fiducial merger orientation angles (orbits a and w–z) in this panel to avoid heavily weighting the fit towards the (q, fgas) combinations for which we have varied the orbit.
By examining the mass deficits of these recoiling BHs relative to the stationary BHs in our simulations, we can quantify the interruption to BH growth caused by GW recoil. Note that we are probing only the effect of GW recoil on the growth of kicked BHs; the actual BH mass deficit resulting from BHs merging in a cosmological framework may differ from these values. In particular, the vk = 0.9vesc simulations result in BH wandering times of at least several Gyr, so in reality these BHs will be lost to the galaxy and never contribute to the MBH–σ* relation. If the galaxy experiences subsequent mergers, especially minor mergers, the kicked BH may be replaced with a smaller BH from the incoming galaxy, creating a larger mass deficit. We also stress that the data plotted here are not meant to represent a random sample, as we have varied parameters in a systematic way. Accordingly, the relations shown are not expected to reproduce the observed MBH–σ* relation. In particular, the scatter in the vk = 0 relation (0.13 dex) is smaller than the observed scatter (∼0.25–0.3 dex, e.g. Tremaine et al. 2002; Aller & Richstone 2007; Hopkins et al. 2007b), owing to a number of factors, such as the single merger orbital configuration used in this plot, the constant feedback efficiency assumed for all simulations and small range of total galaxy mass spanned by the majority of our models. Furthermore, although scaling each kick to 0.9vesc essentially allows us to capture the maximum effect of a single recoil event on BH growth, this does not reflect the fact that in galaxies with lower escape speeds, large vk/vesc recoil events will occur more frequently than in massive galaxies. Such a trend with galaxy mass would steepen the slope of the MBH–σ* relation, though this would also depend on the distribution of other merger parameters, such as q and fgas, with varying galaxy mass. Despite these limitations in comparing our simulations to the observed MBH–σ* relation, our results provide substantial insight into how GW recoil affects the growth of kicked BHs and the inherent complexities of this process. The relative differences in Fig. 15 between the simulations with and without recoil are robust and can easily be understood in physical terms.
It can be seen by visual inspection of Fig. 15 that the effect of high-velocity recoils is an overall downward offset in normalization and an increase in scatter. Because the LOS-averaged values of σ* do not change significantly in the presence of a recoil kick, we focus our discussion on variation in the final BH mass.
The downward shift in Fig. 15 reflects the fact that GW recoil always reduces the final mass of a merged BH. We quantify the ‘mass deficit’ of recoiling BHs relative to stationary BHs by calculating the fractional mass difference for each recoil simulation relative to its vk = 0 counterpart, [Mfin(0) −Mfin(vk)]/Mfin(0), where Mfin is the final BH mass in each case. Fig. 16 shows the mass deficit plotted versus vk for the six merger models discussed previously. The deficits generally increase with vk, though the details vary substantially between models. In the q0.5fg0.1a model, recoiling BHs with vk=vesc have approximately one-third lower final mass than do stationary BHs, while in the q1fg0.3a model, a vk=vesc recoil results in a BH almost five times smaller than its stationary counterpart. Note also that the curves in Fig. 16 all level out at vk/vesc = 0.7–0.9; this marks the kick speed in each model above which recoiling BHs are unable to accrete more gas after their ejected disc has been exhausted. The models with lower fgas achieve this limiting value at lower vk/vesc, owing to the lower gas drag and shallower potentials in these remnants. The blue diamonds in Fig. 15 mark the positions of these varying-vk simulations on the MBH–σ* relation. The points corresponding to lower vk/vesc lie close to the no-recoil (black) line; those with higher vk have lower BH masses.
In Fig. 15, the normalization offset between the vk = 0 BH population and the vk/vesc = 0.9 population is 0.44 dex, and the offset between the vk = 0 BH population and the total population is 0.22 dex. However, it has been demonstrated that the normalization of the black hole–galaxy bulge correlations depends on the efficiency of AGN feedback (Di Matteo et al. 2005; Robertson et al. 2006c; Hopkins et al. 2007a). In our simulations, we assume a constant fraction (5 per cent) of the BH’s bolometric luminosity couples to the surrounding gas as thermal energy. Stationary accreting BHs with lower feedback efficiencies would grow larger before heating the surrounding gas enough to slow down or halt accretion. Depending on how much accretion occurs prior to versus after the BH merger, higher feedback efficiencies could reduce the normalization offset between our recoiling and stationary BH populations. In general, however, stationary BHs will accrete more after the BH merger than recoiling BHs, so there is likely to be a non-zero downward shift in MBH–σ* normalization caused by GW recoil.
As previously stated, the disproportionate effect of recoil events on galaxies with small vesc could have some effect on the MBH–σ* slope. In our simulations, where we have scaled kick speeds to vesc, there is no significant difference between the MBH–σ* slope for recoiling versus stationary BHs. This is to be expected, because the depth of the central potential well scales more or less self-similarly to the total galaxy mass. We do, however, see a trend towards larger BH mass deficits for higher q mergers, which arises because gas is driven more efficiently to the galactic centres during coalescence, and recoiling BHs ‘miss’ a larger accretion phase when they are kicked out of this central dense region.
The variation in BH mass deficits between merger models results in larger intrinsic MBH–σ* scatter for BH populations that have undergone recoils. The correlation for the no-recoil sample has a scatter of 0.13 dex and the vk/vesc = 0.9 correlation has a scatter of 0.24 dex. In other words, the high-recoil sample has almost twice as much intrinsic scatter as the no-recoil sample. This increase in scatter is a universal consequence of GW recoil; recoil events essentially add an extra variable to the determination of final BH masses. As opposed to stationary BHs, the final mass of recoiling BHs depends not only on how much gas is available to be accreted, but also on when this accretion occurs relative to the time of the kick. In mergers with no recoil, MBH is a strong function of the amount of gas driven into the central region (i.e. the depth of the central potential well). The details and timing of the cold gas flow in galaxy mergers depend non-trivially on factors such as q, fgas, the galaxy orbits and the SFR, so it is natural that larger intrinsic scatter is found in the final masses of recoiling BHs. If the kick speed is lower, such that the BH can settle back to the centre while an ample gas supply is still present, the BH mass deficit will be reduced, but the details of this late-phase accretion add yet another element of unpredictability to the final value of the BH mass.
In addition to the larger intrinsic MBH–σ* scatter for a BH population with fixed vk/vesc, we also anticipate an increase in overall scatter caused by the range of kick speeds that would occur in a realistic population of BHs and by BHs that have undergone multiple recoils through their formation history. The contribution of GW recoil to the total MBH–σ* scatter is difficult to predict from our simulations; quantitative predictions would require following detailed BH formation histories in a cosmological framework, which is beyond the scope of this paper. However, we at least gain a sense of how a range of kick speeds affects the MBH–σ* relation by plotting our varied-vk simulations along with the vk = 0 and vk/vesc = 0.9 populations. The right-hand panel in Fig. 15 shows the same data as the left-hand panel but also includes data for the simulations in which we have varied vk/vesc (blue diamonds). The blue line is a fit to these points and the green line is a fit to all data shown. The fit to all the data has a scatter of 0.25 dex, which as expected is larger than the scatter for the vk = 0 population.
The most robust of our conclusions from this analysis is the increase in intrinsic MBH–σ* scatter for BHs that have undergone a single recoil event, by a factor of ∼2 in our simulations. This finding suggests that GW recoil events may be a non-negligible contribution to the scatter in the observed MBH–σ* relation, especially because BHs at z = 0 may have undergone multiple recoil events throughout their formation history.
5.5 Evolution of the central galactic region
Fig. 17 compares the SFRs for simulations with no GW recoil (black dashed curve), with vk/vesc = 0.9 (red solid curve), and with no BHs (blue dotted curve) in a gas-rich, equal-mass merger model (q1fg0.4a). It is clear that after the merger, the simulation with a large recoil kick has a higher SFR than the simulation with vk = 0. The post-merger SFRs are similar in the high-recoil case and in the simulation with no BH at all. This same behaviour can be seen in the high-recoil simulations of model q0.5fg0.3a shown in Fig. 10, though to a lesser degree. These enhanced SFRs occur because when the BH is removed from the centre of the merger remnant, the AGN feedback is also displaced. The lack of central energy input allows star formation to continue unhindered in this region for a longer time. We can actually see the effect of the recoiling AGN feedback on the instantaneous SFR in Fig. 17. The SFR curve for the vk/vesc = 0.9 simulation has a vaguely periodic shape; in fact, each small trough in the SFR coincides with a pericentric passage of the BH through the galactic centre, which imparts a small amount of feedback energy to the dense region. This supports the connection between the recoiling BH dynamics and the galactic SFR, but the main consequence of the recoil event is that star formation is never as strongly quenched as when the BH does not recoil. In the example shown in Fig. 17, the merger remnant with a stationary BH has almost twice the amount of gas at the end of the simulation as the merger with vk/vesc = 0.9. This corresponds to an increase in total star formation of 3.4 × 109 M⊙ in the recoil simulation or ∼3 per cent increase in the total stellar mass, solely due to the displacement of the BH and its feedback.
In Fig. 18, we can clearly see the effect of this enhanced star formation on the remnant stellar density profile. The figure shows the remnant stellar density profiles for the same three simulations as shown in Fig. 17, with vk = 0 (black curves), vk/vesc = 0.9 (red curves) and without BHs (blue curves). The results are further broken down into the contribution from stars formed prior to the galaxy merger (‘old’ stars, thin curves) and those formed during merger-driven starbursts (‘new’ stars, thick curves). In all cases, the ‘new’ stellar population dominates in the central kpc, but in the remnant with a recoiling BH the central stellar density is almost twice that of the remnant with a stationary BH. In fact, this stellar profile very closely resembles that of a remnant without a BH. Therefore, recoil events in gas-rich mergers can prolong starburst phases, creating denser and more massive stellar cusps. We have also examined the (unattenuated) U−B colour evolution of the stellar component in these no-recoil, large-recoil and no-BH simulations. The colours in all three cases are similar; by the end of the simulation, all have U−B∼ 0, that is, these merger remnants are entering the ‘green valley’ in transition from blue to red. However, in the vk/vesc = 0.9 and no-BH simulations, which have prolonged star formation, the stars are slightly bluer than in the no-recoil simulation. While this difference is small, ∼0.03 mag, it further demonstrates that in merger remnants with rapidly recoiling BHs, the (blue) starburst phase may be prolonged and hence the transition to red elliptical slightly delayed, solely due to the effect of GW recoil. In Section 6.5, we discuss the implications of these findings for observable properties of ULIRGs and elliptical galaxies.
We note that this is the opposite effect of one that has been proposed for gas-poor mergers, where repeated BH oscillations through the central region of the merger remnant may scour out a stellar core (Boylan-Kolchin et al. 2004; Gualandris & Merritt 2008). Such core-scouring is expected to occur on scales of less than or equal to tens of pc, which are below our resolution limit. However, we do comment that when even a small amount of gas is present in our merger simulations, it is efficiently driven to the centre of the merger remnant.
Using the SPH/N-body code gadget-3, we have generated a suite of galaxy merger simulations both with and without a recoil kick applied to the central SMBH at the time of BH merger. Owing to the large range of parameters sampled, we were able to analyse systematic trends in recoiling BH behaviour with variation in galaxy mass ratio, total galaxy mass, initial gas fraction, orbital configuration and BH merger time. In addition, we have used a range of recoil kick velocities, which are scaled to the galactic central escape speed in each case. We followed the trajectories and accretion history of the recoiling BHs, as well as evolution of host galaxy properties, such as the SFR and depth of the central potential well. For the BH accretion, the Bondi–Hoyle accretion rates were used, but we also included a time-dependent model for accretion from a gas disc carried along with the recoiling BH.
6.1 Recoil trajectories and kick probabilities
The recoiling BH trajectories in our simulations are characterized by low-angular-momentum orbits (i.e. ‘oscillations’ about the galactic centre) that are also highly three-dimensional. In other words, the baryonic component of the merger remnant dominates the BH orbits even if they extend well into the halo, but the unrelaxed nature of this remnant creates complicated BH trajectories. We find that the amplitude and duration of these oscillations vary widely between different galaxy models and kick velocities; some BHs settle back to equilibrium at the galactic centre within a few Myr, while others may remain on large orbits through the halo for a Hubble time.
The amplitude and duration of BH oscillations in a given recoil event clearly depend on the central escape speed, that is, the depth of the central potential well. However, by normalizing kick speeds in our simulations to vesc at the time of the kick, we are able to compare different galaxy models regardless of their escape speed. We find substantial variation in recoil trajectories between models even for fixed vk/vesc. Recoil oscillations are more readily suppressed in mergers with higher q and fgas. Such merger remnants are more centrally concentrated and gas rich, resulting in steeper central potential wells and increased gas drag and dynamical friction.
In the models for which we have varied vk/vesc, we find that, with the exception of our collisionless simulations, recoiling BHs with vk/vesc = 0.4–0.7 reach galactocentric distances ≲1 kpc and settle back to the centre in <0.5 Gyr. The corresponding range of kick speeds is 310–880 km s−1.
For large recoil velocities, GW recoil events can produce BH oscillations that persist through the end of our simulations, 2.9 Gyr after the kick. The largest orbits can extend well into the galactic halo, with galactocentric distances up to ∼100 kpc. In our models, these long-lived oscillations occur for kicks ≳0.6–0.9vesc, again depending on the mass ratio and gas content of the merger.
Furthermore, as illustrated in Section 4.2, in gas-rich mergers with q∼ 1, the central escape speed may increase rapidly around the time of BH binary coalescence. In this case, the trajectory of the recoiling BH may also depend on the time of the BH merger relative to the formation time of the central density cusp. Of all the factors influencing the BH trajectories, this is the most unpredictable, as it depends on the relative timing of two short-time-scale processes. This adds an element of uncertainty to any attempt, including this study, to model populations of recoiling BHs in gas-rich galaxies. Our assumption of rapid BH mergers (tmrg=tcoal) means that our BH oscillation amplitudes may be upper limits in mergers with strong vesc evolution. However, the effect of this uncertainty on the results presented here is limited, due to the relatively small number of merger models affected and our use of the normalization vk/vesc.
We have also tested whether recoil trajectories are affected by the direction of the kick. Our generic merger remnants (those that are not formed from highly aligned galaxy orbits) are highly disturbed at the time of the BH recoil and therefore do not have disc-like structures, that is, there is no ‘special’ direction in which a recoiling BH would clearly experience much greater drag. Even with the clumpy, irregular density profiles of our merger remnants, we find that the trajectories of recoiling BHs have little dependence on the kick direction. Much more important in determining the amplitude and duration of recoiling BH oscillations is the magnitude of the kick velocity relative to the escape speed, along with the aforementioned global galaxy properties, such as the initial gas fraction and mass ratio. We note that there should indeed be a dependence of the recoil trajectory on the kick direction with respect to the plane of the accretion disc; recoils directed into the disc plane likely experience more gas drag and have shorter oscillation time-scales. In addition, Rossi et al. (2010) have shown that the ejected disc mass and AGN luminosity may depend strongly on the direction of the recoil kick with respect to the accretion disc plane. We do not resolve this small-scale region in our simulations; however, recoil kicks ≳500 km s−1 will always be nearly perpendicular to the orbital plane of the binary (Campanelli et al. 2007a,b; Lousto & Zlochower 2009), which is unlikely to be perpendicular to the orbital plane of the accretion disc. Therefore, for the recoil velocities we consider (≳300 km s−1), we can assume that the BH is unlikely to be kicked directly into the accretion disc plane.
We mention in Section 2 that according to the results of NR simulations, the recoil kick speeds in our fiducial-mass merger models, which we have scaled to the escape speeds of ∼700–1250 km s−1, should occur with fairly high probabilities. Even if the BH spin magnitudes are randomly distributed between 0 and 1, 41.6 per cent of major mergers should have kicks above 500 km s−1 and 13.5 per cent should have kicks above 1000 km s−1. Furthermore, Lousto et al. (2010b) have conducted statistical analysis of GW recoil spins and kick speeds using post-Newtonian approximations calibrated to NR results. They provide kick probability distributions for several mass ratio bins, which we can use to estimate the relative probability of ‘high’ (vk∼vesc) or ‘moderate’ (vk∼ 0.6vesc) kick speeds in various merger models. Coincidentally, for mass ratios 0.3 ≤q≤ 1.0 and 300 ≲vk≲ 1000 km s−1, which represent a large fraction of our fiducial-mass simulations, the kick probability distribution is quite flat. For these models, the relative probability of a given kick speed between models, as well as the relative probability of high versus moderate vk in a given model, is roughly unity. In our highest mass (q1fg0.3M20x) simulation, the probability of a kick near vesc = 2500 km s−1 is about 10 times smaller than the probability of a kick with 0.6vesc and is also ∼10 times smaller than the kick probability for vk∼vesc in the q1fg0.3a model. For mergers with q≲ 0.25, which we do not model in our study, the kick distribution becomes heavily skewed towards velocities <500 km s−1.
An important caveat to the interpretation of recoil kick probabilities is that the relevant velocity for observations is the LOS velocity, not the actual velocity. For random lines of sight, the probability of observing a given kick vk as a kinematic offset >Δvmin is therefore reduced by a factor P(ΔvLOS > Δvmin) = 1 −Δvmin/vk. For example, about 20 per cent of recoils with vk = 1000 km s−1 will have ΔvLOS > 800 km s−1. Lousto et al. (2010b) find that the overall probability distribution for ΔvLOS is lower than the vk distribution by a factor of between ∼2 (for kicks >500 km s−1) and ∼20 (for kicks >2500 km s−1). Similarly, the probability of observing a recoiling BH with spatial offset ΔRLOS > ΔRmin is ; that is, a majority of spatial offsets ≳1.2Rmin will have LOS offsets exceeding the minimum observable offset.
6.2 Recoiling AGN light curves
The final coalescence of two merging, gas-rich galaxies typically triggers both a burst of star formation and a phase of rapid BH accretion. Assuming the BHs themselves coalesce on a short time-scale, the resulting recoil event coincides with a period of high activity in the central region of the merger remnant. Thus, we expect recoiling BHs in gas-rich mergers to have substantially different accretion histories from their stationary counterparts. This is indeed the behaviour seen in our simulations, with greater disruption to the BH accretion history for higher kick velocities.
Recoil events that produce short-period BH oscillations (vk/vesc∼ 0.4–0.8, depending on the merger model) cause the BH accretion rate to drop at the time of the kick, thereby reducing the peak luminosity of the BH. However, the BH remains in the central few kpc of the merger remnant where gas densities are large enough that Bondi–Hoyle accretion is efficient. Because the BH’s feedback energy is deposited in a much larger volume than for stationary BHs, it is unable to heat the gas sufficiently to completely terminate the AGN phase. The recoiling BH may therefore experience a longer active phase by a factor of a few than its stationary counterpart, although at lower luminosities. This phase continues in some cases until the end of the simulation, 2.9 Gyr after the BH merger. Owing to the lower luminosities of these extended AGN phases, we find that they do not increase the final BH mass relative to stationary AGNs. However, this type of accretion episode is distinct from other merger-triggered AGN phases in that it is not accompanied by a simultaneous starburst. The copious amounts of dust produced by starbursts can easily obscure an active central source and the starburst luminosity may overwhelm the AGN luminosity. Recoiling AGNs that continue accreting long after the starburst is complete therefore may be more easily observable.
When recoil kicks are large enough to produce long-period (≳100 Myr) BH oscillations (vk/vesc≳ 0.6–0.9 depending on the merger model), the Bondi–Hoyle accretion rate drops precipitously as the BH is rapidly ejected to a lower density region. Even though the BH remains on a low-angular-momentum orbit and returns several times to the dense central region, its velocity is generally too high during these pericentric passages to allow for substantial gas accretion and the gas is depleted by star formation between subsequent passages.
In such high-velocity recoil events, a portion of the inner accretion disc may remain bound to the BH. We account for this by assigning to each recoiling BH a time-dependent, isolated accretion disc with a mass and initial accretion rate based on the BH mass, recoil velocity and the accretion rate at the time of the kick. The mass of these ejected discs is only a few per cent of the BH mass, so they do not contribute significantly to the final BH mass. However, the accretion from these discs does greatly extend the active lifetime of rapidly recoiling BHs in the regime where Bondi–Hoyle accretion is inefficient. We find that the AGN lifetime in our models may increase from ≪1 Myr to ∼20 Myr when ejected-disc accretion is included, though these are generally low-luminosity AGNs. The peak lifetimes in the Bondi–Hoyle regime are generally longer and correspond to brighter AGNs, but for a fairly narrow range of kick speeds. tAGN in the Bondi–Hoyle regime may be up to ∼300 Myr for low-luminosity AGNs and between ∼10–100 Myr for AGNs accreting at 3–10 per cent LEdd.
6.3 Offset AGN lifetimes
From an observational standpoint, our primary interest is not the total active lifetime, but rather the time for which the BH is active and can be distinguished as a recoiling BH. We find that kinematically offset AGNs – those that are luminous enough to be detected and have velocity offsets large enough to be spectrally resolved – occur in two distinct physical scenarios. In the first scenario, the recoiling AGNs travel far from the galactic centre and are powered by accretion from the ejected disc. This occurs for recoil velocities near or above the central escape speed. In our simulations, AGNs with velocity offsets >800 km s−1 have lifetimes ≲10 Myr in most cases. In most of our galaxy models, recoiling AGNs can only meet this velocity criterion for vk > vesc. For reference, the SDSS quasar study conducted by Bonning et al. (2007), which yielded a null result, had a spectral offset limit of 800 km s−1.
The other physical mechanism for producing kinematically offset recoiling AGNs is via multiple pericentric passages through a central dense region. Note that marginally bound recoiling BHs may also experience close passages, but these will be rapid and few in number; the cumulative probability of observing such an event is low. If a recoiling BH is able to settle back to the galactic centre within ∼0.5–1 Gyr, it will undergo many short-period oscillations through the central region. In this case, the BH may be able to accrete efficiently from the ambient gas during high-velocity close passages, producing a kinematically offset AGN. This scenario favours massive, gas-rich galaxies with high central densities: the larger the central escape speed and the steeper the potential well, the higher the velocity with which a BH can be kicked and still remain tightly bound to the galaxy. In our q = 1, gas-rich mergers, offset AGNs with Δv > 800 km s−1 can have lifetimes up to 100 Myr and may have lifetimes ≳10 Myr at bright (>10 per cent LEdd) luminosities. Especially for these AGNs, an increase in spectral resolution to achieve a limit of 500 km s−1 could make a substantial difference; in our models, the lifetimes increase by a factor of ∼5–10 with this lower value of Δvmin.
The maximum spatial offset achieved by a recoiling AGN in our high-mass simulations is ΔR = 112 kpc by our most generous AGN definition (L > 1.3 × 1043 erg s−1) and 3.2 kpc by our strictest definition (>10 per cent LEdd), but for our fiducial-mass simulations, these maximum offsets are only 15 and 1.3 kpc, respectively. To calculate spatially offset AGN lifetimes, we require a minimum offset ΔRmin = 1 kpc. Such an offset could be resolved by the HST or JWST at z≲ 0.1 or by the SDSS or Chandra at z≲ 0.03. In our simulations, recoiling BHs achieve Rmax > 1 kpc for vk/vesc≳ 0.5–0.7. AGN lifetimes may be quite long near this threshold kick speed (∼100 Myr), though only at low luminosities. At higher kick speeds, lifetimes are generally ∼1–10 Myr and slowly decrease with higher vk, as less gas is bound to the recoiling BHs.
Based on our calculated AGN lifetimes, we conclude that for escaping BHs (vk > vesc), probabilities for observing recoils via spatial offsets and kinematic offsets are similar. For recoiling BHs on bound trajectories, velocity-offset AGNs with high luminosities are favoured for massive, gas-rich galaxies (q = 1, fgas = 0.3–0.6 in our models). Spatially offset AGNs with lower luminosities have similar lifetimes for all models. Recoiling BHs in merger remnants that are both gas poor (fgas≲ 0.1 initially) and result from lower q mergers (q≲ 1/3) are not expected to produce long-lived offset AGNs.
6.4 MBH–σ* scatter and offset
We have demonstrated that recoiling BHs have lower final masses than their stationary counterparts, though the mass deficit varies considerably with kick speed and merger remnant properties. This may affect the formation of the observed correlations between SMBHs and the bulges of their host galaxies. In particular, we consider the possible effects of these mass deficits on the observed MBH–σ* relation by comparing the MBH–σ* correlations of our sets of vk = 0 and vk = 0.9vesc simulations, as well as our simulations with varying vk. We do not find any evidence for a change in MBH–σ* slope caused by recoils with fixed vk/vesc, but in a realistic population of galaxy mergers, galaxies with smaller vesc would be disproportionately affected by recoils with vk near or above vesc. This could contribute to some steepening of the MBH–σ* slope. In our simulations, we find that GW recoil events may contribute to a downward shift in normalization, an increase in overall scatter and an increase in intrinsic scatter of the MBH–σ* relation (see Fig. 15). The magnitude of these effects, particularly the normalization offset and total scatter, is sensitive to the galaxy population considered, and the detailed merger and recoil histories of the central BHs. Therefore, we cannot make quantitative predictions about the contribution of GW recoil to the observed normalization and scatter of the MBH–σ* relation; this would in fact be difficult to predict with any method. A more realistic approach involving BH and galaxy populations derived from a cosmological framework would still be plagued with uncertainties in the recoiling BH masses. We have shown that these masses are sensitive to numerous factors that do not lend themselves to semi-analytical modelling, including the (possibly evolving) central potential depth, the BH merger time-scale and the detailed BH accretion history.
Despite these uncertainties, the existence of a recoiling-BH contribution to the MBH–σ* correlation and its scatter is a necessary consequence of GW recoil events. In particular, we demonstrate that a universal effect of GW recoil is to increase the intrinsic scatter in BH masses. The disruption to BH accretion caused by a recoil event introduces an extra variable into the determination of the final BH mass, namely, the timing of the recoil kick relative to the merger-induced burst of BH accretion (i.e. quasar phase). In our recoil simulations, the intrinsic MBH–σ* scatter of our rapidly recoiling BH population (vk/vesc = 0.9) is almost two times the scatter in our stationary-BH population. This is the contribution to scatter caused solely by interruption of BH growth by a recoil kick; other factors not accounted for here could easily increase the contribution of GW recoil to MBH–σ* scatter. While a more realistic population of BHs will have smaller average mass deficits resulting from a range of kick speeds, the BHs are also likely to have undergone numerous mergers and recoil events in the past, which would increase the overall scatter. Subsequent merger events, in which an incoming BH may replace an ejected BH, could also increase scatter. The relatively small scatter of the empirical relation, 0.25–0.3 dex, places constraints on the contribution of GW recoil and hence on the frequency of large vk/vesc recoil events. Similarly, a non-zero contribution of GW recoil to the observed scatter places constraints on the contribution from other sources, such as redshift evolution of the MBH–σ* relation (e.g. Haehnelt & Kauffmann 2000; Shields et al. 2003; Robertson et al. 2006c).
6.5 Broader implications: starbursts, ULIRGs and elliptical galaxies
Previously, a model for galaxy evolution has been outlined in which mergers play a central role (e.g. Hopkins et al. 2006a, 2008b; Somerville et al. 2008). In this picture, galaxies form originally as discs and the population of ellipticals is built up through time from mergers. Gas-rich mergers pass through intermediate phases in which they host significant star formation and black hole activity, first as ULIRGs, then as quasars and finally as dormant elliptical galaxies.
There is a great deal of observational support for this scenario. Nearby ULIRGs are invariably associated with gas-rich mergers (Sanders et al. 1988a; Sanders & Mirabel 1996) and contain large quantities of molecular gas in their centres (Sanders, Scoville & Soifer 1991), driving nuclear starbursts. Many display evidence for central AGNs through their ‘warm’ spectra (Sanders et al. 1988b), leading to the interpretation that quasars are the descendants of ULIRGs. As the starbursts and quasar activity fade, the merger remnant evolves into a passive elliptical (Toomre & Toomre 1972; Toomre 1977), making the transition from blue to red. The relics of this process can be seen in the central light concentrations in elliptical galaxies that were put into place during the starburst phase of their creation (e.g. Mihos & Hernquist 1994a; Hopkins et al. 2008c), and correlations between SMBHs and, for example, the central potentials of their hosts (Aller & Richstone 2007; Hopkins et al. 2007b).
These fossil signatures, and the timing of the events that led to their formation, are sensitive to the interplay between star formation and black hole growth, and the feedback processes that accompany each of these. In principle, black hole recoil could impact all of these phenomena to varying degrees, providing constraints on the modelling we have presented here as well on the interpretation of the observations more generally.
For example, some ULIRGs display evidence for significant AGN activity in the central starburst region. This is strong evidence that the BHs in these systems did not undergo rapid recoil events and that any recoil motion was quickly damped out. However, merger-triggered central AGNs, particularly those in ULIRGs, may be dust obscured or dwarfed by the luminous starburst. AGNs that are instead kicked out of the central region by a large recoil may be more easily detected, provided a large enough region of the galaxy is observed. Additionally, we have shown that in some cases, bound recoiling AGNs may have longer lifetimes than stationary AGNs. Such AGNs may continue shining after the starburst is complete, which may also make them easier to observe.
Moreover, in the simulations, at least, feedback from black hole growth plays a significant role in quenching star formation (Di Matteo, Springel & Hernquist 2005; Springel, Di Matteo & Hernquist 2005a), allowing the remnant to evolve from a blue, actively star-forming galaxy to a red, quiescent object. During this transition, the remnant can, for some period of time, be classified as a K+A or E+A galaxy (Snyder 2010). If the central BH were ejected, star formation would linger at higher rates than for a stationary BH, extending the duration of the blue to red transition, modifying the stellar populations and ultimately yielding an elliptical with denser, more massive central stellar cusps. Were the BH not to return eventually to the nucleus, then subsequent gas-poor mergers would be less strongly affected by BH scouring, obscuring the relationship between core and cusp ellipticals of the same mass (Hopkins & Hernquist 2010b).
Fig. 18 shows that if the BHs recoil at sufficiently high speed, the remnant relaxes to a state similar to what would have happened had the galaxies not had black holes. Compared to the case of no BH recoil, the remnants with recoiling BHs (or no BHs) have more massive central stellar cusps. These more massive stellar cusps reflect a more extended period of active star formation following final coalescence of the progenitor galaxies, lengthening the period of time to make the blue to red transition. Furthermore, the reduced efficiency of AGN feedback means that the central regions would remain more heavily dust obscured than if there were no BH recoil, making the remnant less likely to be visible as a K+A galaxy (Snyder 2010) and potentially significantly reducing the number of K+A galaxies expected to be observable as merger remnants.
Using our suite of merger models, we can make some qualitative predictions about the morphologies of recoiling BH hosts. Fig. 3 shows that our merger remnants have some common morphological features. In general, at the time of BH merger (and recoil), the merging galaxies have undergone final coalescence (i.e. they no longer have two distinct cores), but are still highly disturbed and have prominent tidal features. This is true even when the progenitor galaxies have low initial gas fractions (4 per cent) or relatively small ratios (0.25). For minor mergers with q < 0.25, the disruption to the primary galaxy would be less pronounced, but such mergers would also yield small GW recoil kicks. Highly aligned (i.e. coplanar) galaxy merger orbits may result in remnants with regular disc structures (see Fig. 3), but these orbits should comprise only a small fraction of major galaxy mergers. Therefore, we conclude that the majority of galaxy merger remnants that produce large GW recoil events have single cores but are still tidally disturbed at the time of the recoil kick. Furthermore, we have shown that the offset-AGN lifetimes for recoiling BHs in our simulations are <100 Myr and are in most cases ∼1–10 Myr. Tidal structures in these major merger remnants generally persist for at least a few hundred Myr, so it follows that most observable recoiling AGNs should be found in unrelaxed remnants that may still have visible rings or tidal tails.
6.6 Comparison to recent work
Shortly before this paper was completed, Guedes et al. (2010) and Sijacki et al. (2010) completed independent studies of recoiling BHs in gaseous mergers. However, both these papers are exploratory in nature, in that only a few examples of galaxy mergers are studied. While there is agreement between these studies on some fundamental conclusions, our large parameter study allows us to expand upon these findings, connect them into a more coherent picture and present some entirely new results. Here we give a brief comparison of this work to the work of Guedes et al. (2010) and Sijacki et al. (2010).
Guedes et al. (2010) use one equal-mass and two minor merger remnants as initial conditions for their simulations. The mass and spatial resolution used are comparable to ours, though in their equal-mass merger, higher resolution is obtained at late stages via particle splitting. The merger orbits were coplanar and prograde, a configuration that results in highly disc-like remnants (see our Fig. 3). Recoils with a range of kick speeds were applied to the BH and after a brief initial phase, the BH trajectories were followed semi-analytically. We note that they use vk/vesc as low as 0.14, while we use only vk > 0.4vesc, such that the BHs do not spend most of their time at radii ∼Rsoft.
Sijacki et al. (2010) use isolated galaxy models as well as one galaxy merger model for their recoiling BH simulations. Like our simulations, their simulations were performed using gadget-3, with similar spatial resolution. However, they use only one merger model, which is massive and gas-rich, resulting in a remnant with a very high central escape speed (3500 km s−1) and steep central potential. This model is therefore fairly extreme, in that recoils with vk∼vesc are required to eject the BH beyond the central few kpc.
Despite the substantially different approaches of these two studies and our own, a common finding is that gas inflow during galaxy mergers can create centrally-dense remnants that impart significant gas drag to recoiling BHs, thereby reducing the amplitude and return time of their trajectories. That this effect remains important even for modest gas fractions (≲0.1) illustrates the limitations of modelling GW recoil dynamics in purely collisionless systems.
Guedes et al. (2010) find that in their minor mergers, recoiling BHs exceed 600 km s−1 only for vk > vesc and they do not calculate AGN lifetimes for escaping BHs. Kinematically offset AGNs have lifetimes up to a few Myr in their q = 1 merger. Sijacki et al. (2010) note that velocity offsets up to 500 km s−1 may occur during the bright AGN phase in their isolated disc model, but do not discuss offset AGNs in their merger model.
Although they use similar values of Lmin and ΔRmin, Guedes et al. (2010) calculate longer spatially offset AGN lifetimes than we find in our models. The disparity in offset AGN lifetimes may arise partly from the different dynamical treatments; Guedes et al. (2010) calculate the recoil trajectories semi-analytically, as opposed to our N-body approach. Further, they assume for the ejected-disc accretion rather than a time-dependent, isolated-disc model. They argue that spatially offset AGNs are naturally longer lived than kinematically offset AGNs, because the former occur near the apocentre and the latter near the pericentre. However, we show that in some cases, kinematically offset AGNs may have longer lifetimes, because the BH is much more likely to be actively accreting during passages though the central dense region.
As we have found in our recoil simulations, Sijacki et al. (2010) show that recoil events interrupt a phase of rapid BH accretion in their gas-rich merger, causing a BH mass deficit relative to a stationary BH. They note that BH mass deficits may introduce scatter into BH mass–host galaxy relationships, which we demonstrate robustly in Section 5.4.
Furthermore, both our results and those of Sijacki et al. (2010) show that the central SFR is higher after a large recoil, due to the displacement of AGN feedback. Sijacki et al. (2010) note that this further depletes the gas supply available to the BH, once it returns to the centre. We demonstrate that SFR enhancement due to recoil is strongest in gas-rich, nearly-equal-mass mergers, though it is present in all of our simulations to some degree, and we have outlined several possible observational consequences in Section 6.5.
In addition to connecting and expanding upon the results of Guedes et al. (2010) and Sijacki et al. (2010), many of our results are entirely new to this work. As we have discussed these previously, we list them here only briefly. First, because we allow the BH merger time to be a free parameter, we are able to demonstrate that recoil trajectories may depend quite sensitively on the BH merger time in q∼ 1, gas-rich mergers. We eliminate much of this uncertainty in our simulations by scaling vk to vesc(tmrg). This approach also enables a much clearer comparison between models. Also unique to our work is that we test a range of recoil kick directions in generic merger remnants, demonstrating that, in general, the kick direction is less important for recoil trajectories than the kick magnitude.
Another novel feature of our models is the inclusion of a time-dependent accretion disc around the recoiling BH. We also calculate recoiling AGN lifetimes using three different luminosity limits, which allows us to differentiate between bright and faint AGNs. We learn from this analysis that some recoiling BHs may actually have longer AGN lifetimes than stationary BHs, though at low luminosities. We also show that kinematically offset AGNs occur via two different mechanisms: repeated pericentric passages in massive, gas-rich remnants and ejected-disc accretion in high-velocity recoils (vk > vesc in our models). We demonstrate that spatially offset AGNs can occur for a wide range of kick speeds and merger models.
Finally, we use our suite of simulations to construct MBH–σ* relations for stationary and rapidly recoiling BH populations. We find that GW recoil introduces a normalization offset, and larger overall and intrinsic scatter into this correlation.
We have shown that recoiling BHs may be observable as offset AGNs; discoveries of these objects would both inform models of hierarchical galaxy formation and constrain event rates for future GW observatories. In the meantime, the richly varied effects of GW recoil on SMBHs and their host galaxies are not to be discounted as important components of merger-driven galaxy evolution.
This work was supported in part by NSF grant AST0907890, and NASA grants NNX08AL43G and NNA09DB30A (for AL). TJC thanks the Keck Foundation for their support.