The many lives of active galactic nuclei-II: The formation and evolution of radio jets and their impact on galaxy evolution

TL;DR

This paper presents a new semi-analytic model integrating radio jet physics into galaxy evolution simulations, enabling detailed predictions of radio AGN properties and their effects on galaxy development.

Contribution

The authors develop and calibrate a novel model that tracks radio jet evolution and its impact within a cosmological galaxy formation framework, improving predictive capabilities.

Findings

Model successfully reproduces the stellar mass--radio luminosity relation at z~0.

Accurately predicts the radio luminosity function up to z~1.

Enables creation of customized AGN mock survey catalogs.

Abstract

We describe new efforts to model radio active galactic nuclei (AGN) in a cosmological context using the SAGE semi-analytic galaxy model. Our new method tracks the physical properties of radio jets in massive galaxies, including the evolution of radio lobes and their impact on the surrounding gas. This model also self consistently follows the gas cooling-heating cycle that significantly shapes star formation and the life and death of many galaxy types. Adding jet physics to SAGE adds new physical properties to the model output, which in turn allows us to make more detailed predictions for the radio AGN population. After calibrating the model to a set of core observations we analyse predictions for jet power, radio cocoon size, radio luminosity, and stellar mass. We find that the model is able to match the stellar mass--radio luminosity relation at $z\sim0$, and the radio luminosity function out to $z\sim1$. This updated model will make possible the construction of customised AGN-focused mock survey catalogues to be used for large-scale observing programs.

Tables

Table 1. Table 1: The new AGN model parameters, their best values, plausible ranges, equation number(s), section(s), and a brief description, as used throughout this paper.

Parameter	Best value	Plausible range	Equation(s)	Section(s)	Description
				Environment
${\beta }_{\mathrm{eff}}$	0.9 $b$	$\le$0.66	8	4.1	Effective $\beta$ in Makino et al. (1998)
$\beta$	3(${\beta }_{\mathrm{eff}}$)	$\le$ 2	13, 4.2.1, 15,16, 33,34	4.2.1,4.2.2, 4.2.3,5.2	Slope in the $\beta$ density profile
$r_0$	0.22 $r_s$	fixed	8,13,17,33, 35	4.1,4.2.1,4.2.3,5.2	Core radius in the Makino et al. (1998) density profile
				Lobe dynamics
${\Gamma }_x$	5/3	fixed	4.2.1	4.2.1	ICM adiabatic index
${\Gamma }_l$	4/3	fixed	4.2.1	4.2.1	Radio lobe adiabatic index
$\eta$	0.35	0.001-1	11	4.2	Jet generation efficiency
${\kappa }_R$	1.0	0.0-1.0	10	4.2	Radio mode feedback efficiency (Croton et al., 2006)
$a_D$	1.0	0.8–1.8	13,4.2.1,34	4.2.1,5.2	Dimensionless constant related to radio source size (Alexander, 2002)
$R_T$	6	1.3–6	4.2.1,34,35	4.2.1,5.2	Cocoon axial ratio
				Radio emission
$\nu$	1.4  GHz	–	35	5.2	Rest-frame observing frequency
$p$	2.6	2.1-2.7	31,32,5.2, 35	5.2	Power-law electron injection index
${\gamma }_1$	1	1–${10}^3$	30	5.2	Minimum Lorentz factor
${\gamma }_2$	${10}^6$	${10}^5$–${10}^6$	30	5.2	Maximum Lorentz factor
$C_2\left(p\right)$	$2.04\times {10}^{-3}$	fixed	5.2	5.2	Constant derived from $p$
$A1$	$2.9\times {10}^{-9}$	fixed	32	5.2	Constant derived from $p$

Equations

$$T_{\mathrm{vir}}=\frac{1}{2}\frac{\mu m_p}{k_B}{V}_{\mathrm{circ}}^2,$$

$$\rho (r)=\frac{m_{\mathrm{hot}}}{4\pi R_{\mathrm{vir}}{r}^2}.$$

$$\rho (r)=\frac{{\rho }_0}{[1+(r/r_c{)}^2{]}^{3\beta /2}},$$

$$\rho (r)={\rho }_0{e}^{-27b/2}(1+r/r_s{)}^{27b/(2r/r_s)}.$$

$$b=\frac{2C}{9\gamma }{\left[\ln \left(1+C\right)-\frac{C}{1+C}\right]}^{-1},$$

$$C=\frac{4}{1+z}{\left(\frac{M_{\mathrm{vir}}}{1.4\times 10^{14}{\mathrm{M}}_{\odot }}\right)}^{-0.13},$$

$${\rho }_0=\frac{m_{\mathrm{hot}}}{4\pi r_s^3}{e}^{27b/2}{\left[{\int }_0^C{x}^2(1+x{)}^{27b/2x}dx\right]}^{-1},$$

$$\rho (r)=\frac{{\rho }_{0}A(b)}{[1+(r/r_{\mathrm{c},\mathrm{eff}}{)}^2{]}^{3{\beta }_{\mathrm{eff}}/2}}.$$

$${\dot{m}}_{\mathrm{BH}}=\frac{2.5\pi G^2{m}_{\mathrm{BH}}^2{\rho }_0}{c_s^3},$$

$${\dot{M}}_{\mathrm{BH}}={\kappa }_R\frac{15}{16}\pi G\mu m_p\frac{kT}{\lambda }{M}_{\mathrm{BH}}.$$

$$Q_{\mathrm{jet},\mathrm{hot}}=\eta {\dot{M}}_{\mathrm{BH}}{c}^2,$$

$$Q_{\mathrm{jet},\mathrm{cold}}=6.3\times 10^{-5}{\left(\frac{M_{\mathrm{BH}}}{10^9{\mathrm{M}}_{\odot }\text{/}}\right)}^{-0.1}{\dot{m}}_{\mathrm{BH}}^{0.2}({\dot{M}}_{\mathrm{BH}}{c}^2)[W],$$

$$r_{\mathrm{cocoon}}(t)=a_D{r}_0{\left(\frac{t}{\tau }\right)}^{3/(5-\beta )}.$$

$$a_D={\left[\frac{R_T^4({\Gamma }_{\mathrm{x}}+1)({\Gamma }_{\mathrm{l}}-1)(5-\beta {)}^3}{18\pi \left(9\left[{\Gamma }_{\mathrm{l}}+({\Gamma }_{\mathrm{l}}-1)R_{\mathrm{T}}^2/2\right]-4-\beta \right)}\right]}^{1/(5-\beta )}.$$

$$r_{\mathrm{shock}}=\frac{1}{\lambda }{r}_{\mathrm{cocoon}},$$

$$1-{\lambda }^3=\frac{15}{4(11-\beta )}.$$

$$T_{\mathrm{shock}}=\frac{15}{16}\frac{3-\beta }{11-\beta }\left(\frac{\mu m_H}{k_B}\right){\dot{r}}_{\mathrm{shock}}^2,$$

$$m_{\mathrm{shock}}(r_{\mathrm{shock}})=16\pi {\rho }_{g0}{r}_0^3\left[\ln \left(1+\frac{r_{\mathrm{shock}}}{r_0}\right)-\frac{r_{\mathrm{shock}}}{r_{\mathrm{shock}}+r_0}\right].$$

$$T_{\mathrm{new}-\mathrm{hot}}=T_{\mathrm{hot}}+\frac{m_{\mathrm{shock}}}{m_{\mathrm{hot}}}(T_{\mathrm{shock}}-T_{\mathrm{hot}}).$$

$$t_{\mathrm{cool}}(r)=\frac{2}{3}\frac{\mu m_p{k}_B{T}_{\mathrm{new}-\mathrm{hot}}}{\rho (r)\Lambda (T_{\mathrm{new}-\mathrm{hot}},Z)},$$

$${\dot{m}}_{\mathrm{cool}}=4\pi {\rho }_g(r_{\mathrm{cool}})r_{\mathrm{cool}}^2{\dot{r}}_{\mathrm{cool}},$$

$${\dot{m}}_{\mathrm{cool}}=4\pi {\rho }_0\left(\frac{r_{\mathrm{cool}}^3}{t_{cool}}\right)\left(\frac{r_{\mathrm{cool}}}{r_s}\right)\left(\frac{e^{13.5b}}{13.5b}\right)$$

$$\times \left(\frac{{\left(1+\frac{r_{\mathrm{cool}}}{r_s}\right)}^{13.5b\left(\frac{r_s}{r_{\mathrm{cool}}}\right)}}{\ln \left(1+\frac{r_{\mathrm{cool}}}{r_s}\right)-\frac{r_{\mathrm{cool}}}{r_s+r_{\mathrm{cool}}}}\right).$$

$${\dot{m}}_{\mathrm{cool}}=\frac{m_{\mathrm{hot}}}{R_{\mathrm{vir}}/V_{\mathrm{vir}}}.$$

$$t_{\mathrm{on}}=120{\left[\frac{m_{\ast }}{10^{11}{\mathrm{M}}_{\odot }}\right]}^{0.7}\mathrm{Myr}.$$

$$t_{\mathrm{return}}=2\frac{r_{\mathrm{shock}}}{c_s}\approx 2\frac{r_{\mathrm{shock}}}{V_{\mathrm{vir}}}.$$

$$\delta =0.05{\left[\frac{m_{\ast }}{10^{11}{\mathrm{M}}_{\odot }}\right]}^{1.5},$$

$$t_{\mathrm{off}}(m_{\ast },t_{\mathrm{on}})=t_{\mathrm{on}}\left[\frac{1}{\delta }-1\right].$$

$$t_{\mathrm{quiet}}=t_{\mathrm{off}}-t_{\mathrm{return}}.$$

$$f_{\mathrm{cool}}=\frac{t_{\mathrm{quiet}}}{t_{\mathrm{on}}+t_{\mathrm{off}}}+\left[\frac{t_{\mathrm{on}}+t_{\mathrm{return}}}{t_{\mathrm{on}}+t_{\mathrm{off}}}\right]\left(\frac{m(r_{\mathrm{cool}})-m(r_{\mathrm{shock}})}{m(r_{\mathrm{cool}})}\right).$$

$${\dot{m^{\prime }}}_{\mathrm{cool}}=f_{\mathrm{cool}}{\dot{m}}_{\mathrm{cool}}.$$

Full text

††pagerange: The many lives of active galactic nuclei-II: The formation and evolution of radio jets and their impact on galaxy evolution–References††pubyear: 2012

1 Introduction

The paradigm of hierarchical galaxy formation provides a remarkably successful framework to explain many properties of the observed galaxy population, particularly at intermediate masses (Press & Schechter, 1974; White & Rees, 1978; White & Frenk, 1991; Lacey & Cole, 1994). With early models over-predicting galaxy counts at both the bright and faint end of the luminosity function, it was quickly realised how important feedback processes were for reproducing the observations. At the low mass end, heating due to the reionization of the Universe at early times, and supernovae ejecta from the most massive stars, became critical to explain the deficit of faint galaxies with respect to the number of low mass halos (Efstathiou, 1992). However in the most massive haloes such feedback modes proved insufficient to offset rapid cooling of the hot gas, leading to a subsequent overproduction of massive elliptical galaxies by many orders-of-magnitude.111The ratio of galaxy stellar mass to dark matter halo mass decreases with increasing mass. Hence, even at a constant specific star formation rate the efficiency with which supernova can eject gas out of the gravitational potential decreases with mass. In reality, massive ellipticals have lower specific star formation rates than spirals, and the impact of the supernova is even smaller. Furthermore, optical observations showed these massive ellipticals to have formed only a small fraction of their stars over the last half of the Hubble time (Bender & Saglia, 1999), in sharp contrast to the models which predicted that these objects should be rapidly forming stars at the present epoch.

A closely related problem was presented by galaxy clusters. Observations of X-ray emission suggested that there should be runaway cooling of the hot ($\sim 10^{7}$K) gas (the so-called “cooling catastrophe”; Cowie & Binney, 1977; Fabian & Nulsen, 1977), however no such rapidly cooling gas at temperatures below about a third of the virial temperature was found (Tamura et al., 2001; Peterson et al., 2003). A new mechanism was required to offset gas cooling on scales of tens of kpc, and suppress excessive galaxy growth from star formation.

Silk & Rees (1998) pointed out that the presence of a supermassive black hole (SMBH) at the centre of every massive galaxy provides such a mechanism, through conversion of a fraction of accreted mass to thermal or kinetic energy (Penrose & Floyd, 1971; Blandford & Znajek, 1977; Blandford & Payne, 1982); this is the mechanism invoked to explain the presence of Active Galactic Nuclei (AGN). By coupling a fraction of the AGN energy output to the rapidly cooling hot halo gas, a new generation of galaxy formation models (Granato et al., 2004; Monaco et al., 2006; Croton et al., 2006; Bower et al., 2006) succeeded in reproducing the quenching of late epoch star formation in the most massive galaxies. While successful in the sense of reproducing the observed galaxy properties, these models do not make detailed predictions for the mechanisms through which AGN energy couples to the gas.

Broadly speaking, AGN feedback can be either mechanical (through radio jets) or radiative (see reviews by Cattaneo et al., 2009; Somerville & Davé, 2015), with radiative feedback affecting the ISM of the host galaxy, and radio jets doing feedback on the larger halo scales. For this reason, Croton et al. (2006) referred to the mode of AGN feedback responsible for suppressing star formation since $z\sim 1$ as “radio” or “maintenance” mode feedback. Dramatic evidence for this mode of feedback in action was found in the Perseus (Böhringer et al., 1991; Fabian et al., 2003) and Virgo (Churazov et al., 2001; Forman et al., 2005) clusters, where lobes of radio emitting plasma were observed to displace the X-ray emitting gas away from the rapidly cooling central regions.

Observationally, the radio loud AGN fraction is found to be a strong function of host galaxy properties: Sadler et al. (1989) found more than 40% of massive ellipticals to host low-luminosity radio sources. Best et al. (2005) showed that this fraction scales strongly with stellar mass, and this scaling is consistent with a picture of “maintenance mode” feedback, in which massive galaxies with rapidly cooling hot haloes are more likely to host AGN, which can in turn provide the feedback needed to offset the otherwise imminent cooling catastrophe. In clusters, the state of the cooling gas appears to be intimately connected to the probability of finding an AGN in the cooling gas: Burns (1990) reported that over 70% of cool core clusters hosted radio sources, compared to only 23% of non-cool core clusters. Other authors (e.g. Mittal et al., 2011) have found similar results.

While many of these sources are compact (Shabala et al., 2008; Sadler et al., 2014), when the jets are powerful enough they can break through the galaxy disk and inflate pairs of lobes of radio emitting synchrotron plasma, with the largest lobes up to a Mpc in size (Laing et al., 1983; Saripalli et al., 1986, see Banfield et al. (2016) for an example of a recent discovery). Morphologically, extended radio lobes are observed to be either core or edge-brightened; such sources are classified as Fanaroff-Riley type I and type II, respectively (Fanaroff & Riley al., 1974).

Large radio sources can do feedback well outside the host galaxy disk via two channels. First, powerful radio sources can drive strong shocks through the surrounding gas (Schoenmakers et al., 2000; Rawlings & Jarvis, 2004; Shabala et al., 2011), heating and uplifting it to large radii. Second, lower power sources in galaxy clusters are often observed to give rise to buoyant bubbles of radio plasma, which displace the hot X-ray emitting gas as they rise through the cluster (Böhringer et al., 1991; Churazov et al., 2001; Fabian et al., 2003; Forman et al., 2005).

Multiple lines of evidence suggest that radio AGN activity must be episodic. Theoretically, AGN intermittency is expected in a self-regulating feedback process (e.g. Bahcall et al., 1997; Kawata & Gibson, 2005; Novak et al., 2011; Gaspari et al., 2015). Observationally, the presence of multiple shocks and ripples around the X-ray cavities in the Perseus (Fabian et al., 2003) and Virgo (Forman et al., 2005) clusters have been interpreted as remnants of multiple AGN outbursts. Perhaps the most dramatic evidence for recurrent radio AGN activity is provided by sources with multiple pairs of radio lobes (e.g.; Schoenmakers et al., 2000). In these “double-double” radio sources, the inner pair of lobes is interpreted to correspond to the most recent AGN outburst, while the outer pair of lobes is due to a previous episode of AGN activity. In a number of sources (e.g. Centaurus A; Israel, 1998; Feain et al., 2009, a similar suggestion has also been been made for Cygnus A, Chon et al. 2012), the jet directions are often misaligned between outbursts, providing a mechanism for isotropizing the coupling of jet energy to cluster gas.

The interaction between radio jets / lobes and the hot atmospheres into which they expand sets both the properties of the observed AGN, such as their sizes, morphologies, luminosities and radio spectra; and, through feedback, those of their host galaxies. The effect of environment on both AGN properties and the amount of feedback they do has been extensively studied through radio source dynamical models (Scheuer, 1974; Begelman & Cioffi, 1989; Kaiser & Alexander, 1997; Alexander, 2002; Shabala et al., 2008; Turner & Shabala, 2015), and numerical simulations (Norman et al., 1982; Reynolds et al., 2002; Basson & Alexander, 2003; Krause, 2005; Mendygral et al., 2012; Hardcastle & Krause, 2013, 2014). Recently, Godfrey & Shabala (2016) pointed out that the often quoted relationship between jet kinetic power and radio luminosity (derived from observations of X-ray deficient cavities, e.g.; Birzan et al., 2004, 2008; Cavagnolo et al., 2010) is in fact driven by strong selection effects, and the true relationship is more complicated due to environmental effects (e.g. Barthel & Arnaud, 1996; Kaiser & Alexander, 1997; Hardcastle & Krause, 2013, 2014). Observationally, Khosroshahi et al. (2017) found that quantities such as the radio luminosity of the brightest group galaxies strongly depend on their environment, such that the brightest group galaxies in dynamically young (evolving) groups are an order-of-magnitude more luminous in the radio than those with a similar stellar mass but residing in dynamically old (evolved) groups (for a definition of old and young groups see Raouf et al., 2014). This finding is consistent with results of hydrodynamical simulations (Raouf et al., 2016), which suggest that the IGM in dynamically evolved groups is hotter for a given halo mass than that in evolving groups.

Despite the importance of feedback from radio jets, to date no galaxy formation model has attempted to predict simultaneously the properties of both galaxies and AGN in detail. Fanidakis et al. (2012) produced an AGN radio luminosity function, however they employed both an arbitrary scaling between jet power and radio luminosity, and an arbitrary normalisation of the radio luminosity function; they therefore did not properly account for the intermittency of the feedback process. Shabala & Alexander (2009) used a dynamical radio source model to quantify the feedback from jets on the hot gas, however they made no predictions for the resulting AGN properties.

A major difficulty in connecting jet models and simulations with observations lies in quantifying the environments into which the jet expands. Density and temperature profiles derived from X-ray observations are the gold standard, however these are typically biased towards low-redshift dense environments. Turner & Shabala (2015) introduced a model which connects radio AGN dynamical models with semi-analytic models of galaxy formation. In this approach, the semi-analytic model was used to estimate the total hot gas mass contained with the virial radius of a given dark matter halo; this gas was then assumed to follow a density profile consistent with observations of local clusters (e.g.; Vikhlinin et al., 2006). Turner & Shabala (2015) used this approach to derive the physical properties (jet powers and ages) of radio AGN in a volume-limited low-redshift sample and found it reproduces well a number of key AGN observables. Here, we adopt their techniques to develop a new prescription for intermittent AGN feedback in the SAGE (Semi-Analytic Galaxy Evolution) galaxy formation model (Croton et al., 2016), updating the more simplistic “radio mode” model introduced in Croton et al. (2006). Our code is publicly available as a fork of the original SAGE repository.222https://github.com/mojtabaraouf/sage

This paper is organised as follows: In Sections 2 and 3 we describe our N-body and semi-analytic framework, respectively. In Section 4 we describe the major features of our model for AGN feedback. Model constraints and predictions are discussed in Section 5. We present the summary of our results in Section 6.

2 The dark matter simulation: Millennium

In this work we use the Millennium Simulation (Springel et al., 2005) N-body dark matter halo merger trees as input into the “Semi-Analytic Galaxy Evolution” (SAGE) galaxy formation model (Croton et al., 2016), which we update with new AGN physics (described below). This simulation is important as it provides the structural backbone onto which galaxies (and hence black holes) can be evolved.

The Millennium Simulation was run using the popular GADGET-2 code and adopted a cosmological model consistent with the first year Wilkinson Microwave Anisotropy Probe data (Spergel et al., 2003, WMAP-1, with parameters $\Omega_{m}=0.25$, $\Omega_{\Lambda}=0.75$ and ${\rm H_{0}}=100{\rm h}$ km $s^{-1}$ $\rm Mpc^{-1}$ where $\rm h=0.73$). The simulation box of $(500\rm h^{-1}\rm Mpc)^{3}$ contained $2160^{3}$ particles and had a mass resolution of $8.6\times 10^{8}\rm h^{-1}M_{\odot}$ per particle. 64 snapshots of the particle evolution were written to disk, spaced approximately logarithmically in scale factor between ${z=127}$ and ${z=0}$. Dark matter halos were then found amongst the mass distribution in each output, characterised as associations of 20 or more bound particles. For this, a combination of the Friends-of-Friends (Davis, 1985) and SUBFIND (Springel et al., 2001) halo finding algorithms were applied. After all dark matter halos had been identified they were then linked into their respective merger trees using the L-HALOTREE code. From this the full growth history of each ${z=0}$ object in the box could be inferred. For more information on the Millennium Simulation see Springel et al. (2005).

Note that the exact simulation and cosmology used is somewhat secondary for the goals of this paper. SAGE, and our new AGN model coupled to it, can be run on any simulation, and the model parameters allow it to be calibrated to match key observations in physically sensible ways.

3 The galaxy formation model: SAGE

We only give a brief introduction to the base SAGE galaxy formation model here and refer the interested reader to Croton et al. (2016) for a full description. Beyond this, the rest of the paper will focus on the model changes, primarily related to supermassive black hole growth and outflows, that lead to a more physically motivated coupling between AGN feedback and galaxy evolution. This includes the key observables those working in this area of research might want such a model to predict.

SAGE is an updated version of the semi-analytic model first introduced in Croton et al. (2006). It analytically follows the movement of baryons through different mass reservoirs, computed on top of the numerically determined evolving dark matter halo mass distribution, characterised by a set of N-body simulation halo merger trees. Baryonic reservoirs include the hot halo gas, cold disk gas, stars in the disk, black holes, and gas ejected from the halo due to feedback events. How this mass moves between the reservoirs is determined by a series of coupled differential equations that describe each physical process believed important. These include hot gas cooling into the disk (hot$\rightarrow$cold), star formation (cold$\rightarrow$ stars), gas heating from supernova or AGN feedback (cold$\rightarrow$ hot), ejection (cold/hot$\rightarrow$ ejected) and later reincorporation (ejected$\rightarrow$ hot), the effects of reionization, and so on. Most processes have one or more efficiency parameters that allow us to individually control their relative importance, although in reality the different components of galaxy and AGN evolution are highly intertwined and not so easily broken up (see e.g. Mutch et al., 2013).

To ensure SAGE produces a galaxy population akin to that observed around us it is calibrated by hand to statistically match a set of key observables. Our primary observable is the local stellar mass function, with a set of secondary observables being the star formation rate density history, net cooling rate–temperature relation and the AGN radio luminosity function (part of the new model extension), shown and discussed below in Section 5.1. All model results assume a Universe where ${h=0.73}$, and when relevant, a Chabrier initial mass function (Chabrier, 2003) to compare model to observed stellar masses.

4 A new model of AGN feedback in SAGE

The impact of an AGN jet on the surrounding gas outside of a galaxy was captured by the model of Shabala & Alexander (2009), who calculated the dynamics of jet-inflated radio lobes that propagate supersonically and shock heat the surrounding gas. Their work is applicable both for “hot mode” and “cold mode” black hole accretion that produces AGN jets with different efficiencies, and hence lends itself nicely for incorporation into SAGE. In the hot mode, also called the radio mode (Croton et al., 2006), accretion occurs at low Eddington rates and the accretion disk is geometrically thick and optically thin. This results in longer cooling times and more powerful jets. In contrast, the cold mode is associated with high Eddington rates and is described via the standard Shakura-Sunyaev geometrically thin, optically thick disk. The cold mode is radiatively efficient and produces a quasi-black-body spectrum, including optical / UV continuum and narrow-line AGN emission.

Both radio lobe dynamics and associated feedback are sensitive to the environment into which the lobes expand (Kaiser & Alexander, 1997; Willott et al., 1999; Turner & Shabala, 2015). For self-consistently between the new AGN and existing galaxy model, we therefore first need to update the SAGE hot halo gas density profile into which the jets propagate. This in turn has an effect on the hot gas cooling rates into the galaxy. We then analytically describe the AGN jet evolution itself and how this replaces the existing more simplistic SAGE radio mode feedback.

We note that all figures and results presented in this paper arise from this new model using the parameters given in Table 1. These parameters are defined in the various subsections below when they are introduced. Control over the model behaviours  is achieved through  parameterizing its different components, essentially to  set  their  relative efficiency, as (currently)  no model of galaxies  can be described from first principles alone. Such “tuning” also accounts for the many finer details that are ill-understood and can not be easily modeled explicitly. Parameter values must be consistent with the physical process they are meant to describe and are set so that the overall model agrees with a key set of galaxy and AGN observations. Beyond a certain point, adding more parameters typically makes the model harder to calibrate.

4.1 The hot gas density profile and cooling

In SAGE every dark matter halo is assumed to carry its cosmic share of baryons, taken to be $f_{b}=0.17$, consistent with the WMAP1 results of Spergel et al. (2003). These baryons begin their life as diffuse hot gas with primordial composition around the galaxy. However, as mentioned above, with time the gas transforms under the action of the many and varied physical processes of galaxy evolution to populate gas reservoirs of several different phases, as well as stars, black holes, and the heavy elements (Baugh et al., 2007).

Considering first this hot gas in virial equilibrium, its temperature can be described by (Sutherland & Dopita, 1993)

[TABLE]

where $m_{p}$ is the mass of proton, $\mu$ is the mean molecular weight of gas, $k_{B}$ is the Boltzmann constant, and $V_{\rm circ}$ is the circular velocity at the virial radius of halo. In what follows we approximate $V_{\rm circ}$ by $V_{\rm vir}$ of the halo. Of course, after an AGN jet ploughs through this gas it is unlikely to be in equilibrium, at least for a time. The new hot temperature of such gas is discussed in the following sections by taking into account the size and properties of an expanding shocked gas bubble due to an AGN outflow.

The hot gas density profile is similarly important as it directly determines the cooling (i.e. feeding) rate of gas into the galactic disk, gas which ultimately leads to the formation of new stars. The original SAGE model, like many other semi-analytic models before it, assumed that this hot gas can be represented as a simple isothermal sphere with temperature given above, and having a density profile (White & Frenk, 1991)

[TABLE]

Here $m_{\rm hot}$ is the total hot gas mass in the halo within the virial radius, $R_{\rm vir}$.

However this approximation is clearly an oversimplification, especially when compared to X-ray observations of nearby cluster systems which show that the density profile is better described by an isothermal $\beta$-model (Cavaliere & Fusco-Femiano, 1978; Fukazawa et al., 2004; Pointecouteau et al., 2004; Vikhlinin et al., 2006),

[TABLE]

where $r_{c}$ defines a characteristic core radius, and $\beta$ describes the changing inner and outer slopes of the profile.

For the purpose of the present study we adopt the density profile described by Makino et al. (1998). The Makino profile closely matches the observationally fit $\beta$-model given by Equation 3, albeit based on a better theoretical foundation. It characterises the mass distribution that arises from isothermal gas in hydrostatic equilibrium analytically embedded in a universal “NFW” (Navarro, Frenk, & White, 1996) dark matter halo. Again we assume that the gas is at the virial temperature of the halo. While a more realistic model would include a clustercentric radius-dependent temperature (e.g. Arnaud et al., 2010), we adopt a single-temperature model to minimize the number of free parameters. Using this model, Makino et al. predict the core density, core radius, $\beta$-parameter, and X-ray luminosity of clusters as function of halo mass and temperature. The Makino profile is given by

[TABLE]

Here $\rho_{\rm 0}$ is the core density of hot gas, while $r_{s}$ is the usual NFW scale radius. If one assumes the hot gas is at the virial temperature of the halo (Equation 1), the $b$ parameter can be written as

[TABLE]

with $\gamma$ being a parameter of order unity that determines the efficiency of shock heating; we adopt $\gamma=1.5$. The concentration, $C=R_{\rm vir}/r_{s}$, was empirically determined by Bullock et al. (2001) to be

[TABLE]

with $z$ the halo’s redshift. Finally, the central density of the hot gas, $\rho_{\rm 0}$, can be determined by integrating the Makino density profile out to the halo virial radius and requiring the total mass equal $m_{\rm hot}$:

[TABLE]

where we have used the change of variable $x\equiv r/r_{s}$.

We can combine the above Makino profile representations for $b$ (Equation 5), $r_{s}=R_{\rm vir}/C$ (using Equation 6), and $\rho_{0}$ (Equation 7) to obtain a $\beta$-model-type expression,

[TABLE]

The parameters completing the profile are: $A(b)=-0.178b+0.982$, $r_{0}\equiv r_{\rm c,eff}=0.22r_{s}$ and $\beta_{\rm eff}=0.9b$, all valid over scales $0.01r_{s}<r<10r_{s}$.

In Figure 1 we compare the density profile of Makino et al. (1998) (Equation 4, red solid line) to that of the $\beta$-model (Equation 3, blue long-dashed line) as function of radius (in units of $R_{\rm vir}$) for the same total hot gas mass ($m_{\rm hot}=10^{13}\rm h^{-1}M_{\odot}$). Also over-plotted is the original isothermal profile ($\propto 1/r^{2}$) used in SAGE (Equation 2, gray short-dashed line). The Makino and $\beta$-model profiles agree comfortably well across the entire range plotted, while the isothermal profile is seen to deviate significantly from the other two on scales less than $\sim 0.1R_{\rm vir}$, where cooling is most significant.

4.2 Jet generation and propagation

Following the original model described in Croton et al. (2006), the accretion rate of gas feeding the black hole is approximated by the Bondi-Hoyle formula (Bondi, 1952),

[TABLE]

where $m_{\mathrm{BH}}$ is the black hole mass and $\rho_{0}$ is the density of accreting hot gas around the black hole. $c_{s}\equiv V_{\rm vir}$ and $G$ are the speed of sound in the gas and the gravitational constant, respectively.

The key unknown above is the central gas density, which can be approximated by assuming a maximal cooling flow in this region, as explained in section 9.1 of Croton et al. (2016). This leads to an expression for $\rho_{0}$ that can be inserted back into Equation 9, after which the accretion rate simplifies to

[TABLE]

Equation 10 is now written in terms of quantities that SAGE naturally produces. Note that the efficiency parameter, $\kappa_{R}$, was introduced by Croton et al. (2016) to modulate the effectiveness of the feedback within the cooling–heating cycle.

Moving beyond this original accretion implementation, we look to build a physical model of the outflows produced as a result of gas accretion onto a black hole. There are two different accretion states that can both happen when removing angular momentum from the accretion disk by viscosity (Narayan et al., 1998; Meier, 2001; Körding et al., 2006; Fender et al., 2004). The “hot mode” occurs for low accretion rates with respect to the Eddington rate ($\dot{M}_{\rm BH}<\alpha_{\rm crit}\dot{M}_{\rm edd}$). Here, the flow is geometrically thick and optically thin, which results in longer cooling times and more powerful jets (for a given black hole mass and accretion rate). The thermal energy of the inflowing gas is advected inward, known as an advection dominated accretion flow (ADAF; Narayan & Yi, 1995).

We relate “hot mode” jet power to the accretion rate via the following relation:

[TABLE]

where $c$ is the speed of light and $\eta$ is the jet efficiency, estimated using relativistic magnetohydrodynamical simulations of jet generation to have typical values between 0.002 and 1; the exact value is largely set by black hole spin (Benson & Babul, 2009). In Section 5.3 we use observations to constrain the best-fitting value of $\eta$ to 0.35.

In contrast, the “cold mode” of AGN accretion occurs when inflow rates are high compared with the Eddington rate ($\dot{M}_{\rm BH}>\alpha_{\rm crit}\dot{M}_{\rm edd}$). Here, the accretion flow can be described by the standard thin disk solution of Shakura & Sunyaev (1973), where the disk is geometrically thin and optically thick. This produces outflows that are radiatively efficient and hence the AGN has weaker jets (for a given black hole mass and accretion rate) in comparison to the ADAF case (Narayan et al., 2002). The cold mode jet power is given by (Shakura & Sunyaev, 1973; Meier, 2001, Equation 5):

[TABLE]

where $\dot{m}_{\rm BH}\equiv\frac{\dot{M}_{\rm BH}}{\dot{M}_{\rm edd}}$.

Note that the parameter $\alpha_{\rm crit}$ defines the transition from hot to cold accretion in units of the Eddington rate and is set at 0.03 throughout this work (see also Shabala & Alexander, 2009; Merloni & Heinz, 2008). In our model we generate both ADAF and thin disk jets, with different jet generation efficiencies.

4.2.1 Radio source expansion

The accretion disk outflows just described drive jets into the surrounding gas, creating expanding cavities and producing synchrotron emission seen as radio lobes. The point at which the jet impacts the far side of the cocoon is seen observationally as a hotspot. The plasma swept up from the jet injection causes these radio lobes to expand supersonically; this characterises the cocoon expansion, the evolution of which we now quantify.

Using the gas density profile described in Section 4.1, and assuming a cocoon axial ratio $R_{T}$ for a jet with kinetic power $Q_{\rm jet}/2$, the dynamical model of Kaiser & Alexander (1997) characterises the evolution of the cocoon radius as

[TABLE]

Here time, $t$, is normalised to a characteristic time-scale, $\tau=\left(\frac{r_{0}^{5}\rho_{0}}{Q_{\rm jet}}\right)^{1/3}$, while the dimensionless constant, $a_{D}$, follows from the work of Kaiser & Alexander (1997) and Shabala & Alexander (2009):

[TABLE]

Assuming a power law density profile of $\rho=\rho_{0}(r/r_{0})^{-\beta}$ for the cocoon expansion where $0<\beta<2$, and axial ratio between $1.3<R_{T}<6$, the value for $a_{D}$ is in the range of 0.8-1.8 (see Table 1). We adopt $a_{D}=1$, a value typical for radio sources in the inner regions of clusters.

4.2.2 The cocoon shock radius

Following Kaiser & Alexander (1997), the relation between the maximum radius of the shocked gas and the radius of radio emitting cocoon plasma is given by

[TABLE]

where the maximum radius is calculated using the AGN active time ($t_{on}$, describe below in Section 4.2.4) and the cocoon radius (Equation 13). $\lambda$ is dependent on the density profile, and we have assumed a self-similar solution for the shape of the cocoon,

[TABLE]

Observations show the existence of somewhat rare, extremely powerful AGN outbursts in massive galaxies (Rawlings & Jarvis, 2004; Miley & De Breuck, 2008; Shabala et al., 2011; Fabian, 2012). Their size extends many hundreds of kpc, and their life-times can be up to $10^{8}$ yr (Shabala et al., 2008). Radio sources sizes, denoted $D$, have a typical range of between 1 and $\sim$1000 kpc (Willott et al., 1999; Sadler et al., 2007). In our model the shocked gas expands up to $\sim 900$ kpc, as can be seen in Figure 2, approximately consistent with such observations. Similarly, in Figure 3 we show the distribution of shock radii as function of stellar mass. Together, these figures indicate that AGN feedback can clear hot and cooling gas out to large radii. In our model, larger radio sources should preferentially be found in more dense environments. On the one hand, denser environments imply slower expansion speeds (through Equation 13); on the other hand, more massive black holes found in such environments have both higher jet powers (Equations 11, 12) and are longer lived (Equation 23). As shown in Figure 2, the net effect is for larger radio sources to preferentially inhabit dense environments in our model. Observationally, many large double-lobed radio galaxies are found in poor environments; on the other hand, many small-scale jets are also associated with Seyfert galaxies (e.g. Ulvestad & Wilson, 1984; Kaviraj et al., 2015). The apparent compactness of many radio AGN may at least in part be due to the lack of sensitivity to diffuse low-surface brightness features in FR-I radio galaxies (Shabala et al., 2017). Our simple model cannot address this selection effect, which is deferred to future work. Given the simplicity of our model we resist the temptation to make more quantitative predictions than this, but it is an area for future work. The feedback is done in two ways, through heating of the gas, and uplifting of the gas, which we describe in the following sections.

4.2.3 Gas heating and cooling

According to Alexander (2002), and following the study by Shabala & Alexander (2009), the difference between the mean isothermal temperature in the shocked gas and the post-shock temperature can be described by

[TABLE]

where the shock velocity is given by $\dot{r}_{\rm shock}={3\over 5-\beta}{r_{\rm shock}\over t_{\rm on}}$, where $t_{\rm on}$ is the active time of the AGN and is defined below in Section 4.2.4. We can then obtain the following integrated hot gas mass using the Makino et al. (1998) profile for a spherical shocked mass distribution, given by

[TABLE]

The shocked mass can never exceed that of the mass of hot gas. When the gas is shock heated we assume an instantaneous mixing of this gas with the hot gas of the halo, resulting in a new hot halo temperature of

[TABLE]

Figure 4 shows the distribution of this new hot gas temperature in units of the original virial temperature, given as a function of stellar mass. In general, AGN in our model act to raise the temperature of the gas by between a few to $20\%$, depending on the stellar mass of the central galaxy in the halo. In this model, we only change the temperature if $T_{\rm shock}>T_{\rm vir}$.

Hot halo gas is not static but cools over time. The cooling time of this gas at each radius is taken as the ratio of its specific thermal energy to the cooling rate per unit volume,

[TABLE]

where $\Lambda(T,Z)$ is the cooling function and depends on the temperature, $T$, and metallicity, $Z$. In such models a cooling radius $r_{\rm cool}$ is then defined as the radius at which $t_{\rm cool}$ is equal to a characteristic time-scale of the system, such as its age, time since last major merger, or dynamical time. Following the standard implementation of SAGE we take the last of these.

For the so-called hot accretion mode, when $r_{\rm cool}<R_{\rm vir}$, the cooling rate of the hot gas is calculated using

[TABLE]

where $\rho_{g}(r_{\rm cool})$ is the Makino et al. (1998) profile gas density at the cooling radius. With this, the cooling rate can be rewritten as

[TABLE]

Alternatively, in the so-called cold accretion mode, when $r_{\rm cool}>R_{\rm vir}$, the hot halo never forms and any infalling gas captured in the dark matter halo potential falls towards the centre on a free fall time-scale, $R_{\rm vir}/V_{\rm vir}$, given by

[TABLE]

4.2.4 The AGN duty cycle and its effect on cooling

Left unchecked, massive dark matter halos in particular can collect prodigious amounts of hot gas which will lead to runaway cooling at rates that are unsupported by the observations. Our AGN model provides an energy counterbalance to such cooling through the heating resulting from an AGN jet. Over long time-scales (100’s of Myr to many Gyr), such energy injection can be approximated as uniform and constant, as assumed in Croton et al. (2006). However, a more realistic model will attempt to describe the intermittent nature of black hole accretion and their resulting outflows and properties.

To this end, suppose that we have a black hole of a given mass being fed gas from its surrounding medium through an accretion disk. We simplify the otherwise complicated resulting AGN duty cycle by breaking it into four primary parts: the time that the jet is on, $t_{\rm on}$, inflating the cocoon as described above; the switching off of the AGN and hence jet when accretion stops; the time it takes for the cocoon to disappear once the jet pressure has been removed, $t_{\rm return}$; and finally the subsequent period of quiescence after the cocoon has dissipated but before the next episode of accretion reignites the AGN and jet, $t_{\rm quiet}$. This cycle is illustrated in Figure 5, and we now look to quantify each phase using the previously established relations.

Observationally, the AGN fraction as a function of stellar mass was used by Turner & Shabala (2015) to estimate an expression for $t_{\rm on}$, approximately given by

[TABLE]

Characteristic time-scales for $t_{\rm on}$ range from between $2\times 10^{7}$ to $5\times 10^{8}$ years for stellar masses between $10^{10}$ and $10^{12}\rm M_{\odot}$. For very low mass galaxies this expression implies radio sources that are so small as to do no feedback.

Once support from the jet for the cocoon has been removed the time for it to collapse back down can be estimated using the free fall time from the shock radius, which we parametrise as

[TABLE]

Here, $c_{s}$ is the sound speed in the gas which we approximate using the virial velocity of the parent halo, $V_{\rm vir}$. The factor of two is an empirically chosen value we use to prevent overcooling in our model. Physically, the slow (slower than the sound speed) return of the shocked gas is due to processes we do not model here, such as the movement of buoyant bubbles once the jet switches off, and the slow mixing of the heated gas throughout the cluster (Basson & Alexander, 2003; Soker, 2015; , 2016).

The total time the AGN is off, $t_{\rm off}$, can be calculated from $t_{\rm on}$ using the observed duty cycle of AGN in the local Universe. Best et al. (2005) and Shabala et al. (2008) both measured the radio loud AGN fraction as a function of galaxy mass; using Figure 3 of Best et al. (2005) the duty cycle can crudely be approximated as

[TABLE]

from which $t_{\rm off}$ is simply

[TABLE]

It then follows that $t_{\rm quiet}$ is just the difference between the total AGN off time and the cocoon return time,

[TABLE]

The presence of inflated cocoons from the AGN is expected to modify the cooling rate in the hot halo in multiple ways. For our model we make a simple approximation and suppress cooling whenever a cocoon is present, with the degree of suppression proportional to the hot gas mass fraction between the shock radius and cooling radius. In other words, only gas that has not been overrun by the shock can cool. This approximation is reasonable, as the active phase corresponds to at most a third of the total cycle time (see Equation 25), and significantly less for lower mass galaxies. Furthermore, jet driven bubbles uplift the shocked gas further away from the galaxy centre, thus stopping it from cooling. Once the cocoon expansion stops, gas can re-fill the evacuated cavity. If this “return” time is less than the “off” time the system will experience a period of quiescence and cooling will be restored. Across the entire cycle the cooling rate may be further modified through the altered hot gas temperature (see Figure 4).

In summary, the time-averaged fraction of cooling in the presence of the AGN duty cycle is

[TABLE]

We include only the first term in the case of modified cooling due to jet heating that leads to higher hot gas temperatures, only the second term in the case of feedback and altered cooling from gas uplifting, and both terms in the full model. The new cooling rate is given by

[TABLE]

5 Results

5.1 Model constraints

Our AGN jet model makes three important changes to the previous SAGE radio-mode prescription. Namely, it uses a more realistic hot gas density profile to calculate cooling, the AGN jet can heat the hot gas to higher temperatures than before, and the combined effect of this new temperature and the jet inflated cavities act to alter the cooling rate in a more realistic way. To calibrate this model we use a set of observables that the output must reasonably compare with, and in a physically sensible way. Only then can we extend our analysis to explore predictions and consequences that can be compared with future observations.

To start, in Figure 6 we show the cooling rates from our model at z=0 as function of their hot halo gas temperature (Equation 18). We compare with a number of X-ray observations of gas surrounding local galaxy clusters and groups: Ponman et al. (1996) and Bharadwaj et al. (2015) who evaluate the the bolometric luminosity within $r_{500}$; Peres et al. (1998) who measure the luminosity inside the cooling radius, and Anderson et al. (2015) who use stacked observations of X-ray emission around locally brightest cluster galaxies. This range of observations serve as a good general indicator of how our new model compares. As the gray model data points and solid median line reveal, our cooling–temperature relation provides a reasonably good fit to the observations, except perhaps for the hottest cluster halos, where we tend to under-predict cooling somewhat. The jet model fares considerably better than the original Croton et al. (2006) and Croton et al. (2016) models, which over–and under–predict these observations by a significant amount, respectively (Croton et al., 2016, see e.g. their figure 7).

The consequence of cooling is a pooling of cold gas in the galactic disk which then leads to star formation (see e.g. Croton et al., 2006). A key motivator to include AGN feedback in semi-analytic models is to understand how the suppression of this infalling gas then leads to a suppression of the formation of new stars. The effectiveness of this feedback, and at the right mass scale, is usually quantified using the galaxy stellar mass function, the shape of which is moulded by both AGN and supernovae. As in previous works, we use the high mass end and knee of the local stellar mass function to constrain our jet model (the low mass end is constrained by other feedback processes and not the focus of the present work).

Figure 7 shows the effect of the jet model on the formation of high mass galaxies compared to the SDSS observations of Baldry, Glazebrook & Driver (2008). First we present the stellar mass function for the original SAGE model without any AGN feedback (gray dashed line), then systematically add the changes described in the previous sections until we arrive at our default version. This version includes the updated hot gas density profile (blue short-dashed line), and both gas heating by the AGN jet (red long-dashed line) and also mass transport through uplifting (blue solid line, our final model).

Lastly, we consider the star formation rate density evolution of the galaxy population with time. In Figure 8 we plot the final model by the blue line and two bounding cases: the original SAGE model (magenta long-dashed line) and our model but with AGN feedback turned off (gray short-dashed line). As can be seen, both SAGE and our new jet model matches the observational compilation of Somerville et al. (2001) reasonably well. This is in contrast to when AGN feedback is turned off, which overproduces star formation at all redshifts except the highest.

5.2 Deriving radio luminosities

We now consider the AGN properties of our model galaxies, a new feature of this work, with a focus on comparing with observables. In particular, we look at the cocoon shock size, jet power, and 1.4GHz radio luminosity.

In our model the jet is assumed to expand into a hot atmosphere with a Makino density profile, characterised by core radius $r_{0}$, core density $\rho_{0}$, and power-law outer slope $\beta_{\rm eff}$ (see Section 4.1). The resulting radio lobes can be seen by their synchrotron emission. This depends on the cocoon magnetic field energy density, $u_{B}$, the distribution of Lorentz factors $\gamma$ of the emitting electrons, $n(\gamma)$, and cocoon volume, $V$. The energy due to the emission of synchrotron radiation is mostly injected at the hotspot between a minimum and maximum Lorentz factor, $\gamma_{1}$ and $\gamma_{2}$, with the initial energy distribution given by

[TABLE]

The electrons inflate the radio-emitting cocoon and lose energy through adiabatic expansion, synchrotron emission, and inverse Compton upscattering of Cosmic Microwave Background photons.

Following Kaiser & Best (2007) and Shabala & Godfrey (2013), under the assumption of equipartition between particle and magnetic field energy density the cocoon radio luminosity at frequency $\nu$ can be written

[TABLE]

where $p$ is the power-law index of the electron energy distribution, and $\alpha=(1-p)/2$. Using the definition of the spectral index for synchrotron emission, $F_{\nu}\propto\nu^{\alpha}$, a reasonable hotspot value for radio galaxies with $\alpha$ between $-0.55$ and $-0.85$ is $p=2.1-2.7$.

$A_{1}$ is a numerical constant given by

[TABLE]

where $r_{e}$ is the classical electron radius, $q$ and $m_{e}$ are the electron charge and mass respectively, $\mu_{0}$ is the permeability of free space, and $f(p,\gamma_{1},\gamma_{2})=\int_{\gamma_{1}}^{\gamma_{2}}\gamma^{1-p}d\gamma$. The function $C_{2}(p)$ can be expressed as

[TABLE]

Here, $\Gamma_{\rm fn}(z)=\int_{0}^{\infty}t^{z-1}e^{-t}dt$ is the Gamma function (, 2009), and for $p=2.6$, $C_{2}(p)=2.04\times 10^{-3}$. The results explored in this paper assume the best values described in Table 1.

For self-similar expansion, the cocoon volume is $V=\pi R_{\rm T}^{2}D^{3}$, where $R_{\rm T}$ is the axial ratio of the source taking values between 1.3 and 6 (Leahy & Williams, 1984; Leahy et al., 1989), and $D=2r_{\rm shock}$ is the diameter of the cocoon. It can be shown that the magnetic field energy density, which is proportional to the cocoon pressure, can be written (Kaiser & Alexander, 1997; Kaiser & Best, 2007, 2008; Shabala & Godfrey, 2013) as

[TABLE]

where $r=\frac{p+1}{4}$, and the constant $f_{p}$ is taken from Kaiser & Best (2007),

[TABLE]

The radio luminosity is then

[TABLE]

Here, $Q_{\rm jet}$ is the total kinetic power from both radio jets.

For a given $t_{\rm on}$, higher jet powers yield both larger sizes (Equation 13 and Equation 14) and radio luminosity (Equation 35). In Figure 9 we show the predicted luminosity–size relation, weighted by AGN active time, $t_{on}$. This figure illustrates that most sources have predicted sizes between 30-100 kpc, and luminosities $<10^{23}W/Hz$. In calculating the radio luminosities below, we adopt the final luminosity at age $t_{\rm on}$, rather than sampling different points along the size-luminosity track. We adopt this simple treatment due to the limited time resolution (260 Myr) of the SAGE output. In future work, we plan to include a more accurate sampling of the luminosity–size tracks. Furthermore, we have made the standard assumption that the synchrotron emitting electrons in the radio lobes are in the minimum energy state, i.e. there is approximate equipartition between the magnetic field and particle energy densities. Observational evidence (e.g. Croston et al., 2005; Shelton et al., 2011; Hardcastle & Krause, 2014) suggests that most sources have sub-equipartition magnetic fields, and the departure from the equipartition value depends on a number of factors, including radio source morphology. Croston et al. (2005) quote median magnetic field values that are marginally sub-equipartition, $0.7B_{\rm eq}$. Adopting this in our work would decrease the radio luminosities by a factor $u_{B}^{(5+p)/4}\sim 0.5$, i.e. a factor of two lower.

5.3 Radio AGN properties

Figure 10 shows the relation between jet power, $Q_{\rm jet}$, and 1.4 GHz radio luminosity. Our results are consistent with that found by Cavagnolo et al. (2010) and Heckman & Best (2014) given by the solid and dashed black lines, respectively. We note that Heckman & Best (2014) estimate the jet mechanical energy of radio sources from the observed cavities and bubbles in the X-ray gas.

In Figure 11 we present the AGN radio luminosity function (RLF) for our jet model in two redshift intervals: z=0 and z=1. This figure indicates that our jet model, using the parameters given in Table 1, produces a positive evolution of the RLF with redshift, and compares well with the observed data of Best & Heckman (2012) and Best et al. (2014) across the redshift range probed, including the increase in the number density of the most luminous sources. To make this figure we have assumed that all sources have the luminosity given at the end of their lifetimes, which is an underestimate (sources typically lose luminosity as they age). Our model has a number of free parameters but is calibrated only at z=0. The fact that it produces roughly the right abundances at z=1 is encouraging. Of course, the model is build from parts which approximate highly complex physical phenomena, and in the future we plan to use better ones - for example a more holistic density profile for cooling and AGN expansion. Such sophisticated approaches to predicting the radio AGN properties (e.g. Turner & Shabala, 2015) are deferred to future work.

In Figure 12 we investigate the observed radio luminosity as a function of stellar mass, again comparing our model against the data of Best & Heckman (2012). Here the gray symbols show the model AGN distribution, with the line and error giving the median and scatter. As can be seen, our model is more aligned with the observed properties of Low Excitation Radio Galaxies (LERGs) across most of the stellar mass range plotted, which are fuelled by cooling flows, although at the lowest masses our model galaxies could be made up of a mix with High Excitation Radio Galaxies (HERGs). The Best & Heckman (2012) sample extends to z=0.3, which corresponds to a completeness of approximately $10^{24}W/Hz$ at 1.4 GHz in radio luminosity (using the 99 percent NVSS completeness cutoff of 3.4 mJy), and $10^{11}M_{\odot}$ in stellar mass (using a SDSS r-band magnitude limit of 17.77 and the relations of Bell et al., 2003). These observations are therefore likely to be complete in stellar mass for the bulk of the bright LERG population, however not for HERGs where the observed median masses should be treated as upper limits.

6 Discussion and Conclusions

In this work we present a new prescription for AGN feedback that models in a more realistic way the intermittent nature of black hole accretion in the galaxy population and its resulting outflows. This update is built upon the SAGE semi-analytic galaxy formation model described in Croton et al. (2016), which itself was an update to the model that introduced the original “radio mode” semi-analytic prescription in Croton et al. (2006). Our enhancements extend to gas cooling from the hot halo to make the cooling–heating cycle self consistent, and include a new hot gas density profile that is better aligned with that observed in X-ray clusters. Overall, and after some minor parameter retuning, our additions retain the good fit SAGE has to a number of key galaxy population statistics, like the stellar mass function (for which AGN heating is critical at the high-mass end) and star formation rate density evolution, while expanding the number of observables it can make predictions for, in particular for AGN and radio galaxies.

Unlike most semi-analytic approaches, the AGN energy injection in our model is spatially distributed. Omma & Binney (2004) pointed out that, for large radio sources, much of the AGN kinetic power can be wasted by coupling to hot gas at large cluster-centric radii, which has long cooling times. By modelling the dynamical evolution of jet-inflated structures explicitly, we address this issue.

In this new model AGN and its resulting feedback are followed through a number of key phases with time: (1) black hole accretion acts to turn on a jet, which quickly expands into the surrounding hot halo gas and inflates a cocoon; (2) after a period of time the accretion stops, and hence also the jet, which removes pressure support from the cocoon leading to its deflation (return), and ultimately dissipation; and finally (3) the galaxy then undergoes a period of quiescence, which lasts until accretion onto the black hole is re-established, after which the cycle begins again. Importantly, only during the quiet phase do we allow cooling of gas onto the galaxy, with the duration of this time set by our prescriptions for the AGN jet on, off, and cocoon return times.

Within our model gas heating from an AGN jet acts to raise the mean temperature of the hot halo by up to $\sim 20\%$ above the virial temperature (Figure 4). When compared to the previous SAGE model, we obtain a superior match to the observed cooling luminosity–X-ray hot gas temperature relation, although we tend to under-predict the observations for the most massive clusters (Figure  6). Our more realistic cooling rates follow from the combined effect of the modified hot gas temperature and the influence that the jet inflated cavities have on the movement of this gas.

New to our model is its ability to make predictions for the properties of radio galaxies. We explore the relationship between jet power, shock radius, and radio luminosity, and compare with current observations (Figures  9 and 10). The model produces a good match to both the observed radio luminosity functions from $z=0$ to $z\sim 1$ (Figure  11), and the local stellar mass–radio luminosity relation (Figure 12). For this later result, our AGN more closely align with the Low-Excitation Radio Galaxy (LERG) population, rather than the High-Excitation (HERG) population, when compared with the observations of Best & Heckman (2012). Indeed, we can match the stellar mass function in a number of ways, but to match the radio luminosity function as well we find our free parameters, as given in Table 1, must be chosen to produce AGN with LERG-like properties.

However our model is not without its simplifications and areas for improvement. For example, the model only follows the active phase of jet inflated cavities and not the rising buoyant bubbles that can influence the hot gas long after the AGN has switched off (as seen from X-ray cavities with radio emission). Thus, our sources may better represent a population of FR-IIs, and our model may underestimate the total amount of feedback done by the AGN. In addition, we ignore the possibility of sound waves distributing heating to large radii. On the other hand, we have also assumed that all cavities have simple morphology with a covering factor of 1. In reality cavities show complicated shapes, with material constantly falling back in behind. As such, these assumptions will tend to overestimate the efficiency of feedback. Alexander (2002) considered the evolution of the swept-up shell of shocked gas, and showed that the cooling time of the swept up gas depends sensitively on the thermal conductivity of the cluster. In the case that the shocked gas evolves isothermally (as would be the case for clusters in which thermal conductivity is not suppressed), radiative cooling was shown to only be significant in the centres of cool core clusters. On the other hand, the shocked gas can cool rapidly if it evolves adiabatically. Our results suggest that thermal conductivity in clusters may provide an important feedback mechanism, as previously suggested by other authors (e.g. Narayan, 2002; McCourt et al., 2013).

Finally, it is important to note that although there is a consistency between our model predictions and observations in the $Q_{\rm jet}-L_{\rm radio}$ relation (Figure 10), this relation derived from observations of X-ray cavities suffers strongly from selection effects (Godfrey & Shabala, 2016). A relation based on dynamical models of radio AGN (Turner & Shabala, 2015) is not subject to such selection effects, but most likely provides a lower limit on the jet power for a given radio luminosity because of the assumption of an equipartition magnetic field in the radio cocoon. Inverse-Compton observations (e.g.; Croston et al., 2005; Shelton et al., 2011; Hardcastle & Krause, 2014) suggest that this assumption may overestimate the magnetic field strength by as much as a factor of four. According to Equation 33, a lower magnetic field would give a higher $Q_{\rm jet}$ for the same radio luminosity.

Despite these shortcomings this paper presents a new effort to incorporate the effects of AGN jet and cavity production into a cosmological model of galaxy formation. Before we can make predictions for future radio  surveys  we first need to have a well constructed cosmological radio AGN model that produces a population which matches a baseline set of key statistics. The majority of our figures are calibration plots, demonstrating that our  model construction and parameter  choices by-and-large  achieve this  (for the properties considered at least).

This model is able to produce catalogues of both normal and radio galaxies across a wide range of mass scales and environments, describing their evolution from high redshift down to the present day. Such catalogues offer a valuable resource for survey teams as they analyse current data and plan for the future. Extensions to our current work will focus on the geometry of bubble evolution (e.g. how do bubbles “pancake” and sweep out the gas), how gas refills behind the bubble as it evolves, and how the local magnetic field helps stabilise the bubble and delays gas from returning. To achieve such detail we will employ more sophisticated high resolution numerical simulations, with the results then generalised and folded into our semi-analytic model, expanding its predictive power and usefulness to the community.

Acknowledgements

We thank the anonymous referee for their constructive comments and suggestions which helped to improve the paper. MR is grateful to Swinburne University and its Centre for Astrophysics and Supercomputing for hosting his stay during which much of this work was carried out. SS thanks the Australian Research Council for an Early Career Fellowship (DE130101399).

The Semi-Analytic Galaxy Evolution (SAGE) model that served as a basis for this work is a publicly available codebase that runs on the dark matter halo trees of a cosmological N-body simulation. It is available for download at https://github.com/darrencroton/sage. The Millennium Simulation on which the semi-analytic model was run was carried out by the Virgo Supercomputing Consortium at the Computing Centre of the Max Plank Society in Garching. It is publicly available online at http://www.mpa-garching.mpg.de/Millennium/ through the German Astrophysical Virtual Observatory.