Dark Sage: Next-generation semi-analytic galaxy evolution with multidimensional structure and minimal free parameters

Abstract

After more than five years of development, we present a new version of Dark Sage, a semi-analytic model (SAM) of galaxy formation that breaks the mould for models of its kind. Included among the major changes is an overhauled treatment of stellar feedback that is derived from energy conservation, operates on local scales, affects gas gradually over time rather than instantaneously, and predicts a mass loading factor for every galaxy. Building on the model’s resolved angular momentum structure of galaxies, we now consider the heating of stellar discs, delivering predictions for disc structure both radially and vertically. We add a further dimension to stellar discs by tracking the distribution of stellar ages in each annulus. Each annulus–age bin has its own velocity dispersion and metallicity evolved in the model. This allows Dark Sage to make structural predictions for galaxies that previously only hydrodynamic simulations could. We present the model as run on the merger trees of the highest-resolution gravity-only simulation of the MillenniumTNG suite. Despite its additional complexity relative to other SAMs, Dark Sage only has three free parameters, the least of any SAM, which we calibrate exclusively against the cosmic star formation history and the $z = 0$ stellar and H i mass functions using a particle-swarm optimisation method. The Dark Sage codebase, written in C and python, is publicly available at https://github.com/arhstevens/DarkSage.

1. Context

Research into galaxy formation remains among the top priorities in the field of astronomy (National Academies of Sciences, Engineering, and Medicine 2021). With the likes of the James Webb Space Telescope in full operation and construction of the Square Kilometre Array underway, this is a promising decade for the field. The endgame of these multi-billion-dollar enterprises is to develop a comprehensive working theory for how the Universe works, of which galaxy formation is a major, crucial part. But the flood of incoming high-quality observations is only as powerful as the models we use to interpret those data. Galaxy formation is a highly complex problem to model, which requires us to ‘solve’ it with numerical simulations. While every model is flawed, simulations are the closest that extragalactic astrophysicists will ever come to conducting a controlled experiment. Having a variety of galaxy formation models, frameworks, and simulations — alongside observations — is paramount to test and develop our theoretical understanding of not just galaxies themselves, but cosmology too.

Two methods for simulating the formation of galaxies in a cosmological context have stood out in recent decades: semi-analytic models (SAMs) and hydrodynamic simulations. Hydrodynamic simulations simultaneously and self-consistently account for gravity, fluid dynamics, and all astrophysical processes deemed relevant to the formation of galaxies. SAMs separate the formation of the universe’s large-scale structure from that of galaxies by first constructing a universe with only gravity, and subsequently evolving galaxies in the gravitationally bound structures in that universe — referred to as ‘haloes’ — with relatively macroscopic descriptions of baryonic astrophysics. In truth, hydrodynamic simulations are still ‘semi-analytic’ at the level of a particle, which is generally comparable in mass to a giant molecular cloud or globular cluster of stars. Or, to frame this in reverse, all baryonic processes are ‘subgrid’ at the level of the (sub)halo in a SAM. A holistic summary of these methods of galaxy formation modelling can be found in a number of review articles (e.g. Somerville & Davé 2015; Vogelsberger et al. 2020).

In this paper, we present a new version of the SAM known as ‘Dark Sage’. Dark Sage has been in development since 2015, with two main versions previously released (Stevens et al. 2016, 2018). The model started as a modified version of SAGE (Croton et al. 2016). While the code architecture is still based on SAGE, the physical models implemented in Dark Sage are almost entirely unique. SAGE was itself built from the Croton et al. (2006) version of the Munich SAM — from which the likes of Gaea (Xie et al. 2020, and references therein) and L-galaxies (Henriques et al. 2020) also stem — with many versions preceding that too (Kauffmann et al. 1999; Springel et al. 2001; De Lucia, Kauffmann & White 2004). While the new Dark Sage that we present here brings together five years of development since its last release, it builds on over 25 years of research in the community.

Outside of the main model papers (Stevens et al. 2016, 2018), Dark Sage has been used to investigate the effects of environment and feedback on the gas content of galaxies (Stevens & Brown 2017), research the angular momentum of galaxies at both high redshifts (Okamura, Shimasaku & Kawamata 2018) and high stellar masses (Porras-Valverde et al. 2021), explain the origin of H i-excess galaxies (Lutz et al. 2018), confirm the origin of the H i size–mass relation (Stevens et al. 2019b), predict the results of H i intensity-mapping experiments (Wolz et al. 2022), develop an H i-dependent halo occupation distribution model for cosmology applications (Qin et al. 2022), and more.a The significant updates that we present in this paper will not only expand the range of use cases for Dark Sage, but also improve the predictive power for the use cases it was originally designed for.

We have endeavoured to adopt a new philosophy with how SAMs should be designed with this new version of Dark Sage. Our main goal is to re-derive the core aspects of galaxy evolution physics from as near to first principles as possible; there are fewer ad hoc descriptions left, and the ones we do use are updated to modern measurements where possible, often over broader mass and redshift ranges. This has meant increasing the coupling between physical prescriptions, which results in the elimination of a large number of free parameters. Simultaneously, we have increased the detail and complexity of the galaxies we evolve, building their structure numerically in more than one dimension. Feedback from stars and AGNs — arguably the most important processes in galaxy evolution — are derived from the law of energy conservation. While the model is certainly not free of phenomenology, our philosophy has meant we can capture the uncertainty in our model almost completely with only three free parameters (a factor 3–6 lower than comparable models tabulated in Table 2), each with well defined boundaries, which we successfully calibrate to observational statistics of galaxies with an objective, automated method.

Dark Sage is a publicly available, version-controlled code (Stevens et al. 2017b), open for the community to use and modify as they desire. The aim of this paper is to describe how the model works, and should serve especially well for those using the code itself. While we present some predictions of the model here, we deliberately do not linger on them. Science questions and specific results will be better served as the subject of future, focussed papers.

2. ‘The next generation’ N-body simulations

The results of Dark Sage in this paper are produced by running the model on the merger trees of the main gravity-only MillenniumTNG N-body simulation. The MillenniumTNG simulationsb (Barrera et al. 2023; Bose et al. 2023; Contreras et al. 2023; Delgado et al. 2023; Ferlito et al. 2023; Hadzhiyska et al. 2023a,b; Hernández-Aguayo et al. 2023; Kannan et al. 2023; Pakmor et al. 2023) are the latest addition to ‘The Next Generation’ (TNG) simulation suite, which began with IllustrisTNGc (Weinberger et al. 2017; Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019b; Pillepich et al. 2018a,b; Springel et al. 2018). Each simulation in the TNG suite has paired N-body and magnetohydrodynamic versions with identical initial conditions, run with the same code [arepo (Springel 2010), although N-body runs have also been completed with gadget-4 (Springel et al. 2021)], and all carry the same $\Lambda$ CDM cosmology: $\Omega_\mathrm{m} = 0.3089$ , $\Omega_\Lambda = 0.6911$ , $\Omega_\mathrm{b} = 0.0486$ , $h = 0.6774$ , $\sigma_8 = 0.8159$ , $n_\mathrm{s} = 0.9667$ (Plank Collaboration 2016). We assume this cosmology for all results in this paper (including observational data that we compare to, where we have adjusted published data accordingly). The MillenniumTNG run used in this paper (hereafter ‘MTNG’, referred to as ‘MTNG740-DM-1’ in other papers) has a comoving box length of $500\,h^{-1}$ cMpc (nearly 740 cMpc) with $4320^3$ particles of mass $1.329 \times 10^8\, h^{-1}\, \mathrm{M}_\odot$ ( $ = 1.962 \times 10^8\, \mathrm{M}_\odot$ ), and was run with arepo.

There are many advantages to pairing a SAM like Dark Sage with MTNG:

• • Both the volume and resolution are high enough to simultaneously capture a large number of galaxy clusters with small enough galaxies to study environmental effects on galaxy evolution.

• • Simulations with different mass resolutions have been run, allowing for future resolution studies.

• • Objects from the SAM can be directly compared against the galaxies formed in the hydrodynamical run (Pakmor et al. 2023) in future studies.

• • Galaxy catalogues from at least one other SAM have already been produced (Barrera et al. 2023, another descendant of the Munich model), allowing for additional future model comparison.

• • With an equivalent volume and eight times the number of particles, MTNG is a natural successor to the original Millennium simulation (Springel et al. 2005), which many SAM results in the literature are based on, including previous iterations of Dark Sage (Stevens et al. 2016, 2018).

• • If a galaxy of a given stellar mass being ‘resolved’ in a simulation requires it to have a minimum number of particles, then by virtue of the fact that only a fraction of baryon particles in a hydrodynamic-simulation halo will become star particles, semi-analytic galaxies built in a simulation of fixed volume and dark-matter particle number are, in principle, resolved (and therefore complete) to lower stellar masses than a hydrodynamic simulation with those same specifications.

Haloes and their substructure (‘subhaloes’) are identified in each snapshot of MTNG with subfind-hbt, described in section 7 of Springel et al. (2021). A minimum of 20 particles must be connected in a friends-of-friends (FoF) group to form a halo. When using the terms ‘halo mass’ or ‘virial mass’, we mean the mass enclosed by a sphere whose mean density is 200 times the critical density of the (simulated) universe (denoted by subscript ‘200c’). This has a one-to-one relationship with both virial radius and virial velocity:

(1) \begin{equation}M_\mathrm{200c} = \frac{100\, H^2(z)\, R_\mathrm{200c}^3}{G} = \frac{V_\mathrm{200c}^3}{10\, G\, H(z)},\end{equation}

where H(z) is the redshift-dependent fractional expansion rate of the universe or ‘Hubble parameter’.

The merger trees for MTNG were originally constructed with the new gadget-4 tree builder (see section 7.4 of Springel et al. 2021). This is similar to the format of lhalotree, originally presented in Springel et al. (2001) and used for the Millennium simulation (see the supplementary material of Springel et al. 2005) and is now in HDF5 format. Dark Sage requires input merger trees in the original binary format of lhalotree. We thus converted the MTNG trees back to the old lhalotree format for this project.

While the final results in this paper are presented with the full MTNG box, all preliminary results and model calibration were performed with a smaller simulation for practical purposes. ‘Mini-MTNG’, as we refer to it, had identical parameters, including mass resolution, but only included   $540^3$ particles in a $62.5\,h^{-1}$ cMpc box (MTNG’s equivalent of ‘Mini-Millennium’).

3. Definitions and design of baryonic reservoirs

Baryonic matter [gas, stars, and black holes (BHs)] in Dark Sage galaxies and haloes is broken into discrete components. Indeed, the primary function of a SAM is to calculate how mass and metals move between these components. For basic context and to familiarise the reader with the nomenclature used throughout this paper, we briefly describe each component in this section.

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 2022). 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. (2016) 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.

4. Baryon accretion, gas cooling, and disc formation

When a halo grows in mass — that is, when stepping between snapshots in the underlying merger trees — mass is added to the CGM (and/or occasionally the IHS; see below) of the central to ensure the baryonic mass fraction of the halo remains in line with the cosmic baryon fraction, $f_\mathrm{bary}^\mathrm{cosmic}$ , modulo a suppression factor, $f_\gamma^+$ , associated with photoionisation heating (Efstathiou 1992; Gnedin 2000). That is, we enforce

(5) \begin{align}f_\gamma^+\, f_\mathrm{bary}^\mathrm{cosmic}\, M_\mathrm{200c} & = m_\mathrm{ejec}^\mathrm{cen} + \sum_\mathrm{gal}^{R_\mathrm{gal}<R_\mathrm{200c}} \bigg( m_*^\mathrm{gal} + m_\mathrm{cold}^\mathrm{gal} \nonumber\\& \quad + m_\mathrm{CGM}^\mathrm{gal} + m_\mathrm{IHS}^\mathrm{gal} + m_\mathrm{BH}^\mathrm{gal} + m_\mathrm{IHBH}^\mathrm{gal} \bigg),\end{align}

where the sum over ‘gal’ includes the central galaxy and all satellites within the virial radius. There is an important subtlety in how the ejected reservoir of the central is handled in equation (5). The entire ejected reservoir is defined to lie outside $R_\mathrm{200c}$ . It is included in the right-hand side of equation (5) because that gas was inside $R_\mathrm{200c}$ ; this gas is ejected by stellar and AGN feedback (described in Sections 8 and 9), and there is no physical reason why that gas should be automatically replaced by other baryons (which would happen if it were not in the equation). In principle, this means that, even in the absence of photoionisation heating, the cosmic baryon fraction only provides an upper limit on the halo baryon fraction in Dark Sage. We can straightforwardly calculate the baryon fraction inside the virial radius of a halo as

(6) \begin{equation}f_\mathrm{bary}^\mathrm{200c} = f_\gamma^+\, f_\mathrm{bary}^\mathrm{cosmic} - \frac{m_\mathrm{ejec}^\mathrm{cen}}{M_\mathrm{200c}}.\end{equation}

The $f_\gamma^+$ factor is adopted from Dark Sage’s predecessors (Croton et al. 2006, 2016), which follow the ‘filtering mass’ fitting function of Kravtsov, Gnedin & Klypin (2004, see their appendix B). This depends on both halo mass and redshift and is primarily relevant during the epoch of reionisation ( $z \gtrsim 6$ , for all haloes) and for haloes with $M_\mathrm{200c} < 10^{12}\,\mathrm{M}_\odot$ (at all epochs). $f_\mathrm{bary}^\mathrm{cosmic}$ is set by the cosmology of the underlying N-body simulation ( $f_\mathrm{bary}^\mathrm{cosmic} = \Omega_\mathrm{b} / \Omega_\mathrm{m} = 0.1573$ in this work).

In previous versions of the model, when a halo grew, pristine cosmological gas (i.e. with close to zero metallicity) was added to the hot-gas reservoir. This still happens as the norm, but with our new framework, preference is first given to mass in the local intergalactic reservoirs to be accreted into the halo. Should a halo have a non-zero LIGM and/or LIGS, mass is proportionally and respectively transferred from these to the hot-gas reservoir and the IHS. Should this be insufficient to raise the halo baryon fraction to $f_\gamma^+\, f_\mathrm{bary}^\mathrm{cosmic}$ , then pristine cosmological gas is added to the hot gas to top it up. See Sections 11 and 12 for how the LIGM and LIGS are built up in the first place.

Gas must first pass through the CGM before reaching the ISM. The process for calculating how much mass is transferred from the CGM to the ISM in a given sub-time-step is similar in concept to previous incarnations of the model (which harken back to White & Frenk 1991) but is not identical.

Given that we assume hot gas to be spherically symmetric and homogenised, the cooling time for gas at a given radius in the CGM is calculated as

(7) \begin{equation}t_\mathrm{cool}(R) \equiv \frac{e_\mathrm{hot}^\mathrm{therm}}{\dot{e}_\mathrm{hot}^\mathrm{rad}(R)} = \frac{3 \, T_\mathrm{hot} \, \bar{\mu} m_\mathrm{p} \, k_\mathrm{B} }{2\, \rho_\mathrm{hot}(R)\, \Lambda(T_\mathrm{hot}, Z_\mathrm{hot})},\end{equation}

where $e_\mathrm{hot}^\mathrm{therm}$ is the specific thermal energy of the CGM, $\dot{e}_\mathrm{hot}^\mathrm{rad}(R)$ is the rate of specific energy loss through radiation, $T_\mathrm{hot}$ is the CGM temperature, $m_\mathrm{p}$ is the mass of a proton, $\bar{\mu} = 0.59$ under the assumption that hot gas is fully ionised, $k_\mathrm{B}$ is the Boltzmann constant, $\rho_\mathrm{hot}(R)$ is the hot-gas density profile, $\Lambda(T_\mathrm{hot}, Z_\mathrm{hot})$ is the tabulated cooling function from Sutherland & Dopita (1993), and $Z_\mathrm{hot}$ is the hot-gas metallicity. We maintain the assumption that the CGM temperature is the same as the virial temperature:d

(8) \begin{equation}T_\mathrm{hot} = T_\mathrm{CGM} = T_\mathrm{vir} \equiv \frac{\bar{\mu} m_\mathrm{p} \, V_\mathrm{200c}^2}{2\, k_\mathrm{B}}.\end{equation}

This formula for $T_\mathrm{CGM}$ is defined for centrals; after infall, satellites are assumed to have fixed $T_\mathrm{CGM}$ . We further maintain the definition of the ‘cooling radius’ as being that where the cooling time equals the halo dynamical time, that is,

(9) \begin{equation}t_\mathrm{cool}(R_\mathrm{cool}) = t_\mathrm{dyn} \equiv \frac{R_\mathrm{200c}}{V_\mathrm{200c}} = \frac{1}{10\,H(z)}.\end{equation}

However, what we have changed is the analytic profile for hot gas.

While previous iterations of SAGE and Dark Sage assumed the CGM to follow the density profile of a singular isothermal sphere, we instead now adopt the so-called ‘beta’ profile. Nominally, a beta profile follows the expression

(10) \begin{equation}\rho_{\beta}(R) = \rho_0 \left[1 + \left(\frac{R}{R_c}\right)^2 \right]^{-3\beta/2},\end{equation}

where $\rho_0 \equiv \rho_{\beta}(0)$ and $R_c$ is the ‘core’ radius. For the application of this profile to the galform SAM, modulo changes to the profile from feedback, Benson et al. (2003) assume that $\beta = 2/3$ and $R_c$ is a fixed fraction of the virial radius of a halo. We follow this idea, but instead of allowing for the profile shape to be modified by galaxy evolution processes, we set an explicit redshift dependence on the latter, that is, $R_c = c_\beta(z)\,R_\mathrm{200c}$ . Finally, with recognition that the CGM is defined in this paper to be within $R_\mathrm{200c}$ [i.e. a volume integral of $\rho_\mathrm{CGM}(R)$ out to $R_\mathrm{200c}$ must return $m_\mathrm{CGM}$ ], and that each component of the CGM is assumed to follow the same density profile, we implement

(11a) \begin{equation}\rho_\mathrm{hot}(R) = \frac{m_\mathrm{hot}}{m_\mathrm{CGM}}\, \rho_\mathrm{CGM}(R),\end{equation}

(11b) \begin{equation}\rho_\mathrm{CGM}(R) = \frac{m_\mathrm{CGM}\, \mathcal{C}_\beta(z)}{4 \pi \, c_\beta^2(z) \, R_\mathrm{200c}^3} \left[1 + \left(\frac{R}{c_\beta(z)\,R_\mathrm{200c}} \right)^2 \right]^{-1},\end{equation}

(11c) \begin{equation}\mathcal{C}_\beta(z) \equiv \left\{1 - c_\beta(z)\, \tan^{-1} \left[c_\beta^{-1}(z) \right] \right\}^{-1}.\end{equation}

For the functional form of $c_\beta(z)$ , we use the fitting function of Stevens et al. (2017a) to gas haloes from the EAGLEe simulations:

(12) \begin{equation}c_\beta(z) = \mathrm{max}\left[0.05, \,0.2\,\exp({-}1.5\,z) - 0.039\, z + 0.28 \right].\end{equation}

The lower limit of 0.05 is imposed by hand; this represents the 16th percentile of haloes at the highest redshift ( $z \simeq 4$ ) studied in Stevens et al. (2017a). As a point of comparison, in their models, Benson et al. (2003) and Font et al. (2008) adopt default (but not fixed) $c_\beta$ values of 0.07 and 0.1, respectively.

For the CGM of satellite galaxies, we still assume a truncation radius. While it does not make complete physical sense to define $R_\mathrm{200c}$ for a satellite subhalo, we adopt $M_\mathrm{200c}$ as the smaller of either the present total subhalo mass or its pre-infall virial mass, then solve for an equivalent $R_\mathrm{200c}$ through equation (1). This ensures that as satellites are tidally stripped, their CGM truncation radius decreases too, which is physically expected. By contrast, as noted above, $T_\mathrm{vir}$ is deliberately not recomputed after infall; if it were, we would be implicitly assuming that tidal stripping causes a decrease in the average CGM temperature, which seems non-physical.

Combining equations (7–12), one can solve for $R_\mathrm{cool}$ . Because of the update in the hot-gas density profile, this solution is different from earlier versions of Dark Sage.

Two modes of gas accretion onto a galaxy can take place; the ‘cold mode’ of accretion takes place when $R_\mathrm{cool} > R_\mathrm{200c}$ , otherwise the galaxy is in the ‘hot mode’ (White & Frenk 1991). In concept, these are unchanged from earlier models. That is, in the cold mode, gas in the CGM that is en route to the ISM is limited more by the time it takes to fall onto the ISM than any time required to cool. In this instance, it becomes somewhat of a misnomer to call the CGM ‘hot gas’, but we nevertheless maintain the nomenclature for consistency with the literature. For a sub-time-step $\Delta t$ , the mass deposited onto the ISM in the cold mode is taken as

(13) \begin{equation}\Delta m_\mathrm{hot \rightarrow cold}^{\operatorname{cold-mode}} = m_\mathrm{hot} \frac{\Delta t}{t_\mathrm{dyn}} - \Delta m_\mathrm{heat}^\mathrm{radio}.\end{equation}

The term $\Delta m_\mathrm{heat}^\mathrm{radio}$ represents an offset in accretion onto the ISM from radio-mode AGN feedback (defined and discussed in Section 9). While this term applies for both the cold and hot mode, it is normally only important in practice for the hot mode.

For the hot mode, the flux of mass cooling through $R_\mathrm{cool}$ is taken as equal to the mass deposition rate onto the galaxy (based on the results of Bertschinger 1989). To calculate this, we solve

(14) \begin{align}\Delta m_\mathrm{hot \rightarrow cold}^{\operatorname{hot-mode}} & = 4\pi \, \rho_\mathrm{hot}(R_\mathrm{cool})\, R_\mathrm{cool}^2 \left(\left. \frac{\mathrm{d} t_\mathrm{cool}}{\mathrm{d} R} \right|_{R \rightarrow R_\mathrm{cool}} \right)^{-1}\, \Delta t \nonumber\\& \quad - \Delta m_\mathrm{heat}^\mathrm{radio}.\end{align}

Assuming equation (11), the solution can be written as

(15) \begin{equation}\Delta m_\mathrm{hot \rightarrow cold}^{\operatorname{hot-mode}} = \frac{m_\mathrm{hot}\, R_\mathrm{cool}\, \mathcal{C}_\beta(z)}{2\, R_\mathrm{200c}\, t_\mathrm{dyn}}\, \Delta t - \Delta m_\mathrm{heat}^\mathrm{radio}.\end{equation}

With $\Delta m_\mathrm{hot \rightarrow cold}$ determined, the final thing to consider is the fraction of that mass that each annulus of the gas disc receives. As each annulus is bound by fixed values of specific angular momentum, this means defining a probability distribution function (PDF) of j for the cooling gas. We base this on the PDF of j for an exponential disc with a constant rotational velocity, that is,

(16) \begin{equation}\mathrm{PDF}_\mathrm{cool}(j) = \frac{4\, j}{j_\mathrm{cool}^2} \exp \left(\frac{-2\, j}{j_\mathrm{cool}} \right),\end{equation}

where $\int_0^\infty \mathrm{PDF}_\mathrm{cool}(j)\, \mathrm{d} j = 1$ by definition, and $j_\mathrm{cool}$ is the mean of the PDF. Equation (16) here is simply a more correct and flexible way of writing equation 6 of Stevens et al. (2016). We emphasise that even though equation (16) is the same as an exponential disc with a constant rotational velocity, that does not mean it is the only type of disc consistent with this PDF. That is, for any arbitrary surface density profile $\Sigma(r)$ , one can find a rotational velocity profile $v_\mathrm{circ}(r)$ that gives the same PDF of j as in equation (16). We are therefore not making any explicit assumptions about the form of $\Sigma(r)$ or $v_\mathrm{circ}(r)$ as gas cools/accretes onto a galaxy; in fact, we solve for these profiles (using additional information per Section 5). Even if one were to make assumptions about the surface density profile or velocity profile of freshly cooled gas, one would still end up with a predictive model for the resulting disc structure after several time-steps, as both the value of $j_\mathrm{cool}$ and the vector along which accretion takes places change between time-steps.

What then sets $j_\mathrm{cool}$ ? The answer to this also depends on the accretion mode. In the cold mode, gas is assumed to be efficiently transported from the halo outskirts to the galaxy via cold streams (as shown to occur in hydrodynamic simulations, e.g. Kereš et al. 2005; Dekel & Birnboim 2006; van de Voort et al. 2011). In the absence of evidence to the contrary, we set

(17) \begin{equation}j_\mathrm{cool}^{\operatorname{cold-mode}} = j_\mathrm{halo},\end{equation}

maintaining the status quo from several SAMs (e.g. Lagos et al. 2018; Henriques et al. 2020). $j_\mathrm{halo}$ is the specific angular momentum of the halo inside the virial radius, which we take from the input merger trees. For the hot mode, we modify the fitting function of Stevens et al. (2017a) — which was first implemented in Dark Sage in Stevens et al. (2018) — by rewriting it in terms of j:

(18) \begin{equation}\log_{10} \left( \frac{j_\mathrm{cool}^{\operatorname{hot-mode}}}{2\, R_\mathrm{200c}\, V_\mathrm{200c}} \right) = 0.23\, \log_{10}(\lambda) - k_d,\end{equation}

where we adopt the Bullock et al. (2001) approximation for halo spin,

(19) \begin{equation}\lambda = \frac{j_\mathrm{halo}}{\sqrt{2}\, R_\mathrm{200c}\, V_\mathrm{200c}}.\end{equation}

While the slope of equation (18) originates from equation 19 of Stevens et al. (2017a), the normalisation is set by hand to recover the result of Stevens et al. (2017a) that $\langle j_\mathrm{cool}^{\operatorname{hot-mode}} / j_\mathrm{halo} \rangle \simeq 1.4 \simeq \sqrt{2}$ when one only considers gas inside the virial radius. Rearranging the above, one can then solve for $k_d$ to force this outcome. Doing so gives

(20) \begin{equation}k_d \simeq \log_{10} \left\langle \lambda^{-0.77} \right\rangle.\end{equation}

To match the low-redshift sample used in Stevens et al. (2017a), we calculate the average of $\lambda^{-0.77}$ for haloes with $10^{11.5} \leq M_\mathrm{200c}/\mathrm{M}_\odot \leq 10^{12.5}$ at $z = 0$ in Mini-MTNG. This returns $k_d = 1.15$ .f Because of the explicit dependence on halo spin in equation (18) — which is known to have a scatter of $\sim$ 0.26 dex and only systematically depend on halo mass to second order in N-body simulations (Bullock et al. 2001; Knebe & Power 2008) — we naturally have a scatter in $j_\mathrm{cool}^{\operatorname{hot-mode}} / j_\mathrm{halo}$ in Dark Sage of $\sim$ 0.2 dex. This is shown as a function of halo mass in Figure 1. As can be seen by comparing the curves and shaded regions in the top panel for all centrals with those exclusively in the hot mode, $M_\mathrm{200c} \simeq 10^{12}\,\mathrm{ M}_\odot$ marks the transition between haloes predominantly in the cold mode and those in the hot mode. There are, nevertheless, some haloes in the hot mode all the way down to the resolution limit of the simulation.

Figure 1. Top panel: Host halo mass versus the ratio of the specific angular momentum of gas cooling from the CGM onto the ISM of the central galaxy in Dark Sage to that of the halo from MillenniumTNG at $z = 0$ . For haloes accreting in the ‘hot mode’, this ratio is determined by equation (18). The scatter and mild slope in this relation is driven entirely by the distribution of halo spins, as seen in the bottom panel. For cold-mode haloes, it is assumed that $j_\mathrm{cool} / j_\mathrm{halo} = 1$ . Running percentiles and means are shown for the full population (thinner lines in the foreground) and when exclusively selecting those in the hot mode (thicker, behind). The dotted line in the top panel marks where $j_\mathrm{cool} / j_\mathrm{halo} = 1.4$ , the average value predicted for hot-mode haloes by Stevens et al. (2017a). Satellite subhaloes are excluded from this plot.

With the above pieces in place, we can calculate the amount of cooling gas that goes into each annulus of the new gas disc. For the ith annulus, the mass received is

(21) \begin{equation}\Delta m_\mathrm{hot \rightarrow cold}^{(i)} = \Delta m_\mathrm{hot \rightarrow cold} \int_{j_\mathrm{inner}^{(i)}}^{j_\mathrm{outer}^{(i)}} \mathrm{PDF}_\mathrm{cool}(j)\, \mathrm{d} j.\end{equation}

If this is the first cooling episode in a halo or $m_\mathrm{cold}$ was 0 at the start of the sub-time-step, this completes the relevant description of how a gas disc is formed. If this new gas is adding to a pre-existing disc in the galaxy (as is the case for most galaxies at most snapshots), there is one more step to consider.

While the magnitudes of the cooling gas and halo’s specific angular momenta are not always assumed to be the same, their directions are. The orientation of a galaxy’s disc in Dark Sage is always tracked by its angular momentum vector. When a new gas disc forms, the direction of its angular momentum vector is set parallel to that of the halo, that is, $\hat{J}_\mathrm{cold} = \hat{J}_\mathrm{cool} = \hat{J}_\mathrm{halo}$ . However, $\hat{J}_\mathrm{halo}$ changes over time for any given halo, meaning sequential accretion episodes are unlikely to be parallel. When a pre-existing gas disc is present, a vector sum of $\vec{J}_\mathrm{cold}$ and $\vec{J}_\mathrm{cool}$ is done to define the new orientation of the resulting gas disc. Both the pre-existing and freshly cooled gas discs are then projected onto this new plane, which uses the same discretised j structure per equation (3a). This means that the mass in the first annulus of the resulting disc will be the sum of the first annuli of the pre-existing and freshly cooled discs plus some fraction of their second (and possibly their third, fourth, etc.) annuli, depending on the angles involved. This aspect of Dark Sage remains unchanged from Stevens et al. (2016).

It is worth mentioning that $\vec{J}_\mathrm{halo}$ is typically a noisy quantity in an N-body simulation (e.g. Benson 2017; Contreras, Padilla & Lagos 2017). While this induces some randomness to the precise orientation and magnitude of the angular momentum of cooling gas in a Dark Sage at a particular instant, this randomness should be averaged out after several time-steps with our method.

5. Computing potentials and converting from specific angular momentum to radius

Recall that disc annuli in Dark Sage galaxies are fixed in terms of specific angular momentum. While this design choice was deliberate, it nevertheless remains crucial to compute the radii of annuli for several of the model’s modules. Annulus radii must be (and are) regularly updated in the model, at least once per sub-time-step per galaxy. To convert j to r, one must first compute the potential of the galaxy and its (sub)halo in the plane of the disc, $\Phi(r)$ . In Stevens et al. (2016), this was achieved by approximating the potential as spherically symmetric and the disc to have precisely centripetal motion at all radii, then iteratively solving for $r_\mathrm{outer}^{(i)}$ for each annulus through

(22) \begin{equation}j_\mathrm{outer}^{(i)} = r_\mathrm{outer}^{(i)}\,v_\mathrm{circ}( r_\mathrm{outer}^{(i)}) .\end{equation}

Stevens et al. (2018) later dropped the assumption of perfect centripetal motion by adding a radially variant factor for the fraction of gravity balanced by circular motion (as opposed to random radial motion), $f_\mathrm{circ}(r)$ :

(23) \begin{equation}v_\mathrm{circ}^2(r) = f_\mathrm{circ}(r)\,v_\mathrm{cent}^2(r) = f_\mathrm{circ}(r)\, r\, \frac{\mathrm{d}\Phi}{\mathrm{d}r}, \end{equation}

where subscripts ‘circ’ and ‘cent’ are respective shorthands for ‘circular’ and ‘centripetal’. As noted in Stevens et al. (2018), the spherical-potential approximation is accurate in the outskirts of a galaxy, but becomes exponentially inaccurate towards the centre. Moreover, solving equation (22) iteratively was a bottleneck in the code. For these reasons we have overhauled the j-to-r conversion.

With the new version of the model, we treat the shape of the potential as a combination of spherically symmetric components (the dark-matter halo, stellar bulges, intracluster stars, and hot gas; shorthand ‘sph’) and two axisymmetric, infinitesimally thin discs (one each for stars and cold gas). This can be written in terms of contributions to the centripetal velocity:

(24) \begin{equation}v_\mathrm{cent}^2(r) = v_\mathrm{sph}^2(r) + v_\mathrm{*,disc}^2(r) + v_\mathrm{gas,disc}^2(r) .\end{equation}

The challenge comes in the implementation of the two lattermost terms. While one can numerically solve for a disc’s potential as a function of radius from first principles, the only way to self-consistently do this in the framework of Dark Sage is by maintaining the iterative j-to-r process and increasing the computational demand of each iteration. For practical reasons, we therefore employ an approximation of a different kind, but one that remains faithful to the general shape of a disc potential. Namely, we assume the analytic approximation of Obreschkow et al. (2009, their equation 37):

(25) \begin{equation}v_\mathrm{disc}^2(r) = \frac{G\, m_\mathrm{disc}}{r_\mathrm{s}} \Biggl[ 1 + 4.8 \exp \left({-}0.35\,\frac{r}{r_\mathrm{s}} - 3.5\,\frac{r_\mathrm{s}}{r} \right) \Biggr]\\\times \Biggl[ \frac{r}{r_\mathrm{s}} + \left(\frac{r_\mathrm{s}}{r}\right)^2 + 2\sqrt{\frac{r_\mathrm{s}}{r}}\Biggr]^{-1}.\end{equation}

This approximation works well for exponential discs with scale radius $r_\mathrm{s}$ . Clearly, it is less accurate for discs with non-parametric disc structure as in Dark Sage. Nevertheless, after testing, we found this to be the best compromise between accuracy and practicality.

Because there is no requirement for Dark Sage discs to be precisely exponential, we need to calculate an analogue for $r_\mathrm{s}$ to use in equation (25). For exponential discs, there is a one-to-one relation between the scale radius and any radius that encloses x per cent of the disc’s mass, $r_x$ . For each of $x = 10, 20, ..., 80, 90$ , we calculate the effective exponential scale radius that corresponds to the disc’s $r_x$ . We then set $r_\mathrm{s}$ to be the average of these nine scale radii. This is a small shift from Stevens et al. (2018), who only used $r_{50}$ and $r_{90}$ to inform $r_\mathrm{s}$ , which relied on an otherwise unnecessary and now non-existent free parameter. The calculation of $r_\mathrm{s}$ and subsequent application of equation (25) are separately done for each of the stellar and gas discs, fulfilling the right-hand side of equation (24).

The contribution to the potential from all other baryons and dark matter is then straightforwardly calculated as:

(26a) \begin{equation}v_\mathrm{sph}^2(r) = \frac{G\,m_\mathrm{sph}({<}r)}{r},\end{equation}

(26b) \begin{align}m_\mathrm{sph}({<}r) = 4 \pi \int_0^r \biggl(\rho_\mathrm{DM}(R) + \rho_\mathrm{CGM}(R) + \sum_{\otimes} \rho_{\otimes}(R) \biggl) R^2\, \mathrm{d}R, \end{align}

where $\otimes = \operatorname{i-bulge}, \operatorname{m-bulge}, \mathrm{IHS}$ . Each of the profiles that contributes to $M_\mathrm{sph}({<} r)$ is analytic. Per Stevens et al. (2016, appendix B), dark matter is assumed to follow an NFW profile (Navarro, Frenk & White 1996) with a baryon-influenced mass–concentration relation (Di Cintio et al. 2014; Dutton & Macciò 2014):

(27a) \begin{equation}\rho_\mathrm{DM}(R) = \frac{M_\mathrm{DM}^\prime\, c^2_\mathrm{NFW}(z)}{4 \pi\, R_\mathrm{200c}^2\, R} \biggl[1 + \frac{c_\mathrm{NFW}(z)\, R}{R_\mathrm{200c}} \biggl]^{-2} \\[4pt]\times \biggl[ \ln\bigl(1+c_\mathrm{NFW}(z)\bigl) - \frac{c_\mathrm{NFW}(z)}{1 + c_\mathrm{NFW}(z)} \biggl]^{-1},\end{equation}

(27b) \begin{align}c_\mathrm{NFW}(z) & = \Biggl(1 \, + 3 \times 10^{-5}\, \exp\biggl\{ 3.4 \biggl[ \log_{10}\left(\frac{m_*}{M_\mathrm{200c}}\right)\nonumber\\& \qquad\qquad\qquad +\, 4.5 \biggr] \biggr\} \Biggr)\, c_\mathrm{DMO}(z),\end{align}

(27c) \begin{equation}c_\mathrm{DMO}(z) = 10^{a_\mathrm{NFW}(z)} + \left(\frac{M_\mathrm{200c}\, h}{10^{12}\, \mathrm{M}_{\odot}}\right)^{b_\mathrm{NFW}(z)},\end{equation}

(27d) \begin{equation}a_\mathrm{NFW}(z) = 0.520 + 0.385\, \exp\left({-}0.617\, z^{1.21}\right),\end{equation}

(27e) \begin{equation}b_\mathrm{NFW}(z) = -1.01 + 0.026\, z.\end{equation}

In principle, the dark-matter mass of the halo should be readily calculated as $M_\mathrm{DM} = \left(1 - f_\mathrm{bary}^\mathrm{200c}\right) M_\mathrm{200c}$ (equation 5). But, again for the sake of computation efficiency, we approximate the matter bound in subhaloes to also be distributed according to the same NFW profile as the dark matter (despite knowing that this is not the most accurate description of how satellites are distributed in Dark Sage — see Qin et al. 2022). We therefore effect

(28) \begin{align}M_\mathrm{DM}^\prime \equiv M_\mathrm{200c} - \big( m_* + m_\mathrm{cold} + m_\mathrm{CGM} + m_\mathrm{IHS} + m_\mathrm{BH} + m_\mathrm{IHBH} \big).\end{align}

Stellar spheroids (the instability-driven bulge, merger-driven bulge, and intrahalo stars) all follow Hernquist (1990) profiles that are truncated at the virial radius:

(29) \begin{equation}\rho_{\otimes}(R) = \frac{m_{\otimes}\, a_{\otimes}\, (R_\mathrm{200c} + a_{\otimes})^2}{2\pi\, R\, (R + a_{\otimes})^3\, R_\mathrm{200c}^2}.\end{equation}

While we have not changed the functional form of the stellar-spheroid profiles, we now calculate each of their scale radii, $a_{\otimes}$ , rather than assuming a scaling relation as done previously. These are discussed in the appropriate sections for each component later in the paper (see equations 71 & 121).

For the final piece of equation (23), we assume a function of the form

(30) \begin{equation}f_\mathrm{circ}(r) = 1 - \exp \left({-}r / r_{\lambda *} \right).\end{equation}

In Stevens et al. (2018), the scale radius $r_{\lambda *}$ took the form $r_{\lambda *} = r_\mathrm{s,*} / \left\{ 3 \left[1 - m_\mathrm{*,bulge}/\left(m_* + m_\mathrm{cold} \right) \right] \right\}$ ,g based loosely on the observational results of Bellstedt et al. (2017). Instead, we now calculate $r_{\lambda *}$ in a more self-consistent way, but maintain the assumption from Stevens et al. (2018) that $f_\mathrm{circ}$ should scale with the local stellar spin parameter, which is defined for an annulus as

(31) \begin{equation}\lambda_*^{(i)} \equiv \frac{ \big\langle v_\mathrm{circ}^{(i)} \big\rangle }{ \sqrt{\big\langle v_\mathrm{circ}^{(i)} \big\rangle^2 + \sigma_*^{(i)\,2}} }.\end{equation}

First, we initialise $r_{\lambda *}$ as $ r_\mathrm{s,*} / 3$ . Then, for every subsequent sub-time-step, we update $r_{\lambda *}$ to

(32) \begin{equation}r_{\lambda *} = \frac{- r\{\lambda_* = 0.5\}}{\ln(0.5)},\end{equation}

where $r\{\lambda_* = 0.5\}$ is smallest radius where $\lambda_*$ rises to 0.5 (interpolated from the annular $\lambda_*^{(i)}$ values). The factor of $-1/\ln(0.5) \simeq 1.44$ translates the factor-of-two scale radius into an exponential scale radius appropriate for equation (30).

With all of the above, we can analytically calculate $v_\mathrm{circ}(r)$ — and therefore j — for any r. To go in the opposite direction, that is, from j to r for annuli, we calculate j(r) for an array of r values. Using this array, we interpolate the j boundaries of the disc annuli to find their corresponding r.

We approximate the surface density of a given disc component to be constant within an annulus, that is,

(33) \begin{equation}\Sigma^{(i)}_\mathrm{disc} \equiv \frac{m^{(i)}_\mathrm{disc}}{\pi \left( r_\mathrm{outer}^{(i)\,2} - r_\mathrm{inner}^{(i)\,2} \right)}.\end{equation}

With this, it follows that the mass in an annulus has an average radius of

(34) \begin{equation}\langle r^{(i)} \rangle \equiv \sqrt{ \frac{ r_\mathrm{inner}^{(i)\,2} + r_\mathrm{outer}^{(i)\,2} }{2} } .\end{equation}

To be consistent with equation (4), we define the average circular velocity of an annulus as

(35) \begin{equation}\big\langle v_\mathrm{circ}^{(i)} \big\rangle \equiv \frac{\left\langle j^{(i)} \right\rangle}{\langle r^{(i)} \rangle} .\end{equation}

Finally, for simplicity and convenience, we define an approximate average potential of each annulus as

(36) \begin{equation}\left\langle \Phi^{(i)} \right\rangle \equiv \frac{\Phi\!\left(r_\mathrm{inner}^{(i)}\right) + \Phi\!\left(r_\mathrm{outer}^{(i)}\right)}{2} .\end{equation}

To demonstrate the resolved kinematics of galaxy discs in Dark Sage, we plot circular velocity profiles and stellar velocity dispersion profiles for random example galaxies across a range of stellar masses in Figure 2. A description of how $\sigma_\mathrm{*,disc}$ builds up over time can be found in Section 7.4.

Figure 2. Circular velocity profiles (top) and stellar velocity dispersion profiles (bottom) of Dark Sage discs at $z = 0$ . We sample one random central galaxy with $\mathrm{sSFR} = 10^{-10.5}\,\mathrm{yr}^{-1}$ and $m_\mathrm{*,bulge} < 0.3\,m_*$ for each stellar mass in intervals of 0.1 dex between $10^{8.5}$ and $10^{11}\,\mathrm{M}_\odot$ with a 0.04-dex search window. The colour of the curves transitions from dark blue to bright yellow as we move up in mass.

One relatively fundamental scaling relation of galaxies is the baryonic Tully–Fisher relation (McGaugh et al. 2000), a tighter relation than in its original stellar form of Tully & Fisher (1977). This relates the peak rotational velocity of a galaxy (a proxy for its potential) with its baryonic mass (in principle, another proxy for potential, if baryonic mass and dark-matter mass scale) for rotation-supported (disc-dominated) systems. Given that we have built the velocity curves of Dark Sage galaxies from their potential and the stellar and H i masses of galaxies are constrained by their observed mass functions at $z = 0$ (Section 13), the baryonic Tully–Fisher relation should fall out of Dark Sage as a natural prediction. As we show in Figure 3 by comparing to data from the SPARCh sample of galaxies (Lelli et al. 2019), it does! Though, we caution that our sample selection to match the SPARC is crude, intended only as a proof of concept here. The galaxies masses of the SPARC sample also have increasingly larger systematic uncertainties due to the adoption of a morphology-independent mass-to-light ratio (see Lelli et al. 2016, 2019). For these reasons, we should not expect a perfect match. For example, the scatter at low masses for Dark Sage of $\sim$ 0.02 dex appears low compared to the intrinsic scatter of $0.040 \pm 0.006$ dex in SPARC (cf. Lelli et al. 2019), but we found in testing other cuts that the scatter in Dark Sage varied up to 0.03 dex.

Figure 3. The Baryonic Tully–Fisher relation for Dark Sage galaxies: that is, the maximum rotational velocity of each galaxy as a function of its stellar + neutral-gas mass. Hex bins show the number density of Dark Sage galaxies with $m_\mathrm{*,bulge} < 0.1\,m_*$ , $m_{\rm H\,\small{I}} \geq 10^8\,\mathrm{M}_\odot$ , and $m_\mathrm{neutral}/m_* \geq 10^{-2}$ . The running solid lines are the median (thick) and 16th and 84th percentiles (thin) of the Dark Sage selection, calculated in bins along the y-axis. The Dark Sage cuts are a crude attempt at a comparable selection to the SPARC sample of galaxies (Lelli et al. 2019), compared here.

6. Gas chemistry and vertical structure of discs

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. (2019):

(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. (2019). Results from Pillepich et al. (2019) 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. 2018). 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. 2019). 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. (2021, 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 (2011) — 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. (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. 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 (2020). 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. (2016), we use the equations from Blitz & Rosolowsky (2004) and Elmegreen (1989):

(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 changei from Stevens et al. (2016) 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 (2010) and Gnedin & Draine (2014). 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. 2017; Catinella et al. 2018). 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. (2019a, 2021).

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. 2018), xCOLD GASS for H $_{2}$ (Saintonge et al. 2017), and SDSS for sSFR (the volume-limited sample used in Brown et al. 2017). 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. Disc instabilities and star formation

The one physical process that arguably is the linch-pin of the entire Dark Sage model is gravitational disc instabilities. The handling of said instabilities drives the majority of star formation, radial redistribution of disc mass, disc heating, bulge growth, and BH growth. While the triggering of an instability in the new Dark Sage is very similar to that of Stevens et al. (2016), some important changes have been made to how instabilities are resolved.

Previously in Dark Sage, stars could form through one of three possible channels. The ‘passive’ channel prescribed each annulus to form stars at a rate proportional to its local H $_{2}$ content. The ‘instability’ channel represented stars formed in bursts driven by local gravitational instabilities. The ‘merger’ channel accounted for extra star formation triggered by gas compression during a galaxy merger. While the former two were physically well defined to operate in a single annulus, the prescription for merger-driven starbursts was an ad hoc, purely phenomenological formula, designed originally to work on global galaxy scales instead (Somerville, Primack & Faber 2001): see section 3.9.3 of Stevens et al. (2016). Because this goes against the ethos of the new model, we have removed the merger channel for star formation. We maintain the passive and instability channels. Elimination of the merger-burst channel does not prevent elevated star formation from taking place. With how mergers are resolved (Section 12), the sudden addition of gas to the descendant’s ISM lowers its stability, making it likely that additional star formation occurs through instabilities directly after the merger.

For more than the above reason, instabilities play a more central role in Dark Sage star formation than they used to (despite their already having been the most important channel for many galaxies, as shown in figure 1 of Stevens & Brown 2017). While passive star formation still occurs, we make an important change in calculating this after resolving instabilities, rather than before. This means that passive star formation is now exclusive to the gas that is stable on the scale of an annulus but still has molecular clouds present. In hindsight, Dark Sage should have been originally set up this way.

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 1964; Binney & Tremaine 1987; Pringle & King 2007). 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 2011). 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 2003). 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 2004, 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. 2016). 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 2013). 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. (2013). 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. 2011) 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. (2016), $\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. (2011). 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. (2011) as well as data from the HERACLES and VERTICO surveys (Leroy et al. 2009; Brown et al. 2021, 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. 2023).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. New stellar feedback model

We introduce a brand-new model for stellar feedback in Dark Sage. Gone are the days of explicitly prescribing a mass loading factor that may or may not be modulated by local surface density. Instead, we now use a purely energy-based argument to calculate the amount of reheated (and ejected) gas from each episode of star formation, in each annulus. Gone also are the assumptions of instantaneous recycling, where we have now implemented a delayed, self-consistent scheme for stellar mass loss, metal enrichment, and associated feedback, which we outline below.

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 (2003) 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 1988). 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 (2014) 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 (2014). 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 1931). 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 ISMk 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 (2003) IMF combined with a Conroy, Gunn & White (2009) simple stellar population model (Lagos et al. 2018). 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. (1995) 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. (2016) 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 2009). 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. 2013), 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. 2023). 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. 2021), 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).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 2013). 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. 2019a 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. (2017), as they are for $z < 0.2$ galaxies from Heckman et al. (2015), while the points from Sugahara et al. (2017) 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. 2015; Chisholm et al. 2017; Sugahara et al. 2017). Based on the captions of tables 1 and 2 of Heckman et al. (2015), 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. 2011; Bower, Benson & Crain 2012; Croton et al. 2016; Lagos et al. 2018). 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.m The assumed functions used in Bower et al. (2012) and Lagos et al. (2018) are similar in slope to what Dark Sage predicts, but systematically higher and lower, respectively. Meanwhile, the asymptotic function of Guo et al. (2011) and the constant used by Croton et al. (2016) appear both qualitatively and quantitatively inconsistent with our energy-conserving framework. The comparison with Lagos et al. (2018) is particularly interesting, as their fitting function was based on the predictive model of Lagos et al. (2013). There are many differences between Dark Sage and the Lagos et al. (2013) 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 2018), 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. (2004) and Bellstedt et al. (2021).

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. (2004) at $z \simeq 0.1$ , covering their 16–84th interpercentile range, and the running median from Bellstedt et al. (2021) 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. (2021), 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. (2004), 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 2008) 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. 2009) (we have corrected the Tremonti et al. 2004 relation for a Chabrier 2003 IMF and our cosmology). For similar reasons, in their comparison of the mass–metallicity relation of the IllustrisTNG simulation with observations, Torrey et al. (2019) ignore the normalisation of the Tremonti et al. (2004) 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. 2021). 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. 2014), the ‘Mapping nearby Galaxies at Apache Point Observatory’ survey (MaNGA; Belfiore et al. 2017), and the ‘Sydney-AAO Multi-object Integral-field spectrograph’ survey (SAMI; Poetrodjojo et al. 2021), which have been homogenized (Sharda et al. 2021; 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.),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. 2016). This is not entirely by coincidence; part of the inspiration for our model design came from Mutch et al. (2016). Similar designs of delayed feedback and chemical enrichment are also incorporated into the galacticus model (Benson 2012) and in De Lucia et al. (2014), 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. Active galactic nuclei

We maintain a two-mode model for feedback from active galactic nuclei (AGNs) in Dark Sage, but have changed the way both modes operate.

9.1. Radio mode

In earlier versions of Dark Sage, the model for radio feedback followed Croton et al. (2016). In this mode, hot gas is passively accreted directly onto the BH, based on the BH’s Bondi (1952) accretion rate. Combined with the ‘maximal cooling flow’ model of Nulsen & Fabian (2000), 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. 2006), where $m_\mathrm{BH}$ is the BH’s mass. In Croton et al. (2016), 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 (1926) 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 1969; Bardeen 1970; Bardeen, Press & Teukolsky 1972; Novikov & Thorne 1973; Shakura & Sunyaev 1973; Stevens 2015). 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. 1972). While $\epsilon_\mathrm{AGN} = 0.0572$ for a non-rotating Schwarzschild (1916) BH, $\epsilon_\mathrm{AGN}$ can theoretically range from 0.0378 to 0.3210, given the anti-parallel and parallel cases for a Kerr (1963) BH rotating at the Thorne (1974) 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. 2011). 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. (2006) and Bower et al. (2006). 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. 2017, 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. (2016), 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. 2018).

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. (2016) 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 2000; Croton et al. 2016) 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 (1926) 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. (2016) 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. The galactic fountain and gas reincorporation

We introduce a new treatment for the galactic fountain, which affects how gas is reincorporated into the hot-gas reservoir to become available for cooling again after being heated by feedback. Our method is new not only for Dark Sage but for SAMs in general.

Early in their development, SAMs assumed the reincorporation timescale to scale linearly with the halo dynamical time (equation 9) (Springel et al. 2001), with the proportionality constant often treated as a free parameter (per De Lucia et al. 2004). Motivated by the difficulty that this simple model posed in recovering observations at multiple epochs (Henriques et al. 2013; Mutch, Poole & Croton 2013), it later became commonplace for reincorporation timescales to have a direct dependence on virial mass or velocity (with various functional forms — e.g. Henriques et al. 2013; White, Somerville & Ferguson 2015; Croton et al. 2016). However, none of these approaches are in line with our philosophy of having a self-consistent feedback framework that minimises the number of free parameters. Instead, we introduce a framework where we update timescales for mass transfer between the CGM and ejected reservoir self-consistently with the energetics of feedback.

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. 2022; also see Bellstedt et al. 2020).

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. Environmental stripping of gas and stars

When a subhalo becomes the satellite of a host halo, several rules for its evolution change. Most notably, baryon reservoirs become prone to environmental effects that can deplete their contents. In Dark Sage, we now account for the effects of starvation, strangulation, tidal stripping, and ram pressure stripping.

In principle, central haloes can also be affected by the environment of nearby large-scale structure. For example, a halo falling into a cosmic filament from a void should experience tidal and ram pressure effects from that filament. In this section, we describe how environmental effects are handled for both central and satellites galaxies in Dark Sage.

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. 2013; Srisawat et al. 2013; Elahi et al. 2019a,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. 2022).

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,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. (2016).

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 (2017). Other criteria for the same phenomenon have been used in other SAMs (e.g. Cora et al. 2018).

We approximate the tidal radius of a satellite galaxy based on the formula of Tormen, Diaferio & Syer (1998, 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}$ ).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 1972). 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. Galaxy mergers

Similar to Croton et al. (2016), a galaxy merger or disruption event occurs in Dark Sage when the baryon fraction of a subhalo rises above a threshold value, $f^\mathrm{ thresh}_\mathrm{bary}$ . In previous versions of SAGE and Dark Sage, the subhalo baryon fraction only considered the cold-gas and stellar mass of the satellite galaxy it housed. We now further include the CGM, IHS, and BH masses of the satellite in this calculation. Where previously a threshold fraction of 1.0 was used, we now impose a threshold fraction of 0.5. Recall that the subhalo mass is the total sum of all particles associated with that object from the underlying N-body simulation. Further recall that that mass represents the sum of both dark matter and baryons, but that all matter in the N-body simulation behaves like dark matter. As a subhalo is tidally stripped, both its dark-matter and hot-gas mass decline. But the stellar and cold-gas content of a satellite typically take longer to be reduced by environmental effects (Section 11). This means the baryon fraction of a subhalo tends to get larger the longer it is a satellite. In principle, it could continue to rise to arbitrarily large numbers greater than 1 — for example, when a subhalo composed of thousands of particles at infall is continually stripped down to the minimum of 20 particles — unless we impose an extra limit. Our physical argument for triggering a merger or satellite disruption event once a subhalo’s baryon fraction rises to 0.5 is that the implicit assumption in the N-body simulation — that the satellite galaxy behaves like a dark-matter-dominated object — no longer holds. It presumably must have reached the point where either the un-modelled influence of baryons on the satellite’s orbit should have caused it to merge with the central galaxy or the stellar and cold-gas content of the galaxy should have been tidally stripped.

The above means that satellite galaxies are generally merged or disrupted before their subhalo is ‘lost’ in the merger trees of the N-body simulation (their being lost in the merger trees triggers a merger or disruption in the same fashion, but these occasions are rare). This also means there are no ‘orphan’ galaxies in Dark Sage.

Arguably, it would be better to include orphan galaxies (i.e. detach satellite galaxies from their host subhalo and let them continue to evolve until different criteria are met for merging or disruption), but this requires changes to Dark Sage’s framework beyond the scope of this version of the model. In practice, this means that the intrahalo and local intergalactic stellar and gas components of Dark Sage centrals include the content that other models might have otherwise associated with orphan galaxies.

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. (2016), this predicted timescale was based on the dynamical friction model of Binney & Tremaine (1987). However, Poulton et al. (2020) 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. (2021) 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. (2021) 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 (1990) 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. (2016), 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. (2016).

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.

13. Calibration of the model

Our new version of Dark Sage stands out from most — possibly all — other SAMs of galaxy formation in its small number of free parameters. We only allow three parameters to be calibrated in the model:

• • $f_\mathrm{move}^\mathrm{gas}$ , the fraction of Toomre-unstable gas that moves to adjacent annuli, as opposed to being consumed in a starburst with feedback (Section 7). This strongly influences star formation rates, the structure and morphology of galaxies, and the rate at which BHs grow. To our knowledge, no other SAM has a parameter that is directly analogous to $f_\mathrm{move}^\mathrm{gas}$ .

• • $\mathcal{E}_\mathrm{SN}$ , the average energy per supernova that is coupled to the ISM/CGM (Section 8). This sets the normalisation for the strength of stellar feedback, which other models usually rely on additional free parameters for.

• • $\epsilon_\mathrm{AGN}$ , the radiative efficiency of BHs, which controls the strength of AGN feedback (Section 9). Again, most other models have more than one parameter to control the overall strength of AGN feedback.

Our three free parameters are calibrated against three observables:

• • the $z = 0$ stellar mass function (SMF) from Driver et al. (2022);q

• • the $z = 0$ H i mass function (HIMF) from Jones et al. (2018);

• • and the CSFH from D’Silva et al. (2023).

The best-fitting values for our parameters, following the method we describe below, are provided in Table 1.

Table 1. Search range of Dark Sage’s free parameters for the PSO calibration and the output best-fitting values to two significant figures. The latter are used in all results throughout this paper. There is no prior expectation on the value of $f_\mathrm{move}^\mathrm{gas}$ ; as it is defined to have a value in (0,1), its search range conservatively covers almost all possible values. $\mathcal{E}_\mathrm{SN}$ is expected to have a value close to $10^{44}$ J, based on literature canon, and the search range around this value is deliberately conservatively wide. The search range of $\epsilon_\mathrm{AGN}$ matches the maximum possible range for a black-hole accretion model where mass moves quasi-statically through a centripetal disc until reaching the innermost stable circular orbit, given the range of possible black-hole spin values.

We invoke a particle-swarm optimisation (PSO) method to calibrate Dark Sage. PSO was introduced by Kennedy & Eberhart (1995) and first adopted to calibrate a SAM by Ruiz et al. (2015). Using PSO, we find the parameter set that provides the minimum multiplicative reduced $\chi^2$ across the three observables. We deliberately multiply these together instead of adding, as this circumvents issues of the quoted errors carrying different meaning across the constraints, and it ensures that an improvement in one reduced $\chi^2$ of a given factor is treated the same as any other constraint. That is, we minimise

(129) \begin{equation}\chi^2_\mathrm{joint} \equiv \prod_\mathrm{var} \left(\frac{1}{N_{p,\mathrm{var}} - 3} \sum_{p = 1}^{N_{p,\mathrm{var}}} \left[ \frac{y_{p,\mathrm{var}}^\mathrm{obs} - y_{p,\mathrm{var}}^\mathrm{mod}}{\sigma_{p,\mathrm{var}}^\mathrm{obs}} \right]^2 \right), \end{equation}

where var = SMF, HIMF, CSFH. Each y (evaluation) and $\sigma$ (uncertainty in y) term for point p is calculated in log-space. We match bin sizes in Dark Sage to the observational data where possible (SMF, HIMF), and interpolate Dark Sage’s inherent bins otherwise (CSFH), to ensure the fairest comparison of every point. Note that we do not include an uncertainty for Dark Sage points in equation (129). In principle, we could have applied a Poisson statistic here (as done by Ruiz et al. 2015), but this would lead to down-weighting points near and after the knee of the SMF and HIMF, which are arguably the most important — and, in practice, the hardest — parts of the mass functions to get right. Counting uncertainties are also already incorporated into the observational uncertainties.r

The PSO code makes use of the pyswarm package for python and acts as a wrapper around Dark Sage. pyswarm implements a PSO technique similar to that discussed by Ruiz et al. (2015). A core difference, however, is that the initial particle positions are randomly drawn from a uniform distribution corresponding to the search space boundary conditions (central column of Table 1), rather than a Latin Hypercube configuration. The PSO code we use was originally developed for the Shark SAM, with a full description due to be published elsewhere. Minimal edits were made from the PSO codebase in the Shark repositorys to allow the code to interface with Dark Sage.

Once the initial conditions are set, the particle positions at the subsequent iteration are computed according to

(130a) \begin{equation}\vec{\mathcal{P}}_m^{[k+1]} = \vec{\mathcal{P}}_m^{[k]} + \vec{\mathcal{V}}_m^{[k]},\end{equation}

(130b) \begin{equation}\vec{\mathcal{V}}_m^{[k]} = \frac{1}{2}\vec{\mathcal{V}}_m^{[k-1]} + \xi_1 \left( \vec{\mathcal{P}}_m^{[b]} - \vec{\mathcal{P}}_m^{[k]} \right) + \xi_2 \left( \vec{\mathcal{P}}_q^{[l]} - \vec{\mathcal{P}}_m^{[k]} \right).\end{equation}

$\vec{\mathcal{P}}_m^{[k]}$ represents the mth particle’s position in the $(f_\mathrm{move}^\mathrm{gas},\,\epsilon_\mathrm{AGN},\,\mathcal{E}_\mathrm{SN})$ parameter space at iteration k. $\vec{\mathcal{V}}_m^{[k]}$ is the particle’s corresponding ‘velocity’. $\vec{\mathcal{P}}_m^{[b]}$ represents particle m’s position at the iteration of best fit (lowest $\chi^2_\mathrm{joint}$ ) across all previous iterations ( $b \leq k$ ). $\vec{\mathcal{P}}_q^{[l]}$ represents the position for the global best fit for all particles across all previous iterations (that of particle q at iteration l). $\xi_1$ and $\xi_2$ are random numbers, each drawn from a uniform distribution in $[0,0.5]$ . The initial velocity of each particle, $\vec{\mathcal{V}}_m^{[0]}$ , is set to a random fraction ( $< 1$ ) of the length of the parameter space in each dimension, and also in a random direction. Should the velocity of a particle move it to outside the search range, it is manually shifted to the edge of the search range.

We must also specify the range in the observational data to which the model is fitted. We use the simulation’s resolution to inform this decision. Our loose criterion is that if galaxies of a given stellar or H i mass exclusively occupy (sub)haloes that have had at least 200 particles for one or more snapshots throughout their primary progenitor branch in the merger trees, then we are ‘complete’ at that mass, and can fit observations from there upwards. 200 particles is a semi-arbitrary number to require completeness to, but in the context of MTNG haloes having a minimum of 20 particles, our criterion strikes a balance between using as much information (as many haloes) as possible without giving too much trust towards MTNG’s resolution limit. The stellar or H i mass to which we are complete to 200 particles depends on the parameters of the model. After exploring the parameter space, we found we could comfortably rely on being complete for $m_* \geq 10^9\,\mathrm{M}_\odot$ and $m_{\rm H\,\small{i}} \geq 10^9\,\mathrm{M}_\odot$ (erring on the side of caution, especially in the case of the stellar mass). These are the limits we adopt for the SMF and HIMF when calculating $\chi^2_\mathrm{joint}$ . In Figure 13, we show stellar and H i mass distributions of Dark Sage galaxies as a function of their historical-maximum (sub)halo mass at $z = 0$ ,t and the 200-particle completeness as a function of each mass.

Figure 13. Historical maximum of the FoF or subhalo mass for Dark Sage galaxies versus their stellar mass (top panel) or H i mass (bottom panel) at $z = 0$ . We define a galaxy as having 200 times the particle mass of MTNG on the x-axis as being resolved. Calibration to the SMF and HIMF only occurs (approximately) above the respective masses where Dark Sage is complete to this mass resolution. The long-dashed, red curves show the completeness as a function of stellar and H i mass. The shaded region on the left is where galaxies are not reliably resolved by this definition. The thin, dashed, horizontal lines are the masses above which we calibrate the SMF and HIMF. As in other figures, hexbins and purple lines show the density of Dark Sage galaxies and the running 16th, 50th, and 84th percentiles.

In total, the PSO code went through two rounds of 16 averageu iterations with 16 particles, meaning 512 parameter sets of Dark Sage were fully simulated with mini-MTNG. Each instance of Dark Sage on mini-MTNG required $\sim$ 2.2 CPUh of walltime. It was important to do more than one run of the PSO code from start to finish to mitigate the possibility of finding a local minimum instead of the global minimum within the defined parameter space. Figure 14 shows the output of each run in the PSO chains against the calibration data, the relative goodness of fit of each run, and the particle positions in the model’s parameter space. We rounded the best-fitting parameters to two significant figures for the final run of the model on MTNG.

Table 2. Stand-out features of Dark Sage compared to its previous version and other semi-analytic models in the literature. Each row represents a different semi-analytic model. Columns represent: (i) the number of free parameters that model has; (ii) the model’s treatment of radial disc structure; (iii) treatment of radial metallicity variation in galaxy discs, if any; (iv) treatment of stellar-age distributions as function of disc radius; (v) treatment of mass loading from stellar feedback. Only if an entry is ‘non-parametric’ can that model make predictions for that property. Note that the total number of nominal parameters a model may have is not what we equate as free parameters. For a parameter to count as ‘free’ by our definition, it has to have been described as being varied in a calibration procedure in the model paper. The most relevant reference for each model is: Dark Sage 2018 (Stevens et al. 2018), L-galaxies (Henriques et al. 2020), Meraxes (Mutch et al. 2016), SAG (Cora et al. 2018), SAGE (Croton et al. 2016), Shark (Lagos et al. 2018).

Figure 14. Three observed relations used to constrain the three free parameters in Dark Sage. Top left: stellar mass function at $z = 0$ with observational data at $z \simeq 0.072$ from Driver et al. (2022). Top middle: H i mass function at $z = 0$ with observational data from Jones et al. (2018). Top right: cosmic star formation history with observational data from D’Silva et al. (2023). Each line in each of the panels corresponds to a run of Dark Sage with a different $(f_\mathrm{move}^\mathrm{gas},\,\epsilon_\mathrm{AGN},\,\mathcal{E}_\mathrm{SN})$ parameter set. The colour indicates the combined reduced $\chi^2$ of the three constraints (equation 129). The bottom panels show where each run lies in Dark Sage’s parameter space. The best fit was found by applying a particle-swarm optimisation code. See Section 13 for further details.

We note that Dark Sage can reproduce any of these three observables to higher accuracy if constrained exclusively by that dataset. In principle, we could also have weighted the constraints in the joint $\chi^2$ calculation had we felt one was more important to reproduce. But in the absence of an obvious reason why one constraint should have greater inherent value over another, we defaulted to having equal weights. Evidently, there is some tension in the model with respect to simultaneously reproducing the SMF and HIMF. This could reflect (i) the simplicity of the molecular-fraction calculation in the model (Section 6.4); (ii) a deeper limitation of the model’s design; (iii) uncertainties in the observational data, bearing in mind that no forward-modelling of any observational uncertainties has been done here; or a combination thereof. We leave a deeper assessment of this for future studies.

To highlight the significance of having so few free parameters, especially in the context of Dark Sage galaxies having multidimensional structure, we compare the number of free parameters and details of galaxy structure with other popular SAMs in the literature in Table 2. We count 15 free parameters for L-galaxies based on the parameters listed that are not fixed in table 1 of Henriques et al. (2020), that is, used in their Monte Carlo Markov Chain calibration procedure; 10 for Meraxes based on table 1 of Mutch et al. (2016), discounting the one parameter listed with a value of zero; 9 for SAG listed in table 1 of Cora et al. (2018) in their PSO calibration; 8 for SAGE based on inside knowledge and what is described as explicitly being calibrated/free in Croton et al. (2016); and 18 for Shark after cross-referencing the parameters listed in table 3 of Lagos et al. (2018) as explicitly described as free in the text and/or departing from their default value. There are many more SAMs in the literature that we could have added to Table 2, in principle. However, we found some papers to be opaque in their description of how their model is parametrised, and we therefore felt that we could not make a fair comparison with them. We expect those models fall in the same range as the others in Table 2.

For a deeper dive into the complexity of calibrating (other) SAMs, we refer the reader to Knebe et al. (2018).

14. Closing remarks

We have presented an overview of the salient features of the brand-new version of the Dark Sage SAM of galaxy formation. Dark Sage bucks the trend of most SAMs by numerically evolving the multidimensional structure of galaxy discs and by only having three free parameters that require calibration against observations.

The predictions from Dark Sage presented in this paper are but the tip of the iceberg as far as what the model can produce. The code codebase is open to the community, which we hope will serve as a useful resource for research into galaxy formation and cosmology.