LINKING THE SPIN EVOLUTION OF MASSIVE BLACK HOLES TO GALAXY KINEMATICS

Observational evidence for the existence of massive black holes (MBHs) in the center of galaxies is circumstantial but convincing. Stellar orbits near the Galactic center point at the presence of a dark and massive object, which is too heavy and compact to be a cluster of low-luminosity bodies (Maoz 1998) or a fermion star (Schödel et al. 2002). Also, near-infrared observations of this object show that it is likely to have an event horizon because if it did not, it would emit thermal radiation and be more luminous than observed (Broderick et al. 2009; Broderick & Narayan 2006). This disfavors exotic horizonless alternatives to general relativity's black holes such as, e.g., boson stars and gravastars. Similarly, MBHs are expected to be present in the center of most galaxies (and not just ours), as their presence is required to explain quasars (QSOs), active galactic nuclei (AGNs; Soltan 1982), and the cosmic downsizing (Scannapieco et al. 2005; Croton et al. 2006; Bower et al. 2006). Indeed, observations of nuclear stellar and gas dynamics, reverberation mapping, spectroscopic single-epoch measurements, and spectral energy distribution (SED) fitting (see, e.g., Ghez et al. 2008; Gillessen et al. 2009; Kormendy & Richstone 1995; Miyoshi et al. 1995; Hicks & Malkan 2008; Blandford & McKee 1982; Kaspi et al. 2005; Calderone et al. 2013; Castignani et al. 2013 and references therein for a discussion of the different methods) corroborate this expectation and allow the mass of MBHs to be estimated. Furthermore, these MBH mass measurements have in turn resulted in the discovery of the MBH–galaxy relations (e.g., Gebhardt et al. 2000; Ferrarese & Merritt 2000; Marconi & Hunt 2003; Häring & Rix 2004; Gültekin et al. 2009), hinting at a symbiotic evolution of the central compact object and its host galaxy.

MBHs are also expected to possess a spin-angular momentum. The combination of mass and spin completely characterizes these objects if they are described by the black hole solutions of general relativity (Kerr 1963).⁷ If instead the gravity theory describing our universe is not exactly given by general relativity, other charges may be needed to describe MBHs, but the spin will still be one of them (see Barausse & Sotiriou 2013; Yagi et al. 2012; Pani et al. 2011; Pani & Cardoso 2009, for a few examples of black holes in alternative gravity theories). Indeed, future space-based gravitational wave interferometers such as ESA's L3 experiment eLISA (The eLISA Consortium et al. 2013) will measure MBH masses and spins with fantastic accuracy (∼0.1% and ∼1%, respectively) and also put general relativity to the test by, e.g., measuring possible additional black hole charges. More importantly from an astrophysical point of view, measurements of MBH spins could unveil still unknown links between the MBHs and their galactic hosts. As an example, knowledge of the spins could constrain the properties of the gas accreting onto the MBHs, e.g., whether it spirals toward it on a preferential plane, on completely isotropic directions, or with other more complex dynamics (Dotti et al. 2013, hereinafter D13).

Until eLISA or a similar space-based gravitational wave interferometer is launched, however, measuring the spins of MBHs is considerably more difficult than measuring their masses. This is because the spin affects the dynamics of stars and gas surrounding the MBH only at distances that are orders of magnitude smaller than the MBH influence radius (e.g., Bardeen & Petterson 1975). Currently, the most accurate way of measuring the spins of MBHs is through the spectra of relativistically broadened Kα iron lines (see Reynolds 2014; Brenneman 2013 for recent reviews on the topic). Indeed, using current X-ray spectrographs (in particular, XMM and Suzaku), a significant sample of MBHs with spin estimates is being built up. Also, more stringent constraints on the spin estimates are becoming available from hard X-ray data from NuSTAR (Risaliti et al. 2013; Marinucci et al. 2014a, 2014b). Finally, spin measurements with Kα iron lines are one of the main scientific drivers of ESA's upcoming L2 mission, ATHENA+ (Nandra et al. 2013).

This situation calls for a theoretical framework capable of interpreting available data for MBH spins and making testable predictions for future measurements. However, detailed modeling of the MBH spin evolution is still missing. Since the spin is a vector, any theoretical model should take the evolution of its direction into account, because that also affects the evolution of the magnitude. This is particularly clear when one considers MBH mergers, which are expected in the late stages of galaxy mergers (e.g., Begelman et al. 1980). The final spin, mass, and kick velocity of the resulting MBH remnant critically depend on the orientations of the spins of the progenitors (Rezzolla et al. 2008; Barausse & Rezzolla 2009; Buonanno et al. 2008; Barausse et al. 2012; Lousto & Zlochower 2009; van Meter et al. 2010).

The main driver of the MBH evolution, however, is not expected to be given by mergers but rather by radiative efficient gas accretion (Soltan 1982; Merloni et al. 2004; Shankar et al. 2004, 2013; Wang et al. 2009; Li et al. 2012). In this case, the spin magnitude evolution depends on whether the gas in the very central regions of the accretion disk is corotating or counterrotating relative to the MBH spin. This, in turn, depends on the larger-scale geometry of the fueling process (i.e., if the gas keeps falling toward the MBH on some preferential plane), as well as on the possible realignment of the MBH–disk system due to the Bardeen–Petterson effect (Bardeen & Petterson 1975). The latter is the interplay between the frame dragging from the rotating spacetime geometry and the viscous stresses in the disk. The details of the effect of the spacetime's geometry on the properties of accretion disks are still not completely understood (see, e.g., Sorathia et al. 2013a, 2013b and references therein). All the studies that follow the coupled evolution of the magnitude and direction of the MBH spin have been performed under some simplifying hypotheses: i.e., assuming that the disk has isotropic viscosity (see the discussion in Sorathia et al. 2013a, 2013b), using order of magnitude estimates for the effect of the alignment (e.g., King & Pringle 2006), or under the small-deformation approximation for the warped disk (e.g., Perego et al. 2009).

These simplifying hypotheses allowed for simple recipes for the spin evolution that have been used in semianalytical models of structure formation and evolution (see, e.g., Volonteri et al. 2005, 2007, 2013; Berti & Volonteri 2008; Lagos et al. 2009; Fanidakis et al. 2011, 2012; Barausse 2012). However, most of these studies (Volonteri et al. 2005; Berti & Volonteri 2008; Lagos et al. 2009; Fanidakis et al. 2011, 2012) assume that the MBH fueling is either always coherent (i.e., accretion always occurs exactly on the equatorial plane) or statistically isotropic (i.e., the MBH accretes small gas clouds with randomly oriented angular momenta). These assumptions were partly relaxed in Barausse (2012) and Volonteri et al. (2007, 2013), who noted that the Bardeen–Petterson effect is likely to align the MBH spin with the angular momentum of the accreting gas in gas-rich nuclear environments, thus making the accretion flow effectively coherent in such situations, as opposed to isotropic accretion, which is more likely in gas-poor nuclear environments. Even with this improvement, however, the possible parameter space of the accretion flow is still limited to two points (coherent accretion in gas-rich environments versus isotropic accretion in gas-poor environments).

Recently, D13 demonstrated that the spin evolution can be significantly different if a more realistic description of the accretion flow is adopted, namely, one that allows regimes intermediate between a perfectly coherent and a perfectly isotropic flow. Dubois et al. (2014a) studied the evolution of the spins of a population of MBHs, taking into consideration the dynamical properties of the fueling gas. They extracted the direction of the angular momentum of the accreting matter on the smallest scales resolved by hydrodynamical simulations (down to ≈10 pc in their highest-resolution run). In a similar spirit, in this paper we couple the semianalytic model of Barausse (2012, hereinafter B12) to the spin evolution model of D13. In particular, we link the degree of coherency/isotropy of the accretion flow to the dynamical parameters of the gas and stars observed in galactic nuclei. We stress that also in our case, most of the observations that we use have a spatial resolution ≳ 100 pc, significantly larger than typical accretion disk scales, and that the dynamical properties of the gas could change at smaller unresolved scales (e.g., Hopkins et al. 2012; Maio et al. 2013; Dubois et al. 2014b). In this sense, therefore, a comparison of the predictions of our model for the MBH spin evolution with actual spin measurements provides a way of testing the dynamics and coherency of the nuclear gas at yet unresolved scales. We attempt such a test with the MBH spin measurements from Kα iron lines available to date. We note, however, that our model provides a flexible framework to also interpret future, more accurate measurements from experiments such as eLISA and ATHENA+.

The paper is organized as follows. In Section 2 we present the details of our model for the evolution of MBHs in their galactic hosts. Improving on B12, we adopt a different star formation law (see Section 2.1 and Appendix A) and, more importantly, we update the prescription for the evolution of the MBH spins under accretion, following D13 (see Section 2.2). The degree of coherence of the accretion flow is then linked to the kinematic properties of the host galaxy in Section 3. We consider, in particular, three models: one where the degree of coherence is related to the gas kinematics, one that links it to the stellar kinematics, and a hybrid model where information from both the gas and stellar kinematics is used. In Section 4, we calibrate the free parameters of our semianalytical galaxy formation model against observations. We present our predictions for the evolution of MBH spins in Section 5, and in Section 6 we compare them to measurements from Kα iron line profiles. In Section 7 we draw our conclusions.

Our model is based on B12, who studied the cosmic evolution of MBH masses and spins using a state-of-the-art semianalytical galaxy formation model. More specifically, B12 describes the dark-matter evolution by merger trees produced with the extended Press–Schechter formalism and suitably modified to reproduce the results of N-body dark-matter simulations (Cole et al. 2008). The baryonic content of galaxies is evolved along the branches of the merger tree and includes several components: a hot-gas phase, a cold-gas disk, a stellar disk, a cold-gas bulge, a stellar bulge, a MBH, and a reservoir of cold gas fueling accretion onto the MBH. Gravitational interactions among the various components are accounted for (although in simplified ways), allowing the computation of the density and velocity profiles of the various components inside each galaxy. Also modeled are a variety of nongravitational interactions, which are summarized schematically in Figure 1.

Zoom In Zoom Out Reset image size

For more details, we refer the reader to B12. Here, we focus instead on our improvements to that model: the star formation law (Section 2.1) and the description of the cosmic MBH spin evolution (Section 2.2).

2.1. Star Formation Model

The star formation model adopted here differs from and improves on the one used in B12. In this section, we describe the general framework and main equations and confine a more detailed presentation to Appendix A.

The first new ingredient is that the star formation is now an explicit function of metallicity, which allows for better modeling of the buildup of stars in different galaxies at various redshifts and environments (see the star formation history of Figure 9, to be compared with Figure 8 of B12). The second is a different prescription for classical bulges and pseudobulges. The former are the result of major galactic mergers, which cause bursts of star formation (Daddi et al. 2010; Genzel et al. 2010), while the latter result from bar instabilities in galactic disks and undergo a slower, disk-like star formation (see, e.g., Kormendy 2013). In B12, both bar instabilities and major mergers were modeled, but star formation in the bulge component was the same irrespective of the bulge's formation mechanism. In practice, our new prescription is based on the work by Krumholz et al. (2009), extended to include the low metallicity (less than 1% solar; Krumholz & Gnedin 2011; Krumholz 2012; Forbes et al. 2014; Kuhlen et al. 2013) and starburst regimes (Daddi et al. 2010; Genzel et al. 2010).

In galactic disks, the star formation density rate is expressed in terms of the fraction f_c of the background gas surface density (Σ_g) that is converted into stars on a timescale t_SF,

$$\begin{equation} \dot{\Sigma }_* = \frac{f_{\rm c} \Sigma _{\rm g}}{t_{\rm SF}}\,. \end{equation} \tag{ 1 }$$

The metallicity dependence enters explicitly in the fraction f_c (see Equation (A4) in Appendix A): down to 2% solar, the lower the metallicity, the smaller the amount of gas available for star-forming. The timescale t_SF (see Equation (A2) in Appendix A) depends on the properties of giant molecular clouds (GMCs), which vary with the gas richness of the environment. This timescale is generally longer than the local free fall timescale. The same star formation law applies to pseudobulges, but it is calculated per unit volume,

$$\begin{equation} \dot{\rho }_* = \frac{f_{\rm c} \rho _{\rm g}}{t_{\rm SF}}. \end{equation} \tag{ 2 }$$

Finally, in classical bulges, we describe the burst of star formation as the sudden consumption of gas in a local dynamical timescale,

$$\begin{equation} \dot{\rho }_{*}= \frac{\rho _{\rm g}}{t_{\rm ff}}, \end{equation} \tag{ 3 }$$

where $t_{\rm ff}=\sqrt{3\pi /(32 G \rho _{\rm g})}$.

The total star formation rates are then performed by integrating Equations (1)–(3) over the whole disk/pseudobulge/bulge.

2.2. Spin Evolution

Two main processes drive the cosmic evolution of MBH spins: gas accretion from the interstellar medium (ISM) or the intergalactic medium (IGM) and mergers between MBHs.

2.2.1. Gas Accretion

Accretion is believed to be the main driver of the MBH spin evolution (Berti & Volonteri 2008; B12), with the possible exception of the most MBHs in massive, gas-poor, low-redshift ellipticals (e.g., Fanidakis et al. 2011). Because of the spin's vectorial nature, its evolution depends on the angular momentum magnitude and direction of the accreting gas. In our implementation, based on D13, the circumnuclear gas reservoir surrounding the MBH is accreted in a sequence of clouds, with a nontrivial angular momentum distribution. Unlike in D13, however, this distribution is not arbitrary but is linked to the galactic properties at larger scales and therefore intimately connected to the galaxy morphology (see Section 3).

Following B12, we distinguish two main astrophysical modes of gas accretion, which prompt different spin evolution channels as illustrated in Figure 2.

• 1.  
  A hot "radio" mode. In this regime, the MBH accretes IGM gas that cools on a timescale longer than its free fall time. The gas is approximately spherically symmetric, so no net angular momentum is acquired by the MBH, and the spin parameter a_bh changes (to first approximation) only because of the increase of the MBH mass:

  $$\begin{equation} \dot{a}_{\rm bh,radio} = -2 a_{\rm bh} \frac{ \dot{M}_{\rm bh,radio}}{{M}_{\rm bh}}, \end{equation} \tag{ 4 }$$

  where $\dot{M}_{\rm bh,radio}$ is the Bondi accretion rate of the hot IGM component (see Equation (44) in B12). This mode has a negligible direct impact on the mass and spin evolution of MBHs because it is typically characterized by very low accretion rates but plays a major role in activating jets capable of exerting an "AGN feedback" on the growth of large cosmic structures (see B12 for more details on the implementation of AGN feedback).
• 2.  
  A cold "QSO" mode. It is associated with the star formation in the bulge, which drives cold ISM gas into a low-angular-momentum reservoir of mass M_res, available for accretion onto the MBH. In practice, we assume that the influx of gas into the reservoir is proportional to the bulge star formation rate (see B12; Granato et al. 2004; Haiman et al. 2004; Lapi et al. 2014). As a result, in our current framework for the star formation (Section 2.1), major mergers tend to produce larger reservoirs and MBH accretion rates than disk instabilities. Once the cold gas has settled in the reservoir, we follow Granato et al. (2004) and assume that it accretes onto the MBH with instantaneous rate,

  $$\begin{equation} \dot{M}_{\rm QSO}=\min (\dot{M}_{\rm visc}, A_{\rm Edd}\dot{M}_{\rm Edd}). \end{equation} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

The viscously driven accretion rate reads $\dot{M}_{\rm visc}= k_{\rm accr}\sigma ^3 M_{\rm res}/(G M_{\rm bh})$, where⁸ k_accr = 10⁻³, σ = 0.65V_vir and V_vir is the halo's virial velocity (Granato et al. 2004; G. L. Granato, 2009 private communication). The Eddington rate $\dot{M}_{\rm Edd}=L_{\rm Edd}/[\eta (a_{\rm bh}) c^2]$ is derived from the Eddington luminosity L_Edd, assuming a radiation efficiency η. As shown below (see Equation (12)), η can be computed as a function of the MBH spin and the fraction of prograde vs. retrograde accretion events. The free parameter A_Edd ⩾ 1 sets the amount of super-Eddington accretion allowed in our model. This prescription for $\dot{M}_{\rm QSO}$ replaces the simpler recipe of B12, which assumed $\dot{M}_{\rm QSO}=M_{\rm res}/t_{\rm accr}$, with t_accr a free parameter on the order of 5 × 10⁸ yr.

The main difference with respect to B12 consists in how the gas in the reservoir is accreted. In B12, depending on the ratio M_res/M_bh, M_res is either accreted coherently or as a collection of small randomly oriented clouds with vanishing total angular momentum. Here, we employ a model that allows for a broader range of configurations for the accreting flow, as we proceed to explain. This model is based on results by Perego et al. (2009) and D13, to which we refer for more details.

Accretion is described as a series of transient accretion disks, each one formed following the tidal disruption of a cloud. The mass M_disk of each individual disk is the minimum between the typical mass of a molecular cloud—M_cloud = 3 × 10⁴ M_☉ is our fiducial value—and the mass M_sg of a disk truncated at the self-gravity radius, where the Toomre parameter Q reaches a value of 1. This is because outside this radius, the accretion disk that forms from the infalling cloud will fragment under its self-gravity and be consumed by star formation. The explicit expression for M_sg is given by (D13)

$$\begin{equation} M_{\rm sg} \approx 2 \times 10^{4} {\alpha }_{0.1}^{-1/45} \left(\frac{f_{\rm Edd}}{\eta _{0.1}}\right)^{4/45} \left(\frac{M_{\rm bh}}{10^6 \, M_\odot }\right)^{34/45} \, M_\odot, \end{equation} \tag{ 6 }$$

where f_Edd is the Eddington ratio, η_0.1 is the disk radiative efficiency normalized to 0.1, and α_0.1 is the disk viscosity parameter in units of 0.1. Throughout this paper we assume α_0.1 = 1. Each accreting lump gets tidally disrupted and settles in a standard thin accretion disk with an associated angular momentum J_disk obtained from its mass, M_disk = min (M_cloud, M_sg):

$$\begin{eqnarray} J_{\rm disk}=\;&\frac{8}{21}\frac{\dot{M}_{\rm bh,QSO} \sqrt{G M_{\rm bh}}}{A_{\nu }}R_{\rm disk}^{7/4}\, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} A_{\nu }=\; & 9.14 \times 10^{6}\alpha _{0.1}^{4/5}\left(\frac{M_{\rm bh}}{10^6 \, M_\odot }\right)^{1/20}\nonumber \\ &\times \left(\frac{f_{\rm Edd}}{\eta _{0.1}} \right)^{3/10}{\rm cm^{5/4}\ s^{-1}}\, \end{eqnarray} \tag{ 8 }$$

(see D13). The extension of the disk R_disk depends on its mass. If M_sg < M_cloud, it is given by (D13)

$$\begin{equation} {R_{\rm disk}}\approx 10^5 R_{\rm g}\, \alpha _{0.1}^{28/45}\left(\frac{M_{\rm bh}}{10^6 \, M_\odot }\right)^{-52/45}\left({f_{\rm Edd}\over \eta _{0.1}}\right)^{-22/45}, \end{equation} \tag{ 9 }$$

where R_g = 2GM_bh/c² is the Schwarzschild radius. If instead M_sg > M_cloud, one has (D13)

$$\begin{eqnarray} {R_{\rm disk}}\approx \; & 4\times 10^4\, R_{\rm g} \left({M_{\rm disk}\over 10^4\,M_\odot }\right)^{4/5}\alpha _{0.1}^{16/25}\nonumber \\ &\times \left(\frac{M_{\rm bh}}{10^6 \, M_\odot }\right)^{-44/25} \left({f_{\rm Edd}\over \eta _{0.1}}\right)^{-14/25}. \end{eqnarray} \tag{ 10 }$$

According to the initial direction of the lump's angular momentum, the Bardeen–Petterson effect will either align or antialign the disk and the MBH spin directions (see Equation (13)), resulting in prograde or retrograde accretion, respectively. As the total mass in the reservoir is consumed by several lumps accreting onto the MBH, the time-averaged spin evolution is determined by the fraction w (1 − w) of lumps that accrete on prograde (retrograde) orbits:

$$\begin{eqnarray} \dot{a}_{\rm bh,QSO} &=& \left \lbrace \left [w \, {L}^{\rm pro}_{\rm _{\rm ISCO}}(a_{\rm bh})+(1-w)\,{L}^{\rm retro}_{\rm _{\rm ISCO}}(a_{\rm bh})\right ]\right. \nonumber\\ &&\left. -2 a_{\rm bh} \left [w \,{E}^{\rm pro}_{\rm _{\rm ISCO}}(a_{\rm bh}) +(1-w)\,{E}^{\rm retro}_{\rm _{\rm ISCO}}(a_{\rm bh})\right ]\right \rbrace \frac{ \dot{M}_{\rm QSO}}{{M}_{\rm bh}}, \nonumber\\ \end{eqnarray} \tag{ 11 }$$

where $L_{\rm _{\rm ISCO}}^{\rm pro}(a_{\rm bh})$ and $E^{\rm pro}_{\rm _{\rm ISCO}}(a_{\rm bh})$ are the specific angular momentum and specific energy at the prograde innermost stable circular orbit (ISCO), respectively, while $L_{\rm _{\rm ISCO}}^{\rm retro}(a_{\rm bh})$ and $E^{\rm retro}_{\rm _{\rm ISCO}}(a_{\rm bh})$ are the same quantities for the retrograde ISCO. Likewise, the accretion efficiency is calculated by weighing the prograde and retrograde thin-disk efficiencies η_pro(a_bh) and η_retro(a_bh):

$$\begin{equation} \eta (a_{\rm bh}) = w \eta _{\rm pro}(a_{\rm bh})+ (1-w) \eta _{\rm retro}(a_{\rm bh}). \end{equation} \tag{ 12 }$$

The fraction w of prograde accretion disks is determined by the Bardeen–Petterson effect (Bardeen & Petterson 1975; King et al. 2005; King & Pringle 2006; Perego et al. 2009; D13) and is given by

$$\begin{equation} \begin{array}{lcc}w=1 & {\rm if} & J_{\rm disk}> 2J_{\rm bh},\\ \displaystyle w=F + \frac{J_{\rm disk}}{2 J_{\rm bh}}(1-F)& {\rm if} & J_{\rm disk} < 2J_{\rm bh}. \end{array} \end{equation} \tag{ 13 }$$

Here, the "isotropy parameter" F measures the fraction of accretion events with angular momentum initially pointing in the hemisphere of the MBH spin (i.e., with angular momentum initially tilted by an angle θ_out < π/2 relative to the MBH spin).⁹ When J_disk > 2J_bh, the MBH spin aligns with the angular momentum of the outer regions of the accretion disk that forms when the lump of matter falls into the MBH, and the alignment happens on a timescale smaller than the accretion timescale (Perego et al. 2009; D13). Most of the accreted gas then has angular momentum aligned with the MBH spin, hence w = 1, and the MBH spins up. In this limiting case, Equation (11) tends to Equation (38) of B12. If instead J_disk < 2J_bh, King et al. (2005) and King & Pringle (2006) showed that on a similarly short timescale, the inner accretion flow's angular momentum antialigns with the MBH spin if the angle between the outer disk angular momentum and the MBH spin satisfies θ_out > π/2 and cos θ_out < −J_disk/(2J_bh). Assuming the fraction 1 − F of accretion events with θ_out > π/2 is distributed isotropically, one then obtains the second expression for w in Equation (13).

Note that in the limit of J_disk/2J_bh ≪ 1 (always valid for very massive MBHs, M_bh ≳ 10⁹ M_☉), our model becomes very simple, as it predicts that the MBHs will move toward an equilibrium spin parameter a_eq function of F, obtained by setting $\dot{a}_{\rm bh,QSO}=0$ in Equation (11) after substituting w with F as given by Equation (13). The value of a_eq as a function of F is shown in Figure 3.

Zoom In Zoom Out Reset image size

At this point, however, we should stress that Equations (11) and (13) (1) only provide an "averaged" version of the stochastic scenario proposed by D13: i.e., our prescription does not follow each individual accreting cloud but only approximates the model of D13 over many such clouds and (2) do not take into account the detailed evolution of the accretion disk during the (short) alignment process. respectively¹⁰ A more accurate treatment of the spin evolution addressing these two shortfalls is beyond the scope of the paper, but we expect the impact of these simplifications on the spin evolution to be small. Regarding (1), the mass M_res of the reservoir available for accretion is typically much larger than M_disk; thus accretion takes place in several "lumps," and the stochastic character of the model of D13 averages out. An exception is in the "tails" of the accretion process, when M_res is small before and after a QSO event. In those tails, however, accretion rates and spin changes are small due to the scarcity of gas. We have indeed tested that replacing M_cloud → min (M_res, M_cloud), which corresponds to changing our model prescription in the tails of the accretion process, does not alter our results significantly. As for (2), neglecting the short alignment process is typically justified, as the alignment takes place on timescales shorter than accretion. Nevertheless, even a small amount of retrograde accretion during the alignment can potentially lower the spin significantly near the a_bh = 1 limit; more details on this point are given in Appendix C.

2.2.2. MBH Mergers

Our prescription for MBH mergers closely follows that of B12, to which we refer for more details. In summary, we use the phenomenological formulae of Barausse & Rezzolla (2009) and van Meter et al. (2010) to predict the final spin and recoil velocity of the MBH resulting from the merger and check whether the recoil ejects the MBH from the galaxy. As for the mass of the MBH produced by the merger, we update the prescription of B12 by adopting the phenomenological formula of Barausse et al. (2012). We stress that these formulae produce results in accurate agreement with fully general relativistic simulations. The only other change from the implementation of B12 regards the criterion adopted to discriminate "wet" (i.e., gas-rich) mergers from "dry" (i.e., gas-poor) mergers. This distinction is important because in wet mergers the spins of the two MBHs are expected to align to each other—we follow here Dotti et al. (2010) and assume that the MBH spins are almost aligned, to within ∼10°, by the circumbinary disk via the Bardeen–Petterson effect—whereas in dry mergers such alignment does not take place, and the spins are isotropically distributed. Improving on B12, who compared M_res to the mass of the MBH binary to distinguish the two regimes, here we assume that a merger happens with almost aligned spins if J_res > 2(S₁ + S₂) (where J_res is the angular momentum of the cold-gas reservoir, and S_{1, 2} are the spins of the two MBHs), while otherwise we assume that the spins are randomly oriented. This prescription is preferable to that of B12 because the Bardeen–Petterson effect is sensitive to the angular momenta of the circumbinary disk and the MBHs, rather than to their masses.

In the previous section, we described how to physically relate the degree of anisotropy of the accretion flow to the MBH spin changes. The details of this relation may still need to be understood, but there is a general consensus in the community on the overall physical picture (Bardeen & Petterson 1975; King et al. 2005; King & Pringle 2006; Perego et al. 2009; D13). The outstanding question instead is what are the actual conditions of the flow at subparsec scales and what determines these conditions. This problem is considerably harder to solve because (1) there is hardly any direct observational evidence at those scales and (2) those scales are extremely challenging to explore theoretically with galaxy formation simulations due to lack of resolution and the unknown role of dissipative/nonlinear processes.

In an attempt to contribute some insight into the problem, we link the dynamical properties of the host galaxies, which can be easily observed, to the properties of the environment in the immediate vicinity of the MBH. As we will argue later in this paper, existing MBH spin measurements allow for testing and placing constraints on this conjectured link. The basic assumption behind our attempt is that the mechanism fueling the MBH—either disk instabilities or galaxy mergers—bridges large and small scales, ensuring that the degree of "disorder" of each galactic component (gas and stars) does not change much from galactic to circumnuclear scales. At kiloparsec scales, a routinely observed measure of "disorder" is the v/σ ratio, where v is the bulk rotation velocity of the system, and σ is its velocity dispersion. Therefore, the remaining issue is whether the gas that is brought down to the MBH and gets accreted has a v/σ more similar to that of the surrounding stars in the bulge or to that of the gas in the galactic disk. In the absence of solid theoretical predictions and observational evidence, we explore a few models, further differentiating between elliptical and spiral galaxies. As a consequence, we are able to predict different spin evolution in different galaxy types, which allows us to constrain our models (and therefore the physics of the fueling mechanism) by comparing them to the observed MBH spin distribution in the local universe.

In the following, we first discuss observations of the v/σ ratio in galaxies (Section 3.1). These data will guide us in constructing models to assign a v/σ ratio to the accreting reservoir, according to the morphology of the host galaxy (Section 3.2). Finally, in Section 3.3 we will describe how to geometrically translate a given a value of the v/σ ratio for the reservoir into the isotropy parameter F used in the previous section to characterize the MBH spin evolution; see Equations (13) and (11).

3.1. Measurements of v/σ

The ratio v/σ for the stellar and/or gaseous components of a galaxy has been measured for a variety of galaxy samples. We collected these data to construct v/σ distributions to implement in our galaxy evolution model. In the following we describe these measurements in detail, discriminating between gas measurements in galactic disks and stellar measurements in bulges.

3.1.1. Gaseous Disks

Measurements of v/σ for the gas component in spiral galaxies have been obtained with integral field emission line spectroscopy with subarcsec resolution, tracing the gas from 10 kpc down to subkiloparsec scales, depending on the redshift of the system. More specifically, Kassin et al. (2012) explicitly discuss the dependence of v/σ on redshift and stellar mass for disk galaxies at z < 1.2 using a sample of several hundred objects. Here, we complement their data with measurement at higher redshift. Wisnioski et al. (2011), Epinat et al. (2012), and Swinbank et al. (2012) report v/σ measurement for star-forming galaxies in the redshift range 1.2 < z < 1.7; their data are combined in a catalog of 30 objects that we take as a representative sample at z ≈ 1.5. Förster Schreiber et al. (2009), Law et al. (2009), and Newman et al. (2013) explore the range 1.5 < z < 2.5; we combine their data in a catalog of 56 objects that we take as a representative sample at z ≈ 2. Gnerucci et al. (2011) measure v/σ of 33 star-forming galaxies at high redshift, in the range 3 < z < 4.8. They find a fraction of rotating systems (defined as systems having v/σ > 1) of ∼1/3, claiming that "the comparison between the SINS analysis at z ∼ 2 and the AMAZE analysis at z ∼ 3.3 suggests that the fraction of rotating objects does not evolve within this redshift interval" (Gnerucci et al. 2011, p. 19). Unfortunately, they do not provide a full set of mass and v/σ measurements for all the objects in their sample. For fitting purposes, we then just assume another set of measurements at z ∼ 3.3, mirroring the ones at z ≈ 2. Each data sample is divided into galaxy stellar mass (M_*) bins, and for each mass bin we compute the mean v/σ and the dispersion around the mean. We therefore end up with a series of measurements of v/σ in the z–M_* plane, which we then fit with a "Nuker's law" of the form

$$\begin{equation} \frac{v}{\sigma }=a\left(\frac{x}{x_0}\right)^{\alpha }\left[1+\left(\frac{x}{x_0}\right)^\beta \right]^{-\alpha /\beta }, \end{equation} \tag{ 14 }$$

where x = log(1 + z), and x₀ = b(logM_* + c).

Equation (14) describes a broken power law where the slope at x > x₀ is forced to zero and the location of the break, x₀, depends on the stellar mass of the system. The parameter β describes the "sharpness" of the transition between the two power laws. The best least-square fit to the data gives the following values for the five parameters: a = 0.6949; b = 0.4548; c = −8.5255; α = −0.8680; and β = 3.2126.¹¹ Data and best fit are shown in Figure 4. To account for the data dispersion around the best fit, we assume that the density distribution for v/σ is log-normal with standard deviation of 0.34 dex around Equation (14).

Zoom In Zoom Out Reset image size

3.1.2. In Stellar Bulges

The kinematics of stellar bulges has also been extensively studied in the literature. The SAURON project (Cappellari et al. 2007) focuses on a sample of 48 ellipticals and lenticular galaxies, completed with 18 objects from the literature, for a total of 66 objects. Kinematic measurements are also available for stellar bulges of spiral galaxies; Fabricius et al. (2012) report measurements of v/σ for bulges in 43 spiral galaxies of all classes, from S0 to Sc, also making a distinction between bulges and "pseudobulges." Using the total sample of 109 objects, we note that, in general, a distinction can be made between ellipticals and lenticular/spirals. The former have lower v/σ, clustering at zero values, whereas the latter have generally higher v/σ, extending to values greater than 1 and with an average greater than 0.5. Moreover, "pseudobulge" spirals tend to rotate faster than classical bulges, as shown in Figure 5. The v/σ distribution is analytically fitted with the following functions:

• 1.  
  ellipticals

  $$\begin{equation} f(x)=ae^{bx^c}, \end{equation} \tag{ 15 }$$

  with a = 14.59, b = −3.98, c = 0.96;
• 2.  
  bulges and pseudobulges in spirals

  $$\begin{equation} f(x)=ax^be^{cx^d}, \end{equation} \tag{ 16 }$$

  with a = 42.24, b = 1.29, c = −2.95, d = 2.81;
• 3.  
  bulges in spirals

  $$\begin{equation} f(x)=\frac{a}{b}x^be^{-0.5[(x-c)/b]^2}, \end{equation} \tag{ 17 }$$

  with a = 2.10, b = 0.24, c = 0.47; and
• 4.  
  pseudobulges in spirals

  $$\begin{equation} f(x)=\frac{a(x-b)}{x^2+cx+d}, \end{equation} \tag{ 18 }$$

  with a = 0.083, b = −0.53, c = −1.38, d = 0.50.

Zoom In Zoom Out Reset image size

The distribution of v/σ for the different galaxy samples, together with the corresponding f(x) fitting formulae given by Equations (15)–(18), is shown in Figure 5. When implemented in our code, we normalize each f(x) so that $\int _0^\infty f(x)dx=1$ to get a probability density function (PDF). The v/σ of each individual galaxy is then drawn randomly from the appropriate PDF according to its morphological type, as described in detail in the next section.

3.2. Implementation in Our Galaxy Formation Model

3.3. Linking v/σ to the Isotropy Parameter F

When the dispersion of the low-angular-momentum reservoir available for accretion is known, one can compute the amount of prograde/retrograde accretion and derive the MBH spin evolution (see Section 2.2). The relation between the v/σ ratio of the accretion flow and the geometrical parameter F introduced earlier to describe the anisotropy of the angular momentum distribution of the accreting clouds can be derived as follows. Suppose that the net rotation is along the z axis. Each lump of matter falling into the MBH has an angular momentum with z-component L_z = mrsin θ(v + σ_φ), where θ and φ are the lump's polar and azimuthal angles in a spherical coordinate system, and σ_φ is the projection of $\boldsymbol{\sigma }$ on the tangential direction $\boldsymbol{e}_\varphi$. Let us now write σ_φ = σcos α, where α is the angle between $\boldsymbol{\sigma }$ and $\boldsymbol{e}_\varphi$. If $\boldsymbol{\sigma }$ is isotropically distributed, cos α has a uniform probability distribution, i.e., dP/dcos α = 1/2. If v/σ > 1, L_z is always positive and F = 1. If 0 ⩽ v/σ < 1, the probability of having L_z > 0 (or L_z < 0) is proportional to 1 + v/σ (or 1 − v/σ), hence F = (1 + v/σ)/2. In summary, we assume

$$\begin{equation} \begin{array}{lll}v/\sigma >1 &\rightarrow &F=1\,;\\ v/\sigma =0 &\rightarrow &F=0.5\,;\\ 0<v/\sigma <1 &\rightarrow &F=(1+v/\sigma)/2\,. \end{array} \end{equation} \tag{ 19 }$$

One subtlety should be noted about Equation (19) and its derivation. When using the parameter F resulting from this equation in Equation (13), one is implicitly identifying the reference system used above to derive Equation (19) (i.e., one where the z axis is along the total angular momentum of the gas reservoir, out of which the accreting lumps are drawn) with the reference system used to derive Equation (13) (i.e., one where the polar angle θ_out is defined relative to a z axis parallel to the MBH spin). Therefore, we are implicitly assuming that the MBH spin is exactly aligned with the average angular momentum of the clouds already when the first cloud hits the MBH. This may not be the case under many circumstances, but the impact of this assumption on our predictions for the MBH spin magnitudes is expected to be small for the following two reasons: (1) When the MBH is light enough (M_bh ≲ 10⁷ M_☉) to be in the J_disk/2J_bh > 1 regime, the MBH spin realigns with each cloud-angular momentum on a timescale shorter than the accretion time. Therefore, assuming the MBH spin to be always aligned with the clouds average rotation axis does not significantly affect the evolution of the spin magnitude. (2) The high-mass regime of our model is exact for the two limiting cases F = 0.5 and F = 1 and also preserves the correct monotonic trend for 0.5 < F < 1 (i.e., the higher F, the higher the MBH spin magnitude). It would indeed be possible to estimate the detailed effect of the initial MBH orientation for 0.5 < F < 1 by comparing Figure 3 with Figures 6 and 7 in D13 for M_bh ≳ 10⁹ M_☉, but we postpone a more detailed analysis to a future work.

As in B12, we anchor our model to a number of observables, with particular attention to quantities characterizing the MBH population and their accretion properties in the local z ∼ 0 universe, where all the currently available MBH spin measurements are performed.

More precisely, we calibrate our free parameters (see Table 1) against observational data for the QSO luminosity function at z = 0.1 (Figure 6), the MBH mass function at z = 0 (Figure 7), the local "Magorrian" relation between the MBH mass and the bulge dynamical mass in ellipticals (Figure 8), the star formation density for 0 ⩽ z ≲ 5 (Figure 9), the stellar mass function at z = 0 (Figure 10), and the observed fraction of galaxies with given morphology at z = 0 (Figure 11). We stress, however, that the results of our model (and the overall conclusions of this paper) do not depend strongly on the values of these free parameters. Also, the results shown in this section are for the pseudobulge model (see Section 3.2), but the disk and hybrid models produce similar results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. The Calibrated Values of the Free Parameters of the Model (See B12 and Section 2.2 for Their Meaning)

	Light Seeds	Heavy Seeds
Mcloud	3 × 104 M☉	3 × 104 M☉
SN, b	0.4	0.4
SN, d	0.1	0.1
fjet	10	10
Ares	6 × 10−3	5.75 × 10−3
AEdd	2.2	1
kaccr	10−3	10−3

Download table as:  ASCIITypeset image

As can be seen from Figure 6, our model gives reasonable predictions for the z = 0.1 QSO luminosity function when compared to the compilation of observational data by Hopkins et al. (2007), despite slightly underestimating it at low luminosities (the disagreement at the high-luminosity end is less statistically significant). The vertical error bars on the data are from Hopkins et al. (2007), while for the horizontal ones we have assumed −0.5 dex and +0.1 dex to account for possible overestimation of the bolometric correction (Lusso et al. 2012). Similarly, Figure 7 shows that our predictions for the z = 0 black hole mass function are in agreement with recent estimates by Shankar (2013). While the reconstruction of the MBH mass function performed by Shankar (2013) is still preliminary and based only on bulge-dominated galaxies, from a comparison to the observational estimates of Shankar et al. (2009; also shown in Figure 7), it is clear that systematic errors still dominate the determination of the mass function at the high-mass end and that our model's prediction is clearly within the observational uncertainties. Figure 8 shows the agreement between a compilation of dynamical bulge masses and black hole masses in elliptical galaxies form McConnell & Ma (2013; see also Häring & Rix 2004; Magorrian et al. 1998) and our predictions for ellipticals.

The improved star formation implementation described in Section 2.1 allows our model to better reproduce data for the star formation history than the original prescription of B12. Figure 9 compares our predictions to reconstructions based either on observations of the stellar mass density and assuming an evolving initial mass function (IMF; Wilkins et al. 2008) or to fits to instantaneous star formation density observations, assuming either a universal or an evolving IMF (Wilkins et al. 2008). Our predictions for the local galaxy mass function (Figure 10) are instead slightly higher than the observational determinations by Bernardi et al. (2010, 2013) at the low- and high-mass ends. The same problem is often seen in galaxy formation models, which tend to overproduce stars (see, e.g., De Lucia & Blaizot 2007; Guo et al. 2011; Khandai et al. 2014), possibly because of overly simplified treatments of AGNs and supernova feedback. More sophisticated feedback schemes might improve the match with observations, and we plan to investigate this issue further in future work. However, one should also keep in mind that determinations of the mass function are typically obtained by estimating galaxy masses through their luminosity, which results in large systematic uncertainties (especially at the high-mass end), due to the assumptions on the stellar mass-to-light ratio, as well as the different possible ways in which one can fit the light profiles (Bernardi et al. 2013). Finally, Figure 11 shows the predicted fraction of ellipticals, spirals, and irregulars as a function of stellar mass, together with the observational estimates of Conselice (2006). More specifically, following Guo et al. (2011), we use the fraction of total baryonic mass in the bulge, B/T = M_b/(M_b + M_d) (with M_b and M_d the bulge and disk masses), to discriminate the various morphologies; that is, we classify a galaxy as an elliptical when B/T > 0.7, as a spiral when 0.03 < B/T < 0.7, and as an irregular when B/T < 0.03.¹⁵ In spite of this simplistic classification, our model reproduces, at least qualitatively, the observational results. Because in Section 6 we will be dealing with accreting systems, in Figure 12 we compare the Eddington ratio, f_Edd, distribution predicted by our model to observations. In the top panel, we compare the local cumulative distribution of f_Edd as measured by Heckman et al. (2004), who inspected 23,000 active galaxies in a complete sample of 123,000 galaxies from the Sloan Digital Sky Survey in the redshift range 0.05 < z < 0.2 to the predictions of our pseudobulge model at z = 0.1 (other models yield similar results). The model reproduces the f_Edd observed trend satisfactorily, even though it slightly underpredicts the fraction of accreting systems at f_Edd > 0.01. The lower panels highlight the difference between the model and the data at f_Edd ≈ 0.1. The figure shows that our model catches the main trends in the low-redshift Eddington ratio distributions, although not perfectly. A better calibration of the model against this quantity may be considered in future work.

Zoom In Zoom Out Reset image size

The spin evolution of each MBH is dictated by a complex interplay of (1) the variable supply of mass available for accretion, (2) the time-changing host galaxy properties and morphology, and (3) MBH binary coalescence events. We will first inspect some examples of individual MBH–galaxy evolutionary tracks to illustrate the effect of those ingredients (Section 5.1). We stress that the few examples we consider do not represent the whole variety of evolutionary paths predicted by our model. Predictions on statistically significant samples of objects will be presented in Section 5.2.

5.1. Spin Evolution of Individual MBHs

Figure 13 gives three examples of evolution of z = 0 spirals for the pseudobulge model. More specifically, we show the evolution of the MBH mass (with MBH mergers shown with filled circles) along the main MBH progenitor history (i.e., at each MBH merger we follow the more massive progenitor back in time). Also shown are the MBH spin, the mass of the reservoir available for accretion onto the MBH, as well as the host galaxy morphology (with the redshift intervals in which the galaxy is elliptical, i.e., B/T > 0.7, marked by a lavender-shaded area).

Zoom In Zoom Out Reset image size

The two right panels show two MBHs accreting at f_Edd ≈ 1.6 × 10⁻² at z ≈ 0. Despite the similarity of the accretion rates and galaxy morphologies at z ≈ 0, the final MBH spins are quite different. In both cases, the spin grows rapidly early on and remains almost maximal as long as the MBH mass is ≲ 10⁶ M_☉. This is because, as discussed earlier, for such low masses the MBH spin is smaller than the typical angular momentum of the accretion disk arising from the capture of a cloud, and as a result the Bardeen–Petterson effect aligns the MBH spin to the disk's angular momentum, thus making accretion effectively coherent. However, the later evolution of the two cases is quite different and reflects the wide kinematic range of bulges and pseudobulges shown in Figure 5, with v/σ above (below) 1 resulting in high (moderate) MBH spins. For instance, in the bottom right panel the spin first drops to ≈0.25 and then rises again to ≈0.8 following a series of MBH mergers and accretion events (triggered by disk instabilities and major mergers). A major galactic merger at z ≈ 0 then triggers a MBH merger and an accretion event that result in a final spin ≈0.75. In the top right panel, instead, accretion and mergers also cause the spin to drop to ≈0.4 at z ≈ 4, but the degree of anisotropy of the MBH fueling (i.e., in the pseudobulge model, the value of v/σ of the host galaxy's stellar population) is never large enough to allow for very large spins in the subsequent evolution. In the left panel, we show a highly sub-Eddington system, which has a final MBH spin of ≈0.5. In this case, a major merger at z ≈ 2 changes the morphology of the galaxy from spiral to elliptical, and the following highly "incoherent" (i.e., with v/σ ≪ 1) accretion causes a dramatic spin down of the MBH. After that episode, the galaxy reacquires a substantial disk because of minor mergers and accretion/cooling of gas, thus becoming a spiral again, but the accretion rate onto the central MBH never becomes high enough to substantially change its spin.

These trends can be compared to those found in the disk model, shown in Figure 14. In the top right panel we consider a galaxy containing an accreting MBH with M_bh ∼ 10⁶ M_☉ at z = 0. For such a small mass, the distinction between ellipticals and spirals and between the various accretion models (e.g., disk versus pseudobulge) is unimportant because the Bardeen–Petterson effect always aligns the MBH spin to the angular momentum of the incoming cloud. This makes accretion essentially coherent and results in almost maximal spins. The bottom left panel shows again an accreting MBH in a spiral galaxy, but with larger final mass (M_bh ∼ 10⁸ M_☉ at z = 0). Again, accretion is effectively coherent (and the spin close to maximal) until the MBH grows larger than ∼10⁶ M_☉. After that, the spin evolution during accretion events depends on the galaxy morphology. In particular, at z ≈ 2 a major galactic merger (not accompanied by a MBH merger) triggers a morphology change to elliptical, as well as a QSO event that pushes the MBH spin to a_bh ≈ 0.3. The galactic disk then regrows via minor mergers and accretion/cooling of gas, and accretion becomes "more coherent" (because the gas in spirals typically has v/σ > 1 at low redshift, see Figure 4). As a result, the MBH spin starts growing again, and at z ≈ 0 an ongoing accretion event pushes it up to almost maximal values. The bottom right panel shows instead a nonaccreting (at z ≈ 0) MBH hosted in an elliptical galaxy. Again, the spin remains close to maximal until M_bh ∼ 10⁶ M_☉, after which the spin evolution is dominated by accretion, which tends to spin the MBH down when the host is elliptical (see the QSO event at z ≈ 2). When the host is a spiral, instead, the spin tends to mildly grow during accretion events.

Zoom In Zoom Out Reset image size

From the above examples we can draw a couple of conclusions. Nonaccreting systems (which we arbitrarily define as MBHs with f_Edd < 0.01) typically retain memory of the last "violent" event that changed the central MBH spin. Since this is likely to be a major merger/disk instability event that temporarily changes the galaxy morphology to bulge-dominated, the MBH spin will likely be small (nonmaximal, in any case), whether the galaxy later retains a bulge-like morphology or acquires a disk again. Highly accreting systems (i.e., with f_Edd > 0.01) will likely have their spin affected by the ongoing accretion episode. In this case, the spin at z = 0 depends on the typical v/σ ratio of the host galaxy. Spirals in the pseudobulge model have a range of possible spins mirroring the wide spread in v/σ. Conversely, spirals in the disk model are almost certain to have nearly maximally spinning MBHs. On the other hand, MBHs in accreting ellipticals tend to have low spins, because of their typically low v/σ ratios, as shown in Figure 5.

5.2. Spin Evolution of the Global MBH Population

Having discussed how accretion, mergers, and galaxy morphology contribute to the MBH spin buildup in individual galaxies, we now turn to the investigation of the cosmic spin evolution of the whole MBH population. In the three panels of Figure 15 we show the spin distributions predicted by the pseudobulge (top plot), disk (central plot), and hybrid (bottom plot) models for selected subsamples of galaxies at different redshifts.

Zoom In Zoom Out Reset image size

As mentioned earlier, all three models adopt the same prescription for the spin evolution in elliptical galaxies (right column in each panel), which results in very similar distributions that do not evolve much from z = 5 to the present epoch. For M_bh ≲ 10⁶ M_☉, J_disk/2J_bh > 1 and MBHs tend to have maximal spins. At large masses, J_disk/2J_bh ≪ 1, and a_bh tends to the equilibrium value given by Figure 3. The median v/σ of elliptical galaxies is ≈0.2, corresponding to F ≈ 0.6, which translates to a median equilibrium spin of 0.3, in good agreement with Figure 15. The evolution of the overall population of spiral galaxies (central column) also shows similarities in all three models. This class of systems experiences more pronounced redshift evolution, and the median of the spin distribution shows a dip between 10⁷ M_☉ and 10⁸ M_☉ appearing at z ∼ 3. The overall distributions are not very different from those of the ellipticals, but spins are on average larger by virtue of the typically higher v/σ, as shown in Figure 5.

The subset of accreting MBHs (i.e., with f_Edd > 0.01) in spirals (right column) shows instead a remarkably different evolution between the pseudobulge and disk models. In the pseudobulge model, these galaxies maintain, on average, higher spins at all redshifts compared to other systems (see discussion in Section 5.1), but the trend of decreasing spin with increasing mass is preserved. Conversely, in the disk model, MBH spins show a stronger redshift evolution, and at z < 2 they tend to be maximal independently of the MBH mass. This is because MBHs are efficiently accreting gas with v/σ > 1 (see Figure 4). As expected, the hybrid model shows an evolution which is halfway between the disk and pseudobulge models.

Our predicted spin distributions are different for different classes of galaxies, which is a testable prediction. In general, MBHs in ellipticals have lower spins than their spiral counterparts, but still with an average value of a_bh ≈ 0.4 and a long tail extending to much higher values. We note that our findings do not rule out spin-powered jet models (see, e.g., Sikora et al. 2007). Indeed, in such models the power of the jet is usually $\propto a_{\rm bh}^2$ (Blandford & Znajek 1977; Tchekhovskoy et al. 2010; Narayan et al. 2014), implying a difference of just a factor of a few in luminosity between maximally and mildly spinning MBHs. It is also clear that in our model, accreting MBHs tend to be biased toward higher spins than nonaccreting ones, especially if located in spiral galaxies. In addition, current measurements of the MBH spins require large X-ray fluxes, which are only possible if the MBH shines at a significant fraction of the Eddington luminosity. Therefore, the current observed sample of MBH spins may not be an unbiased indicator of the overall cosmic MBH population. In the following we will perform a statistical analysis comparing our predictions to the MBH spin measurements available today at z < 0.1, selecting the appropriate sample in our simulations in order to account for both these selection effects.

We carry out here a quantitative comparison between MBH spin measurements and the output of our theoretical models, which link the MBH spin evolution to the properties of the galaxy host. However, existing spin measurements are sparse and are affected by several statistical and systematic uncertainties that are often difficult to quantify. The results of the following analysis should therefore be taken as indicative. Nonetheless, we will show that even with the current data, we can rule out some mass and spin growth scenarios and gain some qualitative insights about the connection between spin evolution and the properties of the galaxy host supplying gas to the accretion flow.

6.1. Observed Sample

MBH spins have been measured via Kα reflection broad line modeling for about 20 objects. We took data from the sample compiled by Reynolds (2014), integrated with objects from Brenneman (2013). We noted some (minor) discrepancies in the numbers reported by the two authors; when necessary we inspected the original measurement papers (see references in Reynolds 2014 and Brenneman 2013) and obtained the values directly from there. All the relevant properties of the sample are shown in Table 2. A meaningful comparison to a theoretical model is possible only if (1) observations have meaningful errors and (2) the general properties of the observed subsample can be isolated, and selection effects are understood. Both issues are somewhat tricky here.

Table 2. Sample of MBHs with Spin Measurements from the Kα Reflection Line

Object Name	Galaxy Type	z	LX(erg s−1)	fEdd	log(Mbh(M☉))	Spin	Adopted PDF
1H0707−495	⋅⋅⋅	0.0411	3.7 × 1043	1.0	6.70 ± 0.4	>0.97	Flat [0.97,0.998]
Mrk1018	S0	0.043	9.0 × 1043	0.01	8.15	$0.58^{+0.36}_{-0.74}$	Flat [0,0.94]
NGC 4051	SAB(rs)bc	0.0023	3.0 × 1042	0.03	6.28	>0.99	Flat [0.99,0.998]
NGC 3783	SB(r)ab	0.0097	1.8 × 1044	0.06	7.47 ± 0.08	>0.88	Flat [0.88,0.998]
1H0419−577	⋅⋅⋅	0.104	1.8 × 1044	0.04	8.18 ± 0.05	>0.89	Flat [0.85,0.998]
3C120	S0	0.033	2.0 × 1044	0.31	$7.74^{+0.20}_{-0.22}$	>0.95	Flat [0.95,0.998]
MCG-6-30-15	E/S0	0.008	1.0 × 1043	0.4	6.65 ± 0.17	>0.98	hGauss [0.998,0.01]
Ark564	SB	0.0247	1.4 × 1044	0.11	<6.90	$0.96^{+0.01}_{-0.06}$	hGauss [0.96,0.04]
TonS180	⋅⋅⋅	0.062	3.0 × 1044	2.15	$7.30^{+0.60}_{-0.40}$	$0.91^{+0.02}_{-0.09}$	hGauss [0.94,0.067]
RBS1124	⋅⋅⋅	0.208	1.0 × 1045	0.15	8.26	>0.97	hGauss [0.998,0.02]
Mrk110	⋅⋅⋅	0.0355	1.8 × 1044	0.16	7.40 ± 0.09	>0.89	Gauss [0.945,0.033]
Mrk841	E	0.0365	8.0 × 1043	0.44	7.90	>0.52	Gauss [0.80,0.17]
Fairall9	Sc	0.047	3.0 × 1044	0.05	8.41 ± 0.11	$0.52^{+0.19}_{-0.15}$	Gauss [0.6,0.1]
SWIFTJ2127.4+5654	SB0/a(s)	0.0147	1.2 × 1043	0.18	7.18 ± 0.07	0.6 ± 0.2	Gauss [0.6,0.1]
Mrk79	SBb	0.0022	4.7 × 1043	0.05	7.72 ± 0.14	0.7 ± 0.1	Gauss [0.7,0.1]
Mrk335	S0a	0.026	5.0 × 1043	0.25	7.15 ± 0.13	$0.83^{+0.09}_{-0.13}$	Gauss [0.81,0.067,<0.92]
Ark120	Sb/pec	0.0327	3.0 × 1045	1.27	8.18 ± 0.12	$0.64^{+0.19}_{-0.11}$	Gauss [0.68,0.093]
Mrk359	pec	0.0174	6.0 × 1042	0.25	6.04	$0.66^{+0.30}_{-0.54}$	Gauss [0.66,0.33,<0.96]
IRAS13224−3809	⋅⋅⋅	0.0667	7.0 × 1043	0.71	7.00	>0.987	Gauss [0.989,0.002]
NGC 1365	SB(s)b	0.0054	2.7 × 1042	0.06	$6.60^{+1.40}_{-0.30}$	$0.97^{+0.01}_{-0.04}$	Gauss [0.97,0.03,<0.98]

Notes. The table is compiled using objects from Reynolds (2014) and Brenneman (2013); galaxy redshifts are taken from the NED database (http://ned.ipac.caltech.edu/). Quoted errors in the mass and spin measurements correspond to the 68% and 90% confidence levels, respectively. The spin PDF we use in our statistical analysis is given in the last column, which identifies three different functional forms: flat (min and max value given in [ ]), half-Gaussian (maximum and σ given in [ ]), and Gaussian (maximum and σ given in [ ]). In the latter case, a third number in [ ], when present, defines a sharp cutoff in the PDF.

Download table as:  ASCIITypeset image

Errors quoted in Table 2 represent the 68% confidence level in the mass and the 90% confidence level in the spin measurements. A complete knowledge of the PDF of those quantities would be desirable, but such information is often missing in the literature, and we can at best put forward educated guesses. For the (log of the) MBH mass, we consider (1) a Gaussian PDF with σ equal to the quoted error when the latter is symmetric with respect to the best measured value, (2) a flat PDF within the given errors when those are asymmetric, and (3) a Gaussian PDF with an arbitrary σ = 0.3 when errors are absent. Spin PDFs are instead derived by visually inspecting the function χ²(a_bh), which is always given in the relevant papers. We identify three families of measurements: (1) sometimes χ²(a_bh) has an approximately flat minimum, and in this case we take a flat PDF in the corresponding range; (2) for a few objects, χ²(a_bh) is extremely skewed on the left of the minimum, and in such cases we take as PDF the left half of a Gaussian distribution; (3) more often, χ²(a_bh) is quite symmetric (at least in the 99% confidence region), calling for a Gaussian model of the PDF. Details of the PDFs are given in the last column of Table 2. Being aware of the arbitrary nature of this procedure, we also consider an alternative model in which we take all errors in masses and spins to have a flat PDF within the range quoted in Column 7 of Table 2. We show in Appendix C that our results are largely independent of the adopted shape of the PDFs.

The obvious selection effect for a Kα line measurement is that the source has to be bright in hard X (see also discussion in Brenneman et al. 2011); in fact, most of the MBHs in the sample reside in Seyfert 1 or narrow-line Seyfert 1 (NLS1) galaxies. Hosts are usually spirals or lenticular galaxies, with the exception of Mrk841. The observed sample is neither flux- nor volume-limited. X-ray luminosities and associated MBH masses have distributions consistent with being log-flat in the ranges 10⁴²–10⁴⁵ erg s⁻¹ and 10⁶–3 × 10⁸ M_☉, respectively. Eddington ratios f_Edd are also evenly distributed in the range 0.01–1. The redshift distribution is also approximately flat in the range 0 < z < 0.07, with two outliers at z > 0.1. Given these facts, the best-matching sample in our model is provided by spiral galaxies (defined by B/T < 0.7) with f_Edd > 0.01 at low redshift. For the purpose of the analysis, we will therefore exclude Mrk841 (which is an elliptical) and include all galaxies in the observed sample with unknown morphology (thus implicitly assuming they are spirals/lenticulars). Discarding the latter would weaken our results (because there would only be 12 objects in the observed sample) but not qualitatively change them.

6.2. Theoretical Sample

Theoretical distributions are computed with our model on a grid in the M_bh − a_bh parameter space. The MBH mass range 10⁶ M_☉ < M_bh < 5 × 10⁸ M_☉ is divided in five equally log-spaced bins, and for each bin, spin distributions are computed on 20 linear bins covering the range 0 < a_bh < 1. We applied Gaussian kernel smoothing to each measured spin with σ = 0.05 to get smoother distributions¹⁶. To check the robustness of our results against spin binning, we also constructed two theoretical distributions using 10 and 30 linear a_bh bins and one distribution considering 20 bins both in log mass and spin. As a sanity check, a comparison among the results obtained with different binning choices is performed in Appendix C.

As mentioned above, all (but one) MBHs with spin measurements are accreting at f_Edd > 0.01 and reside in spiral or lenticular galaxies at z < 0.1. Those systems, although present in fair number, are necessarily rare in our runs (because accreting systems become much rarer than nonaccreting ones at low redshift, see Figure 12), which makes it difficult to construct a smooth two-dimensional (2D) mass-spin distribution to compare to observations, even by adding together our runs for light and heavy seeds. Figure 15, however, shows that the overall spin evolution of the galaxy population is small at z < 1, as a consequence of the reduced galactic activity both in terms of gas cooling/star formation and mergers. Therefore, to improve the statistics of our analysis, we compare observations to z = 0.1 to a fiducial theoretical sample taken at z = 1. For the sake of completeness, effects due to the evolution of the galaxy population at z < 1 are investigated in Appendix C, where we show that they do not significantly affect the conclusions of our analysis.

The predicted spin distributions of accreting MBH hosted in spirals at low redshift are shown in Figure 16, together with the measurements of Table 2. The pseudobulge model qualitatively reproduces the observed trends, with most of the measurements and relative error bars falling in the 95% confidence region predicted by the models. The disk model displays a slightly poorer match, with many events outside the 68% confidence region predicted by the model and just a minor hint of the general spindown observed with increasing MBH mass. The hybrid model nicely matches the data, reproducing both the observed high-spin subsample and the objects with moderate spin.

Zoom In Zoom Out Reset image size

6.3. The 2D Kolmogorov–Smirnov Test

6.4. Bayesian Model Selection

The 2D KS test provides useful indications about the consistency of the data with the predictions of a given model. It is, however, not very practical in assessing which one is the best model among a set of different options. This latter issue can be tackled within the framework of Bayesian model selection by computing the posterior probability of the parameters of a certain model, given a set of data. Because our models are nonparametric, this boils down to the evaluation of the odds ratio for different pairs of models, which we now describe.

Following Sesana et al. (2011), we assume that the observations are uncorrelated, so that the number of objects, n_i, measured in a particular ΔM_bhΔa_bh bin in parameter space can be drawn from a Poisson probability distribution with parameter r_i equal to the bin-integrated rate:

$$\begin{equation} p(n_i) = \frac{(r_i)^{n_i} {\rm e}^{-r_i}}{n_i!}. \end{equation} \tag{ 20 }$$

If we divide the parameter space up into a certain number K of bins, then the information that comes from the data (D) is the number of events in each bin (n_i). The overall likelihood p(D|X) of seeing this data under the model X is the product over all K bins of the Poisson probabilities to see n_i events in the ith bin, given the rate r_i(X) predicted by model X:

$$\begin{equation} p(D|X) = \prod _{i=1}^K \frac{(r_i(X))^{n_i} {\rm e}^{-r_i(X)}}{n_i!}. \end{equation} \tag{ 21 }$$

It is straightforward to take the limit of this expression as the bin sizes tend to zero to derive a continuum version of this equation (Gair et al. 2010), but in this analysis we will stick to binned distributions for simplicity. In the presence of measurement errors, the jth event can then be assigned to different bins, with a probability given by its PDF ρ_j(M_bh, a_bh). We can therefore construct the entire set of possible arrays of events falling in each bin, n_i, with their relative probabilities. Each array is then analyzed and weighted according to its probability to get the overall likelihood of the sample, given the model.

If we want to compare two competitive models, A and B, Bayes' theorem allows us to assign to model A a probability

$$\begin{equation} p(A|D) = \frac{p(D|A)P(A)}{p(D|A)P(A)+p(D|B)P(B)}, \end{equation} \tag{ 22 }$$

whereas the probability of model B is just obtained by swapping A and B in Equation (22). Here, p(D|X) is the likelihood given by Equation (21), and P(X) is the prior probability assigned to model X. The odds ratio of model A over model B is

$$\begin{equation} \Lambda _{AB} = \frac{p(D|A)}{p(D|B)} \frac{P(A)}{P(B)}. \end{equation} \tag{ 23 }$$

If we do not have any reason to prefer a priori model A to model B, then P(A) = P(B) = 0.5, and the odds ratio becomes the likelihood ratio. Moreover, Equation (22) implies that in a two-model comparison we can simply assign a "relative probability" p_A = p(D|A)/(p(D|A) + p(D|B)) to model A and p_B = 1 − p_A to model B.

We compute the odds ratios of our three fiducial models: pseudobulge (p), disk (d), and hybrid (h) according to Equation (23), where we assume P(A) = P(B) = 0.5. Likelihoods are computed according to Equation (21), folding in measurement errors. The shapes of the error PDFs of each single observation are described in Section 6.1, and we factorize ρ_j(M_bh, a_bh) = g(M_bh)h(a_bh), assuming uncorrelated mass and spin measurements.

Results are shown in Table 4. We get log Λ_hd ≈ 2.5, implying that the hybrid model provides a better description of the data than the disk model, at the >99.5% confidence level. This strengthens the results of the 2D KS test, providing compelling evidence in favor of the hybrid model. The pseudobulge model sits somewhat in the middle; it is preferred at the ≈95% level over the disk model, but it is disfavored at about the same level with respect to the hybrid model. To summarize, the odds ratio analysis provides moderate (decisive) evidence that the hybrid model is a better description of the data than the pseudobulge (disk) model.

Table 4. Model Selection Results: Pair Comparisons between Our Three Fiducial Spin Evolution Models, Pseudobulge (p), Disk (d), and Hybrid (h)

Model Pairs	Odds Ratios & Probabilities
Hybrid/pseudobulge	logΛhp = 1.0804	ph = 0.9233	pp = 0.0767
Hybrid/disk	logΛhd = 2.4749	ph = 0.9966	pd = 0.0034
Pseudobulge/disk	logΛpd = 1.3944	pp = 0.9612	pd = 0.0388

Notes. For each two-model comparison we report the log of the likelihood ratio Λ_AB, i.e., the ratio between the probability of model A (p_A) and model B (p_B).

Download table as:  ASCIITypeset image

In this paper, we presented the results of a semianalytical model for the evolution of galaxies and MBHs. The model keeps track (although in simplified ways) of the morphology of the galaxies, as well as the MBH masses and spins. For the first time we link the dynamical properties of the gas fueling the MBHs to the host galaxy properties through different observationally based prescriptions and derive predictions testable with existing observations. We stress that all the other similar investigations predicting MBH spin distributions either assumed "chaotic" accretion (MBHs accreting small gas clouds with isotropically oriented angular momenta) or "coherent" accretion (MBHs accreting all the time on a fixed plane).

Our model predicts different spin distributions for different types of galaxies. To date, only ≈20 MBH spins have been directly measured through Kα iron line fitting (Reynolds 2014; Brenneman 2013). All these MBHs (but one) are hosted in low-redshift late-type galaxies, preventing us from testing our results for different galaxy types. We therefore tested our model by selecting low-redshift accreting MBHs hosted in spirals. For this galaxy class, we have assumed three different gas dynamics prescriptions. In the pseudobulge model we assumed that the dynamics of gas fueling the central MBHs is similar to that of the stars in both galaxy bulges or pseudobulges, according to the nature of the host's spheroidal component. In the disk model, instead, the fueling gas has the more coherent dynamics of the large-scale gaseous disk. The hybrid model shares the same gas dynamics as the disk model for isolated galaxies but assumes a gas dynamics similar to that of stars in spheroids during merger-driven accretion events (see Section 3 for details).

Our models predict interesting features in the MBH spin distribution, and a statistical comparison to observations yields a number of interesting results.

• 1.  
  Different galaxy morphologies result in different MBH spin distributions. Because the geometrical properties of the accretion flow are related to the large-scale kinematics of the galaxy, MBHs hosted in ellipticals tend to have lower spins than those hosted in spirals.
• 2.  
  In general, MBH spins are a decreasing function of the MBH mass: ∼10⁶ M_☉ MBHs tend to be maximally spinning, whereas for masses >10⁸ M_☉ a wide range of values is possible for the spin, depending on the host morphology and MBH accretion rate.
• 3.  
  Accreting MBHs in spirals, i.e., those matching the observational sample, tend to spin fast and do not provide an unbiased indicator of the underlying spin distribution of the overall MBH population.
• 4.  
  The 2D KS tests show that, in general, the hybrid and pseudobulge models are consistent with observations, while the disk model is disfavored. A quantitative assessment of the compatibility of the disk model with the data, however, strongly depends on the maximum value to which MBHs can be spun up by accretion. This issue is discussed in detail in Appendix C, where we show that if this value is a_bh = 0.998 (Thorne 1974), the disk model can be safely ruled out at the >99% confidence level.
• 5.  
  A likelihood odds ratio study shows that the hybrid model—which takes into account the possible angular momentum reshuffling due to merger-driven phases of nuclear activity—provides the best match to the observed data. In a two-model comparison, (1) it is strongly favored (at the >99% confidence level) over the disk model; and (2) it is moderately favored (only at about the 90% confidence level) over the pseudobulge model.
• 6.  
  For the sake of comparison, in Appendix B we have also run three separate sets of simulations adopting the standard coherent-accretion model, as well as two flavors of the chaotic-accretion model. These models provide a much worse description of the data; the coherent model shows features similar to the disk model, but with a poorer agreement to the data, while the chaotic models are unambiguously ruled out by observations.

We emphasize that the study of the evolution of the MBH spin distribution in a cosmological framework is still in its infancy. Because our model is idealized from many points of view, here we highlight some still-needed theoretical and observational improvements. (1) The treatment of the Bardeen–Petterson effect, which plays a major role in determining the fraction of prograde versus retrograde accretion events, is based on simple isotropic-α-viscosity-driven accretion disks, and this assumption has been criticized by Sorathia et al. (2013a, 2013b). The impact of this strong and debated assumption is unknown, since no working models of the MBH spin evolution relaxing it are available to date. (2) We assumed that every gas cloud fueling the MBHs weighs M_cloud ∼ 3 × 10⁴ M_☉. GMCs are observed to have a broad mass spectrum. Observational studies of local galaxies reported typical molecular cloud mass functions with cutoffs at about 10⁵–10⁶ M_☉ (e.g., Williams & McKee 1997; Engargiola et al. 2003; Fukui et al. 2008), making the average cloud mass comparable to the mass that we assume in this paper. We note, however, that the cutoff mass seems to depend on the host galaxy properties (e.g., Hopkins 2012 and references therein). A simple implementation of a "characteristic" cloud mass function is not trivial and beyond the scope of this paper. The effect of different average cloud masses can anyway be easily understood: the cloud mass contributes to setting the maximum MBH mass at which a complete alignment of the spin occurs at every accretion event, effectively making the fueling dynamics unimportant. To first approximation (see D13 for a more detailed discussion), an order-of-magnitude variation of the average cloud mass results in an order-of-magnitude shift along the MBH mass axis of all the spin distributions. The results not showing any significant feature in the spin evolution (e.g., the f_Edd > 0.01 disk model at low redshift, but see also the coherent and chaotic I models discussed in Appendix B) would thus remain mostly unchanged under a variation of the cloud mass. On the other hand, the pseudobulge and hybrid models would predict higher spins for higher masses (lower spins for lower masses) if the cloud mass were higher (lower). We checked that a shift of the spin distributions by one order of magnitude toward larger masses does not affect the compatibility of those two models with the data, yielding a KS probability p_KS > 0.3. A similar shift of the distribution toward lower masses results in a slightly worse description of the data, but still with a KS probability ≳ 0.1. (3) Our dynamical prescriptions for the nuclear gas dynamics are mostly based on gas and/or stellar dynamics observations on scales much larger than the accretion disk or the MBH sphere of influence. Our model would greatly benefit from a description of the gas dynamics on scales ≲ 100 pc, already achievable in some cases with high-resolution ALMA spectroscopy (e.g., Combes et al. 2013, 2014). Finally, (4) in our model we completely neglected the effect of rotational energy extraction from the MBH on the spin. This could in principle power relativistic jets through the Blandford–Znajek mechanism, limiting the maximum spin to lower values in radio galaxies.

From an observational point of view, the main improvement would consist in enlarging the currently small sample of measured MBH spins and extending it to higher redshift and different galaxy types¹⁷. In this regard, the forthcoming Astro-H satellite will have exquisite combined spectral resolution and sensitivity, making Kα-based spin measurements possible perhaps up to z ∼ 1 (Takahashi et al. 2012). Along the same lines, Athena+, to be launched in 2028, will push such measurements to z ∼ 2 for extremely X-ray bright systems, significantly expanding the observed spin sample (Dovciak et al. 2013). On a shorter timescale, the eROSITA satellite will observe >10⁵ AGNs. Image stacking in luminosity and redshift bins will reveal the "average" shape of the Kα line of these systems, making it possible to study trends in the typical spin of AGNs across a large mass range and for redshifts perhaps up to z ∼ 1 (Merloni et al. 2012), thus providing an important benchmark for comparison to theoretical models. On a different note, the eLISA mission, now selected by ESA for the L3 launch slot, will measure the spins of merging MBH binaries across the cosmic history to high precision, providing a sample that will not be biased toward high X-ray luminosities.

Even though the future of spin measurements looks literally bright, we should also keep in mind that the uncertainties on the value of measured spins depend on the technique used. Different methods have been used to constrain the MBH spin distribution at different redshifts and for different masses. Some methods estimate the spin value from the radio properties of the AGNs (e.g., Daly & Sprinkle 2014 and references therein), from the ionizing flux required to produce the observed broad emission lines (Netzer & Trakhtenbrot 2014), or from fitting the continuum from the accretion disk (Laor & Davis 2011b). Some of these methods are still a matter of debate (e.g., Gallo et al. 2012; Raimundo et al. 2012; Laor & Davis 2011a), and all of them require an independent measure of the MBH mass to estimate its spin. In this paper we compared our predictions with the mass-independent measurements obtained through relativistically broadened iron Kα line fitting to prevent any possible spurious mass-spin correlation. We refer to Reynolds (2014) and Brenneman (2013) for a detailed discussion of the uncertainties of these measurements. Since our models predict a significant dependence of the MBH spins on the MBH masses, the uncertainties on the masses are fundamental as well. For example, Mrk 359 has an estimated mass of ≈10⁶ M_☉, obtained through a single-epoch measurement of the continuum and broad Hβ FWHM (Zhou & Wang 2005). The width of the line is only 480 km s⁻¹, significantly smaller than both the typical value for type 1 AGNs and the threshold for being classified as a NLS1. Since the MBH is hosted in a spiral galaxy, the line is peculiarly large to be associated with a single narrow line, and a broad component is expected because of the point-like emission clearly visible in the Hubble Space Telescope imaging (Malkan et al. 1998). The observed line width could be so narrow because of orientation effects, already suggested to be significant for NLS1s (Decarli et al. 2008), causing an underestimation of the MBH mass. (To be on the safe side, we checked that the contribution of Mrk 359 to our analysis is negligible. Its MBH spin is so poorly constrained that removing it from the observational sample has no effect on our results.) More independent measurements of the masses of MBHs with spin estimates would help better constrain their values. The detection of gravitational waves will also greatly improve the situation, giving high-precision measurements of both the MBH masses and spins, which could be compared with predictions for merging systems.

Regardless of the large uncertainties on which we commented above, it is remarkable that two of the models (pseudobulge and hybrid) that we investigated generate samples of MBH masses and spins consistent with observations. Since these models assume that the accreting gas has less coherent dynamics than the larger-scale gas structure (although, in the hybrid case, only in recent merger remnants), our analysis seems to suggest the existence of a physical process that, while decreasing the gas-angular momentum magnitude (triggering accretion in the first place), reshuffles the gas-angular momentum direction as well. Such a reshuffling could be caused by local processes, i.e., star formation and the torques exerted by the gas self-gravity (e.g., Maio et al. 2013), as well as larger-scale gravitational instabilities and violent gas inflows (e.g., Hopkins et al. 2012; Dubois et al. 2014b). On the other hand, in order to achieve such a good agreement between the predictions of the model and the observed spins, we have to require the accreting gas to have, on average, a nonzero angular momentum, i.e., the gas must not accrete isotropically onto the MBH. We have in fact shown that a perfectly isotropic accretion flow would result in significantly lower spins for intermediate to large MBH masses, which is inconsistent with observations (see Appendix B).

We further note that the pseudobulge and the hybrid models, while providing the best descriptions of the data, hint at different (although possibly coexisting) evolutionary scenarios. The hybrid scenario has the highest statistical significance and a simple physical interpretation, with a clear candidate for the source of turbulence. The gas dynamics remains quite coherent for most of the time, which is reasonable because gas can efficiently dissipate turbulent motions on few local orbital times (see the discussion in Maio et al. 2013 and references therein). Only during galaxy mergers the violent reshuffling of the angular momentum allows for more isotropic gas inflows. We note that, strictly speaking, in the hybrid model we adopt a "bulge-like" dynamics for the accreting gas during the whole of the merger-driven accretion events, even though the AGN activity lasts for longer than the galaxy merger itself. Thus, during the tails of the nuclear activity, accretion may be more coherent than in our model's assumptions. However, during a merger-driven AGN phase, most of the gas is accreted in a short, almost-Eddington burst, lasting about a Salpeter time ≲ 0.2 Gyr. The long-lived accretion tails, still observable at low redshifts, present instead a lower Eddington ratio f_Edd ∼ 0.01 and cannot therefore significantly change the MBH spin magnitude.

Although the pseudobulge model is disfavored by our statistical analysis, it still provides a good description of the observational data and is consistent with the evidence that low-redshift AGNs with moderate-to-high Eddington ratios are preferentially hosted in pseudobulges (Heckman & Best 2014). Pseudobulges have recently been associated with stellar bars, i.e., structures able to drive gas inflows toward the galaxy nucleus (see, e.g., Sellwood 2014; Kormendy 2013, for up to date reviews). After the initial instability responsible for the bar formation, the bar itself generally undergoes a thickening that can be either caused by a second "buckling" instability or by a vertical resonance (Sellwood 2014), possibly responsible for the pseudobulge formation (Kormendy 2013). In this scenario, the gas could accrete with a more disk-like dynamics until the bar thickens and fall toward the MBH with a higher-velocity dispersion after the pseudobulge forms if the buckling instability and/or any vertical resonances significantly affect the gas motion.

The paucity of data and the many uncertainties about the different fueling processes do not allow us to firmly reject either of these models. If either model is confirmed by more data, it could point toward a physically motivated galaxy–MBH coevolution scenario for low-redshift disk galaxies, in which the dynamics of gas in nuclear bars and/or mergers drive the MBH mass and spin evolution.

We are indebted to M. R. Krumholz and D. Calzetti for invaluable insights into the metallicity dependency of star formation. We also thank M. Colpi, J. Gair, M. Kesden, A. Merloni, C. Reynolds, A. Tchekhovskoy, and M. Volonteri for insightful comments and discussions. E.B. acknowledges support from the European Union's Seventh Framework Program (FP7/PEOPLE-2011-CIG) through the Marie Curie Career Integration Grant GALFORMBHS PCIG11-GA-2012-321608 and support from the Deutsches Zentrum fur Luft- und Raumfahrt (DLR) through the DFG Grant SFB/TR 7 Gravitational Wave Astronomy. Computations were performed on the GPC supercomputer at the SciNet HPC Consortium, as well as on the "Projet Horizon Cluster" at the Institut d'astrophysique de Paris.

Observations of our and nearby galaxies have established that star formation is ultimately associated with GMCs, giant because they can be very massive (M ∼ 10⁶ M_☉) and extended (∼100 pc). There the star formation rate is determined by the mass fraction f_c in cold gas (generally, but not only, in molecular form) and the timescale t_SF needed to convert it into stars. This latter depends both on the cloud properties, such as density (the higher the density, the faster is the process) and on the presence of feedback-driven turbulence, which slows the collapse. One of the main (observationally supported) assumptions is that the GMC properties are environment-independent, as long as the surrounding ISM has lower pressure than the GMC itself. In particular, the cloud surface density Σ_cl is set by internal processes to be always around Σ_cl ≈ 85 M_☉ pc⁻² in more tenuous environments, where the gas surface density is Σ_g ⩽ 85 M_☉ yr⁻¹ ≡ Σ_th (e.g., Bolatto et al. 2011). Above this threshold, pressure equilibrium between the ISM and the GMC sets Σ_cl ≃ Σ_g.

A.1. Star Formation in the Galactic Disk

In this framework, the local star formation in the galactic disk can be described by

$$\begin{equation} \dot{\Sigma }_* = \frac{f_{\rm c} \Sigma _{\rm g}}{t_{\rm SF}}, \end{equation} \tag{ A1 }$$

where the timescale t_SF has two regimes, according to whether the cloud density Σ_cl is set by the interstellar pressure,

$$\begin{equation} t_{\rm SF}^{-1} = \frac{M_6^{-0.33}}{0.8 \ {\rm Gyr}} \max \left[1,\left(\frac{\Sigma _{\rm g}}{\Sigma _{\rm th}}\right)^{0.67}\right], \end{equation} \tag{ A2 }$$

(Krumholz et al. 2009), where M₆ = M/(10⁶ M_☉) is the GMC mass. The local Jeans mass $M_{\rm j} \approx \sigma _{\rm g}^2/G^2 \Sigma _{\rm g}$ gives a good estimate of the cloud mass, where σ_g is the gas velocity dispersion. In a galactic disk environment, we may assume a condition of marginal gravitational stability and obtain a function of the background gas surface density only,

$$\begin{equation} t_{\rm SF}^{-1} = (\rm 2.6\ Gyr )^{-1} \times \left\lbrace \begin{array}{@{}l l}\left(\frac{\Sigma _{\rm g}}{\Sigma _{\rm th}}\right)^{-0.33}, & \Sigma _{\rm g} < \Sigma _{\rm th},\\ \left(\frac{\Sigma _{\rm g}}{\Sigma _{\rm th}}\right)^{0.34}, & \Sigma _{\rm g} > \Sigma _{\rm th}. \end{array} \right. \end{equation} \tag{ A3 }$$

The cold-gas mass fraction f_c is equivalent to the molecular mass fraction at metallicity greater than ∼1% solar, when H₂ has time to form before collapsing and forming stars. At lower metallicities, instead, star formation will occur in a cold-atomic-gas phase rather than a molecular phase (Krumholz 2012). In general, as the metallicity decreases, the cold gas available for forming stars decreases as well. However, recent observations of nearby spirals and dwarfs (Bigiel et al. 2010) and of the Small Magellanic Cloud (Bolatto et al. 2011) suggest that f_c levels off at around 2% rather than dropping to zero. These considerations led us to adopt the following prescription,

$$\begin{equation} f_{\rm c} = \left\lbrace \begin{array}{@{}l l}1-\left[1+\left(\frac{3}{4} \frac{s}{1+\delta }\right)^{-5}\right]^{-1/5}, & {\rm if}\ f_{\rm c} > 2\%,\\ 2\%, & {\rm otherwise} \end{array} \right. \end{equation} \tag{ A4 }$$

(M. R. Krumholz, 2013 private communication), with

$$\begin{equation*} s = \ln {(1+0.6 \chi)}/(0.04 \Sigma _{\rm 1} Z^{{\prime }}), \end{equation*}$$

$$\begin{equation*} \chi = 0.77 (1+ 3.1 Z^{{\prime }0.365}), \end{equation*}$$

$$\begin{equation*} \delta = 0.0712 (0.1 s^{-1} + 0.675)^{-2.8}, \end{equation*}$$

where Σ₁ = Σ_g/(M_☉ pc⁻²) and Z' is the metallicity in Solar units.

Integrating Equation (A1) over the entire disk surface, we obtain the total star formation rate in the disk.

A.2. Star Formation in the Galactic Bulges

The very same prescription can be applied to the bulge during the periods of "quiescent star formation" in the galaxy, i.e., when violent star-bursting events are not triggered. Practically, we use a volumetric star formation law,

$$\begin{equation} \dot{\rho }_* = \frac{f_{\rm c} \rho _{\rm g}}{t_{\rm SF}}, \end{equation} \tag{ A5 }$$

which we derive directly from Equation (A1), by expressing Σ_g and M as a function of the volumetric gas density ρ_g and the local isothermal sound speed c_s.

Let us start by considering the star formation timescale (Equation (A2)). In regions with densities below threshold, $t_{\rm SF} \propto M_{\rm j}^{0.33}$, where we can simply write the Jeans mass in the more familiar way,

$$\begin{equation} M_{\rm j} = \frac{\pi }{6} \frac{c_{\rm s}^{3}}{G^{3/2} \rho _{\rm g}^{1/2}} \approx 10^6 \, M_{\odot } \left(\frac{c_{\rm s}}{8.3 \ \rm km\ s^{-1}}\right)^{3} \rho _1^{-1/2}, \end{equation} \tag{ A6 }$$

with ρ₁ = ρ_g/(M_☉ pc⁻³). In denser regions of the bulge, one has an additional dependence on Σ_cl, i.e., $t_{\rm SF} \propto M_{\rm j}^{0.33} \Sigma _{\rm cl}^{-0.67}$ (see Equation (A2)). To relate the surface density to the volume density, let us note that for a spherical cloud of mass M and characteristic size L, ρ_cl ≈ M/L³ and L ≈ Σ_cl/ρ_cl. Eliminating L in these two expressions, we obtain $\rho _{\rm cl} \approx M_{\odot }\ {\rm pc}^{-3} \left(\Sigma _{\rm cl}/\Sigma _{\rm th}\right)^{3/2} M_6^{-1/2}$. Inserting Equation (A6) and recalling the pressure equilibrium condition, ρ_g ≈ ρ_cl (Krumholz & McKee 2005), we finally get

$$\begin{equation} \left(\frac{\Sigma _{\rm g}}{\Sigma _{\rm th}}\right) = \rho _1^{1/2} \left(\frac{c_{\rm s}}{8.3 \ \rm km\ s^{-1}}\right). \end{equation} \tag{ A7 }$$

Thus, the overall timescale expression becomes

$$\begin{equation} t_{\rm SF}^{-1} = (\rm 0.8 \ Gyr)^{-1} \times \left\lbrace \begin{array}{@{}l l}\rho _1^{0.165} \left(\frac{c_{\rm s}}{8.3\ \rm km\ s^{-1}}\right)^{-0.99} & \Sigma _{\rm g} < \Sigma _{\rm th},\\ \rho _1^{1/2} \left(\frac{c_{\rm s}}{8.3\ \rm km\ s^{-1}}\right)^{-0.32} & \Sigma _{\rm g} > \Sigma _{\rm th}. \end{array} \right. \end{equation} \tag{ A8 }$$

We now turn our attention to the fraction f_c of mass in cold gas available for star formation (Equation (A4)). Unlike t_SF, f_c depends on the gas density and not on the GMC density. In more tenuous regions, out of pressure equilibrium, the cloud can be denser than the background gas, i.e., η ≡ ρ_cl/ρ_g ⩾ 1. Therefore, in principle, in Equation (A4) we should use a modified version of Equation (A7), e.g., an expression $\Sigma _{\rm g} \propto \eta ^{-1/3} \rho _{\rm g}^{1/2} c_{\rm s}$, which could be used both above (η ≈ 1) and below threshold. However, since f_c has a floor of 2% (see Equation (A4)), which kicks in when η takes values of a few to several, the modification factor (η^−1/3 ≈ 1–2) is negligible relative to the other uncertainties in the derivation of Equation (A7). We will therefore always directly use Equation (A7) in Equation (A4), whether above or below threshold.

A.2.1. Starburst in Merging Galaxies

Observations suggest a link between starbursts and mergers. The main evidence is that the strongest starbursts (ultraluminous and hyperluminous infrared galaxies) are predominantly merging systems at all redshifts (e.g., Elbaz & Cesarsky 2003). More recently, CO observations of galaxy populations hinted that starburst/merging galaxies have a different, more efficient star formation law (Daddi et al. 2010; Genzel et al. 2010). We therefore assume that in merging systems the star formation is driven by different dynamical processes in the bulge which induce star formation over a dynamical time. In practice, we use the same procedure as in B12: in the gaseous bulge that forms after the merger, star formation is regulated by

$$\begin{equation} \dot{\rho }_{*}= \frac{\rho _{\rm g}}{t_{\rm ff}}, \end{equation} \tag{ A9 }$$

where $t_{\rm ff}=\sqrt{3\pi /(32 G \rho _{\rm g})}$ is the local dynamical time for the gas. Again, this expression can be integrated over the bulge volume to yield the total star formation rate.

Alongside our fiducial models, for the sake of comparison we also investigate three models implementing the coherent and chaotic scenarios often used in the literature.

• 1.  
  Coherent model. Each accretion event takes place in a well-defined plane, persisting for the duration of the episode and efficiently spinning-up the hole (Thorne 1974); see for instance, the coherent model of Berti & Volonteri (2008).
• 2.  
  Chaotic I model. For each accretion event, we always take F = w = 0.5, independently of the accreted mass (see for instance, the chaotic model of Berti & Volonteri 2008)¹⁸.
• 3.  
  Chaotic II model. The model of this paper, but with F = 1/2, i.e., with isotropic distribution for the angular momentum of the gas clouds (King & Pringle 2006).

The overall redshift evolution of the MBH spin distributions for the coherent, chaotic I, and chaotic II models is shown in Figure 17 (to be compared to Figure 15). In the coherent model, MBHs tend to be maximally spinning irrespective of redshift, accretion rate, and galaxy host morphology. The opposite is true for the chaotic I model, where spins tend to be small (a_bh < 0.5) and cluster around zero. In the chaotic II model, there is a transition between maximally spinning MBHs at low masses and nonspinning MBHs at high masses. The difference from the chaotic I model is that each accreting lump of matter has a defined angular momentum (although with random direction); light MBHs thus align with it and are efficiently spun up, whereas massive MBHs do not and effectively experience zero-angular momentum accretion averaged over many episodes.

Zoom In Zoom Out Reset image size

The features in the spin distribution discussed above already hint at the difficulty of reconciling these three models with observations. We do not show here a visual comparison, but we report the results of the same statistical analysis performed in Sections 6.3 and 6.4. In this case we summed up accreting MBH in spirals and in ellipticals when building spin distributions. The procedure does not affect the results because in these simplistic models the spin evolution is not connected to the nature of the host galaxy, and it allows us to increase the size of the theoretical sample. The coherent model shows properties similar to the disk one, but with a stronger clustering toward maximally spinning MBHs, resulting in a poorer match to the data, quantified by a KS probability p_KS < 10⁻². Both the chaotic I and chaotic II models result in very different distributions favoring low spins (as shown in the bottom right panels of Figure 17), which are impossible to reconcile with observations (p_KS < 10⁻⁴ in both cases; see numbers in Table 5). This is confirmed by the computation of the odds ratios against, e.g., the hybrid model, which completely discard all three scenarios.

As explained in the main text, to carry out our analysis we used our model to construct discrete theoretical distributions on a specific grid in the M_bh–a_bh parameter space. We then compared spin measurements at z < 0.1 to theoretical predictions at z = 1. In this appendix, we show that our conclusions do not depend significantly on these specific choices.

As mentioned, we choose to compare observations to theoretical distributions at z = 1 rather than at z = 0.1 in order to have a larger statistical sample. The idea behind this choice is that theoretical distributions should not change significantly at z < 1, where QSO and AGN activity and star formation are fading. This is shown to indeed be the case in the left panel of Figure 18, where we plot the redshift evolution of the population of accreting MBHs in spirals for different mass bins and models. The figure shows that distributions generated with the pseudobulge and hybrid models do not change much at z ⩽ 1, but the effect of the small statistical sample at low redshifts is clear, especially at z = 0.1 (thin blue histograms) for the pseudobulge model. Distributions generated with the disk model are slightly more redshift-dependent. At z = 0.1–0.5 (thin blue and medium green histograms) they show more clustering at almost maximal spins, and the moderate spin tails present at z = 1 (thick red histograms) are suppressed. (As we will show below, this has an impact when estimating the model's compatibility with the data).

Zoom In Zoom Out Reset image size

For the analysis of the main text, we evaluated theoretical distribution on a grid made of five equally log-spaced mass bins and 20 linear bins covering the spin range 0 < a_bh < 1. We construct here distributions using instead 10 and 30 linear a_bh bins (by keeping 5 mass bins) and also considering 20 bins both in mass and spin. In the right panel of Figure 18, we show the dependence of the distributions on the a_bh binning. The distributions for the disk and hybrid models both depend on the binning at the high-spin end. This is not surprising because the two models adopt similar accretion prescriptions for MBHs in spirals that recently experienced galactic-disk instabilities. However, as we show below, in the hybrid model this binning dependence does not have a large impact on the comparison with the data. This is because a significant tail at moderate spins is guaranteed by the fact that the model is also partly pseudobulge in nature (see the model description in Section 3.2).

Being aware of these issues, we assess here the impact of the theoretical sample's binning and redshift on the results of our statistical analysis. Results for the 2D KS tests are reported in Table 5. The hybrid and pseudobulge models yield probabilities that are largely independent on the choice of binning and redshift. We just notice a significant discrepancy at z = 0.1 in the subsample of accreting MBHs in spirals in the pseudobulge case. This is entirely due to the small number of objects found in our simulations that cause a very noisy distribution. The chaotic I and chaotic II models are also unaffected by the specific redshift and binning choice and are always ruled out at high significance. Conversely, the match between the disk and coherent models and observations is binning-dependent. Although apparently problematic, this dependence has a physical origin. In our thin accretion disk model, the maximum MBH spin is a_bh = 0.998 (Thorne 1974). Basically, all the simulated MBHs falling in the highest a_bh bin indeed have this spin value. By binning the spin distribution, we are "spreading" those systems on the width of the bin; the larger the bin, the larger the spread, making it easier to reconcile the theoretical distribution with the large number of spins measured in the [0.85,1] range (see Table 2). If the maximum MBH spin is indeed a_bh = 0.998, then the finest spin binning should provide the most trustworthy results. We caution, however, that such a narrow spin distribution, peaked close to the maximal spin a_bh = 1, is affected by some simplifications that we made in our spin evolution model. As described in Sections 2.2 and 3.3, we assumed (1) that at the beginning of each accretion episode the MBH spin is parallel to the total angular momentum of the reservoir and (2) that whenever J_disk/2J_bh > 1, the alignment process is very fast, and accretion is effectively coherent. Relaxing these two assumptions would affect the highly spinning MBHs, because when a_bh ≈ 1, a small amount of retrograde accretion can efficiently spin the black hole down. The effect of relaxing assumption (1) has already been discussed in Section 3.3, so here we will simply estimate the effect of relaxing prescription (2). For a broad range of MBH masses, J_disk/2J_bh can be larger than, but still close to, 1. In these cases, if the accretion event is initially misaligned with respect to $\boldsymbol{J}_{\rm bh}$ by more than π/2, accretion would be retrograde during the first part of the spin realignment process. If a fraction of 10% of the gas mass were accreted on retrograde orbits, the equilibrium spin would be the same as in the case F ≈ 0.9 when neglecting the spin realignment, i.e., a_eq ≈ 0.9 (see Figure 3). Such an effect has been observed and discussed in studies where the spin direction was evolved during each single-accretion event (Dotti et al. 2010; D13). Furthermore, different accretion disk models allow MBH spins only up to a_bh = 0.9–0.95 (e.g., Gammie et al. 2004), and efficient angular momentum extraction via jets might set a maximum equilibrium spin lower than a_bh = 0.9 (see, e.g., Moderski & Sikora 1996; Moderski et al. 1998). It is therefore impossible to completely rule out the disk and coherent models.

We also checked the outcome of the Bayesian model comparison for the different assumptions; results are shown in Table 6. As expected from the 2D KS tests, the consistency of the disk model with observations is highly dependent on the binning size of the distribution. This is confirmed by our model selection exercise, which shows that both log Λ_pd and log Λ_hd generally increase with the number of a_bh bins. Note that for any choice of binning and sample redshift, log Λ_hd > 2, implying that the hybrid model always provides a better description of the data than the disk model, at least at the 99% confidence level. This is particularly interesting, because even though the disk model computed on 10 a_bh bins at z = 1 passes the 2D KS test with flying colors (see Table 5), the odds ratio test provides compelling evidence in favor of the hybrid model. The pseudobulge model sits somewhat in the middle; it is generally preferred at the 90%–99% level over the disk model, but is disfavored at about the same level with respect to the hybrid model (with the anomaly of the z = 0.5, 30 a_bh bins). To summarize, the odds ratio analysis always provides moderate (decisive) evidence that the hybrid model is a better description of the data than the pseudobulge (disk) model, thus proving that our main results are independent of the particular choice of binning and redshift of the theoretical samples.