Mon. Not. R. Astron. Soc. 000, 1–22 (2007) Printed 24 October 2018 (MN LATEX style file v2.2)
The impact of radio feedback from active galactic nuclei in
cosmological simulations: Formation of disk galaxies
Takashi Okamoto1⋆, Rodrigo S. Nemmen2, Richard G. Bower1
1Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham, DH1 3LE
2Instituto de Fı́sica, Universidade Federal do Rio Grande do Sul, Campus do Vale, Porto Alegre, RS, Brazil
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 en-
ergy. 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 pro-
duction in accretion disks, we assume that only the RIAF is responsible for the AGN feedback.
The focus of this paper is to investigate the interplay between galaxies and their central black
holes during the formation of a disc galaxy. 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. This suppression of the star formation regulates mass accretion onto
the black hole and associated AGN feedback. As a result, the nucleus becomes a stochas-
tically fuelled low-luminosity AGN (Seyfert galaxy) with recurrent short-lived episodes of
activity after the star bursts. During the “on” events the AGN produces reasonably power-
ful jets (radio-loud state) and is less luminous than the host galaxy, while in the “off” phase
the nucleus is inactive and “radio-quiet”. 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 ratios between the central black hole
mass and the mass of the host spheroid at z = 0 are ∼ 10−3 regardless of the strength of
either supernova feedback or AGN feedback because the radiation drag model directly relates
the star formation activity in the galactic centre and the mass accretion rate onto the central
black hole.
Key words: black hole physics – galaxies: active – galaxies: starburst – galaxies: Seyfert –
galaxies: evolution – galaxies: formation
1 INTRODUCTION
Recent observations have suggested a fundamental connec-
tion between active galactic nuclei (AGN) and the formation
of galaxies. Firstly, there is now a well established correla-
tion between the properties of galaxies and the masses of
the black holes (BH) at their centres. The mass of the cen-
tral black holes tightly correlates with the mass of the host
galactic bulges (with a median BH-to-bulge mass ratio of
≃ 0.001, e.g. Kormendy & Richstone 1995; Magorrian et al.
1998; Merritt & Ferrarese 2001; McLure & Dunlop 2002;
⋆ E-mail: takashi.okamoto@durham.ac.uk
Marconi & Hunt 2003), as well as with the stellar velocity disper-
sion 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
c© 2007 RAS
http://arxiv.org/abs/0704.1218v3
2 T. Okamoto, R. S. Nemmen and R. G. Bower
sufficient energy to offset the cooling radiation from the cluster
and potentially explain why so little cool (T < 3 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 im-
provements 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 dynami-
cal 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 ex-
ists 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 sur-
rounding 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 cosmolog-
ical 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 proper-
ties 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 sim-
ulations 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 distin-
guish two fundamentally different modes of AGN accretion: radia-
tively efficient geometrically thin accretion flows (standard disks;
Shakura & Sunyaev 1973) and geometrically thick, radiatively in-
efficient 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 ra-
diation drag model directly connects the star formation rate around
BHs and mass accretion rate onto the BHs. Given the limited nu-
merical 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 cen-
tral supermassive BHs in their bulges. In this paper we hence fo-
cus 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 tar-
geted 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 ∼ 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 sam-
ples of Seyfert and LINER galaxies indicate that radio outflows
with extent & 1 kpc are a common feature in Seyfert galax-
ies, and presumably are driven by the AGN jets (Colbert et al.
1996; Gallimore et al. 2006). Gallimore et al. (2006) in particu-
lar find that & 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 ki-
netic power carried by the Seyfert jets as estimated from obser-
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 3
vations can be as high as 1043erg 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 galax-
ies and the evolution of the nuclear activity during the formation of
the disk galaxy are open problems, which we attack using our cos-
mological 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 micro-
scopic 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 re-
sults 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 for-
mation, we construct a physically motivated, well-controlled
model of the star-forming gas. Based on the simplest star for-
mation and feedback recipes used in the first generation of hy-
drodynamic 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 inter-
actions between the phases (Yepes et al. 1997; Springel & Hern-
quist 2003 hereafter SH03; Samland & Gerhard 2003; OEFJ05;
Scannapieco et al. 2006 but see Booth et al. 2007 where they ex-
plicitly model ‘clouds’ by means of sticky particles in SPH sim-
ulations). Our model is based on the one used in OEFJ05 but has
been significantly modified. We describe the details in the follow-
ing subsections.
2.1 Gas cooling
We compute the non-equilibrium cooling/heating rate and ionisa-
tion 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 ≃ 6. Inverse
Compton cooling is also included.
Metals are carried by particles and once assigned to a parti-
cle, 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 par-
ticle that consists of two distinct phases, i.e. hot ambient gas and
cold clouds, once its gas density, ρ, exceeds a density threshold,
We assume that the cold clouds form and grow through ther-
mal instability, that is
Λnet(ρh, uh, Z)
uh − uc
, (1)
where Mc and Mh are mass in cold clouds and mass in the hot
phase associated with the particle, respectively, Λnet is the cooling
function for gas of metallicity Z, and uh and uc represent specific
internal energy of the hot phase and cold clouds, respectively. This
implies that we assume that the gas is thermally unstable when ρ >
ρth. Cold clouds remain at a fixed temperature, Tc = 100K hence
uc is a constant as well. Note, however, as far as Tc is much lower
than 104 K our results are hardly affected by the choice of Tc.
We allow cold clouds to have a mass spectrum
φ(m) ≡
. (2)
Observationally α is in a range between α = 1.5 to 1.9
(Solomon & Rivolo 1989; Fukui et al. 2001; Heyer et al. 2001).
Since the effect of changing the value of α is completely degen-
erate with changing other parameters we will introduce such as star
formation efficiency and thermal conduction efficiency, we simply
fix α = 1.7 and the mass range of clouds, 102–106M⊙, throughout
this paper.
Following Samland & Gerhard (2003), we assume spherical
clouds which follow the mass-radius relation of Elmegreen (1989)
= 190
4 , (3)
where P4 =
104K cm3
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 104M⊙
to 106M⊙ 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
ṁ∗ =
, (4)
where m∗ is stellar mass, tsf = tdyn/c∗ is the star formation
timescale for a cloud of mass m, c∗ is a dimensionless star for-
mation efficiency parameter that must be less than unity, and tdyn
is the dynamical timescale of the cloud. For the mass-radius rela-
tion in equation (3), the dynamical timescale of a cloud of mass m
is given as
tdyn =
32Gρ(m)
≃ 0.32 P
Myr, (5)
where ρ(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 1M⊙ with the mass spectrum φ(m) is ob-
tained by integrating ṁ∗ over the range of masses that are eligible
to form stars,
Ṁ∗ =
Z 106M⊙
104M⊙
ṁ∗(m)φ(m)dm, (6)
c© 2007 RAS, MNRAS 000, 1–22
4 T. Okamoto, R. S. Nemmen and R. G. Bower
where the normalisation
Z 106M⊙
102M⊙
mφ(m)dm = 1
is applied.
2.4 Cloud evaporation
The different phases of the ISM exchange mass by evapora-
tion. As in Samland & Gerhard (2003), we use the model of
Cowie & McKee (1977) for the cloud evaporation, that is,
ṁEVP =
= ηEVP
16πµκr
≃ 4.4× 10
M⊙ Myr
, (7)
where ηEVP, µ, κ, r, and kB are the dimensionless conduction ef-
ficiency 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 1M⊙ then becomes
ṀEVP =
Z 106M⊙
102M⊙
ṁEVP(m)φ(m)dm. (8)
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)
φ∗(m). Throughout this paper we assume that the IMF is always
the Salpeter IMF (Salpeter 1955) whose lower and upper mass lim-
its are 0.1 and 100 M⊙, 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 ηSN10
51 erg of
energy to surrounding gas, where η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 succes-
sor 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 de-
scribe how the physical processes described above are implemented
in our code.
For ρ > ρ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 ∆t as
1− exp
, (9)
where the average star formation timescale t∗ is obtained by using
equation (6):
t∗(P ) =
Ṁ∗(P )
Myr. (10)
We use m∗ = morig/3 as in OEFJ05, where morig 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⊙ explode
as SNe II. The probability pII of a star particle having an event of
SN II explosion during a time-step ∆t is given by
pII =
R t+∆t
rII(t
′)dt′
R t(8M⊙)
rII(t′)dt′
, (11)
where t is the age of the SSP, rII(t) is the SN II rate for the SSP
of age t and t(8M⊙) is the lifetime of the star of original mass
8M⊙. 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⊙) 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 ∆tIa, 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 ∆t as
pIa =
R t+∆t
rIa(t
′)dt′
R t0+∆tIa
rIa(t′)dt′
, (12)
where t0 is the age at which the previous SN Ia time-step finished
and rIa(t) is the SN Ia rate for a SSP of age t. We adopt ∆tIa =
100 Myr in this paper. The probabilistic method described above
statistically relaxes the IRA.
Finally, we update the mass of the hot phase Mh and the spe-
cific internal energy of the hot phase uh. The new mass of the hot
phase Mh
′ is given by solving the thermal energy equation implic-
itly using equation (1),
= Mh +∆MEVP −
Λnet(ρh
′, uh, Z)
uh − uc
∆t, (13)
where ∆MEVP is the mass which evaporates from cold clouds dur-
ing the time-step, and the new density of the hot phase ρh
′ is given
by the mass conservation, MSPH = Mh+Mc, and by the assump-
tion that the volume occupied by a SPH particle remains the same
during the time-step, i.e.
= Vol c +
where Volc 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 evap-
oration also change the specific energy of the hot phase. The new
specific energy of the hot phase uh
′ due to the non-hydrodynamic
processes is also calculated implicitly as
= uh +
∆Q+ (uc − uh
′)∆MEVP
, (14)
where ∆Q is the thermal energy received by the SPH particle dur-
ing 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
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 5
which they showed that the ISM is supported by the turbulent mo-
tion originated in the self-gravity of the gas and galactic rotation.
They found that energy injection by SNe hardly changes the struc-
ture 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 cos-
mological simulation performed with an adaptive mesh refinement
code.
Lacking the heating by the ∆Q term, the ISM loses its pres-
sure quickly and would be fragmented. We therefore introduce min-
imal heating by imposing the minimum pressure for the star form-
ing gas as a function of density,
Pmin(ρ) ∝ ρ
γeff for ρ > ρth, (15)
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
γ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 nor-
malise the minimum pressure as Pmin(ρth) = (γ−1)u4ρth, where
u4 denotes the specific internal energy of the gas with temperature
104 K. Once pressure from the hot phase drops below Pmin, we set
the temperature of the hot phase to Thot and recalculate the mass
in cold clouds so that the ambient pressure becomes Pmin. We set
Thot to 10
6 K in this paper but our results are hardly affected by the
value of Thot as far as Thot >> 10
4 K. In quiescent disks, most of
the star forming gas has Pmin by construction. However, for exam-
ple during a violent merger, the pressure can be significantly higher
than Pmin 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 re-
gions. 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, vw , when it leaves
the disk. Then, we assign the probability pw with which the gas
particle is added to the wind as
MSPHvw
, (16)
where MSPH is the mass of the gas particle and ∆Q is the feedback
energy received during the time-step. By doing this, our model be-
comes insensitive to the numerical resolution, time-step in the sim-
ulations, and the number of neighbours we use to distribute the
feedback energy. When a particle is added to the wind, we modify
the velocity v of the particle according to
= v + vwn. (17)
For the unit vector n, we choose random orientation along the di-
rection v × agrav , where agrav is the gravitational acceleration.
Therefore wind particles are preferentially ejected along the rota-
tion 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 interac-
tions are enabled again once the density of the particle has fallen
to 0.1 ρth, or once a time of ∆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, ρth = 5 ×
10−25 g cm−3, the dimensionless star formation efficiency, c∗ =
0.005, the dimensionless conductivity efficiency, ηEVP = 0.1, the
SN feedback efficiency, ηSN = 1, and the wind velocity vw = 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 ρ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 Ṁw/Ṁ∗ ≃ 3. There will be further
contribution from SNe Ia at later times.
In this subsection, we show the behaviour of our model in ide-
alised simulations of disk formation from virialised gas in the ex-
ternal NFW potential (Navarro et al. 1996). Here, Mvir 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, λ = J |E|1/2/(GM
vir ), 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
j(r) ∝
In Fig. 1 we show the star formation rates in the simulations
for Mvir = 10
12M⊙ with three different resolutions: Ngas =
10000, 50000, and 250000. The star formation rates quickly con-
verge. The simulation with Ngas = 50000, which is a compara-
ble resolution that we will use in cosmological simulations, has
very similar star formation history to that in the simulation with
Ngas = 250000. Hereafter we thus only show the simulations with
Ngas = 50000.
Now we show how star formation rates depend on the values
of feedback parameters, ηSN and vw . We here run the simulations
with two additional models: the slow wind model (ηSN = 1 and
vw = 250 km s
−1) and the weak feedback model (ηSN = 0.5 and
vw = 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 ηw ≃ 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 Mvir = 10
12M⊙ and 10
10M⊙.
Since vw = 500 km s
−1 exceeds the escape velocity of the halo
withMvir = 10
12 M⊙ and vw = 250 km s
−1 does not, the star for-
mation 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, vw = 250 km s
−1, will return to
the disk in ∼ 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.
c© 2007 RAS, MNRAS 000, 1–22
6 T. Okamoto, R. S. Nemmen and R. G. Bower
Figure 1. Star formation rates in simulations of disks growing within ha-
los of Mvir = 10
12M⊙, with different resolutions. The simulations with
Ngas = 10000, 50000, and 250000 are indicated by the solid, dotted, and
dashed lines, respectively.
However situation changes for the halo of Mvir = 10
10 M⊙. 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 ∼ 10−3 M⊙ 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 cos-
mological 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 den-
sity of stars younger than 3 × 107 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 Mvir = 10
10M⊙ 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 dif-
ferent epochs. We find that the global relation in the 1010M⊙ 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 0.
When we find a halo that contains more than Nth dark matter parti-
cles and no BH particles, we put a seed BH on a density peak of the
Figure 2. Star formation rates in simulations of disks growing within halos
of Mvir = 10
12M⊙ (upper panel) and Mvir = 10
10M⊙ (lower panel).
The solid, dotted, and dashed lines indicate the fiducial, slow wind, and
weak feedback models, respectively.
stellar distribution in the halo. We assume 260 M⊙ as the mass of
the seed black hole. This mass is consistent with the theory of stel-
lar evolution, which shows that the nuclear burning in very massive
stars above 260 M⊙ is unable to halt gravitational collapse (e.g.
Heger et al. 2003). Note that as far as the seed BH mass is suffi-
ciently small, our results hardly change because we do not impose
an upper limit on the mass accretion rate. Throughout this paper
we use Nth = 128. The simulation employing Nth = 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 τ̄ = χdρcrc ≃ χdmc/r
where ρc, mc, and rc are the density, mass, and a size of a cloud.
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 7
Figure 3. Star formation rate per unit area versus gas surface density in an
idealised simulation of disk formation in a halo with Mvir = 10
12M⊙ (top
and middle panels). The surface star formation rate densities are computed
using the surface density of stars younger than 3 × 107 yr in cylindrical
bins. The symbols represent star formation rates at t = 0.5, 1, 2, 3 Gyr,
indicated by plus signs, diamonds, triangles, and squares respectively. The
top and middle panels respectively show the fiducial and weak feedback
models. The galaxy in the halo with Mvir = 10
10M⊙ is shown in the
bottom panel, and there both the surface star formation rate density and
the surface gas density are averaged over the galaxy. Points represent the
evolution of the relation from 0.25 Gyr to 3 Gyr with a time-step of 0.25
Gyr. The solid and dotted lines indicate the target relation and the cut-off
surface density respectively.
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 χd is given by
χd = 300 cm
2g−1(ad/0.1µm
−1)(ρs/g cm
−3)(Z/0.3Z⊙),
where ad is the grain radius, ρs is the density of solid material
within the grain (e.g. Spitzer 1978), and Z is the metallicity of
gas. We use ad = 0.1 µm, ρ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
τSFR(t) ≃
Mc(t)
RSFR(t)2
, (18)
where RSFR and Mc 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
Ṁdrag = ηdrag
LSFR(t)
(1− e
−τSFR(t)), (19)
where LSFR is the total luminosity of stars in the star form-
ing region. We calculate LSFR by summing up the bolometric
luminosities of all stars, each of which is obtained as a func-
tion of its age from a look-up table generated by PÉGASE2
(Fioc & Rocca-Volmerange 1997). We here ignore the metallic-
ity 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 ηdrag is maxi-
mally 0.34 in the optically thick regime assuming a spherical star
forming region with a uniform density. In principle ηdrag can be
much larger or smaller depending on the geometry and density
structure of the star forming region (Kawakatu private communica-
tion). We thus change the efficiency and adopt ηdrag = 1 in this pa-
per, 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, RSFR, we obtain χd, Mc, and LSFR in equations (18) and (19)
directly in simulations, and hence we can compute τSFR and Ṁdrag
self-consistently; the effects of outflows are also implicitly incorpo-
rated in the estimate of these variables because we calculate them
self-consistently at each time-step of the simulation. In order to de-
fine the size of the star forming region, RSFR we start from a sphere
around a BH that contains 40 SPH particles, and then we search the
radius of the sphere that maximises Ṁ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, ρth, we consider that
there is not enough gas around the BH to fuel accretion and we
simply set Ṁ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 forma-
tion in the disc and bulge correlate with the mass accretion rate,
Ṁdrag, in our simulations later.
Strictly speaking, one should distinguish Ṁdrag from the ac-
cretion rate onto the BH, Ṁ , because the radiation drag is not likely
to remove the angular momentum thoroughly; some residual an-
gular momentum will terminate the radial contraction of the ac-
creted 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
Ṁ = Ṁdrag for simplicity. The time delay in mass fuelling intro-
duced 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 di-
rectly heat diffuse hot gas confined in dark halos through the pro-
duction 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 ac-
cretion 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
c© 2007 RAS, MNRAS 000, 1–22
8 T. Okamoto, R. S. Nemmen and R. G. Bower
basic features are supported by recent relativistic numerical simu-
lations 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 cos-
mological simulations of the spin evolution of BHs through merg-
ers 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 ≃ 0.7 − 1. We thus only consider spinning BHs in this pa-
per. Meier (2001) parameterised the jet luminosity for standard thin
discs as
jet ≈ 10
erg s
”−0.1
× (1 + 1.1j + 0.29j
), (20)
and for RIAFs as
jet ≈ 10
erg s
αRIAF
× (0.55f
+ 1.5fj + j
), (21)
where αSD and αRIAF are the viscosity parameters for standard
thin discs and for RIAFs, respectively, m9 is the black hole mass
in units of 109M⊙, ṁ ≡ Ṁ/ṀEdd is the accretion rate scaled
to the Eddington rate assuming the 10% energy conversion effi-
ciency (ṀEdd ≡ 22 m9 M⊙ 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 ∆EjetFB = ηAGNLjet∆t is delivered as thermal en-
ergy to the neighbouring 40 diffuse gas particles that have density,
ρ < 0.1ρth. In our simulations, the size of the region heated by
the AGN feedback is typically ∼ 10 kpc. Unless otherwise stated,
we use αSD = αRIAF = 0.1, j = 0.5, and η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 ṁ, and the critical value of the accre-
tion rate which sets the division between the two regimes is given
by ṁcrit ≈ α
2 (e.g. Narayan et al. 1998). Since the RIAF solu-
tion ceases to be valid above ṁcrit, for ṁ 6 ṁcrit a RIAF exists
and for ṁ > ṁ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 ṁ = ṁcrit. For the parameters we chose, we have the following
jet production efficiency in RIAFs,
jet ≈ 2.6× 10
for ṁ 6 ṁcrit, (22)
which is much higher than that in standard discs, given by
jet ≈ 8.1× 10
for ṁ > ṁcrit, (23)
where in equation (23) we adopt m9 = 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 en-
ergy” 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 ΛCDM with the mean matter density Ω0 = 0.3, Hub-
ble parameter, h ≡ H0/100 km s
−1Mpc−1 = 0.7, cosmological
constant term, ΩΛ ≡ Λ0/(3H
0 ) = 0.7, amplitude of mass fluc-
tuations, σ8 = 0.9, and mean baryon density, Ωb = 0.04. The
size of the periodic simulation box is 35.325h−1Mpc and the ini-
tial redshift is z = 49. We put the high-resolution dark matter
particles and gas particles in the region where a halo with mass
Mvir = 1.2 × 10
12h−1M⊙ forms at z = 0. The region external
to this is populated with high-mass dark matter particles (bound-
ary particles), the function of which is to reproduce the appro-
priate tidal fields. The circular velocity, spin parameter and col-
lapse redshift of the selected halo are vc(rvir) = 155 km s
λ ≡ J |E|1/2/(GM5/2) = 0.038 and zc ≃ 1.5, respectively,
where zc 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 parti-
cles are ∼ 2.6× 106 and ∼ 1.7× 107h−1M⊙, respectively. Thus
the Nth = 128 implies that halos larger than ∼ 2.6× 10
9h−1M⊙
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
ǫ = 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ǫ.
4.1 Black hole particles
When we find a halo with Ndm > Nth 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 ac-
cording 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:
MBH −Mp
, (24)
where MBH and Mp are, respectively, the BH mass and the mass
of the BH particle, NSFR is the number of gas particles in the star-
forming region, and Mn 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 impor-
tant 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
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 9
BHs are within a gravitational softening radius and they are grav-
itationally bound. Therefore the BH merger rates obtained in our
simulations are likely to be the maximum possible rate. Techni-
cally, 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 galac-
tic 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 evolu-
tion. As a result, BH particles tend to oscillate around galactic cen-
tres with the velocity dispersions of the stellar systems. To circum-
vent this problem, we allow BH particles to behave as peak tracers
at galactic centres. At each time-step, we compute local stellar den-
sity fields within the gravitational softening radii of BH particles,
ǫ. We then displace the BH particles by a small distance, ∆l, along
the local density gradients. After several tests, we chose
∆l = MIN(0.01ǫ, 0.03 |v| ∆t), (25)
where v is the velocity of a BH particle and ∆t is the time-step.
This ∆l correction sufficiently suppresses the oscillation of the BH
particles at galactic centres.
We adopt a maximum time step for BHs, ∆tmax = 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 ∆tBH = min(∆tmax,∆tdyn), where ∆tdyn 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 par-
ticles defined by their densities, ρ < 0.1ρth in the form of thermal
energy using SPH smoothing. Since only 1% of jet energy is avail-
able 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 red-
shift (z > 1). The disk formation starts at z ∼ 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 parti-
cles 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
Figure 5. Projected stellar distributions within 50 h−1 kpc boxes centred
on the galaxies at z = 0. The left, middle, and right panels show the ‘no-fb’,
‘weak-fb’, and ‘reference’ simulations, respectively. The viewing angles are
chosen to be edge-on for upper panels and face-on for lower panels. The
brightness indicates the projected stellar mass density, and the same scaling
is used for each simulation.
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 in-
teresting 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 gradi-
ents caused by the thermal feedback to blow the gas away. Unfortu-
nately, 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 sur-
face 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
c© 2007 RAS, MNRAS 000, 1–22
10 T. Okamoto, R. S. Nemmen and R. G. Bower
Figure 4. Snapshots of the ‘reference’ model (x–y projections). The rows correspond to z = 6, 3, 2, 1, and 0 from top to bottom. The first and the second
columns show gas and stellar distributions in the comoving 3 h−1 Mpc cubes centred on the main progenitor, respectively, where the brightness represents
three dimensional density. The third and the fourth columns show gas and stellar distributions in the physical 50 h−1 kpc cubes, respectively. The circles in
each panel indicate the position of the BH particles and their radii are proportional to the log of the BH masses.
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 11
Figure 6. Surface brightness profiles in the I band for the ‘no-fb’ (pluses),
‘weak-fb’ (asterisks), ‘reference’ (diamonds) galaxies, and for the galaxy
with the ‘shock-burst’ model of Okamoto et al. (2005; dotted-line). The
solid lines are exponential fits for r < 15 kpc with the scale lengths 5.4,
5.6, and 4.0 kpc for the ‘no-fb’, ‘weak-fb’, and ‘reference’ galaxies, respec-
tively.
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 com-
ponents.
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 for-
mation rates at z ∼ 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 re-
sponsible 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 be-
cause 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 for-
mation 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 ∼ 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 ∼ 10% of jet energy has to be given to the halo
Figure 7. The formation history (star formation rate as a function of time)
of the stars that lie within 25h−1 kpc from the galactic centre at z = 0. The
solid and dotted lines correspond to the galaxies in the simulation without
and with BHs, respectively.
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 re-
lated 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 Ṁ is different from the
mass growth rate, ṀBH; the latter includes the mass increase due to
BH mergers. In all simulations, the central BHs rapidly evolve un-
til t ∼ 5.5 Gyr at which the starburst ceases. The final BH masses
reach ∼ 107M⊙. 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 ex-
plains the decline of the accretion rate onto the central BH. At
t ≃ 5.5 Gyr ṁ reaches the critical accretion rate ṁcrit, below
which the accretion disk becomes a RIAF. Our model has the max-
imum 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 ∼ 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 sup-
c© 2007 RAS, MNRAS 000, 1–22
12 T. Okamoto, R. S. Nemmen and R. G. Bower
Figure 8. Upper panel: the mass growth of the central BH of the main progenitor. Lower panel: the mass accretion rate onto the central BH of the main
progenitor in units of the Eddington rate (ṁ = Ṁ/ṀEdd). The horizontal dashed line indicate the critical mass accretion rate ṁcrit below which the
accretion flow becomes radiatively inefficient (RIAF). The ‘no-fb’, ‘weak-fb’, and ‘reference’ simulations are shown from left to right panels.
press the star formation because it injects energy not into the star
forming dense gas but into the diffuse halo gas (ρ < 0.1ρ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 pro-
duction of jets in the RIAF. The fact that the central surface bright-
ness is almost the same in the three models in spite that the galaxy
with AGN feedback hardly form stars at its centre after z ∼ 1 sug-
gests that the vast majority of stars at the centre have already been
formed at z ∼ 1 in all galaxies.
We now return to the separation of the disk and bulge compo-
nents. As noted earlier, it is best to separate the bulge and disk com-
ponents using a dynamical decomposition (Abadi et al. 2003). For
this, we first compute the angular momentum Jz of each star par-
ticle parallel to the net angular momentum of stars within 10h−1
kpc, and the angular momentum of the corotating circular orbit,
Jc(E), where E is the energy of each particle. The ratio Jz/Jc
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
Jz/Jc(E) ≃ 1 and we see that such a component clearly dom-
inates 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
zlmm > 2) using a different feedback recipe from ours. To iden-
tify 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 Jz/Jc(E) 6 0 are identified
as a half of the spheroidal component. Another half of the spheroid
with Jz/Jc(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
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 13
Figure 9. Mass-weighted probability distributions of the orbital circular-
ity, Jz/Jc(E), of stars within 25h−1 kpc from the galactic centres. The
solid, dotted, and dashed lines indicate the ‘no-fb’, ‘weak-fb’, and ‘refer-
ence’ galaxies, respectively.
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 feed-
back.
In Table 1 we show the total and spheroid’s masses, total lu-
minosity, 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 PÉ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 histo-
ries 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 af-
ter z ∼ 1, it is safe to compare the star formation rates and the mass
accretion rates at low redshift. We find that there are good correla-
tions between the star formation rates of the spheroid and the mass
accretion rates onto the central BHs at low redshift, in particular
between t ∼ 5.5 and 8 Gyr. After the starburst, the star forma-
tion 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 for-
Figure 10. Star formation histories for the disk and spheroid stars within
25h−1 kpc from the centre of each galaxy at z = 0. The ‘no-fb’, ‘weak-fb’,
and ‘reference’ galaxies are shown in the top, middle, and bottom panels,
respectively. The star formation histories for the disk components are repre-
sented by the solid lines and those for the spheroidal components are indi-
cated by the dotted lines. We also show the mass accretion rates, which are
multiplied by 1000, onto the central BHs in the main progenitors (dashed
lines). Note that while star formation not only in the main progenitors but
also in all other building blocks is counted in the star formation histories,
the mass accretion rates only represent the rates onto the central BHs of the
main progenitors.
mation in the disc component (t ≃ 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 ∼ 8 Gyr
there are many short lived star formation episodes in the spheroids,
which do not correspond to gas accretion events. These star forma-
tion 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 be-
tween 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.1rvir contain the stellar systems defined as main pro-
genitors 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 mas-
sive satellites, both of which ruin the dynamical decomposition.
c© 2007 RAS, MNRAS 000, 1–22
14 T. Okamoto, R. S. Nemmen and R. G. Bower
Table 1. Total and spheroid’s mass, g and r band total luminosities, and disk-to-total mass and luminosity ratios for the galaxies at z = 0.
Mtot/ M⊙ MSP/ M⊙ Mg Mr
no-fb 6.01× 1010 1.98× 1010 -20.8 -21.3 0.67 0.83 0.79
weak-fb 5.98× 1010 1.44× 1010 -20.9 -21.3 0.76 0.88 0.85
reference 4.32× 1010 1.21× 1010 -20.6 -21.0 0.72 0.86 0.83
Figure 11. The relation between mass of the central BHs of the main pro-
genitors and the stellar mass of the host spheroids at z = 3, 2, 1.5, 1, 0.5,
0.2, and 0. The relations for the ‘no-fb’, ‘weak-fb’, and ‘reference’ simula-
tion are shown by the plus signs, diamonds, and triangles respectively. The
solid line indicates MBH = 0.001MSP.
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 res-
olution 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 MBH-MSP relations evolve similarly in all simulations be-
cause 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,
MBH ∼ 10
−3MSP, and then it rapidly moves towards the local
relation during the starburst (z = 2 ∼ 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 opti-
cal depth becomes ∼ 1. After the starburst, the mass of BHs stays
almost constant as seen in Fig. 8. The relation is imprinted in the ra-
diation 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.
Figure 12. Colour-magnitude diagrams for the ‘no-fb’ (left), ‘weak-
fb’ (middle), and ‘reference’ (right) galaxies. The galaxies at z =
3, 2, 1.5, 1, 0.5, 0.2 and 0 are indicated by the plus signs connected by
the dotted line. The colours of their spheroidal components are also indi-
cated by the diamonds. The upper and lower thin solid lines represent mean
colours of the red and blue populations taken from Baldry et al. (2004).
The u − r colour of simulated galaxies are converted to the SDSS colours
by using (u − r)(SDSS) = (u − r)(AB) + 0.05 (Abazajian et al. 2003) to
compare with the data by Baldry et al. (2004).
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 galax-
ies (Baldry et al. 2004) are indicated by solid lines. All galaxies
show similar evolution on the colour-magnitude plane. This con-
firms 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 for-
mation, 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 llink = 0.02 〈l〉, where 〈l〉 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 (LAGN/(Ṁc
2)) is ηSD = 10% for a
standard thin disk and for a RIAF we use the equation
ηRIAF = 0.1
ṁcrit
+ δ + 10
where δ represents the fraction of the turbulent energy which heats
the electrons (Quataert 2001). While the value of δ is quite uncer-
tain, we employ δ = 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
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 15
Figure 13. Upper panels: Bolometric luminosities radiated by the stars in the main progenitors (solid lines) and by the AGNs of the main progenitors (dotted
lines). A 10% radiative efficiency is assumed for standard thin disks and the efficiency given by eq. (26) is assumed for RIAFs. Lower panel: Raw feedback
luminosities ejected by SNe (II and Ia) in the main progenitors (solid lines) and by jets from the central BHs (AGNs) of the main progenitors (dotted lines).
From left to right, the ‘no-fb’, ‘weak-fb’, and ‘reference’ simulations are shown.
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 galax-
ies. Since the radiation drag is more efficient for more compact star
forming region, we expect that the significance of the AGN lumi-
nosity 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 ∼ 1. While all the AGNs
have much lower luminosities than their host galaxies at low red-
shift because of the low accretion rates and the low radiative ef-
ficiency 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 suppres-
sion of star formation around the BH by AGN feedback, makes the
AGN invisible most of the time. In the ‘no-fb’ and ‘weak-fb’ simu-
lations, 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 incorpo-
rate 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 lu-
minosities 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
c© 2007 RAS, MNRAS 000, 1–22
16 T. Okamoto, R. S. Nemmen and R. G. Bower
RIAF modes seen in the ‘weak-fb’ simulation. The jet luminosi-
ties suddenly exceed the SN luminosities when the accretion flow
becomes radiatively inefficient. Since in the ‘no-fb’ and ‘weak-fb’
simulations, the efficiencies, ηAGN by which the AGN jet energy
is given to the surrounding halo gas are small (η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, ηAGN = 0.1 is sufficiently high to suppress the star for-
mation around the BH therefore quenching the mass accretion to
the BH. This is an important feature of our AGN feedback. The im-
pulsive 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 ∼ 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 artifi-
cial 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 ef-
ficiencies. In this sections, we compare simulations which all have
strong AGN feedback, but which differ in other aspects. We com-
pare 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 ηAGN = 1 so that all the jet en-
ergy is used to heat the surrounding ISM because a model with
η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, η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 stel-
lar 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 oth-
ers. 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
Figure 14. Same as Fig. 5 but for the ‘ism-fb’ (left), ‘max-fb’ (middle),
and ‘slow-wind’ (left) models.
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 en-
ergy 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 red-
shift (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 ∼ 6 Gyr than the ‘reference’ galaxy. The
‘slow-wind’ galaxy has a completely different star formation his-
tory from others. Owing to the slow wind speed, it has a higher star
formation rate at high redshift including two burst events at t ∼ 2
and 3 Gyr. The AGN in this galaxy becomes underfed earlier than
in other models (at t ∼ 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 particu-
larly at high redshift (t < 6 Gyr) where the accretion disc is thin
(ṁ > 10−2). In these two models, the accretion flows become RI-
AFs 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 en-
ergy is distributed to the halo gas. In the max-fb model, there is no
accretion event after t ≃ 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 ow-
ing to the slow wind speed, the AGN also becomes underfed earlier
(t ∼ 4 Gyr) than in other models. Contrary to other models, the
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 17
Figure 15. Upper panel: star formation histories for the ‘ism-fb’ (left), ‘max-fb’ (middle), and ‘slow-wind’ (right) models. The star formation history for the
reference model (dashed line in Fig. refSFH) is indicated by the dotted line. In the right panel, we also show the ‘slow-nofb’ model in which we use the same
wind speed as in the ‘slow-wind’ model and do not include AGN feedback (dashed line) for comparison. Lower panel: same as the lower panel of Fig. 8 but
for the ‘ism-fb’ (left), ‘max-fb’ (middle), and ‘slow-wind’ (right) models.
mass accretion continues at t 6 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 pre-
sented 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 ra-
diation drag model, that is, the MBH −MSP 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 ∼ 0.1%
of the mass of their host spheroids, MBH ≃ 10
−3MSP, when the
accretion becomes radiatively inefficient.
6 COMPARISON WITH OBSERVATIONS OF AGN IN
SPIRAL GALAXIES
As we have shown in §5.2, for z . 1 the bolometric luminos-
ity of the AGN is feeble compared to the host galaxy and the nu-
cleus undergoes several recurring episodes of activity. During these
episodes when the nucleus is active, its characteristic accretion rate
is in the range ṁ ≈ 10−4 − 10−3, while the typical nuclear bolo-
metric luminosities and jet powers are Lbol ≈ 10
40−1042 erg s−1
and Ljet ≈ 10
42 − 1043 erg s−1 respectively. Given these prop-
erties, at z . 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
c© 2007 RAS, MNRAS 000, 1–22
18 T. Okamoto, R. S. Nemmen and R. G. Bower
Figure 16. Same as Fig. 11 but for the ‘ism-fb’ (pluses), ‘max-fb’ (dia-
monds), and ‘slow-wind’ (triangles) models.
the time is “off” or inactive with ṁ ≈ 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 cen-
tral 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 . 1 agrees well with the median
observational values for LINERs (median Lbol ≈ 10
41 erg s−1)
and Seyferts (median Lbol ≈ 10
42 erg s−1) inferred from the
Palomar spectroscopic survey of nearby galaxies (Ho 2004). Since
the Seyfert episodes of the central AGN coincide with jet out-
burst events during which the AGN would be seen as radio-loud
in our simulations, this suggests that observed Seyfert/LINER nu-
clei would be mostly radio-loud objects rather than radio-quiet
given their formation history, contrary to the generally held no-
tion 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 nu-
clei are radio-loud, a finding that corroborates our results. In par-
ticular, Ho & Ulvestad (2001) obtained that & 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 ≈
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 asso-
ciated 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, un-
collimated 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 demon-
strate 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 en-
train and heat gas as they flow out of the nucleus, or through heat-
ing of the surrounding gas via thermalisation of the kinetic power
of the jet, generating a wind. Either way, these processes bear re-
semblance to the generation of cavities or bubbles of X-ray emit-
ting 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 there-
fore provide evidence for the possible importance of feedback due
to the radio jets in Seyfert galaxies, which could heat the circum-
nuclear 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 ∼ 1042 − 1043erg 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 . 1 in the
simulation, Ljet ≈ 10
42 − 1043 erg s−1 (Fig. 13). It is worth not-
ing 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 η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 ηAGN. For instance, based on ram pressure and terminal ve-
locity 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 ηAGN
may be considerably smaller than 0.1. The other simulations we
present encompass the uncertainty in the value of ηAGN. In this re-
c© 2007 RAS, MNRAS 000, 1–22
AGN feedback in cosmological simulations 19
spect, the simulation adopting ηAGN = 0.01 does not show any
significant difference from the simulation without AGN feedback
(η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 re-
sults 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 forma-
tion 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 stochas-
tic accretion events. He arrived at this conclusion by consider-
ing 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 cen-
tre. 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 . 1 after the galactic disk is formed, the AGN
is active only during a fraction ≈ 9% of the its life in the ‘reference’
model, ≈ 25% in the ‘weak-fb’ model and ≈ 1% in the ‘ism-
fb’ model. This implies that the lifetime of Seyfert activity in the
spiral galaxy is ≈ 109 yr for the ‘reference’ or ‘weak-fb’ model
and ≈ 108 yr for the ‘ism-fb’ model. This is in agreement with
statistical Seyfert lifetime ∼ 108−109 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 106 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 . 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 ac-
cretion, and AGN feedback in cosmological simulations of galaxy
formation. Our treatment of black hole accretion is based on a the-
oretical 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, mo-
tivated by the recent semi-analytic models of Croton et al. (2006)
and Bower et al. (2006), we have distinguished two distinct accre-
tion modes of BH fuelling, namely standard (optically thick, geo-
metrically 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 surround-
ing diffuse gas, and assume RIAFs are responsible for production
of powerful jets. As the first application of this model, we have per-
formed 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 & 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 ∼ 1 because most of the gas at the galactic cen-
tre is consumed by star formation or blown out by the associated
SN feedback. Afterwards, the AGN fuelling rate drops and the ac-
cretion 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 ∼ 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 gen-
erated, which correspond to the radio-loud state. In-between these
Seyfert episodes the nucleus is turned off and the galaxy is inactive.
Below z ∼ 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 feed-
back 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 sup-
press 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 galax-
ies 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 pre-
dicted AGN lifetime and duration of individual episodes of activ-
ity. 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 ∼ 10−3 regard-
less 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 ac-
cretion 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), al-
though the efficiency of the radiation drag we assumed to ob-
tain 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 mod-
c© 2007 RAS, MNRAS 000, 1–22
20 T. Okamoto, R. S. Nemmen and R. G. Bower
els (Baugh et al. 2005; Nagashima et al. 2005a,b) and cosmologi-
cal 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 popula-
tions 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 in-
terplay 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 re-
produce 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 dis-
cussions. TO acknowledges financial support from the Japan Soci-
ety 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 Eu-
ropean Commissions ALFA-II programme through its funding of
the Latin-american European Network for Astrophysics and Cos-
mology (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).
REFERENCES
Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ,
597, 21
Abazajian K., et al., 2003, AJ, 126, 2081
Allen S. W., Dunn R. J. H., Fabian A. C., Taylor G. B., Reynolds
C. S., 2006, MNRAS, 372, 21
Antonucci R., 1993, ARA&A, 31, 473
Baugh C. M., Lacey C. G., Frenk C. S., Granato G. L., Silva L.,
Bressan A., Benson A. J., Cole S., 2005, MNRAS, 356, 1191
Bicknell G. V., Dopita M. A., Tsvetanov Z. I., Sutherland R. S.,
1998, ApJ, 495, 680
Bı̂rzan L., Rafferty D. A., McNamara B. R., Wise M. W., Nulsen
P. E. J., 2004, ApJ, 607, 800
Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
Bondi H., 1952, MNRAS, 112, 195
Bondi H., Hoyle F., 1944, MNRAS, 104, 273
Booth C. M., Theuns T., Okamoto T., 2007, MNRAS accepted
(astro-ph/0701877)
Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S.,
Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645
Baldry I. K., Glazebrook K., Brinkmann J., Ivezić Ž., Lupton
R. H., Nichol R. C., Szalay A. S., 2004, ApJ, 600, 681
Capetti A., Axon D. J., Macchetto F. D., Marconi A., Winge C.,
1999, ApJ, 516, 187
Cattaneo A., Blaizot J., Devriendt J., Guiderdoni B., 2005, MN-
RAS, 364, 407
Cecil G., et al., 2000, ApJ, 536, 675
Cen R., Ostriker J. P., 1992, ApJL, 399, L113
Churazov E., Sazonov S., Sunyaev R., Forman W., Jones C.,
Böhringer H., 2005, MNRAS, 363, L91
Churazov E., Sunyaev R., Forman W., Böhringer H., 2002, MN-
RAS, 332, 729
Colbert E. J. M., Baum S. A., O’Dea C. P., Veilleux S., 1998, ApJ,
496, 786
Colbert E. J. M., Baum S. A., Gallimore J. F., O’Dea C. P., Chris-
tensen J. A., 1996, ApJ, 467, 551
Cowie L. L., McKee C. F., 1977, ApJ, 211, 135
Croton D. J., Springel V., White S. D. M., De Lucia G., Frenk
C. S., Gao L., Jenkins A., Kauffmann G., Navarro J. F., Yoshida
N., 2006, MNRAS, 365, 11
Dalla Vecchia C., 2005, PhD. thesis
Dalla Vecchia C., Bower R. G., Theuns T., Balogh M. L., Maz-
zotta P., Frenk C. S., 2004, MNRAS, 355, 995
Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ,
292, 371
De Lucia G., Springel V., White S. D. M., Croton D., Kauffmann
G., 2006, MNRAS, 366, 499
Di Matteo T., Croft R. A. C., Springel V., Hernquist L., 2003, ApJ,
593, 56
Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., Sijacki, D.,
2007, MNRAS, submitted (arXiv:0705.2269)
Elmegreen B. G., 1989, ApJ, 338, 178
Enoki M., Nagashima M., Gouda N., 2003, PASJ, 55, 133
Escala A., Larson R. B., Coppi P. S., Mardones D., 2004, ApJ,
607, 765
Fabian A. C., Sanders J. S., Taylor G. B., Allen S. W., Crawford
C. S., Johnstone R. M., Iwasawa K., 2006, MNRAS, 366, 417
Falcke H., Wilson A. S., Simpson C., 1998, ApJ, 502, 199
Ferrarese L., Merritt D., 2000, ApJL, 539, L9
Fioc M., Rocca-Volmerange B., 1997, A&A, 326, 950
Fukui Y., Mizuno N., Yamaguchi R., Mizuno A., Onishi T., 2001,
PASJ, 53, L41
Gallimore J. F., Axon D. J., O’Dea C. P., Baum S. A., Pedlar A.,
2006, AJ, 132, 546
Gebhardt K., Bender R., Bower G., Dressler A., Faber S. M., Fil-
ippenko A. V., Green R., Grillmair C., Ho L. C., Kormendy J.,
Lauer T. R., Magorrian J., Pinkney J., Richstone D., Tremaine
S., 2000, ApJL, 539, L13
Governato F., Willman B., Mayer L., Brooks A., Stinson G.,
Valenzuela O., Wadsley J., Quinn T., 2007, MNRAS, 374, 1479
Granato G. L., De Zotti G., Silva L., Bressan A., Danese L., 2004,
ApJ, 600, 580
Haardt F., Madau P., 1996, ApJ, 461, 20
Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103
Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H.,
2003, ApJ, 591, 288
Heyer M. H., Carpenter J. M., Snell R. L., 2001, ApJ, 551, 852
Ho L. C., Fillipenko A. V., Sargent W. L. W. 1997, ApJ, 487, 568
Ho L. C. 2004, in Coevolution of Black Holes and Galaxies, ed.
L. C. Ho (Cambridge: Cambridge Univ. Press), 293
Ho L. C., Ulvestad J. S. 2001, ApJ, 555, 650
c© 2007 RAS, MNRAS 000, 1–22
http://arxiv.org/abs/astro-ph/0701877
http://arxiv.org/abs/0705.2269
AGN feedback in cosmological simulations 21
Hopkins P. F., Hernquist L., 2006, ApJS, 166, 1
Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., Robertson
B., Springel V., 2006, ApJS, 163, 1
Hoyle F., Lyttleton R. A., 1939, in Proceedings of the Cambridge
Philisophical Society The effect of interstellar matter on climatic
variation. p. 405
Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
Kauffmann G., Haehnelt M., 2000, MNRAS, 311, 576
Kaufmann T., Mayer L., Wadsley J., Stadel J., Moore B., 2007,
MNRAS, 375, 53
Kawakatu N., Saitoh T. R., Wada K., 2005, ApJ, 628, 129
Kawakatu N., Umemura M., 2002, MNRAS, 329, 572
Kawakatu N., Umemura M., 2004, ApJ, 601, L21
Kawata D., Gibson B. K., 2005, MNRAS, 358, L16
Keel W. C., White R. E., III, Owen F. N., Ledlow M. J., 2006, AJ,
132, 2233
Kennicutt Jr. R. C., 1998, ApJ, 498, 541
Kinney A. L., Schmitt H. R., Clarke C. J., Pringle J. E., Ulvestad
J. S., Antonucci R. R. J., 2000, ApJ, 537, 152
Kharb P., O’Dea C. P., Baum S. A., Colbert E. J. M., Xu C., 2006,
ApJ, 652, 177
Kitzbichler M. G., White S. D. M., 2007, MNRAS, pp 75–+
Kormendy J., Richstone D., 1995, ARA&A, 33, 581
Kravtsov A. V., 2003, ApJL, 590, L1
Krolik J. H., 1999, in Krolik J. H., ed., Active Galactic Nu-
clei: From the Central Black Hole to the Galactic Environment.
Princeton Univ. Press, Princeton, NJ, c1999
Kukula M. J., Ghosh T., Pedlar A., Schilizzi R. T., Miley G. K.,
de Bruyn A. G., Saikia D. J., 1993, MNRAS, 264, 893
Kukula M. J., Holloway A. J., Pedlar A., Meaburn J., Lopez J. A.,
Axon D. J., Schilizzi R. T., Baum S. A., 1996, MNRAS, 280,
Laor A., 2000, ApJ, 543, L111
Livio M., Ogilvie G. I., Pringle J. E., 1999, ApJ, 512, 100
Maccarone T. J., Gallo E., Fender R., 2003, MNRAS, 345, L19
Magorrian J., Tremaine S., Richstone D., Bender R., Bower G.,
Dressler A., Faber S. M., Gebhardt K., Green R., Grillmair C.,
Kormendy J., Lauer T., 1998, AJ, 115, 2285
Makino J., Funato Y., 2004, ApJ, 602, 93
Malbon R. K., Baugh C. M., Frenk C. S., Lacey C. G., 2007, MN-
RAS (in press)
Marconi A., Hunt L. K., 2003, ApJL, 589, L21
McLure R. J., Dunlop J. S., 2002, MNRAS, 331, 795
Meier D. L., 2001, ApJL, 548, L9
Merritt D., Ferrarese L., 2001, ApJ, 547, 140
Middelberg E., Agudo I., Roy A. L., Krichbaum T. P., 2007, MN-
RAS, 377, 731
Monaco P., Fontanot F., 2005, MNRAS, 359, 283
Morganti R., Oosterloo T., Tsvetanov Z., 1998, AJ, 115, 915
Mundell C. G., Wrobel J. M., Pedlar A., Gallimore J. F., 2003,
ApJ, 583, 192
Nagamine K., Springel V., Hernquist L., 2004, MNRAS, 348, 421
Nagashima M., Lacey C. G., Baugh C. M., Frenk C. S., Cole S.,
2005a, MNRAS, 358, 1247
Nagashima M., Lacey C. G., Okamoto T., Baugh C. M., Frenk
C. S., Cole S., 2005b, MNRAS, 363, L31
Narayan R., 2005, Ap&SS, 300, 177
Narayan R., Mahadevan R., Quataert E., 1998, in Abramowicz
M. A., Bjornsson G., Pringle J. E., eds, Theory of Black Hole
Accretion Disks Advection-dominated accretion around black
holess. p. 148
Narayan R., Yi I., 1995, ApJ, 444, 231
Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
Navarro J. F., White S. D. M., 1993, MNRAS, 265, 271
Nemmen R. S., Bower R. G., Babul A., Storchi-Bergmann T.,
2007, MNRAS, 377, 1652
Nemmen R. S., Storchi-Bergmann T., Yuan F., Eracleous M.,
Terashima Y., Wilson A. S., 2006, ApJ, 643, 652
Novikov I. D., Thorne K. S., 1973, in DeWitt C., DeWitt B. S.,
eds, Black Holes. New York: Gordon & Breach, p. 344
Okamoto T., Eke V. R., Frenk C. S., Jenkins A., 2005, MNRAS,
363, 1299 (OEFJ05)
Okamoto T., Jenkins A., Eke V. R., Quilis V., Frenk C. S., 2003,
MNRAS, 345, 429
Omma H., Binney J., Bryan G., Slyz A., 2004, MNRAS, 348,
Quataert E., 2001, in Peterson B. M., Pogge R. W., Polidan R. S.,
eds, ASP Conf. Ser. 224: Probing the Physics of Active Galactic
Nuclei Low Radiative-Efficiency Accretion Flows. p. 71
Quilis V., Bower R. G., Balogh M. L., 2001, MNRAS, 328, 1091
Rafferty D. A., McNamara B. R., Nulsen P. E. J., Wise M. W.,
2006, ApJ, 652, 216
Rees M. J., Phinney E. S., Begelman M. C., Blandford R. D.,
1982, Nature, 295, 17
Riffel R. A., Storchi-Bergmann T., Winge C., Barbosa F. K. B.,
2006, MNRAS, 373, 2
Rupke D. S., Veilleux S., Sanders D. B., 2005, ApJ, 632, 751
Salpeter E. E., 1955, ApJ, 121, 161
Samland M., Gerhard O. E., 2003, A&A, 399, 961
Sanders R. H., 1984, A&A, 140, 52
Sato J., Umemura M., Sawada K., Matsuyama S., 2004, MNRAS,
354, 176
Scannapieco C., Tissera P. B., White S. D. M., Springel V., 2006,
MNRAS, 371, 1125
Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
Schmitt H. R., Ulvestad J. S., Antonucci R. R. J., Kinney A. L.,
2001, ApJS, 132, 199
Sijacki D., Springel V., 2006, MNRAS, 366, 397
Sijacki D., Springel V., di Matteo T., Hernquist L., 2007, MNRAS,
380, 877
Solomon P. M., Rivolo A. R., 1989, ApJ, 339, 919
Soltan A., 1982, MNRAS, 200, 115
Spitzer L. J., 1962, Physics of Fully Ionized Gases. New York:
Interscience
Spitzer L. J., 1978, Physical Processes in the Interstellar Medium.
New York: Wiley-Interscience
Springel V., 2005, MNRAS, 364, 1105
Springel V., Di Matteo T., Hernquist L., 2005a, ApJL, 620, L79
Springel V., Di Matteo T., Hernquist L., 2005b, NNRAS, 361, 776
Springel V., Hernquist L., 2003, MNRAS, 339, 289 (SH03)
Springel V., Yoshida N., White S. D. M., 2001, NewA, 6, 79
Steinmetz M., Mueller E., 1994, A&A, 281, L97
Taylor G. B., Sanders J. S., Fabian A. C., Allen S. W., 2006, MN-
RAS, 365, 705
Terashima Y., Wilson A. S., 2003, ApJ, 583, 145
Theuns T., Schaye J., Zaroubi S., Kim T.-S., Tzanavaris P., Car-
swell B., 2002, ApJL, 567, L103
Tremaine S., Gebhardt K., Bender R., Bower G., Dressler A.,
Faber S. M., Filippenko A. V., Green R., Grillmair C., Ho L. C.,
Kormendy J., Lauer T. R., Magorrian J., Pinkney J., Richstone
D., 2002, ApJ, 574, 740
Umemura M., 2001, ApJL, 560, L29
Volonteri M., Madau P., Quataert E., Rees M. J., 2005, ApJ, 620,
c© 2007 RAS, MNRAS 000, 1–22
22 T. Okamoto, R. S. Nemmen and R. G. Bower
Wada K., 2001, ApJL, 559, L41
Wada K., Norman C. A., 2007, ApJ, 660, 276
Wada K., Norman C. A., 2001, ApJ, 547, 172
Wang J.-M., Chen Y.-M., Ho L. C., McLure R. J., 2006, ApJL,
642, L111
Whittle M., Wilson A. S., 2004, AJ, 127, 606
Wilson A. S., Ulvestad J. S., 1983, ApJ, 275, 8
Woltjer L., 1968, AJ, 73, 914
Yepes G., Kates R., Khokhlov A., Klypin A., 1997, MNRAS, 284,
Yu Q., Tremaine S., 2002, MNRAS, 335, 965
Yuan F., 2007, in ASP Conf. Ser. 373, The Central Engine of Ac-
tive Galactic Nuclei, ed. L. C. Ho & J.-M. Wang (San Francisco:
ASP), 95
c© 2007 RAS, MNRAS 000, 1–22
	Introduction
	The model of star-forming gas and stellar feedback
	Gas cooling
	Cloud formation and the cloud model
	Star formation
	Cloud evaporation
	Supernova feedback and chemical enrichment
	Numerical implementation
	Winds
	Parameter setting and test simulations
	Black hole accretion and AGN feedback
	Seed black holes
	Mass accretion by the radiation drag
	AGN feedback
	Cosmological simulations of disc galaxy formation with AGN feedback
	Black hole particles
	Models
	Results
	Simulations with different AGN feedback efficiencies
	Relative importance of the AGN
	Dependence on feedback parameters and implementations of feedback
	Comparison with observations of AGN in spiral galaxies
	Summary and Conclusions

  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)

Introduction
	The model of star-forming gas and stellar feedback
	Gas cooling
	Cloud formation and the cloud model
	Star formation
	Cloud evaporation
	Supernova feedback and chemical enrichment
	Numerical implementation
	Winds
	Parameter setting and test simulations
	Black hole accretion and AGN feedback
	Seed black holes
	Mass accretion by the radiation drag
	AGN feedback
	Cosmological simulations of disc galaxy formation with AGN feedback
	Black hole particles
	Models
	Results
	Simulations with different AGN feedback efficiencies
	Relative importance of the AGN
	Dependence on feedback parameters and implementations of feedback
	Comparison with observations of AGN in spiral galaxies
	Summary and Conclusions
