DRIVING OUTFLOWS WITH RELATIVISTIC JETS AND THE DEPENDENCE OF AGN FEEDBACK ...

26 Pages·2012·6.37 MB·English

Submitted to the Astrophysical Journal on April 23, 2012.
PreprinttypesetusingLATEXstyleemulateapjv.11/10/09
DRIVING OUTFLOWS WITH RELATIVISTIC JETS AND THE DEPENDENCE OF AGN FEEDBACK
EFFICIENCY ON ISM INHOMOGENEITY
A. Y. Wagner1, G. V. Bicknell2, M. Umemura1
Submitted to the Astrophysical Journal on April 23, 2012.
ABSTRACT
We examine the detailed physics of the feedback mechanism by relativistic AGN jets interacting
with a two-phase fractal interstellar medium (ISM) in the simulations by Wagner & Bicknell (2011).
We extend the parameter space study of feedback efficiencies to 29 3D grid-based hydrodynamical
simulations with 15 new simulations of systems with low volume factors and varying maximum cloud
sizes on a grid of up to 10243 cells with a resolution down to 1pc per cell. We compare the expan-
sion of the bubbles driven by the jet plasma to that of ideal energy-driven bubbles, and reveal the
mechanism by which dense clouds are accelerated to bubble expansions speeds within the dynamical
time of the bubble. The (negative) feedback efficiency, as measured by the amount of cloud-dispersal
generated by the jet-ISM interactions, is sensitive to the maximum size of clouds in the fractal cloud
distribution;foragivenfillingfactoranddensity,distributionsofsmallercloudsleadtohigheroutflow
velocities. Feedback ceases to be efficient for Eddington ratios P /L (cid:46) 10−4, although systems
jet edd
with large cloud complexes (cid:38) 50pc require jets of Eddington ratio in excess of 10−2 to disperse the
clouds appreciably. Compared to the ISM density and maximum cloud size, the feedback efficiency
dependsweaklyonvolumefillingfactor. Basedonmeasurementsofthebubbleexpansionratesinour
simulations we argue that sub-grid AGN prescriptions resulting in negative feedback in cosmological
simulationswithoutamulti-phasetreatmentoftheISMaregoodapproximationsifthevolumefilling
factor of warm phase material is less than 0.1 and the cloud complexes are smaller than ∼25pc. We
find that the acceleration of the dense embedded clouds is provided by the ram pressure (rather than
the thermal pressure) of the high velocity flow through the porous channels of the warm phase, flow
that has fully entrained the shocked hot-phase gas it has swept up, and is additionally mass-loaded
byablatedcloudmaterial. Thismechanism,reminiscentofatwo-stagefeedbackscenarioproposedby
Hopkins&Elvis(2010),transfers10%to40%ofthejetenergytothecoldandwarmgas,accelerating
it to several 100 to several 1000 kms−1 within a few 10 to 100 Myr. Our predicted velocities match
thoseobservedinarangeofhighandlowredshiftradiogalaxieshostingpowerfulradiojets. Weusea
synthetic radio map generated from one of our models to explain some of the morphological features
of the radio loud quasar and ULIRG 3C 48.
Subject headings: galaxies: evolution – galaxies: formation – galaxies: jets – hydrodynamics – ISM:
jets and outflows – methods: numerical
1. INTRODUCTION trolstheaccretionrateofmatterintothecentralregions
(Umemura 2001; Kawakatu & Umemura 2002). Cosmo-
The formation of galaxies is a non-linear, but to some
logical SPH and semi-analytic models routinely include
degree self-regulatory process; the star-formation effi-
feedback by powerful radio jets or quasar winds, albeit,
cienciesofgalaxiesandthegrowthrateofthecentralsu-
of necessity, using highly simplified models for the feed-
permassive black-holes (SMBH) are thought to be mod-
back. Observationally, ionization diagnostics may not
ified by feedback processes from active galactic nuclei
conclusively distinguish the contributions of radiatively
(AGN) resulting in a tight correlation between SMBH
driven feedback and feedback driven by jet-ISM interac-
mass and the bulge stellar velocity dispersion (the M–σ
tions Holt et al. (2009), although in some cases jet-ISM
relation,Ferrarese&Merritt2000;Gebhardtetal.2000;
interactions are strongly favoured (Dopita et al. 1997;
Tremaine et al. 2002).
Nesvadbaetal.2010). Severalstudiesfindstatisticalcor-
Itisunclear,however,whichtypesofAGNactivityare
relations between AGN activity, outflows, and the sup-
relevantinregulatingbulgeandSMBHgrowth. TheSilk
pression of star-formation (e.g. Schawinski et al. 2007;
& Rees (1998) model invokes an energy-driven quasar
Farrah et al. 2012), but the connection between AGN
windoflowEddingtonratio,whilethemodelsbyFabian
jetsandstar-formationremainsambiguous(Dickenetal.
(1999); King (2003), and Murray et al. (2005) consider
2012).
opacity-regulated momentum-driven outflows requiring
In cosmological SPH simulations (e.g. Okamoto et al.
Eddington ratios of a few percent. Another possibility
2008; Di Matteo et al. 2008; Schaye et al. 2010), grid-
is that the radiation field in the bulge of galaxies con-
based simulations (e.g. Springel 2011; Dubois et al.
2012),andsemi-analyticmodels(e.g.Crotonetal.2006;
[email protected]
1Center for Computational Sciences, University of Tsukuba, Fanidakis et al. 2012) AGN feedback is found to be a
1-1-1Tennodai,Tsukuba,Ibaraki,305-8577,Japan necessary ingredient in order to reproduce the observed
2Research School of Astronomy and Astrophysics, The Aus- galaxy luminosity function and its evolution with red-
tralianNationalUniversity,ACT2611,Australia
2 Wagner, Bicknell, & Umemura
shift, but the relevant range of powers varies between also showed that isotropization of the injected energy
models. Cosmological SPH simulations require energy is even more effective when the jet encounters inhomo-
injection rates described by Eddington ratios η (cid:38) 10−2 geneities because, by virtue of its lightness – the jet par-
while some semi-analytic models find that low-powered ticledensityistypically6to8ordersofmagnitudelower
injection of energy with Eddington ratios of η (cid:38)10−5 is thanthatofISMclouds–thejetisstronglydeflectedby
sufficient. In both methods there exist a variety of “sub- the inhomogeneities. Additional effects, e.g., jet insta-
grid”prescriptionstodepositenergyyieldingdifferentre- bilities and jet precession increase the isotropy of energy
sults. Neither method resolves or treats the galaxy-scale deposition, but are not essential.
physics ofthe interaction ofthe outflows andinterstellar In previous work (WB11) we used grid-based hydro-
medium(ISM)adequately,andoneofouraimsistopro- dynamic simulations to model jet-ISM interactions and
vide a a robust description of sub-grid feedback physics quantifiedthefeedbackefficiencyprovidedbyrelativistic
thatcanbeusedinfuturesemi-analyticandcosmological AGN jets in the core of young, gas-rich radio galaxies.
models. The simulated galaxies typically represent either Com-
Feedback involving mechanical energy input by an pact Steep Spectrum (CSS) or Gigahertz Peaked Spec-
AGN jet, often termed “radio-mode” feedback, has been trum (GPS) sources (Bicknell et al. 1997), which in our
identified as a key mechanism to heat the IGM of the viewareaclassofobjectsexperiencinganearlyphaseof
cluster and prevent a runaway build-up of galaxy mass powerful jet-mediated feedback. In these objects, radio
through further accretion of cooling gas (see Best et al. sourceexpansionisimpededbythedensemulti-phaseen-
2006,andreferencestherein). Well-studiednearbyexam- vironment of the galaxy core in the early phase of their
ples include the Hydra A (Wise et al. 2007), Perseus A evolution. WeconcludedthatAGNjetfeedbackinthese
(Fabian et al. 2006), and M87 in the Virgo cluster (Mil- systems is effective in all galaxies for jets with powers
lionetal.2010),andthephenomenoniswellreproduced 1043 – 1046ergs−1 if the ratio of jet power to Eddington
in cluster-scale grid-based hydrodynamic simulations by luminosity η (cid:38)10−4.
Dubois et al. (2010, 2011) and Teyssier et al. (2011). A unique feature of these simulations is the treatment
Galaxy-scale jet-regulated star-formation (“positive” of the galaxy ISM with a two-point fractal, single point
feedback)maybeveryrelevantathigherredshiftsingas- log-normal warm-phase distribution (clouds) embedded
richgalaxiesandproto-galacticenvironments(DeYoung in a hot atmosphere. We determined feedback efficien-
1989; Bicknell et al. 2000; Reuland et al. 2003; Klamer cies as a function of some of the parameters describing
et al. 2004; Miley et al. 2006; Villar-Mart´ın 2007; Vene- this distribution, e.g. the density and filling factor of
mansetal.2007;Miley&DeBreuck2008). Thecasefor the warm gas. The ISM properties in HzRG are uncer-
jet-induced star-formation in disc galaxies was made in tain; while large reservoirs of molecular gas and HI are
numericalworkasearlyasWoodward(1976)andinmore known to exist, the volume filling factors of the cold and
recent simulations by Fragile et al. (2005) and Gaibler warm gas and the typical sizes of clouds are not known.
et al. (2011). WB11 restricted themselves to volume filling factors of
In some nearby and high-redshift radio galaxies the clouds of 0.42 and 0.13, which are probably at the
(HzRG) neutral and line-emitting gas is seen outflow- higher end of the range of “typical” values. Further-
ingatseveral100kms−1 toseveral1000kms−1 (Gelder- more, we did not investigate the dependence on maxi-
man&Whittle1994;Tadhunteretal.2001;O’Deaetal. mum cloud sizes.
2002;Emontsetal.2005;Holtetal.2008,2011;Morganti With the 15 new simulations presented in this paper,
etal.2005,2007,2010;Nesvadbaetal.2006,2007,2008, wehavenowsubstantiallyextendedthisparameterspace
2010; Lehnert et al. 2011; Guillard et al. 2012; Torresi study to lower filling factors and a variety of maximum
etal.2012). Thealignmentofthejetwithoutflowinggas cloud sizes in the fractal distribution. We also examine
(Pentericcietal.2001;Privonetal.2008),andmatching the acceleration mechanism in more detail, providing an
energetics (Nesvadba et al. 2006, 2007) suggest that the explanation for the high mechanical advantage observed
outflowsaredrivenbythetransferofenergyandmomen- byWB11. Wedescribeourmethodsofcomputationand
tum from the jet to to the dense ISM. This hypothesis parameterspacenextin§2and§3,andpresentourmain
is supported by our previous 3D hydrodynamical simu- results in detail in §4. In §5, we compare our simulation
lations of AGN-jet driven outflows (Wagner & Bicknell results with data from a sample of radio galaxies with
2011, WB11 henceforth). observed outflows. We also discuss other feedback cri-
The question of how a collimated jet may impart en- teria and the review the difficulties in modelling cloud
ergy and momentum isotropically, e.g., to affect the en- ablation. We conclude with a summary of the paper in
tire volume defining the bulge of a galaxy, is frequently §6.
mentioned (De Young 2010; Ostriker et al. 2010). A re-
2. EQUATIONSANDCODE
lated problem is the momentum budget associated with
The system of equations describing the relativistic jet
thedispersionorexpulsionofcloudsinagalaxy. Anim-
plasma, hot atmosphere, and warm clouds in the one
portant feature of AGN jets is that the jet is extremely
fluid approximation is (Landau & Lifshitz 1987):
light and that jet and cocoon are highly overpressured
(underexpanded) with respect to the ambient environ- ∂D ∂Dui
ment (Begelman & Cioffi 1989). Simulations by Sax- + =0; D =Γρ;
∂t ∂xi
ton et al. (2005), Sutherland & Bicknell (2007), Gaibler
etal.(2009),andWB11ofAGNjetsshowhowthelight, ∂Fi + ∂Fiuj + ∂p =0; Fi =ρwΓ2ui/c2 ; (1)
overpressured jet inflates a cocoon that drives a quasi- ∂t ∂xj ∂xi
spherical energy bubble into the ISM. These simulations ∂E ∂Fic2
+ =−ρ2Λ(T); E =ρwΓ2−p.
∂t ∂xi
Relativistic Jet Feedback 3
Theconservedquantities,D,Fi,andE arethelabora- where
tory frame fluid density, components of the momentum (cid:115)
density,totalenergydensity(includingtherestmassen- µ2 (cid:18)σ2 (cid:19)
m=ln , s= ln +1 , (3)
ergy density). The variables p, T, Λ, and ui are pres- (cid:112)σ2+µ2 µ2
sure, temperature, cooling rate, and the components of
the three velocity, respectively. The bulk Lorentz factor andµandσ2 arethemeanandvarianceofthelognormal
is Γ=(cid:0)1−u ui/c2(cid:1)−1/2. The proper rest frame density distribution.
i
Let F(k) be the Fourier transform of the spatial den-
is ρ and w =c2+pγ/ρ(γ−1) the proper rest frame spe-
sity distribution, ρ(r), with k and r as wave vector and
cific enthalpy for an ideal polytropic equation of state,
positionvector,respectively. Thetwo-pointfractalprop-
with index γ.
erty is characterized in Fourier space by a power spec-
We integrate these equations using the publicly avail-
trum, D(k), in wave number, k, that obeys a power-law
able, open-source Eulerian Godunov-type code FLASH
with index −5/3 for a Kolmogorov-type spectrum
(Fryxell et al. 2000) version 3.2 to which we have added
code to incorporate radiative cooling of thermal gas and (cid:90)
code to advance advected scalars in the relativistic hy- D(k)= k2F(k)F∗(k)dΩ∝k−5/3 , (4)
drodynamic solver.
We exploit the adaptive mesh capabilities of FLASH, where the integral of the spectral density, F(k)F∗(k), is
utilizing up to seven levels of grid refinement in a cubi- over all solid angle, Ω.
cal simulation domain of 1kpc3 in physical dimensions, A cube of random numbers that simultaneously satis-
consistingof10243 cellsatamaximumspatialresolution fies Eqn. (2) and Eqn. (4) is generated by the method
1pc. This is twice the resolution of the simulations by outlined in Lewis & Austin (2002). First, a cube with
WB11 and is necessary in order to capture the fractal cell values from a Gaussian distribution with mean m
outlinesofcloudsforsmallfillingfactorsandcloudsizes. and standard deviation s is Fourier-transformed and
Note that a restricted one parameter scaling of physical apodized by a Kolmogorov power law spectrum in
dimensions is possible (Sutherland & Bicknell 2007). wavenumber with index −5/3 and minimum sampling
Tracer variables distinguish jet material and warm wavenumber kmin. The minimum sampling wavenum-
phase gas from each other and from the hot phase ber is, effectively, the average number of clouds per di-
background. We include non-equilibrium, optically thin mension divided by 2, and it determines the scale of the
atomic cooling for T > 104K (Sutherland & Dopita largest fractal structures in the cube relative to the size
1993) and updated solar abundances (Asplund et al. ofthecube. Forexample,ifkmin =20foracubemapped
2005), for which the mean mass per particle, µ = to a domain of extent 1kpc, then k = 20kpc−1 and
m min
0.6165. Thermal conduction, photo-evaporation, self- the largest structures (clouds) would have extents of
gravity, and magnetic fields are not included. R = 1/(2k ) = 25pc. The cube is then trans-
c,max min
WeranoursimulationsontheNationalComputational formed back into real space and exponentiated. Be-
Infrastructure National Facility (NCI NF) Oracle/Sun cause the last step alters the power-law structure in
Constellation Cluster, a high-density integrated system Fourier space, the cube is iteratively transformed be-
of 1492 nodes of Sun X6275 blades, each containing two tween Fourier space and real space until successive cor-
quad-core 2.93GHz Intel Nehalem cpus, and four inde- rections produce a power-law convergence within 1%.
pendent SUN DS648 Infiniband switches.3 We typically To place the fractal cube into the simulation domain
used 256 to 1024 cpus with 3GB of memory per core to it is apodized (in real space) by a spherically symmetric
complete one simulation within two weeks. mean density profile which in the simulations presented
here is flat with mean warm phase density (cid:104)n (cid:105). The
w
3. MODELPARAMETERS,INITIALCONDITIONS,AND porosity of the warm phase arises by imposing an upper
BOUNDARYCONDITIONS
temperature cutoff for the existence of clouds at T =
crit
A crucial ingredient in the simulations described in 3 × 104K, beyond which clouds are deemed thermally
§3 is the two-phase ISM, which consists of a warm unstable. No lower temperature limit is enforced, and
(T ∼ 104K) phase and a hot (T ∼ 107K) phase. In temperatures in the core of clouds may initially be less
particular, we are concerned with the effect of the jet than 100K. The upper temperature cutoff corresponds
plasma on the state and dynamics of the warm phase directly to a lower density cutoff, ρ =µ p/(kT ), if
crit m crit
material. We have, therefore, extended our studies of the pressure, p, is defined. Here, µ is the mean mass
m
parameters related to the warm-phase and identified the per particle of the hot phase. In our simulations the
correctphysicalmechanismthatleadstotheacceleration clouds are in pressure equilibrium with the surrounding
of the clouds. hot phase, thus ρ = µ n T /T , where n and T
crit m h h crit h h
ThewarmphaseISMdensityisinitializedfromacube are the hot phase number density and temperature, re-
of random numbers that simultaneously satisfies single- spectively. The filling factor of the warm phase, within
pointlognormalstatisticsandtwo-pointfractalstatistics. the hemispherical region of radius 0.5kpc, in which it is
Let P(ρ) be the lognormal probability density function distributed, is:
of the random variable ρ, representing density:
(cid:90) ∞
1 (cid:18)−(lnρ−m)2(cid:19) fV = P(ρ)dρ
P(ρ)= s√2πρexp 2s2 , (2) ρcrit  (cid:110) (cid:112) (cid:111)
ln (ρ /µ) σ2/µ2+1
3 For details of the system specifications see 1 crit
http://nf.nci.org.au/facilities/vayu/hardware.php. =21+erf (cid:112)2ln(σ2/µ2+1)  (5)
4 Wagner, Bicknell, & Umemura
Theoriginalfractalcubewasconstructedwithµ=1and ton ratio) to be η = P /L , and φ , ρ, and v as
jet edd w r
σ = 5, and after apodization with a spatially uniform the warm phase tracer (mass fraction in a cell), density,
mean warm phase density distribution, the single point and radial velocity, respectively. A convenient measure
densitydistributionremainslognormal,butwithamean of the efficiency of feedback is the density-averaged ra-
µwh=ose(cid:104)ntwem(cid:105).peFraotruraeninisotthhiesrmwaolrkhoist fipxheadseatdiTsthrib=ut1io0n7,, dreialaltoivuetfltoowthveelvoecliotcyi,ty(cid:104)vdr,iwsp(cid:105)e=rsi(cid:80)onNl=o1fφawgρalvarx/y(cid:80)’sNlb=u1lφgewaρs,
ρcrit/µ = (Th/Tcrit)(nh/(cid:104)nw(cid:105)) is constant everywhere. predicted by the M–σ relation (Silk & Rees 1998; King
The filling factor is, therefore, directly defined by the 2005). Defining the black hole mass in terms of the Ed-
ratioofhotphasedensityandmeanwarmphasedensity, dington ratio M = 4πGm P c/ησ and using the
BH p jet T
nh/(cid:104)nw(cid:105). For a more detailed description of the method M–σ relationfoundbyTremaineetal.(2002)weexpress
togeneratethefractalcubeandadiscussionofthechoice the velocity dispersion as
of statistical parameters for the lognormal probability
distribution and wavenumber power law index we refer σ =1.0η−1/4P1/4 (6)
the reader to the manuscript by Lewis & Austin (2002) 100 jet,45
and the relevant sections and appendixes in Sutherland whereP isthejetpowerinunitsof45ergs−1. When
jet,45
& Bicknell (2007).
(cid:104)v (cid:105) > σ, the jet-ISM interactions result in sufficient
r,w
The general setup, initial conditions, and boundary
feedback of momentum and energy to establish a highly
conditions used here are identical to those of WB11.
dispersed distribution of cold and warm gas within the
WB11 performed 14 simulations of AGN jets with pow-
core of the galaxy.
ers in the range 43 < log(P /ergs−1) < 46. The
jet Themaximumvaluesof(cid:104)v (cid:105)duringtherunasafunc-
r,w
choiceofwarmphasefillingfactors,f ,of0.42and0.13,
V tion of jet power for all 28 simulations that include the
was relatively high and the maximum cloud size fixed at
warmphaseareshowninFig.3, whichupdatesFig.5in
Rc,max ∼25pc (kmin =20kpc−1). WB11. Points of constant hot phase density and filling
In the 15 new simulations presented here, we explore factorareconnectedalongincreasingjetpowerwithlines
new regions of parameter space with filling factors, fV, of specific colors and style according to the legend. The
of0.052and0.027,correspondingtoaveragewarmphase slanted grey lines in Fig. 3 represent the loci of the ve-
densitiesof150cm−3and100cm−3,andk of40kpc−1 locitydispersionalonglinesofconstantη (asdetermined
min
and 10kpc−1, corresponding to maximum cloud sizes of byEqn.6). Thedashedgreylinesinthefigurerepresent
R = 12.5pc and R ∼ 50.0pc. The range in a different M–σ relation of M ∝ σ5, found more re-
c,max c,max BH
jet power and other parameters defining the jet plasma cently by Graham et al. (2011). The locus in this case
remain the same. These jets typically have a density is σ =1.2η−1/5P1/5 . Using either relation, one may
contrast of 10−4 with respect to the ambient hot phase 100 jet,45
then compare the points of the simulations with values
and 10−7 with respect to the embedded clouds. The
forthevelocitydispersionpredictedbytheM–σ relation
pressurecontrastbetweenthejetandthehotphaseISM
for a given value of η. If a point lies above an isoline for
is typically 102 – 103. AGN jets are extremely light,
η, then feedback by a jet of that power P in a galaxy
underexpanded (overpressured) jets. jet
with Eddington limit P /η is effective. Conversely if
The complete list of 29 simulations, including those jet
a point lies below an isoline for η, then feedback is not
from WB11 are given in Table 1. New runs are marked
effective. Equivalently, the point itself marks a critical
with “(cid:73)”.
value of η, η = (P /L ) , below which feedback
crit jet edd crit
ceases to be effective in galaxy with a SMBH of mass
4. RESULTS
M =4πGm P c/ησ .
4.1. Velocities of accelerated clouds and feedback BH p jet T
As observed in previous simulations, the velocities at-
efficiency
tained by clouds match those observed of outflows in
We conducted 15 new simulations to study new re- radio galaxies (Morganti et al. 2005; Holt et al. 2008;
gions in the space spanned by parameters that describe Nesvadba et al. 2006, 2008, 2010; Lehnert et al. 2011;
the distribution of warm phase material in our simula- Guillard et al. 2012; Torresi et al. 2012). The dense
tions as described in §3. Figure 1 shows density maps cores of the clouds in our simulations are accelerated to
ofthreeselectednewsimulationswithlowerfillingfactor a few 100kms−1, while the diffuse ablated material is
and differing maximum cloud sizes to those of previous accelerated to several 1000kms−1. We discovered that
simulations. Toobtainathreedimensionalimpressionof thefeedbackefficiencyoftherelativisticjetonthewarm
the interactions between the jet and the clouds we show phase ISM increases with increasing jet power, decreas-
a volume render of the density of both components from ing mean ISM density, and increasing filling factor, al-
one of our simulations in Figure. 2. The jet plasma is though only two values for the filling factor, f = 0.42,
V
textured in bluish green and the clouds in purple. The and f =0.13, were studied.
V
forward shock outlining the jet-blown energy bubble is Within the new range of parameter space, the main
seen in a translucent grey. An oval excavation is made conclusions reached in WB11 remain valid; feedback is
inthevisualizationofthecloudsinordertoshowthejet effectiveinsystemsinwhichthejetpowerisintherange
plasma flow within. P =1043 –1046ergs−1 andη >η . Furthermore,we
jet crit
Let M , m , c, σ , and σ be the black hole mass, find that, the maximum density weighted radial outflow
BH p T 100
the proton mass, the speed of light, the Thomson elec- velocityofclouds,(cid:104)v (cid:105),orequivalently,thecriticalEd-
r,w
tron scattering cross section, and the velocity dispersion dington ratio of the jets, η , depends weakly on filling
crit
in units of 100kms−1, respectively. We also define the factor, but strongly on the maximum size of clouds in
ratio of jet power to Eddington luminosity (the Edding- the galaxy bulge. The overall lower limit η (cid:38) 10−4
crit
Relativistic Jet Feedback 5
TABLE 1
Simulation parameters
New Simulation logPjet(a) nh(b) pISM/k(c) (cid:104)nw(cid:105)(d) fV(e) kmin(f) Rc,max(g) Mw,tot(h)
(erg) (cm−3) (cm−3K) (cm−3) (kpc−1) (pc) (109M(cid:12))
A........ 45 0.1 106 ··· ··· ··· ··· ···
B......... 46 1.0 107 1000 0.42 20 25.0 16
B(cid:48)........ 46 1.0 107 300 0.13 20 25.0 3.2
(cid:73) B(cid:48)(cid:48)....... 46 1.0 107 150 0.052 20 25.0 0.29
(cid:73) B(cid:48)(cid:48)(cid:48)....... 46 1.0 107 100 0.027 20 25.0 0.15
C......... 46 0.1 106 100 0.42 20 25.0 1.6
C(cid:48)........ 46 0.1 106 30 0.13 20 25.0 0.32
D........ 45 1.0 107 1000 0.42 20 25.0 16
(cid:73) D(cid:48)10...... 45 1.0 107 300 0.13 10 50.0 3.2
D(cid:48)........ 45 1.0 107 300 0.13 20 25.0 3.2
(cid:73) D(cid:48)(cid:48) ...... 45 1.0 107 150 0.052 10 50.0 0.29
10
(cid:73) D(cid:48)(cid:48)....... 45 1.0 107 150 0.052 20 25.0 0.29
(cid:73) D(cid:48)(cid:48) ...... 45 1.0 107 150 0.052 40 12.5 0.29
40
(cid:73) D(cid:48)(cid:48)(cid:48) ...... 45 1.0 107 100 0.027 10 50.0 0.15
10
(cid:73) D(cid:48)(cid:48)(cid:48)....... 45 1.0 107 100 0.027 20 25.0 0.15
E......... 45 0.1 106 100 0.42 20 25.0 1.6
E(cid:48)........ 45 0.1 106 30 0.13 20 25.0 0.32
F......... 44 0.1 106 100 0.42 20 25.0 1.6
F(cid:48)........ 44 0.1 106 30 0.13 20 25.0 0.32
G........ 44 1.0 107 1000 0.42 20 25.0 16
G(cid:48)........ 44 1.0 107 300 0.13 20 25.0 3.2
(cid:73) G(cid:48)(cid:48) ...... 44 1.0 107 150 0.052 10 50.0 0.29
10
(cid:73) G(cid:48)(cid:48)....... 44 1.0 107 150 0.052 20 25.0 0.29
(cid:73) G(cid:48)(cid:48) ...... 44 1.0 107 150 0.052 40 12.5 0.29
40
(cid:73) G(cid:48)(cid:48)(cid:48) ...... 44 1.0 107 100 0.027 20 25.0 0.15
H........ 43 0.1 106 100 0.42 20 25.0 1.6
(cid:73) H(cid:48)........ 43 0.1 106 30 0.13 20 25.0 1.6
(cid:73) I(cid:48)(cid:48)........ 43 1.0 107 150 0.052 20 25.0 0.29
(cid:73) I(cid:48)(cid:48)(cid:48)....... 43 1.0 107 100 0.027 20 25.0 0.15
Note. — Runs with run labels containing the same letter are runs with the same jet power, Pjet, and hot phase density, nh. Runs
labeledwithsingle, double, ortripleprimed(“(cid:48)”)lettersdenotelowerfillingfactorcounterpartstorunswithlessnumberofprimes. All
runs,otherthanthosewhoserunlabelcontainsthevalueofkmin inthesubscript,wereperformedwithkmin=20kpc−1.
(a) Jetpower.
(b) Densityofhotphase.
(c) p/k ofbothhotandwarmphases.
(d) Averagedensityofwarmphase.
(e) Volumefillingfactorofwarmphase.
(f) Minimumsamplingwavenumber.
(g) Maximumcloudsize.
(h) Totalmassinwarmphase.
for efficient feedback found by WB11 is only slightly of the warm phase. In Fig. 5 we show the evolution of
reduced for galaxies containing small cloud complexes the ratio of kinetic energy in clouds to injected jet en-
(Rc,max (cid:46)10pc, kmin =40kpc−1) but jets with Edding- ergy,Ekin,w/Pjett,asafunctionoftforall28simulations
ton ratios of η = 10−2 – 10−1 are required if cloud includingawarmphase. Inallcasesthatfractionishigh,
crit
complexes are large (R (cid:38)50pc, k =10kpc−1). reaching∼0.1–0.4,withdetailsdependingonjetpower
c,max min and ISM properties. The details of the dependence on
WB11 found that the jet-ISM interactions, despite
feedback efficiency on ISM properties are given in the
the porosity of clouds and the radiative losses of shock-
next two sections and the physics of how the high me-
impacted clouds, exhibit a high mechanical advantage,
chanical advantage is sustained and the energy transfer
meaning that substantial momentum transfer from the
occurs are investigated in §4.6.
jet to the clouds occurred through the energy injected
by the jet. We define the mechanical advantage in our
simulations at a given time as the ratio of the total ra- 4.2. Dependence on filling factor
dial outward momentum carried by clouds to the total Figure 6 shows the maximum values of (cid:104)v (cid:105) reached
r,w
momentum delivered by the jet up to that time. Fig- in the simulations as a function of f . The markers de-
V
ure4showsthecurvesforthemechanicaladvantageasa note simulations with equal values of P /n , as indi-
jet h
function of time for all 28 simulations including a warm cated in the legend. The lines of a given color connect
phase. For all simulations the mechanical advantage is simulations of equal power, also indicated by the label
much greater than unity. Most curves fall closely on top letter, and the line color indicates the hot phase density.
of each other along a narrow band up to at least 1Myr. Apart from the cases of different k in the D-series of
min
The high mechanical advantage generally leads to a runs, simulations grouped by connected lines, therefore,
high fraction of jet energy transferred to kinetic energy also indicate runs with equal values of P /n .
jet h
6 Wagner, Bicknell, & Umemura
4
3
2
1
0
−1
−2
−3
−4
4
3
)
2
3
−
1 m
c
0 (
n
−1
g
−2 o
l
−3
−4
4
3
2
1
0
−1
−2
−3
−4
Fig. 1.—Logarithmicdensitymaps(inunitsof cm−3)ofselectednewsimulations. Thedomainextentsineachpanelare1kpc×1kpc.
Theleftcolumnofpanelsshowafaceonviewofinitialthewarmgasdistribution. Thecenterandrightcolumnsofpanelsshowmidplane
slices at an advanced stage of the simulations for z =0 (reflected about x=0) and y =0, respectively. Top row: Run D(cid:48)(cid:48)(cid:48), a very low
filling factor run (fV =0.027); Middle row: Run D(cid:48)10, maximum cloud sizes of Rc,max =50pc; Bottom row: Run D(cid:48)4(cid:48)0, maximum cloud
sizesofRc,max=10pc.
In general, the dependence of (cid:104)v (cid:105) on filling factor of ablated material relative to the total mass of a cloud,
r,w
is weak and non-monotonic. In the B-series, D-series, because the ablation rate is proportional to the cloud
D -series, G-series, and I-series of the simulations (see surface area and the mass is proportional to the cloud
10
table 1 for nomenclature), we observe that for f (cid:38)0.1, volume, which decreases faster than the former for de-
V
lowerfillingfactorsdecreasethefeedbackefficiency,while creasing filling factor. When lowering the filling factor
for f (cid:46) 0.1, lower filling factors increase the feedback in the range f (cid:38) 0.1 the effect of reduced plasma con-
V V
efficiency. The reason for the weak dependence and the finementtimedominatesovertheeffectofincreasedfrac-
non-monotonicity is the competing effects of the cloud tional cloud ablation and results in lower mass-averaged
ablation rate and jet plasma confinement time. On the outflow velocities. In the range f (cid:46) 0.1, the increased
V
onehand,smallerfillingfactorsincreasethevolumeavail- cloud ablation rate dominates over the reduced plasma
able for the jet plasma to flood through, and thereby re- confinement time when reducing f , leading to higher
V
duce the confinement time, which reduces the impulse mass-averaged outflow velocities. Over the range of val-
delivered to the clouds over the confinement time. On ues for the filling factor studied here, these two effects
the other hand, smaller filling factors increase the mass counteract one another, and the dependence of (cid:104)v (cid:105) on
r,w
Relativistic Jet Feedback 7
48 kyr 67 kyr
89.2 pc 296 pc
96 kyr 125 kyr
400 pc 505 pc
158 kyr 187 kyr
624 pc 693 pc
Fig. 2.— Volume render of the density of the jet plasma and clouds for run D(cid:48). The jet plasma is textured in bluish green and the
cloudsinpurple. Theforwardshockoutliningthejet-blownenergybubbleisseenintranslucentgrey. Anovalexcavationismadeinthe
visualizationofthecloudsinordertoshowthejetplasmaflowwithin. Theviewmovesoutwardfromthecoreofthegalaxiesasthebubble
ofjetplasmaexpands,andthephysical(projected)sizeisindicatedbyscalebarsonthebottomrightineachpanel. Thesimulationdata
isreflectedaboutx=0andtheleftsideisrotatedby180◦ aboutthejetaxistoshowabackviewofthesimulation.
8 Wagner, Bicknell, & Umemura
(n,f)
h V
3.2 (1.0,0.42)
(0.1,0.42)
C
(1.0,0.13)
(0.1,0.13) D40′′ C′
3.0 (1.0,0.052)
B′′′
(1.0,0.027)
σ100=1.0η−1/4Pj1e/t,445 E BB′′
2.8 σ100=1.2η−1/5Pj1e/t,545 G40′′ ED′′′′ B′
) D
1s− 5.0 F D′′
km − F′ D′
/w(cid:3)2.6
g(vr,(cid:0) H G′′′G′′
o
l 4.0 H′ G′ D10′′′
2.4 − G
Fig. 5.— Fraction of jet energy going into kinetic energy of the
D10′′ warmphaseasafunctionoftimeforall28simulationscontaininga
I′′′ warmphase. Forallsimulations,0.4(cid:38)Ew,kin/Pjet(cid:38)0.1although
D10′ the maxima and the time taken to reach the maxima depend on
2.2 3.0 I′′ thejetpowerandISMparameters.
−
2.0 2−.0 1−.0 log(η)=0.0 1600 0.0270.052 0.13 0.42 l[oergg(Pcjmet/3nsh−)1]
1400 47
G10′′ C 46
43 44 45 46 45
log(Pjet/ergs−1) 1200 44
C′ 43
Fig. 3.—Maximummeanradialvelocityofclouds,(cid:104)vr,w(cid:105),against D40′′ nh[cm−3]
jet power for the simulations B – I(cid:48)(cid:48)(cid:48) of Table 1. The solid and 1.0
dtoasEhdeddinggretoynlilnuemsinaoresitlyo,cfioorfMco–nσstraenlattiηo,nsthweitrhatpioowoefrsje4tapnodwe5r, 1s]− 1000 B′′′B′′ B 0.1
respectively. The line colors indicate different hot-phase densities m
E
andthelinestylesrepresentdifferentfillingfactors,asindicatedin k 800
thelegend. [w(cid:3) G40′′ B′
vr,(cid:0)600 D′′′ D′′ E′ D
D′
F
400 F′
G′′′ G′′ G′ H
200 D10′′′D10′′ H′ G
I′′′ I′′ D10′
G10′′
0.0 0.1 0.2 0.3 0.4 0.5 0.6
f
V
Fig. 6.— Maximum mean radial velocity of clouds, (cid:104)vr,w(cid:105), ver-
sus cloud volume filling factor for the simulations B – I(cid:48)(cid:48)(cid:48) of Ta-
ble1. Thelinecolorsindicatedifferenthot-phasedensitiesandthe
marker styles group simulations with equal values of Pjet/nh, as
indicatedinthelegend.
f for constant P /n remains weak.
V jet h
The mechanical advantage (Fig.4) is slightly reduced
for systems with lower filling factor down to f =0.027,
V
but the dependence of the efficiency of transfer of jet
energy to kinetic energy of the warm phase (Fig.5)
Fig. 4.— Mechanical advantage versus time for all 28 runs in- on warm-phase filling factor parallels the weak (non-
cluding clouds. The mechanical advantage here is defined as the monotonic)dependenceofthemaximumoutflowvelocity
totaloutwardradialmomentumofcloudsattimetdividedbythe on filling factor.
total momentum delivered by the jet up to time t. The mechani-
Note that, by reducing the filling factor, we are also
caladvantageinallsimulations(cid:29)1indicatingstrongmomentum
couplingintheenergy-drivenregime. reducing the total mass of the warm phase. in contrast
tothis, wemaykeepthetotalmassandfillingfactorthe
same but change the maximum size of clouds by varying
Relativistic Jet Feedback 9
R [pc] between k =10kpc−1 and k =20kpc−1.
c min min
50 25 16.7 12.5 The reason for the strong cloud-size dependence and
linear scaling is that changing the cloud sizes at con-
1200 stant filling factor changes the rate of ablation relative
D40′′ to the total cloud mass without changing the jet plasma
confinement time. This is because only the amount of
surface area exposed to ablation relative to the volume
1000
of a cloud changes. Since k ∝ R−1, where R is the
min c c
Pjet=1045 ergs−1 cloud radius, the ratio of surface area to volume of a
cloud scales linearly with k . For higher k , the rate
min min
] 800 of ablation relative to the total mass of clouds increases,
1
−
s while the confinement time of the jet plasma does not
m change compared to runs with different k but iden-
min
[k 600 G40′′ tical fV. This allows for a far higher fraction of warm
w(cid:3) phasemasstobeacceleratedtohighervelocities,increas-
vr,(cid:0) D′′ ing the maximum value of (cid:104)vr,w(cid:105) reached in the run.
Anequivalentstatementtotheaboveexplanationuses
400 Pjet=1044 ergs−1 the concept of a jet-cloud “interaction depth”, τ , for a
jc
given distribution of clouds with varying k . In anal-
min
G′′ ogy to optical depth, the effective interaction depth may
be written τ = n πR2 R , where n is the num-
200 jc c c,max bulge c
D10′′ ber of clouds per unit volume (the number density of
clouds), and R is the radius of the region in the
bulge
G10′′ bulgewhichcontainsclouds. Thecloudsmaybethought
0 10 20 30 40 of as N scattering centers with a cross-section πRc2,max
kmin [kpc−1] randomly distributed in the volume (4π/3)Rb3ulge, and
the interaction depth may be thought of as a measure
Fig. 7.— Maximum mean radial velocity of clouds, (cid:104)vr,w(cid:105), of the average number of jet-cloud interactions any jet
versus minimum sampling wave number, kmin (and correspond- stream starting from the origin, including trajectories
ing maximum cloud size, Rc,max), for runs in the D-series (dia-
mondpoints,Pjet =1045ergs−1)andG-series(triangularpoints, along secondary streams, will experience. This formu-
Pjet = 1044ergs−1). The runs with filling factor fV = 0.052 are lationofaninteractiondepthisindeedrelevantbecause,
markedwithblackmarkers,labeled,andconnectedwithlinesand as we demonstrate in §4.6, the jet streams carrying en-
showthevariationof(cid:104)vr,w(cid:105)withkmin. Thegreymarkersclustering trainedhotandwarmphasematerialaredirectlyrespon-
aroundablackmarkerarerunsdifferingonlyinfillingfactor.
sible for the acceleration of clouds through their ram
pressure. The total number of clouds in the bulge is
k . The results for this are shown next.
min N = f R3 /R3 = n R3 . Therefore, the num-
V bulge c,max c bulge
4.3. Dependence on maximum cloud sizes ber density of clouds is nc = fV/Rc3,max, and the inter-
action depth is τ = πf (R /R ) = πf k .
Let us look at the D-series of runs, for which we have jc V bulge c,max V min
Hence, for a fixed f , τ ∝k .
varied the maximum size of clouds by varying k . The V jc min
min
The linear relation between the ratio of surface area
values of k are denoted by the subscript of the run
min
to volume of a cloud and k , or equivalently, the lin-
labels in Fig. 3. min
In Fig. 7, we plot the sequences D(cid:48)(cid:48) , D(cid:48)(cid:48), D(cid:48)(cid:48) , and ear relation between τjc and kmin, leads to a linear re-
G(cid:48)(cid:48) , G(cid:48)(cid:48), G(cid:48)(cid:48) against k with black10markers,40labels, lation between kmin and (cid:104)vr,w(cid:105), which is seen between
an1d0 connect4e0d by a linem. inThe grey, unlabelled markers k = 20kpc−1 and k = 40kpc−1. It is, how-
min min
are other runs with varying filling factors but otherwise ever, not clear whether one may extrapolate this rela-
identical parameters. The sequences in this figure, but tion to larger cloud complexes with sizes characteristic
also those in both Figs. 3 and 6, show that the depen- of giant molecular clouds (GMC), say of order several
dence of mean velocity on the maximum cloud size in a 100pc, given that the scaling between k = 10kpc−1
min
simulationisverystrong,andthatismuchstrongerthan and k = 20kpc−1 is steeper than linear. A possi-
min
the dependence on filling factor. ble reason for the steepening at larger R is that the
c,max
By halving the size-scale of clouds from kmin = larger inter-cloud voids cause a decollimation of the jet
20kpc−1 (D(cid:48)(cid:48)) to kmin = 40kpc−1 (D(cid:48)4(cid:48)0), the feedback streamsleadingtolessefficientmomentumtransfer. The
provided by the jet accelerates the clouds to a velocity a scaling may also be affected by resolution limitations to
factor of two greater, from 600kms−1 to 1200kms−1. capturingthefractalsurfaceofclouds, andbystatistical
Doubling the cloud sizes from k = 20kpc−1 to variations in the decreased number of jet-cloud interac-
min
kmin = 10kpc−1 decreases the maximum cloud veloc- tions for small kmin. It is difficult to predict the feed-
ities reached in the simulation by a factor of 3, from back efficiency with respect to GMCwith scales of order
200kms−1 to 600kms−1. η is therefore more sensi- 100pc from our simulations because these are generally
crit
tive to the maximum sizes of clouds than the volume not spherical and the effective interaction cross-section
filling factor of clouds. Moreover, the scaling between depends on orientation with respect to the jet streams.
k and (cid:104)v (cid:105) is nearly linear between k =20kpc−1 The simulations by Sutherland & Bicknell (2007) and
min r,w min
Gaibler et al. (2011) show that, if the molecular gas is
and k =40kpc−1, and somewhat steeper than linear
min
10 Wagner, Bicknell, & Umemura
distributed as a large coherent complex in a disc-like ge- slight deceleration can even be seen in some runs as the
ometry, the coupling between jet and ISM in terms of bubbleisincreasinglymass-loadedbywarmphasemate-
negative feedback through gas expulsion is weak. Ob- rial. In a few runs a return toward the theoretical line is
servations of some gas-rich radio galaxies indicate that visible, with gradients steeper than the theoretical limit
themoleculargasisnotcoupledasstronglyintooutflows for an energy-driven wind of given injection power. This
withthejetastheneutralorionizedmaterial(Ogleetal. happens because the medium inside the pressure bubble
2010; Guillard et al. 2012). while the jet plasma is confined by clouds is at a higher
Theexplanationsgivenherealsoapplytotheinfluence pressurethanthatofabubblethatisexpandinginaho-
of cloud sizes on the mechanical advantage and energy mogeneousatmospherewithambientdensityn . During
h
transfer efficiency from jet to warm phase in these sys- the breakout phase of the jet from the region filled with
tems. Both the mechanical advantage (Fig.4) and the clouds,thejetplasmaburstsoutoftheoutermostporous
energytransferefficiency(Fig.5)aresignificantlygreater channels and momentarily fills volumes at a faster rate
in systems with smaller cloud sizes. than a bubble that was not impeded and confined for
The sense in which the fractional cloud dispersal rate some duration by a porous, dense distribution of clouds.
depends on cloud sizes is the same as the dependence of The first panel (a) shows four runs of differing filling
theconditionsforstarformationinacloudonitssize,in factor from the D-series, for which k = 20kpc−1. A
min
that, the larger a cloud the more likely it is to collapse bubble evolving in a system with larger filling factor ex-
due to an external pressure trigger. Thus, whether jet pands more slowly. As we decrease the filling factor, the
mediated feedback induces or inhibits star-formation is deviation from an energy-driven bubble become smaller.
a sensitive function of the statistics of the warm phase We see the same behaviour for the runs in the G-series
distribution, in particular its size distribution. (panel d). The larger volume of channels available for
thejetplasmatofloodthrough,andtheresultingsmaller
4.4. The expansion rate of the quasi-spherical bubble confinementtime,isthedominantfactorthatdefinesthe
bubble expansion rate. The same trend is visible in the
In this section, we determine the departure of the out-
secondpanelforsimulationsofdifferingfillingfactor,for
flow energetics from that of an energy-driven bubble as
which k = 10kpc−1 (panel c), although the effect is
functions of warm phase parameters. Because our simu- min
muchweaker. Thisistheresultoftheconfinementtimes
lations include radiative cooling and a porous two-phase
andmassloadingfromhydrodynamicablationceasingto
ISM, we expect the energetics of the bubble that sweeps
vary much as the maximum cloud sizes increases.
up the ISM imparting momentum and energy to the
Theexpansionrateprofilesforrunswithdifferingmax-
clouds to lie between the energy-driven and momentum-
imum cloud sizes, but equal filling factor at f = 0.052
drivenlimits. Wediscussmomentum-drivenandenergy- V
is shown in panels (b) and (e) for the D-series and G-
driven outflows in relation to work in the literature sep-
series, respectively. The expansion rate of the bubble
arately in §4.5.
deviates increasingly from the theoretical rate with de-
Figure 8 contains six panels showing the evolution of
creasingmaximumcloudsizes. Thereasonforthisisthe
the bubble radius with time for different runs in the G-
increased mass ablation rate relative to the total mass
series, D-series, and B-series. We defined the bubble ra-
of clouds, as described in §4.3. Smaller maximum cloud
dius to be the radius of a hemisphere whose volume is
sizes for the same porosity lead to higher mass-loading
equivalenttothatsweptupbythepressurebubbleinthe
rates and decreased expansion rates of the bubble.
simulation. In each panel, the solid black line and the
While feedback efficiencies are mainly sensitive to the
dotted black line represent the theoretical, self-similar,
maximumcloudsizes,thedeviationofthebubbleexpan-
spherically symmetric evolution of the forward shock ra-
sion rate from a theoretical energy-driven rate is sensi-
diusandcontactdiscontinuity,respectively,ofanenergy-
tive to both the maximum cloud sizes and filling factor.
drivenbubble(wind)inauniformmediuminStage1,as
Within the parameter range studied here the deviation
defined by Weaver et al. (1977) (see also §6 Bicknell &
depends more strongly on filling factor than maximum
Begelman 1996). That stage represents an adiabatically
cloud size. For systems containing large clouds with
expanding bubble with constant injection power, which
smallfillingfactors,thebubbleevolutionapproachesthat
in our case is P . The solutions of the first stage are
jet of an ideal energy-driven bubble. Thus, it would seem
applicable here because radiative losses, although they
thatinthislimit,theseresultsencourageasub-gridAGN
improve the structural integrity of clouds and their sur-
feedback prescription in cosmological models, in which
vival time (Cooper et al. 2009), are energetically unim-
energy is injected isotropically into a small region, even
portant in our simulations. The forward shock and dis-
if the multiphase ISM conditions in the cores of gravita-
continuity evolve according to R =0.88(P /n )1/5t3/5
2 jet h tional potentials are not adequately resolved. However,
and R =0.86R , respectively. The location of the thin
C 2 thislimitisnotthesameasthatwhichleadstothemost
shell in the momentum-driven limit of a bubble expand-
efficient cases of negative AGN feedback. The latter is
ing in a uniform medium of mass density ρ is delin-
h attained in the limit of small filling factors and small
eatedbythedashedblackline,andgivenbytheequation
(cid:112) cloudsizes. Adistributionoflargerclouds, instead, may
R = 3/2(p˙/ρ )1/4t1/2, where p˙ is the momentum
shell h lead to positive feedback, e.g., pressure-triggered star-
injection rate (Dyson 1984).
formation.
Without focusing on a particular panel in Fig. 8, we
Cosmological SPH models commonly invoke negative
note that the bubble radius in some runs follows that
AGN feedback in a single-phase ISM, which essentially
of the theoretical prediction for an energy-driven bub-
corresponds to the hot phase in our simulations. The
ble closely, while in others the bubble radius initially in-
heating rate of the hot phase will therefore likely al-
creasesmoreslowlythantheratepredictedbytheory. A