© ESO, 2010

1. Introduction

Elliptical galaxies invariably contain central supermassive black holes (SMBHs), and there exists a tight relationship between the characteristic stellar velocity dispersion σ (or stellar mass M_∗) of the host system and the SMBH mass M_BH (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000; Tremaine et al. 2002; Yu & Tremaine 2002). These relations clearly indicate a co-evolution of the SMBHs and their host spheroids. Several investigations have been dedicated to this subject, either by using hydrodynamical simulations (e.g., Ciotti & Ostriker 1997, 2001, 2007; Ciotti et al. 2009, 2010) or one-zone models (e.g., Sazonov et al. 2005,hereafter SOCS; Ballero et al. 2008, Matteucci 2008). In some work, the effect of galaxy merging has also been taken into account (e.g., Hopkins et al. 2005, 2006). The main advantage of hydrodynamical models is that complex physical phenomena effects (such as shock waves, jets, radiation transport, etc.) can be taken into account. However, the computational time of the simulations force us to search for faster methods that allow a more systematic exploration of the parameter space, which is prohibitive with hydrodynamical simulations. In this framework, hydro-simulations can be used to set the acceptable range for parameters to be adopted in toy models (e.g., the duty cycle value, see Sect. 2.1), while simpler models, such as that used here, are useful toidentify the most interesting cases that can be simulated in detail with hydrodynamical codes.

The general idea behind one-zone models is to work with “average” equations that capture some aspect of a more complicated situation. In practice, some of the equations are exact (such as for example the mass and energy balance equation). On the other hand, some of the physical variables are volume or mass averaged (e.g., the mean gas and temperature of the interstellar medium, respectively), or finally computed at some fiducial radius of the assumed gas distribution. From this point of view, the specific galaxy and dark matter halo profiles do not enter directly into the code, but are needed to obtain realistic mean values to be used in the equations.

A preliminary investigation of a physically based one-zone toy model has been done in SOCS, where it is shown that for a typical quasar spectral energy distribution (SED) the final M_BH produced by feedback clearly reproduces the observed M_BH-σ relation. The new models discussed here contain important improvements with respect to SOCS. First of all, the modelization of accretion onto SMBH. In this paper we examine how the effects of radiation pressure due to the Thomson electron scattering modify in a self-consistent way the spherical Bondi flow. Another improvement is the treatment of type Ia supernova rates and mass losses due to the evolution of the galactic stellar population, as we now adopt the Kroupa (2005) initial mass function (IMF) coupled with the evolutionary prescription of Maraston (2005). Finally, the dark matter halo is now described by a finite-mass Jaffe (1983) distribution instead of a singular isothermal sphere.

We focus on a scenario in which the masses of the central SMBH and the host galaxy grow in a dark matter halo, which is replenished by accretion of gas of cosmological origin. We follow star formation and we also consider the mass return from the evolving stellar populations. The combined effect of SNIa heating and radiative feedback, during episodes when the luminosity from the central SMBH approaches its Eddington limit, heats and drives much of the remaining gas out of the galaxy, limiting both future growth of the SMBH and future star formation to low levels. We do not consider the merging phenomenon, and we restrict ourselves to the evolution of an isolated galaxy: it is well known that significant AGN activity may also be present in isolated systems, because of the large amounts of gas produced by passively evolving stellar populations (e.g., Mathews 1983; Ciotti et al. 1991; Bregman & Parriott 2009; Kaviraj et al. 2010; Cisternas et al. 2010). More specifically, the one-zone models discussed in this paper are developed as an alternative approach to studying the “feedback modulated accretion flows” of Ciotti et al. (2009, 2010). There are two major differences between the present approach and hydro-simulations. The first, and most obvious, is the impossibility to resolve physical phenomena associated with specific length and timescales. The second is that we consider heating feedback only, while in the current version of the hydrodynamical code by Ciotti et al. (2009, 2010), both mechanical feedback, and radiation pressure are also considered. However, despite these shortcomings, with the present models we can attempt to simulate the process of galaxy formation, which is beyond the current possibilities of the hydro-simulations and at the same time we can explore the parameter space.

The evolutionary scenario that we consider here addresses several key observational findings. First, that giant ellipticals are old – they end their period of vigorous star formation early in cosmic time, since the radiative output from the central SMBHs limits (in cooperation with the energetic input due to star formation) the gas content to be at levels for which ongoing star formation is minimal. Secondly, gas-rich, actively star-forming galaxies at redshift z ~ 3, including Lyman break galaxies and bright submillimeter SCUBA galaxies, generally exhibit AGN activity (Steidel et al. 2002; Alexander et al. 2003; Alexander 2009; Lehmer et al. 2005; Donley et al. 2010), indicating that their central SMBHs continue to grow. This suggests that the formation of a spheroid probably closely preceeds a quasar shining phase, as verified by spectroscopic observations indicating that quasars occupy metal-enriched environments (e.g., Hamann & Ferland 1999; Schawinski et al. 2009; Wild et al. 2010). The redshift evolution of the quasar emissivity and the star formation history of spheroids is thus expected to have evolved roughly in parallel since z ~ 3, which is also consistent with observations (e.g., Haiman et al. 2004; Heckman et al. 2004). Among the most important observational predictions of the model is the length of the so-called “obscured accretion phase” (e.g., Comastri 2004). This phase is defined as the period of time when a high column density is associated with a high accretion rate onto the central SMBH. We also study the relation between the duration of the obscured phase and the corresponding “cold phase” (defined by a low mass-weighted gas temperature), and how they depend on the adopted parameters. Finally, in a large set of models we explore how the final SMBH mass is related to the final stellar mass, as a function of the dark matter halo mass, the amount of cosmological infall and the accretion rate of the gas, the SMBH accretion radiative efficiency, and finally the initial SMBH mass. We find that the models are in close agreement with the observed Magorrian relation when we assume high efficiencies (ϵ ~ 0.1) and relatively low initial SMBH masses (M_BH ≲ 10⁶   M_⊙).

This paper is organized as follows. In Sect. 2, we present the physics adopted in the simulations, with special emphasis on the differences from and improvements on SOCS. Section 3 is devoted to the description of the adopted self-regulated accretion model, while in Sect. 4 we present our findings. The main results are summarized and discussed in Sect. 5, while several useful interpolating functions that we use to compute the stellar mass losses and the SNe Ia rate are presented in the Appendix.

2. The model

The galaxy model and the input physics adopted for the simulations have been improved with respect to SOCS, in particular the treatment of SNIa heating, the mass return rate from evolving stars, and the accretion onto the central SMBH. Several aspects of this model were described in SOCS and Ciotti & Ostriker (2007), and we present here only the relevant modifications.

2.1. The unchanged physics

For completeness, we briefly summarize the aspects of the input physics that remain unchanged with respect to SOCS. The differential equation for the gas mass balance is $\mathit{Ṁ}\mathrm{g}\mathrm{=}\mathit{Ṁ}\inf \mathrm{-}\mathit{Ṁ}\mathrm{\ast }\mathrm{-}{\mathit{f}}_{\mathrm{Edd}}\mathit{Ṁ}\mathrm{BH}\mathrm{-}\mathit{Ṁ}\mathrm{esc}\mathit{.}$(1)The first source term on the right-hand side (r.h.s.) describes the cosmological infall in a pre-existing (and time-independent) dark matter halo $\mathit{Ṁ}\inf \mathrm{=}\frac{{\mathit{M}}_{\inf }}{{\mathit{\tau }}_{\inf }}{\mathrm{e}}^{\mathrm{-}\mathit{t}\mathit{/}{\mathit{\tau }}_{\inf }}\mathit{,}$(2)where M_inf is the total gas mass accreted during the simulation (but in general not equal to the final stellar galaxy mass M_∗), and τ_inf is the characteristic infall timescale. A more appropriate description of the cosmological gas infall could be obtained by multiplying the r.h.s. of Eq. (2) by a factor t/τ_inf (e.g., Johansson et al. 2009), so that the total infalling mass is unaffected, while the infall starts in a more gentle way. For completeness, we performed some simulation with this different description, and we found (as expected) that for reasonable values of τ_inf the results are not much affected. Therefore, for consistency with SOCS we retain Eq. (2). The second source term is the net star-formation rate $\mathit{Ṁ}\mathrm{\ast }\mathrm{=}\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{-}\mathit{Ṁ}\begin{array}{}\mathrm{w}\\ \mathrm{\ast }\end{array}\mathit{,}$(3)where $\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{\ast }}{\mathit{M}}_{\mathrm{g}}}{{\mathit{\tau }}_{\mathrm{\ast }}}$(4)is the instantaneous star-formation rate, and $\mathit{Ṁ}\begin{array}{}\mathrm{w}\\ \mathrm{\ast }\end{array}$ is the mass return by the evolving stellar population (Appendix A.1); following SOCS, in the simulations α_∗ is fixed to be 0.3. The characteristic star-formation timescale τ_∗ is defined as ${\mathit{\tau }}_{\mathrm{\ast }}\mathrm{=}\begin{array}{}\max \end{array}\mathrm{\left(}{\mathit{\tau }}_{\mathrm{dyn}}\mathit{,}{\mathit{\tau }}_{\mathrm{cool}}\mathrm{\right)}\mathit{,}$(5)where the dynamical timescale τ_dyn is given by Eq. (18), and the mean gas cooling time τ_cool is estimated to be ${\mathit{\tau }}_{\mathrm{cool}}\mathrm{=}\frac{\mathit{E}}{\mathit{Ė}\mathrm{C}}\mathrm{·}$(6)In the equation above, Ė_C is the fiducial cooling rate given by $\mathit{Ė}\mathrm{C}\mathrm{=}{\mathit{n}}_{\mathrm{e}}{\mathit{n}}_{\mathrm{p}}\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\left(\frac{{\overline{\mathit{\rho }}}_{\mathrm{g}}}{{\mathit{m}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\frac{\mathit{X}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{X}\mathrm{)}}{\mathrm{2}}\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{,}$(7)where ${\overline{\mathit{\rho }}}_{\mathrm{g}}$ is the instantaneous mean gas density (see Sect. 2.3), $\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2.18}\mathrm{\times }{\mathrm{10}}^{-18}}{{\mathit{T}}^{\mathrm{0.6826}}}\mathrm{+}\mathrm{2.706}\mathrm{\times }{\mathrm{10}}^{-47}{\mathit{T}}^{\mathrm{2.976}}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1} \mathrm{cm}$(8)is the cooling function (Mathews & Bregman 1978; see also Ciotti & Ostriker 2001), and X = 0.7 is the hydrogen mass abundance (for simplicity we assume complete ionization). In agreement with Eqs. (6) and (7), the mean gas internal energy is $\mathit{E}\mathrm{=}\frac{\mathrm{3}{\overline{\mathit{\rho }}}_{\mathrm{g}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathrm{2}\mathit{\mu }{\mathit{m}}_{\mathrm{p}}}\mathit{,}$(9)where μ = (0.25 + 1.5X + 0.25Y)^-1 ≃ 0.62 is the mean atomic weight. It follows that ${\mathit{\tau }}_{\mathrm{cool}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }{\mathit{r}}_{\mathrm{g}}^{\mathrm{3}}{\mathit{k}}_{\mathrm{B}}{\mathit{m}}_{\mathrm{p}}}{\mathit{\mu X}\mathrm{(}\mathit{X}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{M}}_{\mathrm{g}}}\frac{\mathit{T}}{\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathit{,}$(10)where M_g is the instantaneous value of the total gas mass and r_g is the gas distribution scale radius (Sect. 2.3). The third source term is the total accretion rate onto the SMBH $\mathit{Ṁ}\mathrm{BH}\mathrm{=}\mathit{Ṁ}\mathrm{BH}\mathit{,}\mathrm{acc}\mathrm{+}{\mathit{\beta }}_{\mathrm{BH}\mathit{,}\mathrm{\ast }}\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathit{.}$(11)The first term describes gaseous accretion (Sect. 3), whereas the second term represents the contribution by the coalescence of stellar remnants of massive stars, as discussed in SOCS. The numerical coefficient f_Edd ≈ 10^-2−10^-3 needs to be implemented in the one-zone code (which by definition is unable to model the different spatial scales of the problem) to represent the observed time variation of quasars. In practice, f_Edd represents the “duty-cycle”, and its value is constrained by both observations (e.g. Yu & Tremaine 2002; Haiman et al. 2004) and simulations (Ciotti & Ostriker 2007; Ciotti et al. 2009, 2010). When the thermal energy of the interstellar medium (ISM) of the galaxy is high enough, the gas is able to escape from the dark matter potential well, at a fiducial escape rate computed as $\mathit{Ṁ}\mathrm{esc}\mathrm{=}\{\begin{array}{c}\\ \frac{{\mathit{M}}_{\mathrm{g}}}{{\mathit{\tau }}_{\mathrm{esc}}}\mathit{T}⩾{\mathit{\eta }}_{\mathrm{esc}}{\mathit{T}}_{\mathrm{vir}}\mathit{,}\\ \\ \mathrm{0}\mathit{T}\mathit{<}{\mathit{\eta }}_{\mathrm{esc}}{\mathit{T}}_{\mathrm{vir}}\mathit{.}\end{array}$(12)The parameter η_esc is of the order of unity, while the expression for T_vir is given in the next section. Finally, the escape characteristic timescale is ${\mathit{\tau }}_{\mathrm{esc}}\mathrm{=}\frac{\mathrm{2}{\mathit{r}}_{\mathrm{g}}}{{\mathit{c}}_{\mathrm{s}}}\mathit{,}$(13)where c_s is the speed of sound and r_g is the scale length of the gas distribution (Sect. 2.3). Energy input into the ISM is provided by the thermalization of supernova ejecta (both SNII and SNIa). The treatment of SNII is the same as in SOCS, and the new numerical treatment of SNIa is described in Appendix A.2. Additional contributions to the ISM energetics come from the thermalization of red giant winds and radiative feedback due to accretion onto the SMBH. As in SOCS, we adopt the average quasar SED obtained from the cosmic background radiation field supplemented by information from individual objects. We recall that the UV and high energy radiation from a typical quasar can photoionize and heat a low density gas up to an equilibrium Compton temperature (T_C ≈ 2 × 10⁷ K) that exceeds the virial temperatures of giant ellipticals. Following SOCS, we also consider adiabatic cooling in the case of gas escaping and heating/cooling due to inflow/outflowing galactic gas. The gas temperature is therefore determined at each time-step by integrating the equation of the internal energy per unit volume $\mathit{Ė}\mathrm{=}\mathit{Ė}\mathrm{H}\mathit{,}\mathrm{SN}\mathrm{+}\mathit{Ė}\mathrm{H}\mathit{,}\mathrm{w}\mathrm{+}\mathit{Ė}\mathrm{H}\mathit{,}\mathrm{AGN}\mathrm{-}\mathit{Ė}\mathrm{C}\mathrm{+}\mathit{Ė}\mathrm{ad}\mathrm{+}\mathrm{3}\frac{\mathit{Ṁ}\inf \mathit{\lambda }{\mathit{\upsilon }}_{\mathrm{esc}}^{\mathrm{2}}\mathrm{-}\mathit{Ṁ}\mathrm{esc}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{16}\mathit{\pi }{\mathit{r}}_{\mathrm{g}}^{\mathrm{3}}}\mathit{,}$(14)where Ė_H,SN is the energy due to SNIa and SNII, Ė_H,w describes the thermalization of red giant winds, Ė_H,AGN is the AGN heating, Ė_ad is the adiabatic cooling in the case of galactic winds, and λ is a dimensionless parameter (). Finally, we force the gas to remain above 10⁴ K.

2.2. The dark matter halo

For simplicity, the dark matter (DM) halo in our model is the only contributor to the gravity of the galaxy. In addition, the gas and the stellar density distributions in the simulations are assumed to be proportional to the local (unevolving) DM density distribution, modeled as a Jaffe (1983) profile (Sect. 2.3) ${\mathit{\rho }}_{\mathrm{h}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{h}}}{\mathrm{4}\mathit{\pi }}\frac{{\mathit{r}}_{\mathrm{h}}}{{\mathit{r}}^{\mathrm{2}}\mathrm{(}{\mathit{r}}_{\mathrm{h}}\mathrm{+}\mathit{r}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(15)For a DM halo total mass M_h, the scale length r_h is fixed following Lanzoni et al. (2004). In that paper, 13 massive DM halos obtained from high-resolution cosmological N-body simulations were carefully analyzed. In particular, for each halo the “overdensity radius” r_Δ (sometimes called virial radius) was determined, and it was shown that r_Δ does not differ significantly from the true virial radius r_vir of the system (with a scatter  ≲ 20%). Averaging over the 13 clusters reported in Table 2 in Lanzoni et al. (2004), we found that ${\mathit{M}}_{\mathrm{h}}\cong \mathrm{0.03}{\mathit{r}}_{\mathrm{\Delta }}^{\mathrm{3}}\mathit{,}$(16)where M_h is in 10⁹   M_⊙ and r_Δ is in kpc units. For the Jaffe profile r_vir = 2r_h, and from Eq. (16) we can link the scale-length r_h to M_h by assuming r_vir = r_Δ. We finally recall that the virial (3D) velocity dispersion associated with Eq. (15), is given by ${\mathit{\sigma }}_{\mathrm{vir}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{h}}}{\mathrm{2}{\mathit{r}}_{\mathrm{h}}}\mathrm{·}$(17)The halo mean circular velocity is estimated as ${\mathit{\upsilon }}_{\mathrm{c}}^{\mathrm{2}}\mathrm{=}\mathrm{2}{\mathit{\sigma }}_{\mathrm{vir}}^{\mathrm{2}}$, and the dynamical time required in Eq. (5) is therefore given by ${\mathit{\tau }}_{\mathrm{dyn}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{h}}}{{\mathit{\upsilon }}_{\mathrm{c}}}\mathrm{·}$(18)

2.3. The galaxy model

One of the new aspects of the model evolution investigated in this paper is the time extent of the Compton thin and Compton-thick phases, as defined at the end of this section. In SOCS, the adopted gas density distribution was a singular isothermal sphere and the major drawback of this choice for the evaluation of average quantities is the divergence of the total mass. For this reason, we now use the more realistic Jaffe density profile, i.e., we assume that at each time the gas density distribution is ${\mathit{\rho }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{g}}}{\mathrm{4}\mathit{\pi }}\frac{{\mathit{r}}_{\mathrm{g}}}{{\mathit{r}}^{\mathrm{2}}\mathrm{(}{\mathit{r}}_{\mathrm{g}}\mathrm{+}\mathit{r}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(19)where r_g is the gas scale radius and M_g is the instantaneous value of the total gas mass, obtained by integrating Eq. (1). For simplicity, r_g = r_h during the whole simulation. However, it is possible to generalize the present approach and consider a time-dependent r_g.

The mean gas density, which is needed in Eqs. (7) and (9), is evaluated to be the instantaneous mean value within the half-mass radius r_g${\overline{\mathit{\rho }}}_{\mathrm{g}}\mathrm{=}\frac{\mathrm{3}{\mathit{M}}_{\mathrm{g}}}{\mathrm{8}\mathit{\pi }{\mathit{r}}_{\mathrm{g}}^{\mathrm{3}}}\mathit{,}$(20)while the virial gas temperature T_vir is obtained from the hydrostatic equilibrium and the Jeans equations solved in the gravitational potential of the DM halo neglecting the self-gravity of the gas. Simple algebra then shows that ${\mathit{T}}_{\mathrm{vir}}\mathrm{=}\frac{\mathit{\mu }{\mathit{m}}_{\mathrm{p}}\mathit{G}{\mathit{M}}_{\mathrm{h}}}{\mathrm{6}{\mathit{k}}_{\mathrm{B}}{\mathit{r}}_{\mathrm{g}}}\mathrm{·}$(21)To quantify the relative importance of the optically thin and thick phases, the code computes the gas column density as $\mathrm{⟨}{\mathit{N}}_{\mathrm{H}}\mathrm{⟩}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{Pg}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathrm{2}\mathit{\pi }{\mathit{R}}^{\mathrm{2}}}\mathit{,}$(22)where M_Pg(R) is the projected gas mass enclosed in a circle of radius R, and the factor of 2 at the denominator takes into account that only one side of the gas column actually obscures the center. Therefore, the fiducial column density depends not only on the total gas mass, but also on the aperture radius adopted. To simulate observational work (e.g. Risaliti et al. 1999), we decided to define R to be one hundredth of the effective radius R_e¹, i.e., in Eq. (22) R ≅ 0.0074r_g . From Eq. (A.17) $\mathrm{⟨}{\mathit{N}}_{\mathrm{H}}\mathrm{⟩}\cong \mathrm{32.7}\frac{{\mathit{M}}_{\mathrm{g}}}{{\mathit{r}}_{\mathrm{g}}^{\mathrm{2}}}\mathit{,}$(23)and r_g is fixed during the simulation,  ⟨ N_H ⟩  depends only on M_g. We refer to the “thick” phase if  ⟨ N_H ⟩  ≥ 10²⁴ cm^-2, while galaxies with column density  ⟨ N_H ⟩  ≥ 10²² cm^-2 are considered obscured.

2.4. The code

The simulations are performed by numerical integration of the previous evolutionary differential equations using a forward Euler scheme. The numerical code is based on the code developed in SOCS. The time step is defined as the minimum among the characteristic times associated with the different physical processes $\mathit{\delta t}\mathrm{=}\mathit{B}\begin{array}{}\min \end{array}\mathrm{\left(}{\mathit{t}}_{\inf }\mathit{,}{\mathit{t}}_{\mathrm{\ast }}\mathit{,}{\mathit{t}}_{\mathrm{g}}\mathit{,}{\mathit{t}}_{\mathrm{BH}}\mathit{,}{\mathit{t}}_{\mathrm{E}}\mathrm{\right)}\mathit{,}$(24)where the different subscripts indicate the specific aspect of the physics involved and, in general, for a quantity X_i, t_i ≡ X_i/ | Ẋ_i | . The dimensionless coefficient B ≤ 1 is used to improve the code accuracy. We performed several test simulations and we found that B ≅ 0.1 leads to rapid convergence and excellent agreement with results obtained in SOCS, when used with their model galaxies and input physics.

3. Self-regulated Bondi accretion

As often assumed in the literature, the SMBH accretion rate is determined as the minimum between the Bondi accretion rate Ṁ_B and Eddington limit Ṁ_Edd$\mathit{Ṁ}\mathrm{BH}\mathit{,}\mathrm{acc}\mathrm{=}\begin{array}{}\min \end{array}\mathrm{\left(}\mathit{Ṁ}\mathrm{Edd}\mathit{,}\mathit{Ṁ}\mathrm{B}\mathrm{\right)}\mathit{.}$(25)In the previous formula, the Eddington accretion rate is $\mathit{Ṁ}\mathrm{Edd}\mathrm{\equiv }\frac{{\mathit{L}}_{\mathrm{Edd}}}{\mathit{ϵ}{\mathit{c}}^{\mathrm{2}}}\mathit{,}$(26)where 0.001 ≲ ϵ ≲ 0.1 is the accretion efficiency, and ${\mathit{L}}_{\mathrm{Edd}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi cG}{\mathit{M}}_{\mathrm{BH}}\mathit{\mu }{\mathit{m}}_{\mathrm{p}}}{{\mathit{\sigma }}_{\mathrm{T}}}$(27)is the Eddington luminosity and σ_T is the Thomson cross-section. The classical Bondi accretion rate is $\mathit{Ṁ}\mathrm{B}\mathrm{=}{\mathit{\lambda }}_{\mathrm{c}}\mathrm{4}\mathit{\pi } {\mathit{G}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{BH}}^{\mathrm{2}} {\mathit{\rho }}_{\mathrm{\infty }}{\mathit{c}}_{\mathrm{s}\mathit{,}\mathrm{\infty }}^{-3}\mathit{,}$(28)where ρ_∞ and c_s,∞ are the gas density and the speed of sound at infinity, respectively, and λ_c is a dimensionless coefficient of the order of unity (e.g., Bondi 1952; Krolik 1999). In Eq (28), radiative effects are not taken into account, so that in Eq. (25) there is a sharp transition between the pure hydrodynamical and radiation-dominated regimes. Fortunately, in the optical thin regime dominated by electron scattering, it is possible to extend the classical Bondi accretion solution and take into account the radiation pressure effects, so that the transition between Bondi-limited and Eddington-limited accretion can be described in a more consistent way. After solving this problem, we discovered that is had already been fully described by Taam et al. (1991) and Fukue (2001), and here we just summarize the main features of the modified accretion model. The result is obtained basically, by noting that in the optically thin regime the radiation pressure scales as 1/r², thus reducing the effective gravitational force at each radius by the same amount. By imposing self-consistentcy, i.e. by requiring that the effective accretion rate $\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}$ is determined by the effective gravity, one finally obtains the following expression for the modified coefficient ${\mathit{\lambda }}_{\mathrm{c}\mathit{,}\mathrm{eff}}\mathrm{=}{\mathit{\lambda }}_{\mathrm{c}}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}}{\mathit{Ṁ}\mathrm{Edd}}\right)}^{\mathrm{2}}\mathrm{·}$(29)Therefore, the self-consistent Bondi accretion rate satisfies the equation $\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}\mathrm{=}\mathit{Ṁ}\mathrm{B}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}}{\mathit{Ṁ}\mathrm{Edd}}\right)}^{\mathrm{2}}\mathit{,}$(30)that can be solved for $\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}$. After discarding the unphysical solution, we have $\mathit{ṁ}\mathrm{\equiv }\frac{\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}}{\mathit{Ṁ}\mathrm{Edd}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left[\mathrm{2}\mathrm{+}\mathit{r}\mathrm{-}\sqrt{\mathrm{4}\mathit{r}\mathrm{+}{\mathit{r}}^{\mathrm{2}}}\right]\mathit{,} \mathit{r}\mathrm{\equiv }\frac{\mathit{Ṁ}\mathrm{Edd}}{\mathit{Ṁ}\mathrm{B}}\mathrm{·}$(31)For r → 0 (high accretion rates) ṁ → 1, so that $\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}$ tends to the Eddington accretion, while for r → ∞ (low accretion rates) ṁ → 0. The full solution is indicated in Fig. 1 by the heavy solid line: we note how Eq. (25) overestimates accretion onto SMBH in the range .

4. Results

In the following sections, we present the main properties of a set of simulations, focusing on three important issues and comparing the new results with those in SOCS. However, before discussing the whole set of the new models, in the next section we present in detail the evolution of a representative model and three of its possible variants.

4.1. The evolution of the reference model and some of its variants

Figure 2 shows the time evolution of important quantities of our galaxy reference model (RM). This model is characterized by a DM halo of total mass M_h = 4 × 10¹¹      M_⊙, which corresponds to a halo scale length r_h = 11.86 kpc and R_e = 8.83 kpc, a fiducial circular velocity calculated according to Eq. (17) of  ~270 km s^-1, and a characteristic infall time 2 Gyr. The assumed total mass of the cosmological gas infall is M_inf = 10¹¹      M_⊙, corresponding to a dark-to-total mass ratio of 80% (provided that all the infalling gas forms stars). Other simulation parameters are α_∗ = 0.3, β_BH, ∗  = 1.5 × 10^-4, ϵ = 0.1, η_SN = 0.85, η_esc = 2, and a stellar mass-to-light ratio of 5 (see Eq. (A.11)). As in SOCS, the initial SMBH mass is 10⁸      M_⊙; the duty circle, f_Edd, is fixed at 0.005.

As can be seen from the comparison of Fig. 2 with Figs. 7 and 8 in SOCS, the global evolution of the new models is qualitatively very similar to those in SOCS. In particular, the time evolution of the model, from the beginning up to  ~9 Gyr, is characterized by a cold phase (we assume that a cold phase corresponds to T_gas ≤ 10⁵ K) of high density and low temperature. The gas density at the Bondi radius remains between 1 and 10² particles per cm³, while the mean gas density is  ~10^-1−10^-2 particles per cm³. At the beginning of the cold phase, about 2.17 Gyr are spent in the Compton-thick phase (see Fig. 2d). The remaining part of the cold phase is obscured, while the accretion luminosity remains high, at L_BH ≳ 10⁴⁶ erg s^-1. As soon as the mean gas density decreases (Fig. 2c green dotted line) the cooling becomes inefficient, the gas heating dominates, and the temperature increases to the virial temperature. The total duration of the cold phase is  ~9 Gyr, which corresponds to a decrease in  ⟨ N_H ⟩  of about 2 orders of magnitude. The obscured phase lasts  ~9.6 Gyr and the average column density is then around  ⟨ N_H ⟩  ~ 10²⁰ cm^-2. The durations of the cold and the obscured phases are very similar, and indeed, the two phases are related. Until the gas density is high, the gas can radiate efficiently the energy input due to the AGN, and its temperature remains below 10⁵ K. Because of this large amount of cold gas, the star formation proceeds at high rates, until most of the gas is consumed. At this point, the radiative cooling of the gas becomes inefficient and the heating due to the AGN causes the gas temperature to increase. The combined effects of the gas consumption and the increase temperature stop the star formation. A specific feature of the cold phase can be seen in Fig. 2, where characteristic temperature oscillations are apparent. These oscillations, already presented and discussed in SOCS, are due to the combined effect of gas cooling and AGN feedback. They terminate when the gas density falls below some threshold determined by the cooling function and because less and less gas is produced by stars and accreted by the galaxy, while SNIa heating declines less strongly. When the density is high (at early times), the cooling time is instead, very short, and AGN heating is radiated efficiently. The two competitive effects produce the temperature (and density) oscillations. In the hydro-simulations, the spatial and temporal structures of these oscillations is quite complicated, as feedback and cooling act on several different spatial and temporal scales (from a month to 10 Myr). In the one-zone models, the oscillations are instead dependent on the “duty-cycle” parameter f_Edd in Eq. (1), which is constrained by observational and theoretical studies (see Sect. 2.1).

To summarize, the observational properties of the RM would correspond to a system initially Compton-thick, that then switches to an obscured phase, and in the past 3.2 Gyr the galaxy is unobscured. The SMBH accretion history is shown in Fig. 3, where it is compared to that predicted by the SOCS formulation in a model otherwise identical to the RM. The major difference between the two models is in the final value of M_BH: in particular, Eq. (25) would predict a present-day M_BH a factor  ~1.5 higher than that for modified accretion. This is because the accretion in SOCS remains for almost all of the galaxy life, at the Eddington limit, while in the new treatment the self-regulation maintains the accretion at lower rates.

Before presenting the results of a global exploration of the parameter space, we focus on a few obvious questions. For example, what happens if we keep all the model parameters fixed and reduce the radiative accretion efficiency? What happens if we increase or decrease the total mass of infalling gas? Or if we double the gas infalling time? Of course, these simple examples do not cover all the possible cases. However, these models will provide a guide as we investigate the parameter space. The relevant properties of the additional “reference” models are listed in Table 1.

The first variant of the RM model is model RM₁ obtained by reducing the radiative accretion efficiency from ϵ = 0.1 to ϵ = 0.001. Overall, this reduced efficiency model passes through the same evolutionary phases as model RM: an initial Compton-thick phase, followed by an obscured phase, and finally a low-density unobscured phase. The initial cold high density phase is also very similar to that of RM in Fig. 2. The most important and expected difference is in the final mass of SMBH, which is higher by a factor of  ≈3.5 than the RM. A reduction in the radiative efficiency also produces a slightly larger amount of escaped gas, because of the shorter cold phase. As a consequence, a larger amount of gas produced by the evolving stars is lost from the galaxy during the low-density, hot phase.

The effect of a reduction in M_inf is explored in model RM₂, while all the others parameters are the same as in RM. Qualitatively, the evolution is again very similar to that of RM in Fig. 2. The main difference is the absence of the initial Compton-thick phase because of the lower gas density. Both the final stellar mass and final SMBH mass are lower than in model RM, as can be seen from Table 1. Overall, model RM₂ does not display remarkable or unexpectedproperties. The only noticeable aspect is that the escaped mass is higher than in model RM, but this is again due to the shorter cold phase.

A complementary model to RM₂ is RM₃, where we double the value of M_inf while maintaining all the other parameters identical to those of model RM. The qualitative evolution is again similar to that of RM in Fig. 2, but, at variance with RM₂, the Compton-thick phase is now present. Not surprisingly, the Compton-thick phase in model RM₃ extends for a longer period (~5 Gyr) than in the other variants of RM (see Table 1), and its total mass of new stars is also the highest. However, the total mass that escapes is not as high as one would expect, the infall mass being a factor of 2 higher than in the RM yet the escaped mass being almost the same. The main reason for this behaviour is the very massive star formation, which is almost double that of RM. We note that model RM₃ is the closest, in the RM family, to the galaxy models studied in the hydrodynamical simulations of Ciotti et al. (2009, 2010), as far as the final stellar mass and the final M_BH are concerned.

We finally discuss model RM₄, which is identical to model RM, but has twice as long an infall time. The main effects are the longest cold phase in RM family (see Table 1), and the absence of the initial Compton-thick phase. This latter characteristic is due to the time dilution of the infalling gas density, which prevents the possibility of reaching very high column density values. Overall, the final SMBH mass is not affected significantly by the extended infall phase, but its escaped gas masses is quite high because the total mass in new stars is (as expected) lower than in the RM model.

To summarize, these preliminary experiments have revealed that, at fixed dark matter halo, sensible variations in the input parameters do not produce remarkable different results. The relevant differences are found mainly in the amount of star formation and the different durations of the cold and obscured phases, while the final SMBH mass appears to be mainly affected by the value of the radiative efficiency ϵ.

4.2. Exploring the parameter space

We now present the general results (summarized in Tables 2 and 3) of our examination of the parameter space. The SMBH initial mass in the models in Table 2 is 10⁸   M_⊙, while in Table 3 it is 10⁵   M_⊙. From the astrophysical point of view, the first choice mimics a scenario in which the central SMBHs are already quite massive at the epoch of galaxy formation, while in the second case the main growth occurs with galaxy formation. Each of the 5 main families of models (M₁, ..., M₅) consists of galaxy models characterized by the same dark halo mass, M_h, ranging from 10¹⁰   M_⊙ (the M₁ family) up to 10¹²   M_⊙ (the M₅ family). In each family, we have investigated 6 models with different values of infalling gas-to-dark matter ratio α = M_inf/M_h, and finally different radiative efficiencies ϵ in the self-regulated Bondi accretion. In practice, in each family we explore the effects of different total infalling gas mass and different efficiencies (spanning the commonly accepted range of values) at fixed M_h.

Almost independently of M_h and the initial value of M_BH, we can recognize some general trends. For example in each group of 3 models characterized by identical parameters but decreasing ϵ, it is apparent how the total amount of stars formed decreases at decreasing ϵ, while the final M_BH increases. In the total mass budget of the galaxy, an important quantity is the total mass of gas ejected, M_esc. As expected, we find that massive galaxies, with larger infall gas masses also eject more mass. However, in all cases, the escaped gas mass is much lower than the final stellar mass, i.e.,  ~10% of it.

Finally, higher final gas masses are obtained, at fixed M_h, for higher M_inf and lower efficiencies, because of the less effective feedback. We note however that the product of ϵ and the mass accreted by the SMBH decreases with decreasing ϵ, i.e., even though the accreted mass is higher, the integrated energy output is lower, and so is the feedback effect².

A qualitative illustration of the effect of the reduction in ϵ on the SMBH accretion history is shown in Fig. 4. When the efficiency decreases, the final mass of the SMBH increases and, we indicate by the green line, the corresponding evolution of identical models but with the SOCS treatment, i.e., where the accretion is determined by the minimum of Ṁ_B and Ṁ_Edd. As expected, the differences in M_BH reduce with decreasing ϵ, because the Eddington accretion regime in the SOCS models (when the major differences from the Bondi-modified case are established), becomes less and less important, and $\mathit{Ṁ}\begin{array}{}\mathrm{e}\\ \mathrm{B}\end{array}$ approches Ṁ_B.

All models have a transition phase from obscured to unobscured. The Compton-thick phase is present in almost all models that have a cold phase (18 of 30 simulations); however, in the less massive set of models (M₁) the Compton-thick phase is absent. This is consistent with the scenario in which the majority of the most massive black holes spend a significant amount of time growing in an earlier obscured phase (see Kelly et al. 2010; Treister et al. 2010).

Finally, the mean SMBH-to-star ratio is not far from the value inferred from the present-day Magorrian et al. (1998) relation, though on the high mass side. This is not surprising, because of our use of a one-zone model. In any case, we emphasize that the SMBH feedback in the models was able to remove (in combination with star formation and/or gas escape) most of the infalling gas (which, if accreted onto the central SMBH in a cooling-flow like solution, would lead to the final SMBH masses being  ~2 orders of magnitude higher than the observed ones).

Independent of the particular characteristics of the single models, Fig. 5 clearly illustrates that the final SMBH masses, in models with a quite high initial M_BH (Table 2), are higher than implied by the observed Magorrian relation. For this reason, we explored other families of models, in which the initial mass of the SMBH was reduced to 10⁶   M_⊙, 10⁵   M_⊙, and 10³   M_⊙. From Fig. 5, it is apparent how galaxy models with initial SMBH masses  ≲ 10⁶   M_⊙ (circles) closely agree to within 1σ of the observed dispersion with observations, but only when the radiative efficiency is high, i.e. ϵ = 0.1. It is even more remarkable that the final SMBH mass is proportional to the final M_∗, as all the models started with the same initial SMBH mass. This strongly indicates that co-evolution is the most plausible explanation of the proportionality between M_∗ and M_BH, in line with other observational evidence (e.g., see Haiman et al. 2004, and references therein).

To compare our present findings with those of a complementary hydrodynamical approach, in Fig. 5 we plot the results of the high-resolution hydro-simulations developed by Ciotti et al. (2010, Table 1, Cols. 5 and 6) for a representative galaxy with stellar mass of  ~3 × 10¹¹   M_⊙. The final M_BH are represented by crosses, and the different (luminosity weighted) radiative efficiencies are in the range 0.003 ≤ ϵ ≤ 0.133. We note that the crosses indicate almost all the hydrodynamical models that have been studied in detail so far, and this shows the importance of one-zone models as a complementary approach to exploring the parameter space. We also note thatthe final mass of the SMBHs, as computed in the hydrodynamical simulations, lead to an accurate reproduction of the Magorrian relation but the comparison with the one-zone models is delicate. In the hydrodynamical simulations (at variance with the one-zone models), the initial phases of galaxy formation are not simulated, and the focus is on the maintaining low SMBH masses in the presence of stellar mass losses that, if accreted onto the central SMBH as in an undisturbed cooling flow, would produce a final BH mass of about a factor of  ~100 higher than the observed ones. For this reason, in the hydrodynamical models the galaxy is already assumed to have formed, and the initial mass of the SMBH is just slightly lower than the mass predicted by the Magorrian relation.

5. Conclusions

This paper is a natural extension of a previous paper by Sazonov et al. (2005). From a technical point of view, several aspects of the input physics have been improved. In particular, we have adopted a different description of the accretion rate, which now follows the self-consistent modified Bondi theory (Taam et al. 1991; Fukue 2001), instead of the minimum between the Eddington and the (classical) Bondi rate. Moreover, the SNIa rate is now computed using a high-precision multi-exponential approximation for the time-kernel in the convolution integral of the star formation rate, instead of the standard power-law. This avoids the need to store the entire star formation history and permits more rapid numerical simulations. The time-dependent mass return rate from the evolving stellar population is now computed using the Kroupa (2005) initial mass function and the mass return rate as a function of the stellar mass given by Maraston (2005), instead of the Ciotti et al. (1991) formulae. For the stellar mass-return rate, we also numerically implemented a multi-exponential fit. Finally, the dark-matter potential well (and the galaxy stellar distribution) are now described as Jaffe (1983) models.

From the astrophysical point of view, in addition to the standard outputs considered in SOCS (e.g. final SMBH mass, final stellar mass, etc), we now also focus on the duration of the “cold phase”, of the “obscured phase”, and the “Compton thick phase” as fiducially computed using the gas temperature and the column density of cold gas. The main parameters adopted to fix a model are the dark-matter halo mass (distributed to reproduce cosmological expectations), the amount of gas deemed to flow onto the dark matter halo, and finally the radiative accretion efficiency. As a separate parameter, we also consider the initial mass of the central SMBH. After some preliminary model exploration and verifying that the new treatment, when applied to the SOCS galaxy model reproduces the previous results, we performed a significant exploration of the parameter space.

Our main results can be summarized as follows:

• 1.

  Almost independently of the dark matter halo mass and theradiative efficiency, the computed galaxy models have an initialphase that extends for some Gyrs, in which the gas temperature islow and the gas density is high. These galaxies would be definedobscured quasars, were we to measure their gas column densitywithin an aperture radius of the order of R_e/100.

• 2.

  At late times, all the SMBHs are found to be in a low accretion state without much feedback, and the ISM is overall optically thin. The galactic ISM is at about the virial temperature of the dark matter potential well and should be emitting in X-rays, as found for elliptical galaxies in the local universe.

• 3.

  Interestingly, we have found that only a specific class of models, are found to agree with the observed Magorrian relation, i.e. only those with low initial SMBH masses (≲10⁶   M_⊙) and high radiative accretion efficiency, ϵ ~ 0.1. Higher initial SMBH masses, or lower radiative efficiencies lead to final SMBH masses that are too high. Therefore, this result implies that the seed SMBHs should be quite small, but that their mass accretion should occurmainly with high radiative efficiency, in agreement with observational findings.

• 4.

  In addition, we have also showed how the self-regulated Bondi accretion recipe can be easily implemented in numerical codes, and how it leads to a lower SMBH mass accretion than the more common “on-off” Eddington regulation.

A series of final comments are in order. The first is that the presentedone-zone models seem to be able to reproduce, without much fine tuning, two different phases of galaxy evolution, namely a first obscured phase where much of the SMBH accretion and star formation occurs, followed by a hot phase in which SMBH accretion is highly sub-Eddington and star formation is (almost) entirely absent. The robustness of the two-phase evolution, characterized both by one-zone and hydrodynamical simulations is mainly due to the combined effect of (1) the secular decrease in the mass return rate from the evolving stellar populations; (2) the time dependence of the SNIa heating (after the first Gyrs of evolution); and (3) the cooling function, which are identical in the hydrodynamical and the one-zone models. The increase in the specific ISM heating with increasing time, and that a substantial degassing occurs only when the gas temperature is roughly higher by a factor of  ~4 than the virial temperature, leads to the appearance of the two phases in the two types of model.

The second is that to reproduce the present-day Magorrian relation, the optimal combination of parameters is a quite low initial SMBH mass and a quite high radiative efficiency. We propose that this represents a useful constrain of semi-analytical investigations (also in the context of high-z galaxy merging). As a side-product of this work, we also anticipate that the presented multi-exponential fit for stellar evolution and SNIa heating will be useful in both hydrodynamical and semi-analytical works.

As a final comment, we point out a major improvement in the code. In a future study, we will relax the assumption of a fixed gas density distribution, imposing a time dependence of r_g(t) as a function of the gas thermal content at that given time. We expect that this additional ingredient will cause the models to become more sensitive to the adopted parameters on the one hand by increasing the gas density during the cold phases (hence producing a stronger feedback by decreasing the value of the ionization parameter and increasing the mass accretion rate), and on the other hand by decreasing the gas density during the hot phase, so reducing further the Eddington ratios at late times. We naturally, also expect that both the star-formation and accretion histories, as well as the durations of the obscured and Compton-thick phases to be affected.

---

Acknowledgments

The anonymous referee is acknowledged for his/her comments that greatly improved the paper. We thank A. Comastri and J. P. Ostriker for useful discussions, and F. Marinacci for suggestions on the numerical code. The project is supported by ASI-INAF under grants I/023/05/00.

References

Appendix A: Interpolating functions for stellar evolution

A.1. Stellar mass losses

The general expression for the stellar mass-return rate from an evolving stellar population is $\mathit{Ṁ}\begin{array}{}\mathrm{w}\\ \mathrm{\ast }\end{array}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{W}}_{\mathrm{\ast }}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}\mathrm{)}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathit{,}$(A.1)where $\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}$ is the instantaneous star-formation rate and W_∗(t − t′) is the normalized stellar death rate for a stellar population of age t − t′. For an initial mass function (IMF) of unitary total mass and age t, ${\mathit{W}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{IMF}\mathrm{\left[}{\mathit{M}}_{\mathrm{TO}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right]}\mathrm{\left|}\mathit{Ṁ}\mathrm{TO}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right|}\mathrm{\Delta }{\mathit{M}}^{\mathrm{w}}\mathrm{\left[}{\mathit{M}}_{\mathrm{TO}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right]}\mathit{,}$(A.2)where M_TO(t) is the mass of stars entering the turn-off at time t, and ΔM^w their mass loss. Following Ciotti & Ostriker (2007), we assume that $\mathrm{\Delta }{\mathit{M}}^{\mathrm{w}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{0.923}\mathit{M}\mathrm{-}\mathrm{0.48}\mathit{,} \mathrm{0.08}⩽\mathit{M}⩽\mathrm{8.5}\\ \mathit{M}\mathrm{-}\mathrm{1.4}\mathit{,} \hspace{0.17em}\mathrm{8.5}\mathit{<}\mathit{M}⩽\mathrm{40}\\ \mathit{M}\mathit{/}\mathrm{2}\mathit{,} \mathit{M}\mathit{>}\mathrm{40.}\end{array}$(A.3)In our models, we adopt a Kroupa (2001) IMF, with a minimum mass of 0.08    M_⊙ and a maximum mass of 100    M_⊙, while Ṁ_TO(t) is taken from Maraston (2005) $\begin{array}{c}\\ \log {\mathit{M}}_{\mathrm{TO}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\cong \mathrm{2.982}\mathrm{+}\mathrm{0.213}\log \mathit{t}\mathrm{-}\mathrm{0.108}\mathrm{(}\log \mathit{t}{\mathrm{)}}^{\mathrm{2}}\\ \hspace{0.17em}\mathrm{+}\mathrm{0.006}\mathrm{(}\log \mathit{t}{\mathrm{)}}^{\mathrm{3}}\mathit{.}\end{array}$(A.4)In the formula above, t is in yrs and M_TO in solar masses. The main computational problem posed by the evaluation of the integral in Eq. (A.1) is that, in principle, the entire history of $\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}$ must be stored, which requires a prohibitively large amount of memory and computational time. In the special case of an exponential time dependence of W_∗, this problem can be fortunately avoided, as shown in Ciotti & Ostriker (2001, 2007). In the present case, W_∗ is not an exponential function, being more similar to a power law. However, we assume that a fit in terms of exponentials has been obtained, i.e. ${\mathit{W}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}}\frac{{\mathrm{e}}^{\mathrm{-}\frac{\mathit{t}}{{\mathit{b}}_{\mathit{i}}}}}{{\mathit{a}}_{\mathit{i}}}\mathit{,}$(A.5)where the timescales a_i and b_i are known. From the above equation and Eq. (A.1), it follows that $\mathit{Ṁ}\begin{array}{}\mathrm{w}\\ \mathrm{\ast }\end{array}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}}\frac{{\mathit{F}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{a}}_{\mathit{i}}}\mathit{,}$(A.6)where ${\mathit{F}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}}{{\mathit{b}}_{\mathit{i}}}}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathit{.}$(A.7)We divide the integral into ${\mathit{F}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{t}\mathrm{-}\mathit{\delta t}}^{\mathit{t}}\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}}{{\mathit{b}}_{\mathit{i}}}}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}\mathrm{-}\mathit{\delta t}}\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}}{{\mathit{b}}_{\mathit{i}}}}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathit{,}$(A.8)where δt represents the last time-step. The evaluation of the first integral can be done by using a simple trapezoidal rule and only the values of $\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}$ at t and t − δt are required. The second integral is transformed by adding and subtracting δt in the exponential term, and simple algebra finally shows that ${\mathit{F}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathit{\delta t}}{\mathrm{2}}\left[\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{\delta t}}{{\mathit{b}}_{\mathit{i}}}}\mathit{Ṁ}\begin{array}{}\mathrm{+}\\ \mathrm{\ast }\end{array}\mathrm{(}\mathit{t}\mathrm{-}\mathit{\delta t}\mathrm{)}\right]\mathrm{+}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{\delta t}}{{\mathit{b}}_{\mathit{i}}}}{\mathit{F}}_{\mathit{i}}\mathrm{(}\mathit{t}\mathrm{-}\mathit{\delta t}\mathrm{)}\mathit{.}$(A.9)Therefore, the introduction of the multi-exponential fit allow us to compute the instantaneous mass return rate at time t just by storing the values of the F_i functions at the previous time-step. We performed a non-linear fit of the tabulated values of Eq. (A.2) (obtained with the exact formulae). We found that an acceptable approximation over the entire Hubble time is obtained with the coefficients in the first row of Table A.1 (with a percentage error of at most of 10% and only at very late times, when the mass return rate of a stellar population is negligible).

A.2. SN Ia

In a given simple stellar population, the total (i.e. volume-integrated) energy release by SNIa is ${\mathit{L}}_{\mathrm{SN}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{E}}_{\mathrm{SN}}{\mathit{R}}_{\mathrm{SN}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(A.10)where E_SN = 10⁵¹ erg is the fiducial energy released by SNIa event and R_SN is the SNIa rate. From Eqs. (11), (12) in Ciotti & Ostriker (2007), the SNIa rate for a stellar population of mass M_∗ can be written as ${\mathit{R}}_{\mathrm{SN}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{-13}}{{\mathrm{{\rm Y}}}_{\mathrm{B}\hspace{0.17em}\mathrm{\odot }}}\frac{{\mathit{M}}_{\mathrm{\ast }}}{ {\mathit{M}}_{\mathrm{\odot }}}{\left(\frac{\mathit{t}}{{\mathit{t}}_{\mathrm{H}}}\right)}^{\mathrm{-}\mathit{s}} \mathrm{\left[}\mathrm{yr}{\mathrm{s}}^{-1}\mathrm{\right]}\mathit{,}$(A.11)where we adopt a stellar mass-to-light ratio Υ_{B    ⊙}  = 5, θ_SN = 1, h = 0.7, t_H = 13.7 Gyr, and s = 1.1. In our model with a star formation rate of Ṁ_∗, it follows that in each time-step an amount Ṁ_∗δt of stars is added to the galaxy body. Therefore, the instantaneous total rate of SNIa explosion at time t is ${\mathit{R}}_{\mathrm{SN}}^{\mathrm{T}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{-13}}{{\mathrm{{\rm Y}}}_{\mathrm{B}\hspace{0.17em}\mathrm{\odot }} {\mathit{M}}_{\mathrm{\odot }}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathit{Ṁ}\mathrm{\ast }\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}{\left(\frac{\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\prime }}}{{\mathit{t}}_{\mathrm{H}}}\right)}^{-1.1}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathit{.}$(A.12)We expanded the dimensionless power-law kernel in the integral above again using a (dimensionless) multi-exponential function, so that all the considerations in Sect. A.1 hold and in particular Eq. (A.9). The functions F_i now contain new values of b_i (in units of Gyr^-1). The parameters of the multi-exponential fit (with percentual errors  <3% over 13.7 Gyr) are listed in the second row of Table A.1. Therefore, ${\mathit{L}}_{\mathrm{SN}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{38}}}{{\mathrm{{\rm Y}}}_{\mathrm{B}\hspace{0.17em}\mathrm{\odot }} {\mathit{M}}_{\mathrm{\odot }}}\sum _{\mathit{i}}\frac{{\mathit{F}}_{\mathit{i}}}{{\mathit{a}}_{\mathit{i}}} \mathrm{erg}{\mathrm{s}}^{-1}$(A.13)and finally, the average SNIa heating per unit time and unit volume, needed in the gas energy equation is $\mathit{Ė}\mathrm{H}\mathit{,}\mathrm{SNIa}\mathrm{=}\frac{\mathrm{3}{\mathit{L}}_{\mathrm{SN}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{8}\mathit{\pi }{\mathit{r}}_{\mathrm{g}}^{\mathrm{3}}}\mathit{,}$(A.14)where again we use half-mass values.

A.3. Column density for the Jaffe model

In our model the gas density at each computation time is assumed to be proportional to ρ_h and given by Eq. (19). The corresponding projected gas distribution is $\begin{array}{c}\\ {\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{\equiv }\mathrm{2}{\mathrm{\int }}_{\mathit{R}}^{\mathrm{\infty }}\frac{{\mathit{\rho }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{r}\mathrm{d}\mathit{r}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}}}\mathrm{·}\end{array}$(A.15)We note that the projected density formally diverges for R → 0, so that, to estimate a realistic value of the column density we need to compute a mean value of the projected density within a fiducial aperture. The explicit formulae for Σ_g(R), and the associated projected mass encircled with R, M_Pg(R), are given in Jaffe (1983). In particular, to estimate  ⟨ N_H ⟩  in Eq. (22) we consider M_Pg(R) $\begin{array}{c}\\ {\mathit{M}}_{\mathrm{Pg}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}{\mathit{R}}^{\mathrm{\prime }}\hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}\\ \mathrm{=}{\mathit{M}}_{\mathrm{g}}\mathrm{-}\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathit{R}}^{\mathrm{\infty }}\mathit{r}{\mathit{\rho }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}}\mathrm{d}\mathit{r}\mathit{.}\end{array}$(A.16)Integration of Eq. (A.16) yields $\begin{array}{c}\end{array}\frac{{\mathit{M}}_{\mathrm{Pg}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{{\mathit{M}}_{\mathrm{g}}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\{\begin{array}{c}\\ \frac{\mathit{\pi }}{\mathrm{2}}\mathrm{-}\frac{\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array} \mathrm{arcsech}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}}{\sqrt{\mathrm{1}\mathrm{-}\begin{array}{}\mathrm{˜}\\ {\mathit{R}}^{\mathrm{2}}\end{array}}}\mathit{,}\mathrm{0}\mathit{<}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{<}\mathrm{1}\mathrm{;}\\ \frac{\mathit{\pi }}{\mathrm{2}}\mathrm{-}\frac{\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array} \mathrm{arcsec}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}}{\sqrt{\begin{array}{}\mathrm{˜}\\ {\mathit{R}}^{\mathrm{2}}\end{array}\mathrm{-}\mathrm{1}}}\mathit{,}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{>}\mathrm{1}\mathrm{;}\end{array}$(A.17)where $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathrm{\equiv }\mathit{R}\mathit{/}{\mathit{r}}_{\mathrm{g}}$ and M_Pg(r_g) = (π/2−1)M_g.

All Tables

All Figures