THE ASSEMBLY OF SUPERMASSIVE BLACK HOLES AT HIGH REDSHIFTS

We construct DM merger trees based on the algorithm by Volonteri et al. (2003), which allows only binary mergers. Similar numerical algorithms (e.g. Somerville & Kolatt 1999; Cole et al. 2000) give somewhat different results, as have been discussed recently by Zhang et al. (2008). We will not reproduce the Volonteri et al. (2003) recipe here in full; instead we present a brief review. We take the extended Press–Schechter (EPS) mass function (Press & Schechter 1974; Lacey & Cole 1992),

$$\begin{eqnarray} \frac{dN}{dM}&=&\frac{1}{\sqrt{2\pi }}\frac{M_{0}}{M}\frac{d\sigma ^{2}_{M}}{dM} \frac{\delta _{c}(z)-\delta _{c}(z_{0})}{\left(\sigma ^{2}_{M}-\sigma ^{2}_{M_{0}}\right)^{3/2}} \nonumber\\ &&\times \exp \left\lbrace -\frac{1}{2}\frac{\left[\delta _{c}(z)-\delta _{c}(z_{0})\right]^{2}}{\sigma ^{2}_{M}-\sigma ^{2}_{M_{0}}}\right\rbrace, \end{eqnarray} \tag{ 2 }$$

which gives for a parent halo of mass M₀ at redshift z₀ the number of progenitor halos in the range M ± dM/2 at a higher redshift z. Here, σ_M is the amplitude of the linear matter fluctuations at redshift z = 0, smoothed by a top hat window function whose scale is such that the enclosed mass at the mean density is M (computed using the fitting formula for the transfer function provided in Eisenstein & Hu 1999), and δ_c is the redshift-dependent critical overdensity for collapse. We take the limit as Δz ≡ z − z₀ → 0:

$$\begin{equation} \frac{dN}{dM}=\frac{1}{\sqrt{2\pi }}\frac{M_{0}}{M}\frac{d\sigma ^{2}_{M}}{dM} \;\left(\sigma ^{2}_{M}-\sigma ^{2}_{M_{0}}\right)^{-3/2} \;\frac{d\delta _{c}}{dz} \;\Delta z. \end{equation} \tag{ 3 }$$

The two advantages of taking this limit are that (1) by linearizing the expression, the z and M dependences separate, allowing a tabulation as functions of the parent and progenitor halo mass, and (2) separating Δz allows for a simple algorithm for adaptive time steps to make sure that fragmentations produce binaries at most (no triplets).

For a parent halo of mass M_p after a small step Δz, the mean number of "minor" fragments in the mass range M_lo < m < M_p/2 is given by

$$\begin{equation} N_{p}=\int _{M_{\rm lo}}^{M_{p}/2} \frac{dN}{dM}dM \propto \Delta z. \end{equation} \tag{ 4 }$$

We choose Δz adaptively such that N_p ≪ 1, which ensures that multiple fragmentations are unlikely in a given single time step. We place a lower limit of 10⁻³ (in redshift units) for the time step to keep computation times manageable. The integral in Equation (4) diverges as M_lo → 0, making it computationally prohibitive to compute the merger history of arbitrarily small halos. To avoid this problem, all progenitors below a fixed mass resolution M_res are considered jointly as accreted mass and not tracked individually. Any progenitor with M < M_res(z) is thus discarded from the tree and its prior history is disregarded. For our calculations, we choose M_res(z) to be the mass corresponding to a virial temperature of 1200 K, corresponding to M_res ∼ 4.7 × 10⁵ M_☉ at z = 40 and 3.4 × 10⁶ M_☉ at z = 10. Theoretical studies have concluded that Pop-III stars can start forming at lower virial temperatures,² but numerical considerations have forced us to adopt a somewhat larger threshold. We do not impose an explicit upper redshift limit, and we run the tree until our halos at z = 6 are entirely broken up into progenitors with M < M_res. As a check on our Monte Carlo merger tree algorithm, in Figure 1, we present the progenitor mass functions of a 10¹² M_☉, z₀ = 6 parent halo at redshifts of z = 8, 13, 21, and 34, together with the Poisson errors of the merger tree output and the predictions from the EPS conditional mass function (Equation (2)). Our merger tree results are consistent with the EPS conditional mass function up to redshift z ≈ 40, with agreement within a factor of 2 for most mass bins and redshift values. Our merger tree results are consistent with the EPS conditional mass function up to redshift z ≈ 40, with agreement within a factor of 2 for most mass bins and redshift values. In particular, the numerical mass function agrees well with the EPS prediction for the low-mass progenitors, even at very high redshift, but the higher-mass progenitors are underpredicted by a factor of up to 2. We note that Cole et al. (2000) used a numerical algorithm similar to the one we adopted, to construct a halo merger-tree. As discussed in Zhang et al. (2008), that algorithm results in a similar inaccuracy.

Zoom In Zoom Out Reset image size

The conditions under which the first BHs form are highly uncertain, though numerical simulations (Abel et al. 2002; O'Shea & Norman 2007) do provide useful indications. We parameterize our ignorance in terms of a seeding fraction, such that a fraction f_seed of all halos reaching the critical virial temperature T_seed form a seed BH. There are physical mechanisms that make a low seeding fraction plausible: the first stars may form only in rare, baryon-rich overdense regions with unusually low angular momentum, and seed remnants may receive ejecting kicks from collapse asymmetry mechanisms similar to those responsible for high-velocity pulsars. Furthermore, radiative and other feedback processes may prohibit H₂ formation, cooling, and star formation in the majority of low-mass minihalos at high redshift (e.g., Haiman et al. 1997; Mesinger et al. 2006). Since the LISA event rate, especially at the earliest epochs, will depend primarily on the abundance of BHs present, it is highly sensitive to the seeding function.

We choose two fiducial seeding models, the first with T_seed = 1200 K (the minimum value required for metal-free molecular line cooling and star formation) for a Pop-III remnant seed BH with m_seed = 100 M_☉. The second model has T_seed = 1.5 × 10⁴ K and m_seed = 10⁵ M_☉, inspired by the "direct-collapse" models of more massive BHs from the central gas in halos with a deep enough potential to allow atomic cooling (Oh & Haiman 2002; Bromm & Loeb 2003; Volonteri & Rees 2005; Begelman et al. 2006; Spaans & Silk 2006; Lodato & Natarajan 2006). If Eddington accretion is the main mode of growth, then we do not expect the choice of seed mass for each type of model to qualitatively affect our results, other than the obvious linear scaling of the overall BH mass function with the initial seed mass. Only the binary mass ratios affect recoil magnitudes, and the subsequent orbital dynamics depend minimally on the BH mass.

The DM profile for the earliest halos is found to be similar to the NFW (Navarro et al. 1997) profile of lower redshift, more massive DM halos (Abel et al. 2000; Bromm et al. 2002; Yoshida et al. 2003). However, the composition and spatial distribution of the baryons, at the time when the seed BH appears in these halos, is poorly understood, and is unconstrained by observations. This is unfortunate, since these quantities play a pivotal role in determining the orbital dynamics and growth rate of a recoiling BH.

A steep profile with a cusp will retain BHs more effectively, owing both to a deeper gravitational potential well and a larger dynamical drag force at the halo center. The baryon distribution will also determine the accretion history of the central BH by determining the accretion rate as the BH rests near the halo's potential center, or as it oscillates in a damped orbit through the halo following a recoil displacement event.

In addition, whether the baryons are gaseous or stellar has nontrivial consequences, owing to the difference in the dynamical friction force between the two cases. A collisional medium provides a greater dynamical friction force than a stellar or DM medium with the same density profile (Ostriker 1999; Escala et al. 2004). Because of the difference in the drag force, an environment dominated by gas, and not by stars (or DM), has several possible consequences on BHs: (1) binaries coalesce more rapidly; (2) a BH recoiling in gas has a higher likelihood of being retained in its parent halo; and (3) any "vagrant" BH that is displaced from the baryon-rich center of the gravitational potential of its host halo takes less time to return there, reducing episodes of suppressed accretion. In three-dimensional simulations of star formation in metal-free minihalos suggest that star formation is inefficient, with either a single star, or at most a few stars, forming at the center of the halo (Abel et al. 2000; Bromm et al. 2002; Yoshida et al. 2003, 2008). Since in the context of this paper we are concerned with the pre-re-ionization universe, we work with the assumption that stars are rare before z ≳ 6 and that the baryons in our halos are mostly gaseous.

We model each galaxy as a spherically symmetric mass distribution with two components: a DM halo with an NFW profile, and a superimposed baryonic component. Previous studies on this subject (see, e.g., Volonteri et al. 2003; Madau & Quataert 2004) have often assumed a noncollisional singular isothermal sphere (SIS) profile for the mass distribution. This is justified if the gas does not cool significantly below the virial temperature of the DM halo, and if it has little angular momentum (so that it is supported thermally, rather than by rotation). In most halos whose virial temperature is above 10⁴ K, this assumption is less justified, and a disk may form at the core of the DM halo (Oh & Haiman 2002). The direct-collapse scenario in Begelman et al. (2006) and also VR06 adopt such a "fat disk" configuration. However, the central densities of such disks are within the range of our adopted spherical profiles. For simplicity, here we only consider three different prescriptions for spherical gas distribution. Our fiducial gas profile is a cuspy, ρ ∝ r^−2.2 power law, where we have taken the power-law index of 2.2 as suggested by numerical simulations of the first star-forming minihalos (Machacek et al. 2001). This profile is established in halos that are able to cool their gas via H₂, and describes the gas distribution at the time of the first star formation.

It is possible, however, that the typical seed BHs are surrounded by a very different gas distribution, at the time of their formation. First, the progenitor Pop-III stars of the first seed BHs are here assumed to form in DM halos of mass ∼10^5–6 M_☉. The UV radiation from the star will photoheat, and easily blow out most of the gas from low-mass minihalos, even before the star collapses to leave behind a seed BH (e.g., Whalen et al. 2004). In this case, the remnant BH will find itself in a DM halo devoid of gas, and can only start accreting once a merger with another, gas-rich halo has taken place, or until the parent halo has accreted enough mass to replenish its gas (e.g., Alvarez et al. 2008). We therefore make the simple assumption that no gas is present, until the minihalo containing the newly formed seed BH merges with another halo, or grows sufficiently—assumed here to be a factor of 10—in mass. However, we will examine the consequences of this assumption below, by performing runs without such a blow-out phase.

Second, as mentioned above, feedback processes may prohibit H₂ formation and cooling in the majority of the low-mass minihalos (e.g., Haiman et al. 1997; Mesinger et al. 2006). The gas in such minihalos remains nearly adiabatic, and cannot contract to high densities. To allow for this possibility, we will study a variant for the effective gas profile. Specifically, we adopt the gas distribution in these halos by the truncated isothermal sphere (TIS) profile proposed in Shapiro et al. (1999), which has an r⁻² profile at large radii, but has a flat core at the center owing to the central gas pressure. The density profile is normalized (here, and also in our fiducial gas profile above) such that the total baryon-to-DM mass ratio inside the virial radius r₂₀₀ equals the cosmological value. Both the DM and the baryonic components are assumed to extend out to 10r_vir, at which point the density falls to the background value. This is consistent with infall-collapse models of Barkana (2004).

We next have to make important assumptions about the treatment of mergers between DM halos and their resident BHs.

First, we consider the merger between two DM halos, with the more massive halo referred to as the "host" and the less massive as the "satellite" halo. The Press–Schechter formalism and our merger tree consider as "merged" two halos that become closely gravitationally bound to each other. However, if the mass ratio of a halo pair is large, then in reality the smaller halo can end up as a satellite halo, and its central BH will never merge with that of the more massive halo. We therefore require in our models that for BHs in such halo pairs to coalesce, the halo merger timescale must be shorter than the Hubble time. We take the standard parameterization of the merger time:

$$\tau _{\rm merge}\approx x\frac{M_{1}}{M_{2}}\tau _{\rm dyn}\approx 0.1 x\frac{M_{1}}{M_{2}}t_{\rm Hub},$$

where M₁/M₂ > 1 is the ratio of the halo masses, x ≲ 1 is some dimensionless factor that encodes the orbit geometry (e.g., circularity) and dynamical friction, and τ_dyn and t_Hub are the dynamical and Hubble times, respectively. It has been suggested (Boylan-Kolchin et al. 2006) that radial infall along filaments may be preferred in the mergers of elliptical galaxies. Boylan-Kolchin et al. (2008) give a fitting formula for the merger time based on numerical simulations. Their Equation (5) reduces approximately to τ_merger/τ_dyn ≈ 0.45 for a moderately noncircular orbit with circularity (defined as the ratio of the orbit's specific angular momentum to the angular momentum of the circular orbit with the same energy) of 0.5. We take a moderate value of x ≈ 0.5. That is, if M₁/M₂ < 20, then the BH in the smaller halo is considered to be "stuck" out in the orbiting satellite halo and never merges with the central BH of the more massive halo. This choice also ensures that the vast majority of BH binaries in our simulation do not have extreme mass ratios, as the BH masses co-evolve with the host halos.

We next make assumptions regarding the timescales involving BH dynamics in their host halos, as follows: (1) if the two merging halos each contain a central BH, the two holes are assumed to form a binary efficiently, i.e., we assume there is no delay, in addition to the time taken by the DM halos to complete their merger; (2) the binary is then assumed to coalesce in a timely fashion, prior to the interaction with a third hole; and (3) the binary coalescence is assumed to take place at the center of the potential of the newly merged halo. The first assumption has been addressed by Mayer et al. (2007), who report that the increased drag force of gas in wet mergers allows the timely formation of SMBH binaries. The second assumption is valid for binaries in extremely gas-rich environments (see Milosavljevic & Merritt 2003 for a review), or in triaxial galaxies (Berczik et al. 2006). As for the third assumption: given that the timescales of orbital damping are comparable to the intramerger timescale for BH velocities of ≳σ_SIS, unperturbed BHs free-falling during the halo merger process are likely to end up near the center of the potential within the merging timescales of their hosts. We do not include triple-BH interactions in our analysis.

The assumption that BH binaries merge efficiently following the mergers of their host halos may be unjustified in our models in which initially, the DM halo is devoid of gas, since gas is generally believed to be necessary for prompt coalescence. However, this inconsistency will not affect our conclusions, for the following reasons. First, we find that a BH merger in a gasless environment is a rare event, as it occurs only if both parent halos are low-mass halos that had formed seed BHs relatively recently (i.e., neither halos have yet grown in mass by a factor of 10). Second, members of such binaries will have equal (or, in actuality, similar) masses, since they have not been able to add to their seed mass by accretion. If the BHs can merge efficiently without gas, the coalesced product will likely be ejected, owing to a shallow halo potential and a relatively large kick of an equal-mass merger. If the binary does not coalesce efficiently, it will coalesce once the parent halo merges with a gas-rich, BH-free halo, or once the parent halo accretes enough gas to facilitate the merger. Such belated mergers will also presumably take place with a mass ratio of close to unity and will likely result in ejection, regardless of whether significant gas accretion takes place before coalescence. Now, consider the case of a stalled binary encountering a third BH before gas enrichment of the halo. If the third BH is much more massive, it will not be ejected by gravitational interaction or recoil. There will be a massive BH in the center of the host halo, and whether the two smaller seed BHs were ejected or swallowed by the larger BH is of little consequence to our analysis, especially given the rarity of double-seed binaries. Thus, inefficient binary coalescence is of concern only when a double-seed binary encounters a third BH of comparable mass before the host halo is gas-enriched. Such triple-seed systems are likely to be extremely rare indeed, and unlikely to affect the overall mass function at z = 6. We anticipate that the main effects of a gas-depleted host halo will be increasing the number of similar-mass mergers following the initial epoch of seed formation, and slightly reducing the time available for gas accretion.

2.5.1. Probability Distribution of Kick Velocities

To calculate the recoil velocities of coalesced BHs, we employ the formulae provided in Baker et al. (2008), which are used to fit their numerical results,

$$\begin{eqnarray} v_{\rm recoil}^{2}&=v_{m}^{2}+v_{\perp }^{2}+v_{\parallel }^{2} +2v_{\perp }\left(v_{m}\cos \xi +v_{\parallel }\sin \xi \right), \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} v_{m}&=A\eta ^{2}\sqrt{1\hbox{--}4\eta }\left(1+B\eta \right), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} v_{a^{\parallel }}&=H\frac{\eta ^{2}}{1+q}\left(a_{1}\cos \theta _{1}-qa_{2}\cos \theta _{2}\right), \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} v_{a^{\perp }}&=&-K\frac{\eta ^{3}}{1+q} [a_{1}\sin \theta _{1}\cos \left(\phi _{1}-\Phi _{1}\right) \nonumber\\ &&-qa_{2}\sin \theta _{2}\cos \left(\phi _{2}-\Phi _{2}\right)], \end{eqnarray} \tag{ 8 }$$

where q ≡ m₂/m₁ ⩽ 1 is the binary mass ratio, η ≡ q/(1 + q)² is the symmetric mass ratio, and θ_1,2 are the angles between the BH spin vectors $\vec{a}_{1,2}=\vec{S}_{1,2}/m_{1,2}$ and the binary orbital angular momentum vector. The angles ϕ_1,2 denote the projection of the spin vectors onto the orbital plane, measured with respect to a fixed reference angle. As seen from the equations themselves, v_m is the kick component that depends only on the symmetric mass ratio; $v_{a^{\parallel }}$ and $v_{a^{\perp }}$ depend on the mass ratio and the projection of the binary spins parallel and perpendicular, respectively, to the orbital angular momentum.³ Φ₁(q) = Φ₂(1/q) are constants for a given value of q that encode the dependence of the kick and orbital precession on the initial spin configuration. We use the mean values given in Baker et al. (2008) for the fitting parameters: A = 1.35 × 10⁴ km s⁻¹, B = −1.48, H = 7540 km s⁻¹, K = 2.4 × 10⁵ km s⁻¹, and ξ = 215°. We assume spherical symmetry in our host DM halos, so we are concerned only with the recoil magnitudes and not with the kick orientations.

Following Schnittman & Buonanno (2007), for every recoil event, we assign to both members of the binary spin magnitudes in the range 0.0 ⩽ a_1,2 ⩽ 0.9, randomly selected from a uniform distribution. We consider two scenarios for the spin orientation: a case where the orientation is completely random, with 0 ⩽ θ_1,2 ⩽ π and 0 ⩽ ϕ_1,2 ⩽ 2π chosen randomly from a uniform distribution; and one where the spins are completely aligned with the orbital angular momentum vector.⁴ The latter scenario is motivated by Bogdanovic et al. (2007), who argued that external torques (such as those provided by a circumbinary accretion flow) may help align the binary spins with the orbital angular momentum prior to coalescence, making kicks of ≳200 km s⁻¹ physically unfavored. While the argument was originally used to explain the lack of quasars recoiling along the line of sight at lower redshifts (Bonning et al. 2007; although see Komossa et al. 2008 for a candidate recoil detection), the same spin–orbit alignment will also impact the early assembly history of SMBHs, by allowing less massive halos to retain recoiling BHs.

Berti & Volonteri (2008) have provided a merger-tree model to follow the spin evolution of BHs through gas accretion and binary merger events. In this work, we opt not to track the spins of individual BHs due to the uncertainties involved. For instance, if circumbinary disks can act to align the spins of each binary member, this would present a scenario significantly different from the scenario presented by Berti & Volonteri. As we will show in later sections, the spin prescription does not appear to play a primary role in determining the mass function of z = 6 BHs.

We show the recoil velocity distribution for both orientation scenarios as a function of the mass ratio q in Figure 2. The figure shows the mean, 1σ, and maximum values for the recoil velocity magnitude from 10⁶ random realizations at a given value of 0.01 ⩽ q ⩽ 1. For q > 0.1, if spins are randomly oriented then kicks for similar-mass mergers are in the 100–1000 km s⁻¹ range, with a handful of kicks above 1000 km s⁻¹ and a maximum possible kick of ≈3000 km s⁻¹; for spins aligned with the orbital angular momentum, kicks are typically below 200 km s⁻¹, with the maximum allowed kick no more than 300 km s⁻¹. For q ≳ 0.1 and for random spin orientations, the maximum kick v_max(q) is achieved close to where v_∥ ≫ v_⊥ is maximized; this occurs when a_1,2 = 0.9, cos(ϕ_1,2 − Φ_1,2) = 1, and sin θ₁ ≈ −sin θ₂ ≈ 1. At these spin parameter values, v_max(q) is well approximated by $\sqrt{v_{m}^{2}+v_{\parallel }^{2}}$, and is a monotonically increasing function of q. If the spins are aligned with the orbital angular momentum vector, then v_∥ = 0 and the maximum kick occurs when a₁ = 0.0, a₂ = 0.9. Also, in this aligned case, the spin-independent component v_m is the dominant term for q ≲ 0.6 and peaks at q ≈ 0.345, resulting in the maximum and mean values for the recoil speed peaking between these q values.

Zoom In Zoom Out Reset image size

2.5.2. Trajectories of Kicked Black Holes

Given the mass distribution of the host halo and a recoil speed generated from the method detailed in the previous subsection, we numerically integrate the equation for the radial motion of a BH with mass m,

$$\begin{equation} \frac{dv}{dt} = -\frac{GM(r)}{r^{2}} + a_{\rm DF} - v\frac{\dot{m}}{m}, \end{equation} \tag{ 9 }$$

where r(t) is the radial displacement of the BH from the center of the host halo and v(t) is the BH's radial velocity. The first term on the right-hand side is the acceleration due to Newtonian gravity with M(r) the total (DM + baryon) mass enclosed inside spherical radius r; the second is the drag deceleration due to dynamical friction; and the third is the deceleration due to mass accretion. A similar calculation of the kicked BH's trajectory has been performed by Madau & Quataert (2004)—the main difference from our prescription is that they assumed the halo to have a collisionless SIS profile, and adopted parameters describing galactic bulges in the local universe, whereas we use the hybrid DM+gas profile described above, and adopt parameters relevant to low-mass halos at high redshifts.

For a noncollisional medium (in our case, for DM), the dynamical friction is described by the standard Chandrasekhar formula (see, e.g., Binney & Tremaine 1988),

$$\begin{eqnarray} a_{\rm DF}^{\rm star}(r, v)&=&-4\pi G^{2}m\rho (r)\;\frac{1}{v} \;\ln \Lambda \nonumber\\ &&\times\left[{\rm erf}(X)-\frac{2}{\sqrt{\pi }}X\exp (-X^{2})\right], \end{eqnarray} \tag{ 10 }$$

where ln Λ is the Coulomb logarithm and $X=v/(\sqrt{2}\sigma _{\rm DM})$, with σ_DM the velocity dispersion of the DM halo. In a collisional medium, the density wave in the wake of an object traveling at near or above the sound speed is enhanced via resonance, an effect that has no counterpart in collisionless media. This results in an enhancement of the dynamical friction force, for which Ostriker (1999) has derived an analytical formula. However, the Ostriker prescription is known to overpredict the drag force at slightly supersonic velocities when compared with numerical simulations. While Escala et al. (2004) provides a fitting formula that better fits the numerical results at low speeds, their formula suffers from the opposite problem, and overpredicts the drag for highly supersonic motion. We therefore adopt a hybrid prescription, adopting the Escala et al. (2004) formula for motion with $\mathcal {M}<\mathcal {M}_{\rm eq}$ and the Ostriker formula for $\mathcal {M}>\mathcal {M}_{\rm eq}$, where $\mathcal {M}_{\rm eq}$ is the Mach number (= |v|/c_s) where the two prescriptions predict equal drag. The entire prescription is described by

$$\begin{equation} a_{\rm DF}^{\rm gas}(r, v)=-4\pi G^{2}m\rho (r)\frac{1}{v} \;\times \; f(\mathcal {M}),\qquad {\rm where} \end{equation} \tag{ 11 }$$

$$\begin{equation} f(\mathcal {M})= \left\lbrace \begin{array}{@{}ll@{}} 0.5\ln \Lambda \left[{\rm erf}\left(\frac{\mathcal {M}}{\sqrt{2}} \right)-\sqrt{\frac{2}{\pi }} \mathcal {M}\exp (-\mathcal {M}^{2}/2)\right] &\\ \quad \quad \qquad {\rm if} \;\;\;0\le \mathcal {M} \le 0.8;\\ 1.5\ln \Lambda \left[{\rm erf}\left(\frac{\mathcal {M}}{\sqrt{2}} \right)-\sqrt{\frac{2}{\pi }} \mathcal {M}\exp (-\mathcal {M}^{2}/2)\right] & \\ \quad \qquad \quad {\rm if}\;0.8<\mathcal {M} \le \mathcal {M}_{\rm eq};\\ \frac{1}{2}\ln \left(\frac{\mathcal {M} +1}{\mathcal {M} -1}\right)-\ln \Lambda & \\ \quad \quad \qquad {\rm if}\;\;\;\;\;\;\;\;\mathcal {M}>\mathcal {M}_{\rm eq}.\end{array} \right. \end{equation} \tag{ 12 }$$

The Coulomb logarithm ln Λ is not a precisely known parameter, but is generally agreed to be ≳1 for both the stellar and the gaseous cases. We adopt the value ln Λ = 3.1 used in Escala et al. (2004), which yields $\mathcal {M}_{\rm eq}\approx 1.5$.

The drag force depends on the local gas sound speed. Instead of attempting to compute a temperature profile explicitly, we make the approximation that the gas sound speed is constant and given by the isothermal sound speed of the halo virial temperature. We believe this to be justified from numerical simulation results that show the gas temperature to vary by at most a factor of ≈3 within the virial radius despite a steeper-than-isothermal (∝r^−2.2) density profile (see, e.g., Machacek et al. 2001). While local variations in the sound speed may have significant effects when v ∼ c_s, the recoil events of interest here are for the most part highly supersonic. Recoil events with v ∼ c_s will result in quick damping of the BH orbit and for the purposes of this paper will in all likelihood be indistinguishable from the zero-recoil calculation in terms of their accretion history.

The virial temperature is given by (e.g., Barkana & Loeb 2001)

$$\begin{equation} T_{\rm vir}\approx 370(1+z)\left(\frac{\mu }{0.6}\right) \left(\frac{M}{10^{7}\,M_{\odot }}\right)^{2/3}\left(\frac{\Omega _{0}h^{2}}{0.14}\right)^{1/3} \;\mathrm{K}, \end{equation} \tag{ 13 }$$

and the isothermal sound speed is

$$\begin{equation} c_{s}\approx 1.8\;(1+z)^{1/2}\left(\frac{M}{10^{7}\,M_{\odot }}\right)^{1/3} \left(\frac{\Omega _{0}h^{2}}{0.14}\right)^{1/6}\mathrm{km}\;{\rm s}^{-1}. \end{equation} \tag{ 14 }$$

We also employ a simplified prescription for the velocity dispersion of noncollisional matter, given by the SIS expression evaluated at the virial radius: $\sigma _{\rm SIS}=\sqrt{GM/2r_{200}}$. The simplified expression agrees with the exact velocity dispersion for the NFW profile within ∼20% inside the virial radius.

Because matter that is bound to the BH does not contribute to dynamical friction, we follow Madau & Quataert (2004) and truncate the density profiles at the BH radius of influence, r_BH ≈ Gm/σ²_SIS. The density is furthermore assumed to be constant inside this radius. Although the BH will drag with it the surrounding gravitationally bound matter, effectively increasing its mass, the additional mass is small owing to the large initial recoil velocity (e.g., Lippai et al. 2008), and we have ignored this mass enhancement here.

Figure 3 illustrates the effect of the halo matter distribution on the BH orbit. For each of the three orbits shown in the figure, we adopt M = 10⁸ M_☉, m = 10⁵ M_☉, z = 20, and v_kick = 100 km s⁻¹. The black curve shows, for reference, the radial orbit for a pure NFW profile. The red curve corresponds to the case which includes an NFW DM component and a power-law gas profile with ρ ∝ r^−2.2. The blue curve is for a halo with an NFW DM component and a TIS gas profile.

Zoom In Zoom Out Reset image size

In our calculations, we are mainly interested in whether the kicked BH is ejected and lost (i.e., cannot contribute to the final SMBH mass at z = 6), or is retained (i.e., eventually returns to the nucleus, and can be incorporated into the z = 6 SMBH). We place the following retention condition for recoiling BHs: the BH must return to within 1/10 of the virial radius of the newly merged host halo within 1/10 of the Hubble time. The fate of a (SM)BH placed in an orbit extending to the outskirts of its host halo is uncertain: even if it is not lost through tidal interactions with an incoming merging halo, it is not likely to form a binary that hardens efficiently. We therefore impose the above conservative cutoff, in order to avoid tracking these vagrant BHs. Our retention threshold velocity, v_ret, above which recoiling BHs do not return within the prescribed time limit and are considered lost, is a numerically calculable function of M, m, z and halo composition. In order to minimize computation time, we tabulate v_ret in the range 5 < z < 40, 10⁵ M_☉ < M < 10¹⁵ M_☉ and 10⁻⁶ < m/M < 1 and approximate the result with a fitting formula that accurately reproduces the exact numerical results within 5% in the tabulated range. In principle, we should compute the retention velocity as a function of the time to the next merger experienced by the halo. However, we find the time dependence to be weak. The distribution of the time intervals between halo mergers in a given merger tree has a sharp peak at ≈10⁻¹t_Hub, with far fewer mergers occurring at ∼10⁻²t_Hub and ∼t_Hub. At these tails of the distribution, v_ret varies by ≲10% from the value for 10⁻¹t_Hub. We find that v_ret ∼ 5–8 × σ_SIS. This is comparable to the escape velocity for a nondissipative pure SIS profile that is truncated at the BH radius of influence, v_esc ≈ 5σ_SIS if m = 10⁻³M (YM04).

The weak dependence on v_ret on the return time limit is a counterintuitive result, but it can be understood as follows. There is a minimum kick speed that is required to displace the hole beyond 0.1r₂₀₀, which represents lim_t→0v_ret; and there is a maximum kick, lim_t→∞v_ret, beyond which the BH remains completely unbound from the halo, even in the presence of drag. v_ret(t) then is a function of time that is always between these two extreme values. However, owing to the high central density of our gas-dominated halos, the difference between these two limits is small, ∼10%. Since this difference is smaller than other model uncertainties, such as those stemming from discrepancies from the actual density profile of high-redshift DM halos (e.g., clumping or triaxiality) and the merger tree prescription, we simply use the retention velocity computed for the approximate median time limit, 0.1t_Hub.

Finally, for simplicity, we treat the halo as a static mass distribution during each recoil event. That is, we ignore the cosmological evolution of the DM potential, and we assume that the recoiling BH does not affect the medium by clumping or heating it. Note, however that the latter feedback may play a nontrivial role in real systems, since the kinetic energy of a recoiling hole can be comparable to the gravitational binding energy of the entire host halo and can be expected to cause significant disruption of the surrounding matter.

We turn now to our prescription for the accretion rate of the BHs in our model. Of particular interest is the possibility that the gravitational recoil effect will significantly limit the ability of kicked BHs to accrete gas, by displacing them into low-density regions for prolonged periods, and/or by limiting through high relative velocities the amount of gas that can be gravitationally captured. One can imagine a scenario in which an SMBH whose progenitors have survived numerous kicks but have spent long episodes in underdense regions may have a final mass much less than that predicted by simple Eddington growth. We therefore follow the accretion rate self-consistently, as the recoiling holes proceed along their radial orbits. Specifically, in our models a BH embedded in gas is assumed to undergo standard Bondi–Hoyle–Littleton (BHL) accretion,

$$\begin{equation} \dot{m}(r,v) = \frac{4\pi G^{2}\rho _{b}(r)}{(c_{s}^{2}+v^{2})^{3/2}}m^{2}. \end{equation} \tag{ 15 }$$

The accretion rate is capped at the Eddington rate,

$$\begin{equation} \dot{m}=\frac{1-\epsilon }{\epsilon }\frac{m}{t_{\rm Edd}}, \end{equation} \tag{ 16 }$$

where t_Edd = 44 Myr and is the radiative efficiency, for which we assume = 0.10.

If the gas density is too low, or if the sound speed or its velocity with respect to the gas disk is too high, a BH may not be able to accrete at the Eddington rate even if it is close to the center of a halo. Because the BHL rate is proportional to m², an underfed BH with initial mass m₀ will eventually reach the Eddington accretion rate (∝m) at a threshold mass

$$\begin{equation} m_{\rm Edd}=\frac{1-\epsilon }{\epsilon }\frac{c_{s}^{3}}{4\pi G^{2}\rho _{b\;} t_{\rm Edd}} \approx 3500\left(\frac{c_{s}}{4\;\mathrm{km}\;{\rm s}^{-1}}\right)^{3}\! \left(\frac{\rho _{b}}{\,M_{\odot }^{-3}}\!\right)^{-1} \!\! \!\!M_{\odot }, \end{equation} \tag{ 17 }$$

where 4 km s⁻¹ is the isothermal sound speed for an ionized hydrogen gas at 1200 K. The typical central density of a 1200 K halo with a TIS gas profile is ∼5 × 10⁻³(1 + z)^7/3 M⁻³_☉. Sub-Eddington accretion rates are not an issue for BHs in a power-law gas profile, as the steep profile provides a sufficient central density for immediate Eddington accretion.

The time it takes for a BH with m₀ < m_Edd to reach this threshold mass, assuming that the local density remains constant, is

$$\begin{eqnarray} t_{\rm crit} &=& \left(\frac{m_{\rm Edd}}{m}-1\right) \frac{\epsilon }{1-\epsilon } t_{\rm Edd} \approx 2.7\left(\frac{m_{0}}{100\,M_{\odot }}\right)^{-1} \nonumber\\ &\times&\left(\frac{c_{s}}{4\;\mathrm{km}\;{\rm s}^{-1}}\right)^{3} \left(\frac{\rho _{b}}{\,M_{\odot }^{-3}}\right)^{-1} {\rm Gyr}; \end{eqnarray} \tag{ 18 }$$

here, t_crit ∼ a few 100 Myr for 100 M_☉ halos embedded in 1200 K TIS halos at z ≳ 20, and t_crit ≳  Gyr for z ≲ 14. This means that for TIS gas profiles, seed holes will spend a significant fraction or all of the available time prior to z ≈ 6 accreting below the Eddington rate.

The difference between BHL and Eddington accretion rates as they relate to BH growth is also discussed in VR06. However, in that paper the context is for the direct formation of m > 10⁴ M_☉ intermediate-mass BHs through super-Eddington BHL accretion. We here adopt the opposite extreme assumption, i.e., that the BH radiates efficiently at all times, and its accretion rate obeys the Eddington limit. The BHL rate can then initially be sub-Eddington in TIS halos, owing to the low BH mass and low gas density (the baryon density required to fuel BHL accretion at the Eddington rate was also discussed by Turner 1991).

To illustrate the impact of extended sub-Eddington growth phases, we perform a simple analytical calculation. Suppose that a BH is formed with mass m_seed at redshift z in a halo with virial temperature 1200 K, and that the local gas density is held constant at the value when the BH was formed. Note that for this exercise, we assume that gas density is constant even as the halo around the BH is growing by merging with other halos—in other words, we assume these mergers deliver gas to the nucleus containing the original BH, roughly maintaining a constant density at the Bondi radius around the BH. Figure 4 shows the maximum possible mass that can be attained by such a BH growing in isolation through gas accretion alone, assuming that the accretion rate is determined solely by Equations (15) and (16) and that accretion is not supply limited. If the host halo has a steep power-law cusp, the accretion rate is Eddington throughout regardless of when the seed BH is formed. However, if the central fuel density is low, then it is possible for the local BHL accretion rate to be significantly sub-Eddington initially. In such a scenario, the earliest-forming seeds are the only ones able to reach the Eddington rate; the late-forming seeds are unable to reach the Eddington rate before z = 6. In this scenario, the late-forming seeds, which are easily identifiable by the drastically shallower growth slope in the figure, cannot grow rapidly enough to contribute to the SMBH population. Note that assuming a constant gas density in this calculation gives an optimistic accretion rate for the TIS case, as in those profiles the central gas density generally decreases with Hubble expansion and significantly reduce the BHL rate.

Zoom In Zoom Out Reset image size

It is computationally expensive to numerically integrate the individual orbits and accretion histories of every recoiling BH in our simulations. We therefore tabulate the accretion growth in Eddington units during the first 0.1t_Hub of the orbit, and describe the results in a fitting formula in the same manner as we have done for the retention velocity. Given a specific prescription for the baryon distribution, we tabulate across four relevant variables in the following ranges: 10⁵ M_☉ < M < 10¹⁵ M_☉, 10⁻⁶ < m/M < 1, 5 < z < 40, and 0 < v_kick/v_ret < 1.

In the absence of a kick, and if the accretion rate were always at the Eddington limit, the SMBH mass in a given halo at z ≈ 6 is easily approximated by

$$\begin{equation} m_{\rm SMBH}\approx N_{\rm seed}\; \exp \left(f_{\rm duty}\frac{\epsilon }{1-\epsilon }\frac{t_{{\rm seed},6}}{t_{\rm Edd}}\right)m_{\rm seed}, \end{equation} \tag{ 19 }$$

where m_seed is the seed BH mass, N_seed is the number of seeds in the merger history of the halo, t_seed,6 is the available time between the typical seed formation time, and z ≈ 6 and f_duty is the time-averaged duty cycle for accretion. Equation (19) represents the ideal, maximally efficient scenario for SMBH assembly, and we can use it to effectively measure the cumulative impact of underfed accretion, recoil-induced ejection, and other factors that limit the assembly efficiency.

The measurements of clustering of quasars in the SDSS suggest that the duty cycle of active (luminous) accretion increases steeply with redshift at 3 ≲ z ≲ 6, with the most active quasar BHs at z ≈ 6 showing 0.6 ≲ f_duty ≲ 0.9 (Shen et al. 2007; Shankar et al. 2008). We therefore adopt duty cycles of ⩾0.6. Although it is likely that SMBHs regulate their own growth through feedback mechanisms, we do not address such scenarios a priori. To keep our models as simple as possible, we will only impose the loosest upper limit on the SMBH accretion rate: they must not accrete more mass than there the total mass of baryons in the host halo. However, we will discuss an alternative ad hoc model below, which is able to reproduce the well-known relation between SMBHs and the velocity dispersions of the bulges of host galaxies (the m–σ relation).

Explicitly constructing the SMBH mass function at z = 6 over a wide mass range is computationally intractable. The host halo mass inferred from the observed quasar space density from the z ∼ 6 quasars is several ×10¹² M_☉ (e.g., Fan 2006). For every halo with this mass, there are ∼10⁷ halos with 10⁸–10⁹ M_☉. Blindly calculating the SMBH mass for every halo with M > 10⁸ M_☉ using our trees-plus-orbits algorithm would be prohibitively expensive computationally.

We therefore carry out a piecewise calculation of the SMBH mass function that is computationally tractable. The procedure is as follows: (1) we group the halo population into logarithmic mass bins of size x < log₁₀M < x + Δx; (2) for each bin, we select ≳10²–10⁴ individual Monte-Carlo-generated halos with masses randomly generated from the Press–Schechter distribution at z = 6; (3) we simulate the BH population for each such halo using our trees-plus-orbits algorithm, and assume the resulting sample is representative of all z = 6 halos in the mass bin; (4) for each bin, we multiply the sample by the appropriate weight to construct the entire Press–Schechter halo mass function, ∫^{x + Δx}_xdN/dln Mdln M; and finally (5) sum the contributions from each bin. The result is a Press–Schechter distribution of halos with M > 10⁸ M_☉ at z = 6, with a statistical representation of the corresponding SMBH mass function. The bins used and their relevant properties, including the number of Monte Carlo halos that were cloned to populate the full mass function, are listed in Table 1. This numerical shortcut is not used for the most massive halos. Forty halos are expected above 10^12.85 M_☉, and these are simulated individually.

This method allows a fast calculation of the BH mass function, with the caveat that there must be enough BHs in each bin to keep statistical uncertainties to a minimum. Unfortunately, this is not always preventable for models with extremely low f_seed, and the reader will notice statistical noise in the results of such models. In some cases, we increase the halo sample by a factor of 10 in an attempt to reduce the errors.

In all, our simulations represent the Press–Schechter population of DM halos in a comoving volume of ≈280 Gpc³, roughly equal to the comoving volume that was probed by the SDSS survey between z ≈ 5.7 and 6.4. We have simulated a total of ≈1.17 × 10⁵ DM halos (see the column N_sim in Table 1), and through the procedure described above extrapolated the results to represent the ≈3 × 10¹² Press–Schechter halos (with M > 10⁸ M_☉) expected to be present in the 280 Gpc³ volume.

As stated in the Introduction, our primary goals are (1) to understand the possible ways in which the >10⁹ M_☉ SMBHs may have been assembled at redshift z > 6, and (2) whether the LISA event rates are sufficiently different in competing models so that one can disentangle the physical assembly scenario from the LISA data stream. As a first step toward these goals, we would like to survey all feasible combinations of the physical assembly parameters, understand the impact of each model parameter, and look for corresponding give-away features in the predicted LISA stream.

Although a broad simulation survey is beyond the scope of this paper, we have undertaken a cursory tour of the basic parameter space. We begin by considering two basic seed models: 100 M_☉ seeds forming as remnants of Pop-III stars in minihalos when they first reach the virial temperature T_vir ⩾ T_seed = 1200 K and 10⁵ M_☉ seeds forming as a result of direct collapse of gas in halos when they first reach T_vir ⩾ 1.5 × 10⁴ K. For each case, we consider halos with an NFW DM component and a gaseous component with either a cuspy ρ_gas ∝ r^−2.2 power law or a corey TIS profile. The TIS models consistently failed to produce SMBHs by z = 6, with typical maximum BH masses of ≲10³ M_☉. In those models, the central gas densities are too low to allow for prolonged episodes of accretion near the Eddington rate, as noted in Section 2.5.

For each type of seed, we vary three sets of parameters: (1) the seeding fraction 10⁻³ ⩽ f_seed ⩽ 1, i.e., the probability that a halo reaching the critical temperature will form a seed BH; (2) the time-averaged accretion rate in Eddington units, characterized by the duty cycle f_duty, for which we use f_duty = 0.6, 0.8, and 1.0 (note that f_duty and the radiative efficiency are degenerate in our prescription; see below); and (3) whether at the time of their merger, the BH spins are randomly oriented or aligned with the orbital angular momentum of the binary. The spin magnitudes are chosen from a uniform random distribution 0 < a_1,2 < 0.9 (Schnittman & Buonanno 2007). We do not track the evolution of BH spins in our models. While more rapidly spinning BHs are capable of accreting at higher efficiency, we neglect this effect and assume a global efficiency coefficient in our models. Of the four main ingredients of the SMBH assembly introduced in Section 1, we here therefore vary three: the seed rarity, the accretion rate, and the recoil dynamics. For the fourth, the merger rate, we have simply assumed that BHs merge when their parent halos do, as extreme-mass BH binaries are rare in our simulation, given the threshold we have imposed on the mass ratio for halo mergers. We will call the above our fiducial set of parameters.

For each realization, we compute the mass function of BHs at z = 6 with m > 10⁵ M_☉; the rate of events in the LISA observational mass range, ∼10^4–7/(1 + z) M_☉; the fraction of DM halos hosting a central massive BH; and m/M, the ratio between the mass of the SMBH and its host DM halo, which serves as a proxy for the m–σ relation.

We begin with Figure 5, showing the mass function for the m_seed = 100 M_☉ Pop-III seed model. This and the companion figures are organized with different values for f_seed in different columns, the two spin prescriptions in different rows, and the different duty cycles in different line styles (and different colors, in the online version). This will be the default organization of our eight-panel figures unless noted otherwise.

Zoom In Zoom Out Reset image size

The numbers in the upper right-hand corner of each panel represent the total SMBH mass density, log₁₀[ρ_•/( M_☉ Mpc⁻³)], for all BHs with m > 10⁵ M_☉. Because of statistical fluctuations, for multiple model realizations with identical parameters this value can vary by ≲10% for f_seed ≳ 10⁻¹, and as much by a factor of a few for f_seed ≲ 10⁻³. For each of these simple models, this density is exceedingly high compared to the SMBH density of the local universe. For the present discussion, we will overlook this point as we address the effects of the various model factors; we will introduce the additional constraint from ρ_• in subsequent sections.

The most stringent observational requirement for the high-mass end of the z ≈ 6 SMBH mass function comes from the SDSS observation of the z ≈ 6.4 quasar J1054+1024, which has an inferred mass of ∼4 × 10⁹ M_☉. Since this object was detected in a region ∼10% of the sky, we estimate that ≳10 similar objects may exist at z ∼ 6. We have adopted this as a rough indication of the lower limit of the SMBH mass function at redshift z = 6, and represent this limit by the upper right quadrangle in each of the panels, delineated by the red dashed lines. Note that these lower limits are conservative, since they do not require the presence of any additional SMBHs with comparable mass that are "off." Since high values for the duty cycle—near unity—are suggested by quasar clustering measurements (as discussed above), and are also required for growing the most massive SMBHs (as we find below), this correction is only of order a factor of ∼2. As seen in the figure, the high-mass end of the SMBH mass function is so steep that increasing the lower limit on the required SMBH space density by even ∼2 orders of magnitude would make little difference to our conclusions.

The first conclusion one can draw from this figure is that several combinations of parameters can produce SMBHs massive enough to power the quasar J1054+1024. We will return to this and other observational constraints in the following subsection, but let us for now focus on understanding the effects of our various parameters and their combinations. In general, the effect of varying each of the parameters is relatively straightforward to understand. Increasing the accretion rate, increasing the seed occupation fraction, and aligning the spins all tend to result in more massive BHs. Note that the accretion rate in Eddington units, f_duty, is degenerate with the accretion efficiency, as $\dot{m}\propto f_{\rm duty}\times (1-\epsilon)/\epsilon$. For example, a model with = 0.15 and f_duty = 1.0 is equivalent to = 0.10 and f_duty = 0.63. We have used = 0.10 throughout the results presented here. With this value for the efficiency, it is possible to build the most massive SDSS quasar SMBHs by z ≈ 6, starting with 100 M_☉ seeds. If the efficiency is ≈0.15, however, it is only possible to build the SMBHs in question with the most optimistic assumptions: f_duty ≈ 1, f_seed ≳ 0.1, and spin alignment would all be required.

We note two nontrivial observations to be made from Figure 5. First, if the seed fraction is low, the spin orientation has a minimal effect on the BH mass function. This is because if seeds are extremely rare, they are likely to grow in isolation for much of their existence along with their host halos. The first mergers are not likely to occur until the gravitational potentials of their host halos are deep enough to retain them from any recoil event. Second, the increase in the SMBH abundance from increasing the seeding fraction has a tendency to plateau. This is the reverse situation compared to the low f_seed limit: if seeds are very common, they are likely to experience multiple BH mergers very early in the merger tree, when their masses are still comparable, and ejections become very common. As f_seed increases, assembly becomes increasingly inefficient at early times. In models where the seeding halo temperatures are lower, the efficiency saturates at lower values of f_seed. Furthermore, in models with the shallower TIS gas profiles, we find that increasing the occupation fraction beyond a certain "sweet spot" value rapidly decreases the final SMBH masses. The reasons for this are that (1) the seed BHs in these models barely grow, making near-equal-mass (q ∼ 1) BH mergers, and therefore large kicks, much more common and (2) the gas drag is reduced, making it easier to kick holes out of these halos. We conclude that arbitrarily increasing the number of seed holes contributing to the assembly process is not necessarily an efficient solution to the SMBH assembly problem. In fact, excessive seeding can lead to a different conflict with observations—overproducing the mass density in lower mass SMBHs—to which we will return in the following subsection.

Let us turn now to Figure 6, which shows the BH occupation fraction in M > 10⁸ M_☉ DM halos at z = 6, and the BH-to-halo mass ratio (which here serves as a proxy for the m–σ relation; see, e.g., Ferrarese 2002) for each of the models. Here, we have arbitrarily defined our occupation fraction to mean the fraction of DM halos that host a BH with a minimum mass of 10⁴ M_☉, as we are interested in the population of nuclear BHs and not stellar remnants. Menou et al. (2001) have shown that, in the low-redshift universe, the fraction of DM halos hosting an SMBH will approach unity, regardless of the initial occupation fraction at earlier times. We find that at z ≈ 6, the occupation fraction in halos with M ≲ 10¹¹ M_☉ can still be significantly below unity, depending primarily on f_seed. In principle, a future survey that can determine the quasar luminosity function (LF) at z = 6 to several magnitudes deeper than the SDSS could look for this drop in the occupation fraction, since it will produce a flattening of the LF at magnitudes below some threshold. The spin prescription has a noticeable effect on the z = 6 occupation fraction for f_seed ≳ 10⁻². On the other hand, the duty cycle essentially only affects the mass of the BHs, and not their presence or absence, and so has a minimal effect on the occupation fraction, within the parameter range shown. Note that we do not expect to reproduce the m–σ relation in our simulations, as we employ simple and z-independent accretion and seeding prescriptions, and our models have no feedback to enforce the relation. Instead, the m/M relation should be taken as a sanity check that we are producing a physically viable BH population. Note that in some of the models shown in Figure 5, the BHs grow much larger than the m/M relation observed in nearby galaxies. In particular, as the figure reveals, SMBHs in the lower mass (∼10⁸ M_☉) halos in the models with f_duty = 1 tend to consume most of the gas in their parent halos, clearly an improbable result.

Zoom In Zoom Out Reset image size

For an explicit comparison, alongside the m/M relation produced by our models, we have plotted the value expected based on measurements of the m–σ relation in local galaxies. Specifically, we show the expression from Wyithe & Loeb (2003),

$$\begin{eqnarray} m&=&10^{12}\epsilon _{0}\left(\frac{M}{10^{12}\,M_{\odot }}\right)^{\gamma /3} \left(\Omega _{0}h^{2}\right)^{\gamma /6}(1+z)^{\gamma /2} \nonumber\\ &&\approx 1.16\times 10^{7}\left(\frac{M}{10^{12}\,M_{\odot }}\right)^{\gamma /3} (1+z)^{\gamma /2}. \end{eqnarray} \tag{ 20 }$$

We adopt their parameter choices of ₀ = 10^−5.4 and γ = 5. This expression satisfies the SDSS constraints at the high-mass end of the mass function. It also agrees well with the relation suggested by Ferrarese (2002) for SMBHs in the local universe, m ∼ 10⁷(M/10¹² M_☉)^1.65. As the figure shows, our predicted m/M relation tends to have the opposite slope than the one inferred from the observed m–σ relation: in our results m/M decreases with mass or stays roughly constant as M increases, but this is the opposite of the empirical trend. This is due mainly to the difference in the growth rates of halos and holes: our simple prescription for steady, exponential accretion for the BHs can significantly exceed the growth rate of DM halo masses due to the accretion of unresolved low-mass halos in the EPS merger tree. As a result, in some cases, the host halo mass predicted for the z = 6 quasars is as low as 10¹¹ M_☉, which is an order of magnitude lower than would be predicted from the extrapolation of the locally measured m/M relation, or from the inferred space density of the host halos (e.g., Haiman & Loeb 2001). However, any extrapolation of the m–σ or m/M relation to high redshift, and the masses of halos that host the brightest z > 6 quasars, at present, have large uncertainties, and do not robustly preclude such low values. As will be discussed below, the overly rapid growth of SMBHs in this suite of models motivates modifications to the modeling, including a model in which the extrapolated m/M relation holds by assumption.

Figure 7 shows estimates for the LISA event rate, calculated for all binary mergers satisfying 10⁴ M_☉ < (m₁ + m₂)(1 + z) < 10⁷ M_☉ as (see, e.g., Menou et al. 2001),

$$\begin{equation} \frac{d^{2}N}{dz\; dt}=4\pi c \; d_{\rm com}^{2}(z) \frac{\Delta N}{\Delta z \; \Delta V}, \end{equation} \tag{ 21 }$$

where ΔN is the number of SMBH merger events in the tree in a time step Δz and a simulated comoving volume ΔV, and d_com is the comoving distance. Although there is a mild dependence on the duty cycle/accretion rate and the kick prescription, it is evident that for our assembly models f_seed has the greatest effect in setting the rate of SMBH binary mergers detectable by LISA. Because the initial merger rates scale as f²_seed, the measured event rate is extremely sensitive to the BH number population. Since the merger activity peaks near z ≲ 10, LISA should be able to measure the masses of most of these SMBH binaries to high prevision (Hughes 2002; Arun et al. 2008). Although the raw detection rate—without any information on BH spin or mass ratio—alone will be richly informative, the ability to determine the source masses with high fidelity should be very useful in further constraining the growth rate of the first SMBHs, and breaking degeneracies between the allowed parameter combinations. The observed distribution of binary masses as a function of redshift will provide direct snapshots of the mass function and shed light on its evolution, independently from quasar luminosity surveys (e.g., Yu & Tremaine 2002; Willott et al. 2006).

Zoom In Zoom Out Reset image size

Another family of assembly models that has been frequently discussed in the literature is the so-called "direct-collapse" model, wherein BHs with m ∼ 10^4–6 M_☉ are formed rapidly from gas cooling via atomic H in halos with virial temperature T ≳ 10⁴ K (Oh & Haiman 2002; Bromm & Loeb 2003; Volonteri & Rees 2005; Begelman et al. 2006; Spaans & Silk 2006; Lodato & Natarajan 2006). We simulate such a family of models, for the same set of parameters f_seed and f_duty and the same spin alignments. We choose m_seed = 10⁵ M_☉ and T_seed = 1.5 × 10⁴ K. We show the mass functions, the hole–halo relations, and the LISA rates (the same information as in Figures 5, 6, and 7) in Figures 8, 9, and 10. Although the main differences are all fairly intuitive, it is instructive to address how the direct-collapse models differ from the Pop-III seed models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, from Figure 8 we see that it is much easier to construct more massive SMBHs owing to the larger seed masses. In fact, the models with f_seed ≳ 0.1 become unphysical, since the SMBHs in these models would exceed the baryon mass (Ω_b/Ω_m)M of their parent halos. The second point is that the mass function is very differently distributed between the Pop-III and direct-collapse scenarios. The reader will note that for each set of parameters, the overall SMBH density does not differ significantly between the corresponding seed models. However, this is deceiving as the mass function is clearly different, with the direct-collapse seeds resulting in a more "top heavy" SMBH population. For the range of parameters surveyed, the Pop-III model has ≳90% of the SMBH mass in the m < 10⁷ M_☉ range if D = 0.6, compared to ≲30% for the D = 0.6 direct-collapse seed models. For D = 1.0 and random spin alignment, ∼2%–5% of the SMBH density resides in the most massive m > 10⁹ M_☉ BHs if the seeds are Pop-III; for the same parameters, ∼60%–70% of the mass is in the billion-plus solar mass BHs in the corresponding direct-collapse models. Note that the total BH mass density remains extremely high; again, we will address this point shortly.

Third, there is a slightly weaker dependence on the spin orientation of the BH binaries. This is most apparent by comparing the z = 6 occupation fractions in Figures 6 and 9, and is due to the deeper potentials of the host halos in the direct-collapse case. Fourth, even though the m/M relation continues to have a slope opposite to the locally observed trend, there appears to be a break in the relation at M ∼ 10⁹ M_☉, accompanied by a drop in the occupation fraction. This is due to the simple fact that halos below this mass threshold cannot have many T > 1.5 × 10⁴ K progenitors and the model similarly prohibits intermediate-mass BHs. This cutoff contributes to the top-heaviness of the mass function for these models. Fifth, LISA rates are lower by 1–2 orders of magnitude than in the Pop-III seed models, because there are fewer 1.5 × 10⁴ K halos than 1200 K halos (another factor is that the seed BHs are already born with a mass near the middle of LISA's logarithmic mass range, so they spend only about half the time in the LISA band, compared to the Pop-III seeds). It is worth noting, in particular, that it is possible to build the SDSS quasar BHs in ways that produce no detectable GW events for observation with LISA beyond z ∼ 6 (in contrast, in the successful Pop-III models, a minimum of a few events are predicted). In such scenarios, SMBHs are extremely rare until z ∼ 6, at which point the SMBH occupation fraction will evolve toward unity fairly rapidly, as described by Menou et al. (2001).

Finally, in Figure 11 we present the mass function and the LISA detection rate in five variants of a fiducial f_duty = 0.65, f_seed = 1, aligned spin model with Pop-III BH remnant seeds. We choose these values because they just barely are able to match the abundance of the most massive SDSS quasar BHs, and because they produce the most BH mergers, and thus the effects of varying the other recoil-related parameters will be the most discernible. We show the fiducial model in bold. In the modified models, we (1) allow the gas density profile to be shallower, with a r⁻² power law—this is to allow for the possibility that the gas surrounding the BH has not cooled and condensed to high density (dark blue curves); (2) require the BH spins to be near-maximal at a_1,2 = 0.9, instead of choosing them randomly from the range 0.0 to 0.9—this is to allow for the possibility that all SMBHs at high z are rapidly spinning, e.g., due to coherent accretion (green curves); (3) allow halos of all mass ratios to merge, where previously we had considered a halo with mass less than 1/20th that of its merger companion to become a satellite (yellow curves); (4) assume that the Pop-III seed progenitors are not able to blow out the gas in the host minihalo (red curves); and (5) ignore the effects of accretion suppression due to episodes of recoil-induced wandering, and instead assume that all BHs accrete at f_duty unless ejected (light blue curves).

Zoom In Zoom Out Reset image size

We also ran models with a TIS profile for the gas component. We found that these models always failed to produce any SMBHs above 10⁶ M_☉ by redshift z = 6 due to the initial phase of sub-Eddington accretion of seed BHs, and their results are not shown in Figure 11. This implies that SMBHs must be continuously surrounded by dense cores of gas that was able to cool at the centers of DM halos—feeding holes with the low-density gas in DM halos whose gas was unable to cool does not allow for high enough accretion rates (Turner 1991 emphasized the same issue for the growth of z ∼ 4 quasar BHs).

The results shown in Figure 11 give insight to the importance of the assumptions that went into our models. In particular, of all the parameters that we have varied for fixed f_seed and f_duty, the most significant for the SMBH mass function, by far, are the spin orientation and the limit on the halo mass ratio for timely mergers. Both of these have a similar effect of increasing the number of BHs, especially at the massive end. The former result—that maximizing the spin increases the SMBH mass function—may seem surprising at first, since generally large spins imply larger kicks. However, this is not always the case, as can be understood from Equations (5)–(8). Under the assumption that both spins are aligned with the orbital angular momentum, v_∥ = 0, and v_⊥ ∝ (a₁ − qa₂) is maximized for unequal spins (i.e., a₁ = 1 and a₂ = 0 for a typical q ∼ 1); setting a₁ = a₂ = 0.9 therefore eliminates the largest kicks, and allows more BHs to be retained. Comparatively smaller differences are visible in the mass function for the other model variants. Figure 11 also shows that adding in the unequal-mass halo mergers and increasing the spins affect the LISA event rates differently: the former adds new events mostly at z ≲ 10, where unequal-mass mergers are more common, whereas increasing the spin mostly adds new events at z ≳ 10. Interestingly, ignoring the blow-out has little impact on the SMBH mass function at z = 6, but it does shift the LISA events to higher redshifts, especially at z ≳ 10. Figure 11 suggests that the LISA event rate can be useful in disentangling these three effects.

Perhaps the most surprising inconsequential variation is ignoring episodes of reduced accretion due to the BH wandering in low-density outskirts following recoil. The reason for this is simply that lengthy oscillating orbits are relatively rare if the central gas density is high; the BHs either return quickly, or are ejected (by assumption). Recall that our definition for a "retained" BH is that the recoiled hole must return to within 1/10 of the halo's virial radius within 1/10 of the Hubble time. For low kick speeds, the BH does not get very far, because the gas provides both a steep gravitational potential barrier and a strong dissipative sink. The orbit will thus be rapidly damped, with only very brief periods of underfeeding. If a BH is kicked hard enough to reach the outskirts of the halo, there is very little dissipation there, to tug it toward the center or to further damp its oscillation. Since the radial velocity at this point is low, it is likely to remain outside 1/10 of the virial radius for a significant period. The bottom line is that the range of initial kicks that would take the BH to the low-density outskirts to cause significant underfeeding, while still allowing it to return quickly enough to be "retained" is simply negligibly small.

Blecha & Loeb (2008) recently performed a more detailed analysis of the orbits of recoiling SMBHs that include a multi-component halo mass distribution and three-dimensional orbits. They report that recoil velocities of between 100 km s⁻¹ and the escape velocity lead to significant suppression of the accretion rate, with SMBHs accreting only ∼10% of its initial mass over 10⁶–10⁹ years of wandering through the halo. We note that (1) in our simulations we ignore the longest-wandering BHs through our prescribed retention threshold; (2) typical kick magnitudes for SMBHs are lower than 100 km s⁻¹ in our simulations, a point we explain below; and (3) their prescription of the baryon distribution results in a lower central density than our models, as in that paper they are concerned with typical galaxies at low redshift, and not with minihalos and protogalaxies. We conclude that prolonged periods of wandering and underfeeding are unlikely to have a significant effect on the mass function of high-z SMBHs as a whole. Oscillations may play a more prominent role in the growth history of SMBHs if large kicks are more common (and the retained holes tend to be marginally retained), if the halo gravitational potential is more shallow, or if the effect of dynamical friction on BH orbits is less than what we have considered in this paper. Also, halo triaxiality (Blecha & Loeb 2008) and/or a clumpy mass distribution (Guedes et al. 2008) could increase the time that kicked holes take to return to the halo center (or they may never return). In principle, this may increase the impact of these oscillations. However, in practice, the inner, baryon-dominated regions of the galaxies are likely close to smooth and spherical. The BHs that are not ejected in our models typically do not leave these regions and so will not be subject to large changes in their orbits from these effects.

The results presented thus far seem to paint a relatively simple set of relevant parameters for SMBH assembly. There is the seeding function f_seed, which governs the BH merger rate and therefore to a large extent the LISA event rate. The event rate also depends on the time-averaged accretion rate, parameterized here by the duty cycle f_duty, and the initial seed mass M_seed, as these set the evolution of the observable mass spectrum. The seeding prescription and f_duty determine the mass function almost entirely if f_seed ≪ 1. If f_seed ≳ 0.1, then the spins of the BH binary play a role in setting the recoil speeds and the subsequent evolution.

Once typical spin and mass parameters of merging SMBHs are determined by LISA, either by direct measurement, or perhaps by extrapolating from detections at lower redshift, combined with the event rate and what is known about the upper end of the mass function, this information is likely sufficient to give at least a strong indication on the typical SMBH growth rate and initial seed mass.

Our simulations above also confirm a result reported by VR06, namely that SMBHs form primarily through repeated mergers of the most massive SMBHs merging with less-massive (SM)BHs. This is because the gravitational rocket speeds decrease rapidly as the mass ratio q decreases. The first few BHs that "outweigh" their neighbors—be it through being endowed with more mass at birth, accreting faster or being fortunate enough to survive the first mergers—will be less likely to be ejected from their host halos. This survival advantage becomes a runaway effect, as each subsequent merger will result in a lower value of q for the next merger.

Now that we have a first-glance grasp of the assembly parameters and their most basic observational characteristics, we can turn to identifying actual candidate models for the formation of m ≳ 10⁹ M_☉ SMBHs before z ≈ 6.

The suite of models discussed above has demonstrated that there are several feasible ways to build the SMBHs. These models have so far focused only on the number density of m ≳ 10⁹ M_☉ SMBHs, and ignored indirect constraints that exist on the mass function at lower masses. In particular, the total mass density in SMBHs with masses in the range 10⁶ M_☉ ≲ m ≲ 10⁹ M_☉ in the local universe has been estimated by several authors, who find ∼4 × 10⁵ M_☉ Mpc⁻³ (to within a factor of ∼2; see, e.g., the recent paper by Shankar et al. 2009a and references therein). Furthermore, a comparison of the locally observed mass density with the mass density inferred from accretion by the evolving quasar population suggests that ∼90% of the total local SMBH mass density is attributable to quasar accretion. This implies that the total SMBH mass density increased by a factor of ∼10 between z ∼ 6 and the local universe (Yu & Tremaine 2002; Haiman et al. 2004; Shankar et al. 2009a), which then places an indirect constraint on the SMBH mass function, down to m ∼ 10⁵ M_☉, at z = 6.

To be specific, we set the following upper limit on the expected total z ≈ 6 SMBH mass density in m ≳ 10⁵ M_☉ SMBHs:

$$\begin{equation} \rho _{\rm SMBH, 5+}(z\,{=}\,6)\,{\sim}\,0.1\,{\times}\,\rho _{\rm SMBH, 6+}(z\,{\approx}\,0)\,{\sim}\,4\,{\times}\,10^{4}\,M_{\odot}\;\mathrm{Mpc}^{-3}, \end{equation} \tag{ 22 }$$

that is, we assume that the total mass density of SMBHs with mass m > 10⁵ M_☉ at z = 6 is at most 10% of the total inferred mass density of SMBHs with mass m > 10⁶ M_☉ in the local universe. The major caveat to making such an expectation is that we assume that the low-mass end of the BH mass function has grown by a factor of 10, and that high-redshift quasar luminosity surveys have sufficiently accounted for selection effects of any SMBHs that may be hidden by inactivity or by being too dim.

The analysis that follows below is similar to that of Bromley et al. (2004), who considered the upper limit to the z ≈ 6 SMBH mass density in weighing the plausibility of various z = 6 quasar BH assembly models. The main difference is that Bromley et al. (2004) did not consider the gravitational recoil effect. Adding in this recoil makes the problem significantly worse. This is because the recoil necessitates more optimistic assumptions in order to build up the ∼10⁹ M_☉ SMBHs, which tends to increase, by a large factor, the predicted number of lower mass ∼10⁶ M_☉ SMBHs, which arise later, and whose growth is therefore less sensitive to the kicks. In other words, the kicks preferentially suppress the abundance of the most massive SMBHs which arise from the earliest seeds; as a consequence, we predict steeper SMBH mass functions than the models considered in Bromley et al. (2004). We performed a series of assembly calculations explicitly without recoil, and find that these indeed produce mass functions with a flatter slope and a higher normalization, owing to greater SMBH masses at the high-mass end, and a higher occupation fractions at all masses. If there are no kicks, less optimistic assumptions are required to produce the SDSS quasars, and the overall mass density is lower in no-kick scenarios that produce the minimum number of ≳10⁹ M_☉ SMBHs.

The basic result we find is that among the models presented thus far, all of those that match the SDSS abundance of the most massive SMBHs overshoot the above value by 2 or more orders of magnitude. One possible solution to this apparent problem is that there is no problem at all: we have simply set our constraint too severely. Perhaps not all 10⁵ M_☉ SMBHs at z ≈ 6 have evolved to become m > 10⁶ M_☉ BHs in galactic nuclei in the local universe. Still, there is no obvious way to "hide" these low-mass SMBHs locally, and the overprediction in the SMBH densities in our simulations are very large: we find that in our fiducial models, we typically need to reduce the population of SMBHs in the mass range 10^5–7 M_☉ by a factor of ≳100, and those in the 10^7–9 M_☉ range by a factor of ≳10.

One logical solution, and one that has been suggested by Bromley et al. (2004), is to simply stop making seeds below some cutoff redshift. The earliest seeds contribute the majority of the total mass of the most massive SMBHs, and late-forming seeds tend to end up in lower mass SMBHs. In principle, then, by introducing a cutoff redshift for seed formation, one can suppress the low-mass end of the mass function, while still allowing for the formation of the most massive SMBHs. Such a model would also be in line with our physical understanding of seed BHs. We know Pop-III stars stopped forming relatively early on in the universe, with the halt in production being due to trace metal contamination (e.g., Omukai et al. 2008), radiative feedback from the UV and/or X-ray background during the early stages of re-ionization (e.g., Haiman et al. 2000), the higher turbulence of gas in the centers of later halos, or a combination of these factors.

This proposed solution amounts to keeping f_duty a constant while allowing f_seed to evolve with z; in the simplest case, as a step function dropping to f_seed = 0 below some redshift. Figure 12 shows the fractional contribution (dM/dz_seed)/M from seeds forming at different redshifts to the final mass at z = 6 in three different bins of the z = 6 SMBH mass function, for two different models that are marginally able to assemble the SMBHs powering J1054+1024. The model in the left panel has f_seed = 10⁻³, while the model on the right has f_seed = 1. Note that contributions to each mass bin are normalized to integrate to unity, but the two lower mass bins ρ(10^5–7 M_☉) and ρ(10^7–9 M_☉) make up the vast majority of the total mass density. The formation epochs of the seeds contributing the majority of the mass of the different SMBHs are very distinct in the model with lower f_seed, but overlap significantly for f_seed = 1. What accounts for this qualitative difference in the assembly epochs? There are two ways to build 10⁹ M_☉ holes: accretion onto the earliest-forming holes that undergo few mergers and grow mostly in isolation, and the mashed-together product of many seeds that may have formed later. If the occupation fraction is high, one expects both populations to contribute, and therefore a wide distribution for dM/dz_seed. If the occupation fraction is low, SMBHs with many progenitors contributing become rare, and we find that the most massive SMBHs can only be formed from isolated seeds in the earliest minihalos, through few mergers, and so the dM/dz curves in these models are sharply peaked.

Zoom In Zoom Out Reset image size

Our task is to eliminate ≈99% of the BH mass in the lower mass bins, while leaving most of the m ≳ 10⁹ M_☉ holes unaffected. By simply examining the normalized dM/dz curves in Figure 12, one can see that simply cutting off seed formation at an arbitrary redshift for the f_seed = 1 model is not a viable solution to the overproduction problem, as there is no way to eliminate the lower mass SMBHs without eliminating a significant fraction of the 10⁹ M_☉ holes. Cutting off the seed production can produce successful mass functions for models with f_seed ≪ 1, but only if the seed cutoff occurs at very high redshifts, typically z ∼ 30. Essentially, the solution calls for a very brief and early period of seed BH formation, and very rare seed formation there after. An unfortunate generic consequence of this early cutoff redshift is that it quickly chokes off the LISA event rates.

We also simulated models where seed production continues beyond the cutoff redshift—with the same probability f_seed as in the minihalos at z > z_cut—in halos with T_vir > 10⁵ K. Such halos could continue to form BHs if heating by the UV background is the primary mechanism for seed suppression, as they are able to shield their central gas from the UV radiation (Dijkstra et al. 2004). We find that models with such a partial cutoff overproduce the SMBH mass density if f_seed ≳ 10⁻⁴. This result implies that either the initial occupation fraction is very low (it is still possible to make the SDSS BHs with this low f_seed; see Table 2 below), or else some other feedback beyond UV radiation, such as metal enrichment, stops seed BHs from forming in all halos at z < z_cut, even in the rare, more massive ones.

In short, the requirement in models in which the duty cycle is a constant, but seeds stop forming suddenly below some redshift, can be simply summarized as follows: the only way to build SMBH mass functions that satisfy observational constraints and indications at both the low-mass and the high-mass ends is from extremely rare seeds that form during a brief and very early epoch. We also note that these models are attractive because (1) there are physical reasons for the seeds to stop forming below some redshift, and (2) there is independent empirical evidence, from constraints on the re-ionization history from WMAP measurements of the cosmic microwave background polarization anisotropies, that the ionizing luminosity in high-redshift minihalos was suppressed by a factor of ≳10 (Haiman & Bryan 2006).

We present in Table 2 the parameters for four such successful models. While it is not computationally feasible to search the entire parameter space for such models, we present two typical examples for both the Pop-III-remnant and direct-collapse seed BH scenarios. f_seed is low in each of these examples, as we have argued above that they must be. Although we have listed the spin prescriptions, they are relatively unimportant because seeds are rare (see Section 3.1). Given a particular value for f_seed, the only free parameters are the accretion rate f_duty and the seed cutoff redshift z_cut. For the Pop-III models, we find in the range f_seed ⩾ 10⁻⁴ that seeds must typically stop forming at z_cut > 30, with a lower cutoff z_cut ∼ 20 for the direct-collapse models. An important conclusion is that in each of these models, GW events are too rare (<10⁻³/dz/yr for z < 30) for LISA detections beyond z ∼ 6 to occur within the mission's lifetime.

The mass density of ejected BHs is exceedingly low in each of the successful scenarios, <10⁻³ M_☉ Mpc⁻³. Because the total number of seeds is small, so are the number of ejected holes. We only give an upper limit here, because the ejected holes are too rare for us to give a robust value given the statistical limitations of our "halo cloning" method for populating the entire halo population at z = 6 (the mass density in ejected BHs can be large in models with large f_seed; see below).

These models represent the simplest scenarios for SMBH formation, requiring a very brief period of seed formation and a prolonged period of accretion at rates comparable to the Eddington rate, and consequentially they represent the most pessimistic predictions for LISA's observational prospects.

While the suppression of BH seed formation is an attractive possibility that fits constraints on the z = 6 SMBH mass function, it is clearly not unique. One alternative solution to the overproduction problem is to simply reduce the accretion rate of lower mass BHs at lower redshifts. This is again physically plausible: accretion could be choked off as a result of the baryonic gas being churned into stars, being heated and dispersed by re-ionization feedback, or through self-induced negative feedback where the BH's own accretion-powered radiation stops the gas supply. Rather than try to model such a time-dependent (and probably mass-dependent) mass accretion scenario, we examined several model variants, in which BHs are allowed to accrete just enough mass to match the value inferred by the m–σ relation between BH mass and host halo velocity dispersion. That is, at each time step t → t + Δt, all BHs were assumed to grow in mass by m → m + Δm such that the m–M relation is satisfied at the new host halo mass and redshift. However, whenever this requires super-Eddington growth, i.e., if (m + Δm)/m > exp(Δt/t_Edd), then Eddington growth is applied instead. The main additional assumption here is that the m–σ relation remains valid at all redshifts (which is at least consistent with a comparison between the evolution of quasars and early-type galaxies at 0 < z ≲ 6; Haiman et al. 2007; Shankar et al. 2009b). As above, we adopt Equation (20) as our extrapolated m/M relation.

The BHs in these models form as 100 M_☉ seeds, and, given their host halo mass and redshift, accrete to match this relation as closely as possible without exceeding the Eddington accretion rate. If the simulation completes with the mean BH accretion rate well below the Eddington rate at all redshifts, then it is consistent with satisfying the Eddington accretion rate and the m–M relation inferred by Equation (20). As we shall find below, in our models the maintenance of the m–σ relation does not typically require that the Eddington accretion rate is saturated (see Figure 15 below) as long as the f_seed ≳ 10⁻². If both the seed mass and the seeding fraction are low, it is increasingly difficult to satisfy Equation (20) at higher redshifts while simultaneously satisfying the Eddington upper limit. We find that for f_seed ≲ 10⁻², accretion must saturate at the Eddington rate for much of z ≳ 15 until the extrapolated m–σ relation is satisfied, with the mass function falling below this relation at earlier stages of growth.

Recoil velocities are calculated with the spin magnitudes chosen uniformly between 0.0 and 0.9. As with our previous models, we run simulations where the spins are either randomly oriented or completely aligned with the angular momentum vector of the binary orbit. This class of models in effect represents the most optimistic LISA expectations, as it allows us to keep numerous seeds, while simply adjusting the accretion rate, as described above, to keep the mass function within bounds. Note that the recoil speeds are also minimized by our choice for the spin alignment.

We show the mass functions and occupation fractions for this model in Figure 13, for three different values of the cutoff redshift below which new seeds are not formed, z_cut = 0, 12 and 18. Note that it is still possible to form the SDSS quasar-SMBHs via Eddington-limited accretion by z ≈ 6 even if seeds form in just 0.1% of all 1200 K halos and only before z = 18. We do not plot the m/M relation, as it is satisfied in the form of Equation (20) in all cases shown here, by construction. These models also satisfy, by construction, the upper limit for the SMBH mass density. In all of the models shown in Figure 13, ρ_SMBH,5+(z = 6) ≳ 1.3 × 10⁴ M_☉ Mpc⁻³ if f_seed ⩾ 10⁻². In the f_seed = 10⁻³ models, we find ρ_SMBH,5+(z = 6) ∼ 1.0 × 10⁴ M_☉ Mpc⁻³. The difference in ρ_SMBH is due to varying occupation fractions at the low end of the halo mass function. The mass functions in Figure 13 have a much shallower slope overall than those in Figures 5 and 8. For the mass functions shown earlier, the steeper slopes were due to the ratio m/M increasing with decreasing host halo mass; the observed m/M relation has the opposite trend.

Zoom In Zoom Out Reset image size

We show the LISA event rates in these alternative models in Figure 14. Most significantly, we note that in the f_seed = 1 versions of these successful models, the LISA event rate can be as high as 30 yr⁻¹. (Note that this number can be even higher if seeds can form in minihalos down to a virial temperature that is significantly lower than our assumed fiducial value of 1200 K.) The rate is somewhat suppressed when compared to the earlier Pop-III seed models (Figure 7) that exceeded realistic indications on the SMBH mass density, because the massive BHs in the LISA band are rarer due to the more modest growth rates. We draw the attention of the reader to the apparent independence of the detection rate on the seed fraction and seed formation cutoff in the cases where f_seed ≳ 10⁻¹ and z_cut ≲ 12. Because BH ejections probabilities are lower in these models when compared to the constant-accretion scenarios of Figures 7 and 10, and because the m/M ratios are the same function of M in all models shown, the LISA rates saturate and converge once the occupation fraction in the most massive halos approach unity. Because the T = 1200 K halos form in greatest abundance at z ∼ 20, the z_cut = 12 case hardly differs from the case with no seed cutoff.

Zoom In Zoom Out Reset image size

A key characteristic of any SMBH assembly scenario is the balance between growth through BH mergers and growth through gas accretion. As discussed above, the two must strike a balance such that they are able to account for the most massive observed quasar-SMBHs at z ∼ 6, while also not exceeding the total observed SMBH mass density. In Table 3, we illustrate the relative importance of mergers versus growth in the models presented in this paper: the four successful constant-accretion models from Table 2; two of the unsuccessful constant-accretion models, also in Table 2, which overproduce the universal SMBH mass density; and four of the models that explicitly follow the extrapolated m–σ relation via Equation (20). The values shown in the table (in log₁₀) are the sum of the initial masses of all the seed BHs that enter our merger trees; the total mass of galactic BHs at z = 6;⁵ and the total mass of BHs ejected before z = 6. We also calculate the ratio of the total (galactic and ejected) BH mass at z = 6 to the total initial seed BH mass, which gives a simple measure of the growth through gas accretion. We see immediately the contrast between the two types of successful models: the constant-accretion scenario relies on gas accretion for much of the growth, typically several orders of magnitude in the total BH mass; where the self-regulating models essentially describe the most heavily merger-driven scenarios possible, requiring accretion-driven growth of as little as a factor of a few.

These models also produce a significant population of ejected BHs. Even though ejection rates are lower on the whole than our constant-accretion models (compared to the unsuccessful constant accretion Model X in Tables 2 and 3), two factors contribute to the ejected BH mass being comparable to the galactic BH mass at z = 6. First, seed BHs are allowed to be very common, especially in contrast to the successful constant-accretion Models A through D; this results in a far greater number of total merger events, and a high total number of ejections despite the lower ejection probabilities. Second, the surviving BHs do not grow nearly as rapidly in these models as in the constant-accretion scenarios, so that the ratio between the total galactic BH mass and the total ejected mass can remain low (whereas in the constant-accretion scenarios, retained holes can rapidly outgrow their escaped counterparts). In the f_seed = 1 models, the ejected holes can outnumber and outweigh their retained galactic counterparts with mass densities of ∼3 × 10⁴ M_☉ Mpc⁻³ and number densities of ∼100 Mpc⁻³. The mass function of the ejected holes is peaked slightly above the seed mass (because holes are most likely to be ejected at the earliest stages of their evolution, when their host halos are the least massive). All of the ejected BHs are ≲10⁴ M_☉ if spins are aligned, but in rare instances, SMBHs as massive as ∼10⁸ M_☉ are ejected in our models with randomly oriented spins (the ejected SMBHs with masses above m > 10⁶ M_☉ have a very low number density, SMBHs of ∼4 × 10⁻⁵ Mpc⁻³, even in the model with the most ejections (random orientation, no cutoff redshift, f_seed = 1).

These self-regulating accretion models work by adjusting the BH accretion rates according to the mass growth of their host halos. We plot the accretion rates in units of the Eddington rate for this new set of models in Figure 15. We do so for the z_cut = 12 case with four combinations for f_seed and spin alignment, and for three different BH mass ranges: 10³ M_☉ ⩽ m ⩽ 10⁶ M_☉, 10⁶ M_☉ ⩽ m ⩽ 10⁸ M_☉, and m ⩾ 10⁸ M_☉. Note that the accretion rates must be slightly higher if BH binary spins are randomly oriented, in order to compensate for the higher ejection rates. Similarly, accretion rates are higher if seeds are less common, in order to compensate for the reduced merger-driven growth. For the models shown, the duty cycle (the mean accretion rate in Eddington units) for the most massive SMBHs converge to ∼0.2 at z ≈ 6, though it can be as high as ≳0.5 at z ≳ 8 if merger-driven growth is hindered by low occupation fraction or recoil-induced ejections.

Zoom In Zoom Out Reset image size

We note that similar merger-tree models tracking SMBH growth have been published. For example, Koushiappas et al. (2004) presented a similar model where the SMBHs are assembled primarily through mergers of directly collapsed halo cores. Bromley et al. (2004) also presented SMBH assembly model wherein gas accretion activity was triggered by major mergers of the BHs' host halos; in their model, a set fraction of the baryonic mass of the host was fed to the BH at each major merger. Their prescription (albeit without gravitational recoil) successfully produced the most massive SDSS SMBHs before redshift z = 6 without overproducing the mass density. In general, this type of assembly model is fairly easily tuned to broadly reproduce the m–σ relation, as the parallel mass growth of BHs and their host halos is built in.