3.1. Gas components

The gas reservoirs in Dark Sage galaxies and haloes include:

• • The interstellar medium (ISM) or ‘cold gas’, which refers to gas in a galaxy’s disc. Each disc is broken into $N_\mathrm{ann}$ annular segments, each fixed a priori in terms of specific angular momentum. Each segment has fields for its mass, metallicity, and radial boundaries. The ISM is assumed to be axisymmetric and thin, but its radial structure is built numerically over time in the model. We use the subscript ‘cold’ to refer to ISM properties in equations, even though this medium includes some ionised gas.

• • The ‘hot gas’ reservoir, which represents the bulk of the circumgalactic medium (CGM), that is, the gas in the halo that is distinct from that in the galaxy. This reservoir is always assumed to be homogenous in terms of its temperature and metallicity, and therefore is described entirely by two fields: its total mass and total metal content. It is assumed to follow an analytic density profile. ‘Hot gas’ represents that which is available to be accreted onto the galaxy, subject to cooling (i.e. be it through the ‘hot mode’ or ‘cold mode’).

• • The ‘fountain’ reservoir, which accounts for gas in the CGM that has recently been reheated out of the ISM and is in the process of mixing with the hot gas. It is assumed to follow the same density profile and have the same specific energy as the hot-gas reservoir. We label this reservoir ‘fount’ in equations.

• • The ‘outflowing’ reservoir, abbreviated to ‘outfl’ in equations. Similar to the fountain reservoir, this represents gas in the CGM that has been reheated due to feedback. But its specific energy is higher, such that it will escape the halo after an appropriate time.

• • The ‘ejected’ reservoir, denoted by subscript ‘ejec’, which represents gas that has been expelled from the halo’s virial radius due to feedback (after transiting the outflow reservoir). It is reincorporated into the CGM over time.

• • The ‘local intergalactic medium’ (LIGM), which represents the remaining gas in the halo that is outside the virial radius of all its subhaloes. This is a single field per halo; in practice, it is associated with the central subhalo. The LIGM primarily provides a means of tracking metal enrichment in cosmological gas that can eventually be accreted by the central galaxy.

Whenever we generically refer to the ‘CGM’ in text and equations, we mean the sum of the hot, outflowing, and fountain reservoirs:

(2) \begin{equation}m_\mathrm{CGM} \equiv m_\mathrm{hot} + m_\mathrm{fount} + m_\mathrm{outfl}.\end{equation}

The temperature and density profile of all 3 CGM components is assumed to be the same.

3.2. Stellar components

The stellar components in Dark Sage galaxies and haloes include:

• • The stellar disc, which — like the ISM — is broken into a series of annuli. These annuli are fixed in terms of specific angular momentum in the same way as the ISM, but the stellar disc need not always be coplanar with the ISM. We denote fields concerning the stellar disc with the subscript ‘^*,disc’.

• • The instability-driven bulge, which — as its name suggests — is built from disc instabilities driving stellar mass into a galaxy’s centre. This component has zero angular momentum by definition, and is therefore supported against gravity entirely by random motions.

• • The merger-driven bulge, which is more akin to a classical galaxy bulge or ellipsoid. This constitutes stars that enter the bulge as a direct result of galaxy mergers. We calculate a nominal specific angular momentum for the merger-driven bulge.

• • Intrahalo stars, otherwise referred to as the stellar halo or ‘IHS’ for short. This is built up through tidal stripping and the disruption of satellite galaxies before they have the chance to merge. By definition, this exclusively tracks stars inside the virial radius of a subhalo. Satellites can have a non-zero IHS.

• • Local intergalactic stars (LIGS), the stellar analogue for the LIGM. It is built up similarly to IHS, but accounts exclusively for tidally stripped stellar mass in a halo that falls outside the virial radius. This field is always zero for satellites by construction.

One major infrastructural change we have made to the stellar components is that each one is now broken into $N_\mathrm{age}$ age bins. This means there are $(N_\mathrm{ann} + 4) \times N_\mathrm{age}$ stellar-mass and stellar-metallicity fields for each galaxy (the ‘ $+4$ ’ accounts for the spheroidal components and the LIGS).

The introduction of stellar-age bins serves several purposes that expand Dark Sage’s capabilities. For one, it means it is possible to reconstruct the formation history of any disc segment or bulge component of a galaxy (or the galaxy as a whole) from a single-snapshot output. This means Dark Sage can now predict age gradients in stellar discs and even three-dimensional age–metallicity–radius maps of populations of galaxies, which can be tested against modern observations (e.g. with survey data from integral field spectrographs). In this sense, Dark Sage galaxies have ‘multidimensional’ structure, á la the title of this paper. Component-wise star formation histories of Dark Sage galaxies are trivially reconstructed from a single-snapshot output, allowing for easy comparison with the outcome of models applied to the spectral energy distributions of survey galaxies by component (notably, e.g. Robotham, Bellstedt & Driver Reference Robotham, Bellstedt and Driver2022). The age bins also allow us to add detail to our treatment of stellar evolution and feedback, where we now model the distribution of supernovae as a function of time since birth for stellar populations (Section 8).

$N_\mathrm{age}$ is a user-defined integer. In this work, we take $N_\mathrm{age} = 30$ , which is high enough to ensure that the delayed feedback scheme is sufficiently converged with $N_\mathrm{age}$ (not shown here) and to recover the shape of galaxies’ star formation histories, but not so high that it prohibitively slows the code. The age bracket of each age bin is interpolated from the snapshot times of the underlying N-body simulation (see Section 3.5); $N_\mathrm{age}$ is therefore capped by the number of time steps in the merger trees. If $N_\mathrm{age} = N_\mathrm{snap} - 1$ , then each age bin has a one-to-one correspondence with snapshot intervals in the mergers trees. Later, in Figure 7, we show the age bins used in this work.

Whenever we refer to ‘the bulge’ without further specificity (or use subscript ‘^*,bulge’), we mean the combination of the instability- and merger-driven bulges. Whenever we refer to the total stellar content of a galaxy (or shorthand ‘^*’), we mean the combination of the stellar disc and both bulges, but not intracluster stars (nor the LIGS).

3.3. Disc annuli

We maintain the design of Stevens et al. (Reference Stevens, Croton and Mutch2016) in keeping the annular boundaries for Dark Sage discs as

(3a) \begin{equation}\frac{j_\mathrm{outer}^{(i)}}{\mathrm{kpc\,km\,s}^{-1}} = 1.4^{i-1}\, h^{-1},\end{equation}

(3b) \begin{equation}j_\mathrm{inner}^{(i)} =\left\{\begin{array}{l@{\quad}r}0, & i=1\\j_\mathrm{outer}^{(i-1)}, & i \geq 2\\\end{array}\right.,\end{equation}

(3c) \begin{equation}i \in \left[1, 2, ..., N_\mathrm{ann}-1, N_\mathrm{ann}\right],\end{equation}

with $N_\mathrm{ann} = 30$ . Throughout this paper, we use the superscript ‘(i)’ to refer specifically to the ith annulus. When no such superscript appears for a disc property, we are referring to that property for the whole disc. For convenience, in parts of this paper, we will also refer to the mean properties of an annulus, for example,

(4) \begin{equation}\left\langle j^{(i)} \right\rangle \equiv \frac{j_\mathrm{inner}^{(i)} + j_\mathrm{outer}^{(i)}}{2} .\end{equation}

The above applies for both stellar and gas discs. Whenever we refer to ‘the disc’, we mean the combination of the ISM and stellar disc.

3.4. BH components

We expand the consideration of where massive BHs can exist in Dark Sage. Each subhalo now is now granted several BH reservoirs:

• • The central BH of the galaxy maintains the definition as in previous versions of Dark Sage. It is what powers AGN feedback. Each galaxy is effectively assumed to hold a single central BH; when two galaxies merge, their BHs are assumed to merge too.

• • The ‘intrahalo black holes’ reservoir (IHBH) tracks the total number and summed mass of BHs that end up inside the halo but outside the galaxy itself. This reservoir grows through the disruption of satellite galaxies inside $R_\mathrm{200c}$ of their parent halo. This reservoir has no effect on galaxy evolution.

• • The reservoir for ‘local intergalactic black holes’ (LIGBH) similarly tracks the number and summed mass of BHs stripped from satellites outside the virial radius. Satellite galaxies themselves always have zero LIGBHs by construction.

3.5. Time-steps

In the vernacular of this paper, a full ‘time-step’ refers to the interval between snapshots in the merger trees of the N-body simulation that Dark Sage is run on. MTNG has 264 time-steps, starting at $z \simeq 29.3$ and finishing at $z = 0$ (the initial conditions were produced at $z = 63$ ). The separation between each pair of the first 31 snapshots is $\Delta \log_{10}(a) = 0.014$ , where a is the cosmological expansion factor. The interval in $\Delta \log_{10}(a)$ is halved for the next 33 snapshots and halved again for the remainder.

The movement of mass, metals, angular momentum, et cetera, between reservoirs in Dark Sage is ultimately governed by a series of coupled differential equations. These equations are solved numerically and sequentially. Rather than setting the ‘unit’ of time used in the numerical solution to that of the time-step in the merger trees, we break each time-step into three ‘sub-time-steps’ in Dark Sage.

The use of sub-time-steps helps mitigate the fact that we must chose an order in which we apply our astrophysical processes; in principle, most processes should happen simultaneously (though, some have a natural order in which they should be done). This way, we loop through all physical processes multiple times over a time-step, rather than once. Each sub-time-step within the same full time-step covers the same interval in cosmic time. We denote the cosmic-time interval of a sub-time-step as $\Delta t$ throughout this paper.

In earlier versions of Dark Sage, there were 10 sub-time-steps between each time-step. Those versions were based on the original Millennium simulation, which only had 64 snapshots in its merger trees, versus 265 in MTNG. The total number of sub-time-steps has therefore moderately increased from those versions (640 then versus 792 now).

Rather than explicitly writing out the network of coupled differential equations in Dark Sage altogether, we present equations in this paper in the form that they are solved within a sub-time-step. This best reflects how the equations are written in the codebase.

6.1. Velocity dispersion and disc height

Previously in Dark Sage, it was assumed that the velocity dispersion of gas in galaxy discs, $\sigma_\mathrm{cold}$ , was a universal constant for all cosmic time. We now introduce a redshift dependence to $\sigma_\mathrm{cold}$ , but maintain the (false but simple) assumption that for a given snapshot this is constant across all annuli of all galaxies. The function we employ for this is empirical, based on the average redshift evolution for observations of neutral gas shown in figure 6 of Übler et al. (Reference Übler2019):

(37) \begin{equation}\frac{\sigma_\mathrm{cold}(z)}{\mathrm{km\,s}^{-1}} = 11.0 + 11.3\,z.\end{equation}

Our redshift-zero value of $11\,\mathrm{km\, s}^{-1}$ is what had been previously assumed at all redshifts. The redshift coefficient is taken from equation 1 of Übler et al. (Reference Übler2019). Results from Pillepich et al. (Reference Pillepich2019) suggest it might be sensible to add a stellar-mass dependence to equation (37), but we leave such modifications for future model developments.

We assume that the vertical structure of discs follows a square-hypersecant profile, typical for those that are self-gravitating and isothermal (see Benítez-Llambay et al. Reference Benítez-Llambay, Navarro, Frenk and Ludlow2018). For a given annulus, we invoke

(38a) \begin{equation}\rho_\mathrm{cold}^{(i)}(\zeta) = \frac{\Sigma_\mathrm{cold}^{(i)}}{2\, \zeta_s^{(i)}}\, \mathrm{sech}^2\! \left(\frac{\zeta}{\zeta_s^{(i)}} \right), \end{equation}

(38b) \begin{equation}\zeta_s^{(i)} = \frac{\sigma_\mathrm{cold}^2(z)}{G\, \Sigma_\mathrm{cold}^{(i)}}.\end{equation}

Here, $\zeta$ represents height from the mid-plane of the disc, where $\zeta_s^{(i)}$ is the local scale height. By definition, the surface density of an annulus must be $\Sigma_\mathrm{cold}^{(i)} = \int_{-\infty}^{\infty} \rho_\mathrm{cold}^{(i)}(\zeta)\, \mathrm{d}\zeta$ . The average three-dimensional density of an annulus is therefore

(39) \begin{equation}\big\langle \rho_\mathrm{cold}^{(i)} \big\rangle = \frac{1}{\Sigma_\mathrm{cold}^{(i)}} \int_{-\infty}^{\infty} \rho_\mathrm{cold}^{(i)\,2} (\zeta)\, \mathrm{d}\zeta = \frac{G}{3} \left( \frac{\Sigma_\mathrm{cold}^{(i)}}{\sigma_\mathrm{cold}(z)} \right)^2 .\end{equation}

We do not differentiate between the dispersions (or scale heights) of ionised and neutral gas in Dark Sage gas discs, despite the theoretical expectation and empirical evidence that suggest these should differ (Übler et al. Reference Übler2019). Relatively speaking, this is not a great source of uncertainty in the model; the observed systematic difference between the neutral and ionised gas velocity dispersions of galaxies is much less than, for example, the observed scatter in either dispersion at fixed redshift, which we do not model.

When stars are formed (see below), they inherit the velocity dispersion of the gas that formed them. We record a stellar velocity dispersion for each stellar-age bin within each annulus. By construction, this means that old stars carry a higher velocity dispersion than young stars in galaxy discs. As the stellar-age bins are typically wider than the average snapshot interval in the merger trees, each episode of star formation requires updating the velocity dispersion of that age bin. This is done by taking a mass-weighted average of the squares of the velocity dispersions of the new and pre-existing stars.

6.2. Hydrogen fraction

In another update to the model, we adjust the mass fraction of hydrogen, X, in a given gas disc annulus. Naturally, where metallicities are high, one expects to find a lower hydrogen fraction. We adopt the fitting function of Stevens et al. (Reference Stevens2021, their equation B2), which accurately recovers the hydrogen fraction of gas cells in IllustrisTNG from metallicity alone (i.e. without needing prior knowledge of the helium mass fraction):

(40) \begin{equation}X =\left\{\begin{array}{l@{\quad}r}0.753 - 1.26\, Z, & Z \geq 0.025\\\mathrm{min} \left[0.76,\,0.762 - 2.3\, Z + 24.2\, Z^2 \right], & Z < 0.025\\\end{array}\right..\end{equation}

The helium mass fraction is then straightforwardly recovered as

(41) \begin{equation}Y = 1 - X - Z .\end{equation}

We impose a metallicity floor of $Z = 10^{-10}$ on all baryon reservoirs, which accounts for the trace amounts of lithium formed shortly after the Big Bang. Due to the ‘primordial lithium problem’ — see Fields (Reference Fields2011) — this is lower than what one would expect the initial metallicity of the Universe to be based on Big Bang nucleosynthesis alone. Metallicity is then built up from stellar enrichment (see below). In practice, the precise value of the metallicity floor is irrelevant, provided it is negligible relative to the yield of metals produced by a stellar population, which is of order $10^{-2}$ .

6.3. Neutral fraction

We introduce a new method to Dark Sage for calculating the fraction of the ISM that is neutral/ionised. Previously, this had been assumed to be a constant for all annuli and all galaxies. We now instead invoke a model of photoionisation equilibrium within each annulus, providing a predictive framework for the neutral fraction radial profiles of galaxies. This is done in two stages, considering the situations when (i) local sources of ionising radiation are dominant and (ii) the universal ionising background is instead dominant.

In what follows, we seek a solution for

(42) \begin{equation}f_\mathrm{neutral}^{(i)} \equiv \frac{ m^{(i)}_\mathrm{neutral} }{ m^{(i)}_\mathrm{cold} } = \frac{ m_{\rm H\,\small{I}}^{(i)} + m_\mathrm{H_2}^{(i)} }{ X^{(i)}\, m_\mathrm{cold}^{(i)} }.\end{equation}

In practice, we calculate this as

(43) \begin{equation}f_\mathrm{neutral}^{(i)} = \mathrm{min} \left[ f_\mathrm{neutral}^{\star(i)},\,f_\mathrm{neutral}^{\mathrm{bg}(i)} \right].\end{equation}

6.3.1. Local photoionisation

Let us start with the approximation demonstrated by Rahmati et al. (Reference Rahmati, Schaye, Pawlik and Raičević2013) that the local photoionisation rate of hydrogen in a galaxy scales linearly with its star formation rate. Writing this in terms of ionisation and star formation rate densities, we have

(44) \begin{equation}\dot{n}_\mathrm{ion} \simeq q_\gamma\, \rho_\mathrm{SFR}, \end{equation}

where $q_\gamma = 2 \times 10^{53}\, \mathrm{s^{-1} (M_\odot\, yr^{-1})^{-1}} = 6.31 \times 10^{60}\, \mathrm{M}_\odot^{-1}$ (Rahmati et al. Reference Rahmati, Schaye, Pawlik and Raičević2013). The question then becomes: how do we implement this on an annulus-by-annulus basis?

Similar to equation (39), we can calculate the average three-dimensional SFR density in terms of an annulus’s SFR surface density, given $\rho_\mathrm{SFR}^{(i)} (\zeta) \propto \rho_\mathrm{cold}^{(i)} (\zeta)$ :

(45) \begin{equation}\big\langle \rho_\mathrm{SFR}^{(i)} \big\rangle = \frac{1}{\Sigma_\mathrm{cold}^{(i)}} \int_{-\infty}^{\infty} \rho_\mathrm{cold}^{(i)} (\zeta)\, \rho_\mathrm{SFR}^{(i)} (\zeta)\, \mathrm{d}\zeta.\end{equation}

Approximating $\rho_\mathrm{SFR}$ in equation (44) as $\langle \rho_\mathrm{SFR}^{(i)} \rangle$ above and solving the integral, we take

(46) \begin{equation}\dot{n}_\mathrm{ion}^\mathrm{(i)} = \frac{q_{\gamma}\, G\, \Sigma_\mathrm{cold}^{(i)}\, \Sigma_\mathrm{SFR}^{(i)}}{3\, \sigma_\mathrm{cold}^2(z)}.\end{equation}

To solve for $f_\mathrm{neutral}^{(i)}$ , we need to balance ionisation with the recombination of atoms. The recombination rate density of hydrogen can be written as

(47) \begin{equation}\dot{n}_\mathrm{recom} = n_\mathrm{e}\, n_\mathrm{p}\, \alpha,\end{equation}

where $n_\mathrm{e}$ and $n_\mathrm{p}$ are the respective number densities of free electrons and free protons, and $\alpha = 4 \times 10^{-13}\,\mathrm{cm^3\,s^{-1}}$ is the recombination coefficient for hydrogen that we adopt. For an annulus, we consider average three-dimensional number densities, consistent with our earlier assumptions that there is no vertical variation in the neutral fraction (or ionisation–recombination balance):

(48a) \begin{equation}\dot{n}_\mathrm{recom}^{(i)} = \langle n_\mathrm{e}^{(i)} \rangle \, \big\langle n_\mathrm{p}^{(i)} \big\rangle\, \alpha,\end{equation}

(48b) \begin{equation}\big\langle n_\mathrm{e}^{(i)} \big\rangle = \frac{4}{3\,X^{(i)} + 1} \big\langle n_\mathrm{p}^{(i)} \big\rangle,\end{equation}

(48c) \begin{equation}\big\langle n_\mathrm{p}^{(i)} \big\rangle = \left( 1 - f_\mathrm{neutral}^{(i)} \right)\, X^{(i)}\, \frac{\big\langle \rho_\mathrm{cold}^{(i)} \big\rangle}{m_\mathrm{p}}\end{equation}

(48d) \begin{equation}= \left( 1 - f_\mathrm{neutral}^{(i)} \right)\, X^{(i)}\, \frac{G}{3\,m_\mathrm{p}} \left(\frac{\Sigma_\mathrm{cold}^{(i)}}{\sigma_\mathrm{cold}(z)} \right)^2 .\end{equation}

The difference between $\big\langle n_\mathrm{p}^{(i)} \big\rangle$ and $\langle n_\mathrm{e}^{(i)} \rangle$ is driven by extra free electrons sourced from ionised helium. For simplicity, we approximate that for every hydrogen atom that is split into a free electron and proton, there is a helium atom split into a free electron and He ii ion. In reality, it takes more energy to ionise a helium atom than a hydrogen atom, and the cross section of the two elements are not identical either, but these minutiae are unimportant for the detail modelled here. We also reasonably neglect any contributed free elections from He iii or ions of other elements, given their low cosmic abundances.

Under photoionisation equilibrium,

(49) \begin{equation}\dot{n}_\mathrm{ion}^{(i)} = \dot{n}_\mathrm{recom}^{(i)}.\end{equation}

Equation (48) serves as the right-hand side of equation (49). If we take equation (46) as the left-hand side, we are implicitly treating $f_\mathrm{neutral}^{(i)}$ as $f_\mathrm{neutral}^{\star (i)}$ . After combining these equations, we can rearrange and solve for the neutral fraction (when local photoionisation dominates) in terms of knowns:

(50) \begin{equation}f_\mathrm{neutral}^{\star (i)} = 1 - \frac{1}{X^{(i)}} \sqrt{ \frac{ 3 \left(3\, X^{(i)} + 1\right) q_\gamma\, m_\mathrm{p}^2\, \sigma_\mathrm{cold}^2(z)\, \Sigma_\mathrm{ SFR}^{(i)} }{ 4\, \alpha\, G\, \Sigma_\mathrm{cold}^{(i)\,3} } }.\end{equation}

For convenience, we note that the constants in this equation can be reduced to

(51) \begin{equation}\frac{3\, q_\gamma\, m_\mathrm{p}^2}{4\, \alpha\, G} = 0.2495\,\mathrm{cm^{-6}\,g^2\,s^3}.\end{equation}

One final decision remains: how do we quantify $\Sigma_\mathrm{SFR}^{(i)}$ for equation (50)? While we could tie this to the sub-time-step’s instantaneous SFR in the annulus, this would require a simultaneous solution for $\Sigma_\mathrm{SFR}^{(i)}$ and $f_\mathrm{ neutral}^{(i)}$ , because, as we will show in Section 7.3, local star formation explicitly depends on local molecular fraction, which requires knowing the neutral fraction first. But a better option exists. What equation (44) really symbolises is that young stars produce the majority of local photoionising radiation in galaxies. The age bins we have in each Dark Sage disc annulus effectively tell us how many young stars are present. That is, to be self-consistent with our treatment of stellar evolution (where any stars formed in the same age bin are treated as if they have an identical formation time — see the next section), we use the annulus’s average SFR for the age bin that the present sub-time-step falls in to calculate $\Sigma_\mathrm{SFR}^{(i)}$ for equation (50). This is both physically motivated and circumvents a circular solution for $f_\mathrm{neutral}^{\star (i)}$ .

We note that AGN should, in principle, add an additional source of local photoionisation. This is not something we have considered for this version of Dark Sage. We leave it for future developments of the model.

6.3.2. Photoionisation background

We adopt a redshift-dependent, universal photoionising radiation background for neutral hydrogen to place a ceiling on the neutral fraction of the ISM of Dark Sage galaxies. When this source of photoionisation dominates over local sources, we impose

(52) \begin{equation}\dot{n}_\mathrm{ion} = \Gamma_{\rm H\,\small{I}}^\mathrm{eff}(z)\, n_\mathrm{H, neutral}.\end{equation}

Here, the effective rate of hydrogen photoionisation, $\Gamma_{\rm H\,\small{I}}^\mathrm{eff}(z)$ , is taken from Table D1 of Faucher-Giguère (Reference Faucher-Giguère2020). In practice, for Dark Sage disc annuli, we effect

(53) \begin{equation}\dot{n}_\mathrm{ion}^{(i)} = \Gamma_{\rm H\,\small{I}}^\mathrm{eff}(z)\, f_\mathrm{neutral}^{\mathrm{bg}(i)}\, X^{(i)}\, m_\mathrm{p}^{-1}\, \big\langle \rho_\mathrm{ cold}^{(i)} \big\rangle .\end{equation}

Combining this with equations (48 & 49), after some omitted algebra, we end up with a quadratic equation for $f_\mathrm{neutral}^{\mathrm{bg}(i)}$ in terms of knowns. The requirement that $f_\mathrm{neutral}^{\mathrm{bg}(i)} \leq 1$ by definition eliminates one of the quadratic solutions, leaving

(54a) \begin{equation}f_\mathrm{neutral}^{\mathrm{bg}(i)} = 1 + \gamma - \sqrt{\gamma^2 + 2\, \gamma},\end{equation}

(54b) \begin{equation}\gamma \equiv \frac{3 \left(3\, X^{(i)} + 1 \right) m_\mathrm{p}\, \Gamma_{\rm H\,\small{I}}^\mathrm{eff}(z)}{8\, \alpha\, G\, X^{(i)}} \left( \frac{\sigma_\mathrm{ cold}(z)}{\Sigma_\mathrm{cold}^{(i)}} \right)^2 .\end{equation}

6.4. Molecular fraction

We return to using a formula based on the mid-plane pressure of each disc annulus, $P_\mathrm{mid}^{(i)}$ , as the default for calculating the molecular fraction of gas. Similar to the implementation in Stevens et al. (Reference Stevens, Croton and Mutch2016), we use the equations from Blitz & Rosolowsky (Reference Blitz and Rosolowsky2004) and Elmegreen (Reference Elmegreen1989):

(55a) \begin{equation}f_\mathrm{H_2}^{(i)} = \left[ \left(\frac{P_\mathrm{mid}^{(i)}}{5.93 \times 10^{-13}\, \mathrm{Pa}}\right)^{-0.92} + 1 \right]^{-1},\end{equation}

(55b) \begin{equation}P^{(i)}_\mathrm{mid} =\left\{\begin{array}{l@{\quad}r}\dfrac{\pi\, G\, \Sigma_\mathrm{cold}^{(i)}}{2} \left(\Sigma_\mathrm{cold}^{(i)} + \dfrac{\sigma_\mathrm{cold}(z)\, \Sigma_*^{(i)}}{\sigma_*^{(i)}} \right), & \theta_\mathrm{ disc} \leq 10^\circ \\[15pt]\dfrac{\pi\, G}{2} \Sigma_\mathrm{cold}^{(i)\,2} \,, & \theta_\mathrm{disc} > 10^\circ\end{array}\right. .\end{equation}

The only practical changeFootnote i from Stevens et al. (Reference Stevens, Croton and Mutch2016) is how $\sigma_\mathrm{cold}$ and $\sigma_*^{(i)}$ are calculated, which is now consistent with that described above in this paper. $\theta_\mathrm{disc}$ is the angular misalignment between the stellar and gas discs.

Other prescriptions for $f_\mathrm{H_2}^{(i)}$ are available in the code, such as those based on the works of McKee & Krumholz (Reference McKee and Krumholz2010) and Gnedin & Draine (Reference Gnedin and Draine2014). We do not adopt them in this paper, though, as we found in testing that extra free parameters were likely required to calibrate the model to observations with these prescriptions.

In the top two panels of Figure 4, we show the H i and H $_{2}$ fractions of Dark Sage galaxies as a function of stellar mass. These relations are predictions of the model, and are compared against the representative xGASS and xCOLD GASS samples of galaxies (Saintonge et al. Reference Saintonge2017; Catinella et al. Reference Catinella2018). We note that this raw comparison is simply for a rough reference; to make a fair, accurate comparison between and any simulations and this dataset requires deeper considerations, as has been explored in detail in Stevens et al. (Reference Stevens2019a, Reference Stevens2021).

Figure 4. The planes of stellar mass versus H i fraction (top panel), H $_{2}$ fraction (middle panel), and specific star formation rate (bottom panel) in Dark Sage galaxies at $z = 0$ . Representative observational samples with $m_* \geq 10^9\,\mathrm{M}_\odot$ are compared from xGASS for H i (Catinella et al. Reference Catinella2018), xCOLD GASS for H $_{2}$ (Saintonge et al. Reference Saintonge2017), and SDSS for sSFR (the volume-limited sample used in Brown et al. Reference Brown2017). Thick lines represent medians. Thin lines are 16th and 84th percentiles. For xGASS and xCOLD GASS, the lines assume non-detections carry their upper-limit value, while the shaded region covers the 16th–84th interpercentile range assuming non-detections have zero mass of that type. SDSS contours approximately encapsulate 38, 68, and 95 per cent of the sample (thicker to thinner). Dark Sage hexbins and solid lines in the bottom panel use the time-step-averaged SFR (almost instantaneous), while the dashed lines represent the average SFR in the last age bin ( $\sim$ 1 Gyr).

7.1. Calculating (in)stability

Each annulus of each disc of each galaxy is routinely checked for its gravitational stability according to Toomre’s Q criterion. We calculate this local stability parameter first for each of the stellar and gas components separately, following the respective equations

(56a) \begin{equation}Q_*^{(i)} = \frac{ \kappa^{(i)}\, \sigma_\mathrm{*,disc}^{(i)} }{ 3.36\, G\, \Sigma_*^{(i)} }, \end{equation}

(56b) \begin{equation}Q_\mathrm{cold}^{(i)} = \frac{ \kappa^{(i)}\, \sigma_\mathrm{cold}(z) }{ \pi\, G\, \Sigma_\mathrm{cold}^{(i)} }, \end{equation}

(56c) \begin{equation}\kappa^{(i)} \equiv \sqrt{ \frac{ 2\, \big\langle v_\mathrm{circ}^{(i)} \big\rangle \left( j_\mathrm{outer}^{(i)} - j_\mathrm{inner}^{(i)} \right) }{ \langle r^{(i)} \rangle^2 \left( r_\mathrm{outer}^{(i)} - r_\mathrm{inner}^{(i)} \right) } }\end{equation}

(Toomre Reference Toomre1964; Binney & Tremaine Reference Binney and Tremaine1987; Pringle & King Reference Pringle and King2007). Provided $\theta_\mathrm{disc} \leq 10^\circ$ , $Q_*^{(i)}$ and $Q_\mathrm{cold}^{(i)}$ are then combined into a total stability, $Q^{(i)}$ , through

(57a) \begin{equation}Q^{(i)} = \left(\frac{1}{\mathrm{min} \left[ Q_\mathrm{cold}^{(i)}, Q_*^{(i)} \right] } + \frac{W^{(i)}}{\mathrm{max} \left[ Q_\mathrm{cold}^{(i)}, Q_*^{(i)} \right] } \right)^{-1}, \end{equation}

(57b) \begin{equation}W^{(i)} \equiv \frac{2\, \sigma_*^{(i)}\, \sigma_\mathrm{cold}(z)}{ \sigma_*^{(i)\,2} + \sigma_\mathrm{cold}^2(z)}\end{equation}

(Romeo & Wiegert Reference Romeo and Wiegert2011). If $Q^{(i)} < 1$ , the annulus is gravitationally unstable.

By definition, it is not physically possible for something to remain in an unstable state for an extended period of time. This demands that something must happen to bring the annulus back to $Q^{(i)} = 1$ . We therefore seek to resolve instabilities as soon as they occur.

Physically, one would expect instabilities to produce star formation from the collapse of unstable gas. Instabilities can also cause a rearrangement of mass in a galaxy. In principle, that rearrangement of mass can break the symmetry of a galaxy’s gravitational potential, for example, through creating a bar, which can create a radially variant torque on the disc, thereby redistributing angular momentum (see, e.g. Athanassoula Reference Athanassoula2003). This would, in turn, establish radial migration in the galaxy’s disc and alter the orbits of stars, which could also affect the stellar velocity dispersion (and disc scale height) of the disc and create a pseudobulge (see Kormendy & Kennicutt Reference Kormendy and Kennicutt2004, and references therein). In Dark Sage, we do not explicitly model bars nor asymmetric gravitational potentials, while pseudobulges are inherently considered to be part of the disc (for relevant discussion, see Stevens et al. Reference Stevens, Croton and Mutch2016). We nevertheless do model the end result that instabilities have on discs, namely in the form of star formation, radial migration, and disc heating. We do this in a way that is agnostic to the details of why each aspect of instability resolution happens the way it does, as it ultimately does not matter for how the model is constructed.

There are two ways to return an unstable disc annulus to stability: either mass must be removed from the annulus or its velocity dispersion must be raised. We have already prescribed $\sigma_\mathrm{cold}$ to be a strict function of z, meaning the former must be how gaseous instabilities are resolved. Cold gas can be removed from an annulus in one of two ways: by transferring it to adjacent annuli, or by induced star formation triggering stellar feedback. Stellar mass can only be removed by transferring it to adjacent annuli. However, we allow the velocity dispersion of stars to also rise. Each phase therefore has two modes of instability resolution. The fractional importance of each mode is controlled by a single free parameter, dubbed $f_\mathrm{move}^\mathrm{gas}$ .

The first step to resolving an instability is deciding which baryon phase the onus is on to solve it. A return to stability ( $Q^{(i)} \geq 1$ ) can be guaranteed once both $Q_\mathrm{cold}^{(i)} \geq 1 + W^{(i)}$ and $Q_*^{(i)} \geq 1 + W^{(i)}$ . If one phase already satisfies this requirement, the onus is entirely on the other phase to resolve the instability. If neither criterion is met, both phases have to do some work.

If $\theta_\mathrm{disc} > 10^\circ$ , the stability of the gas and stellar discs are treated entirely independently. Respective instabilities are then triggered when $Q_\mathrm{cold}^{(i)} < 1$ and $Q_*^{(i)} < 1$ .

7.2. Resolving gas instabilities

Once a gas instability is identified, it is resolved in two steps. Both steps involve removing gas from the annulus. The amount of unstable gas that must be removed from that annulus to return it to stability is

(58) \begin{equation}m_\mathrm{cold, unstable}^{(i)} = m_\mathrm{cold}^{(i)} \left( 1 - \frac{Q_\mathrm{cold}^{(i)}}{Q_\mathrm{cold, min}^{(i)}} \right),\end{equation}

where $Q_\mathrm{cold,min}^{(i)}$ is greater than $Q_\mathrm{cold}^{(i)}$ (for a stable disc, $Q_\mathrm{cold}^{(i)} = Q_\mathrm{cold, min}^{(i)}$ by our definition) and typically equal to 1 or $1+W^{(i)}$ , but can vary based on the circumstance of the instability trigger as described in Section 7.1. The first step in resolving this instability is to redistribute gas in the disc. This is done by moving a controlled fraction of the unstable mass to the two adjacent annuli, in proportion such that angular momentum is conserved. Mathematically, we can write this as

(59) \begin{equation}\Delta m_\mathrm{cold}^{(i) \rightarrow (i \pm 1)} = \frac{ \left\langle j^{(i)} \right\rangle - \left\langle j^{(i \mp 1)} \right\rangle }{ \pm \left\langle j^{(i+1)} \right\rangle \mp \left\langle j^{(i-1)} \right\rangle } f_\mathrm{move}^\mathrm{gas}\, m_\mathrm{cold, unstable}^{(i)} .\end{equation}

Whenever mass is transferred to adjacent annuli in this way, it carries with it metals in proportion to its metallicity. Strictly, equation (59) only applies for $i \in [2,N_\mathrm{ann}-1]$ . For the outermost annulus, mass transfer can only happen inwards, meaning we abandon conservation of angular momentum, and simply adopt

(60) \begin{equation}\Delta m_\mathrm{cold}^{(N_\mathrm{ann}) \rightarrow (N_\mathrm{ann}-1)} = f_\mathrm{move}^\mathrm{gas}\, m_\mathrm{cold, unstable}^{(N_\mathrm{ann})} .\end{equation}

For the innermost annulus, any inward-moving gas is assumed to be accreted by the BH. We approximate the gas accreted onto the BH to have lost all its angular momentum. We still adhere to angular momentum conservation for outward- and inward-moving gas in this case, meaning

(61a) \begin{equation}\Delta m_\mathrm{cold}^{(1) \rightarrow (2)} = \frac{ \left\langle j^{(1)} \right\rangle}{ \left\langle j^{(2)} \right\rangle } f_\mathrm{move}^\mathrm{gas}\, m_\mathrm{cold, unstable}^{(1)}, \end{equation}

(61b) \begin{equation}\Delta m_\mathrm{cold \rightarrow BH}^{(1)} = \frac{ \left\langle j^{(2)} \right\rangle - \left\langle j^{(1)} \right\rangle }{ \left\langle j^{(2)} \right\rangle } f_\mathrm{move}^\mathrm{gas}\, m_\mathrm{cold, unstable}^{(1)} .\end{equation}

Note that the growth of the BH through this channel is independent of its present mass. This is how Dark Sage is able to seed zero-mass BHs. We discuss the consequences of this ‘quasar’ form of BH accretion in Section 9.

The movement of unstable gas through the disc to the innermost annulus and subsequently onto the BH is the primary means by which BH accretion occurs. The only other way for BHs to gain a substantial amount of mass is when they merge together. An important test of this framework is whether Dark Sage can recover the well-known relation between a galaxy’s central BH mass and bulge mass. What makes this particularly exciting is that this is now a prediction of the model; as we show in Section 13, this relation has no explicit bearing on the calibration of Dark Sage’s free parameters. By design, instabilities and mergers can simultaneously grow both the BH and bulge in Dark Sage (if both gas and stars are involved), though this does not guarantee quantitative agreement with the empirical relation between these two properties.

We show in Figure 5 that Dark Sage predicts a BH–bulge mass relation that overlaps favourably with observations (cf. Scott, Graham & Schombert Reference Scott, Graham and Schombert2013). Interestingly, Dark Sage’s relation is not unimodal; this is seen in the intensity of the hexbins and clarified by the misalignment of equivalent running percentiles when the Dark Sage data are binned along each of the two axes. No cuts have been applied to Dark Sage here. A deeper assessment of this relation relative to observations is left for future work.

Figure 5. Black hole–bulge mass relation for Dark Sage galaxies at $z = 0$ . The x-axis is the sum of the merger- and instability-driven bulge masses. Because this is not a simple unimodal relation in Dark Sage, we show running percentiles (thick for median, thin for 16th and 84th) when the model’s data are binned along the x-axis (solid) and y-axis (dotted). Neither median traces the high-number-density ridge of Dark Sage galaxies. Compared is a compilation of observational data from Scott et al. (Reference Scott, Graham and Schombert2013). While previous versions of Dark Sage calibrated to these data, this outcome is now a prediction of the model.

Figure 6. The resolved molecular Kennicutt–Schmidt relation for Dark Sage galaxies at $z = 0$ . The dashed Dark Sage line is the median relation for all annuli. The hexbins are grey-scaled according to the logarithmic number of annuli within them. For a fairer comparison to observations, the hexbins only count annuli with $\Sigma_* > 10\,\mathrm{M}_\odot\,\mathrm{pc}^{-2}$ from galaxies with $10^{8.9} < m_*/\mathrm{M}_{\odot} < 10^{11}$ and $0.02 < m_\mathrm{H_2} / m_{\rm H\,\small{I}} < 1.13$ in (sub)haloes with at least 200 particles at one point in their history. The solid lines represent the median (thick) and 16th and 84th percentiles (thin) for those binned annuli after weighting each annulus by its area. This weighting further improves the comparison to observations, which typically use pixels of fixed length $\simeq$ 1 kpc. The best-fitting relation from observations (Bigiel et al. Reference Bigiel2011) and equivalent scatter is shown by the shaded region. Contours encapsulate 38, 68, and 95 per cent of pixels across the HERACLES and VERTICO surveys that are detected in both axes.

After moving gas, part of the remaining unstable gas is then consumed in star formation, with the last of it reheated out of the disc by the instantaneous component of the stellar feedback associated with those new stars:

(62) \begin{align}\Delta m_\mathrm{cold \rightarrow *}^{(i)\operatorname{i-SF}} + \Delta m_\mathrm{cold \rightarrow fount}^{(i)\operatorname{i-SN}} + \Delta m_\mathrm{cold \rightarrow outfl}^{(i)\operatorname{i-SN}} = \left(1 - f_\mathrm{move}^\mathrm{gas}\right)\, m_\mathrm{cold, unstable}^{(i)} .\end{align}

The balance between each of the terms on the left-hand side of this equation is described in Section 8. Note that $\Delta m_\mathrm{cold \rightarrow *}^{(i)\operatorname{i-SF}}$ represents the net production of stellar mass at the end of the present sub-time-step, and is therefore less than the annulus’s instability-driven star formation rate multiplied by $\Delta t$ (this is also explained in Section 8). The superscripts i-SF and i-SN are intended to differentiate instability-driven star formation and supernovae from the passive channel of SF (Section 7.3).

After a gas instability is resolved, we recalculate $Q_*^{(i)}$ . This is necessary, as the stellar mass in the annulus will have changed, which inevitably (temporarily) raises not only the probability of a stellar instability triggering but also the mass of unstable stellar mass.

7.3. Passive star formation

After resolving instabilities (including what follows in Section 7.4), we allow for a second mode of star formation that depends explicitly on the local molecular gas surface density. This ‘passive’ mode of star formation in Dark Sage obeys

(63) \begin{equation}\Sigma_{\operatorname{p-SFR}}^{(i)} = \epsilon_\mathrm{H_2}\, \Sigma_\mathrm{H_2}^{(i)}, \end{equation}

where $\epsilon_\mathrm{H_2}$ is the star formation efficiency of H $_{2}$ , equivalent to the inverse of a depletion timescale. In a change from Stevens et al. (Reference Stevens, Croton and Mutch2016), $\epsilon_\mathrm{H_2}$ is no longer treated as a free parameter. Instead, this is fixed to $\epsilon_\mathrm{H_2} = (2.35\,\mathrm{Gyr})^{-1}$ , consistent with the empirical best-fitting relation from Bigiel et al. (Reference Bigiel2011). Of course, there must be a limit to the amount of star formation in a given annulus set by the total mass of gas present, accounting for the gas that must be removed from the annulus due to feedback from the largest stars that form in that population and effectively go supernova immediately. In practice, we therefore effect

(64) \begin{equation}\Delta m_\mathrm{cold \rightarrow *}^{(i)\operatorname{p-SF}} = \mathrm{min} \bigg[ \pi \left( r_\mathrm{outer}^{(i)\,2} - r_\mathrm{inner}^{(i)\,2} \right) \left( 1 - f_\star^\mathrm{ return} \right) \Sigma_{\operatorname{p-SFR}}^{(i)}\, \Delta t, \\m_\mathrm{cold}^{(i)} - \Delta m_\mathrm{cold \rightarrow fount}^{(i)\operatorname{p-SN}} - \Delta m_\mathrm{cold \rightarrow outfl}^{(i)\operatorname{p-SN}} \bigg] .\end{equation}

Many of the terms here will be defined in Section 8.

Despite equation (63) explicitly tying star formation and H $_2$ surface density together, Dark Sage predicts a non-trivial resolved molecular Kennicutt–Schmidt (RMKS) relation, that is, $\Sigma_\mathrm{H_2}$ versus $\Sigma_\mathrm{SFR}$ . This is predictive because (i) equation (63) only applies to stars formed in stable annuli, (ii) the majority of star formation happens in unstable annuli, and (iii) there is no explicit tying of the molecular fraction to disc stability. In Figure 6, we show the RMKS relation for (a relevant subset of) Dark Sage galaxies at $z = 0$ . This is compared to the relation fitted to observations by Bigiel et al. (Reference Bigiel2011) as well as data from the HERACLES and VERTICO surveys (Leroy et al. Reference Leroy2009; Brown et al. Reference Brown2021, respecitvely). Note that, despite HERACLES and VERTICO respectively probing field and cluster environments, these two datasets have been demonstrated to have consistent RMKS relations (Jiménez-Donaire et al. Reference Jiménez-Donaire2023).Footnote j

7.4. Resolving stellar instabilities

Similar to gas, there are two steps in resolving a stellar instability. But these steps are different in detail. The mass of stars deemed necessary to move out of the annulus to restore stability is set by

(65) \begin{equation}m_\mathrm{*, unstable}^{(i)} = m_\mathrm{*,disc}^{(i)}\, f_\mathrm{move}^{*(i)} \left( 1 - \frac{Q_\mathrm{*}^{(i)}}{Q_\mathrm{*, min}^{(i)}} \right),\end{equation}

(66) \begin{equation}f_\mathrm{move}^{*(i)} = 1 - \frac{\sigma_\mathrm{cold}(z)}{\sigma_\mathrm{*,disc}^{(i)}} \left( 1 - f_\mathrm{move}^\mathrm{gas} \right) .\end{equation}

By construction, it is always true that $\sigma_\mathrm{cold}(z) \leq \sigma_\mathrm{*,disc}^{(i)}$ , and therefore $f_\mathrm{move}^{*(i)} \in (0,\,1)$ . The unstable mass of stars is moved to adjacent annuli in the same j-conserving fashion as done for gas:

(67) \begin{equation}\Delta m_\mathrm{*,disc}^{(i) \rightarrow (i \pm 1)} = \frac{ \left\langle j^{(i)} \right\rangle - \left\langle j^{(i \mp 1)} \right\rangle }{ \pm \left\langle j^{(i+1)} \right\rangle \mp \left\langle j^{(i-1)} \right\rangle } m_\mathrm{*, unstable}^{(i)}\end{equation}

when $i \in [2, N_\mathrm{ann}-1]$ . For the outermost annulus, we again neglect angular momentum conservation by invoking

(68) \begin{equation}\Delta m_\mathrm{*,disc}^{(N_\mathrm{ann}) \rightarrow (N_\mathrm{ann}-1)} = m_\mathrm{*, unstable}^{(N_\mathrm{ann})} .\end{equation}

Inward-moving stars from the innermost annulus are transferred to the instability-driven bulge. This bulge component has no angular momentum by design. As such, we effect

(69a) \begin{equation}\Delta m_\mathrm{*,disc}^{(1) \rightarrow (2)} = \frac{ \left\langle j^{(1)} \right\rangle}{ \left\langle j^{(2)} \right\rangle }\, m_\mathrm{*, unstable}^{(1)}, \end{equation}

(69b) \begin{equation}\Delta m_\mathrm{*,disc \rightarrow \operatorname{i-bulge}}^{(1)} = \frac{ \left\langle j^{(2)} \right\rangle - \left\langle j^{(1)} \right\rangle }{ \left\langle j^{(2)} \right\rangle }\, m_\mathrm{*, unstable}^{(1)} .\end{equation}

As stars transfer from one annulus to the next (or to the instability-driven bulge), they carry with them the same age–metallicity–velocity dispersion distribution. The metallicity and velocity dispersion squared of each age bin in the receiving stellar component are updated with a mass-weighted sum of the moving stars and that already in the receiving annulus.

By itself, the motion of stellar mass prescribed above does not raise $Q_*^{(i)}$ to its minimum stable value. To complete this, we invoke ‘heating’ of the stellar disc. That is, the stellar velocity dispersion of the disc annulus is raised by a factor equal to

\begin{align*}f_\mathrm{move}^* + \left(1 - f_\mathrm{move}^* \right)\, \frac{Q_\mathrm{*, min}^{(i)}}{Q_\mathrm{*}^{(i)}} \, .\end{align*}

This change applies to each age bin individually. The functional form of equation (66) is something we have put in by hand intentionally so that disc heating does not continue to occur unregulated in a perpetually unstable disc. The increase to the annular stellar velocity dispersion is applied after moving the unstable mass of stars.

It is because stars are collisionless and gas collisional that we have different $f_\mathrm{move}$ values — that is, different radial-migration rates — for gas and stars. The functional form of equation (66) ensures that stellar velocity dispersions in the disc cannot grow unchecked. The more stellar discs are heated, the less they will be heated at the next instability.

The size of an instability-driven bulge is set by the velocity dispersion of the stars it acquires. This is calculated from the Jeans equation, under the assumptions that the instability-driven bulge is both spherically symmetric and isothermal:

(70) \begin{equation}\frac{\mathrm{d} \Phi}{\mathrm{d} R} = -\frac{\sigma_{\operatorname{i-bulge}}^2}{\rho_{\operatorname{i-bulge}}} \frac{\mathrm{d} \rho_{\operatorname{i-bulge}}}{\mathrm{d} R}.\end{equation}

Substituting equation (29) — and its derivative, where $\otimes \rightarrow \operatorname{i-bulge}$ — into (70) and rearranging gives

(71) \begin{equation}a_{\operatorname{i-bulge}} = 3 \left( \frac{1}{\sigma^2_{\operatorname{i-bulge}}} \frac{\mathrm{d} \Phi}{\mathrm{d} R} - \frac{1}{R} \right)^{-1} - R.\end{equation}

In principle, this expression should be solvable at any R. However, because of the approximation of isothermality (i.e. treating $\sigma_{\operatorname{i-bulge}}$ as radially invariant), equation (71) does not return identical answers at any two radii. For that reason, starting at small radii and working outwards incrementally, we continuously solve for $a_{\operatorname{i-bulge}}$ , averaging as we go, until reaching a converged answer.

8.1. Energy available from feedback

Let us begin with the assumption that the birth masses of stars, $\mathcal{M}_{\star}$ , for any population are distributed according to a Chabrier (Reference Chabrier2003) initial mass function (IMF, $\phi$ ), where we impose

(72a) \begin{equation}\phi(\bar{\mathcal{M}}_{\star}) = \left\{\begin{array}{l@{\quad}r}0, & \bar{\mathcal{M}}_{\star} < 0.1 \\[3pt]\dfrac{A_\star}{\bar{\mathcal{M}}_{\star}} \exp \left[-\dfrac{1}{2} \bigg(\dfrac{\log_{10} \big(\frac{\bar{\mathcal{M}}_{\star}}{0.08}\big)}{0.69}\bigg)^2 \right], & 0.1 \leq \bar{\mathcal{M}}_{\star} < 1\\[9pt]k_\star\, \bar{\mathcal{M}}_{\star}^{-2.3}, &1 \leq \bar{\mathcal{M}}_{\star} \leq 100 \\[3pt]0, & \bar{\mathcal{M}}_{\star} > 100\end{array}\right., \end{equation}

(72b) \begin{equation}\bar{\mathcal{M}}_{\star} \equiv \mathcal{M}_{\star} / \mathrm{M}_{\odot},\end{equation}

where $A_\star$ and $k_\star$ are constants. Recognising that (i) $\phi(\bar{\mathcal{M}}_{\star})$ must be continuous at $\bar{\mathcal{M}}_{\star} = 1$ and (ii)

(73) \begin{equation}\int_{0}^{\infty} \bar{\mathcal{M}}_{\star}\, \phi(\bar{\mathcal{M}}_{\star})\, \mathrm{d}\bar{\mathcal{M}}_{\star} = \int_{0.1}^{100} \bar{\mathcal{M}}_{\star}\, \phi(\bar{\mathcal{M}}_{\star})\, \mathrm{d}\bar{\mathcal{M}}_{\star} = 1\end{equation}

leads to the derivation that $k_\star \simeq 0.238$ and $A_\star \simeq 0.843$ .

To calculate the mass fraction of a stellar population returned to the ISM over a specified period of time, one must also know both the lifetimes ( $t_\mathrm{life}$ ) and end-of-life remnant masses ( $\mathcal{M}_\mathrm{rem}$ ) of stars as a function of their birth mass. For the former, we impose the textbook standard

(74) \begin{equation}\frac{t_\mathrm{life}}{\mathrm{Gyr}} = 10\,\bar{\mathcal{M}}_{\star}^{-5/2}\end{equation}

(Harwit Reference Harwit1988). This is readily derived from the assertions that the luminosity of a star scales as $\mathcal{M}_\star^{3.5}$ and that the Sun has an approximate lifetime of 10 Gyr.

For the remnant masses of stars, we roughly follow Madau & Dickinson (Reference Madau and Dickinson2014) in effecting

(75a) \begin{equation}\bar{\mathcal{M}}_\mathrm{rem}(\bar{\mathcal{M}}_{\star}) =\left\{\begin{array}{l@{\quad}r}\bar{\mathcal{M}}_{\star}, & \bar{\mathcal{M}}_{\star} < 1\\0.444 + 0.084\,\bar{\mathcal{M}}_{\star}, & 1 \leq \bar{\mathcal{M}}_{\star} \leq 7\\0.419 + 0.109\,\bar{\mathcal{M}}_{\star}, & 7 < \bar{\mathcal{M}}_{\star} < 8\\1.4, & 8 \leq \bar{\mathcal{M}}_{\star} \leq 50\\\bar{\mathcal{M}}_{\star}, & \bar{\mathcal{M}}_{\star} > 50\end{array}\right.\end{equation}

(75b) \begin{equation}\bar{\mathcal{M}}_\mathrm{rem} \equiv \mathcal{M}_\mathrm{rem} / \mathrm{M}_{\odot}.\end{equation}

The term for $\bar{\mathcal{M}}_\mathrm{rem}(7 < \bar{\mathcal{M}}_{\star} < 8)$ is something we have manually introduced as a means of linearly interpolating between the adjacent mass ranges, which more closely follow Madau & Dickinson (Reference Madau and Dickinson2014). The value of $\mathcal{M}_\mathrm{rem} = 1.4\,\mathrm{M}_{\odot}$ for $8 \leq \bar{\mathcal{M}}_{\star} \leq 50$ assumes that the remnant neutron star or BH has a mass close to the upper limit of white dwarfs (á la Chandrasekhar Reference Chandrasekhar1931). The basis for $\bar{\mathcal{M}}_\mathrm{rem} = \bar{\mathcal{M}}_{\star}$ when $\bar{\mathcal{M}}_{\star} > 50$ is that the collapse of the star is so rapid, its entire mass is committed to the remnant BH (which we still class as ‘stellar mass’ in the model). Given how quickly the IMF diminishes at high stellar masses (equation 72a), variations on the assumptions for $\bar{\mathcal{M}}_\mathrm{ rem}(\bar{\mathcal{M}}_{\star} \gtrsim 50)$ are inconsequential for our treatment of stellar evolution and feedback.

The mass returned to the ISMFootnote k for a given star equals $\mathcal{M}_\star - \mathcal{M}_\mathrm{rem}$ . We approximate that that mass is returned instantaneously at $\Delta t = t_\mathrm{life}$ after the star’s birth. We can therefore calculate the fraction of a stellar population’s initial mass that is returned to the ISM over a cosmic-time interval $[t_0-t_\mathrm{form}, t_1-t_\mathrm{form}]$ as

(76a) \begin{equation}f_{\star}^\mathrm{return} = \int_{\bar{\mathcal{M}}_1}^{\bar{\mathcal{M}}_0} \left[\bar{\mathcal{M}}_{\star} - \bar{\mathcal{M}}_\mathrm{rem}(\bar{\mathcal{M}}_{\star}) \right] \phi(\bar{\mathcal{M}}_{\star})\, \mathrm{d}\bar{\mathcal{M}}_{\star},\end{equation}

(76b) \begin{equation}\bar{\mathcal{M}}_x = \left( \frac{t_x - t_\mathrm{form}}{10\,\mathrm{Gyr}} \right)^{-2/5},\end{equation}

where the latter comes directly from equation (74), with $x = 0,1$ , while $t_\mathrm{form}$ denotes the time at which the stellar population formed. Note that because we define $t_1 > t_0$ and higher-mass stars die sooner, it follows that $\bar{\mathcal{M}}_1 < \bar{\mathcal{M}}_0$ , which is why $\bar{\mathcal{M}}_1$ is the lower bound for the integral. The fraction of the remaining mass of a population (where the returned mass from all stars with $\bar{\mathcal{M}}_{\star} > \bar{\mathcal{M}}_0$ has already been subtracted from the population’s total mass) that is returned to the ISM can likewise be calculated as

(77a) \begin{equation}f_\mathrm{rem}^\mathrm{return} = f_{\star}^\mathrm{return} \frac{m_\mathrm{*,pop}^\mathrm{form}}{m_\mathrm{*,pop}^\mathrm{rem}},\end{equation}

(77b) \begin{align}\frac{m_\mathrm{*,pop}^\mathrm{form}}{m_\mathrm{*,pop}^\mathrm{rem}} & = \bigg[ \int_{0.1}^{\bar{\mathcal{M}}_0} \bar{\mathcal{M}}_{\star}\, \phi(\bar{\mathcal{M}}_{\star})\, \mathrm{d}\bar{\mathcal{M}}_{\star}\nonumber\\[4pt]& \quad + \int_{\bar{\mathcal{M}}_0}^{100} \bar{\mathcal{M}}_\mathrm{rem}(\bar{\mathcal{M}}_{\star}) \, \phi(\bar{\mathcal{M}}_{\star})\, \mathrm{d}\bar{\mathcal{M}}_{\star} \bigg]^{-1},\end{align}

where $m_\mathrm{*,pop}^\mathrm{form}$ and $m_\mathrm{*,pop}^\mathrm{rem}$ represent the formation mass and remaining mass of the stellar population, respectively. This expression is practically the most useful in the design of Dark Sage, as we only explicitly track the remaining mass of a given population.

When mass is returned from a stellar population to the ISM, metals are returned to the ISM in proportion to the metallicity of that stellar population. In addition to this, new metals are added to the ISM, born from fusion in the stars’ cores and during supernova explosions. We maintain that the mass of new metals produced by a stellar population over its lifetime is a fixed fraction of the population’s birth mass (ignoring the metals already present). This fraction is known as the ‘yield’, which is expected to be $\sim$ 0.03, based on a Chabrier (Reference Chabrier2003) IMF combined with a Conroy, Gunn & White (Reference Conroy, Gunn and White2009) simple stellar population model (Lagos et al. Reference Lagos2018). These new metals are deposited into the ISM at a rate proportional to the mass-loss rate of stars in that population. This occurs before any reheating from feedback (Section 8.2). Note that we do not track the individual abundances of metal elements, only their summed content.

Let us hypothesise for a moment that we know the average energy released by a supernova that is used in heating gas out of the ISM and/or ejecting it out of the CGM, which we label $\mathcal{E}_\mathrm{SN}$ (note that we make no assumptions about the total energy released per supernova). Then we must ask: how many supernovae will there be for a given stellar population and when will they explode? Equation (74) already provides us the answer to the latter — any star that produces a supernova does so at the end of its life. For the former, we assume all stars with $8 \leq \bar{\mathcal{M}}_{\star} \leq 50$ produce a (core-collapse/Type-II) supernova. Additionally, we assume some fraction of stars with $1 \leq \bar{\mathcal{M}}_{\star} \leq 8$ produce a (accretion-driven/Type-Ia) supernova. This means we effectively assume that each star with $1 \leq \bar{\mathcal{M}}_{\star} \leq 8$ has a fixed probability, $p_{\operatorname{SN-Ia}}$ , of producing a supernova at the end of its life. The total number of supernovae of each type per unit mass of stars formed in a population, $n_\mathrm{SN}$ , is thus

(78a) \begin{equation}\frac{n_{\operatorname{SN-II}}}{\mathrm{M}_{\odot}^{-1}} = \int_8^{50} \phi(\bar{\mathcal{M}}_{\star})\, \mathrm{d}\bar{\mathcal{M}}_{\star},\end{equation}

(78b) \begin{equation}\frac{n_{\operatorname{SN-Ia}}}{\mathrm{M}_{\odot}^{-1}} = p_{\operatorname{SN-Ia}} \int_1^8 \phi(\bar{\mathcal{M}}_{\star})\, \mathrm{d}\bar{\mathcal{M}}_{\star}.\end{equation}

To determine $p_{\operatorname{SN-Ia}}$ , we look to observations. Broadly speaking, the rate of Type-II supernovae are observed to exceed that of Type Ia. For example, Tsujimoto et al. (Reference Tsujimoto, Nomoto, Yoshii, Hashimoto, Yanagida and Thielemann1995) note that in the Milky Way, the ratio of Type-Ia to Type-II supernovae — equivalent to $n_{\operatorname{SN-Ia}} / n_{\operatorname{SN-II}}$ — is $\sim$ 0.15, but that it is more like 0.2–0.3 in the Magellanic Clouds. By contrast, Maggi et al. (Reference Maggi2016) find $n_{\operatorname{SN-Ia}} / n_{\operatorname{SN-II}} \simeq 0.74$ in the Large Magellanic Cloud (LMC), but they expect this to be above the true average ratio due to the short timescale that their data probe (order $10^4$ yr) and the shape of the LMC’s star formation history. In the absence of tighter observational constraints, we choose to adopt the Milky Way value of $n_{\operatorname{SN-Ia}} / n_{\operatorname{SN-II}} = 0.15$ . Through equation (78), we can therefore solve for $p_{\operatorname{SN-Ia}} \simeq 0.01$ . In line with equation (76), we can calculate the total energy imparted on the local ISM by a stellar population over a given time interval, which we frame as a birth-mass interval, as

(79a) \begin{equation}\Delta E_\mathrm{SN} = \mathcal{E}_\mathrm{SN}\, m_\mathrm{*,pop}^\mathrm{form} \int_{ \bar{\mathcal{M}}_1}^{ \bar{\mathcal{M}}_0} \varphi(\bar{\mathcal{M}}_{\star}) \, \mathrm{ d}\bar{\mathcal{M}}_{\star},\end{equation}

(79b) \begin{equation}\varphi(\bar{\mathcal{M}}_{\star}) \equiv\left\{\begin{array}{l@{\quad}r}p_{\operatorname{SN-Ia}}\, \phi(\bar{\mathcal{M}}_{\star}), & \bar{\mathcal{M}}_{\star} < 8 \\\phi(\bar{\mathcal{M}}_{\star}), & \bar{\mathcal{M}}_{\star} \geq 8 \\\end{array}\right..\end{equation}

We show how our mass-return and supernova models look as a function of star birth mass and time since population birth in Figure 7.

Figure 7. Cumulative mass fraction of a stellar population returned to the ISM (top, equation 76) and cumulative number of supernovae per unit initial mass of that population (bottom) as a function of star birth mass (left) and time since the population’s formation (right) implemented in our new stellar evolution/feedback model in Dark Sage. For reference, the vertical lines in the right-hand panels represent the age bins used in Dark Sage for this work (treating $t_\mathrm{form}$ as 0; see Section 3.2). For stellar populations born at $t_\mathrm{form} > 0$ , the vertical lines would effectively be shifted to the left to describe how the evolution of that population is discretised in Dark Sage.

All stars are assumed to form at the middle of an age bin. Returned mass, the production of new metals, and feedback from stars born within the same age bin occur instantaneously, always assuming half the age bin has elapsed, regardless of where in that bin the current sub-time-step actually is. This is necessary, as once the model has moved into the time window of the next age bin, all stars formed in the previous age bin are assumed to be part of the same population with a single $t_\mathrm{form}$ . This is the nature of discretised time; regardless of how many bins one has, one always has to treat the ‘current’ bin in a special way. This also means that if we were to set the number of age bins to 1, we would return to the classical instantaneous recycling approximation (with the caveat that the recycling fraction is solved for, rather than given as an explicit input). Whenever star formation occurs, we check if it is necessary to rescale the local SFR to ensure that the mass consumed in star formation plus that associated with the instantaneous portion of stellar feedback does not exceed the mass of the gas annulus (in line with equation 64).

Several numbers in our model for stellar feedback could have been treated as free parameters. For example, the ratio of Type-Ia to Type-II supernovae for a population of stars, the parameters controlling the shape of the IMF, and the normalisation and power index of stellar life-times as a function of birth mass. Instead, we have opted to fix those numbers, and allow one stellar feedback parameter to vary freely in our calibration: $\mathcal{E}_\mathrm{SN}$ . Our rationale is that variation in most of the other parameters would be degenerate with changes to $\mathcal{E}_\mathrm{SN}$ to first order for calibration purposes. The exact value of $\mathcal{E}_\mathrm{SN}$ should, therefore, not be interpreted too literally, as it compensates for uncertainty in numerous assumptions we have made. Nevertheless, we should not expect $\mathcal{E}_\mathrm{SN}$ to depart too far from the canonical value of $10^{44}$ J (for a discussion on the energy output of supernovae, see e.g. Smartt Reference Smartt2009). Sure enough, as we show in Section 13, we find the best value of $\mathcal{E}_\mathrm{SN}$ to be ${\sim}0.8 \times 10^{44}$ J.

There are further details one may wish to consider in a stellar feedback model, such as how the metallicity of a stellar population affects the rate of supernovae (e.g. Kistler et al. Reference Kistler, Stanek, Kochanek, Prieto and Thompson2013), how the energetics of Type-Ia and Type-II supernovae differ, and how the IMF may change with metallicity and formation time (e.g. Li et al. Reference Li, Liu, Zhang, Tian, Fu, Li and Yan2023). Implicitly, our ignorance of these details is also folded in to the best-fitting value of $\mathcal{E}_\mathrm{SN}$ .

8.2. Impartation of energy on the ISM

With the energy from stellar feedback determined, the next step is to decide how that energy is used. The standard practice for doing this is to reheat an appropriate mass of gas out of the ISM and/or eject it. There are three major differences between how we do this now versus previous versions of Dark Sage:

1. 1. We do not impose a mass loading factor. Instead, we solve for the amount of gas reheated out of the ISM using energy-conservation arguments, as we outline below.

2. 2. We split this reheated mass into two categories: the fraction with enough energy to escape the halo’s virial radius, and that which will remain bound to the CGM.

3. 3. We do not instantaneously transfer mass to the hot and/or ejected reservoirs, but instead initially transfer them to the new fountain and outflowing reservoirs, thereby enforcing a delay before that gas is available for cooling.

Let us imagine that a supernova outflow begins carrying the mass prescribed by its return fraction of the source stellar population according to its age and the IMF assumptions laid out in Section 8.1. This initial blast then collides with the ISM, picking up mass and driving a wind that removes that mass from the ISM. Once mixed, the gas must also be heated to the temperature of the CGM, which is assumed to be uniform (equation 8). Assuming this wind conserves energy, we have

(80) \begin{equation}\Delta E_\mathrm{SN}^{(i)} = \frac{1}{2}\, \Delta m_\mathrm{cold \rightarrow CGM}^{(i) \rm SN} \left(V_\mathrm{200c}^2 + v_\mathrm{wind}^{(i)\,2}\right).\end{equation}

Here, we have defined $v_\mathrm{wind}^{(i)}$ in the annulus’s rotating frame as the equivalent initial wind velocity were the outflow were to hypothetically reach thermal equilibrium with the CGM immediately. In reality, we would expect the wind to start colder and faster, gradually heating and slowing (i.e. exchanging kinetic energy for thermal energy in addition to gravitational potential energy) before fully mixing with the CGM. Defining $v_\mathrm{wind}$ as we have simply provides a mathematical shortcut to reach our desired outcome (below). While we do model how long it takes an outflow to incorporate into the CGM (or travel beyond; see Section 10), this timescale is not informed by $v_\mathrm{wind}$ .

Dark Sage is designed such that any feedback-affected gas must have its energy raised to at least that of the hot-gas reservoir. As such, we can solve for $v_\mathrm{wind}^{(i)}$ under the condition that

\begin{align*}\Delta E_\mathrm{SN}^{(i)} \rightarrow \left(e_\mathrm{hot} - e_\mathrm{cold}^{(i)}\right) \Delta m_\mathrm{cold \rightarrow CGM}^{(i) \rm SN}, \end{align*}

where each e represents the specific energy of the subscripted reservoir. Combining with equation (80) and rearranging gives

(81) \begin{equation}v_\mathrm{wind}^{(i)} = \sqrt{2 \left(e_\mathrm{fount} - e_\mathrm{cold}^{(i)}\right) - V_\mathrm{200c}^2}.\end{equation}

The specific energy of the hot-gas reservoir can be broken into three components: its specific potential energy, specific thermal energy, and specific kinetic energy (i.e. that which is non-thermal, as noted to be non-negligible by e.g. Lochhaas et al. Reference Lochhaas2021), that is,

(82) \begin{equation}e_\mathrm{hot} = \langle \Phi \rangle_\mathrm{hot} + e_\mathrm{hot}^\mathrm{kin} + e_\mathrm{hot}^\mathrm{therm}.\end{equation}

To calculate the average potential of gas in the CGM, we numerically solve

(83) \begin{equation}\langle \Phi \rangle_\mathrm{hot} = \frac{4 \pi}{m_\mathrm{hot}} \int_0^{R_\mathrm{200c}} \rho_\mathrm{hot}(R)\, \Phi(R)\, R^2\, \mathrm{d}R.\end{equation}

The potential profile, $\Phi(R)$ , is found by numerically integrating equation (23).Footnote l Note that $\Phi(r)$ is defined in the plane of the disc in equation (23). As a necessary approximation for computational efficiency, we assume this functional form works for three-dimensional radii, R, as well. In the limit where the ratio of the disc mass to halo mass goes to zero, this approximation is accurate. However, in general, it is imperfect.

The specific kinetic energy of hot gas is informed entirely by bulk rotational motion. We approximate its bulk rotational velocity based on its specific angular momentum (assumed here to be equal to that of the halo) and the average radius of a gas element within it:

(84a) \begin{equation}e_\mathrm{hot}^\mathrm{kin} = \frac{j_\mathrm{hot}^2}{2\,\langle R \rangle^2_\mathrm{hot}} = \frac{j_\mathrm{halo}^2}{2\,\langle R \rangle^2_\mathrm{hot}}, \end{equation}

(84b) \begin{equation}\langle R \rangle_\mathrm{hot} = \frac{4\pi}{m_\mathrm{hot}} \int_0^{R_\mathrm{200c}} \rho_\mathrm{hot}(R)\, R^3\, \mathrm{d}R \end{equation}

(84c) \begin{align} = \frac{1-\ln(2)}{2}\, R_\mathrm{200c}\, c_\beta^2(z)\, \mathcal{C}_\beta(z).\end{align}

As we assume $T_\mathrm{hot} = T_\mathrm{vir}$ throughout the CGM, the thermal specific energy of the hot-gas reservoir is readily found as

(85) \begin{equation}e_\mathrm{hot}^\mathrm{therm} = \frac{V_\mathrm{200c}^2}{2}. \end{equation}

For satellites, we use $V_\mathrm{200c}$ at infall here, in line with our assumption that $T_\mathrm{hot}$ remains fixed after infall from Section 4.

Similar to equation (82), the specific energy of an annulus in the ISM can be broken down as

(86) \begin{equation}e_\mathrm{cold}^{(i)} = \left\langle \Phi^{(i)} \right\rangle + e_\mathrm{cold}^{(i)\rm kin}.\end{equation}

No thermal term is added for the ISM; we have assumed that the temperature of the ISM is negligible compared to that of the CGM, such that

\begin{align*}T_\mathrm{CGM} - T_\mathrm{cold} \simeq T_\mathrm{vir}.\end{align*}

The potential specific energy of an annulus follows equation (36), while its specific kinetic energy is readily found through contributions from its circular velocity and average vertical motion:

(87) \begin{equation}e_\mathrm{cold}^{(i)\rm kin} = \frac{\big\langle v_\mathrm{circ}^{(i)} \big\rangle ^2 + \sigma_\mathrm{cold}^2(z)}{2}.\end{equation}

While the above provides the necessary equations to calculate $\Delta m_\mathrm{cold \rightarrow CGM}^{(i) \rm SN}$ , we are yet to specify how this feedback-affected gas is handled in detail. As noted above, we distribute this gas into two transitory reservoirs: the fountain reservoir and the outflowing reservoir, meaning

(88) \begin{equation}\Delta m_\mathrm{cold \rightarrow CGM}^{(i) \rm SN} = \Delta m_\mathrm{cold \rightarrow fount}^{(i) \rm SN} + \Delta m_\mathrm{cold \rightarrow outfl}^{(i) \rm SN}.\end{equation}

The fountain reservoir is defined to have the same specific energy as the hot-gas reservoir:

(89) \begin{equation}e_\mathrm{fount} = e_\mathrm{hot}.\end{equation}

We track a timescale for which this gas is reincorporated in the hot-gas reservoir, which we describe in Section 10.1. The outflowing reservoir represents gas with energy to escape the (sub)halo’s virial radius, and is gradually transferred to the ejected reservoir (see Section 10.2). How, then, do we decide the fraction of $\Delta m_\mathrm{cold \rightarrow CGM}^{(i) \rm SN}$ that goes to each of the fountain and outflowing reservoirs?

Consider the ideal that any gas kicked from the ISM with a velocity greater than the halo escape velocity,

(90) \begin{equation}\langle v_\mathrm{esc}^{(i)}\rangle ^2 = 2 \left[ \Phi(R_\mathrm{200c}) - \left\langle \Phi^{(i)} \right\rangle \right], \end{equation}

should be placed in the outflowing reservoir, and that the remainder goes to the fountain reservoir. For an observer external to the galaxy — that is, in the inertial frame of the halo — the appropriate velocity to compare against a given $\langle v_\mathrm{ esc}^{(i)}\rangle$ is not simply $v_\mathrm{wind}^{(i)}$ , but rather the vector sum of it and the local rotation velocity (recall that $v_\mathrm{wind}$ is defined to have already taken care of the fact that the wind slows from heating, which is what allows us to compare it to the escape velocity in the first place). We assume that the energy from a supernova is emitted with spherical symmetry in the rotating frame of the annulus. Any gas on the leading side will therefore be kicked with a heightened speed, and conversely that on the trailing side to a lower speed. One can solve for the angle between these two vectors, where the magnitude of the sum is equal to the escape velocity, through the cosine rule. The fraction of gas expelled from the ISM with a velocity above the escape velocity of the halo is then equivalent to the fraction of possible angles less than this critical angle. Recognising that the range of possible angles follows a $\sin(\theta)$ distribution — where $\theta \in [0,\pi]$ — then provides the ratio of ejected gas to reheated gas (after some omitted mathematics):

(91) \begin{equation} \frac{\Delta m_\mathrm{cold \rightarrow fount}^{(i) \rm SN} }{\Delta m_\mathrm{cold \rightarrow outfl}^{(i) \rm SN} } = \left( \frac{1}{2} - \frac{ \left\langle v_\mathrm{esc}^{(i)}\right\rangle ^2 - \big\langle v_\mathrm{circ}^{(i)} \big\rangle ^2 - v_\mathrm{wind}^{(i)\,2} }{4\, \big\langle v_\mathrm{circ}^{(i)} \big\rangle \, v_\mathrm{wind}^{(i)}} \right)^{-1} -1 .\end{equation}

Built into equation (91) is the implicit assumption that

\begin{align*}\big\langle v_\mathrm{circ}^{(i)} \big\rangle - v_\mathrm{wind}^{(i)} < \left\langle v_\mathrm{esc}^{(i)}\right\rangle < \big\langle v_\mathrm{circ}^{(i)} \big\rangle + v_\mathrm{wind}^{(i)} .\end{align*}

However, there is freedom in the solution for $v_\mathrm{wind}^{(i)}$ to not obey this restricted range. If $\big\langle v_\mathrm{circ}^{(i)} \big\rangle + v_\mathrm{wind}^{(i)} < \langle v_\mathrm{esc}^{(i)}\rangle$ , then none of the gas kicked out of the ISM is capable of escaping the halo, meaning we set $ \Delta m_\mathrm{cold \rightarrow outfl}^{(i) \rm SN} = 0$ . For the inverse reason, if $\big\langle v_\mathrm{circ}^{(i)} \big\rangle - v_\mathrm{wind}^{(i)} > \langle v_\mathrm{esc}^{(i)}\rangle$ , we set $ \Delta m_\mathrm{cold \rightarrow fount}^{(i) \rm SN} = 0$ .

In our framework, energy conservation dictates that

(92) \begin{align}\Delta E_\mathrm{SN}^{(i)} = \Delta m_\mathrm{cold \rightarrow fount}^{(i) \rm SN} \!\left( e_\mathrm{fount} - e_\mathrm{cold}^{(i)} \right) + \Delta m_\mathrm{cold \rightarrow outfl}^{(i) \rm SN} \!\left( e_\mathrm{outfl}^{(i) \rm SN} - e_\mathrm{cold}^{(i)} \right) .\end{align}

The final quantity left to define is $e_\mathrm{outfl}^{(i) \rm SN}$ , which is the specific energy of the material presently being added to the outflow reservoir from a given annulus. In truth, this mass has a distribution of specific energies, for which we know the minimum and maximum, based on the picture presented above:

\begin{align*}& e_\mathrm{outfl}^{(i) \rm SN} \geq \Phi(R_\mathrm{200c}) + e_\mathrm{hot}^\mathrm{therm},\\[5pt]& e_\mathrm{outfl}^{(i) \rm SN} \leq \left\langle \Phi^{(i)} \right\rangle + e_\mathrm{hot}^\mathrm{therm} + \frac{1}{2} \left[ \left( \big\langle v_\mathrm{circ}^{(i)} \big\rangle + v_\mathrm{ wind}^{(i)} \right)^2 + \sigma_\mathrm{gas}^2 \right].\end{align*}

We approximate $e_\mathrm{outfl}^{(i) \rm SN}$ as the average of these two limits:

(93) \begin{align}e_\mathrm{outfl}^{(i) \rm SN} = e_\mathrm{hot}^\mathrm{therm} + \frac{ \Phi(R_\mathrm{200c}) + \left\langle \Phi^{(i)} \right\rangle }{2} + \frac{1}{4} \left[ \left( \big\langle v_\mathrm{circ}^{(i)} \big\rangle + v_\mathrm{wind}^{(i)} \right)^2 + \sigma_\mathrm{gas}^2 \right].\nonumber\\[3pt]\end{align}

Unlike the other reservoirs, the specific energy of the outflow reservoir is updated as mass is added to it:

(94) \begin{equation}e_\mathrm{outfl}(t + \Delta t) = \mathrm{max} \Bigg[ \Phi(R_\mathrm{200c}) + e_\mathrm{hot}^\mathrm{kin} + e_\mathrm{hot}^\mathrm{therm}, \\\frac{ m_\mathrm{outfl}(t)\, e_\mathrm{outfl}(t) + \sum_{i=0}^{N_\mathrm{ann}} \Delta m_\mathrm{cold \rightarrow outfl}^{(i) \rm SN}\, e_\mathrm{outfl}^{(i) \rm SN} }{ m_\mathrm{ outfl}(t) +\sum_{i=0}^{N_\mathrm{ann}} \Delta m_\mathrm{cold \rightarrow outfl}^{(i) \rm SN}} \Bigg].\end{equation}

The first term inside the ‘max’ function ensures the outflowing energy exceeds the fountaining energy, which it must by definition. It must at least carry enough energy to reach $R_\mathrm{200c}$ and have the equivalent of the thermal and kinetic specific energy of the hot gas. $e_\mathrm{outfl}$ is needed in other parts of the model, which we come to later in the paper.

We come back to what happens to the feedback-affected gas next in Section 10.

8.3. Mass loading factor

Our new stellar feedback model means we now have an entirely predictive framework for determining the mass loading factors of galaxies experiencing feedback. The mass loading factor of a galaxy is defined as the ratio of its gas reheating rate from stellar feedback to its star formation rate:

(95) \begin{equation}\eta \equiv \frac{\dot{m}_\mathrm{cold \rightarrow fount}^\mathrm{SN} + \dot{m}_\mathrm{cold \rightarrow outfl}^\mathrm{SN}}{\mathrm{SFR}}, \end{equation}

where this excludes any consideration of gas heating from quasar feedback, but includes all channels of star formation. Here, $\dot{m}$ and SFR refer to the mass transfer rates and star formation rates averaged over the previous full time step. This is physically more informative than the most recent sub-time-step, as it is less susceptible to the order in which we resolve galaxy evolution processes (this is why the sub-time-steps are there in the first place).

Dark Sage is one of few SAMs in the literature to make a true prediction for the mass loading factor from stellar feedback (for another example, see Lagos, Lacey & Baugh Reference Lagos2013). Comparing Dark Sage’s mass loading predictions with observations and other models is not trivial; definitions of mass loading factors in the literature are far from uniform. This is in part because both real and hydrodynamic-simulation galaxies are not so easily broken into discrete reservoirs as SAMs define galaxies and haloes to have. For the former, it is typically necessary to measure the outflow rate (the numerator in equation 95) at some radius. Our definition of $\eta$ is effectively measured at a radius of zero. This means our definition of $\eta$ should be systematically larger than that calculated from observational studies (and hydrodynamic simulations).

While we make a face value comparison with some observational data from the literature in the left panel of Figure 8, for the above reason and others (see discussion in Nelson et al. Reference Nelson2019a on this topic), we caution against over-interpreting the apparent offset Dark Sage has with the data. Points for individual galaxies at $z \simeq 0$ are shown from Chisholm et al. (Reference Chisholm, Tremonti, Leitherer and Chen2017), as they are for $z < 0.2$ galaxies from Heckman et al. (Reference Heckman, Alexandroff, Borthakur, Overzier and Leitherer2015), while the points from Sugahara et al. (Reference Sugahara2017) are bins of galaxies from their $z \simeq 0$ sample.

Figure 8. Mass loading factor predicted by Dark Sage as a function of stellar mass (left) and maximal circular velocity (right-hand panel) for central galaxies with $\mathrm{sSFR} \geq 10^{-11}\,\mathrm{yr}^{-1}$ at $z = 0$ . Greyscale hexbins represent a two-dimensional histogram of Dark Sage galaxies, where the solid, purple lines follow the running median (thick) and 16th and 84 percentiles (thin). Compared are observations in the left-hand panel (Heckman et al. Reference Heckman, Alexandroff, Borthakur, Overzier and Leitherer2015; Chisholm et al. Reference Chisholm, Tremonti, Leitherer and Chen2017; Sugahara et al. Reference Sugahara2017). Based on the captions of tables 1 and 2 of Heckman et al. (Reference Heckman, Alexandroff, Borthakur, Overzier and Leitherer2015), we assume an error of 0.3 dex in stellar mass and $\sqrt{0.2^2 + 0.25^2}$ dex in $\eta$ . Various other models are compared in the right-hand panel (including Guo et al. Reference Guo2011; Bower, Benson & Crain Reference Bower, Benson and Crain2012; Croton et al. Reference Croton2016; Lagos et al. Reference Lagos2018). Note that $V_\mathrm{max}$ is from subfind, which is different to what we derive from the rotation curves we build; this is the fairest quantity to compare to other semi-analytic models.

In the right-hand panel of Figure 8, we compare Dark Sage’s predicted mass loading factors as a function of maximum halo circular velocity against the input functional mass loading factors of other SAMs.Footnote m The assumed functions used in Bower et al. (Reference Bower, Benson and Crain2012) and Lagos et al. (Reference Lagos2018) are similar in slope to what Dark Sage predicts, but systematically higher and lower, respectively. Meanwhile, the asymptotic function of Guo et al. (Reference Guo2011) and the constant used by Croton et al. (Reference Croton2016) appear both qualitatively and quantitatively inconsistent with our energy-conserving framework. The comparison with Lagos et al. (Reference Lagos2018) is particularly interesting, as their fitting function was based on the predictive model of Lagos et al. (Reference Lagos2013). There are many differences between Dark Sage and the Lagos et al. (Reference Lagos2013) version of galform, placing an elaborative discussion and comparison between the models outside the scope of this paper.

The slope and scatter of $\eta$ as a function of either $m_*$ or $V_\mathrm{max}$ are the outcome of our complex network of equations in Dark Sage, although they ultimately boil down to the diversity of gravitational potentials at fixed mass and the increase in average potential depth at higher masses. The normalisation of $\eta$ in Dark Sage is effectively set by the value of $\mathcal{E}_\mathrm{SN}$ . As we show in Section 13, this is calibrated to other observations in a way that is blind to Figure 8.

8.4. Metal loading

In a departure from previous versions of Dark Sage, we do not assume that gas reheated from the ISM carries the same metallicity as the ISM. We imagine a scenario where the energy from a stellar feedback event expands with spherical symmetry (in the ISM’s frame) into an ISM where all atoms (metals and non-metals) are distributed randomly and are isothermal. If the cross section of each atom were the same, we would expect the same amount of energy to be imparted on each atom on average. While heavier elements may have a greater cross section for interacting with an outflow when stationary, this is offset by the lighter elements having higher thermal motion at fixed temperature. In the absence of a detailed model for these competing effects, we approximate the energy imparted on each atom to be equal. Because we have implicitly assumed elsewhere that anything reheated carries a common specific energy, this approximation must mean some metals are left behind. As oxygen is the most abundant metal in the Universe, and carbon the second-most abundant, the average atomic mass of a metal should be a little less than $16\,m_\mathrm{H}$ . Non-metals (hydrogen and helium) should have an average atomic mass of around $7\,m_\mathrm{H}/4$ . This implies that for every unit mass of non-metals reheated from an annulus, we should see $\sim$ $Z^{(i)}_\mathrm{cold}/9$ times that mass of metals reheated. Or, in other words, the metallicity of the locally produced outflow is one tenth that of the local ISM. This is what we enforce.

One potential issue with our treatment of metals in stellar feedback outflows is that some observations find evidence that the metallicity of outflows is typically elevated relative to the ISM (e.g. Chisholm, Tremonti & Leitherer Reference Chisholm, Tremonti and Leitherer2018), which is the opposite of what we have imposed. However, there is once again more nuance than this in comparing our model to observations. As we explain more in Section 10, the outflow that we follow from a galaxy in Dark Sage represents that gas that will eventually either redistribute itself within the CGM, with the same energy as the CGM, and/or be ejected from the halo. We essentially consider all the material in an outflow to carry one of two possible energy states. This means we lose the information of an energy or velocity distribution in the outflow. In practice, we should expect the low-energy portion of that distribution to never fully reincorporate into the quasi-hydrostatic CGM, but instead to return to the ISM whence it came on a relatively short timescale. One might also hypothesise that an outflow has a metallicity distribution, which is anti-correlated with the velocity distribution, by virtue of metals being heavier. It is therefore conceivable that the metallicity of an outflow observed close to the disc of a galaxy might differ significantly from the metallicity of the portion of that outflow that ends up mixing with the CGM. Moreover, the statistics on the metallicity of outflows are weak, and the measurement uncertainties are often dominated by systematics.

One observed relation that is much more enshrined in canon is that between the stellar mass and gas metallicity of galaxies. In previous versions of the model, this relation informed the calibration of the model’s free parameters. As we discuss in Section 13, our calibration method has changed significantly, and the mass–metallicity relation of galaxies is no longer used. This means it is now a prediction of the model. We show this relation for Dark Sage in Figure 9. This is compared to the observation-based relations of Tremonti et al. (Reference Tremonti2004) and Bellstedt et al. (Reference Bellstedt2021).

Figure 9. The relation between stellar mass and ISM metallicity for galaxies at $z = 0$ . We only include Dark Sage galaxies with $m_\mathrm{neutral} \geq 10^8\,\mathrm{M}_{\odot}$ to ensure there is appreciable gas present for a meaningful gas metallicity. Dark Sage data are plotted in the same way as Figure 8. Compared is the observed relation from Tremonti et al. (Reference Tremonti2004) at $z \simeq 0.1$ , covering their 16–84th interpercentile range, and the running median from Bellstedt et al. (Reference Bellstedt2021) at $z = 0.07$ .

We note that there are hidden systematic uncertainties in the compared mass–metallicity relations in Figure 9, which should not make the apparent systematic offset between them and Dark Sage a point of concern. For Bellstedt et al. (Reference Bellstedt2021), assumptions in their method of modelling galaxies’ spectral energy distributions are a natural source of systematic uncertainty, such as their enforced proportionality in the evolution of $Z_\mathrm{gas}$ with stellar-mass growth. There is also $\sim$ 1 Gyr more of cosmic evolution from their $z = 0.07$ result to our Dark Sage result at $z = 0$ . For the comparison with Tremonti et al. (Reference Tremonti2004), whose data represent $z \simeq 0.1$ , the largest sources of systematic uncertainty are from the emission line diagnostics used to measure the oxygen abundance (see Kewley & Ellison Reference Kewley and Ellison2008) and the conversion from oxygen abundance to metallicity. In our comparison, we assume

(96) \begin{equation}12 + \log_{10}(\mathrm{O} / \mathrm{H}) = 8.69 + \log_{10}(Z_\mathrm{gas} / 0.0142),\end{equation}

where 0.0142 is the estimated metallicity of the Sun (Asplund et al. Reference Asplund, Grevesse, Sauval and Scott2009) (we have corrected the Tremonti et al. Reference Tremonti2004 relation for a Chabrier Reference Chabrier2003 IMF and our cosmology). For similar reasons, in their comparison of the mass–metallicity relation of the IllustrisTNG simulation with observations, Torrey et al. (Reference Torrey2019) ignore the normalisation of the Tremonti et al. (Reference Tremonti2004) relation, moving it down in their figure 6. While we could have taken the same approach in our Figure 9, we opted not to for the sake of fidelity.

The change to how we model metal loading in Dark Sage turned out to be crucial to the outcome shown in Figure 9. That is, had the outflows from the ISM carried higher metallicity, $Z_\mathrm{cold}$ would have been systematically lower in Dark Sage galaxies. While we outlined above why some systematic offset between Dark Sage and observations as presented in Figure 9 is expected — in this case, an apparent factor of $<$ 1.4 — larger offsets naturally become increasingly harder to waive away.

With the resolved detail of Dark Sage discs, we can also predict the gas metallicity gradients of galaxies. This is a rare quality for a SAMs, which few others have (as another example, see Yates et al. Reference Yates, Henriques, Fu, Kauffmann, Thomas, Guo, White and Schady2021). As a first prediction from the new Dark Sage, we show the radius-weighted average metallicity gradient of Dark Sage galaxies as a function of stellar mass in Figure 10, where

Figure 10. Radius-weighted average of the metallicity gradient in gas discs in Dark Sage galaxies at $z = 0$ . Data for Dark Sage are presented in the same way as in previous figures. Dotted lines with error bars represent the average relations from three IFS surveys: the ‘Calar Alto Legacy Integral Field spectroscopy Area’ survey (CALIFA; Sánchez et al. Reference Sánchez2014), the ‘Mapping nearby Galaxies at Apache Point Observatory’ survey (MaNGA; Belfiore et al. Reference Belfiore2017), and the ‘Sydney-AAO Multi-object Integral-field spectrograph’ survey (SAMI; Poetrodjojo et al. Reference Poetrodjojo2021), which have been homogenized (Sharda et al. Reference Sharda, Krumholz, Wisnioski, Archaryya, Federrath and Forbes2021; Acharyya et al. in prep.).

(97) \begin{align}\left\langle \frac{ \mathrm{d} \log_{10} (Z_\mathrm{cold})}{ \mathrm{d} r} \right\rangle & \equiv \left[ \sum_{i=1,\, m_\mathrm{cold}^{(i)} \neq 0}^{N_\mathrm{ann}-1} \log_{10} \left(Z_\mathrm{cold}^{(i)}\right) \right] \nonumber\\[5pt]& \quad \times \left[ \sum_{i=1,\, m_\mathrm{cold}^{(i)} \neq 0}^{N_\mathrm{ann}-1} \left(\langle r^{(i+1)} \rangle - \langle r^{i} \rangle \right) \right]^{-1}.\end{align}

We choose to explicitly weight by the radial range of each annulus, as observational gradients are typically implicitly weighted in the same way, that is, data will be binned by radius and a gradient fitted to those bins. Such observations from integral field spectroscopy (IFS) surveys are compared in Figure 10 (to appear in Acharyya et al. in prep.),Footnote n though we caution that this is still a cursory comparison, and should therefore not be interpreted as a means of falsifying Dark Sage. While the vast majority of Dark Sage galaxies have negative metallicity gradients (which is the expected norm), Dark Sage is also able to produce galaxies with positive metallicity gradients. The relation between gas metallicity gradient and stellar mass is relatively flat in Dark Sage, as it is in observations. But this is dependent on how the metallicity gradient is calculated. Though not shown here, reproducing Figure 10 with a mass-weighted average gradient across the annuli in a Dark Sage disc results in an increasingly steep relation towards low stellar masses. Metallicity gradients normalised by the stellar disc half-mass radius look very similar (but for a re-scaled y-axis).

8.5. Direct impartation of energy on the CGM

The instability-driven bulge, merger-driven bulge, and intrahalo stars (collectively ‘spheroid stars’) can produce feedback. While no stars form in spheroids in Dark Sage, the delayed feedback scheme means that spheroid stars still produce Type-Ia supernovae. We assume that spheroid supernovae cannot reheat gas from the disc. Instead, the metals and energy they produce are imparted directly onto the CGM. Summing over all stellar-age bins, the energy from spheroid supernovae, $\Delta E_\mathrm{SN}^\mathrm{sph}$ , is balanced by generating an outflow from the CGM. In practice, this means transferring mass from the hot and fountain reservoirs to the outflowing reservoir. Because we have defined $e_\mathrm{fount} = e_\mathrm{hot}$ , we transfer mass proportionally out of the hot and fountain reservoirs, assuming this gas carries the same specific energy as that pre-computed for the outflowing reservoir:

(98a) \begin{equation}\Delta m_\mathrm{hot \rightarrow outfl}^\mathrm{SN} = \mathrm{min} \left[m_\mathrm{hot}, \,\frac{\Delta E_\mathrm{SN}^\mathrm{sph}}{e_\mathrm{outfl} - e_\mathrm{hot}} \frac{m_\mathrm{hot}}{m_\mathrm{hot} + m_\mathrm{fount}} \right],\end{equation}

(98b) \begin{equation}\Delta m_\mathrm{fount \rightarrow outfl}^\mathrm{SN} = \mathrm{min} \left[m_\mathrm{fount}, \,\frac{\Delta E_\mathrm{SN}^\mathrm{sph}}{e_\mathrm{outfl} - e_\mathrm{hot}} \frac{m_\mathrm{fount}}{m_\mathrm{hot} + m_\mathrm{fount}} \right].\end{equation}

8.6. Comparative remarks

Although vastly different in how the feedback is applied, a clear parallel can be drawn between our model for delayed stellar feedback and that in the meraxes SAM (Mutch et al. Reference Mutch, Geil, Poole, Angel, Duffy, Mesinger and Wyithe2016). This is not entirely by coincidence; part of the inspiration for our model design came from Mutch et al. (Reference Mutch, Geil, Poole, Angel, Duffy, Mesinger and Wyithe2016). Similar designs of delayed feedback and chemical enrichment are also incorporated into the galacticus model (Benson Reference Benson2012) and in De Lucia et al. (Reference De Lucia, Tornatore, Frenk, Helmi, Navarro and White2014), although again different in implementation to Dark Sage. We refer the reader to the other papers for more context.

In previous incarnations of SAGE and Dark Sage, the returning of mass from stars to the ISM and the production of new metals was done after reheating gas from feedback. Implicitly, this meant it was guaranteed that some (local, metal-rich) cold gas would forcibly remain after an episode of star formation. We have deliberately changed this. Now, mass is returned and metals added to the ISM before we apply feedback to the ISM. Provided there is sufficient energy, this means that returned gas can immediately find its way to the CGM.

We explored the possibility of a distributive feedback model in Dark Sage, where the energy from supernovae in one annulus could affect the gas in adjacent annuli. Our default assumption was that energy would be distributed from an annulus centre isotropically. For infinitesimally wide annuli, one could imagine that most of the energy from feedback should be directed towards neighbouring annuli instead of out of the disc. However, with our set-up in Dark Sage, it is almost always true that the width of a disc annulus is much greater than the disc scale height, implying a negligible fraction of feedback energy should transfer through the disc. As such, in the end, we reasonably approximated the feedback taking place in an annulus to be independent from other annuli in a galaxy.

9.1. Radio mode

In earlier versions of Dark Sage, the model for radio feedback followed Croton et al. (Reference Croton2016). In this mode, hot gas is passively accreted directly onto the BH, based on the BH’s Bondi (Reference Bondi1952) accretion rate. Combined with the ‘maximal cooling flow’ model of Nulsen & Fabian (Reference Nulsen and Fabian2000), the accretion rate can be derived as

(99) \begin{equation}\dot{m}_\mathrm{BH}^\mathrm{radio} = \frac{15 \pi \, G \, \bar{\mu} m_\mathrm{p}\, k_\mathrm{B} \, T_\mathrm{vir}}{16\, \Lambda(T_\mathrm{vir}, Z_\mathrm{hot})} m_\mathrm{BH}\end{equation}

(cf. Croton et al. Reference Croton2006), where $m_\mathrm{BH}$ is the BH’s mass. In Croton et al. (Reference Croton2016), an additional free parameter was added to equation (99). However, with the changes we have made to the model, we no longer need this. We note that we explicitly do not allow $\dot{m}_\mathrm{BH}^\mathrm{radio}$ to exceed the BH’s Eddington (Reference Eddington1926) accretion rate in the code, but it is typically orders of magnitude below this anyway.

A half-century-old standard, we assume that mass accreting onto a BH quasi-statically moves through equatorial, circular orbits — all the while losing energy and angular momentum — until reaching the orbit of lowest energy, where it is captured (see Lynden-Bell Reference Lynden-Bell1969; Bardeen Reference Bardeen1970; Bardeen, Press & Teukolsky Reference Bardeen1972; Novikov & Thorne Reference Novikov, Thorne, DeWitt and DeWitt1973; Shakura & Sunyaev Reference Shakura and Sunyaev1973; Stevens Reference Stevens2015). The power emanating from the accretion disc — which drives an AGN — is then given directly by the accretion rate:

(100) \begin{equation}\dot{E}_\mathrm{AGN} = \epsilon_\mathrm{AGN}\, \dot{m}_\mathrm{BH}\, c^2,\end{equation}

where c is the speed of light in a vacuum and $\epsilon_\mathrm{AGN}$ is the ratio of inertial (gravitational) mass lost by that accreted to its rest mass (at infinity), often referred to as the ‘radiative efficiency’. Note that $\dot{m}_\mathrm{BH}$ represents the rest-mass accretion rate of the BH, but $m_\mathrm{BH}$ represents the inertial mass of the BH. Over a sub-time-step, the growth of a BH is therefore

(101) \begin{equation}\Delta m_\mathrm{BH} = (1 - \epsilon_\mathrm{AGN})\, \dot{m}_\mathrm{BH}\, \Delta t .\end{equation}

The difference between rest mass and inertial mass was not properly accounted for in previous versions of the model. This means BH feedback is now even more self-regulatory. That is, for a fixed rest mass of accreted material, a BH can either grow by a lot and produce a weak AGN (low $\epsilon_\mathrm{AGN}$ ), or grow less and produce a stronger AGN (high $\epsilon_\mathrm{AGN}$ ). But the less the BH grows, the weaker any future radio-mode feedback will be. In principle, to conserve inertial mass (energy), that lost by the material accreted onto the BH should be added somewhere. Without doing so, the baryonic mass of the halo within which that BH resides will drop. In practice though, this is inconsequential not only because this comprises a negligible fraction of a halo’s baryonic mass, but also because it is automatically accounted for by Dark Sage’s infall prescription anyway. At the next time step, baryonic mass will be added to the CGM of that halo to top it back up to the cosmic baryon fraction (modulo any suppression from photoionisation feedback, and including ejected gas in the accounting — see Section 4). We implicitly approximate rest mass and inertial mass as equivalent for all other components of galaxies and haloes (as literally all cosmological simulations and SAMs do), which we do not expect to have any adverse effects whatsoever.

Under the above picture, the value of $\epsilon_\mathrm{AGN}$ depends exclusively on the spin of the BH (neglecting charge) and the relative orientation of the accretion disc to the spin vector (Bardeen et al. Reference Bardeen1972). While $\epsilon_\mathrm{AGN} = 0.0572$ for a non-rotating Schwarzschild (Reference Schwarzschild1916) BH, $\epsilon_\mathrm{AGN}$ can theoretically range from 0.0378 to 0.3210, given the anti-parallel and parallel cases for a Kerr (Reference Kerr1963) BH rotating at the Thorne (Reference Thorne1974) limit, and can reach 0.4226 if the Thorne limit is breached. However, we do not explicitly consider BH spin in Dark Sage (for an example of a SAM that does, see Fanidakis et al. Reference Fanidakis and Baugh2011). The addition of such a feature would be an entire project in and of itself. In the meantime, we opt to treat $\epsilon_\mathrm{AGN}$ as a free parameter that is applied uniformly for all BHs.

Equation (100) applies to both modes of accretion, but the way that that energy is handled depends on the mode. For the radio mode, the energy is used as a preventative form of feedback. This has been standard for SAMs since Croton et al. (Reference Croton2006) and Bower et al. (Reference Bower, Benson, Malbon, Helly, Frenk, Baugh, Cole and Lacey2006). After the gross cooling rate of a galaxy is calculated (Section 4), we calculate how much of that cooling would be prevented by the energy released from the BH accretion disc:

(102a) \begin{equation}\Delta m_\mathrm{heat}^\mathrm{radio} = \frac{\dot{E}_\mathrm{AGN}^\mathrm{radio}\, \Delta t}{e_\mathrm{hot} - \langle e_\mathrm{cold} \rangle} = \frac{\Delta E_\mathrm{AGN}^\mathrm{radio}}{e_\mathrm{hot} - \langle e_\mathrm{cold} \rangle}, \end{equation}

(102b) \begin{equation}\langle e_\mathrm{cold} \rangle \equiv \sum_{i=1}^{N_\mathrm{ann}} e_\mathrm{cold}^{(i)} \int_{j_\mathrm{inner}^{(i)}}^{j_\mathrm{outer}^{(i)}} \mathrm{PDF}_\mathrm{cool}(j)\, \mathrm{d} j,\end{equation}

providing the final pieces to equations (13–15).

The term ‘radio mode’ is derived from the presence of radio jets that are generated from significantly sub-Eddington rates of BH accretion. These pervade the CGM and are what keeps the gas hot, or, rather, prevents it from cooling. While we have no explicit consideration of AGN jets in Dark Sage (but see Raouf et al. Reference Raouf, Shabala, Croton, Khosroshahi and Bernyk2017, where there is an explicit consideration of this in the Radio-SAGE model), the radio mode still implicitly relies on there being jets. The influence of these jets persists over extended timescales. For this reason, AGN memory in a SAM is important; without it, a model can be vulnerable to spurious, non-physical increases and decreases in radio-mode heating.

In Croton et al. (Reference Croton2016), radio-mode memory was achieved by introducing a ‘heating radius’, which never decreased with time and effectively prevented gas in the CGM internal to it from cooling. This radius would increase when periods of radio AGN activity exceeded what the galaxy had had previously. This model was adopted in previous versions of Dark Sage, and has been adopted in other SAMs too (e.g. Lagos et al. Reference Lagos2018).

For the new Dark Sage, we achieve radio-mode memory differently. We abandon any consideration of a heating radius and instead move mass from the hot reservoir to the fountain reservoir equal to that offset from cooling:

(103) \begin{equation}\Delta m_\mathrm{hot \rightarrow fount}^\mathrm{radio} = \Delta m_\mathrm{heat}^\mathrm{radio}.\end{equation}

This does not affect energy conservation, as the fountain and hot reservoirs are defined to have the same specific energy. But this does create a lasting impact by effectively disabling the AGN-affected gas from cooling for a dynamical time. This is not only less contrived than the Croton et al. (Reference Croton2016) heating radius model, but it also allows the suppression of cooling from radio-mode feedback in a galaxy to wax and wane moderately over time, rather than effectively forcing it to stay the same or increase.

Provided $m_\mathrm{BH} > 0$ and $m_\mathrm{hot} > 0$ , radio-mode accretion and feedback will occur in a Dark Sage galaxy at every sub-time-step.

9.2. Quasar mode

Whereas the radio mode is coupled to the CGM and is primarily responsible for the quenching of high-mass galaxies, the quasar mode is coupled to the ISM and is the primary growth channel for BHs. In earlier versions of SAGE and Dark Sage, quasar-mode accretion and feedback were triggered in two instances: from disc instabilities and when galaxy mergers occurred. Instability-driven accretion in Dark Sage requires gas to first flow into the innermost annulus; any unstable gas there that is not consumed/removed by a starburst and its associated feedback (Section 7.2) is accreted onto the central BH. The previous treatment of merger-driven BH accretion in Dark Sage was ad hoc and not well motivated theoretically; its phenomenological design (see Kauffmann & Haehnelt Reference Kauffmann and Haehnelt2000; Croton et al. Reference Croton2016) was better suited for SAMs that only treat global galaxy properties. We have now abandoned this channel of quasar-mode accretion entirely, leaving only the instability channel. Though, this is not to say that mergers no longer affect when BHs grow. The additional gas placed in the relevant annuli after a merger inevitably drives an instability that typically cascades all the way into the galaxy centre. These merger-induced instabilities play a major role in BH growth. Effectively, all we have done is retire a redundant prescription; this has no adverse consequences on the model’s predictions and it means the model has fewer free parameters.

$\dot{m}_\mathrm{BH}^\mathrm{quasar}$ is set by the amount of unstable gas in the innermost annulus of a galaxy and one free parameter that controls what fraction of this gas is consumed in a starburst ( $f_\mathrm{move}^\mathrm{gas}$ , see Section 7.2). Unlike the radio mode, there is nothing to explicitly prevent quasar-mode accretion from exceeding the Eddington (Reference Eddington1926) rate. Equation (100) then gives the power that drives quasar-mode feedback.

Quasar-mode feedback is ejective, not preventative. The energy is first dumped into the innermost annulus of the ISM, working outward until all the energy output over a given sub-time-step is used. Quasar feedback on the ISM differs from stellar feedback in that all feedback-affected gas is sent to the outflowing reservoir, with none going to the fountain reservoir. The physical justification for these choices is that the gas nearest the BH should be maximally affected before energy from it starts coupling to the outer ISM. Using the pre-computed specific energy of the outflowing reservoir (Section 8.2), we effect

(104a) \begin{equation}\Delta m_\mathrm{cold \rightarrow outfl}^{(i) \rm quasar} = \mathrm{min} \left[m_\mathrm{cold}^{(i)}, \,\frac{\Delta E_\mathrm{left}^{(i)}}{e_\mathrm{outfl} - e_\mathrm{ cold}^{(i)}} \right],\end{equation}

(104b) \begin{equation}\Delta E_\mathrm{left}^{(i)} =\left\{\begin{array}{l@{\quad}r}\Delta E_\mathrm{AGN}^\mathrm{quasar}, & i=1\\[4pt]\Delta E_\mathrm{AGN}^\mathrm{quasar} - \sum_{k=1}^{i-1} \Delta E_\mathrm{used}^{(k)}, & i \geq 2 \\\end{array}\right.,\end{equation}

(104c) \begin{equation}\Delta E_\mathrm{used}^{(k)} = \Delta m_\mathrm{cold \rightarrow outfl}^{(k) \rm quasar} \left(e_\mathrm{outfl} - e_\mathrm{cold}^{(k)}\right) .\end{equation}

This equation is solved for each annulus i until either reaching $i = N_\mathrm{ann}$ or $\Delta E_\mathrm{left}^{(i)} = 0$ . Whenever the former is the trigger to stop, there must still be some non-zero excess energy

(105) \begin{equation}\Delta E_\mathrm{excess}^\mathrm{quasar} = \Delta E_\mathrm{AGN}^\mathrm{quasar} - \sum_{i=1}^{N_\mathrm{ann}} \Delta E_\mathrm{used}^{(i)}.\end{equation}

Identical to our treatment of excess energy from stellar feedback, we use as much of this energy as possible to eject gas from the hot reservoir:

(106a) \begin{equation}\Delta m_\mathrm{hot \rightarrow outfl}^\mathrm{quasar} = \mathrm{min} \left[m_\mathrm{hot}, \,\frac{\Delta E_\mathrm{excess}^\mathrm{quasar}}{e_\mathrm{outfl} - e_\mathrm{hot}} \frac{m_\mathrm{hot}}{m_\mathrm{hot} + m_\mathrm{fount}} \right],\end{equation}

(106b) \begin{equation}\Delta m_\mathrm{fount \rightarrow outfl}^\mathrm{quasar} = \mathrm{min} \left[m_\mathrm{fount}, \,\frac{\Delta E_\mathrm{excess}^\mathrm{quasar}}{e_\mathrm{outfl} - e_\mathrm{ hot}} \frac{m_\mathrm{fount}}{m_\mathrm{hot} + m_\mathrm{fount}} \right].\end{equation}

The main differences from the Stevens et al. (Reference Stevens, Croton and Mutch2016) implementation of quasar-mode AGN feedback is our far more comprehensive consideration of the specific energy differences between reservoirs, the allowance for an annulus to have its gas partially ejected (as opposed to it being ‘all or nothing’), and the timescale on which this gas reincorporates back into the hot reservoir (described in the next section).

10.1. Fountaining in the CGM

Gas in the fountain reservoir is intended to reincorporate into the hot-gas reservoir to become available for cooling again. We must decide the timescale that that gas remains in the transitory fountain reservoir. In the absence of a detailed consideration of how efficiently gas mixes (beyond the scope of this version of Dark Sage), we make a simple assumption that gas takes a halo dynamical time to transit the fountain reservoir. We track an average ‘fountain time’ whenever mass is added to the reservoir:

(107) \begin{align}t_\mathrm{fount}(t + \Delta t) = \mathrm{max} \bigg[0, \frac{m_\mathrm{fount}(t)\, t_\mathrm{fount}(t) + \Delta m_\mathrm{\Sigma \rightarrow fount}\, t_\mathrm{dyn}}{m_\mathrm{fount}(t) + \Delta m_\mathrm{\Sigma \rightarrow fount}}- \Delta t \bigg],\end{align}

where $\Delta m_\mathrm{\Sigma \rightarrow fount}$ is the sum of contributions to the fountain reservoir for the sub-time-step.

For each sub-time-step, we move a fraction of the fountaining mass to the hot reservoir in proportion to the fraction of $t_\mathrm{fount}$ that has elapsed:

(108) \begin{equation}\Delta m_\mathrm{fount \rightarrow hot} = \mathrm{min} \left[1,\,\frac{\Delta t}{t_\mathrm{fount}}\right]\, m_\mathrm{fount}.\end{equation}

10.2. Ejection of outflowing gas

Gas is in the outflowing reservoir is, by definition, eventually going to escape the halo. Again, we need a timescale to decide how long gas added to the outflowing reservoir will remain in it before being transferred to the ejected reservoir. We estimate the instantaneous outflow timescale based on the specific energy of the gas added to the reservoir at that time. If we assume the thermal energy of the outflowing gas remains constant, we can approximate the average velocity of outflowing gas travelling as that which balances energy at the halo centre and at the virial radius:

(109a) \begin{equation}t_\mathrm{outfl}^\mathrm{inst} = \frac{2\, R_\mathrm{200c}}{v^\mathrm{inst}_\mathrm{outfl}(0) + v_\mathrm{outfl}^\mathrm{inst}(R_\mathrm{200c})},\end{equation}

(109b) \begin{equation}v_\mathrm{outfl}^\mathrm{inst}(R) \equiv \sqrt{2\, \left(e_\mathrm{outfl}^\mathrm{inst} - \Phi(R) - e_\mathrm{hot}^\mathrm{therm}\right)}.\end{equation}

When mass is added to the outflow reservoir from supernova feedback,

\begin{align*}e_\mathrm{outfl}^\mathrm{inst} = e_\mathrm{outfl}^\mathrm{SN}.\end{align*}

Otherwise

\begin{align*}e_\mathrm{outfl}^\mathrm{inst} = e_\mathrm{outfl}.\end{align*}

In the same manner as the fountain reservoir, we track an average outflow time as mass is continually added to the outflowing reservoir, ensuring to reduce this as cosmic time elapses:

(110) \begin{align}t_\mathrm{outfl}(t + \Delta t) = \mathrm{max} \bigg[0, \frac{m_\mathrm{outfl}(t)\, t_\mathrm{outfl}(t) + \Delta m_\mathrm{\Sigma \rightarrow outfl}\, t_\mathrm{outfl}^\mathrm{inst}}{m_\mathrm{outfl}(t) + \Delta m_\mathrm{\Sigma \rightarrow outfl}} - \Delta t \bigg],\end{align}

where $\Delta m_\mathrm{\Sigma \rightarrow outfl}$ represents the additions to the outflowing reservoir from all sources for the sub-time-step. For central galaxies and satellites outside their host’s virial radius, mass is added to their ejected reservoir according to:

(111) \begin{equation}\Delta m_\mathrm{outfl \rightarrow ejec} = \mathrm{min} \left[1,\,\frac{\Delta t}{t_\mathrm{outfl}}\right]\, m_\mathrm{outfl}.\end{equation}

Satellites inside the virial radius of their host do not have an ejected reservoir by construction. For these galaxies, the mass on the left-hand side of equation (111) is instead placed in the fountain reservoir of the central (and its fountain timescale is updated accordingly).

10.3. Effect of halo growth

As a halo grows in mass, the depth of its potential well and thermal energy both increase. This can cause outflowing gas that previously had enough energy to escape the halo to no longer have enough. Whenever a (sub)halo increases in mass, we first check whether any of the outflowing gas can still reach the virial radius. This requires

\begin{align*}\Phi(R_\mathrm{200c}, t) + e_\mathrm{hot}^\mathrm{therm}(t) > e_\mathrm{outfl}(t - \Delta t).\end{align*}

If this inequality is not satisfied, all gas in the outflowing reservoir is transferred to the fountain reservoir. Otherwise, we compute the average specific energy required for the outflowing gas to maintain its outflow timescale.

While we only explicitly track the average outflow timescale and energy of gas in the CGM, one would physically expect the outflowing gas to have a distribution of energy. The lower-energy gas in this distribution must be that which is moved into the fountain reservoir as the result of halo growth. This should, therefore, raise the specific energy of the material that remains in the outflowing reservoir. In principle, this should reduce the outflow time. However, as the halo grows, the assumed temperature of the remaining outflowing gas rises — which means it gets slower, if we assume its total energy is conserved — and the distance it must travel to reach $R_\mathrm{200c}$ also increases. Both of these effects should raise the outflow time. Because there are competing effects that should both raise and lower $t_\mathrm{outfl}$ , we take the agnostic stance of keeping the outflow time constant, and then use energy conservation to calculate the fraction of the outflowing gas that should be moved to the fountain reservoir:

(112a) \begin{equation}\Delta m_\mathrm{outfl \rightarrow fount}\, e_\mathrm{fount}(t) + m_\mathrm{outfl}(t)\, e_\mathrm{outfl}(t)\\ = m_\mathrm{outfl}(t - \Delta t)\, e_\mathrm{outfl}(t - \Delta t),\end{equation}

(112b) \begin{equation}e_\mathrm{outfl}(t) = \Phi(R_\mathrm{200c},t) + e_\mathrm{hot}^\mathrm{therm}(t) + \frac{1}{2} v_\mathrm{outfl}^2(R_\mathrm{200c},t),\end{equation}

(112c) \begin{align}v_\mathrm{outfl}(R_\mathrm{200c},t) = \frac{2\, R_\mathrm{200c}(t)}{t_\mathrm{outfl}} - \sqrt{2 \left[ e_\mathrm{outfl}(t {-} \Delta t) {-} \Phi(0,t) - e_\mathrm{hot}^\mathrm{therm}(t) \right]}.\end{align}

10.4. Reincorporation of ejected gas

When gas is added to the ejected reservoir from the outflowing reservoir, we calculate an expected time for that gas to reincorporate back into the halo. From the specific energy of the outflowing reservoir, we can calculate both the speed of the ejected gas at the virial radius, and the radius at which the gas will cease to have kinetic energy and must turn around. We approximate the average speed the gas moves during this period as half its ejected speed at the virial radius. With the understanding that the gas must travel the distance from $R_\mathrm{200c}$ to $R_\mathrm{turn}$ twice, we effect

(113a) \begin{equation}t_\mathrm{reinc}^\mathrm{inst} = \frac{4 \left( R_\mathrm{turn} - R_\mathrm{200c} \right)}{v_\mathrm{outfl}(R_\mathrm{200c})},\end{equation}

(113b) \begin{equation}\Phi(R_\mathrm{turn}) + e_\mathrm{hot}^\mathrm{therm} = e_\mathrm{outfl}.\end{equation}

To solve these simultaneous equations efficiently in the code, we approximate the shape of the potential outside the virial radius with an analytic NFW profile.

As with the other timescales, we average the reincorporation timescale as more mass is added:

(114) \begin{align}t_\mathrm{reinc}(t + \Delta t) = \mathrm{max} \bigg[0, \frac{m_\mathrm{ejec}(t)\, t_\mathrm{reinc}(t) + \Delta m_\mathrm{outfl \rightarrow ejec}\, t_\mathrm{reinc}^\mathrm{inst}}{m_\mathrm{ejec}(t) + \Delta m_\mathrm{outfl \rightarrow ejec}} - \Delta t \bigg].\end{align}

Ejected mass is then reincorporated directly into the hot-gas reservoir at a rate set by

(115) \begin{equation}\Delta m_\mathrm{eject \rightarrow hot} = \mathrm{min} \left[1,\,\frac{\Delta t}{t_\mathrm{reinc}} \right] m_\mathrm{ejec}.\end{equation}

10.5. Baryon breakdown

With all the baryon reservoirs and the equations for how mass moves between them in Dark Sage now described, it is informative to see how the baryonic mass inside haloes is broken between them at $z = 0$ . We show precisely this in Figure 11.

Figure 11. The breakdown of baryons in Dark Sage haloes at $z = 0$ . We calculate the average mass of each reservoir for all haloes in each halo mass bin of width 0.05 dex. Bars are stacked on top of each other, such that the height of all bars in a column, less the ejected gas, gives the total baryon fraction inside the halo. All matter for satellite galaxies inside the virial radius is included. Black holes are reasonably neglected from this plot. The dotted horizontal line represents the cosmic baryon fraction. Photoionization heating suppresses the total baryon fraction below the cosmic value for $M_\mathrm{200c} \lesssim 10^{13}\,\mathrm{M}_\odot$ . The ejected gas contributes towards the halo baryon fraction for the purposes of cosmic accretion, even though it is outside the halo (as discussed in Section 4). The outflowing-gas reservoir is nigh negligible at all halo masses.

The stellar mass fraction as a function of halo mass is as one would expect, peaking at $M_\mathrm{200c} \simeq 10^{12}\,\mathrm{M}_{\odot}$ . The CGM fraction increases monotonically, dominating the baryon budget at high halo masses. Of the mass in the CGM, the outflowing reservoir contains little mass at all halo mass scales, showing that the time spent in this transitory reservoir is less than the time spent in the ejected reservoir before reincorporation. A significant portion of the CGM is in the fountain reservoir for $M_\mathrm{200c} \lesssim 10^{13}\,\mathrm{M}_{\odot}$ . This highlights that most of the gas reheated from stellar feedback does not have sufficient energy to escape the halo, and that notably less star formation and stellar feedback operates in the group–cluster regime. At these mass scales, after fountaining gas reincorporates into the hot-gas reservoir, it is offset from cooling from radio-mode AGN feedback.

The age bins in the stellar reservoirs allow us also reconstruct the star formation histories of galaxies on a per-component basis. As an example of Dark Sage’s abilities in this area, we use the summed birth mass of stars in common stellar components across galaxies within wide stellar-mass bins to build cosmic star formation density histories (CSFHs) in the top panels of Figure 12. Because each stellar component has its own metallicity, we can also reconstruct the average metallicity of the star-forming gas as a function of cosmic time for each component. The averages of these for the binned galaxies are shown in the bottom panels of Figure 12. This plot shows that high-mass galaxies are dominated in stellar mass by merger-driven bulges, while the stellar mass in low-mass galaxies is predominantly in discs. It also shows that spheroidal stars (be they in either bulge or the intrahalo stars; and even the LIGS) have a notably higher average age than disc stars. Except for the highest-mass galaxies, merger-driven bulges tend to be metal-poor relative to discs. There is also a smooth transition in both the age and metallicity of disc stars, getting younger and more metal-rich from the outskirts to the centre and into the instability-driven bulge. These predictions from Dark Sage can be tested with component-based SED-fitting models applied to large-scale galaxy surveys (e.g. Robotham et al. Reference Robotham, Bellstedt and Driver2022; also see Bellstedt et al. Reference Bellstedt2020).

Figure 12. Top panels: Contributions to the cosmic star formation density history based on the stellar components of galaxies in three wide stellar-mass bins. This is reconstructed from the $z = 0$ output of Dark Sage using the stellar-age bins in each stellar component after multiplying the mass in those components by $m_\mathrm{*,pop}^\mathrm{form} / m_\mathrm{*,pop}^\mathrm{rem}$ . Bottom panels: the average metallicity of those stars, equal to the average metallicity of the star-forming gas at the time of formation. The dashed, thin lines are the per-annulus contribution to the disc stars, smoothly changing colour from teal for outer annuli ( $i \rightarrow N_\mathrm{ann}$ ) to blue for intermediate annuli, to magenta for inner annuli ( $i \rightarrow 1$ , i.e. towards the instability-driven bulge, which functions like a ‘zeroth’ annulus).

11.1. Central galaxies

While we have no explicit treatment of large-scale environmental effects on central haloes, this is considered implicitly from how the galaxies are evolved within the merger trees of the N-body simulation. Where an increase in virial mass leads to the addition of new baryons in the halo in order to maintain a target baryon fraction (Section 4), so too are baryons removed when the virial mass decreases to ideally maintain this same fraction. Only the CGM and IHS of a central are stripped in this circumstance. This means that the baryon fraction can become higher than the target value in rare circumstances where removing the CGM and IHS outright is not enough to bring it back down.

When the CGM and IHS are stripped this way, they are done so proportionally to their mass. Those stripped gas and stars are transferred to the central’s LIGM and LIGS, respectively. This means that those baryons will return to the halo preferentially over fresh primordial gas, should the halo later grow (Section 4). Helpfully, this mitigates against potential halo-finding and/or tree-building issues that can lead to artificial, momentary losses in halo mass (Behroozi et al. Reference Behroozi, Wechsler, Wu, Busha, Klypin and Primack2013; Srisawat et al. Reference Srisawat2013; Elahi et al. Reference Elahi2019a,b); should this happen in Dark Sage, baryons will initially be taken out of the halo, but then the same baryons (i.e. gas with the same metallicity and stars with the same metallicity and age distribution) will return a snapshot later.

Note that in a change from previous versions of the model, we do not remove the ejected gas of central haloes when they are tidally stripped, as this mass is already considered to be beyond the virial radius.

11.2. Satellite galaxies

Much more detail is given to (and needed for) the environmental effects on satellites. To start with, immediately when a galaxy becomes a satellite, all gas in its LIGM is transferred to the host’s equivalent reservoir. This is simply in line with how we have defined that reservoir; it does not make sense to track a LIGM for a satellite. The same applies for LIGS and LIGBHs.

Satellites are shut off from cosmological gas accretion entirely. This starves the CGM from any replenishment. However, gas can continue to cool from the CGM of the satellite onto its ISM. Likewise, while the satellite remains outside its host virial radius, its CGM will continue to reincorporate gas from the satellite’s ejected reservoir. This is physically justifiable from the results from hydrodynamic simulations that have shown that satellites can still continue to accrete new gas, even while they are being affected by environmental effects that remove gas (Wright et al. Reference Wright and Lagos2022).

The CGM of a satellite is gradually stripped, an effect often referred to as ‘strangulation’ in the literature. In a change from earlier versions of Dark Sage, we now also strip the intrahalo stars of a satellite gradually, alongside the CGM (previously, the IHS was instantly transferred to the central upon infall). While this change is inconsequential for galaxy evolution, it allows us to track the stellar haloes of groups after they fall into a cluster. The rate of stripping of the CGM and IHS is based on (as follows) the subhalo mass-loss rate from the N-body simulation, which is caused by gravitational tides. We routinely check the total baryonic mass of the satellite subhalo relative to its expected baryon fraction and total subhalo mass. Should the former be above expectation, we proportionally remove mass (and metals) from the CGM and IHS to bring the baryonic mass as close to expectation as possible,Footnote o transferring that mass to the respective reservoirs of the central. The decision to strip the CGM and IHS of satellites this way is to be consistent with how those same reservoirs are stripped for centrals when their halo mass drops (Section 11.1). This philosophy extends from Croton et al. (Reference Croton2016).

In truth, the above paragraph is only an accurate description of a satellite that has crossed the virial radius of its host. Until that point, we preference removal of the ejected reservoir of the satellite over the CGM. Any stripped gas and stars of satellites outside $R_\mathrm{200c}^\mathrm{host}$ actually go to the respective local intergalactic reservoir. Only if the ejected gas mass of the satellite is zero does Dark Sage start to strip the CGM. This is expedited by the choice that once the satellite crosses $R_\mathrm{200c}^\mathrm{host}$ , its ejected reservoir is transferred to the central’s CGM immediately and in its entirety. The assumption is that the gravitational influence of the host halo is sufficiently larger than the satellite, such that any ejected gas can no longer be reincorporated into the satellite’s CGM. Because this gas is lost inside the virial radius of the host, there should arguably be no delay in its incorporation to the host’s CGM, hence why we do not put it into the central’s ejected reservoir.

Material from satellite galaxies themselves (cold gas and stars in the disc) is also stripped, but follows a different procedure. For this, we calculate the strength of both tides and ram pressure on the satellite galaxy. Previous versions of Dark Sage only accounted for ram pressure, which only affects the gas. However, we only turn these processes on once $m_* + m_\mathrm{cold} > m_\mathrm{hot}$ . The reason for this is that the hot-gas halo should initially shield the cold-gas disc from ram pressure, as should the CGM and IHS experience tidal stripping before the baryons in the galaxy disc that are more tightly gravitationally bound. Our exact criterion for switching this shielding off is not theoretically derived; it is simply maintained from Stevens & Brown (Reference Stevens and Brown2017). Other criteria for the same phenomenon have been used in other SAMs (e.g. Cora et al. Reference Cora2018).

We approximate the tidal radius of a satellite galaxy based on the formula of Tormen, Diaferio & Syer (Reference Tormen, Diaferio and Syer1998, cf. their equation 6):

(116a) \begin{equation}r_\mathrm{tidal} = R_\mathrm{sat}\, M_\mathrm{sat}^{1/3}\, \Upsilon^{-1/3},\end{equation}

(116b) \begin{equation}\Upsilon \equiv\left\{\begin{array}{l@{\quad}r}2\, M_\mathrm{halo}({<}R_\mathrm{sat}) & \\[4pt] - R_\mathrm{sat} \left. \dfrac{\mathrm{d} M_\mathrm{halo}({<}R)}{\mathrm{d} R}\right|_{R=R_\mathrm{sat}}, & R_\mathrm{sat} < R_\mathrm{200c}^\mathrm{host}\\[4pt] 2\, M_\mathrm{200c}^\mathrm{host}, & R_\mathrm{sat} \geq R_\mathrm{200c}^\mathrm{host}\\\end{array}\right..\end{equation}

Strictly, this is only accurate for satellites with circular orbits where $r_\mathrm{tidal} \ll R_\mathrm{sat}$ and $M_\mathrm{sat} \ll M_\mathrm{halo}({<}R_\mathrm{ sat})$ (and therefore $M_\mathrm{sat} \ll M_\mathrm{200c}^\mathrm{host}$ ).Footnote p However, we only enforce the lattermost criterion for tidal stripping to take place, requiring that $M_\mathrm{sat} < 0.1\,M_\mathrm{halo}({<}R_\mathrm{sat})$ (and by extension that $M_\mathrm{sat} < 0.1\,M_\mathrm{200c}^\mathrm{host}$ ). With that criterion passed, we remove all gas and stars in each disc annulus where $r_\mathrm{outer}^{(i)} \geq r_\mathrm{tidal}$ . The stripped gas and stars are respectively transported to the CGM and IHS of the central if $R_\mathrm{sat} \leq R_\mathrm{200c}^\mathrm{host}$ . If the satellite is instead outside the virial radius, the stripped gas and stars go to the respective local intergalactic reservoirs. In the instance that $r_\mathrm{tidal} < r_\mathrm{outer}^{(1)}$ , we trigger a tidal disruption event where all the satellite’s baryons are transferred to the CGM and IHS, and the satellite’s existence is no longer considered thereafter. As we explain in Section 12, this is not the only means by which a disruption event is triggered. With the exception of tidal disruption events, the stellar bulge of a satellite galaxy is left unaffected by environmental effects in Dark Sage.

With tidal effects resolved, any remaining gas in the ISM of a satellite is then subject to ram pressure stripping. Consistent with previous versions of Dark Sage, for each annulus, we calculate the ram pressure exerted on a satellite galaxy as

(117) \begin{equation}P_\mathrm{ram} = \rho_\mathrm{hot}^\mathrm{cen}(R_\mathrm{sat})\, V_\mathrm{sat}^2\end{equation}

and compare this against the gravitational restoring force per unit area for each annulus, approximated as

(118) \begin{equation}P_\mathrm{restore}^{(i)} =\left\{\begin{array}{l@{\quad}r}2\pi\, G \, \Sigma_\mathrm{cold}^{(i)} \left( \Sigma_\mathrm{cold}^{(i)} + \Sigma_{*}^{(i)} \right), & \theta_\mathrm{disc} \leq 10^\circ \\[5pt]2\pi\, G \, \Sigma_\mathrm{cold}^{(i)\,2}, & \theta_\mathrm{disc} > 10^\circ \\\end{array}\right.\end{equation}

(cf. Gunn & Gott Reference Gunn and Gott1972). For each annulus where $P_\mathrm{ram} > P_\mathrm{restore}^{(i)}$ , its gas is transferred to the central’s CGM (or LIGM if it is beyond $R_\mathrm{200c}^\mathrm{host}$ ).

12.1. Merger or disruption event?

In the above situation, how do we then choose between merging the satellite or disrupting it? The answer lies in how long the subhalo has been a satellite. The moment a subhalo becomes a satellite, a predicted merger timescale is calculated. If less than this time has elapsed when the above is triggered, the satellite presumably must have been disrupted. Otherwise, it should have reasonably merged already or be in the process of doing so imminently. In Croton et al. (Reference Croton2016), this predicted timescale was based on the dynamical friction model of Binney & Tremaine (Reference Binney and Tremaine1987). However, Poulton et al. (Reference Poulton, Power, Robotham and Elahi2020) has since shown that this timescale — in addition to several other merger timescale formulae commonly used in the literature — only works for a limited range of infall scenarios, and even then carries a factor of $\sim$ 2 uncertainty. Building on this, Poulton et al. (Reference Poulton, Power, Robotham, Elahi and Lagos2021) built a more sophisticated model by fitting a function to timescales derived from the meticulous tracking of subhalo orbits in cosmological simulations with the OrbWeaver code. We now adopt this timescale in Dark Sage.

The Poulton et al. (Reference Poulton, Power, Robotham, Elahi and Lagos2021) timescale, hereafter $t_\mathrm{merge}$ , can be calculated as follows (cf. their equations 4–7, 11–12).

(119a) \begin{equation}\frac{t_\mathrm{merge}}{5.5\, R_\mathrm{peri}^{0.2}} =\left\{\begin{array}{l@{\quad}r}\sqrt{\dfrac{R_\mathrm{200c}^\mathrm{host}\, R_\mathrm{sat}^{1.6}}{G\, M_\mathrm{halo}({<}R_\mathrm{sat})}}, & R_\mathrm{sat} < R_\mathrm{200c}^\mathrm{host}\\\\\dfrac{R_\mathrm{200c}^\mathrm{host}\, R_\mathrm{sat}^{0.3}}{\sqrt{G\, M_\mathrm{200c}^\mathrm{host}}}, & R_\mathrm{sat} \geq R_\mathrm{200c}^\mathrm{host}\\\end{array}\right.,\end{equation}

(119b) \begin{equation}R_\mathrm{peri} = \frac{J_\mathrm{orb}^2}{(1+\varepsilon)\, G\, M_\mathrm{host}\, M_\mathrm{sat}\, M_\mathrm{redu}}, \end{equation}

(119c) \begin{equation}\varepsilon = \sqrt{ 1 + \frac{2\, E_\mathrm{orb}\, J_\mathrm{orb}^2}{(G\, M_\mathrm{host}\, M_\mathrm{sat})^2\, M_\mathrm{redu}} }, \end{equation}

(119d) \begin{equation}J_\mathrm{orb} = M_\mathrm{redu}\, | \vec{R}_\mathrm{sat} \times \vec{V}_\mathrm{sat} |, \end{equation}

(119e) \begin{equation}E_\mathrm{orb} = \frac{M_\mathrm{redu}\, V_\mathrm{sat}^2}{2} + M_\mathrm{sat}\, \Phi(R_\mathrm{sat}), \end{equation}

(119f) \begin{equation}M_\mathrm{redu} \equiv \frac{M_\mathrm{sat}\, M_\mathrm{host}}{M_\mathrm{sat} + M_\mathrm{host}},\end{equation}

(119g) \begin{equation}M_\mathrm{host} \equiv\left\{\begin{array}{l@{\quad}r}M_\mathrm{halo}({<}R_\mathrm{sat}), & R_\mathrm{sat} < R_\mathrm{200c}^\mathrm{host}\\M_\mathrm{200c}^\mathrm{host}, & R_\mathrm{sat} \geq R_\mathrm{200c}^\mathrm{host}\\\end{array}\right..\end{equation}

Every term in the above equations is known/calculable for each subhalo in Dark Sage, making this straightforward to implement. $R_\mathrm{sat}$ represents the distance between the central galaxy’s centre of mass and the satellite galaxy’s centre of mass. $R_\mathrm{peri}$ represents the expected separation of the galaxies at the pericentre of their orbit with total angular momentum $J_\mathrm{orb}$ , energy $E_\mathrm{orb}$ , and ellipticity $\varepsilon$ . $V_\mathrm{sat}$ is the relative speed between the two galaxies. $M_\mathrm{redu}$ is the ‘reduced’ mass of the two-galaxy system. The mass terms for the central and satellite — $M_\mathrm{host}$ and $M_\mathrm{sat}$ , respectively — account for all matter, including dark matter.

When a disruption event occurs, if the satellite galaxy is inside the host’s virial radius, the remaining entirety of the satellite’s CGM, ISM, and ejected gas are transferred to the central’s hot-gas reservoir. Similarly, the entirety of the satellite’s stellar mass (including its IHS) is transferred to the central’s IHS. The size of the central’s IHS is then updated based on where the satellite was disrupted:

(120) \begin{equation}\langle R \rangle_\mathrm{IHS}^\mathrm{cen} (t + \Delta t) = \frac{m_\mathrm{IHS}^\mathrm{cen}(t)\, \langle R \rangle_\mathrm{IHS}^\mathrm{cen}(t) + \left( m_*^\mathrm{sat} + m_\mathrm{IHS}^\mathrm{sat} \right) R_\mathrm{sat}} {m_\mathrm{IHS}^\mathrm{cen}(t + \Delta t)},\end{equation}

where the average size of the IHS relates to its scale radius, $a_\mathrm{IHS}$ , of a truncated Hernquist (Reference Hernquist1990) sphere via

(121) \begin{equation}\frac{\langle R \rangle_\mathrm{IHS}}{a_\mathrm{IHS}} = 2 \left(\frac{a_\mathrm{IHS}}{R_\mathrm{200c}} + 1\right)^2 \ln \left(\frac{R_\mathrm{200c}}{a_\mathrm{IHS}} + 1 \right) - 2 \frac{a_\mathrm{IHS}}{R_\mathrm{200c}} -3.\end{equation}

This equation can be solved for $a_\mathrm{IHS}$ iteratively. If the disrupted satellite was outside the virial radius, the gas and stars of that satellite are instead transferred to the LIGM and LIGS components of the halo, respectively.

The resolution of a merger is more complicated. How the remaining baryons are handled depends on whether the merger is ‘major’ or ‘minor’.

A minor merger represents when one galaxy vastly outweighs the other. Conversely, a merger is deemed major when the merging systems have comparable mass. In practice, a major merger is defined in Dark Sage when

(122) \begin{align}\mathrm{min}\left[m_\mathrm{*}^\mathrm{cen} + m_\mathrm{cold}^\mathrm{cen}, m_\mathrm{*}^\mathrm{sat} + m_\mathrm{cold}^\mathrm{sat} \right] \\ \geq f_\mathrm{major} \times \mathrm{max}\left[m_\mathrm{*}^\mathrm{cen} + m_\mathrm{cold}^\mathrm{cen}, m_\mathrm{*}^\mathrm{sat} + m_\mathrm{cold}^\mathrm{sat} \right]\end{align}

(while the central always has a greater total mass than its satellites, it is not strictly guaranteed that it has greater galactic baryonic mass). The treatment of minor mergers remains similar to Stevens et al. (Reference Stevens, Croton and Mutch2016), but the handling of the ISM in major mergers has notably changed. To more accurately reflect the implicit assumption in the handling of minor mergers that $M_\mathrm{host} \gg M_\mathrm{sat}$ , we set $f_\mathrm{major} = 0.1$ . Previous versions of the model had this set to 0.3.

12.2. Resolving minor mergers

The ISM of a satellite is assumed to be dispersed uniformly across a subset of the annuli of the central’s disc. The orbital specific angular momentum of the satellite is projected onto the central disc’s rotation axis to find the annulus with equivalent j. Rather than dumping the entire ISM of the satellite into this one annulus of the central, we deposit it across several annuli. This is intended to account for the fact there is actually a distribution of orbital j of the satellite’s ISM relative to the central, given that the satellite’s ISM has its own rotation prior to merging. We approximate this by selecting the maximum number of annuli ( $l-k+1 \geq 1$ ) where

(123a) \begin{equation}j_\mathrm{outer}^{(k) \mathrm{cen}} > \left( \vec{j}_\mathrm{orb} \cdot \hat{j}_\mathrm{cold}^\mathrm{cen} \right) - j_\mathrm{spread},\end{equation}

(123b) \begin{equation}j_\mathrm{inner}^{(l) \mathrm{cen}} < \left( \vec{j}_\mathrm{orb} \cdot \hat{j}_\mathrm{cold}^\mathrm{cen} \right) + j_\mathrm{spread},\end{equation}

(123c) \begin{equation}\vec{j}_\mathrm{orb} \equiv \vec{R}_\mathrm{sat} \times \vec{V}_\mathrm{sat},\end{equation}

(123d) \begin{equation}j_\mathrm{spread} \equiv \left( \vec{R}_\mathrm{sat} \cdot \hat{j}_\mathrm{cold}^\mathrm{cen} \right)\, \left( \hat{j}_\mathrm{cold}^\mathrm{sat} \cdot \hat{j}_\mathrm{cold}^\mathrm{cen} \right)\, V_\mathrm{200c,infall}^\mathrm{sat}.\end{equation}

For each annulus i of the central’s gas disc where $k \leq i \leq l$ , we deposit $m_\mathrm{cold}^\mathrm{sat} / (l-k+1)$ into that annulus. Metals are deposited in proportion to the satellite’s average ISM metallicity. This is identical to the procedure from Stevens et al. (Reference Stevens, Croton and Mutch2016).

In reality, the manner in which we calculate $j_\mathrm{spread}$ or how we distribute the satellite’s ISM across annuli k through l in the central are of little consequence. The sudden addition of mass to those annuli typically drives a local instability, which promptly redistributes some of that mass over a larger number of annuli (and also removes some of that gas from those annuli due to feedback from instability-driven star formation — see Section 7.2).

The above takes place regardless of whether the orbit of the satellite was prograde or retrograde relative to the central’s gas disc. However, if the orbit was retrograde, this would artificially add angular momentum to the descendent galaxy’s disc, when instead there should be a net loss. To rectify this violation of angular momentum conservation, after a retrograde minor merger, we effectively shrink the descendant’s gas disc by proportionally transferring matter inward, such that the resulting $j_\mathrm{cold}^\mathrm{cen}$ is as it ought to be, that is,

(124) \begin{equation}\vec{j}_\mathrm{cold}^\mathrm{cen}(t + \Delta t) = \Bigl[ m_\mathrm{cold}^\mathrm{cen}(t) + m_\mathrm{cold}^\mathrm{sat} \Bigl]^{-1} \Bigl[ m_\mathrm{cold}^\mathrm{cen}(t)\, \vec{j}_\mathrm{cold}^\mathrm{cen}(t) \\+ m_\mathrm{cold}^\mathrm{sat}\, \hat{j}_\mathrm{cold}^\mathrm{cen} \left( \vec{j}_\mathrm{orb} \cdot \hat{j}_\mathrm{cold}^\mathrm{cen} \right)\Bigl].\end{equation}

This equation is imposed before determining whether any local instabilities are triggered.

We maintain the choice that all stars in the satellite galaxy are transferred to the merger-driven bulge of the central from previous versions of the model. While we would expect collisional gas to settle into the disc structure of the ISM, it would be contrived to force the orbits of collisionless stars acquired from a minor merger to match the stellar disc. The merger-driven bulge exists in the model precisely to account for stars whose orbits are less ordered (there is no analogue for gas). While this principle has not changed, we have added several details in how the merger-driven bulge is updated after a minor merger. For example, we now also transfer any remaining IHS of the satellite to the central’s merger-driven bulge too. That is,

(125) \begin{equation}m_{\operatorname{m-bulge}}^\mathrm{cen}(t + \Delta t) = m_{\operatorname{m-bulge}}^\mathrm{cen}(t) + m_*^\mathrm{sat} + m_\mathrm{IHS}^\mathrm{sat}.\end{equation}

In all likelihood, most of the IHS of the satellite would have already been stripped (per Section 11) and any remaining IHS should be relatively nuclear (and therefore small compared to the size of the central’s bulge).

As another new feature in Dark Sage, we now approximately track the angular momentum and velocity dispersion of the merger-driven bulge, which in turn is used to calculate its size. After a minor merger, we update the angular momentum of the central by summing the rotational angular momentum of what the bulge originally had with the total rotational angular momentum of the satellite. That is, we conserve bulge rotation, but neglect any transfer of orbital angular momentum. The assumption is any orbital angular momentum is lost to dynamic friction before stars are added to the bulge. The specific angular momentum of the merger-driven bulge after a minor merger is then:

(126) \begin{align}\vec{j}_{\operatorname{m-bulge}}^\mathrm{cen}(t + \Delta t) & = m_{\operatorname{m-bulge}}^\mathrm{cen\,-1}(t + \Delta t) \left[ m_{\operatorname{m-bulge}}^\mathrm{cen}(t)\, \vec{j}_{\operatorname{m-bulge}}^\mathrm{cen}(t)\right.\nonumber\\[5pt]& \quad \left. +\, m_\mathrm{*,disc}^\mathrm{sat}\, \vec{j}_\mathrm{*,disc}^\mathrm{sat} + m_{\operatorname{m-bulge}}^\mathrm{sat}\, \vec{j}_{\operatorname{m-bulge}}^\mathrm{sat} \right].\end{align}

For the merger-driven bulge’s velocity dispersion, we add the mass-weighted squares of the dispersion of each relevant component and a further term that accounts for the radial motion of the satellite relative to the central. That is,

(127) \begin{align}\sigma_{\operatorname{m-bulge}}^\mathrm{cen\,2}(t + \Delta t) = m_{\operatorname{m-bulge}}^\mathrm{cen\,-1}(t + \Delta t) \Biggl[m_{\operatorname{m-bulge}}^\mathrm{cen}(t)\, \sigma_{\operatorname{m-bulge}}^\mathrm{cen\,2}(t) \\[5pt]+ m_{\operatorname{m-bulge}}^\mathrm{sat}\, \sigma_{\operatorname{m-bulge}}^\mathrm{sat\,2} + m_{\operatorname{i-bulge}}^\mathrm{sat}\, \sigma_{\operatorname{i-bulge}}^\mathrm{sat\,2}+ \left(m_*^\mathrm{sat} + m_\mathrm{IHS}^\mathrm{sat}\right) \\[5pt] \times \left(\vec{V}_\mathrm{sat} \cdot \hat{R}_\mathrm{sat} \right)^2+ \sum_{i=1}^{N_\mathrm{ann}} m_\mathrm{*,disc}^{(i) \mathrm{sat}}\, \sigma_\mathrm{*,disc}^{(i) \mathrm{sat}\,2} \Biggl] .\end{align}

The CGM and ejected reservoirs of the satellite are added to the hot-gas reservoir of the central. In essence, we assume there is no significant delay for that gas to become available for cooling. There is no impact on whether the satellite was inside or outside the central’s virial radius at the time the merger was triggered.

12.3. Resolving major mergers

During a major merger, the structure of the remnant galaxy is changed significantly from either of its progenitors. First, we define the plane for the new gas disc to exist, chosen to be that of the last orbital plane between the two galaxies. We then calculate the projected specific angular momentum component of each gas disc annulus in each galaxy relative to the centre of mass of the two-galaxy system in the direction of this plane. The annulus with the corresponding j in the new disc is where the gas from the old annulus is placed. Note the difference from earlier versions of the model in our choosing the centre-of-mass frame to determine how the new gas disc is built. Not only will the two progenitor discs be mixed, but the surface density profile of the new disc is likely to look entirely different to both progenitors. Similar to the case for minor mergers, if any retrograde gas is added to the new disc, we shrink it to ensure angular momentum conservation is satisfied in the plane. In earlier versions of the model, we triggered a specific starburst when a retrograde merger occurred. We no longer do this, as this is effectively taken care of by the instability prescription.

All stars in the discs and bulges of both progenitors are moved to the merger-driven bulge of the descendant after a major merger, while the two IHS components of the progenitors are summed to form the IHS of the descendant. We set the merger-driven bulge to have a specific angular momentum equal to the orbital specific angular momentum of the two-progenitor system in the centre-of-mass frame. The velocity dispersion of the merger-driven bulge is updated consistently with how it is during minor mergers (where terms for the central that are only present for the satellite in equation 127 are added). The size of the merger-driven bulge is calculated by solving the Jeans equation, in the same manner as for the instability-driven bulge, except this time accounting for the non-zero angular momentum:

(128) \begin{equation}a_{\operatorname{m-bulge}} = 3 \left[ \frac{1}{\sigma^2_{\operatorname{m-bulge}}} \left( \frac{\mathrm{d} \Phi}{\mathrm{d} R} - \frac{j_{\operatorname{m-bulge}}^2}{R^3} \right) - \frac{1}{R} \right]^{-1} - R.\end{equation}

The CGM and the ejected reservoir of the satellite are instead transferred to the fountain reservoir of the central, and the fountain time is updated when a major merger occurs. The idea here is that the collision of the two comparable mass CGMs should drive turbulence in the gas, causing a delay for it become available for cooling onto the remnant galaxy.