LEAVING THE DARK AGES WITH AMIGA

Galaxies develop within dark matter (DM) halos through mergers and gas accretion. This "hierarchical scenario" (Rees & Ostriker 1977; Silk 1977; White & Rees 1978; Blumenthal et al. 1984; White & Frenk 1991) explains the main observed galaxy properties. However, some aspects of the nearby universe resist being satisfactorily recovered (e.g., Benson 2010; Cattaneo et al. 2006), and the increasing amount of data at progressively higher redshifts, z's, is permanently challenging our ideas within this theoretical framework.

Great progress has been achieved in the last decades in this field thanks to the use of hydrodynamic simulations (e.g., Tissera et al. 1997; Steinmetz & Navarro 1999; Springel 2000, 2005; Nagamine et al. 2004; Schaye et al. 2010) and semi-analytic models (SAMs; Kauffmann et al. 1993, 1999; Cole et al. 1994; Somerville & Primack 1999; Cole et al. 2000; Hatton et al. 2003; Benson et al. 2003; Menci et al. 2005; Bower et al. 2006; Monaco et al. 2007; Ricciardelli & Franceschini 2010; Font et al. 2011). SAMs are more flexible and more easily inform on the main properties of objects. However, they have the reputation of describing the baryon physics by means of recipes that are too simple and include too many parameters. Simulations certainly provide more detailed information and are, in principle, based on first principles. However, they involve the same recipes and parameters as SAMs at subresolution scales.

However, all of these tools suffer for the same limitations: the huge amounts of memory and CPU time involved. This is annoying for two reasons. First, galaxy formation is a non-linear process where the feedback of luminous objects to the intergalactic medium (IGM) plays a central role (see, e.g., Manrique & Salvador-Solé 2014, hereafter MSS, and references therein). However, those limitations prevent us from self-consistently treating the coupling of these two baryonic phases. In particular, the ionizing background, with important consequences for dwarf galaxies (e.g., Benson et al. 2002; Hambrick et al. 2011; Mamon et al. 2011), must be treated in an ad hoc fashion. Second, galaxy properties depend on those of their earlier low-mass progenitors. However, the highest redshift z and the minimum halo mass M_H that can be reached in studies of nearby galaxies are about 7 and 10⁹ M_☉, respectively, in both SAMs (e.g., Bower et al. 2006) and simulations (e.g., Schaye et al. 2010).⁵ Of course, when studies focus either on small regions or high z's, the limits are less stringent, yet still too restrictive.

More importantly, the first generation stars formed by molecular cooling from the original pristine gas, the so-called Population III (Pop III) stars, are responsible for the initial metal enrichment and reionization of the IGM, as well as for the seeds of massive black holes (MBHs). Their local effects can be studied in detail by means of high-resolution hydrodynamic simulations (Wise & Abel 2007, 2008; Turk et al. 2009; Stacy et al. 2010; Kim et al. 2011; Prieto et al. 2011). However, the limited dynamic range of simulations prevents us from simultaneously analyzing their cosmological effects. The only attempt to date to account for the feedback of Pop III stars in hydrodynamic simulations comes from Ciardi et al. (2000); meanwhile, in the case of SAMs better adapted, in principle, to the study of galaxy formation on cosmological scales, there is only the work by Choudhury & Ferrara (2005). Unfortunately, in both studies, the baryon physics is dealt with by means of too simple analytic recipes and galaxies are not realistically modeled.

Last but not least, IGM is a composite (several chemical species), multiphase (singly and doubly ionized bubbles and subbubbles embedded in a neutral background), inhomogeneous (density and temperature fluctuations) environment, whose accurate analytic modeling is hard to achieve without important simplifying assumptions.

An improved analytic treatment of IGM has recently been developed by MSS. In the present Paper, we couple this treatment with AMIGA, an Analytic Model of Intergalactic-medium and GAlaxy evolution specifically devised to monitor those cosmic components since the dark ages. AMIGA includes molecular cooling and Pop III stars with the most dramatic feedback. To save memory, AMIGA does not rely on the construction of individual halo merger trees, but instead on the interpolation in grids of halo properties that are progressively built starting from well-known boundary conditions. In addition, it makes the causal connection of physical processes usually dealt with independently from each other. This reduces the number of free parameters and leads to a model which is internally more consistent than previous SAMs. Its application to the study of reionization is provided elsewhere (E. Salvador-Solé & A. Manrique 2015, in preparation). Here, we describe the general model with special emphasis on its novelties.

The outline of the paper is as follows. In Section 2, we describe the general procedure followed in AMIGA. Sections 3, 4, 5, 6, and 7 are, respectively, devoted to the modeling of DM, gas, stars, galaxies, and MBHs. In Section 8, we summarize the mass and metallicity evolution of the different baryonic phases and galactic components, and in Section 9, we describe how the final photometric properties of luminous objects are computed. Finally, in Section 10, we discuss the main achievements and some fundamental results of AMIGA.

The specific results shown throughout the paper correspond to plausible values of the AMIGA parameters in the concordant ΛCDM cosmology characterized by $\Omega _\Lambda =0.712$, Ω_m = 0.288, Ω_b = 0.0472, H₀ = 69.3 km s⁻¹ Mpc⁻¹, n_s = 0.97, and σ₈ = 0.830 (Hinshaw et al. 2013). Whenever possible, they are compared to observational data in order to assess the performance of the models. The reader is referred to E. Salvador-Solé & A. Manrique (2015, in preparation) for detailed information on the source of these data and the parameter values.

To minimize the memory and CPU time requirements of AMIGA, we pay special attention to treating every random process entering the problem in the best-suited way. If every single event has a noticeable, possibly critical effect, then the random process is dealt with in a fully probabilistic fashion. Otherwise, it is dealt with in a deterministic fashion by calculating its secular action in the desired time interval. In both cases, we use either analytic or well-sampled numerical probability distribution functions (PDFs).

AMIGA does not proceed by explicitly constructing Monte Carlo or N-body halo merger trees, but by interpolating the typical properties of halos in neutral and ionized regions in two parallel three-dimensional (3D) grids with n_z log-bins of 1 + z, n_H log-bins of DM masses M_H, and n_a linear bins of halo ages t_a ≡ t(z) − t_H, where t(z) is the cosmic time corresponding to z and t_H is the halo formation time, defined as the time of the last major merger having caused the rearrangement of the system. These grids are progressively built (see Figure 1), starting at a redshift z_max high enough for halos of all masses to have trivial properties down to the redshift z_obs of observations; at every z, from a value of M_H low enough (10 M_☉ at z_max = 60) for halos to also have trivial properties up to a value high enough (10¹⁵ M_☉ at z_obs = 0) for them to be highly improbable; and at every couple of z and M_h values, for a set of halo ages spanning the relevant time interval. In this way, integrating at every z over halo ages for the halo formation time PDF, and over halo masses for the appropriate halo mass function (MF), we determine the instantaneous change induced by luminous sources in the IGM properties at that z. This is a notable improvement compared to ordinary SAMs where the feedback of luminous sources at a given z is only known for the small number of halo masses and ages covered by the discrete branches of the merger tree that is being built.

Zoom In Zoom Out Reset image size

To obtain the typical properties of a halo at a new grid point (z, M_H, a), we chose the masses and formation times of its progenitors according to the corresponding PDFs, and we find their typical properties by interpolation in the pieces of grids previously built covering all possible halo ages, all halo masses less than M_H, and all redshifts higher than z. Based on the properties of the progenitors, we determine those of the halo at formation, and then follow its evolution through continuous accretion until reaching the z of the final object (see Figure 1). The accretion rate is given by the analytic halo growth model (Section 3), and the composition of the accreted matter at any moment is well known: (1) substantial halos whose properties are obtained by interpolation within the grids, (2) tiny halos with trivial properties lying outside the grids, and (3) a well-determined fraction of non-trapped intergalactic gas (Section 4.1). Finally, averaging the properties of the halo obtained from each progenitor configuration, we obtain the quantities to be stored in the new point of the grid.

All halo properties, including the baryonic content, stored at every point of the grids and the main physical processes where they are involved are represented in the notation used throughout the paper in Figure 2. Taking advantage of the fact that satellites are numerous, and hence can be dealt with statistically, we do not store the values of their individual properties as for the central galaxy, but we use their occupation numbers in a multi-dimensional space of galactic properties, with n_f linear bins of formation time, i.e., the last moment the satellite structure was rearranged, n_m log-bins in total mass, n_bm log-bins in baryonic to total mass ratio, 2 n_sb log-bins in disk and spheroid stellar to total baryonic mass ratio, 3 n_Z log-bins in disk and spheroid stellar metallicity and disk gas metallicity, and n_σ log-bins in disk central surface density. This latter property is used to calculate the disk scale radius given its mass, while the spheroid scale radius is calculated using the average dissipative contraction factor (see Section 6.1) of those central spheroids with identical stellar masses and formation times. Finally, the mass of satellite MBHs is calculated using the constant average MBH to stellar mass ratio of central spheroids with identical masses.

Zoom In Zoom Out Reset image size

AMIGA is implemented in an OMP (shared memory) parallelized code with 32 CPUs. The run time mostly depends on the size of the satellite array and the value of z_obs as the non-null occupation numbers filling that array increase with decreasing z. For z_obs = 2 and the minimum sizes of the interpolation grids (n_z = 51, n_H = 91, and n_a = 3) and of the satellite array (n_f = 8, n_m = 38, n_bm = 24, n_sb = 6, n_Z = 4, and n_σ = 4) ensuring convergence, it takes about 76 hr. For such standard dimensions, the typical properties of a halo in the grid arise from 3² distinct progenitor configurations, and the typical properties of galaxies in halos with M_H at z arise from to 2 × 3³ different halo progenitor configurations, with the additional factor of three arising from the different halo ages, and the factor two arising from the fact that the host halo may form either in a neutral region (before it harbors galaxies) or in an ionized one.

At the time of matter–radiation equality, DM begins to cluster in halos that merge with each other and grow from small to large scales. Halos will become the backbone of galaxies, and so it is mandatory to accurately model their mass growth, inner structure, and kinematics. This is achieved in AMIGA within the framework of the excursion set formalism, as in normal non-hybrid SAMs, but in a slightly modified version referred to as the modified extended Press–Schechter (MEPS) formalism (Salvador-Solé et al. 1998; Raig et al. 2001), with important advantages compared to the usual extended Press–Schechter (EPS) model.

The conditional MF in the EPS model diverges in the limit of small M_H. Thus, merger trees are infinitely ramified, which forces one to adopt a finite resolution in mass and time. The finite resolution in mass prevents us from properly dealing with the capture of low-mass halos contributing to accretion, and introduces some uncertainty in the total number and mass of the resolved progenitors (Somerville & Primack 1999). While the finite resolution in time prevents us from accurately dealing with mergers, since the conditional MF ensures only that a halo with M_H at t is found at t' > t within another halo with $M_{\rm H}^{\prime }> M_{\rm H}$, it does not inform us of the exact moment when the incorporation takes place. However, the properties of a halo at a node of the merger tree are inferred from those of its progenitors at the previous node, assuming the merger takes place at that moment, evolved until the time of the final node. In other words, the timing and properties of the evolving objects do not match those of the real merging process (see Figure 3). To minimize the effects of such inaccuracy, a relatively small time step must be adopted, but then the need of storing all the information on the merger tree prevents reaching a very high z.

Zoom In Zoom Out Reset image size

The MEPS formalism making the distinction between minor and major mergers does not have any of these drawbacks. Major mergers are really binary (Raig et al. 2001), so they can be dealt with statistically, in a fully accurate fashion with no need to introduce any limited resolution. Minor mergers can be dealt with fully accurately through their global secular action. Specifically, one can calculate the halo DM mass accretion rate, $\dot{M}_{\rm H}$, as a function of M_H and t. In addition, the MEPS formalism allows one to calculate the PDFs of halo formation times and progenitor masses (Raig et al. 2001), which is not available in the usual excursion set formalism. All of these differences lead to exact merger trees with realistic discrete branching (see Figure 3). Furthermore, the MEPS formalism also allows one to accurately derive the inner structure and kinematics of virialized DM halos (Salvador-Solé et al. 2007, 2012a, 2012b). We stress that, contrary to ordinary SAMs, all of these halo properties are determined in AMIGA directly from the cosmology considered without a single free parameter.

AMIGA monitors the evolution of halos lying in neutral and ionized regions, separately. Such a distinction is important because pristine gas only falls inside halos lying in neutral regions; in ionized regions, the IGM is polluted with metals produced in galaxies. The feedback of luminous sources to the ionized IGM is computed using the halo MF in those ionized regions, which is slightly different from that of halos in neutral regions (see MSS).

Until recombination at z ∼ 1100, radiation pressure prevents the ionized gas from falling into the halo potential wells. Nonetheless, until z = 150(Ω_bh²/0.023)^2/5 − 1 ≃ 150, the abundance of free electrons is high enough for the neutral gas to be kept thermalized with the cosmic background radiation (CMB). At that z, the residual abundance of free electrons ($\bar{x}_{\rm e}\approx 3.1\times 10^{-4}$) freezes out and the gas begins to undergo adiabatic cooling. At the beginning, the gas is too hot to be trapped by the only mini-halos significantly abundant at those z. Only after z ∼ 50 is the gas cold enough to fall into the potential wells of reasonably abundant halos with M_H ∼ 10⁵ M_☉, giving rise to the formation of the first generation stars.

4.1. Unbound IGM

Luminous sources reionize and reheat the diffuse unbound IGM, hereafter simply the IGM. UV photons with short mean free paths ionize the bubbles around them, which then grow and progressively overlap. Inside these bubbles, subbubbles with doubly ionized helium develop due to the smaller fraction of more energetic UV photons. X-ray photons with a much larger mean free path give rise instead to a uniform background also heating the IGM through Compton scattering (and through secondary ionizations, neglected in AMIGA).

Some amount of the diffuse gas in the IGM is accreted by sufficiently massive halos (inflows) or is expelled from them (outflows). The gas mass inflow rate, $\dot{M}_{\rm hg}^{\rm in}$, is proportional to $\dot{M}_{\rm H}$ (Section 3) with a proportionality factor equal to the current baryon mass fraction in the IGM, calculated from the original total baryon fraction taking into account the gas gains and losses into and from halos. Gas outflows are triggered by supernova (SN) or active galactic nucleus (AGN) driven winds (Sections 5.3.3 and 7.2, respectively) when they cause the hot gas in the halo to become unbound. Its typical rate, $\dot{M}_{\rm hg}^{\rm out}$, is taken to be equal to the hot gas mass lost over the wind duration. Outflows from halos also cause the chemical enrichment of the unbound IGM. As this effect takes place only in the vicinity of halos hosting luminous sources, AMIGA assumes that the metal pollution of IGM affects ionized bubbles only.

AMIGA deals with the properties of the neutral and ionized IGM separately, distinguishing between He ii and He iii ionized regions. The evolution with z of the IGM temperature, T_IGM, or, more exactly, the average temperature T_j of the gas in phases j = I, II, and III, for the neutral, singly, and doubly ionized regions, respectively, is according to the differential equation (MSS)

$$\begin{equation} \frac{d \ln T_{\rm j}}{d \ln (1+z)}= 2+\frac{d \ln (\mu _{\rm j}\varepsilon _{\rm j}/n_{\rm j})}{d \ln (1+z)}, \end{equation} \tag{ 1 }$$

where μ_j, ε_j, and n_j are the average mean molecular weight, energy density, and baryon density, respectively, in region j. In Equation (1), the term equal to 2 on the right gives the cosmological adiabatic cooling, and the second term includes Compton heating–cooling from the CMB and X-rays, heating–cooling through ionization–recombination of the various chemical species, collisional cooling of hot neutral regions, the achievement of energy equipartition of newly ionized–recombined material, inflows–outflows from halos, and cooling by collisional ionization and excitation. The heating–cooling by gravitational compression–expansion of density fluctuations vanishes, neglecting non-linear effects, after averaging over each ionization phase (MSS).

Figure 4 illustrates the kind of IGM temperature evolution that can be obtained depending on the values of the free parameters of AMIGA, such as b_cl (see below for the rest). The specific solutions shown correspond to two distinct Pop III star initial mass functions (IMFs): the top-heavier one leads to double H i reionization, while the less top-heavy one leads to single reionization. As can be seen, in the case of double reionization, there is a marked dip in the temperature of the (singly and doubly) ionized IGM over the redshift interval between the two full ionizations (at z = 5.5 and 10.3), which is absent in the case of single reionization. This is due to the drop in flux of the ionizing photons at the first complete ionization (z = 10.3), when Pop III stars stop forming because of the lack of neutral regions. This causes a short period of H i recombination until the UV flux from normal galaxies becomes (at z = 8) high enough for reionization to start again.

Zoom In Zoom Out Reset image size

4.2. Trapped Hot Gas

4.2.1. Structure

The hot gas in equilibrium within halos is not assumed to be isothermal, as often done in SAMs, but with a polytropic equation of state with index Γ = 1.2, consistent with the MEPS formalism for halo growth. The gas that is accreted by a halo is shock-heated and deposited at the instantaneous virial radius of the halo, meaning that its spatial distribution grows from the inside out, similar to that of the DM (Salvador-Solé et al. 2012b), with the only difference being that the gas has a polytropic equation of state as a result of the shock, whereas DM follows a density profile à la NFW (Navarro et al. 1997) set by the rate at which halo accretes the non-collisional DM. It is the preservation of the ratio of total energies between the gaseous and dark components that fixes the value ∼1.2 of Γ (Solanes et al. 2005). Although this reasoning applies to halos grown by pure accretion, the same result holds for halos that have suffered major mergers (Salvador-Solé et al. 2012a, 2012b). Such a gaseous structure is not only expected on theoretical grounds, but it is also supported by observations (Ponman et al. 2003; Pratt et al. 2010) and simulations (Voit et al. 2005; Short et al. 2010). Furthermore, it leads to X-ray scaling relations from galaxy groups to rich clusters that are in very good agreement with observation (Solanes et al. 2005; Bode et al. 2009). Thus, the inner structure of the hot gas within halos is also calculated directly from the cosmology considered with no free parameters.

4.2.2. Cooling

The hot gas radiates, cools, and contracts in a runaway process that leads to the infall of cold gas to the halo center. The treatment of cooling carried out in AMIGA is the same as in conventional SAMs. The cooling radius, r_cool(t), encompassing the gas, having had time to cool since the formation of the halo, is found by equating the age of the halo, t_a, to the characteristic cooling time, given by the ratio between the energy density and the emissivity, ${\cal E}_{\rm hg}(r)/\dot{\cal E}_{\rm hg}(r)$, of the hot gas. The cooling rate is then

$$\begin{equation} \dot{M}_{\rm cool}= 4\pi \,r_{\rm cool}^2\,\mu m_{\rm p}n_{\rm hg}[r_{\rm cool}(t)]\,{d r_{\rm cool}\over d t}, \end{equation} \tag{ 2 }$$

where m_p is the proton mass, μ is the mean molecular weight of the hot gas, and n_hg is its particle number density.

Star formation begins to proceed at a significant rate in mini-halos with T_hg ∼ 10³ K. At these temperatures, the gas is essentially neutral and atomic cooling is not effective. The only way such a primordial gas can radiate is by spontaneous emission of roto-vibrational molecular levels excited by collisions of atoms with H₂ molecules (with a fraction as small as x ∼ 10⁻⁶). In more massive halos with virial temperatures above 10⁴ K, H₂ molecules are dissociated by collisions with atoms, and the spontaneous emission of atomic electronic levels excited by collisions of atoms with free electrons is the dominant cooling mechanism. Then, the higher the metal abundance, the more effective cooling is. Four possible cases are found.

• 1.  
  If the metallicity is higher than a critical value Z_crit and the gas is ionized (T_hg > 10⁴ K), then atomic line cooling takes place. The emission through metallic lines continues to operate when the gas cools below 10⁴ K or when the density increases to the point that the gas becomes shielded to ionizing photons and recombines. Molecular and dust emission then become active (see next item), but the total cooling rate is limited by the initial atomic rate.
• 2.  
  If the metallicity is higher than Z_crit and the gas is not ionized (no stars in the halo), then molecular cooling can proceed by means of H₂ and many other molecules synthesized on dust grains as well as by the emission of dust itself. The rate of this complex cooling process is unknown and, contrary to case (1), it is not limited by the rate of any previous cooling mechanism. However, this case can be neglected because the high metallicity is indicative that star formation has taken place in the halo progenitors.
• 3.  
  If the metallicity is below Z_crit and T_hg > 10⁴ K, then a sequence of two different processes takes place. Since at such temperatures H₂ molecules are dissociated due to collisions with H atoms, the first mechanism to operate is atomic cooling. However, once T_hg drops below 10⁴ K, most H ii recombines and, as there are no free electrons, atomic cooling halts. Then, the gas switches to H₂ cooling for a concentration of this molecule corresponding to equilibrium (Oh & Haiman 2002).
• 4.  
  If the metallicity is lower than Z_crit but T_hg < 10⁴ K, then atomic cooling is not efficient because there are essentially no free electrons that can excite atoms (not only because the gas is mostly recombined, but also because the remaining electrons do not have enough energy to excite hydrogen atoms at those temperatures, the lowest excitation level requiring an energy of 1.2 × 10⁵ K), and the gas cools directly by H₂-molecular emission.

The latter two cases presume of course that the gas is not ionized by luminous sources. Otherwise, molecules could not form and molecular cooling would not be effective. The critical metallicity, Z_crit, is taken to be equal to 10⁻⁴ Z_☉ (Santoro & Shull 2006; Smith et al. 2009; Schneider & Omukai 2010).

In case (1), the emissivity leading to Equation (2) is given by the usual expression for atomic cooling:

$$\begin{equation} \dot{\cal E}_{\rm hg}(r)=n_{\rm hg}^2(r)\Lambda [T_{\rm hg}(r),Z_{\rm hg}]\,, \end{equation} \tag{ 3 }$$

with the cooling function Λ[T_hg(r), Z_hg] drawn from Sutherland & Dopita (1993).

In cases (3) and (4), the cooling rate depends on the abundance of H₂ molecules, and the situation is more complex due to the strong feedback that stars have on the H₂ concentration. In metal-free gases, this molecule forms through reactions catalyzed either by electrons or by protons. See Table 1 for the different possible reactions: the second and third channels for electrons, and the fourth and fifth ones for protons. The former of these channels is the most efficient, although for completeness both channels are included in AMIGA. The first reaction corresponds to recombination. We will come back to recombination later on.

Table 1. Rates for the Reactions Involved in H₂ Formation

Reaction	Rate Coefficient
	(cm3 s−1)
H+ + e− → H + hν	k1 = 8.4 × $10^{-11}T_3^{0.2}/ \left(T^{0.5}\left[ 1+T_6^{0.7} \right]\right)$
H + e− → H− + hν	k2 = 1.4 × $10^{-18}T^{0.928} e^{-T/[1.62\times 10^4]}$
H + H− → H2 + e−	k3 = 1.3 × 10−9
H + H+ → H$_2^+ + h\nu$	k4 = 2.1 × 10−23T1.8e−20/T
H + H$_2^+ \rightarrow$ H2 + H+	k5 = 6.4 × 10−10

Note. T is the gas temperature in K and T_n = T/10ⁿ.

Download table as:  ASCIITypeset image

The concentration of H₂ in the gas of a newly born halo is computed in AMIGA according to the reactions and corresponding rates given in the second column of Table 1 (taken from Hutchings et al. 2002) for the appropriate density and temperature of the gas in the halo and from the initial concentrations and total abundances of all of the initial chemical species, namely, H, H⁺, H⁻, H₂, H$_2^+$, and e⁻, previously calculated for each progenitor, starting from the trivial initial concentrations outside the interpolation grid of AMIGA given by Galli & Palla (1998).

Provided that there is no star in the halo (otherwise molecular cooling is inhibited), the gas cools efficiently and contracts until n_hg and T_hg reach some critical values n_c and T_c, respectively, equal to 10⁴ cm⁻³ and 100 K. Then cooling halts. The H₂ emissivity leading to such a stable state is

$$\begin{equation} \dot{\cal E}_{\rm hg}(r)=f_{\rm H_2}(r)\,n_{\rm hg}^2(r)\,\Lambda _{\rm H_2}[T_{\rm hg}(r)]\,, \end{equation} \tag{ 4 }$$

where $f_{\rm H_2}$ is the number fraction of H₂ molecules, and $\Lambda _{\rm H_2}(T_{\rm hg})$ is the associated cooling function given by Galli & Palla (1998) for the H₂ concentration at equilibrium (case (3)) or calculated in the manner explained above (case (4)). After reaching the minimum temperature T_c, the cold gas accumulates at the halo center until one Bonner-Ebert, or simply one Jeans mass,

$$\begin{equation} M_{\rm J}=\left(\frac{\gamma \pi k_{\rm B}T_{\rm c}}{G \mu m_{\rm p}} \right)^{3/2}\rho _{\rm c}^{-1/2}\,, \end{equation} \tag{ 5 }$$

is reached. In Equation (5), γ is the adiabatic index, k_B is the Boltzmann constant, and ρ_c is the mean inner density of the isothermal sphere with temperature T_c or, more exactly, with temperature equal to the maximum between T_c and the CMB temperature at that z, as the cold gas is heated by the background radiation. Then it fragments and collapses to form a small cluster of metal-free stars of about 1000 M_☉.⁶

All physical processes calculated so far, directly from the cosmology considered with no free parameters, are consistent with both observations and simulations.

5.1. Star Formation

Due to the presence of metals, case (1) leads to the formation of ordinary Population I and II (Pop I and II) stars, whereas cases (3) and (4) lead to Pop III stars.

In the metal-rich case, the cold gas collected at the halo center tends to settle down in a centrifugally supported disk or directly feeds a central spheroid (see Section 6), where it gives rise to star formation.

The star formation rate (SFR), $\dot{M}_{\rm sf, C}$, in the galactic spheroid (C=B) or disk (C=D) is taken according to the usual Schmidt–Kennicut law (Kennicutt 1998),

$$\begin{equation} \dot{M}_{\rm sf, C}=\alpha _{\rm G}\frac{M_{\rm cg,C}}{{\tau _{\rm dyn,C}}}, \end{equation} \tag{ 6 }$$

where M_{cg, C} is the mass of cold gas available, τ_{dyn, C} is the dynamical timescale at the half-mass radius of the component C (see Section 6.1), and α_G is the star formation efficiency, taken as a free parameter.

Pop I and II stars form according to the IMF, ϕ(m), along the zero-age main sequence and evolve along the respective mass- (and metallicity-) dependent evolutionary tracks.⁷ AMIGA is ready for any wanted IMF, but the default one is the modified Salpeter IMF proposed by Wilkins et al. (2008), consistent with observations of the local IMF, which is characterized by a power-law index equal to −1 in the range 0.1 < m/M_☉ <0.5 and equal to the Salpeter value (−2.35) in the range 0.5 < m/M_☉ <100. This does not preclude, of course, that the IMF of metal-poor stars is much more top-heavy (i.e., with a greater lower mass or a less steep logarithmic slope).

Pop III stars are believed, indeed, to reach masses well above 100 M_☉. Those with m ≲ 130 M_☉ would explode as type II SNe and produce metals typically according to the yield p of ordinary Pop I and II stars (see Section 5.3.4), whereas those with 130 ≲ m/M_☉≲ 260 would explode as pair instability SNe (PISNe) and release about half their mass as metals. Finally, those with m ≳ 260M_☉ would collapse into a black hole (BH) and leave no yield at all (Heger & Woosley 2002). Thus, denoting the mass fractions of metal-free stars in these three mass ranges of increasing mass as f₁, f₂, and f₃, we see that the parameters f₂ and f₃ (f₁ = 1 − f₂ − f₃) are enough to characterize the Pop III star IMF and evolution.

5.2. Stellar SEDs

To calculate the emission from normal galaxies or low-mass (m < 130 M_☉) Pop III stars, an SED is assigned to each group of stars in the color–magnitude diagram according to its time-varying spectral type and luminosity class (see Section 9.1 for details). The spectrum of massive (m > 130 M_☉) Pop III stars is approximated by that of a black-body with effective temperature T_eff equal to ∼10⁵ K regardless of their mass (Bromm et al. 2001). The superposition of all these spectra gives rise to the synthetic SED of the whole stellar population of the galactic component under consideration as a function of time.

The contribution F_λ(t) to the galactic flux at wavelength λ at time t of a given stellar population is

$$\begin{equation} F_\lambda (t)= \int ^t_0 d t^{\prime } \int _{m=0}^{\infty } d m \, \dot{M}_{\rm sf}(t-t^{\prime }) \,\phi (m)\, f_\lambda (m,t^{\prime },Z_{\rm s}), \end{equation} \tag{ 7 }$$

where f_λ(m, t', Z_s) is the flux at wavelength λ of one individual star with initial mass m, initial metallicity Z_s, and age t' = t − t_s, t_s being its formation time, with an origin at the zero-age main sequence (i.e., f_λ(m, t', Z_s) = 0 for t' greater than the lifetime of the star), provided by the adopted stellar population synthesis model (Section 9.1).

5.3. Stellar Feedback

Stars affect the surrounding IGM in three different ways: by increasing its metallicity through SNe and stellar mass losses, by reheating it mechanically through SN shocks, and through radiative losses. Since being richer in metals equates to the hot gas cooling more easily, metal enrichment is a positive feedback for star formation. On the contrary, reheating by SNe may cause part of the metal-enriched interstellar medium (ISM) to escape from the halo (outflows) and, consequently, it is a negative feedback like the photo-dissociation of molecules and the photo-ionization and reheating of the IGM. In all of these feedback processes, there are important differences between Pop I and II and Pop III stars.

5.3.1. Photo-dissociation

H₂ is dissociated by photons with energy below the Lyman limit, in the so-called Lyman–Werner (11.28–13.6 eV) bands (Haiman et al. 2000). This effect operates when SFR becomes intense enough for a cosmic soft UV background to build up.

The rate at which dissociation takes place can be approximated by (Abel et al. 1997)

$$\begin{equation} k_{\rm diss}=1.38\times 10^{-12}J_{21}(h\nu =12.87\, {\rm eV}) \end{equation} \tag{ 8 }$$

in s⁻¹ with J₂₁ being the flux in units of 10⁻²¹ erg $\rm {s}^{-1}\ \rm {cm}^{-2}$ $\rm {Hz}^{-1}\ \rm {str}^{-1}$. This flux should be essentially homogeneous and isotropic since the distance traveled by the photons is far larger than the mean separation between halos (absorption by the medium can be neglected). In these conditions, redshifting of the photons must be taken into account. Thus, by integrating the flux of dissociating photons from all Pop III stars in a given volume, one can obtain the emissivity of dissociating photons at each z, j_ν(z), and from this the corresponding flux

$$\begin{equation} J_\nu (z)=\int _z^{z_{\rm max}}dz^{\prime }c\,\frac{d t}{dz^{\prime }}j_{\nu ^{\prime }}(z^{\prime }), \end{equation} \tag{ 9 }$$

where c is the speed of light and z_max = (13.6/11.28)(1 + z) − 1 is the redshift at which photons with an energy of 13.6 eV are redshifted to 11.28 eV. Note that the upper bound in the integral of Equation (9) coming from the fact that any UV photon emitted at a redshift z_em > z_max falling into the soft UV bands after redshifting to z will have been previously absorbed by the neutral IGM. Actually, the flux given by Equation (9) is shielded inside halos due to the molecules produced since virialization. Thus, to calculate the photo-dissociation rate, the emissivity given in Equation (9) must be multiplied by the shielding factor F_sh, estimated through (Draine & Bertoldi 1996)

$$\begin{equation} F_{\rm sh}={\rm min}\left[1,\left(\frac{N_{\rm H_2}}{10^{14}\,\rm {cm}^{-2}}\right)^{-3/4}\right], \end{equation} \tag{ 10 }$$

where $N_{\rm H_2}$ is the column H₂ density.

5.3.2. Photo-ionization

To calculate the flux of H i/He i- and He ii-ionizing photons emitted by Pop III and Pop I and II stars (as well as AGNs), one must consider the different SEDs of the emitting objects. The SED of metal-rich and low-mass metal-poor stars is provided by the stellar population synthesis model, taking into account the star formation and metallicity histories of the emitting populations. The flux of ionizing photons from zero-metallicity Pop III stars and the nebular emission they induce is computed according to Schaerer (2002) for the particular IMF (i.e., the mass fractions f₁, f₂, and f₃) considered.

At each z, we compute the flux of ionizing photons escaping from galaxies in halos with different masses along the whole (relevant) range, and integrate for the halo MF corresponding to that z in the ionized environment. For halos with virial temperatures lower than 10⁴ K, the escape fraction of photons above the Lyman continuum limit, Ly_c, is obtained by subtracting those photons captured by the neutral gas present in it. While if the temperature is higher than 10⁴ K, we assume some escape fraction, f_esc. The possibility that f_esc increases with increasing z has been considered in order to obtain a reionization at high z, as suggested by the analysis of CMB anisotropies (Pawlik et al. 2009; Kuhlen & Faucher-Giguère 2012; Alvarez et al. 2012). However, such behavior is not supported by observations, and so in AMIGA we adopt the usual, more conservative assumption of a constant f_esc taken as a free parameter. Recombination, both inside and outside halos where ionizing photons are produced, is also taken into account as it leads to the absorption of more ionizing photons.

The evolving H ii and He iii volume filling factors, Q_HII and Q_He iii, are governed by the differential equations for trivial initial conditions at the dark ages (MSS):

$$\begin{equation} \frac{d Q_{{\rm S}\, \scriptsize{II}}}{d t}\!=\!\frac{\langle {\dot{N}}_{{\rm S}\, \scriptsize{II}}\rangle }{\langle n_{\rm S}\rangle } \!-\!\left[\!\left\langle \! \frac{\alpha _{{\rm S}\, \scriptsize{I}}(T_{\rm IGM})}{\mu ^{\rm e}}\!\right\rangle _{{\rm \!\!S}\, \scriptsize{II}}\!\!\frac{C\,\langle n_{\rm b}\rangle }{a^3(t)}\!+\!\frac{d \ln \langle n_{\rm S}\rangle }{{\rm d}t}\right]\!Q_{{\rm S}\, \scriptsize{II}}, \end{equation} \tag{ 11 }$$

where the subscripts S, SI, and SII stand for H, H i, and H ii, or He i, He ii, and He iii, and angular brackets indicate averages over the regions denoted by a subscript (lacking any subscript, the average is over the whole IGM). The average of a function f(T_IGM) of the IGM temperature in the region j is taken to be equal to $f(T_{\rm j})+({\rm d}^2 f/{\rm d}T^2){{\sigma _{\rm T}^2}_{\rm j}}/2$, with T_j equal to the mean temperature in that region and ${\sigma _{\rm T}^2}_{\rm j}$ the corresponding variance. In Equation (11), 〈n_b〉 is the comoving cosmic baryon density, a(t) is the cosmic scale factor, μ^e is the electronic contribution to the mean molecular weight, $\dot{N}_{\rm SII}$ is the comoving metagalactic ionizing photon rate density due to luminous sources and recombinations (calculated according to Meiksin 2009) to He ii and He i ground states for H i-ionizing photons (for simplicity, the contribution from He ii Lyα recombinations is neglected), α_SI is the recombination coefficient to the SI species, and C is the ionized clumping factor.

Cosmological N-body simulations allow one to estimate the clumping factor from the observed fluctuations in DM, C_ρ, for which we have practical analytic fits (Iliev et al. 2007; Raicević & Theuns 2011). If baryons traced mass, then C_ρ and C would be identical. Unfortunately, the limits of the diffuse IGM are hard to established in terms of the DM density field. On the other hand, the increased pressure in the ionized gas may largely reduce its density fluctuations (Miralda-Escudé et al. 2000; Pawlik et al. 2009; Finlator et al. 2012). For this reason, we adopt the relation C = b_clC_ρ with the matter clumping factor C_ρ drawn from simulations and the bias factor b_cl taken as a free parameter.

5.3.3. SN Reheating

The secular effects of reheating on the unbound IGM are accounted for through Equations (1), while the extra energy of the non-gravitational origin imparted to the hot gas trapped in halos is transferred jointly with the gas itself (and galaxies) to the descendants of every halo.

X-ray photons produced in SNe (from free-free emission and inverse Compton scattering of CMB photons by relativistic electrons) and, to a lesser extent, emitted from very massive Pop III stars, ordinary binary stars, and AGNs Compton heat the IGM and increase the entropy-floor of the non-trapped gas. The fraction of the SN energy converted to X-rays is about 1% (Oh & Haiman 2003).

SNe also mechanically reheat the ISM in the disks and spheroids of normal galaxies as well as the hot gas in halos with metal-poor galaxies. For a given stellar mass at formation, M_sf, some fraction is in massive fast-evolving stars that quickly explode as SNe. When an SN occurs, some amount of the energy released is imparted to the surrounding ISM, causing part of it to blow off to the halo, in the case of ordinary stars, or to be directly imparted to the hot gas in the halo, which can be ejected from it (usually it is), in the case of very massive metal-free stars.

The condition that the reheated ISM leaves the component C of galaxies and joins the hot gas in the halo leads to the usual expression,

$$\begin{equation} M_{\rm rh,C}=\epsilon _{\rm C}\frac{2\,\eta _{\rm SN}E_{\rm SN}}{V_{\rm C}^2-V_{\rm hg}^2} M_{\rm sf,C}\, \end{equation} \tag{ 12 }$$

for the mass of reheated gas, with its time-derivative leading to the rate, $\dot{M}_{\rm rh,C}$, at which the ISM of component C is reheated and expelled from the galaxy. In Equation (12), E_SN ≈ 10⁵¹ erg is the energy produced in one Type II SN, η_SN = 0.0144 is the number of SN explosions per solar mass unit over the typical duration of a starburst (0.2 Gyr) of stars formed instantaneously with the (modified) Salpeter IMF, V_hg is the thermal velocity of the hot gas at the halo half-mass radius, V_C is the circular velocity at the radius used to set the typical dynamical time of component C, and _C is the corresponding SN reheating efficiency.

Hydrodynamic studies indicate that reheating of ISM by type II SNe triggers galactic winds only in spheroids, which is also in agreement with observation. The reason for this would be that in disks, only a small fraction of the SN energy is directed toward the plane where the ISM lies, which greatly diminishes the heating efficiency. For this reason, we take _C to be equal to one for C=B (e.g., Dekel & Silk 1986; Mac Low & Ferrara 1999; Strickland & Stevens 2000) and equal to zero for C=D. However, when the gas has little angular momentum and the stable disk is found to be smaller than the corresponding spheroid (bulge), we assume that it settles in an oblate pseudo-bulge rather than in a thin disk, and so _C is then also taken to be equal to one. Of course, the effective amount of reheated gas leaving a galactic component depends not only on the reheating efficiency, but also on the gravitational pull of the galaxy, accounted for through the circular velocity V_C.

In the case of Pop III star clusters, Equation (12) also holds, but with V_hg = 0 and the circular velocity V_C at the half-mass radius of the halo, as the reheated gas is then expelled from it out to the unbound IGM. On the other hand, the energy E_SN liberated by one PISN explosion is two orders of magnitude larger than in Type II SNe (Fryer et al. 2001; Heger & Woosley 2002), and the expected number η_SN of SN explosions per solar mass unit of stars formed in a typical starburst is 0.0015 (Schaerer 2002). The reheating efficiency of PISNe explosions is likely also different from that of normal SNe, but the exact value does not matter provided that it is large enough for mini-halos to lose the hot gas in those explosions. This is indeed what happens for equal to unity, as also adopted for simplicity in AMIGA.

The stellar mass loss, M_{loss, C}, going into the ISM of component C from stars with masses from m₁ to m₂ in a stellar population with total mass M_{sf, C} at formation is

$$\begin{equation} M_{\rm loss,C}= M_{\rm sf,C}\int _{m_1}^{m_2} [m - w(m)]\; \phi (m)\;d m, \end{equation} \tag{ 13 }$$

where w(m) is the mass of the remnant after the star with m dies. This expression can be readily extended in order to account for the entire star formation history of a given stellar population. This leads to the following stellar mass loss rate:

$$\begin{equation} \dot{M}_{\rm loss,C}= \int _{0}^{m_{\rm up}} \dot{M}_{\rm sf,C}[t-\tau (m)]\; [m - w(m)]\; \phi (m)\; d m, \end{equation} \tag{ 14 }$$

where τ(m) is the lifetime of stars with m and m_up is the IMF upper mass.

5.3.4. Metal Enrichment

The amount of metals ejected by stars into the ISM over their life and when they die as SN explosions depends on whether they are metal-rich or metal-poor. As mentioned, the yield of Pop III stars depends on their initial mass. According to the definition of f₁, f₂, and f₃ (Section 5.1), the mass fraction in massive (m > 130 M_☉) Pop III stars ending up locked into BHs is

$$\begin{equation} \beta _{\rm III}=f_3, \end{equation} \tag{ 15 }$$

while their yield is

$$\begin{equation} p_{\rm III}=0.5\,f_2. \end{equation} \tag{ 16 }$$

Once the hot gas metallicity reaches the value Z_crit, the H₂ density is no longer relevant as atomic cooling becomes the most efficient cooling mechanism. Then, Pop II stars begin to form in disks and spheroids. Pop I and II stars and the less massive Pop III stars liberate metals by type II SN explosions and, to a lesser extent, throughout their life. In AMIGA, we follow the mass loss of Pop I and II stars according to their specific evolution and compute the mass of metals they eject, which are supposed to mix with the cold ISM under the instantaneous recycling approximation (IRA; Tinsley 1980).

Equation (14) can also be readily adapted for the computation of the metal mass gained by ISM due to stellar evolution by means of the substitution

$$\begin{equation} m - w(m) \longrightarrow p(m)\,m, \end{equation} \tag{ 17 }$$

where p(m) is the yield of stars with m, which we can approximate by the global average value p, taken to be equal to 0.03,⁸ Thus, neglecting, according to IRA, the lifetime of massive stars, i.e., those essentially contributing to the chemical enrichment of the ISM, we arrive at the following metal mass loss rate into the ISM of component C due to stellar evolution:

$$\begin{equation} \dot{(ZM)}_{\rm loss,C}= p\,\dot{M}_{\rm sf,C}\!\!\int ^{m_{\rm up}}_{m_{\rm eff}} \!\!d m\,m\,\phi (m)\equiv p\,C_{\rm eff}\dot{M}_{\rm sf,C}, \end{equation} \tag{ 18 }$$

where m_eff = 10 M_☉ is the effective lower mass of stars producing metals, and C_eff is equal to 0.1 for the (modified) Salpeter IMF.

6.1. Inner Structure

The disk stability condition used in AMIGA is the simple global one provided by van den Bosch (1998). The shape of the disk of the central galaxy is computed self-consistently through the iterative procedure described in Mo et al. (1998) from the specific angular momentum of the gas at the cooling radius in the halo, taken to be equal to that of DM distributed according to the results of N-body simulations (e.g., Catelan & Theuns 1996; Bullock et al. 2001). This completely determines the scale length r_D or, equivalently, the central surface density, $\Sigma (0)=M_{\rm D}/(2\pi \ 0.83\ r_{\rm D}^2)$, of the exponential disk with a total mass equal to M_D. Hence, the disk structure is also set without introducing any free parameter. When a central disk galaxy becomes a satellite (see Section 6.2), it conserves its shape.

If the disk is unstable or its stability radius is smaller than the spheroid radius, then the cold gas coming from the halo directly goes into the spheroid. As the gas reaching the spheroid is the first to cool, it contributes to the spheroid with a very low angular momentum, which is neglected for simplicity. Some cold gas also reaches the central galaxy through captured satellites (see Section 6.2). The orbital momentum of satellites is assumed to be randomly distributed, so that such captures do not to alter (in average) the angular momentum of the disk.

Due to the lack of analytic treatment for violent relaxation, spheroids are the only systems whose inner structure cannot be causally linked to the cosmic properties set by cosmology. AMIGA assumes them using 3D density profiles of the Hernquist (1990) form (whose projection in two dimensions approximates the r^1/4 law) with scale length of r_B = r_e/1.81, where r_e is the effective (half-mass) radius of the two-dimensional profile. Spheroids forming with no gas, and hence suffering no dissipative contraction, are assumed to satisfy the relation $r_{\rm e}\approx A M_{\rm B}^{\gamma _{\rm B}}$ à la Kormendy (1977) between the effective radius r_e and total stellar mass M_B, with constants A and γ_B so as to recover the observed values of r_e of local spheroids with extreme masses (Shen et al. 2003; Guzman et al. 1997).⁹ Meanwhile, those forming with some amount of gas suffer dissipative contraction from the previous initial configuration during star formation. In the Appendix, we derive the following physically motivated differential equation for the dissipative contraction of the scale radius,

$$\begin{equation} r_{\rm B}^2(t)\,\frac{d r_{\rm B}}{d t}= -\frac{Z_{\rm cg,B}^{1/2}(t)M_{\rm cg,B}(t)}{Z_\odot ^{1/2}\rho _{\rm dis}\tilde{\tau }_{\rm acc}}, \end{equation} \tag{ 19 }$$

where $\tilde{\tau }_{\rm acc}$ is the universal time elapsed since the formation of the spheroid to the quenching of star formation due to the action of the AGN, and ρ_dis is a critical dissipation density, taken as a free parameter. When contraction is so marked that the density of the cold gas reaches the typical density (10⁶ particles cm⁻³) of dense molecular cores in local galaxies, AMIGA assumes that the gas cloud fragments to form stars without suffering any further contraction.

6.2. Galaxy Interactions

7.1. MBH Feeding

MBHs are supposed to arise from the BH remnants of very massive (m > 260 M_☉) Pop III stars, which coalesce in one mini-MBH per star cluster. In galaxy mergers or captures, the MBHs of the progenitor galaxies are assumed to migrate, by dynamical friction, to the center of the new spheroid where they form a binary system. Binary MBHs break if the recoil velocity produced by the emission of gravitational waves (in the Newtonian approximation; Fitchett 1983) is larger than the galaxy escape velocity (Blecha et al. 2011), in which case the less massive MBH escapes to the halo. Otherwise, the binary system quickly coalesces (Merritt et al. 2007).

MBHs also grow by accreting part of the gas that reaches their host spheroid. Spheroids collect matter in three different ways: by means of cooling flows of gas with low angular momentum, at wet mergers of similarly massive galaxies, and via disk-to-bulge mass transfers. Part of the gas loses angular momentum and reaches the central region where it feeds the central MBH.

Following Hatziminaoglou et al. (2003), AMIGA assumes that the gas mass accretion curve onto the MBH, $M_{\rm BH}^{\rm g}(t)$, scaled to the total accreted mass, has a bell-shaped universal form with a characteristic timescale of τ_acc (the AGN duty cycle; see below) equal to 0.1 Gyr. The only exception is at the beginning of gas cooling after halo formation if the central galaxy is a naked stellar spheroid. Then, the angular momentum of the falling gas is very small and there is no hindrance for the gas to directly reach the center of the main galaxy, and so the accretion rate onto the central MBH is simply taken to be equal to the cooling rate in the halo.

Apart from gas, MBHs accrete stars lying at the center of the spheroid (Milosavljević & Merritt 2003) at a rate of

$$\begin{equation} \dot{M}_{\rm BH}^{\rm s}=\alpha _{\rm BH}\frac{M_{\rm s,B}}{\tau _{\rm BH}}, \end{equation} \tag{ 20 }$$

where τ_BH is the dynamical time in the region of gravitational influence of the MBH, with a typical radius r_BH equal to $GM_{\rm BH}/\sigma _{\rm B}^2$, with σ_B being the stellar velocity dispersion in the spheroid, and where α_BH is the MBH feeding efficiency. In principle, α_BH should be taken as a free parameter, but for all reasonable values tried, the resulting stellar feeding is insignificant compared to the gas feeding, so we have taken it to simply be equal to 0.01.

As a result of all these feeding mechanisms, MBHs grow at the center of spheroids in such a way that they end up satisfying the observed Magorrian et al. (1998) relation between MBH masses and the stellar masses of the host spheroids.

7.2. AGN Feedback

AMIGA assumes that all AGNs have the same typical intrinsic spectrum independent of redshift. The continuum is described by two power laws crossing each other at a wavelength equal to 1100 Å (big blue bump). The optical spectral index in frequency (we define as "optical" the slope longward of 1100 Å) has a typical value of 0.5, and the UV index shortward of 1100 Å is equal to 1.76 (Wang et al. 1998). The most important emission lines (Ly_α, Ly_β, Mg ii, C iii, C iv, Si iv, H_α, H_β, and H_γ) and the small blue bump centered at ∼3000 Å are added to the above continuum with varying equivalent widths. From such a spectrum and the bolometric luminosity of any given AGN, inferred as explained in Section 9.2, one can readily compute its rest-frame extinction-free flux of ionizing photons and the associated energy.

AGNs contribute to the X-ray background with 0.04 of their bolometric luminosities (Vasudevan & Fabian 2007). However, the most important feedback from AGNs is the mechanical reheating of the gas inside galaxies. As has been noted, when new gas reaches the spheroid, a starburst takes place and the MBH begins to accrete gas. At about half the increasing branch of the MBH accretion curve, the gas reheated by the enlightened AGN begins to be expelled back into the halo, which will ultimately cause the quenching of the ongoing starburst (Springel 2005) and the braking of gas accretion onto the MBH. The total mass increase of the MBH is estimated as the mass of gas remaining in the spheroid at the maximum of the accretion curve minus the mass of gas reheated by the AGN and expelled back into the halo, with the reheating rate being given by

$$\begin{equation} \dot{M}_{\rm rh}^{\rm AGN}=\epsilon _{\rm AGN}\frac{2\,L(t)\,c^{-2}}{V_{\rm B}^2-V_{\rm hg}^2}, \end{equation} \tag{ 21 }$$

where _AGN is the quasar-mode (Bower et al. 2006) AGN heating efficiency, taken as a free parameter, and L(t) is the AGN bolometric luminosity (see Section 9.2). Of course, AGNs radiate at most at the Eddington limit, and so in low-mass MBHs, this reheating rate may not be enough to expel all of the gas remaining in the spheroid after the feeding of the MBH. Then, the starburst continues, with the dynamical timescale τ_{dyn, B} of the final contracted spheroid, until all of the gas is exhausted.

A second AGN feedback is the so-called radio-mode heating of the hot intrahalo gas (Croton et al. 2006) taking place when the MBH lies within a naked spheroid directly fed by cooling flows with small angular momentum (see Section 7.1). In this case, about one-tenth of the bolometric AGN luminosity is transferred mechanically to the hot gas in the halo through relativistic jets (Croton et al. 2006; Allen et al. 2006), which slows down the cooling of the hot gas there, possibly even halting it in the case of sufficiently massive MBHs. Such reheating is completely determined by the amount of cold gas feeding the MBH and the AGN radiation model described above, so it introduces no extra parameter.

Finally, AGNs also ionize and reheat the IGM outside halos. The escape fraction of ionizing photons from AGNs is taken to be equal to the previously mentioned escape fraction f_esc of ionizing UV photons from galaxies.

As a consequence of all the preceding processes, baryons circulate through the different phases and galactic components (see Figure 2).

Specifically, in periods between sudden mass changes due to halo mergers and galaxy captures and mergers, the masses of such phases and components in any given halo evolve smoothly according to the following set of differential equations:

$$\begin{eqnarray} {d M_{\rm hg}\over d t}&=& \sum _{i=1}^N\!\left \lbrace \!\dot{M}_{\rm rh}^{\rm AGN (i)}\!+\!\!\!\sum _{\mathrm{C}=B,D} \!\!\left[\dot{M}_{\rm rh,C}^{(i)}-\dot{M}_{\rm cool,C}^{(i)}\right]\!\!\right \rbrace \nonumber \\ &&+\dot{M}_{\rm hg}^{\rm in}-\dot{M}_{\rm hg}^{\rm out}\\ {d M_{\rm cg,C}^{(i)}\over d t}&\!=\!&\dot{M}_{\rm cool,C}^{(i)}\!-\!\dot{M}_{\rm rh,C}^{(i)}\!-\!\dot{M}_{\rm sf,C}^{(i)}\!+\!\dot{M}_{\rm loss,C}^{(i)}\!-\!\dot{M}_{\rm BH,C}^{\rm g (i)}\\ {d M_{\rm s,C}^{(i)}\over d t}&=&\dot{M}_{\rm sf,C}^{(i)}-\dot{M}_{\rm loss,C}^{(i)}-\dot{M}_{\rm BH,C}^{\rm s (i)}\\ {d M_{\rm BH}^{(i)}\over d t}&=& \dot{M}_{\rm BH}^{\rm g (i)}+\dot{M}_{\rm BH}^{\rm s (i)}, \end{eqnarray} \tag{ 22 }$$

where all of the rates denoted by a dot on the right are known functions of the evolving DM mass, M_H, provided by the EPS formalism, and the hot gas, cold gas, stellar and MBH masses are M_hg, $M_{\rm cg,C}^{(i)}$, $M_{\rm s,C}^{(i)}$, and $M_{\rm BH}^{(i)}$, whose evolution is being followed.

To render the notation in Equations (22)–(25) more compact, we have introduced the following definitions according to whether galaxies are the central object or satellites, and the galactic components are disks or spheroids: $\dot{M}_{\rm BH,D}^{\rm g (i)}=\dot{M}_{\rm BH,D}^{\rm s (i)}=\dot{M}_{\rm rh,D}^{\rm AGN (i)}=0$, while $\dot{M}_{\rm BH,B}^{\rm g (i)}=\dot{M}_{\rm BH}^{\rm g (i)}$ and $\dot{M}_{\rm BH,B}^{\rm s (i)}=\dot{M}_{\rm BH}^{\rm s (i)}$; in addition, $\dot{M}_{\rm cool,C}^{(i\ne 1)}=0$, and either $\dot{M}_{\rm cool,B}^{(1)}=0$, $\dot{M}_{\rm cool,D}^{(1)}=\dot{M}_{\rm cool}$, and $\dot{M}_{\rm sf,B}^{(1)}=\dot{M}_{\rm rh,B}^{(1)}=0$ or $\dot{M}_{\rm cool,B}^{(1)}=\dot{M}_{\rm cool}$, $\dot{M}_{\rm cool,D}^{(1)}=0$, and $\dot{M}_{\rm sf,D}^{(1)}=\dot{M}_{\rm rh,D}^{(1)}=0$, depending on whether or not star formation takes place in a stable disk.

Those mass exchanges between phases are, of course, accompanied by metal exchanges. Consequently, in periods between sudden mass changes due to captures and mergers, the metallicity of the hot gas, cold gas, and stars, Z_hg, $Z_{\rm cg,C}^{(i)}$, and $Z_{\rm s,C}^{(i)}$, in any given halo with mass M_H, and of the ionized IGM associated with it, $Z_{\rm IGM}^{\rm H}$ (the metallicity Z_IGM of the total ionized IGM is the result of the metal losses by all halos lying in ionized regions), evolve according to the following set of equations:

$$\begin{eqnarray} {d \left[Z_{\rm IGM}^{\rm H} M_{\rm IGM}^{\rm H}\right] \over d t}&=&-Z_{\rm IGM}\dot{M}_{\rm hg}^{\rm in}+Z_{\rm hg}\dot{M}_{\rm hg}^{\rm out}\\ {d [Z_{\rm hg}M_{\rm hg}] \over d t}&\!=\!&\sum _{\,i=1}^N \!\!\left \lbrace Z_{\rm cg,B}^{(i)}\dot{M}_{\rm rh,B}^{\rm AGN (i)}\!+\!\!\sum _{\mathrm{C}=\!B\!,D} \!\!\left [Z_{\rm cg,C}^{(i)}\dot{M}_{\rm rh,C}^{(i)}\right.\right.\nonumber \\ &&\left.\left.- Z_{\rm hg}\dot{M}_{\rm cool,C}^{(i)}\right ]\right \rbrace +Z_{\rm IGM}\dot{M}_{\rm hg}^{\rm in}-Z_{\rm hg}\dot{M}_{\rm hg}^{\rm out}\nonumber\\ \\ {d \left[Z_{\rm cg,C}^{(i)}{M_{\rm cg,C}^{(i)}}\right] \over d t}& =& Z_{\rm hg}\dot{M}_{\rm cool,C}^{(i)}-Z_{\rm cg,C}^{(i)}\dot{M}_{\rm rh,C}^{(i)}\nonumber \\ &&+\left[p\,C_{\rm eff}- Z_{\rm cg,C}^{(i)}\right]\,\dot{M}_{\rm sf,C}^{(i)}-Z_{\rm cg,C}^{(i)}\dot{M}_{\rm BH,C}^{\rm g (i)}\nonumber\\ \\ {d \left[Z_{\rm s,C}^{(i)}{M_{\rm sf,C}^{(i)}}\right]\over d t} &=& Z_{\rm cg,C}^{(i)}\,\dot{M}_{\rm sf,C}^{(i)}-Z_{\rm s,C}^{(i)}\dot{M}_{\rm BH,C}^{\rm s (i)}, \end{eqnarray} \tag{ 26 }$$

coupled to the previous set (Equations (22)–(25)), where $M_{\rm IGM}^{\rm H}$ is the mass of that part of the ionized IGM associated with the halo, equal to M_H times the current baryon mass fraction in ionized regions.

Strictly speaking, Equations (22)–(29) hold for halos harboring normal galaxies. In the case of primordial Pop III star clusters, the corresponding equations are somewhat different due to the fact that there is neither cooling in halos (stars photo-dissociate and even photo-ionize the hot gas) nor cold gas in galaxies. Then, the mass and metals lost by Pop III stars go directly into the hot gas, and the gas reheated through PISN leaves the halo, liberating metals into the surrounding IGM. In fact, these are essentially the only outflows from halos, opposed to the inflows mentioned in Section 2 and hence the only vector for the metal enrichment of the IGM (Rollinde et al. 2009; Greif et al. 2010; Wise et al. 2012). Indeed, the gas ejected from normal galaxies through type II SNe- and AGN-driven flows go into the halo where it enriches the metallicity of the hot gas (e.g., Scannapieco & Brüggen 2010; Springel 2005). In principle, it might also leave the halo, but according to the value of V_hg adopted in Equations (12) and (21), the reheated gas leaves the specific energy of the hot gas in the halo essentially unaltered,¹⁰ and so the possibility of those outflows is actually ignored.

AMIGA also follows the detailed exchanges of carbon. The reason for this is that carbon abundance is a more direct observable than metallicity Z, while the carbon mass fraction in the yields p_III and p of metal-poor and metal-rich stars are very different. For such a monitoring, we adopt the carbon mass fraction in the two yields provided by Schaerer (2002) and Ryan-Weber et al. (2009).

When reheated gas (with increased metallicity) returns to the halo, it takes some time to mix with the hot gas present there. In fact, during the smooth evolution of a halo, viscosity causes the gas ejected from the central spheroid (not from satellites) to remain stuck at the cooling front where it will be the next to cool.¹¹ The surviving reheated gas mixes with the hot gas when the halo merges or is accreted. Specifically, in a merger, the gas recently reheated is mixed with the hot gas lying in the inner h_rec fraction. If the halo is accreted, then the gas recently reheated is mixed with the outer 1 − h_rec fraction of the new halo. Only in the next major merger or accretion event is the surviving part of that reheated gas definitely mixed with the hot gas in the new halo. The hot gas recycling fraction h_rec is a free parameter of AMIGA.

Therefore, although this is not reflected in Equations (22)–(25), during periods of smooth evolution between halo captures and mergers, the hot gas in halos is stored in two separate compartments: the outer initial 1 − h_rec fraction, where the gas participates in inflows–outflows with the outer IGM, and the inner initial h_rec fraction, where it participates in cooling–reheating exchanges with galaxies. When the inner compartment is empty, cooling-reheating continues in the outer one.

In Figures 5 and 6, we show the evolution of the main cosmic mass densities and mass-weighted metallicities predicted by AMIGA in the same models as in Figure 4. We will come back to these figures in Section 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the previous sections, we described the modeling of the temperatures, metallicities, and ionized fractions of the various IGM phases, and the structural, chemical, kinematic, and dynamic properties of luminous objects. This is sufficient for some applications of the model (see, e.g., E. Salvador-Solé & A. Manrique 2015, in preparation). However, for other applications, the photometric properties of luminous objects must also be modeled.

9.1. Galaxy Luminosities

9.2. AGN Luminosities

To compute the photometric properties of AGNs, we not only need the masses of the associated MBHs and the rate at which they accrete matter, noted in Section 7.1, but also a radiation model of these objects. In AMIGA, we adopt the simple model developed by Hatziminaoglou et al. (2003) assuming that the radiative pressure onto the infalling gas, opposed to the gravitational pull by the MBH, produces damped oscillations that quickly reach a stationary regime. The AGN bolometric luminosity is then fully determined by the MBH feeding rate, according to

$$\begin{equation} L(t)=\frac{1}{2}\dot{M}_{\rm BH}^{\rm g} V_{\rm last}^2, \end{equation} \tag{ 30 }$$

where $V_{\rm last}^2\approx c^2/K_{\rm last}$ is the squared velocity at the last marginally stable Keplerian orbit around the MBH,¹² with radius equal to K_last = 9.2 times the Schwarzschild radius (Hatziminaoglou et al. 2003).

The bolometric luminosity function of AGNs is computed as the product of the time each source spends in the desired luminosity bin, known from the luminosity curve (Equation (30)), times the specific AGN reactivation rate (Hatziminaoglou et al. 2003). In the case of AGNs activated in galaxy mergers or in direct cooling flows, such a reactivation rate coincides with the formation rate of halos times the typical total number of central galaxy mergers or disk-instability episodes taking place after the halo forms. In the case of AGNs activated in tidal interactions among satellites, the distribution of bolometric AGN luminosities is directly related to the distribution of times elapsed since the corresponding (Poissonian distributed) satellite interactions.

Once the bolometric luminosities of AGNs are known, their spectrum described in Section 7.2 readily leads to their luminosities in any desired observer-frame photometric band, which is then properly corrected for extinction by H i according to Meiksin (2006). AMIGA assumes that the only period AGNs are visible, after correction for dust obscuration according to Gaskell et al. (2004), is after the accretion rate has reached its maximum value and the dispersion of gas is the most marked. Before that moment, they are completely enshrouded within the gas cloud.

AMIGA is a very complete, detailed, analytic model of galaxy formation devised to account fully consistently for the coupled evolution of luminous objects (galaxies and AGNs) and IGM since the dark ages. It incorporates molecular cooling and Pop III stars, the luminous objects with the most dramatic feedback, and accurately accounts for the intertwined evolution of the halo MF and the IGM temperature and ionization state.

AMIGA treats all those aspects of galaxy formation that can be causally linked to each other and to the underlying cosmology (DM clustering, halo and hot gas structure and kinematics, cooling, disk formation, and BH growth) without any free parameter. The only free parameters in the model,

Hot gas and IGM:

    b_cl: H ii clumping bias

    h_rec: hot gas recycling fraction

Pop III stars:

    p_III: yield of massive stars

    β_III: stellar mass fraction ending locked in BHs

Normal galaxies:

    α_G: star formation efficiency

    f_esc: escape fraction of ionizing photons

    ρ_dis: critical dissipation density

Galaxy interactions:

    Δ_i: minimum relative mass for interactions (0.01)

    χ_DB: disk-to-bulge mass transfer efficiency (0.01)

AGN:

    _AGN: quasar-mode heating efficiency

concern poorly known aspects (small-scale ionized gas distribution, stellar and AGN feedback, spheroid structure, and galaxy interactions) that are disconnected from each other.

Contrary to the usual procedure, no free parameters are used to specify the initial conditions (IGM metallicities, temperatures, and ionization state, MBH masses, ionizing UV fluxes, etc.). The modeling starts from trivial initial conditions at the dark ages, and follows the formation of the first generation galaxies with Pop III stars, characterized by only two parameters simultaneously fixing their IMF and feedback (p_III and β_III or, alternatively, f₂ and f₃; Equations (15) and (16)).

The fact that there is no artificial freedom or poorly motivated initial conditions that may spuriously facilitate the fitting of observations renders the predictions of AMIGA particularly reliable.

However, the number of free parameters is still quite large. A first analysis of the parameter space reveals that if the parameters characterizing Pop III stars adopt values outside some "acceptability ranges," then the IGM metallicity never becomes high enough to trigger the formation of normal galaxies and MBHs. The acceptability range of each of these parameters is very robust in the sense that it is independent of the value of its partner.

From the meaning of the different parameters and the way Pop III stars form, it is clear that the properties of these stars are fully determined by the parameters in the corresponding set. No other parameter can influence them. However, the situation is different for the properties of normal galaxies and MBHs. They depend, of course, on the parameters in their respective sets, but they may depend on the properties of Pop III stars as well. Indeed, Pop III stars are responsible for the metal enrichment of the IGM that is incorporated into halos. Since the higher the metallicity of the hot gas, the higher the cooling rate, the properties of Pop III may influence not only galactic metallicities, but also the amount of cold gas falling into galaxies, and hence the structural properties of these objects. On the other hand, Pop III stars reheat the IGM, increasing its entropy floor, and hence the minimum mass of halos able to trap gas. Since the larger the mass of a halo is, the smaller its concentration and the less intense is cooling, the amount of cold gas feeding the most abundant dwarf galaxies will be smaller. The question then arises: do the properties of Pop III stars influence those of normal galaxies?

To answer this fundamental question, next we compare the results of AMIGA drawn from the two plausible models differing only in the values of the Pop III star parameters, used in the previous figures. All of the remaining parameters are taken with identical values,¹³ so any possible difference in the final properties of normal galaxies will be due to the influence of Pop III star properties.

As shown in Figure 6, the higher metallicity of the ionized IGM found in the case of the top-heaviest Pop III star IMF (with the largest f₂ value) affects the metallicities of the hot gas, cold gas, and stars, which are a little higher in this model than in the one with less top-heavy Pop III star IMF (with the smallest f₂ value and, hence, lower IGM metallicity). However, normal galaxies, particularly dwarf ones, eject large amounts of metals in the hot gas, so its metallicity in any given halo quickly increases. Nonetheless, halos keep on accreting low-metallicity IGM and there are always new halos accreting gas for the first time, so the average hot gas metallicity increases very slowly, causing a similar trend in the metallicities of stars and cold gas in galaxies. Due to the permanently renewed memory of the IGM metallicity, the convergence of the metallicities in the two models is delayed until z ∼ 8 (a little later in the case of stars, due to their larger memory concerning their history, the cold gas mass is instead continuously renewed).

The structural properties of normal bright galaxies and MBHs, as traced by the masses of their galactic components and phases, show a similar behavior. As shown in Figure 5, they are even less sensitive to the properties of Pop III stars at very high z's. The reason for this is that atomic cooling is much more sensitive to the density of the hot gas than to its metallicity, while the hot gas density is very similar in the two models because of the very similar temperature of the corresponding ionized IGM (see Figure 4). Indeed, in the model with the top-heaviest Pop III star IMF, ionized bubbles grow more rapidly due to the larger abundance of massive Pop III stars, and hence the larger rate of ionizing photons. However, the temperature reached by IGM in bubbles is essentially the same in the two models because of the mass-independent SED of high-mass Pop III stars.

The structural properties of normal bright galaxies in the two models are so similar that despite the different masses of MBH seeds in the two models (through f₃), no significant difference is found in the MBH mass densities (see Figure 5). The reason for this is that although the masses of the coalesced Pop III BH remnants depend on the Pop III star IMF, as soon as spheroids begin to grow, MBHs accrete such large amounts of gas compared to the mass of their seeds that MBHs rapidly lose the memory of those seeds.

Therefore, one fundamental result of AMIGA is that the structural properties of normal bright galaxies and MBHs essentially decouple from those of Pop III stars; there is only a small coupling in the metallicities. In this sense, SAMs and simulations with non-self-consistent initial conditions should correctly predict the properties of normal bright galaxies and MBHs in the low and moderately high-z universe, or even in the high-z one provided we do not care about metallicities. This justifies previous studies on galaxy formation focusing on the properties of nearby galaxies achieved by means of models with inaccurate initial conditions. The situation is completely different, however, if one is interested in predicting accurate galactic metallicities or accurate properties, at any z, of small luminous objects (Pop III star clusters and normal dwarf galaxies), or even if one is interested in the evolution of the universe at very high z where the effects of Pop III stars are the most marked. Then, the use of a model like AMIGA is mandatory.

The results of the two models analyzed also show that spheroids grow in parallel to MBHs, and disks grow in parallel to spheroids, so that the MBH to spheroid and spheroid to disk mass ratios remain rather constant (see Figure 7). The former effect is the consequence of star formation in spheroids being quenched by ISM reheating by AGNs, whose bolometric luminosities are self-regulated by the MBH feeding. The latter is the consequence of the fraction of cold gas going into disks or spheroids, which depends on the spheroid mass.

Zoom In Zoom Out Reset image size

The constant MBH to spheroid mass ratio is first reached in massive galaxies because, when dwarf galaxies form, MBHs are tiny. This is why the average MBH to spheroid mass ratio shows a small increasing trend with decreasing z. Of course, for the stationary regime in both ratios to be reached, all galactic components must develop freely. At low z's (≲ 7), cooling becomes increasingly inefficient and disks begin to fade, while spheroids keep on growing through galaxy mergers that become increasingly dry. The spheroid to disk mass ratio then begins to increase with decreasing z. However, the MBH to spheroid mass ratio is kept unaltered. The reason for this is that in dry mergers, both MBHs and spheroids continue to grow, the former through the coalescence of the progenitor MBHs and the latter through the addition of the stellar populations of the merging galaxies. Just as the MBH to (stellar) spheroid mass ratio of the two progenitors is the same, so is the ratio in the final object. This explains why the MBH to spheroid mass ratio remains unchanged until z = 0.

Therefore, a second fundamental result of AMIGA is that the growth of all galactic components is ultimately controlled by that of MBHs. In fact, these objects play a crucial role in the evolution not only of normal galaxies, but also, through the feedback of those objects, in the evolution of all cosmic properties. For this reason, any model or simulation of galaxy formation must necessarily deal self-consistently with MBHs. This is true regardless of whether or not the model or simulation includes accurate MBH seeds.

The total number of free parameters in the list above is 10. However, some of them can be removed if one concentrates in the high-z universe (at, say, z ⩾ 2). Galaxy interactions indeed only play a significant role in the detailed morphological appearance of galaxies. As the galactic morphologies are unresolved at very high z, one can study in a first step the formation of the first luminous objects, taking those two parameters with fixed reasonable values (like those quoted in parentheses in the list above), and adjust them in a second step by studying the local universe.

The total number of free parameters in the first step then reduces to eight. Although this number is still considerably large, it is smaller than the number of independent data sets available on the cosmic properties at z ⩾ 2 (see E. Salvador-Solé & A. Manrique 2015, in preparation), so the problem is well constrained. Moreover, the fact that the structural properties of normal bright galaxies, the only observable at those z's, are independent of the properties of Pop III stars translates into the decoupling of their respective parameters. As shown in E. Salvador-Solé & A. Manrique (2015, in preparation), this notably simplifies the adjustment of the free parameters of the model. In fact, a very simple fitting procedure can then be devised where all the parameters are adjusted sequentially through the fit of one independent data set each. This has the added advantage of rendering the complex process of galaxy formation particularly easy to comprehend.

This work was supported by the Spanish DGES grant AYA2012-39168-C03-02 and the Catalan DIUE grant 2009SGR00217. G.B. acknowledges support for this work from the National Autonomous University of México, through grants IA102311 and IB102212-RR182212. We are thankful to Guy Mathez for his inestimable encouragement to start this project.

Dissipative contraction of spheroids during star formation is caused by the loss of orbital energy by the dense nodes of cold gas where stars form, hereafter the gas clouds, that move inside the spheroid. Gas clouds also radiate internal energy as they contract and fragment to form stars, but this energy loss does not alter the cloud orbits. For simplicity, we neglect the adiabatic contraction (or expansion, in the final gas ejection) of the dissipationless component (stars and DM).

The specific kinetic energy of the gas associated with the orbital velocity of clouds supporting the spheroid is

$$\begin{equation} e=-2f\frac{GM_{\rm B}}{r_{\rm B}}=-\frac{8\pi }{3} g G \bar{\rho }_{\rm B}r_{\rm B}^2, \end{equation} \tag{ A1 }$$

where $\bar{\rho }_{\rm B}$ is the mean spheroid density and f and g are two constants dependent on the specific spheroid density profile. The emission power per unit gas mass at the base of the dissipative contraction of the system can be assumed to satisfy the simple equation

$$\begin{equation} \frac{d e}{d t}=-\epsilon _{\rm dis}\frac{3kT_{\rm cg,B}}{2\mu m_{\rm p}\tau _{\rm dis}}, \end{equation} \tag{ A2 }$$

where _dis is the dissipation efficiency, T_{cg, B} is the effective temperature of the gas accounting for the orbital kinetic energy of clouds, and τ_dis is the dissipation timescale.

To determine the expression for τ_dis, we take into account that the orbital energy radiated per unit gas mass is also equal to

$$\begin{equation} \frac{d e}{d t}=- f_{\rm dis}\Lambda [T_{\rm cg,B}(t),Z_{\rm cg,B}(t)] \frac{\bar{n}_{\rm cg,B}(t)}{\mu m_{\rm p}}, \end{equation} \tag{ A3 }$$

where $\bar{n}_{\rm cg,B}$ is the mean particle density in the spheroid and f_dis is the fraction of the radiated energy that can be associated with the orbital motion of clouds. The cooling function Λ(T_{cg, B}, Z_{cg, B}) for a gas at T_{cg, B} of the order of 10⁵ K, corresponding to halos with relevant masses, and Z_{cg, B} spanning from 10⁻² Z_☉ to 1 Z_☉ is proportional to T_{cg, B} with the proportionality factor inversely proportional to the square root of Z_{cg, B} (Sutherland & Dopita 1993), that is, Λ(T_{cg, B}, Z_{cg, B}) ≈ 3kT_{cg, B}/[2n_cτ_c(Z_{cg, B})], where τ_c(Z_{cg, B}) ≈ τ_c(Z_☉)(Z_{cg, B}/Z_☉)^−1/2 for appropriate values of the characteristic number density n_c and time τ_c(Z_☉). We thus have $\tau _{\rm dis}\approx (\epsilon _{\rm dis}/f_{\rm dis})\tau _{\rm c}(Z_{\rm cg,B})n_{\rm c}/\bar{n}_{\rm cg,B}$.

Substituting this expressions for τ_dis into Equation (A2), taking into account the virial relation $3kT_{\rm cg,B}/(\mu m_{\rm p})=fGM_{\rm B}/r_{\rm B}=(4\pi /3) g G \bar{\rho }_{\rm B}r_{\rm B}^2$, we arrive at the following approximate equation for the scale radius at a time t after the beginning of star formation and dissipation,

$$\begin{eqnarray} r_{\rm B}^2\,d r_{\rm B}&=& -\frac{3\epsilon _{\rm dis}\,f}{16\pi f_{\rm dis}\,g} \left[\frac{Z_{\rm cg,B}(t)}{Z_\odot }\right]^{1/2}\frac{M_{\rm cg,B}(t)}{n_{\rm c}}\frac{d t}{\tau _{\rm c}(Z_\odot)}\nonumber\\ &=& -\frac{Z_{\rm cg,B}^{1/2}(t)M_{\rm cg,B}(t)}{Z_\odot ^{1/2}\rho _{\rm dis}}\frac{d t}{\tilde{\tau }_{\rm acc}}, \end{eqnarray} \tag{ A4 }$$

where $\tilde{\tau }_{\rm acc}$ is the universal time interval elapsed between the formation of the spheroid and the moment when star formation is quenched due to the reheating and ejection of the remaining gas by AGNs, and the parameter ρ_dis is a critical dissipation density where all (known and unknown) constant factors are encapsulated. Equations (A4), (23), and (28) for C=B determine the coupled evolution of the contracting scale radius and the cold gas mass and metallicity in the spheroid.