The impact of radio feedback from active galactic nuclei in cosmological simulations: Formation of disk galaxies

TL;DR

This paper introduces a new AGN feedback model in cosmological galaxy formation simulations, focusing on RIAFs, which suppress star formation and produce properties consistent with observed Seyfert galaxies.

Contribution

The novel aspect is the implementation of RIAF-based AGN feedback in simulations, linking accretion modes to galaxy evolution and black hole growth.

Findings

AGN feedback suppresses star formation in bulges

Model reproduces properties of Seyfert galaxies

Black hole to spheroid mass ratio ~10^{-3} at z=0

Abstract

In this paper, we present a new implementation of feedback due to active galactic nuclei (AGN) in cosmological simulations of galaxy formation. We assume that a fraction of jet energy, which is generated by an AGN, is transferred to the surrounding gas as thermal energy. Combining a theoretical model of mass accretion onto black holes with a multiphase description of star-forming gas, we self-consistently follow evolution of both galaxies and their central black holes. The novelty in our model is that we consider two distinct accretion modes: standard radiatively efficient thin accretion disks and radiatively inefficient accretion flows which we will generically refer to as RIAFs; motivated by theoretical models for jet production in accretion disks, we assume that only the RIAF is responsible for the AGN feedback. We find that, after an initial episode of bursting star formation, the accretion rate onto the central black hole drops so that the accretion disk switches to a RIAF structure. At this point, the feedback from the AGN becomes efficient and slightly suppresses star formation in the galactic disk and almost completely halts star formation in the bulge. As a result, the nucleus becomes a stochastically fuelled low-luminosity AGN (Seyfert galaxy) with recurrent short-lived episodes of activity after the star bursts. Our model predicts several properties of the low-luminosity AGN including the bolometric luminosity, jet powers, the effect on kpc-scale of the radio jet and the AGN lifetime, which are in broad agreement with observations of Seyfert galaxies and their radio activity. We also find that the mass ratios between the central black hole and the the host spheroid at z = 0 are ~10^{-3} regardless of the strength of either supernova feedback or AGN feedback. (abridged)

Figures

Tables

Table 1. Table 1: Total and spheroid’s mass, g 𝑔 g and r 𝑟 r band total luminosities, and disk-to-total mass and luminosity ratios for the galaxies at z = 0 𝑧 0 z=0 .

	$M_{\mathrm{tot}}/M_{\odot }$	$M_{\mathrm{SP}}/M_{\odot }$	$M_g$	$M_r$	${\left(\frac{D}{T}\right)}_{\mathrm{Mass}}$	${\left(\frac{D}{T}\right)}_g$	${\left(\frac{D}{T}\right)}_r$
no-fb	$6.01\times {10}^{10}$	$1.98\times {10}^{10}$	-20.8	-21.3	0.67	0.83	0.79
weak-fb	$5.98\times {10}^{10}$	$1.44\times {10}^{10}$	-20.9	-21.3	0.76	0.88	0.85
reference	$4.32\times {10}^{10}$	$1.21\times {10}^{10}$	-20.6	-21.0	0.72	0.86	0.83

Equations

$${\frac{\mathrm{d}{M}_{\mathrm{c}}}{\mathrm{d}t}|}_{\mathrm{TI}}=-{\frac{\mathrm{d}{M}_{\mathrm{h}}}{\mathrm{d}t}|}_{\mathrm{TI}}=\frac{{\Lambda }_{\mathrm{net}}({\rho }_{\mathrm{h}},u_{\mathrm{h}},Z)}{u_{\mathrm{h}}-u_{\mathrm{c}}}\frac{M_{\mathrm{h}}}{{\rho }_{\mathrm{h}}},$$

$$\varphi (m)\equiv \frac{\mathrm{d}{N}_{\mathrm{c}}}{\mathrm{d}m}\propto m^{-\alpha }.$$

$$\frac{m}{r^2}=190\left(\frac{M_{\odot }}{\text{p}{c}^2}\right)P_4^{1/2},$$

$${\dot{m}}_{\ast }=\frac{\mathrm{d}{m}_{\ast }}{\mathrm{d}t}=\frac{m}{t_{\mathrm{sf}}}=c_{\ast }\frac{m}{t_{\mathrm{dyn}}},$$

$$t_{\mathrm{dyn}}={\left[\frac{3\pi }{32G\rho (m)}\right]}^{1/2}\simeq 0.32P_4^{-3/8}{\left(\frac{m}{M_{\odot }}\right)}^{1/4}\mathrm{Myr},$$

$${\dot{M}}_{\ast }={\int }_{10^4{M}_{\odot }}^{10^6{M}_{\odot }}{\dot{m}}_{\ast }(m)\varphi (m)\mathrm{d}m,$$

$${\int }_{10^2{M}_{\odot }}^{10^6{M}_{\odot }}m\varphi (m)\mathrm{d}m=1$$

$${\dot{m}}_{\mathrm{EVP}}$$

$${\dot{M}}_{\mathrm{EVP}}={\int }_{10^2{M}_{\odot }}^{10^6{M}_{\odot }}{\dot{m}}_{\mathrm{EVP}}(m)\varphi (m)\mathrm{d}m.$$

$$p_{\ast }=\frac{M_{\mathrm{c}}}{m_{\ast }}\left[1-\exp \left(-\frac{\Delta t}{t_{\ast }}\right)\right],$$

$$t_{\ast }(P)={\left(\frac{{\dot{M}}_{\ast }(P)}{M_{\odot }}\right)}^{-1}\mathrm{Myr}.$$

$$p_{\mathrm{II}}=\frac{{\int }_t^{t+\Delta t}{r}_{\mathrm{II}}(t^{\prime })\mathrm{d}{t}^{\prime }}{{\int }_t^{t(8M_{\odot })}{r}_{\mathrm{II}}(t^{\prime })\mathrm{d}{t}^{\prime }},$$

$$p_{\mathrm{Ia}}=\frac{{\int }_t^{t+\Delta t}{r}_{\mathrm{Ia}}(t^{\prime })\mathrm{d}{t}^{\prime }}{{\int }_t^{t_0+\Delta t_{\mathrm{Ia}}}{r}_{\mathrm{Ia}}(t^{\prime })\mathrm{d}{t}^{\prime }},$$

$${M_{\mathrm{h}}}^{\prime }=M_{\mathrm{h}}+\Delta M_{\mathrm{EVP}}-\frac{{\Lambda }_{\mathrm{net}}({{\rho }_{\mathrm{h}}}^{\prime },u_{\mathrm{h}},Z)}{u_{\mathrm{h}}-u_{\mathrm{c}}}\frac{M_{\mathrm{h}}^{\prime }}{{{\rho }_{\mathrm{h}}}^{\prime }}\Delta t,$$

$$\frac{M_{\mathrm{SPH}}}{\rho }={Vol}_{\mathrm{c}}+\frac{M_{\mathrm{h}}}{{\rho }_{\mathrm{h}}},$$

$${u_{\mathrm{h}}}^{\prime }=u_{\mathrm{h}}+\frac{\Delta Q+(u_{\mathrm{c}}-{u_{\mathrm{h}}}^{\prime })\Delta M_{\mathrm{EVP}}}{{M_{\mathrm{h}}}^{\prime }},$$

$$P_{\min }(\rho )\propto {\rho }^{{\gamma }_{\mathrm{eff}}}\mathrm{for}\rho >{\rho }_{\mathrm{th}},$$

$$p_{\mathrm{w}}=\frac{\frac{1}{2}{M}_{\mathrm{SPH}}{v_{\mathrm{w}}}^2}{\Delta Q},$$

$${\mathbf{v}}^{\prime }=\mathbf{v}+v_{\mathrm{w}}\mathbf{n}.$$

$$j(r)\propto \frac{M(r)}{M_{\mathrm{vir}}}.$$

$${\tau }_{\mathrm{SFR}}(t)\simeq \frac{3{\chi }_{\mathrm{d}}}{4\pi }\frac{M_{\mathrm{c}}(t)}{R_{\mathrm{SFR}}(t{)}^2},$$

$${\dot{M}}_{\mathrm{drag}}={\eta }_{\mathrm{drag}}\frac{L_{\mathrm{SFR}}(t)}{c^2}(1-e^{-{\tau }_{\mathrm{SFR}}(t)}),$$

$$L_{\mathrm{jet}}^{\mathrm{SD}}$$

$$L_{\mathrm{jet}}^{\mathrm{RIAF}}$$

$$L_{\mathrm{jet}}^{\mathrm{RIAF}}\approx 2.6\times 10^{-1}\dot{M}{c}^2\mathrm{for}\dot{m}\le {\dot{m}}_{\mathrm{crit}},$$

$$L_{\mathrm{jet}}^{\mathrm{SD}}\approx 8.1\times 10^{-5}\dot{M}{c}^2\mathrm{for}\dot{m}>{\dot{m}}_{\mathrm{crit}},$$

$$p_n=\frac{M_{\mathrm{BH}}-M_{\mathrm{p}}}{N_{\mathrm{SFR}}}\frac{1}{M_n},$$

$$\Delta l$ = ${\rm MIN}(0.01\epsilon,0.03\ |\mbox{\boldmath$v$}|\ \Delta t),$$

$${\eta }_{\mathrm{RIAF}}=0.1\left(\frac{\dot{m}}{{\dot{m}}_{\mathrm{crit}}}+\delta +10^{-3}\right)$$

Full text

††pagerange: The impact of radio feedback from active galactic nuclei in cosmological simulations: Formation of disk galaxies–References††pubyear: 2007

1 Introduction

Recent observations have suggested a fundamental connection between active galactic nuclei (AGN) and the formation of galaxies. Firstly, there is now a well established correlation between the properties of galaxies and the masses of the black holes (BH) at their centres. The mass of the central black holes tightly correlates with the mass of the host galactic bulges (with a median BH-to-bulge mass ratio of $\simeq 0.001$, e.g. Kormendy & Richstone, 1995; Magorrian et al., 1998; Merritt & Ferrarese, 2001; McLure & Dunlop, 2002; Marconi & Hunt, 2003), as well as with the stellar velocity dispersion of the bulges (Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002). This discovery points to a fundamental connection between the growth of central BHs and the formation of stellar spheroids in galaxies. Secondly, high resolution X-ray observations of galaxy clusters have revealed large, radio-plasma cavities in intra-cluster medium. These are usually associated with episodic outbursts from a central radio galaxy and indicate that huge amounts of mechanical energy are being deposited into the intracluster medium by powerful AGN-driven jets (e.g. Allen et al., 2006; Bîrzan et al., 2004; Fabian et al., 2006; Rafferty et al., 2006; Taylor et al., 2006). Simulations of the impact of these AGN-driven radio cavities suggest that this powerful feedback provides sufficient energy to offset the cooling radiation from the cluster and potentially explain why so little cool ($T<3\,{\rm keV}$) gas is seen in these systems (Quilis et al., 2001; Churazov et al., 2002; Dalla Vecchia et al., 2004; Omma et al., 2004; Sijacki & Springel, 2006; Sijacki et al., 2007).

Motivated by these discoveries, simple prescriptions for the feedback from AGN (aka. “radio-mode feedback”) have recently been incorporated into semi-analytic galaxy formation models (e.g. Granato et al., 2004; Monaco & Fontanot, 2005; Cattaneo et al., 2005; Croton et al., 2006; De Lucia et al., 2006; Bower et al., 2006). Including radio-mode AGN feedback has resulted in dramatic improvements in the models’ ability to match the sharp decline of the galaxy luminosity function and to explain the “down-sizing” seen in the evolution of the galaxy population. In particular, Bower et al. shows that, by assuming that AGN feedback operates only in quasi-hydrostatic halos where the cooling time is longer than the dynamical time (Dalla Vecchia, 2005; Sijacki & Springel, 2006), the galaxy luminosity functions in local and higher redshifts universe can be matched well. Croton et al. (2006) and Kitzbichler & White (2007) showed that similar results were obtained by modifying the Bondi-Hoyle-Lyttleton accretion rate formula (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952) to account for the multiphase structure of the accreting gas in rapidly cooling halos.

Recent simulations also have begun to track the impact of AGN feedback on the galaxy population (Kawata & Gibson, 2005; Sijacki & Springel, 2006; Di Matteo et al., 2005; Springel et al., 2005b; Sijacki et al., 2007; Di Matteo et al., 2007). One of the first hydrodynamic simulation of galaxy formation which invoked AGN feedback was performed by Kawata & Gibson (2005) by using cosmological simulations with smoothed particle hydrodynamics (SPH). They reproduced observed X-ray and optical properties of elliptical galaxies by injecting thermal energy into the centre of the main progenitor at $z<1$, assuming that sufficient BH mass was present for the AGN to be active when a convergent gas inflow exists at the centre of the main progenitor. A self-consistent model for AGN “quasar-mode” feedback associated with BH accretion in simulations of galaxy mergers was proposed by Di Matteo et al. (2005) and Springel et al. (2005b). They estimated accretion rate onto BHs by using a Bondi-Hoyle-Lyttleton parameterisation and injected some fraction of accreted rest mass energy into the surrounding interstellar medium (ISM) in the form of thermal energy. By using this model and performing a series of merger simulations, Hopkins et al. (2006) simulated evolution of quasar luminosity and predicted luminosity functions of quasars.

In this paper, we introduce a new methodology to incorporate BH growth and AGN radio feedback self-consistently in cosmological simulations of galaxy formation. We are motivated to do this by the recent semi-analytic models that we have discussed above. These indicate that the distinction between “quasar” and “radio” modes is necessary to achieve good matches to the global properties of galaxies. We distinguish these two modes depending on the accretion rate onto the central BHs assuming that the radio-mode exists only when the accretion rate onto the BH is low.

We estimate the mass accretion rate onto the central BHs by combining a radiation drag model by Kawakatu & Umemura (2002) with our description of the ISM and star formation. A virtue of the radiation drag model is that it relates star formation activity in galactic centre, which can be resolved in our simulations, to mass accretion onto a central BHs. It implies that a fixed fraction of the star formation in bulges is converted into BH mass just as assumed in most semi-analytic models (e.g. Croton et al., 2006; Bower et al., 2006). Since our estimate of the accretion rate onto central BHs is based on radiation from stars, the mass growth of BHs in our simulations is resolution independent so long as the star formation is modelled in a resolution independent fashion.

In order to model the AGN radio-mode feedback, we distinguish two fundamentally different modes of AGN accretion: radiatively efficient geometrically thin accretion flows (standard disks; Shakura & Sunyaev, 1973) and geometrically thick, radiatively inefficient accretion flows (which we will generically refer to as RIAFs; Narayan et al., 1998; Narayan, 2005; Nemmen et al., 2006; Yuan, 2007). We assume only the latter are responsible for the radio-mode feedback through production of powerful jets (e.g. Rees et al., 1982; Meier, 2001; Maccarone et al., 2003; Churazov et al., 2005). This assertion is supported by theoretical models for jet production since the jet power is intimately linked to the scale of the magnetic field loops in the vicinity of the last stable orbit around the BH (e.g. Livio et al., 1999; Meier, 2001; Nemmen et al., 2007). Since RIAFs exist only when the accretion rate is much lower than the Eddington rate (Narayan et al., 1998), AGNs can produce strong feedback in objects which have low specific star formation rates because the radiation drag model directly connects the star formation rate around BHs and mass accretion rate onto the BHs. Given the limited numerical resolution inherent in our simulations, we adopt a subgrid model for the impact of the AGN feedback. We simply assume that a certain fraction of the available jet power thermally couples to nearby diffuse halo gas.

Until today, most of simulations of galaxy formation with AGN feedback aimed to form elliptical galaxies because the AGN feedback is the most plausible mechanism that halts the star formation activity in elliptical galaxies (Kawata & Gibson, 2005; Springel et al., 2005b; Di Matteo et al., 2005; Sijacki et al., 2007). Therefore roles of AGN feedback in disk galaxy formation are still largely unknown in spite that disk galaxies also harbour central supermassive BHs in their bulges. In this paper we hence focus on the impact that feedback from an AGN jet might have on the formation of a large disk galaxy, and the resulting coevolution of the AGN and the host galaxy. The majority of AGN hosted by disk galaxies are known to be either Seyfert nuclei or the lowest luminosity AGN, the low-ionisation nuclear emission-line regions (LINERs) (Ho et al. 1997; see Ho 2004 for a review). Several targeted observations of individual objects (e.g., Wilson & Ulvestad 1983; Morganti, Oosterloo, & Tsvetanov 1998; Cecil et al. 2000; Whittle & Wilson 2004; Keel et al. 2006; Kharb et al. 2006; Middelberg et al. 2007) indicate that Seyfert galaxies commonly host compact radio jets spanning $\sim 10-100$ pc which strongly affect the state of the ISM and contribute to its ionisation via shocks (e.g., Kukula et al. 1993, 1996; Falcke, Wilson, & Simpson 1998; Capetti et al. 1999; Bicknell et al. 1998; Mundell et al. 2003; Riffel et al. 2006). It is also common that kpc-scale radio structures are observed, which are thought as extensions of the inner radio jets (e.g., Wilson & Ulvestad 1983; Morganti, Oosterloo, & Tsvetanov 1998; Cecil et al. 2000; Whittle & Wilson 2004; Keel et al. 2006; Kharb et al. 2006; Middelberg et al. 2007), and sometimes blow out of the galactic disk with no preferential direction (Schmitt et al., 2001; Gallimore et al., 2006). Very Large Array surveys of samples of Seyfert and LINER galaxies indicate that radio outflows with extent $\gtrsim 1$ kpc are a common feature in Seyfert galaxies, and presumably are driven by the AGN jets (Colbert et al., 1996; Gallimore et al., 2006). Gallimore et al. (2006) in particular find that $\gtrsim 44\%$ of the galaxies in their complete sample display kpc-scale radio outflows, and they cannot rule out that most Seyfert galaxies produce extended radio outflows. The kinetic power carried by the Seyfert jets as estimated from observations can be as high as $10^{43}{\rm erg}\ {\rm s}^{-1}$ (e.g., Wilson & Ulvestad 1983; Morganti, Oosterloo, & Tsvetanov 1998; Kharb et al. 2006). If an appreciable fraction of this power reaches out of the nuclear region, then the jets may affect the environment of the AGN and be a source of feedback in Seyfert galaxies. The impact of feedback from AGN jets on the nuclear environment of the host disk galaxies and the evolution of the nuclear activity during the formation of the disk galaxy are open problems, which we attack using our cosmological simulations. In this paper, we show that distinguishing the two modes of AGN accretion yields results which are in broad agreement with observations of Seyferts.

The paper is organised as follows. We first describe our microscopic descriptions of the ISM, star formation, and stellar feedback in Section 2 and BH growth and AGN feedback in Section 3. We present the details of our cosmological simulations in Section 4 and the results in Section 5. Finally, we summarise and discuss our results in Section 6. Readers who are observation-oriented and are not interested in the technical details of the simulation may wish to skip directly to Section 6, where we compare our results with observations of nuclear activity in Seyfert galaxies.

2 The model of star-forming gas and stellar feedback

In order to study the effect of AGN feedback on galaxy formation, we construct a physically motivated, well-controlled model of the star-forming gas. Based on the simplest star formation and feedback recipes used in the first generation of hydrodynamic simulations of galaxy formation (Navarro & White, 1993; Steinmetz & Mueller, 1994; Cen & Ostriker, 1992; Katz et al., 1996), several authors have introduced ‘multiphase’ models for star-forming gas in which the ISM is treated as a number of distinct phases by formulating differential equations that model the interactions between the phases (Yepes et al. 1997; Springel & Hernquist 2003 hereafter SH03 ; Samland & Gerhard 2003; OEFJ05 ; Scannapieco et al. 2006 but see Booth et al. 2007 where they explicitly model ‘clouds’ by means of sticky particles in SPH simulations). Our model is based on the one used in OEFJ05 but has been significantly modified. We describe the details in the following subsections.

2.1 Gas cooling

We compute the non-equilibrium cooling/heating rate and ionisation state of each gas particle by using the cooling functions in Theuns et al. (2002) which include cooling by H, He, and metals, in the presence of an evolving but uniform ultraviolet background (Haardt & Madau, 1996) that is switched on at $z\simeq 6$. Inverse Compton cooling is also included.

Metals are carried by particles and once assigned to a particle, they remain with it. Nonetheless, effective mixing takes place because we use the smoothed metallicity (smoothed in the same way as other SPH quantities) when computing the cooling rate and defining the metallicity of a newly-formed star.

2.2 Cloud formation and the cloud model

As in SH03 and OEFJ05 we treat an SPH particle as a hybrid particle that consists of two distinct phases, i.e. hot ambient gas and cold clouds, once its gas density, $\rho$, exceeds a density threshold, $\rho_{\rm th}$.

We assume that the cold clouds form and grow through thermal instability, that is

[TABLE]

where $M_{\rm c}$ and $M_{\rm h}$ are mass in cold clouds and mass in the hot phase associated with the particle, respectively, $\Lambda_{\rm net}$ is the cooling function for gas of metallicity $Z$, and $u_{\rm h}$ and $u_{\rm c}$ represent specific internal energy of the hot phase and cold clouds, respectively. This implies that we assume that the gas is thermally unstable when $\rho>\rho_{\rm th}$. Cold clouds remain at a fixed temperature, $T_{\rm c}=100K$ hence $u_{\rm c}$ is a constant as well. Note, however, as far as $T_{\rm c}$ is much lower than $10^{4}$ K our results are hardly affected by the choice of $T_{\rm c}$.

We allow cold clouds to have a mass spectrum

[TABLE]

Observationally $\alpha$ is in a range between $\alpha=1.5$ to $1.9$ (Solomon & Rivolo, 1989; Fukui et al., 2001; Heyer et al., 2001). Since the effect of changing the value of $\alpha$ is completely degenerate with changing other parameters we will introduce such as star formation efficiency and thermal conduction efficiency, we simply fix $\alpha=1.7$ and the mass range of clouds, $10^{2}$–$10^{6}M_{\odot}$, throughout this paper.

Following Samland & Gerhard (2003), we assume spherical clouds which follow the mass-radius relation of Elmegreen (1989)

[TABLE]

where $P_{4}=\frac{P/k}{10^{4}{\rm K}~{}{\rm cm}^{3}}$ is the pressure in the ambient hot phase which is equal to the effective pressure of the particle.

2.3 Star formation

We assume that giant molecular clouds in mass range of $10^{4}M_{\odot}$ to $10^{6}M_{\odot}$ are eligible to form stars with a star formation timescale that is proportional to the dynamical time of each cloud. Therefore the star formation rate for a cloud of mass $m$ is

[TABLE]

where $m_{*}$ is stellar mass, $t_{\rm sf}=t_{\rm dyn}/c_{*}$ is the star formation timescale for a cloud of mass $m$, $c_{*}$ is a dimensionless star formation efficiency parameter that must be less than unity, and $t_{\rm dyn}$ is the dynamical timescale of the cloud. For the mass-radius relation in equation (3), the dynamical timescale of a cloud of mass $m$ is given as

[TABLE]

where $\rho(m)$ is the mean density of a cloud of mass $m$.

It should be noted that now the star formation timescale is a function of ambient pressure and it could be very short in high-pressure medium. Therefore we naturally have ‘shock-induced’ starbursts that OEFJ05 introduced by changing the star formation efficiency by hand when the rate of change of the entropy by shock heating exceeds a threshold value. The total star formation rate for clouds of total mass of $1\,M_{\odot}$ with the mass spectrum $\phi(m)$ is obtained by integrating $\dot{m}_{*}$ over the range of masses that are eligible to form stars,

[TABLE]

where the normalisation

[TABLE]

is applied.

2.4 Cloud evaporation

The different phases of the ISM exchange mass by evaporation. As in Samland & Gerhard (2003), we use the model of Cowie & McKee (1977) for the cloud evaporation, that is,

[TABLE]

where $\eta_{\rm EVP}$, $\mu$, $\kappa$, $r$, and $k_{\rm B}$ are the dimensionless conduction efficiency parameter, the mean molecular weight, the conductivity where we use the Spitzer value (Spitzer, 1962), the radius of the cloud, and the Boltzmann constant. The total evaporation rate for clouds of total mass $1\,M_{\odot}$ then becomes

[TABLE]

2.5 Supernova feedback and chemical enrichment

We consider each stellar particle as a single stellar population (SSP) having its own age, metallicity and initial mass function (IMF) $\phi_{*}(m)$. Throughout this paper we assume that the IMF is always the Salpeter IMF (Salpeter, 1955) whose lower and upper mass limits are $0.1$ and $100$ $M_{\odot}$, respectively. The recycled mass fraction, type II and type Ia supernova (SN) rates, and metal production rates for some species are calculated according to the age, metallicity, and IMF of the SSP. Each SN is assumed to give $\eta_{\rm SN}10^{51}$ erg of energy to surrounding gas, where $\eta_{\rm SN}$ is the feedback efficiency parameter. We refer the readers to OEFJ05 and references therein for a more detailed description.

2.6 Numerical implementation

We have implemented the physics described above into a publicly available PM-TreeSPH code GADGET2 (Springel, 2005), a successor of the TreeSPH code GADGET (Springel et al., 2001). We do not use the instantaneous recycling approximation (IRA) which is often assumed in studies of galaxy formation. In this section we describe how the physical processes described above are implemented in our code.

For $\rho>\rho_{\rm th}$, an SPH particle is eligible to form cold clouds. Hence it is eligible to form stars. We compute the probability $p_{*}$ of a SPH particle spawning a new star particle of mass $m_{*}$ during a time-step $\Delta t$ as

[TABLE]

where the average star formation timescale $t_{*}$ is obtained by using equation (6):

[TABLE]

We use $m_{*}=m_{\rm orig}/3$ as in OEFJ05 , where $m_{\rm orig}$ denotes the original SPH particle mass. When an SPH particle spawns a star particle, the mass in cold clouds is reduced by $m_{*}$.

We relax the IRA through a probabilistic treatment of SNe II and Ia. First, we assume that stars more massive than $8M_{\odot}$ explode as SNe II. The probability $p_{\rm II}$ of a star particle having an event of SN II explosion during a time-step $\Delta t$ is given by

[TABLE]

where $t$ is the age of the SSP, $r_{\rm II}(t)$ is the SN II rate for the SSP of age $t$ and $t(8M_{\odot})$ is the lifetime of the star of original mass $8M_{\odot}$. When a star particle satisfies the condition to have an event of SN II we distribute mass, metals, and feedback energy, which are expelled by the total SNe II between $t=0$ and $t(8M_{\odot})$ to the surrounding gas particles using the SPH kernel weighting.

Since the SN Ia rate is much lower than the SN II rate, we allow star particles to have only one SN Ia explosion event for each SN Ia time-step $\Delta t_{\rm Ia}$, which is much longer than the time-steps for dynamical calculation. By doing this, we substantially reduce the computational expense that is needed mainly for the neighbour search to distribute energy, mass, and metals. As for SN II, we can compute the SN Ia probability during a time-step $\Delta t$ as

[TABLE]

where $t_{0}$ is the age at which the previous SN Ia time-step finished and $r_{\rm Ia}(t)$ is the SN Ia rate for a SSP of age $t$. We adopt $\Delta t_{\rm Ia}=100$ Myr in this paper. The probabilistic method described above statistically relaxes the IRA.

Finally, we update the mass of the hot phase $M_{\rm h}$ and the specific internal energy of the hot phase $u_{\rm h}$. The new mass of the hot phase ${M_{\rm h}}^{\prime}$ is given by solving the thermal energy equation implicitly using equation (1),

[TABLE]

where $\Delta M_{\rm EVP}$ is the mass which evaporates from cold clouds during the time-step, and the new density of the hot phase ${\rho_{\rm h}}^{\prime}$ is given by the mass conservation, $M_{\rm SPH}=M_{\rm h}+M_{\rm c}$, and by the assumption that the volume occupied by a SPH particle remains the same during the time-step, i.e.

[TABLE]

where ${\it Vol}_{\rm c}$ is the volume occupied by the cold clouds computed using the mass-radius relation (eq. 3).

In addition to the adiabatic heating/cooling and shock heating, non-hydrodynamic processes such as SN feedback and cloud evaporation also change the specific energy of the hot phase. The new specific energy of the hot phase ${u_{\rm h}}^{\prime}$ due to the non-hydrodynamic processes is also calculated implicitly as

[TABLE]

where $\Delta Q$ is the thermal energy received by the SPH particle during the time-step.

Note however that we model the SN feedback not as the thermal heating of the ISM but as the kinetic heating described in the following subsection as winds. This is motivated by the results from high-resolution 2D and 3D simulations of the ISM (Wada & Norman, 2001; Wada, 2001; Wada & Norman, 2007) in which they showed that the ISM is supported by the turbulent motion originated in the self-gravity of the gas and galactic rotation. They found that energy injection by SNe hardly changes the structure of the ISM. Kravtsov (2003) also indicated that the energy feedback from SNe does not change the structure of the ISM but significantly reduces the amount of gas in the galaxy in their cosmological simulation performed with an adaptive mesh refinement code.

Lacking the heating by the $\Delta Q$ term, the ISM loses its pressure quickly and would be fragmented. We therefore introduce minimal heating by imposing the minimum pressure for the star forming gas as a function of density,

[TABLE]

by which we mimic the heating from self-gravity and galactic shear motion claimed by Wada & Norman (2001), Wada (2001), and Wada & Norman (2007). We adopt in this paper the value $\gamma_{\rm eff}=1.4$ for the effective adiabatic index, which is close to the value found by Wada & Norman (2007) for dense gas and is stable against gravitational instability (the critical value is 4/3). We normalise the minimum pressure as $P_{\rm min}(\rho_{\rm th})=(\gamma-1)u_{4}\rho_{\rm th}$, where $u_{4}$ denotes the specific internal energy of the gas with temperature $10^{4}$ K. Once pressure from the hot phase drops below $P_{\rm min}$, we set the temperature of the hot phase to $T_{\rm hot}$ and recalculate the mass in cold clouds so that the ambient pressure becomes $P_{\rm min}$. We set $T_{\rm hot}$ to $10^{6}$ K in this paper but our results are hardly affected by the value of $T_{\rm hot}$ as far as $T_{\rm hot}>>10^{4}$ K. In quiescent disks, most of the star forming gas has $P_{\rm min}$ by construction. However, for example during a violent merger, the pressure can be significantly higher than $P_{\rm min}$ due to shock heating and hence the gas can have very short star formation timescale.

2.7 Winds

As discussed in SH03 , we explicitly assume that cold clouds and hot ambient medium remain tightly coupled in high-density regions. Hence, there is no obvious route for the high entropy gas to escape from the star forming regions. We thus extend our feedback model to account for galactic winds by SN feedback. Our model is based on SH03 and is modified to take the non-instantaneous SN explosions into account.

First, we specify the velocity of the wind, $v_{\rm w}$, when it leaves the disk. Then, we assign the probability $p_{\rm w}$ with which the gas particle is added to the wind as

[TABLE]

where $M_{\rm SPH}$ is the mass of the gas particle and $\Delta Q$ is the feedback energy received during the time-step. By doing this, our model becomes insensitive to the numerical resolution, time-step in the simulations, and the number of neighbours we use to distribute the feedback energy. When a particle is added to the wind, we modify the velocity ${\mathbf{v}}$ of the particle according to

[TABLE]

For the unit vector ${\mathbf{n}}$, we choose random orientation along the direction ${\mathbf{v}}\times{\mathbf{a}}_{\rm grav}$, where ${\mathbf{a}}_{\rm grav}$ is the gravitational acceleration. Therefore wind particles are preferentially ejected along the rotation axis of a spinning object as the ‘axial wind’ in SH03 . We also decouple a new wind particle from hydrodynamic interactions for a brief time such that the winds originate from a region close to the surface of the star forming region. The full hydrodynamic interactions are enabled again once the density of the particle has fallen to 0.1 $\rho_{\rm th}$, or once a time of $\Delta t$ = 20 Myr has elapsed, whichever happens earlier.

2.8 Parameter setting and test simulations

In this paper we use the following parameters for star formation and SN feedback as a fiducial model; the threshold density, $\rho_{\rm th}=5\times 10^{-25}$ g cm*-3*, the dimensionless star formation efficiency, $c_{*}=0.005$, the dimensionless conductivity efficiency, $\eta_{\rm EVP}=0.1$, the SN feedback efficiency, $\eta_{\rm SN}=1$, and the wind velocity $v_{\rm w}=500$ km s*-1*. The relatively low value of $c_{*}$ compared to other papers comes from the fact that we use the dynamical time of each cloud not the dynamical time for SPH particles. We choose the values of $\rho_{\rm th}$ and $c_{*}$ so that our model reproduce the observed relation between surface gas density and surface star formation rate density (Kennicutt, 1998). The most important parameters in our model are the feedback efficiency and the wind velocity. The wind velocity we use is comparable to the strong wind model in Nagamine et al. (2004). For our choice of the parameters and the IMF, the wind production rate by SNe II is $\dot{M}_{\rm w}/\dot{M}_{*}\simeq 3$. There will be further contribution from SNe Ia at later times.

In this subsection, we show the behaviour of our model in idealised simulations of disk formation from virialised gas in the external NFW potential (Navarro et al., 1996). Here, $M_{\rm vir}$ is the virial mass of the system and we always assume that 10 per cent of the mass is in baryonic form. The initial angular momentum in terms of the spin parameter, $\lambda=J|E|^{1/2}/(GM_{\rm vir}^{2/5})$, is set to 0.07 with the assumption that the specific angular momenta $j(r)$ of spherical shells are all aligned and their magnitude is given by

[TABLE]

In Fig. 1 we show the star formation rates in the simulations for $M_{\rm vir}=10^{12}M_{\odot}$ with three different resolutions: $N_{\rm gas}=$ 10000, 50000, and 250000.

The star formation rates quickly converge. The simulation with $N_{\rm gas}=50000$, which is a comparable resolution that we will use in cosmological simulations, has very similar star formation history to that in the simulation with $N_{\rm gas}=250000$. Hereafter we thus only show the simulations with $N_{\rm gas}$ = 50000.

Now we show how star formation rates depend on the values of feedback parameters, $\eta_{\rm SN}$ and $v_{\rm w}$. We here run the simulations with two additional models: the slow wind model ($\eta_{\rm SN}=1$ and $v_{\rm w}=250$ km s*-1*) and the weak feedback model ($\eta_{\rm SN}=0.5$ and $v_{\rm w}=250$ km s*-1*). The wind speed in the slow wind model is half of that in the fiducial model, therefore it has an uncomfortably high mass loading factor $\eta_{\rm w}\simeq 12$. Of course the mass loading in the weak feedback model becomes half.

Fig. 2 shows the star formation rates in galaxies with these three wind models for halos of $M_{\rm vir}=10^{12}M_{\odot}$ and $10^{10}M_{\odot}$. Since $v_{\rm w}=500$ km s*-1* exceeds the escape velocity of the halo with $M_{\rm vir}=10^{12}\ M_{\odot}$ and $v_{\rm w}=250$ km s*-1* does not, the star formation is more strongly suppressed in the fiducial model in spite of the same feedback efficiency. We find that the star formation rates of the slow wind and weak feedback models have regularly time spaced peaks of activity. For the adopted NFW potential, the wind particles with the initial velocity, $v_{\rm w}=250$ km s*-1*, will return to the disk in $\sim 0.13$ Gyr. This time-scale might be responsible for the periodicity seen in the star formation rates. The star formation rate is highest in the weak feedback model for this halo as expected. However situation changes for the halo of $M_{\rm vir}=10^{10}\ M_{\odot}$. Now the mass loading factor is more important than the wind velocity because both wind velocities exceeds the escape velocity of this halo. Hence the star formation rate is highest in the fiducial model and lowest in the slow wind model although all models have star formation rates $\sim 10^{-3}\;M_{\odot}\,{\rm yr}^{-1}$ after 1 Gyr. Since we have to suppress star formation even in relatively large halos to form disk galaxies (OEFJ05, ), we employ the fiducial parameter set in our cosmological simulations.

In Fig. 3 we show the relations between the star formation rate per unit area and the surface gas density for the fiducial (top panel) and weak feedback (middle panel) models. We estimate the surface star formation rate densities by computing the surface density of stars younger than $3\times 10^{7}$ yr in cylindrical bins. Both models show reasonably good agreement with the target relation (Kennicutt, 1998, ; the target relation we show is eq. (25) of SH03 ). The use of different wind parameters do not affect this relation in our simulations, while the ISM can have higher gas density in the weak feedback model. For the galaxy in $M_{\rm vir}=10^{10}M_{\odot}$ halo, winds are so effective that there are only a handful of young stars in the galaxy. We hence only show, in the bottom panel of Fig. 3, the relation between the global surface gas density and the global surface star formation rate , which are averaged over the galaxy. The different points in this panel correspond to the relations at different epochs. We find that the global relation in the $10^{10}M_{\odot}$ halo is also in reasonable agreement with the target relation.

3 Black hole accretion and AGN feedback

3.1 Seed black holes

In the cosmological simulations we run a friends-of-friends halo-finder (Davis et al., 1985) 200 times on the fly from $z=15$ to [math]. When we find a halo that contains more than $N_{\rm th}$ dark matter particles and no BH particles, we put a seed BH on a density peak of the stellar distribution in the halo. We assume 260 $M_{\odot}$ as the mass of the seed black hole. This mass is consistent with the theory of stellar evolution, which shows that the nuclear burning in very massive stars above 260 $M_{\odot}$ is unable to halt gravitational collapse (e.g. Heger et al., 2003). Note that as far as the seed BH mass is sufficiently small, our results hardly change because we do not impose an upper limit on the mass accretion rate. Throughout this paper we use $N_{\rm th}=128$. The simulation employing $N_{\rm th}=256$ yields almost the same results except for the total number of BHs. The insensitivity on the choice of these parameters relies on our AGN feedback model that is only effective when the mass accretion rate onto the BH is much lower than the Eddington rate (see §3.3).

3.2 Mass accretion by the radiation drag

One mechanism yielding the proportionality of the BH mass to its host bulge mass has been proposed by Umemura (2001) and Kawakatu & Umemura (2002). According to Kawakatu & Umemura, in the central regions of galaxies the drag due to stellar radiation may result in a loss of angular momentum of a clumpy ISM and gives rise to mass inflow toward the centre. The optical depth of a gas cloud is $\bar{\tau}=\chi_{\rm d}\rho_{\rm c}r_{\rm c}\simeq\chi_{\rm d}m_{\rm c}/r_{\rm c}^{2}$, where $\rho_{\rm c}$, $m_{\rm c}$, and $r_{\rm c}$ are the density, mass, and a size of a cloud. In principle we can use the cloud model we described in §2.2. We however assume that the clouds are all identical and randomly distributed in a star forming region as in Kawakatu & Umemura for simplicity. The mass extinction coefficient $\chi_{\rm d}$ is given by $\chi_{\rm d}=300\ {\rm cm}^{2}{\rm g}^{-1}(a_{\rm d}/0.1{\rm\mu m}^{-1})(\rho_{\rm s}/{\rm g}\ {\rm cm}^{-3})(Z/0.3Z_{\odot})$, where $a_{\rm d}$ is the grain radius, $\rho_{\rm s}$ is the density of solid material within the grain (e.g. Spitzer, 1978), and $Z$ is the metallicity of gas. We use $a_{\rm d}=0.1\ \mu{\rm m}$, $\rho_{\rm s}=1$ g cm*-3*, and $Z$ is directly obtained in the simulations as the mean metallicity of the gas in a star forming region. The total optical depth of a star forming region can then be approximated as

[TABLE]

where $R_{\rm SFR}$ and $M_{\rm c}$ are the size of a star forming region and the total mass of the clouds in the star forming region.

The mass accretion rate due to the radiation drag mechanism described above is well approximated by

[TABLE]

where $L_{\rm SFR}$ is the total luminosity of stars in the star forming region. We calculate $L_{\rm SFR}$ by summing up the bolometric luminosities of all stars, each of which is obtained as a function of its age from a look-up table generated by P$\acute{\rm E}$GASE2 (Fioc & Rocca-Volmerange, 1997). We here ignore the metallicity dependence of the stellar luminosity and use a table for the solar metallicity because the dependence on metallicity is weak. Kawakatu & Umemura found that the efficiency $\eta_{\rm drag}$ is maximally $0.34$ in the optically thick regime assuming a spherical star forming region with a uniform density. In principle $\eta_{\rm drag}$ can be much larger or smaller depending on the geometry and density structure of the star forming region (Kawakatu private communication). We thus change the efficiency and adopt $\eta_{\rm drag}=1$ in this paper, considering more centrally concentrated density distribution in real star forming regions than assumed in Kawakatu & Umemura.

Once we define the radius of a star forming region around a BH, $R_{\rm SFR}$, we obtain $\chi_{\rm d}$, $M_{\rm c}$, and $L_{\rm SFR}$ in equations (18) and (19) directly in simulations, and hence we can compute $\tau_{\rm SFR}$ and $\dot{M}_{\rm drag}$ self-consistently; the effects of outflows are also implicitly incorporated in the estimate of these variables because we calculate them self-consistently at each time-step of the simulation. In order to define the size of the star forming region, $R_{\rm SFR}$ we start from a sphere around a BH that contains 40 SPH particles, and then we search the radius of the sphere that maximises $\dot{M}_{\rm drag}$ by increasing the radius of the sphere. We employ this radius as the size of the star forming region. When the mean gas density of the initial sphere is lower than the threshold density for star formation, $\rho_{\rm th}$, we consider that there is not enough gas around the BH to fuel accretion and we simply set $\dot{M}_{\rm drag}=0$. Kawakatu & Umemura (2004) argued that radiation from bulge stars contributes to the mass accretion to the central hole more efficiently than that from disc stars. We crudely mimic this effect by using spherically averaged values to define star forming region; by doing this we put less weight to the star forming gas that has a flattened distribution. We will show how star formation in the disc and bulge correlate with the mass accretion rate, $\dot{M}_{\rm drag}$, in our simulations later.

Strictly speaking, one should distinguish $\dot{M}_{\rm drag}$ from the accretion rate onto the BH, $\dot{M}$, because the radiation drag is not likely to remove the angular momentum thoroughly; some residual angular momentum will terminate the radial contraction of the accreted gas (Sato et al., 2004). Granato et al. (2004) introduced the viscous timescale in which the mass accumulated at the galactic centre by the radiation drag is fed to the BH. We however assume $\dot{M}=\dot{M}_{\rm drag}$ for simplicity. The time delay in mass fuelling introduced by the viscous timescale may become important when we compare the luminosity of AGNs and their host galaxies in detail.

3.3 AGN feedback

Motivated by the success of recent semi-analytic models (Croton et al., 2006; Bower et al., 2006), we assume that AGNs directly heat diffuse hot gas confined in dark halos through the production of jets. As the jet production mechanism, we employ the model proposed by Meier (2001) which is an hybrid version of the Blandford-Znajek (Blandford & Znajek, 1977) and the Blandford-Payne (Blandford & Payne, 1982) processes. In this hybrid model the jet power is generated by the magnetic field threading the accretion flow inside and outside the the ergosphere of the BH. This model is able to draw upon the rotational energy of the accretion flow as well as the spinning black hole in order to drive jets, and its basic features are supported by recent relativistic numerical simulations of jet formation (e.g. Hawley & Krolik, 2006b). We refer the readers to Meier (2001) and Nemmen et al. (2007) for more details.

Meier (2001) presented models both for non-spinning (Schwarzschild) and for spinning (Kerr) BHs. Semi-analytic cosmological simulations of the spin evolution of BHs through mergers and gas accretion (Volonteri et al., 2005) and estimates of the radiative efficiencies of global populations of quasars based on Soltan-type arguments (e.g. Soltan, 1982; Yu & Tremaine, 2002; Wang et al., 2006) suggest that most nearby massive BHs are rapidly rotating. Recently Nemmen et al. (2007) also reported that BHs in nearby elliptical galaxies are likely to have high spins $j\simeq 0.7-1$. We thus only consider spinning BHs in this paper. Meier (2001) parameterised the jet luminosity for standard thin discs as

[TABLE]

and for RIAFs as

[TABLE]

where $\alpha_{\rm SD}$ and $\alpha_{\rm RIAF}$ are the viscosity parameters for standard thin discs and for RIAFs, respectively, $m_{9}$ is the black hole mass in units of $10^{9}M_{\odot}$, $\dot{m}\equiv\dot{M}/\dot{M}_{\rm Edd}$ is the accretion rate scaled to the Eddington rate assuming the 10% energy conversion efficiency ($\dot{M}_{\rm Edd}\equiv 22\ m_{9}\ M_{\odot}\ {\rm yr}^{-1}$), and $f$ and $g$ are the ratios of the actual angular velocity and the azimuthal magnetic field to those calculated by Narayan & Yi (1995) respectively. According to Meier (2001) we employ $f=1$ and $g=2.3$. Since the limited numerical resolution of our cosmological simulations prevents us from resolving the details of the propagation of jets and its impact on the surrounding medium, we simply assume that a fraction of the jet energy $\Delta E_{\rm FB}^{\rm jet}=\eta_{\rm AGN}L_{\rm jet}\Delta t$ is delivered as thermal energy to the neighbouring 40 diffuse gas particles that have density, $\rho<0.1\rho_{\rm th}$. In our simulations, the size of the region heated by the AGN feedback is typically $\sim 10$ kpc. Unless otherwise stated, we use $\alpha_{\rm SD}=\alpha_{\rm RIAF}=0.1$, $j=0.5$, and $\eta_{\rm AGN}=0.1$. The adopted value of $j$ corresponds to the intermediate case between Schwarzschild (no-spin) and maximally rotating BHs.

As previously mentioned, we adopt two distinct regimes of accretion flows: standard thin discs (optically thick, geometrically thin, radiatively efficient) and RIAFs (optically thin, geometrically thick, radiatively inefficient). The parameter which controls the state of the accretion flow is $\dot{m}$, and the critical value of the accretion rate which sets the division between the two regimes is given by $\dot{m}_{\rm crit}\approx\alpha^{2}$ (e.g. Narayan et al., 1998). Since the RIAF solution ceases to be valid above $\dot{m}_{\rm crit}$, for $\dot{m}\leq\dot{m}_{\rm crit}$ a RIAF exists and for $\dot{m}>\dot{m}_{\rm crit}$ a thin disc occurs. As a consequence of the adopted theoretical model of jet production, we have strong AGN feedback only for RIAFs and the maximum AGN feedback occurs at $\dot{m}=\dot{m}_{\rm crit}$. For the parameters we chose, we have the following jet production efficiency in RIAFs,

[TABLE]

which is much higher than that in standard discs, given by

[TABLE]

where in equation (23) we adopt $m_{9}=\dot{m}=1$. The equations (22) and (21) correspond to the “radio-loud mode” in our simulations, and the equations (23) and (20) represent the “radio-quiet mode”.

Note that, in principle, we could have even higher efficiencies of jet production, because the energy is drawn from the “spin energy” stored in the BH. We explore the case of a higher efficiency corresponding to a higher BH spin later in §5.3.

4 Cosmological simulations of disc galaxy formation with AGN feedback

To study the role of AGN feedback in disk galaxy formation, we use the initial conditions presented in OEFJ05 , from which they produced an extended disk galaxy. The background cosmology is the so-called $\Lambda$CDM with the mean matter density $\Omega_{0}=0.3$, Hubble parameter, $h\equiv H_{0}/100$ km s*-1Mpc${}^{-1}=0.7$, cosmological constant term, $\Omega_{\Lambda}\equiv\Lambda_{0}/(3H_{0}^{2})=0.7$, amplitude of mass fluctuations, $\sigma_{8}=0.9$, and mean baryon density, $\Omega_{\rm b}=0.04$. The size of the periodic simulation box is $35.325h^{-1}$Mpc and the initial redshift is $z=49$. We put the high-resolution dark matter particles and gas particles in the region where a halo with mass $M_{\rm vir}=1.2\times 10^{12}h^{-1}M_{\odot}$ forms at $z=0$. The region external to this is populated with high-mass dark matter particles (boundary particles), the function of which is to reproduce the appropriate tidal fields. The circular velocity, spin parameter and collapse redshift of the selected halo are $v_{\rm c}(r_{\rm vir})=155$ km s-1*, $\lambda\equiv J|E|^{1/2}/(GM^{5/2})=0.038$ and $z_{\rm c}\simeq 1.5$, respectively, where $z_{\rm c}$ is defined as the redshift at which the main progenitor had half the final halo mass.

The masses of the SPH and high-resolution dark matter particles are $\sim 2.6\times 10^{6}$ and $\sim 1.7\times 10^{7}h^{-1}M_{\odot}$, respectively. Thus the $N_{\rm th}=128$ implies that halos larger than $\sim 2.6\times 10^{9}h^{-1}M_{\odot}$ are allowed to have BHs. The gravitational softening length are kept fixed in comoving coordinates for $z>3$; thereafter they are frozen in physical units at a value (of the equivalent Plummer softening) of $\epsilon=0.5$ and $1$ kpc for the baryonic particles (SPH particles, stars, and BHs) and high-resolution dark matter particles, respectively. The gravitational force obeys the exact $r^{-2}$ law at $r>2.8\epsilon$.

4.1 Black hole particles

When we find a halo with $N_{\rm dm}\geq N_{\rm th}$ that does not contain any BH particles and is not contaminated by the boundary particles, we turn the star particle that has the largest stellar density into a seed BH. We thus discriminate the BH mass from the mass of the BH particle (particle mass). The former is updated every time-step according to the mass accretion rate obtained from the radiation drag model. Once BH mass exceeds its particle mass, the BH particle increases its particle mass by absorbing gas particles $n$ in the star-forming region around it with the probability:

[TABLE]

where $M_{\rm BH}$ and $M_{\rm p}$ are, respectively, the BH mass and the mass of the BH particle, $N_{\rm SFR}$ is the number of gas particles in the star-forming region, and $M_{n}$ is the mass of the gas particle $n$.

BH particles can increase their BH masses and particle masses by BH mergers. However, how long it would take to harden a BH binary until finally gravitational wave emission becomes important is still unclear (Makino & Funato, 2004; Escala et al., 2004). We simply assume that BHs can merge into a single BH if two or more BHs are within a gravitational softening radius and they are gravitationally bound. Therefore the BH merger rates obtained in our simulations are likely to be the maximum possible rate. Technically, we apply the friends-of-friends algorithm to find groups of BHs which may potentially merge originating new BHs.

An interesting problem arises owing to the finite numerical resolution. In reality, supermassive BHs are much heavier than stars and dark matter. Therefore they are quickly brought to the galactic centres and stay at the bottom of the potential wells owing to the dynamical friction. In cosmological simulations, the mass of a BH particle is comparable to those of gas and stellar particles and lighter than the dark matter particles in the early stage of its evolution. As a result, BH particles tend to oscillate around galactic centres with the velocity dispersions of the stellar systems. To circumvent this problem, we allow BH particles to behave as peak tracers at galactic centres. At each time-step, we compute local stellar density fields within the gravitational softening radii of BH particles, $\epsilon$. We then displace the BH particles by a small distance, $\Delta l$, along the local density gradients. After several tests, we chose

[TABLE]

where $v$ is the velocity of a BH particle and $\Delta t$ is the time-step. This $\Delta l$ correction sufficiently suppresses the oscillation of the BH particles at galactic centres.

We adopt a maximum time step for BHs, $\Delta t_{\rm max}=1$ Myr, in order to follow mass accretion and associated feedback with a reasonable time resolution. Hence a time step for a BH particle is defined as $\Delta t_{\rm BH}=\min(\Delta t_{\rm max},\Delta t_{\rm dyn})$, where $\Delta t_{\rm dyn}$ is a time step determined by the dynamics.

4.2 Models

We first compare three simulations in order to show the effects of AGN feedback. The first simulation does not have feedback from AGN. We refer this model as the ‘no-fb’ model. In the second model, we distribute 1% of jet energy to the ambient halo gas particles defined by their densities, $\rho<0.1\rho_{\rm th}$ in the form of thermal energy using SPH smoothing. Since only 1% of jet energy is available for heating, we call this model the ‘weak-fb’ model. In the ‘reference’ model, we adopt the same AGN feedback scheme as in the weak-fb model but assume 10% of jet energy is available to heat the ambient halo gas. We will explore other models in which we change the implementation of AGN feedback and strength of SN feedback in §5.3.

5 Results

To give a visual impression of our simulations, we show some snapshots of the simulation with AGN-radio feedback (‘reference’ model) in Fig. 4. There are a number of merger events at high redshift ($z>1$). The disk formation starts at $z\sim 1$, and finally a disk galaxy forms. A lot of BHs form in the small halos and they are brought to the main progenitor by merging satellites. The BH particles nicely trace the centres of the stellar objects resulting from our correction applied to the positions of BH particles. There are also many BHs which are not associated to any stellar objects. They are scattered into the halo when small satellites, which harbour a small BH at their centres, are destroyed during mergers. Most of the drifting BHs hence have very small mass.

In the following sections,we present the main results of the paper. In §5.1, we compare the results of simulations run with the standard feedback parameters (the ‘reference simulation’) with simulations using weaker AGN feedback (the ‘weak-fb’ model) and no AGN feedback at all (the ‘no-fb’ model). In §5.2 we look at the importance of the AGN feedback for the bolometric luminosity of the forming galaxy, while in §5.3, we explore the dependence of the results on the way in which the feedback is coupled to the surrounding gas and the impact of assuming a lower efficiency of the feedback from SNe.

5.1 Simulations with different AGN feedback efficiencies

Now we compare galaxies formed in three simulations, i.e. the ‘no-fb’, ‘weak-fb’, and ‘reference’ simulations. We first show the stellar distributions within $50h^{-1}$ kpc boxes centred on the galaxies at $z=0$ in Fig. 5. The edge-on and face-on projections are selected to be perpendicular and parallel to the angular momentum vector of the stellar component within the central $10h^{-1}$ kpc sphere. In all of the models, the galaxy has an extended stellar disk. The ‘no-fb’ and ‘weak-fb’ models are very similar. From the face-on view, it is found that the ‘reference’ galaxy has a disk with lower surface stellar density.

Fig. 6 shows the $I$ band surface brightness profiles for these galaxies. We also show the surface brightness profile of a disk galaxy obtained by OEFJ05 using the same initial condition. An interesting difference from the surface brightness profile of the galaxy by OEFJ05 can be seen at the centre. Our galaxies do not show the extreme central concentration. The origin of this difference is the different implementation of SN feedback. While we implement it as kinetic winds, they deposited feedback energy in the form of thermal energy. Therefore they had to calculate the pressure gradients caused by the thermal feedback to blow the gas away. Unfortunately, the SPH does not allow steep gradients in physical quantities by definition unless one uses a sufficiently high resolution.

All our three galaxies show similar central surface brightness. The main difference due to AGN feedback can be seen in the surface brightness of the disks. While ‘no-fb’ and ‘weak-fb’ galaxies show the almost same surface brightness profiles, the ‘reference’ model has a less extended disk. We do not try to determine the bulge-disk decomposition based on the surface brightness profiles since they provide only partial information about the morphology of galaxies. Below, we will show the decomposition based on the dynamics of stars (Abadi et al., 2003). This is a more accurate way to characterise the relative importance of the bulge and disk components.

To investigate effects of AGN feedback quantitatively, we present the star formation histories of the stars that lie within $25h^{-1}$ kpc from the galactic centres at $z=0$ in Fig. 7. All galaxies have qualitatively similar star formation histories and have peak star formation rates at $z\sim 1.5$. During the starburst ($t=2.5-6$ Gyr), the ‘reference’ galaxy shows lower peak star formation rate compared with other galaxies. As we will show later, the AGN is not a RIAF during the starburst. Hence AGN feedback from a RIAF is not responsible for this lower star formation rate. The feedback from standard disk (equation (23)) might suppress the starburst. However we should not overinterpret the differences during the starburst because the system is strongly non-linear during this epoch and hence small changes can result in relatively large difference. In fact, the model which employs stronger AGN feedback has a higher star formation rate during the starburst than in the ‘reference’ model as we will show later in Fig. 15.

After the starburst, the ‘reference’ galaxy has significantly lower star formation rate than other galaxies. This suppression of star formation is likely to be caused by AGN feedback from a RIAF. The lower star formation rate in the ‘reference’ galaxy explains its less extended disk in Fig. 6 because it is the star formation after $z\sim 1$ that builds up the stellar disk (see Fig. 4). The ‘weak-fb’ galaxy has almost the same star formation rate as the ‘no-fb’ galaxy. Hence it can be said that, in order to give a visible impact on galaxy formation, at least $\sim 10\%$ of jet energy has to be given to the halo gas. We will later explore the case in which the jet energy is not given to the diffuse halo gas but delivered to the surrounding ISM.

Since the strength of AGN feedback in our model is closely related to the central BH mass and the accretion rate onto the BH, we plot the evolution of the central BH of the main progenitor and the mass accretion rate in units of the Eddington rate in Fig. 8 in each model. Note that the mass accretion rate $\dot{M}$ is different from the mass growth rate, $\dot{M}_{\rm BH}$; the latter includes the mass increase due to BH mergers. In all simulations, the central BHs rapidly evolve until $t\sim 5.5$ Gyr at which the starburst ceases. The final BH masses reach $\sim 10^{7}M_{\odot}$. The stepwise jumps seen in the evolution of the BH mass indicate the mass growth by mergers.

The consumption of the available gas for star formation explains the decline of the accretion rate onto the central BH. At $t\simeq 5.5$ Gyr $\dot{m}$ reaches the critical accretion rate $\dot{m}_{\rm crit}$, below which the accretion disk becomes a RIAF. Our model has the maximum AGN feedback when the mass accretion rate is equal to the critical accretion rate. This epoch is indicated by the crossing of the solid line and the horizontal dashed line in the lower panel of Fig. 8. Once the accretion disk becomes a RIAF after $z<1$ , in the ‘reference’ model, the accretion rate drops to zero rapidly and there is no mass accretion after the strongest production of AGN feedback at $t\sim 5.5$ Gyr except for some episodic accretion events. This implies that there is almost no star formation activity at the centre of the main progenitor at $z<1$ for this galaxy. On the other hand, lacking the strong AGN feedback, there are not such sudden decrease in accretion rates in the ‘no-fb’ and ‘weak-fb’ models. The gradual decline is caused by consumption of gas due to star formation and associated SN feedback(winds).

Our implementation of AGN feedback cannot directly suppress the star formation because it injects energy not into the star forming dense gas but into the diffuse halo gas ($\rho<0.1\rho_{\rm th}$). Thus it is the suppression of the cooling of the halo gas around the centre of the galaxy which is responsible for the halting of star formation around the BH in the ‘reference’ simulation. The impulsive feature seen in the accretion rate implies the existence of a self-regulated cycle that operates between gas cooling, star formation, and production of jets in the RIAF. The fact that the central surface brightness is almost the same in the three models in spite that the galaxy with AGN feedback hardly form stars at its centre after $z\sim 1$ suggests that the vast majority of stars at the centre have already been formed at $z\sim 1$ in all galaxies.

We now return to the separation of the disk and bulge components. As noted earlier, it is best to separate the bulge and disk components using a dynamical decomposition (Abadi et al., 2003). For this, we first compute the angular momentum $J_{z}$ of each star particle parallel to the net angular momentum of stars within $10h^{-1}$ kpc, and the angular momentum of the corotating circular orbit, $J_{c}(E)$, where $E$ is the energy of each particle. The ratio $J_{z}/J_{c}$ defines an orbital circularity.

In Fig. 9, we show the probability distribution of this orbital circularity for stars within $25h^{-1}$ kpc from the centre of each galaxy. A disk component should have $J_{z}/J_{c}(E)\simeq 1$ and we see that such a component clearly dominates in each galaxy. Comparing to the results by OEFJ05 , the suppression of the formation of spheroidal components seems to rely largely on the strong SN feedback by winds. Governato et al. (2007) also produced disk dominated galaxies in halos with quieter merger histories (in which the epoch of the last major merger is $z_{\rm lmm}>2$) using a different feedback recipe from ours. To identify the stars in the spheroid, we assume a non-rotating spheroid, i.e. that stars in the spheroid are symmetrically distributed around zero in Fig. 9, where all stars having $J_{z}/J_{c}(E)\leq 0$ are identified as a half of the spheroidal component. Another half of the spheroid with $J_{z}/J_{c}(E)>0$ are defined statistically so that the total angular momentum of the spheroid becomes zero. All remaining stars are identified as the disk component. As expected from the similarity in the star formation histories, all galaxies have similar distributions in this plane while the disc component in the ‘reference’ model is slightly less significant because of the relatively strong AGN feedback.

In Table 1 we show the total and spheroid’s masses, total luminosity, and the disk-to-total ratios in mass, $g$ band, and $r$ band. Since colour of a stellar population is sensitive to its metallicity, we here include the metallicity dependence to compute the luminosity at each passband by using P$\acute{\rm E}$GASE2 (Fioc & Rocca-Volmerange, 1997). The masses of the spheroidal components directly reflect amplitude of starbursts and therefore the ‘no-fb’ galaxy has the most massive spheroidal component. Nevertheless, all galaxies have disk-dominated morphology.

It is interesting to see when each component forms and how the formation history of each component is related to the accretion history of the central BH. In Fig. 10, we show star formation histories of the disk and spheroidal components in each galaxy at $z=0$. We also show the rescaled accretion rate onto the central BH of the main progenitor of each galaxy. Note that the star formation rates include star formation in all building blocks (all progenitors) but the mass accretion rates are only the rates onto the central BHs of the main progenitors. Thus at high redshift, the mass accretion rates presented in the figure considerably underestimate the net accretion rates. Since the galaxies do not undergo any significant mergers after $z\sim 1$, it is safe to compare the star formation rates and the mass accretion rates at low redshift. We find that there are good correlations between the star formation rates of the spheroid and the mass accretion rates onto the central BHs at low redshift, in particular between $t\sim 5.5$ and 8 Gyr. After the starburst, the star formation rates of the spheroids rapidly decreases in all simulations. The accretion rates onto the BHs also decrease at the same time. The contribution from disk stars to the mass accretion rate due to the radiation drag effect should be largest during the peak of star formation in the disc component ($t\simeq 5$ Gyr) because at this epoch the disc is gas rich and very compact. We have however confirmed that even at this epoch the contribution from disc stars is less than 20% and most of the time less than 5% in all models. Considering the fact that half of the spheroid stars are defined statistically, the correlation between the star formation rate and the mass accretion rate of the spheroid in each galaxy is significant. After $t\sim 8$ Gyr there are many short lived star formation episodes in the spheroids, which do not correspond to gas accretion events. These star formation events occur in small satellites before they become a part of the spheroidal components by minor mergers.

Next we show in Fig. 11 the evolution of the relations between mass of the central BHs of the main progenitors and the mass of the host spheroids at $z=3$, 2, 1.5, 1, 0.5, 0.2, and 0. We here define the main progenitors as stars within 10% of the virial radii at each redshift. We have confirmed that the spheres with radii of $0.1r_{\rm vir}$ contain the stellar systems defined as main progenitors by the friends-of-friends algorithm in §5.2. We have also checked by eye that at selected redshifts the main progenitors are not undergoing major mergers and the spheres do not contain massive satellites, both of which ruin the dynamical decomposition. We then compute the mass of the spheroidal component of each main progenitor based on the dynamical decomposition, noting that this method of decomposition is less vulnerable to insufficient resolution compared with the decomposition based on surface density profiles. This feature is important at high redshift where the size of the galaxies are much smaller than that at $z=0$. Fig. 11 shows that the $M_{\rm BH}$-$M_{\rm SP}$ relations evolve similarly in all simulations because the BHs and spheroids mainly increase their mass during the starbursts and our AGN feedback is not effective during this phase. Before the starburst, the BH mass is far below the local relation, $M_{\rm BH}\sim 10^{-3}M_{\rm SP}$, and then it rapidly moves towards the local relation during the starburst ($z=2\sim 1$). The delay in the mass evolution of the BHs compared to the spheroids is resulting from the fact that the mass accretion rate due to the radiation drag model depends on the metallicity of the ISM. The mass accretion by this model is not efficient until the ISM is metal enriched and the optical depth becomes $\sim 1$. After the starburst, the mass of BHs stays almost constant as seen in Fig. 8. The relation is imprinted in the radiation drag model because the model provides a good correlation between star formation in a spheroid and the mass accretion onto a BH. It is important to note that the radiation drag model predicts considerably lower mass ratio between central BHs and their host spheroids before the starburst.

In order to illustrate how the colours of our galaxies evolve, we show the evolutions of the main progenitor of each galaxy and its spheroid in a colour-magnitude diagram (Fig. 12). The mean colours of the red and blue populations of the SDSS galaxies (Baldry et al., 2004) are indicated by solid lines. All galaxies show similar evolution on the colour-magnitude plane. This confirms that the AGN feedback has only minor effects on evolution of disk galaxies. The colour of the galaxies and their spheroids at $z=0$ are consistent with the observed colours of local blue and red galaxies, respectively.

5.2 Relative importance of the AGN

In order to show the relative importance of the AGN in galaxy formation, we show the bolometric luminosities generated by stars in the main progenitors and the bolometric luminosities radiated by the AGNs of the main progenitors in the upper panels of Fig. 13 for our three simulated galaxies. A main progenitor is defined by applying the friends-of-friends algorithm for stellar particles with a linking length $l_{\rm link}=0.02\left\langle l\right\rangle$, where $\left\langle l\right\rangle$ is the the mean separation of dark matter particles. To estimate the luminosity of an AGN, we assume that the efficiency of conversion of rest mass energy of the accretion flow into radiation ($L_{\rm AGN}/(\dot{M}c^{2})$) is $\eta_{\rm SD}=10\%$ for a standard thin disk and for a RIAF we use the equation

[TABLE]

where $\delta$ represents the fraction of the turbulent energy which heats the electrons (Quataert, 2001). While the value of $\delta$ is quite uncertain, we employ $\delta=0.1$ as commonly assumed in spectral models of RIAFs (Nemmen et al., 2006). Note that the radiative efficiency for a thin disk can be a factor of a few higher than $0.1$ in the case of a spinning BH, and may reach efficiencies as high as $42\%$ for an extreme Kerr BH (Novikov & Thorne, 1973).

It is found that the main progenitors are always brighter than the AGNs for these galaxies as we expect for disk-dominated galaxies. Since the radiation drag is more efficient for more compact star forming region, we expect that the significance of the AGN luminosity would increase for galaxies which are undergoing gas-rich mergers. There is a notable difference among our models when the accretion flows become RIAFs after $z\sim 1$. While all the AGNs have much lower luminosities than their host galaxies at low redshift because of the low accretion rates and the low radiative efficiency of the RIAF, only the AGN in the ‘reference’ galaxy has almost zero luminosity right after the accretion disc switches to a RIAF. In this simulation, the low accretion rate, due to the suppression of star formation around the BH by AGN feedback, makes the AGN invisible most of the time. In the ‘no-fb’ and ‘weak-fb’ simulations, the AGNs still continuously emit light for a few giga years after the onset of the RIAFs. This result suggests that if 10% of the jet luminosity is used to heat the ambient halo gas AGNs become quite dim and can be hardly observed after the starburst regime. It should be noted that it is beyond the scope of this work to incorporate in the simulations possible inclination effects in the appearance of the AGN, caused by dust obscuration inside the nucleus as in the obscuring dust torus model invoked in unified models of AGN (e.g. Antonucci 1993; Krolik 1999).

In the lower panels of Fig. 13 we show the raw feedback luminosities ejected by the total SNe in the main progenitors and jet luminosities produced by the AGNs computed by using equations (20) and (21). Until the end of the starburst, the feedback energy is dominated by energy produced by SNe except for a few temporary RIAF modes seen in the ‘weak-fb’ simulation. The jet luminosities suddenly exceed the SN luminosities when the accretion flow becomes radiatively inefficient. Since in the ‘no-fb’ and ‘weak-fb’ simulations, the efficiencies, $\eta_{\rm AGN}$ by which the AGN jet energy is given to the surrounding halo gas are small ($\eta_{\rm AGN}=0$ and 0.01 respectively), the effects of the AGN feedback are none or negligible in these simulations. In the ‘reference’ simulation, the efficiency, $\eta_{\rm AGN}=0.1$ is sufficiently high to suppress the star formation around the BH therefore quenching the mass accretion to the BH. This is an important feature of our AGN feedback. The impulsive events in Fig. 13 represent epochs where the AGN would be seen as radio-loud, and they constitute the radio-loud mode of the AGN. The duration in which the AGN is in the radio-loud mode is strongly dependent on the feedback efficiency, since the AGN feedback terminates the radio-loud mode through the suppression of gas cooling and hence star formation. Interestingly enough, in the ‘reference’ simulation, the AGN spends most of the time in the radio-quiet mode and only a fraction 8.8% of the time in the radio-loud mode. Contrarily, in the ‘no-fb’ and ‘weak-fb’ galaxies, the AGNs spend $\sim 25\%$ of the time in the radio-loud mode. Note that we are likely to overestimate the star formation rate at the centre and hence AGN activity particularly in the quiescent star forming regime because of the angular momentum transfer due to the artificial and numerical viscosities which can bring gas particles to the centre (Okamoto et al., 2003; Kaufmann et al., 2007).

5.3 Dependence on feedback parameters and implementations of feedback

In §5.1, we compared simulations with different AGN feedback efficiencies. In this sections, we compare simulations which all have strong AGN feedback, but which differ in other aspects. We compare three additional models. The models investigate the coupling of the AGN feedback, the role of BH spin and the effect of the wind speed generated by SN feedback.

In the first model, AGN feedback energy is not given to the ambient halo gas (diffuse gas) but deposited to the nearest 40 ISM particles (dense gas). We use $\eta_{\rm AGN}=1$ so that all the jet energy is used to heat the surrounding ISM because a model with $\eta_{\rm AGN}=0.1$ does not show any visible difference from the ‘no-fb’ model. We call this model the ‘ism-fb’ model. In the second model, we assume the higher spin, $j=0.9$, and the maximum feedback efficiency, $\eta_{\rm AGN}=1$, as we may need such strong feedback in order to explain the X-ray properties of clusters of galaxies (e.g. Nemmen et al., 2007). As a result, the efficiency of AGN feedback in this model becomes 17 times higher than in the ‘reference’ model . We dub this model the ‘max-fb’ model. In the third model, we use the same AGN feedback as in the ‘reference’ model but employ a slow wind speed, 250 km s*-1* for SN feedback in order to explore effects of the wind speed. We refer this model as the ‘slow-wind’ model.

In Fig. 14, we show the edge-on and face-on views of the stellar distributions of the galaxies produced by the three additional models. The ‘ism-fb’ model is similar to the ‘reference’ model. As expected the disk component of the ‘max-fb’ galaxy is less significant than that of the ‘reference’ model since stronger AGN feedback suppresses the disk formation more strongly. The ‘slow-wind’ model exhibits completely different morphology from others. There is no sign of a disk component.

Now we carry out quantitative comparison by showing the star formation histories of the galaxies and the accretion histories to the central BHs (Fig. 15). In spite of the fact that we employ an order of magnitude higher AGN feedback efficiency in ‘ism-fb’ model, the star formation history of this galaxy is very similar to that of the ‘reference’ model. It suggests that delivering AGN feedback energy to the diffuse hot gas causes the stronger feedback effect than depositing it into the surrounding dense ISM where cooling time is very short. The ‘max-fb’ galaxy shows a resembling star formation history to the ‘reference’ model (upper-middle panel) at high redshift ($t<6$ Gyr) where the AGN harbours a thin accretion disc. Once the accretion flow becomes a RIAF, strong AGN feedback suppresses the star formation dramatically and this galaxy have far less star formation after $t\sim 6$ Gyr than the ‘reference’ galaxy. The ‘slow-wind’ galaxy has a completely different star formation history from others. Owing to the slow wind speed, it has a higher star formation rate at high redshift including two burst events at $t\sim 2$ and $3$ Gyr. The AGN in this galaxy becomes underfed earlier than in other models (at $t\sim 4$ Gyr) and the AGN feedback strongly suppresses the subsequent star formation. We also show the model which employs the same wind speed as the ‘slow-wind’ model but without AGN feedback (dashed line in the upper-right panel). By comparing these two slow-wind models, we can see that the AGN feedback is responsible for the suppression of the star formation at low redshift. The bulge-to-total mass ratio of the ‘slow-wind’ model is 0.9 and the $u-r$ colour of this galaxy is consistent with the mean colour of the red population at $z=0$.

The evolution of the accretion rates in the ‘ism-fb’ and ‘max-fb’ models are similar to that in the ‘reference’ model particularly at high redshift ($t<6$ Gyr) where the accretion disc is thin ($\dot{m}>10^{-2}$). In these two models, the accretion flows become RIAFs after the starbursts as in the ‘reference’ model. After that the mass accretion onto the central BH in the ‘ism-fb’ model is more strongly suppressed than in the ‘reference’ model, while the ‘ism-fb’ model has a higher star formation rate. Since the AGN feedback energy is deposited into the surrounding star forming ISM particles in this model, it can terminate star formation around the BH more directly than in the ‘reference’ model where the AGN feedback energy is distributed to the halo gas. In the max-fb model, there is no accretion event after $t\simeq 6$ Gyr. The strong AGN feedback in this model removes all the low angular momentum gas around the BH and prohibits subsequent star formation around the galactic centre. The ‘slow-wind’ model shows a different accretion history from other models. Since the starbursts occur earlier in this model owing to the slow wind speed, the AGN also becomes underfed earlier ($t\sim 4$ Gyr) than in other models. Contrary to other models, the mass accretion continues at $t\leq 4$ Gyr. This is because the slow wind speed, $250$ km $s^{-1}$, does not allow gas particles to escape from galactic centre and the AGN feedback energy is not given to the star forming gas. The accretion rate might become more episodic if we deposit a significant fraction of the AGN feedback energy to the central star forming gas as in the ‘ism-fb’ model.

In Fig. 16, we show the relations between mass of the central BHs of the main progenitors and the mass of the host spheroids at $z=3$, 2, 1.5, 1, 0.5, 0.2 and 0. The relations for ‘ism-fb’ and ‘max-fb’ models evolve similarly to the relations for the models presented in Fig. 11. The relation for the ‘slow-wind’ model reaches the local relation earlier (at $z=2$) because this galaxy undergoes starbursts earlier than other galaxies. At low redshift, this model has slightly smaller BH mass to spheroid mass ratio than in other models. This is because infalling satellites increase the mass of the spheroid without changing the mass of the central BH very much when these infalling galaxies are disk dominated. Due to the slow wind speed, these infalling systems have more stars than in other models. Fig. 11 and Fig. 16 tell us an important feature of the radiation drag model, that is, the $M_{\rm BH}-M_{\rm SP}$ relation after the last starburst does not depend on either the wind speed or efficiency of AGN feedback. In all models, the central BH masses are $\sim 0.1$% of the mass of their host spheroids, $M_{\rm BH}\simeq 10^{-3}M_{\rm SP}$, when the accretion becomes radiatively inefficient.

6 Comparison with observations of AGN in spiral galaxies

As we have shown in §5.2, for $z\lesssim 1$ the bolometric luminosity of the AGN is feeble compared to the host galaxy and the nucleus undergoes several recurring episodes of activity. During these episodes when the nucleus is active, its characteristic accretion rate is in the range $\dot{m}\approx 10^{-4}-10^{-3}$, while the typical nuclear bolometric luminosities and jet powers are $L_{\rm bol}\approx 10^{40}-10^{42}\;{\rm erg}\ {\rm s}^{-1}$ and $L_{\rm jet}\approx 10^{42}-10^{43}\;{\rm erg}\ {\rm s}^{-1}$ respectively. Given these properties, at $z\lesssim 1$ during which the galactic disk is formed, the AGN can be considered to be in a highly sub-Eddington, low-luminosity or Seyfert state which undergoes a stochastic cycle of activity. The Seyfert nucleus is sometimes turned on, but most of the time is “off” or inactive with $\dot{m}\approx 0$, resembling the nucleus of a normal spiral galaxy. The “on” episodes are induced by the inward motion of clumps of cold gas into the centre of the galaxy, which provide a replenishment of new stars. The radiation drag or Poynting-Robertson effect caused by stellar radiation in the central star-forming region then dissipates the angular momentum of the gas and fuels the AGN. Our model suggests that reasonably powerful jet outbursts accompany each “on” episode of the low-luminosity AGN, suggesting that, in analogy to the observed states of AGNs, during the active phases the Seyfert nucleus would be seen as radio-loud, becoming radio-quiet during the “off” phase. In this section, we compare the properties of our simulated AGN with those of observed ones in disk galaxies.

In the reference model, the range of values obtained for the AGN bolometric luminosity for $z\lesssim 1$ agrees well with the median observational values for LINERs (median $L_{\rm bol}\approx 10^{41}\;{\rm erg}\ {\rm s}^{-1}$) and Seyferts (median $L_{\rm bol}\approx 10^{42}\;{\rm erg}\ {\rm s}^{-1}$) inferred from the Palomar spectroscopic survey of nearby galaxies (Ho, 2004). Since the Seyfert episodes of the central AGN coincide with jet outburst events during which the AGN would be seen as radio-loud in our simulations, this suggests that observed Seyfert/LINER nuclei would be mostly radio-loud objects rather than radio-quiet given their formation history, contrary to the generally held notion existent in the literature (e.g., Krolik 1999; Laor 2000). Ho & Ulvestad (2001) and Terashima & Wilson (2003) however investigated several low-luminosity AGN and their radio-loudness classification, obtaining that a large fraction of Seyfert/LINER nuclei are radio-loud, a finding that corroborates our results. In particular, Ho & Ulvestad (2001) obtained that $\gtrsim 60\%$ of the sources in their sample are radio-loud.

In our simulations, the numerical resolution is not sufficient to resolve the details of the propagation of the jet as it flows out of the AGN and how it transfer its power to the surrounding medium, for instance the spatial resolution is 0.5 kpc. As such, we assume that a fraction of the jet power is injected as thermal energy in the diffuse gas around the AGN on kpc-scales. The size of the region which is affected by the feedback from the AGN may be as large as $\approx 10-15$ kpc, although the jet itself is not resolved in the simulation. Two important questions arise: (a) are the heating effects due to jets generated in Seyfert galaxies able to affect kpc-scales and (b) can they be a source of feedback? Observations provides appealing clues to answering these questions.

Several observational studies of Seyferts suggest that Seyfert jets may indeed reach kpc-scales and beyond. The Very Large Array surveys of Colbert et al. (1996) and Gallimore et al. (2006) find that extended kpc-scale radio outflows most probably associated with jets are common in Seyfert galaxies, although they cannot rule out entirely the possibility that these outflows arise from starburst-driven winds in some sources. These radio outflows have random orientations with respect to the host galaxy, a result which is maintained for jets observed from sub-kpc to kpc scales (Kinney et al., 2000; Schmitt et al., 2001), and have distorted, uncollimated morphologies. Therefore the observed radio outflows could provide a broad effective surface for depositing the jet power in the surrounding gas, as envisioned in our simulations. Several targeted observations of individual Seyfert galaxies also demonstrate that kpc-scale radio jets are not a rare feature in disk galaxies (e.g., Wilson & Ulvestad 1983; Morganti, Oosterloo, & Tsvetanov 1998; Cecil et al. 2000; Whittle & Wilson 2004; Keel et al. 2006; Kharb et al. 2006; Middelberg et al. 2007).

X-ray observations also reveal kpc-scale outflows of X-ray-emitting hot plasma in Seyfert galaxies (Colbert et al., 1998; Rupke, Veilleux, & Sanders, 2005), cospatial with the large scale radio outflows. Two favoured interpretations are that these kpc-scale X-ray outflows are generated either by the AGN jets that entrain and heat gas as they flow out of the nucleus, or through heating of the surrounding gas via thermalisation of the kinetic power of the jet, generating a wind. Either way, these processes bear resemblance to the generation of cavities or bubbles of X-ray emitting plasma in the central elliptical galaxies of clusters of galaxies, which are inflated by the radio jets (e.g., Dalla Vecchia et al. 2004; Bîrzan et al. 2004; Fabian et al. 2006). These observations therefore provide evidence for the possible importance of feedback due to the radio jets in Seyfert galaxies, which could heat the circumnuclear plasma and generate X-ray emitting outflows.

Observational estimates of the kinetic power carried by the jets in Seyferts using different methods indicate that it is in the range $\sim 10^{42}-10^{43}{\rm erg}\ {\rm s}^{-1}$ (e.g., Wilson & Ulvestad 1983; Morganti, Oosterloo, & Tsvetanov 1998; Bicknell et al. 1998; Capetti et al. 1999; Kharb et al. 2006). This observational range is in agreement with the values of the jet power which are generated during the radio-loud episodes for $z\lesssim 1$ in the simulation, $L_{\rm jet}\approx 10^{42}-10^{43}\;{\rm erg}\ {\rm s}^{-1}$ (Fig. 13). It is worth noting that the jet power in our simulation is estimated using a physical model for the energetics of the jet, which depends on fundamental parameters of central engine, especially the accretion rate which is self-consistently calculated.

The fraction of the jet power generated by the AGN which is actually available to heat the gas on kpc-scales is a matter of debate. In the ‘reference’ simulation we adopt that a fraction $\eta_{\rm AGN}=0.1$ of the jet power is deposited as thermal energy of the diffuse gas, although it is possible that this may be an overestimate of the value of $\eta_{\rm AGN}$. For instance, based on ram pressure and terminal velocity considerations of the lobes of kpc-scale radio outflows in Seyferts, Gallimore et al. (2006) argues that the Seyfert jets lose almost all their power within the inner kpc, implying that $\eta_{\rm AGN}$ may be considerably smaller than 0.1. The other simulations we present encompass the uncertainty in the value of $\eta_{\rm AGN}$. In this respect, the simulation adopting $\eta_{\rm AGN}=0.01$ does not show any significant difference from the simulation without AGN feedback ($\eta_{\rm AGN}=0$). On the other hand, when 100% of the jet energy is transferred to the surrounding star forming gas (‘ism-fb’ model), most of the feedback energy is rapidly radiated away, and the results are very similar to the ‘reference’ model in which 10% of jet energy is given to the ambient halo gas. Thus, a better knowledge on the “microscopic” physics of jet-ISM interactions is required in order to improve the understanding of the effects of AGN feedback in Seyfert galaxies.

Our simulations show that during the cosmological formation of a typical spiral galaxy, a low-luminosity AGN (Seyfert or LINER) is formed which has a stochastic cycle of activity, with many short-lived recurrent accretion events accompanied by jet outbursts. This finding agrees with the conclusion of Sanders (1984) that the Seyfert activity is the result of short-lived stochastic accretion events. He arrived at this conclusion by considering plausible mechanisms that may fuel the central black holes of Seyferts, which are stochastic in nature, and our simulations track the evolution of the AGN by implementing one of these mechanisms: the accretion of molecular clouds to the galactic centre. Our findings therefore provide support to the assumption of Hopkins & Hernquist (2006), central to their work, that the central supermassive BHs of Seyferts stochastically increase their mass by accretion of cold gas.

In the period $z\lesssim 1$ after the galactic disk is formed, the AGN is active only during a fraction $\approx 9\%$ of the its life in the ‘reference’ model, $\approx 25\%$ in the ‘weak-fb’ model and $\approx 1\%$ in the ‘ism-fb’ model. This implies that the lifetime of Seyfert activity in the spiral galaxy is $\approx 10^{9}$ yr for the ‘reference’ or ‘weak-fb’ model and $\approx 10^{8}$ yr for the ‘ism-fb’ model. This is in agreement with statistical Seyfert lifetime $\sim 10^{8}-10^{9}$ yr estimated by counting the fraction of spiral galaxies that are observed to be Seyferts (Woltjer, 1968; Sanders, 1984; Ho et al., 1997). Sanders (1984) estimated that each single episode of activity would not last longer than $10^{6}$ yr based on the size of Seyfert jets, a prediction which is in agreement with the detailed study of Kharb et al. (2006). Unfortunately the time-step we used to output snapshot files is 16.8 Myr and therefore we may miss many of such short-lived events. Considering this, the ‘ism-fb’ model may be more comfortable than the ‘reference’ model because there is no long-lived accretion event in the ‘ism-fb’ model for $z\lesssim 1$ (lower-left panel of Fig. 15). It is also possible that the observed short duration of activity is not caused by the large-scale mass accretion to the galactic centres but by the physical processes in the accretion disk (e.g., instabilities) which we do not take into account in this paper.

7 Summary and Conclusions

In this paper, we have introduced a new methodology that enables us to simultaneously model star formation, SN feedback, BH accretion, and AGN feedback in cosmological simulations of galaxy formation. Our treatment of black hole accretion is based on a theoretical model of the mass accretion into the galactic centre due to radiation drag. This connects the kpc-scale star formation activity in the galactic centre with the BH accretion rate. Furthermore, motivated by the recent semi-analytic models of Croton et al. (2006) and Bower et al. (2006), we have distinguished two distinct accretion modes of BH fuelling, namely standard (optically thick, geometrically thin) Shakura-Sunyaev disks and (geometrically thick, optically thin) RIAFs. This allows us to distinguish radio loud and radio quiet AGNs depending on their mass accretion rate. In this model, we only consider AGN radio feedback in which we assume a fraction of the jet power is thermally coupled with the surrounding diffuse gas, and assume RIAFs are responsible for production of powerful jets. As the first application of this model, we have performed cosmological simulations of the formation of a disk galaxy and its coevolution with the central AGN.

Our results reveal a fundamental AGN-starburst connection during the evolution of the galaxy and the central AGN. During the starburst period ($z\gtrsim 1$) there is a lot of gas available to be accreted and the AGN is fuelled at high accretion rates, as a consequence the accretion disk is in the standard/thin state and the AGN is relatively luminous and produces weak jets (radio-quiet state). The starburst phase ceases at $z\sim 1$ because most of the gas at the galactic centre is consumed by star formation or blown out by the associated SN feedback. Afterwards, the AGN fuelling rate drops and the accretion disk then switches to a radiatively inefficient (RIAF) state. Due to a combination of small accretion rates and the low radiative efficiency of the RIAF, the AGN becomes a low-luminosity AGN or Seyfert nucleus which is almost invisible compared to the host galaxy. At this point even a little star formation activity around the central BH triggers intense jet production that regulates further gas cooling to the centre. Thus, in the quiescent star forming regime after $z\sim 1$ the AGN has a stochastic fuelling, undergoing several short-lived recurrent episodes of activity. During these accretion episodes the AGN brightness increases and jet outbursts are generated, which correspond to the radio-loud state. In-between these Seyfert episodes the nucleus is turned off and the galaxy is inactive.

Below $z\sim 1$ AGN feedback becomes important if 10% of jet energy is available to heat the ambient diffuse halo gas (or 100% is fed into the surrounding ISM). This powerful episodic AGN feedback suppresses star formation in the disc and almost completely halts the star formation in the bulge. This effect is minor when the wind speed due to SN feedback is fast compared to the circular velocity of the host halo because the SN feedback alone can suppress the star formation in the bulge reasonably well. When the wind speed is slow compared to the circular velocity of the host halo, the AGN feedback becomes more important. This leads to the formation of an old, red elliptical galaxy in the ‘slow-wind’ model.

Many properties of the simulated low-luminosity AGN and its jets are in broad agreement with observations of Seyfert galaxies and their jets: namely its typical nuclear bolometric luminosity and radio luminosity, the coincidence of the radio-loud state with the Seyfert episodes, the fact that the jet may reach kpc-scales and strongly affect the surrounding gas on these scales, and the predicted AGN lifetime and duration of individual episodes of activity. The stochastic pattern of nuclear activity that emerges from our simulation seems to confirm the prediction of Sanders (1984) for Seyfert galaxies.

The ratios between the central BH mass and the mass of the host spheroid at $z=0$ in our simulations are $\sim 10^{-3}$ regardless of the strength of either SN or AGN feedback. This result is direct outcome of the radiation drag model which relates the star formation rate in the central star-forming region and the mass accretion rate to the central BH. The ratio is in good agreement with the observed ratio (Kormendy & Richstone, 1995; Magorrian et al., 1998; Merritt & Ferrarese, 2001; McLure & Dunlop, 2002), although the efficiency of the radiation drag we assumed to obtain this ratio is somewhat higher than the value suggested by Kawakatu & Umemura (2002). However, as we have discussed, the geometry and density structure of star forming regions can easily change the efficiency. We also note that recent semi-analytic models (Baugh et al., 2005; Nagashima et al., 2005a, b) and cosmological simulations (OEFJ05, ) suggest that stars may be born with top-heavy IMFs (and thus stronger radiation fields and higher radiation drag efficiencies) in starbursts. The accretion rate by radiation drag can be much higher with this assumption because stellar populations with top-heavy IMFs produce more radiation.

The self-consistent treatment of BH accretion and associated AGN feedback in galaxy formation enables us to explore the interplay between galaxies and AGNs in cosmological simulations of galaxy formation. The self-regulation of mass accretion due to feedback from RIAFs produces galaxies whose AGN properties are well matched to the observations. The model is also qualitatively consistent with recent semi-analytic models that successfully reproduce many of the global properties of the galaxy population. The effects of our AGN feedback implementation on formation of groups of galaxies, where more powerful AGN feedback from larger BHs is expected, will be presented in a forthcoming paper.

Acknowledgments

We greatly appreciate the detailed reading and shrewd comments of the anonymous referee that strengthened this paper. We thank Takayuki Saitoh, Nozomu Kawakatu, Yasuyuki Watabe, Carlos Frenk, Tom Theuns and Thaisa Storchi-Bergmann for useful discussions. TO acknowledges financial support from the Japan Society for the Promotion of Science for Young Scientists (1089) under which part of this work was carried out. RSN acknowledges support from the Brazilian institutions CNPq and CAPES, and by the European Commissions ALFA-II programme through its funding of the Latin-american European Network for Astrophysics and Cosmology (LENAC). RGB thanks PPARC for the support of a senior research fellowship. We are grateful to Volker Springel for making the GADGET2 code public. All simulations were performed on the Cosmology Machine at the Institute for Computational Cosmology in Durham University. Part of analyses were carried out on GRAPE system at the Center for Computational Astrophysics, CfCA, of the National Astronomical Observatory of Japan (Project-ID: g06b13).