GALAXY FORMATION WITH SELF-CONSISTENTLY MODELED STARS AND MASSIVE BLACK HOLES. I. FEEDBACK-REGULATED STAR FORMATION AND BLACK HOLE GROWTH

The ZEUS astrophysical hydrodynamics module included in Enzo is employed to solve the Euler equations for the collisional baryon fluid represented by grids (Stone & Norman 1992a, 1992b; Anninos & Norman 1994). While known to introduce spurious effects, this scheme is widely used with AMR because of the stability of its solutions and the acceptable error when combined with high resolution.

Dark matter, stars, and MBHs are treated as collisionless particles which interact only by the gravitational force. To evolve the particle positions and velocities, the gravitational dynamics are solved by an N-body adaptive particle-mesh solver. After particles are gridded onto the mesh by the cloud-in-cell interpolation, the Poisson equation is solved on the discretized density grids via fast Fourier transform and multigrid solvers (Hockney & Eastwood 1988; O'Shea et al. 2004).

Enzo decides whether each parent cell needs to be refined into eight child cells based on the mass of the cell in gas or in particles. The time step is also adaptively determined level by level so that the time step dt satisfies

$$\begin{eqnarray} dt &\le 0.3 \times { \Delta x \over c_{\rm s} } = 0.3 \times {\rm (sound\,\,crossing \,\, time)} \end{eqnarray} \tag{ 2 }$$

for all the cells at that level. Here c_s is the sound speed of the gas, and we choose the Courant–Friedrichs–Lewy safety number of 0.3. As decisions for refinement are made recursively, the resulting data set is a nested grid-patch structure. In our work, the grids are adaptively refined down to 15.2 pc resolution. This value is in accord with the Jeans length for a dense gas clump of n = 125 cm⁻³ at ∼200 K, at which point a corresponding Jeans mass of 16, 000 M_☉ collapses to spawn a molecular cloud particle.

We refine the cells by factors of two in each axis, on gas and particle overdensities of eight. The mass thresholds, M_ref, above which a cell refines are functions of a refinement level l as

$$\begin{eqnarray} M_{\rm ref, gas}^{l} & = & 2^{-0.378 l} M_{\rm ref, gas}^0 = 2^{-0.378 l} \cdot 0.125\Omega _{\rm b} \rho _0 \Delta x^3\quad\ \ \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} M_{\rm ref, part}^{l} & = & 2^{-0.105 l} M_{\rm ref, part}^0 = 2^{-0.105 l} \cdot 0.125 \Omega _{\rm m} \rho _0 \Delta x^3,\quad \end{eqnarray} \tag{ 4 }$$

where factors 0.125 = 8(1/2³)² guarantees to refine all the cells of the first two nested levels (Section 3.1). Δx is the cell size at a root grid and ρ₀ = 3H²/8πG is the critical density. Ω_m = 0.27, Ω_b = 0.044, and H = 71  km  s⁻¹  Mpc⁻¹ are matter density, baryon density, and the Hubble constant, respectively.

For example, at the finest static level, l = 2, a cell is refined if it has more mass than 8.9 × 10⁵ M_☉ in gas or 6.7 × 10⁶ M_☉ in particles. At level l = 11 (Δx = 15.2  pc at z = 3) a cell is refined if more than 8.4 × 10⁴ M_☉ = 5  M_Jeans(125  cm⁻³, 200  K) in gas or 3.5 × 10⁶ M_☉ = 47  M_{DM, smallest} in particles. This way we refine the grids more on small scales, which allows us to focus our computational resources more on the dense star-forming regions, making the simulation super-Lagrangian (O'Shea & Norman 2008).

We use non-equilibrium chemistry model to track six species (H, H⁺, He, He⁺, He⁺⁺, ${{\textit {e}}}^-$) by following six collisional processes among them. At the same time, Enzo's cooling module considers collisional excitation cooling, collisional ionization cooling, recombination cooling, Bremsstrahlung cooling, and cosmic microwave background Compton cooling to compute the radiative loss of internal gas energy (Anninos et al. 1997; Abel et al. 1997). Added to these primordial cooling rates is the metallicity-dependent metal cooling rate ΔΛ(Z) = Λ_net(Z) − Λ_net(0) above 10⁴ K, where Λ_net is the net cooling rate tabulated in Sutherland & Dopita (1993). Cooling below 10⁴ K is also enabled with fine structure metal-line cooling by C, O, and Si (Glover & Jappsen 2007; Wise & Abel 2008). This treatment ensures that a thin galactic disk forms by being cooled below 10⁴ K, the approximate virial temperature of the ISM in a galactic disk.

We further refine our module with photoionization heating at z < 3 by the metagalactic background UV of quasars and galaxies (Haardt & Madau 1996, 2001), which is known to give rise to a warm diffuse ISM and prevent star formation in optically thin gas (Ceverino & Klypin 2009). An approximate self-shielding factor is applied when the heating term is added (Cen et al. 2005). While not introducing a marked difference in overall results analyzed here, inclusion of this additional heating term results in a more realistic ISM.

Our molecular cloud particle formation is based on Cen & Ostriker (1992) formalism with several important modifications. With a fixed formation efficiency of _* = 0.5, the finest cell of physical size Δx = 15.2  pc and gas density ρ_gas produces a molecular cloud particle of mass

$$\begin{eqnarray} &&M_{\rm MC} = \epsilon _* \rho _{\rm gas} \Delta x^3 \end{eqnarray} \tag{ 5 }$$

when the following four criteria are met.

• 1.  
  the proton number density exceeds the threshold n_thres = 125 cm⁻³,
• 2.  
  the velocity flow is converging, i.e., $\nabla \cdot \mathbf {v} < 0$,
• 3.  
  the cooling time t_cool is shorter than the dynamical time t_dyn of the cell: $E_{\rm int} / \dot{E} < [3\pi / (32G\rho _{\rm gas})]^{1/2}$, and
• 4.  
  the particle produced has at least M_thres = 8000 M_☉.

The consequence of our criteria is the following. The gas in the finest cell is converted into a particle as soon as the cell has accumulated more than M_thres/_* = 16, 000 M_☉, the Jeans mass at n = 125 cm⁻³ at ∼200 K. Because 8000–42,000 M_☉ is instantly removed from the cell every time a particle is created, the gas mass in the finest cell never reaches the refinement threshold M^{l = 11}_{ref, gas} = 84, 000 M_☉ described in Section 2.2, ensuring the consistency between the refinement criteria and the particle formation.⁹ The values used here are in good agreement with those corresponding to collapsing giant molecular clouds (GMCs; McKee & Ostriker 2007) where star-forming molecular clumps are enshrouded by cold atomic gas.

As an additional note, differences from more traditional star particle formation criteria such as in Tasker & Bryan (2008) include: (1) the Jeans condition ρ_gasΔx³ > M_Jeans is removed because this condition could have allowed mass greater than M_Jeans to accumulate while not being properly resolved until a particle finally forms, (2) the factor Δt/t_dyn in Equation (1) of Tasker & Bryan (2008) is removed in order to instantly create a particle and not leave any unresolved mass behind, and (3) stochastic star formation is not imposed.

With these modifications, our criteria guarantee that a particle forms before an unphysically large mass begins to accrete onto any unresolved dense clump. It is worth to emphasize the differences between our molecular cloud formation criteria and prior studies. While many previous studies with particle-based codes (e.g., Mihos et al. 1991; Katz 1992; Mihos & Hernquist 1994; Springel & Hernquist 2003; Governato et al. 2007) place a star particle using the Schmidt relation (ρ_SFR ∼ ρ^1.5_gas; Schmidt 1959), we deposit a particle when a gas cell of a typical molecular cloud size actually becomes Jeans unstable. For this reason, the particle in our simulation represents a star-forming molecular cloud that is self-gravitating, is thus decoupled from the gas on the grid.¹⁰ It is tagged with its mass M_MC, dynamical time $t_{\rm dyn} = {\tt max} ([3\pi / (32G\rho _{\rm gas})]^{1/2}, 1.0 \,\,{\rm Myr})$, creation time t_cr, and metallicity. Each molecular cloud particle gradually yields an actual stellar mass, M_star(t), over 12 t_dyn which then contributes to stellar feedback (see Figure 1 and Section 2.5).

Zoom In Zoom Out Reset image size

Observational evidence suggests that only ∼2% of the gas in GMCs is converted into stars per dynamical time (Krumholz & McKee 2005; Krumholz & Tan 2007, and references therein). Numerical studies also indicate that turbulence, magnetic fields, or radiation pressure can make the star formation process surprisingly slow (e.g., Murray et al. 2010; Wang et al. 2010). To reflect these observations in our simulation, only 20% of the molecular cloud particle mass, M_MC, turns into an actual stellar mass, M_star(t), over 12 t_dyn by

$$\begin{eqnarray} M_{\rm star}(t) & = & 0.2 \, M_{\rm MC} \int _0^{\tau } \tau \,^{\prime } e^{-\tau \,^{\prime }}\,d\tau \,^{\prime } \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&\qquad\quad\ = 0.2 \, M_{\rm MC} [ 1- (1+ \tau) e^{-\tau } ], \end{eqnarray} \tag{ 7 }$$

where τ = (t − t_cr)/t_dyn. In this formulation, the production of the stellar mass peaks at t_dyn. As 7.5 × 10⁻⁷ of the rest mass energy of M_star is gradually deposited into the cell in which the particle resides,¹¹ this thermal stellar feedback replenishes the energy loss to radiative cooling. At the same time, the rest of the molecular cloud particle mass, 0.8 M_MC, slowly returns to the gas grid. This again reflects the fact that most of the gas in GMCs does not end up locked in stars in a few dynamical time, but is blown out into the ISM to be recycled. Meanwhile, 2% of the ejected mass is counted as metals, contributing to the metal enrichment of the ISM (see Figure 1).

Overall, our feedback treatment corresponds to the energy of 10⁵¹ erg for every 750 M_☉ of actual stellar mass formed. Although Type II supernovae explosions are its dominant source (Spitzer 1990; Tasker & Bryan 2006, 2008), this feedback also models various other types such as protostellar outflows (Li & Nakamura 2006; Nakamura & Li 2008), photoionization (McKee 1989), and stellar winds (Oey et al. 2001). Therefore no explicit time delay is necessary between the formation of a molecular cloud and the start of stellar feedback. This thermal feedback heats the mass of ∼10⁴ M_☉ in a <30 pc cell up to ∼10⁷  K, but a multiphase medium (McKee & Ostriker 1977) is naturally established without using any sub-resolution model. The so-called overcooling problem (Somerville & Primack 1999; Balogh et al. 2001) is absent in our simulation since the cooling time of these hot cells is much longer than the sound crossing time (Kim et al. 2009).

A 10⁵ M_☉ MBH is put as a seed at the center of each simulated galaxy. It is treated as a collisionless sink particle, but grows in mass by accreting gas from its surroundings. We estimate the rate of accretion by employing the Eddington-limited spherical Bondi–Hoyle formula (Bondi & Hoyle 1944; Bondi 1952):

$$\begin{eqnarray} &&\ \ \ \dot{M}_{\rm BH} = {\tt min} (\dot{M}_{\rm B} \,,\, \dot{M}_{\rm Edd})\qquad\qquad\qquad\quad \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&\qquad\quad = {\tt min} \left({4 \pi G^2 M^2_{\rm BH} \rho _{\rm B} \over { c^3_{\rm s} } } \,\,,\,\, {4 \pi {\it G M}_{\rm BH} m_{\rm p} \over { \epsilon _{\rm r} \sigma _{\rm T} c } } \right), \end{eqnarray} \tag{ 9 }$$

where M_BH is the mass of an MBH, c_s is the sound speed of the gas at the cell the MBH resides in, m_p is the mass of a proton, and σ_T is the Thomson scattering cross-section. Note that, when compared with Equation (1) the nondimensional parameter α is absent. ρ_B is the density at the Bondi radius

$$\begin{eqnarray} &&R_{\,\rm B} = { 2 {\it G M}_{\rm BH} \over c^2_{\rm s}} \, \simeq \, 8.6 \,\,{\rm pc} \left({ M_{\rm BH} \over 10^5\, M_{\odot }} \right) \left({ 10 \,{ \rm km\, s^{-1}} \over c_{\rm s} }\right)^{2},\quad \end{eqnarray} \tag{ 10 }$$

and is extrapolated from the density ρ_gas of the cell of size Δx where the MBH resides by

$$\begin{eqnarray} &&\rho _{\rm B} = \rho _{\rm gas} \,\cdot \, {\tt min} ((\Delta x/R_{\rm B})^{1.5}, 1.0) < \rho _{\rm gas}. \end{eqnarray} \tag{ 11 }$$

Here an r^−3/2 density profile is assumed inside the sphere of R_B (Wang et al. 2010). Adopting a radiative efficiency _r = 0.1 for a non-rotating Schwarzschild black hole (Shakura & Sunyaev 1973; Booth & Schaye 2009), the Eddington rate for a 10⁵ M_☉ black hole is ≃ 0.002 M_☉ yr⁻¹. For more discussion on the accretion rate adopted here, see Appendix A. To minimize any numerical artifacts, the gas mass accreting onto the MBH is uniformly subtracted from grid cells within a Bondi radius. The MBH also inherits the momentum of the accreting gas.

Most importantly, to probe the gas dynamics accreting onto the MBH and to fully incorporate the MBH in a galactic simulation, it is imperative to always reach the resolution close to the Bondi radius around the MBH. To resolve the gas around the MBH with the best resolution available, eight nearby cells close to the MBH are required to successively refine down to 15.2 pc (proper) at all times. In practice, the MBH naturally sits at the densest region most of the time, surrounded by many finest cells. While our spatial resolution is still slightly too large to resolve the Bondi radius of a 10⁵ M_☉ black hole, Equation (10), it is enough to resolve the Bondi radii of more massive MBHs such as in nearby X-ray luminous galaxies (e.g., ∼120 pc for SMBH in M87; Allen et al. 2006). This shows that our simulations are beginning to depict the self-consistent coevolution of both galaxies and MBHs in one comprehensive framework. Admittedly, this resolution is still far from the Schwarzschild radius of any black hole

$$\begin{eqnarray} &&R_{\,\rm Sch} = { 2 {\it G M}_{\rm BH} \over c^2} \, \simeq \, 10^{-8} \,{\rm pc} \left({ M_{\rm BH} \over 10^5\, M_{\odot }} \right), \end{eqnarray} \tag{ 12 }$$

which is needed to thoroughly describe its accretion disk. Due to our resolution limit, an MBH particle in our framework represents not just the black hole itself, but also includes accreting gas and stars deep within the galactic nucleus; in other words, the Bondi–Hoyle accretion estimate does not accurately model the physics below the resolution limit (see Section 5.2).

We now turn our attention to the feedback of an accreting MBH. The gravitational potential energy of the gas accreting onto a black hole is extracted during the gravitational infall. Assuming an infall down to the innermost stable orbit of an accretion disk, the conversion rate from the rest mass energy to feedback energy is 10%, previously defined as the radiative efficiency _r. Hence the bolometric radiation luminosity of an MBH is

$$\begin{eqnarray} &&L_{\rm BH} = \epsilon _{\rm r} \dot{M}_{\rm BH} c^2. \end{eqnarray} \tag{ 13 }$$

As was discussed earlier, for a long time a thermal energy deposition has been the dominant strategy to treat the feedback of an accreting MBH (Springel et al. 2005b; Di Matteo et al. 2005; Colberg & Di Matteo 2008; Booth & Schaye 2009; Callegari et al. 2010; Teyssier et al. 2010, see Sijacki et al. (2008) or Debuhr et al. (2011) for other approaches). Without question, it has been an effective approximation characterizing the impact of an accreting MBH on a resolved scale when sufficient resolution or full radiative transfer is inaccessible (see Figure 2(c); Springel et al. 2005b). Despite its practical efficiency, however, better feedback models are imperative for high-resolution galaxy formation studies where the Bondi radius is starting to be resolved. In the next two sections, we explain the detailed implementations of two modes of MBH feedback: radiative and mechanical. The thermal feedback model previously used can be regarded as an approximation of these two feedback channels combined.

Although the radiation from the MBH in a galaxy was tested in spherically symmetric or axisymmetric models (Ciotti & Ostriker 2007; Proga et al. 2008; Ciotti et al. 2009; Kurosawa et al. 2009; Park & Ricotti 2010), a three-dimensional radiative transfer calculation of the impact of an MBH has never been performed in galactic scale simulations. In what follows, we treat the MBH as a point source of radiation and carry out a three-dimensional transport computation to evolve the radiation fields (see Figure 2(a); note that molecular cloud particles are not treated as radiation sources). Achieving high resolution around the MBH is critical here because, if otherwise, the optical depth of the radiation could be small even at the smallest resolved distance from the MBH (Omma et al. 2004).

Zoom In Zoom Out Reset image size

Enzo's radiative transfer module incorporates the adaptive ray tracing technique (Abel & Wandelt 2002) with the hydrodynamics, energy, and chemistry solvers. It has been applied to problems such as the radiative feedback from Population III stars (Abel et al. 2007; Wise & Abel 2008; Wise et al. 2010) and from Population III black holes (Alvarez et al. 2009). For algorithmic and numerical details of Enzo radiative transfer we refer the readers to Wise & Abel (2011); and, here we briefly describe the machinery relevant to the presented results. First the luminosity of the MBH is assigned by Equation (13). Then 768 (=12 × 4³; Healpix level 3) rays are isotropically cast with a monochromatic energy of E_ph = 2  keV, a characteristic temperature of an averaged quasar spectral energy distribution (SED; Sazonov et al. 2004, 2005; Ciotti & Ostriker 2007).¹² Consequently the number of photons per each initial ray is

$$\begin{eqnarray} &&P_{\,\,\rm init} = { L_{\rm BH}\, dt_{\rm ph} \over E_{\rm ph} \cdot 768} = { \epsilon _{\rm r} (\dot{M}_{\rm BH}dt_{\rm ph}) c^2 \, \over E_{\rm ph} \cdot 768} \end{eqnarray} \tag{ 14 }$$

given the photon time step dt_ph which we set as the light-crossing time of the entire computational domain. This choice is justified because the photons are in a free streaming regime, and the energy deposited by the radiation per time step is relatively small. Each ray is traced at speed c until the ray reaches the edge of the computational domain or most of its photons (99.99995%) are absorbed. It is adaptively split into four child rays whenever the area associated with a ray becomes larger than 0.2 (Δx)² of a local cell.

Photons in the emitted ray then interact with the surrounding gas in three ways: they (1) ionize the gas, (2) heat the gas, and (3) exert momentum onto the gas. First, the ray loses its photons when it photoionizes H, He, and He⁺ with the respective photoionization rates of

$$\begin{eqnarray} k_{\rm ph, H} & = & { P_{\,\,\rm in} (1-e^{-\tau _{\rm H}}) (E_{\rm ph}Y_{k, \rm H} /E_{i, \rm H}) \over n_{\rm H} (\Delta x)^3 dt_{\rm ph} } \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&\ \ \ k_{\rm ph, He} = { P_{\,\,\rm in} (1-e^{-\tau _{\rm He}}) (E_{\rm ph}Y_{k, \rm He} /E_{i, \rm He}) \over n_{\rm He} (\Delta x)^3 dt_{\rm ph} } \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&k_{\rm ph, He^+} = { P_{\,\,\rm in} (1-e^{-\tau _{\rm He^+}}) \over n_{\rm He^+} (\Delta x)^3 dt_{\rm ph} },\qquad\qquad\ \ \end{eqnarray} \tag{ 17 }$$

where P_in is the number of photons coming into the cell, τ_H = n_Hσ_Hdl is the optical depth, n_H is the hydrogen number density, σ_H is the energy-dependent hydrogen photoionization cross-section (Verner et al. 1996), dl is the path length through the cell, and E_i = 13.6,  24.6,  54.4 eV are the ionization thresholds for H, He, He⁺, respectively. The factors Y_k are the energy fractions used for ionization when secondary ionizations are considered (Shull & van Steenberg 1985).¹³

Second, the excess energy above the ionization threshold, E_i, heats each of the species with the photoheating rates of

$$\begin{eqnarray} &&\Gamma _{\rm H} = { P_{\,\,\rm in} (1-e^{-\tau _{\rm H}}) E_{\rm ph} Y_{\Gamma } \over n_{\rm H} (\Delta x)^3 dt_{\rm ph} }, \,\,\,\,{\rm etc.}, \end{eqnarray} \tag{ 18 }$$

where Y_Γ is the fraction of energy deposited as heat when secondary ionizations are taken into account.¹⁴ The 2 keV soft X-ray photon can also scatter off and heat an electron resulting in the Compton heating rate of

$$\begin{eqnarray} &&\Gamma _{\rm C} = { P_{\,\,\rm in} (1-e^{-\tau _{\rm e}}) \Delta E(T_{\rm e}) \over n_{\rm e} (\Delta x)^3 dt_{\rm ph} }, \end{eqnarray} \tag{ 19 }$$

where τ_e = n_eσ_KNdl is the optical depth, n_e is the electron number density, σ_KN is the nonrelativistic Klein–Nishina cross-section (≃ σ_T; Rybicki & Lightman 1979), and ΔE(T_e) = 4k_BT_e · (E_ph/m_ec²) is the nonrelativistically transferred energy to an electron at T_e (Ciotti & Ostriker 2001). It should be noted that, in Compton scattering, a photon loses its energy by a factor of ΔE(T_e)/E_ph, but essentially keeps propagating without being absorbed. However, in order to model this with monochromatic photons, we instead subtract $P_{\,\,\rm in} (1-e^{-\tau _{\rm e}}) \Delta E(T_{\rm e})/E_{\rm ph}$ photons from the ray. This is another way a ray loses its photons while traveling through a cell. Combined, the total heating rate by absorbed and scattered photons becomes

$$\begin{eqnarray} &&\Upsilon = n_{\rm H} \Gamma _{\rm H} + n_{\rm He} \Gamma _{\rm He} + n_{\rm He^+} \Gamma _{\rm He^+} + n_{\rm e} \Gamma _{\rm C}\,. \end{eqnarray} \tag{ 20 }$$

Lastly, photons exert outward momentum to the gas when they are taken out from the ray either by photoionization or by Compton scattering. It was claimed that the radiation pressure from the MBH may markedly alter the environment near the MBH, especially within ∼0.1 kpc in radius (Haehnelt 1995; Debuhr et al. 2011). The large-scale galactic wind driven by deposited photon momentum is also considered as a possible explanation for the M_BH–σ_bulge relation (Murray et al. 2005). The added acceleration onto the cell by the radiation pressure is calculated by

$$\begin{eqnarray} &&{\bf a}_{\rm ph} = {d{\bf p}_{\rm ph} \over m_{\rm cell} dt_{\rm ph}} = { P_{\,\,\rm lost} E_{\rm ph} \over \rho _{\rm gas} (\Delta x)^3 c dt_{\rm ph} } {\hat{\bf r},} \end{eqnarray} \tag{ 21 }$$

where dp_ph is the photon momentum exerted onto the cell in dt_ph, P_lost is the number of photons lost in the cell, and ${\hat{\bf r}}$ is the directional unit vector of the ray. Neglecting the radiation pressure on dust grains is conservative because its inclusion would further enhance the negative feedback effect (see Section 5.2).

Observations find that a significant portion of the energy extracted during the accretion onto an MBH is released as mechanical energy, creating bipolar jets (Bridle & Perley 1984; Pounds et al. 2003) or inflating cavities (Fabian et al. 2002; McNamara et al. 2005) at the sites of active galactic nuclei (AGNs). A number of authors have used a numerical approach to explore the effectiveness of jets in heating up a cooling flow (Fabian et al. 1994; Peterson & Fabian 2006); most of them targeted the gas dynamics in galaxy clusters with ∼kpc resolution excluding detailed galactic scale physics (e.g., Omma et al. 2004; Cattaneo & Teyssier 2007; Antonuccio-Delogu & Silk 2008; Dubois et al. 2010). In the meantime, a numerical analysis on stellar winds from nuclear disk or MBH jets has been carried out in a galactic scale, but only in an one-dimensional context (Ciotti et al. 2009, 2010; Shin et al. 2010) Here, we construct a mechanical feedback model of an MBH applicable in three-dimensional galactic simulations, which creates accretion-rate-dependent subrelativistic bipolar jets launched at the vicinity of the MBH (see Figure 2(b)).

Let us assume that all of the bolometric luminosity of the MBH, L_BH, is converted to the "mechanical" power of jets. Because the ejecta has to climb out of the potential well of the MBH, the "kinetic" power of the jets is less than L_BH by

$$\begin{eqnarray} L_{\rm BH} & = & P_{\rm mech}\qquad\qquad\qquad\qquad\qquad \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} &&\qquad = P_{\rm kin} + {\rm (gravitational \,\, potential \,\, energy)}. \end{eqnarray} \tag{ 23 }$$

Therefore the "kinetic" power of the jets, as we introduce at a scale of R_jet = 2Δx = 30.4  pc, can be written as

$$\begin{eqnarray} &&P_{\rm kin} = \epsilon _{\rm kin} L_{\rm BH} = \epsilon _{\rm kin} \epsilon _{\rm r} \dot{M}_{\rm BH} c^2 = {1\over 2} \dot{M}_{\rm jet} v_{\rm jet}^2, \end{eqnarray} \tag{ 24 }$$

where _kin < 1 is the "kinetic" coupling constant denoting the fractional energy available for the kinetic motion of the jets (see Figure 2(b)). $\dot{M}_{\rm jet}$ is the mass ejection rate of the jets, and v_jet is the jets velocity when introduced in the simulation. Hence _kin encapsulates not only the acceleration of the jets powered by the AGN central engine, but also the gravitational "redshift" from the scale of an accretion disk (∼R_Sch) to a resolved scale of jets in simulations (∼R_jet). Ciotti et al. (2009) provides estimates for an MBH of $l = \dot{M}_{\rm BH} / \dot{M}_{\rm Edd} = 0.005$ as

$$\begin{eqnarray} \epsilon _{\rm kin} & = & {P_{\rm kin} \over \epsilon _{\rm r} \dot{M}_{\rm BH} c^2} = {0.0125 \over \epsilon _{\rm r}(1+400l)^4} \simeq 0.0015 \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \eta _{\rm jet} &\equiv & {\dot{M}_{\rm jet} \over \dot{M}_{\rm BH}} = {0.2 \over (1+100l)^4} \simeq 0.04,\qquad \end{eqnarray} \tag{ 26 }$$

based on which we fiducially adopt conservative values of _kin = 10⁻⁴ and η_jet = 0.05.

With _kin and η_jet now fixed, the kinetic motion of the jets can be fully described. First, as usual, out of a sphere of R_B centered on the MBH the accreting mass is taken out at every finest hydrodynamical time step dt; then 5% of the accreted mass, $\dot{M}_{\rm jet}dt = 0.05\,\,\dot{M}_{\rm BH} dt$, is set aside as a mass of jets. Now Equation (24) yields the initial jet momentum, $(\dot{M}_{\rm jet}dt) v_{\rm jet}$, with

$$\begin{eqnarray} &&v_{\rm jet} = c \left({ 2\epsilon _{\rm kin} \epsilon _{\rm r} \over \eta _{\rm jet} } \right)^{1/2} = 6000\,\, {\rm km\, s^{-1}} \end{eqnarray} \tag{ 27 }$$

for _r = 0.1. This value of v_jet is consistent with numerous observational evidence (e.g., Biretta & Junor 1995; Junor et al. 1999; Homan et al. 2009) and relativistic MHD simulations (e.g., Vlahakis & Königl 2003, 2004) suggesting the existence of at least mildly relativistic AGN jets on scales of 1–10 pc from a central engine. This v_jet is also well matched with the velocity of momentum-driven AGN winds discussed by King (2010). Finally the launch speed of the surrounding cells is found by averaging the momentum of jets and the preexisting gas in those cells (Wang et al. 2010).

One may want to continuously launch the jets at every finest hydrodynamical time step. However, if the injected mass of jets is minuscule compared to the preexisting mass in surrounding cells, the jets make little or no dynamical impact on the surrounding cells after being mass-weighted averaged with them. Since it is unfeasible to resolve all gas cells around the MBH down to $\dot{M}_{\rm jet} dt$, an alternative approach is indispensable. Moreover, there is growing observational evidence of double-lobed radio galaxies (or double–double radio galaxies) implying that the jets have launched in an episodic fashion with jets interruption timescales of 10⁵–10⁸ years (Stawarz 2004; Saikia et al. 2006). These two considerations lead us to adopt the following method: every time the accumulated jet mass, $\Sigma \,\dot{M}_{\rm jet} dt$, exceeds the threshold of 300 M_☉ it is injected in collimated bipolar jets of a width of five finest cells in the vicinity of the MBH. This approach renders jets intermittent (once every 30 Myr if $\dot{M}_{\rm BH} = 10^{-5}\, M_{\odot }\, {\rm yr}^{-1}$) and dynamically important in our calculation.

The jets are injected parallel and anti-parallel to the total angular momentum $\mathbf {L}$ of the accreted gas up to that point. The angular momentum vector $\mathbf {L}$ changes its direction frequently while it asymptotes to the overall galactic rotation axis. This implementation is motivated by the observations of X-shaped radio galaxies where the radio jets rapidly reorient themselves by the interaction with the surrounding gas or by mergers (Merritt & Ekers 2002; Gopal-Krishna et al. 2003). Lastly, since the mechanical or the radiative feedback alone may not describe the whole picture, we include hybrid models in which each of these two channels constitutes half of the MBH bolometric luminosity, L_BH (Sim-RMF; see Table 1).

Note that the mechanical channel has not been a main driver of MBH feedback in the presented calculation because, with highly suppressed mass accretion rate, jets have launched only a few tens of times in 350 Myr (see Section 4.2). We later comment upon its efficiency in Section 5.2.

A three-dimensional cubic volume of 16 comoving Mpc on a side is set up at z = 60 assuming a flat ΛCDM cosmology with dark energy density Ω_Λ = 0.73, matter density Ω_m = 0.27, baryon density Ω_b = 0.044, and Hubble constant h = 0.71 (in the unit of H₀ = 100  km  s⁻¹  Mpc⁻¹). A scale-invariant primordial power spectrum (spectral index n = 1, Eisenstein & Hu 1999) is adopted with σ₈ = 0.81, the rms density fluctuation amplitude in the sphere of 8 h⁻¹  Mpc.

We identify a dark matter halo of ∼10¹² M_☉ at z = 3 by performing a coarse-resolution adiabatic run. Then we recenter the density field around this halo and set up a new initial condition which preserves the same large-scale power yet contains a small-scale power as well, with a 128³ root grid and a series of two nested child grids of twice finer resolution each (160³ cells for level l = 1 and 200³ for l = 2). Therefore the finest nested grid at level l = 2 spans 6.25 comoving Mpc on a side, contains 200³ dark matter particles of 9.6 × 10⁵ M_☉, and manifests the equivalent resolution of a 512³ unigrid. Initially all the cells throughout l = 2 grids are allowed to be further refined; however, the volume in which additional refinement is enabled ($\mathbf {V}_{\rm ref}$ ; in the shape of a rectangular solid) continually shrinks in size in such a way that it encloses only the smallest dark matter particles.¹⁵  An initial metallicity of Z = 0.003 Z_☉ is also set up everywhere to track the metallicity evolution and to facilitate cooling below 10⁴ K.

Our initial condition is first evolved to z = 3 with a low-resolution (121.6 pc) refinement strategy and a particle formation and feedback recipe, without an accreting MBH. At z = 3 we split each dark matter and star particle inside the focused volume ($\mathbf {V}_{\rm foc}$; a rectangular hexahedron of 1.28 comoving Mpc on a side, a subset of $\mathbf {V}_{\rm ref}$) into 13 child particles using the particle refinement technique by Kitsionas & Whitworth (2002). This algorithm places child particles on a hexagonal close packed array and has been applied to many particle-based applications requiring enhanced particle resolution in a resimulated region (e.g., Bromm & Loeb 2003; Kitsionas & Whitworth 2007; Yoshida et al. 2008). After the particle splitting procedure, each dark matter particle in $\mathbf {V}_{\rm foc}$ represents a collective mass of 74,000 M_☉. Across $\mathbf {V}_{\rm foc}$ cells are now allowed to refine up to 11 additional levels, achieving maximum spatial resolution of 15.2 pc at z ∼ 3 (see Section 2.2).

Consequently, this process produces our focused object at z = 3 dubbed a model galaxy, on which a suite of high-resolution simulations is performed (Figure 3). The model galaxy has a mass of M_vir ≃ 9.2 × 10¹¹ M_☉ at z = 3 and a corresponding virial radius of

$$\begin{eqnarray} &&R_{\rm vir} = M_{\rm vir}^{1/3} \left[ {H^2_0h^2 \Omega _{\rm m} \Delta _{\rm c} \over 2 G} \right]^{-{1/3}} \simeq 310 \,\,\,{\rm comoving \,\,\, kpc}\quad\ \end{eqnarray} \tag{ 28 }$$

given h = 0.71, Ω_m = 0.27, and Δ_c = 200. The dark matter halo represented by ∼1.1 × 10⁷ particles constitutes ∼88% of the total mass. About ∼1.0 × 10⁷ particles contain 8.0 × 10¹⁰ M_☉ of stellar mass, whereas the rest, 3.5 × 10¹⁰ M_☉, is in gaseous form available for future star formation, either in the ISM or in the embedding halo. There is no shortage of gas supply, as the gas from outside the halo continuously falls inward either by spherical accretion or by cold accretion along one of the multiple filaments (Dekel et al. 2009; Ceverino et al. 2010). The halo has spin parameters of 0.051 for dark matter and 0.069 for gas. At the center of all lies a 10⁵ M_☉ MBH we plant as a gravitational seed. This choice of the initial MBH mass lies below the Magorrian et al. (1998) relationship assuming 10% of the stellar mass is in the bulge, which may have resulted in a weaker mode of MBH feedback—possibly a "radio-mode" analog—and the negligible gas expulsion by the MBH (see Sections 4.2 and 4.4). Therefore the reported results should not be interpreted as a general picture of MBH feedback. It remains to be seen whether more massive MBHs or fast growing MBHs have different effects. We will come back to this issue in Section 4.4.

Zoom In Zoom Out Reset image size

First we check the validity of our molecular cloud formation criteria (Section 2.4) and stellar feedback (Section 2.5) by comparing star formation rate (SFR) with gas density. Figure 4 displays a relation between the SFR surface density and the gas surface density in Sim-SF at z = 2.75. In a (20  kpc)³ box centered on an MBH, each data point is made by taking the mean values in a (1  kpc)³ bin, the typical aperture size of Σ_SFR − Σ_gas studies for spatially resolved nearby galaxies (e.g., Kennicutt et al. 2007). Here Equation (7) is used to calculate the stellar mass newly spawning in each cell, and the data points below the observation limit, 10⁻⁵ M_☉ yr⁻¹ kpc⁻², are discarded.

Zoom In Zoom Out Reset image size

The molecular cloud formation and stellar feedback, joined with high spatial resolution, work together to self-regulate star formation. However, authors note that our best fit to this particular snapshot of the galaxy is steeper than the observed trend of z ∼ 0 with larger dispersion.

Now we turn to the topic of an accreting MBH and its feedback. We focus on how MBH feedback changes its surrounding ISM and how it locally suppresses molecular cloud formation. For this purpose, we hereafter examine the snapshot of the model galaxy at z = 2.75 (or 2410 Myr after the big bang), about 220 Myr after an accreting MBH is placed at the center of the galaxy. The model galaxy now has a mass of M_vir ≃ 9.0 × 10¹¹ M_☉ and correspondingly, a virial radius of R_vir ≃ 80  kpc (proper; radii are hereafter in proper kpc, not comoving, unless marked otherwise).

When an MBH starts to accrete gas, the gravitational potential energy of the accreting gas is released in the form of radiation and jets. Even in the case of a slowly growing MBH, as in our simulations (a possible "radio-mode" analog; $\dot{M}_{\rm BH} \sim 0.05 \,\dot{M}_{\rm Edd}$; see Section 4.4), the feedback from the MBH is known to play a major role in regulating star formation and its own growth (Croton et al. 2006; McNamara & Nulsen 2007; Sijacki et al. 2007; Puchwein et al. 2008).

The density and temperature structures in the central regions of the galaxies from Sim-SF (left; without MBH feedback) and Sim-RMF (right; with radiative and mechanical MBH feedback) are shown in Figures 5 and 6. In particular, in Figure 6 for Sim-RMF, a hot region of size ∼2 kpc surrounding the MBH is prominent, which remarkably contrasts with a much colder ISM in Sim-SF. This region is heated up to 10⁶ K by ionizing photons heating hydrogen and helium, and scattering off electrons. The latter, i.e., Compton heating, is important especially in a highly ionized region. When the fractions of neutral species (or singly ionized He) are low, the effect of photoheating on H, He, and He⁺ is mild; instead, the contribution of Compton heating on electrons is relatively large. The hot temperature at the center of Sim-RMF is also evident in the radially averaged temperature profile of Figure 7. Note that the 1–2 kpc distance over which the gas is heated is consistent with the characteristic distance found in an analytic study (Figure 4 of Sazonov et al. 2005) out to which the gas is heated by photoionization and Compton scattering.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A hot temperature in the inner core of the galaxy in Sim-RMF leads to a significant deprivation of cold, dense star-forming gas. Figure 8 illustrates how the structure of ISM is changed by MBH feedback, in terms of joint probability distribution functions of gas density and temperature.

• 1.  
  The left figure depicts a typical ISM without MBH feedback but still with stellar feedback (Sim-SF). It features a multiphase ISM that is naturally achieved in adaptively refined mesh, including cold, dense star-forming gas (T <10⁴ K, ρ > 10⁻²⁴  g  cm⁻³), and hot diffuse supernovae bubbles. As expected, above the molecular cloud formation threshold (n_thres = 125 cm⁻³; denoted by a dashed line), gas cells immediately turn into molecular cloud particles, and thus no cell is left behind unresolved.
• 2.  
  On the right, the gas cells of density >10⁻²⁴  g  cm⁻³ are now heated up to 10⁶–10⁷ K, populating zone "A." These cells are mostly located on the disk in relatively close proximity (<2 kpc) to the central MBH. These cells are very stable against fragmentation because they are so hot that they can barely cool down to a typical molecular cloud temperature in a dynamical time (i.e., t_cool ≫ t_dyn). For that reason, the cells have hard time to fulfill the condition (3) of the molecular cloud formation criteria described in Section 2.4. The heating by the X-ray radiation thus increases the amount of dense gas in the vicinity of the MBH which is incapable of condensing into stars.¹⁷

Zoom In Zoom Out Reset image size

Therefore, the feedback from even a slowly growing MBH retains hot dense gas at a galactic center which otherwise could have created strong star formation. Figure 9 dramatically demonstrates the distinct changes in radially averaged density profiles when MBH feedback is included.

• 1.  
  The left figure of gas density profiles displays that the radiation from the MBH keeps a substantial amount of gas at the core of the galaxy and inhibits the gas from all turning into stars. The total gas density of Sim-RMF (the green dashed line) at 0.1 kpc from the MBH is about ∼5 times as high as that of Sim-SF (the red solid line). The thin blue dotted line represents the initial profile of gas at z = 3 when the simulation restarts with 15 pc resolution. However in Sim-RMF, only a minimal fraction of the gas within 0.2 kpc radius is below 10⁴ K (the cyan dot-dashed line) whereas in Sim-SF, 10%–30% of the gas in the same region is considered to be cold, thus potentially star forming (the pink dotted line).
• 2.  
  The deprivation of cold dense gas in Sim-RMF inevitably prompts the suppression of star formation activity in the inner core of the galaxy. The right figure reveals the SFR density (in M_☉ yr⁻¹ kpc⁻³) as a function of distance from the MBH. Equation (7) is again used to calculate the new stellar mass being generated in each cell. The SFR density of Sim-RMF at 0.1 kpc from the MBH is reduced by more than ∼50% when compared with that of Sim-SF.¹⁸   Overall, the X-ray radiation from the MBH severely suppresses the SFR of Sim-RMF inside the 0.2 kpc sphere.

Zoom In Zoom Out Reset image size

The gas masses enclosed within a radius r, M_gas(< r), are plotted in Figure 10, showing the impact of MBH feedback as a powerful energy source to reshape the galactic gas distribution. Sim-RMF harbors ∼2.5 times more gas (1.3 × 10⁸ M_☉) within 0.2 kpc than what Sim-SF does (5.1 × 10⁷ M_☉), but almost none of this gas is below 10⁴ K. This gas in the core is not consumed by star formation, nor is pushed away by any mechanical outflow. Note that in Sim-RMF cold gas mass within 10 kpc is reduced by ∼50% displaying how far the MBH radiation reaches.

Zoom In Zoom Out Reset image size

Readers should also note that the enclosed gas masses at the virial radius (80 kpc proper) are almost identical between Sim-SF and Sim-RMF, implying there has been no massive gas expulsion driven by the MBH. This is one of the key differences from previous numerical studies. Very strong gas expulsions were frequently observed in previous studies in which the MBH feedback energy is only thermally deposited to a few neighboring gas particles (e.g., Springel et al. 2005b). In contrast, most of the gas is still bound to our simulated galaxies because (1) the energy released from our slowly growing MBH is relatively small (a possible "radio-mode" analog),¹⁹ (2) the gas mass to which the MBH energy is coupled is large in our radiative feedback formalism, and (3) the mass accretion rate onto our MBH is not high enough to repeatedly drive jets. Note again that the mechanical channel of MBH feedback is not a main driver of feedback in the presented calculation. The MBH has not doubled its mass at the end of our calculation (after 350 Myrs; see Section 4.4); and, with this highly suppressed mass accretion rate, jets have launched only a few tens of times in 350 Myr. We come back to the efficiency issue of mechanical feedback in Section 5.2.

To summarize, we have shown that MBH feedback, especially its radiation, alters the multiphase ISM of the surrounding gas and thus deprives the galactic inner core of cold, dense star-forming gas. Two consequences arise from the lack of star-forming gas at the galactic center: locally suppressed star formation and the associated change in stellar distribution. We discuss these topics in the following sections.

The inner core of the galaxy in Sim-RMF becomes a sterile environment for molecular cloud formation (star formation) because the gas is hot and turbulent, therefore Toomre stable. As a consequence, star formation is suppressed locally in the inner core, as shown in Figure 9. Figure 11 displays how the stellar mass density profile changes in the inner core region as a result of the locally suppressed star formation in Sim-RMF. Again we use the snapshot at z = 2.75, about 220 Myr after the MBH is placed at the center of the model galaxy. The stellar mass density at 0.1 kpc in Sim-RMF is only less than ∼20% of that of Sim-SF. The mass density of young stars (age <100 Myr) shows the similar drastic reduction.

Zoom In Zoom Out Reset image size

The stellar masses enclosed within a radius r, M_star(< r), are shown in Figure 12. The stellar mass enclosed within the 0.1 kpc sphere of Sim-RMF is almost an order of magnitude smaller than that of Sim-SF. The suppressed star formation replaces the steep inner cusp of stellar density profile with a flattened core ∼0.3 kpc in radius, i.e., stars are less concentrated at the galactic center of Sim-RMF. Considering the relatively small amount of energy released from the slowly growing MBH, the difference between these two lines is quite remarkable. Together with Figure 10, one expects that the stellar to gas mass ratio inside the 0.2 kpc sphere of Sim-RMF will be much smaller than that of Sim-SF. Note also that the total stellar mass enclosed at the virial radius (80 kpc) is almost indistinguishable between Sim-SF and Sim-RMF. In other words, star formation in Sim-RMF is not globally suppressed, but only locally suppressed at the center. This is because our MBH feedback is not strong enough to unbind a large amount of gas (see Section 4.2) or to globally abolish cold, star-forming clumps in the entire ISM.

Zoom In Zoom Out Reset image size

Figure 13 exhibits the evolution of the new stellar mass (molecular cloud particles born after the MBH is placed at 2190 Myr, or at z = 3) inside a 0.2 kpc sphere centered on the MBH. The plot demonstrates that since 2300 Myr the star formation activity in the inner core of Sim-RMF is suppressed. Naturally, alteration in the stellar distribution ensues at the center of the galaxy with an active MBH. Figure 14 strikingly contrasts the distribution of newly formed stars (molecular cloud particles of age <10 Myr).

• 1.  
  In the top row, a three-dimensional rendering of newly formed particles is constructed at a ∼45° angle from the disk plane. Here the light intensities of newly formed particles are integrated along the lines of sight. At the center of the stellar distribution, i.e., the densest peak, lies the MBH particle.
• 2.  
  In the bottom row, newly formed stellar masses are projected along the z-axis of the simulation box which makes a ∼472 angle with the angular momentum vector of the gas within a 5 kpc sphere centered on the MBH. The morphological difference at the inner core of the galaxy is particularly evident. Stars in Sim-SF are highly concentrated at the center, while stars in Sim-RMF are less concentrated but form spiral-like structures at ∼0.2 kpc radius from the MBH.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, it is shown that MBH feedback suppresses star formation locally at the galactic inner core, thus significantly changing the stellar distribution there. This new channel of feedback is particularly interesting because it is dominant only in the local surroundings of the MBH. Unlike stellar feedback, which operates globally, this new suppression mechanism does not require additional star formation and/or extensive mass expulsion out of the galactic potential.

Heating by MBH feedback, which locally suppresses star formation, also makes the MBH to self-regulate its own growth. Figure 15 shows the MBH accretion history versus time for Sim-SF and Sim-RMF.

• 1.  
  The MBH in Sim-SF has grown exponentially to 3 × 10⁶ M_☉ in 350 Myr, already about >30 times more massive than its initial mass of 10⁵ M_☉. The MBH maintained the accretion rate of 0.2–0.6 $\dot{M}_{\rm Edd}$ during this time period, corresponding to the unhindered growth of the MBH when there is no mechanism to self-regulate itself other than stellar feedback.
• 2.  
  Over the course of the same period the MBH in Sim-RMF has grown by only ∼70% to 1.7 × 10⁵ M_☉. The heated and diffused ISM in the vicinity of the MBH considerably suppresses the Bondi–Hoyle accretion estimate to as low as ${\sim} 0.02 \,\, \dot{M}_{\rm Edd}$ in this period. This indicates that the MBH feedback described in previous sections is a possible "radio-mode" analog where the accretion rate is ${\sim} 0.05 \,\dot{M}_{\rm Edd}$ (Croton et al. 2006). It also presents a potential route to the relatively low-mass MBH at the center of the Milky Way (3 × 10⁶ M_☉).²⁰

Zoom In Zoom Out Reset image size

Therefore, the feedback from an MBH is confirmed as an effective mechanism for slowing down the accretion of gas onto itself. Without having to suddenly unbind all the surrounding gas, the MBH self-regulates its growth by heating up the neighborhood and keeping it hot for an extended period of time. This finding is consistent with the work by Debuhr et al. (2011), who claimed that the growth of an MBH could be "self-regulated," rather than "supply-limited" (as in Springel et al. 2005b; Teyssier et al. 2010) where quasar-like MBH feedback drive energetic large-scale outflows to unbind a significant amount of gas.

A state-of-the-art numerical framework which fully incorporates gas, stars, and a central MBH is developed. Our simulation, for the first time, followed the comprehensive evolution of a massive star-forming galaxy with self-consistently modeled stars and an MBH. Our novel framework renders a completely different, yet physically more accurate picture of how a galaxy and its embedded MBH evolve under each others influence, providing a powerful means in understanding the coevolution of galaxies and MBHs. Our main results and new advancements are as follows.

• 1.  
  Molecular cloud formation and feedback. We have included a new model of molecular cloud formation and stellar feedback in our code (Sections 2.4 and 2.5). Unlike previous star formation recipes based on the Schmidt relation, a particle spawns when a gas cell of a typical molecular cloud size, 15.2 pc, actually becomes Jeans unstable. Then the molecular cloud particle gradually produces stellar mass while returning a large fraction of mass back to the gas with thermal feedback energy, modeling the observed slow star formation in molecular clouds. Thermal stellar feedback is shown to self-regulates star formation (Section 4.1).
• 2.  
  Massive black hole accretion and feedback. We have successfully developed a self-consistent model of accretion of gas onto an MBH and its radiative and mechanical feedback effects (Sections 2.6–2.8). Gas accretion onto the MBH is estimated with the Bondi–Hoyle formula, but without any boost factor, as we begin to resolve the Bondi radius. Monochromatic X-ray photons from the MBH are followed through three-dimensional adaptive ray tracing, rendering the radiative feedback of an MBH; here, rays of photons ionize and heat the gas, and exert momentum onto the gas. Finally, the mechanical feedback of the MBH is represented by bipolar jets with velocities of ∼10³  km  s⁻¹ launched from the vicinity of the MBH accretion disk, well resolved in our high-resolution AMR simulations. Our approach is significantly different from the previous recipes for MBH feedback in galactic scale simulations to date (e.g., Springel et al. 2005b; Sijacki et al. 2007; Booth & Schaye 2009; Teyssier et al. 2010; Debuhr et al. 2011), yet more accurately presents the physics of MBHs when properly incorporated with a high dynamic range.
• 3.  
  Locally suppressed star formation. By investigating the coevolution of a 9.2 × 10¹¹ M_☉ galactic halo and its 10⁵ M_☉ embedded MBH at z ∼ 3, we show that MBH feedback, especially its radiation, heats the surrounding ISM up to 10⁶ K through photoionization and Compton heating and thus locally suppresses star formation in the inner core of a galaxy (Section 4.2). The feedback also considerably changes the stellar distribution at the galactic center. This new channel of feedback from a slowly growing MBH is particularly interesting because it is only locally dominant, and does not require the heating of gas globally on the disk, or instigate a massive gas expulsion out of the galactic potential (Section 4.3).
• 4.  
  Self-regulated black hole growth. MBH feedback is also demonstrated to be an effective mechanism for slowing down the accretion of gas onto the MBH itself. Without necessarily unbinding all of its surrounding gas, the MBH self-regulates its growth by keeping the surrounding ISM hot for an extended period of time (Section 4.4). Therefore, our results possibly are consistent with a "radio-mode" analog of MBH feedback.

Our method limits the use of ad hoc formulation and instead more accurately models the physics of galaxy formation. As a result, four key components of galactic scale physics, (1) molecular cloud formation, (2) stellar feedback, (3) MBH accretion, and (4) MBH feedback, work self-consistently in one comprehensive framework. As an example, the radiation and jets from the MBH heat up the surrounding gas and create hot regions, but the thermal couplings of the radiative and mechanical energy are all carried out by the shock-capturing radiation hydrodynamics AMR scheme itself, not by any presupposed thermal deposition model. In our framework, one should also be able to couple small-scale physics (such as molecular cloud formation and feedback) with large-scale physics (such as quasar-driven galactic outflows) without any sub-resolution model. These first results undoubtedly demonstrate that we can now develop an unabridged, self-consistent numerical framework for both galaxies and MBHs.

While proven to be fruitful already in producing robust results, our comprehensive galaxy formation framework is only the first step forward in the right direction. Imminent future projects and improvements are as follows.

• 1.  
  Merging galaxies. What is experimented in this work is a quiescent form of MBH feedback, possibly a "radio-mode" analog. Meanwhile, MBH feedback is known to make a dramatic difference in galaxy mergers, which is frequent in hierarchical structure formation. The gas funneled into galactic centers triggers "quasar-mode" MBH feedback (Hopkins et al. 2008), which subsequently reshapes the relation between black hole masses and bulge stellar distributions (Di Matteo et al. 2005; Johansson et al. 2009). Detailed research on this more intense mode of MBH feedback will be the topic of a future paper.
• 2.  
  Parameter studies. More comprehensive parameter studies should follow, especially in the parameterization of MBH feedback and the efficiency of stellar feedback. The results should be compared and calibrated with observations such as bulge to disk mass ratio and gas to stellar mass ratio (e.g., Gavazzi et al. 2008) or with analytical investigations (e.g., Sazonov et al. 2005). In particular, the disk-bulge decomposition of simulated galaxies will be the subject of subsequent analysis of our simulations.
• 3.  
  Improving mechanical feedback. In the results presented herein, mechanical feedback is energetically secondary to radiative feedback because the mass accreted onto the MBH is not large enough to repeatedly drive jets. These infrequent jets easily penetrate the ISM without necessarily creating sizable shocks or entraining a large amount of gas. However, a few mechanisms will be considered in the future which could have enhanced the effectiveness of jets. Magnetic fields could aid the jets in efficiently depositing outflow momentum onto the infalling gas, as was shown by the studies on the evolution of jets in the presence of magnetic fields (Dubois et al. 2009; Wang et al. 2010). Cosmic rays accelerated by relativistic jets and shock fronts (Skillman et al. 2008) could boost the effectiveness of jets, too.
• 4.  
  Improving radiative feedback. For now, monochromatic X-ray photons are utilized to carry the energy of MBH radiation (Section 2.7); however, a better model will be needed to describe the polychromatic energy distribution of MBH radiation. Ideally one wants to have a large number of spectral energy bins, each of which is separately followed through three-dimensional ray tracing. Given the computationally challenging nature of polychromatic radiative transfer, however, tabulated rates of photoionization and photoheating as functions of optical depth can be a good alternative. Moreover, to accurately quantify the radiative feedback on the gas in the vicinity of an MBH, the pressure force on dust grains needs to be computed. This could have increased the radiation pressure in the presented results, especially in the central <kpc region. For this purpose, dust models in Rocha et al. (2008) will need to be considered.
• 5.  
  Adding supplementary feedback channels. An MBH particle in our work represents not just the black hole itself, but also includes accreting gas and stars deep within the galactic nucleus. Thus, there is a need for other feedback channels, such as stellar winds from a nuclear disk (Ciotti et al. 2009). The nuclear disk winds can be implemented as thermal deposition of energy, working in conjunction with the aforementioned radiative and mechanical feedback. Stellar UV radiation from the nuclear disk can also be incorporated into the radiative feedback of the MBH. Including this supplementary feedback will reveal the multi-faceted nature of the coupling of MBH energy with its surroundings.
• 6.  
  Improving accretion estimate. The accretion estimate using the Bondi–Hoyle formula will need to be improved, especially when the gas disk around the MBH can be resolved down to the Bondi radius. Different estimates such as the ones considering gas angular momentum (Hopkins & Quataert 2010; Levine et al. 2010) are attractive candidates that should be explored.
• 7.  
  Nonthermal pressure sources. Nonthermal pressure sources such as magnetic fields (Wang & Abel 2009), stellar UV radiation, and cosmic rays are missing in this work, but should be included in future simulations.

We also recognize that the results from our experiment can provide the community with better sub-resolution models for MBH physics. For example, the radial profile of heating rates by the MBH in our simulation can be tabulated; in a coarsely resolved particle-based simulation, one can deposit thermal energy according to this radial dependence into a volume larger than a typical smoothing kernel. This can be a useful means for improving the particle-based simulations as well as for speeding up future large-scale AMR calculations, such as the formation of high-redshift quasars and the reionization of intergalactic helium.