# Quenching massive galaxies across cosmic time with the semi-analytic model
Shark v2.0

Claudia del P. Lagos1,2,3, Matías Bravo4, Rodrigo Tobar1, Danail
Obreschkow1,2, Chris Power1,2, Aaron S. G. Robotham1,2, Katy L. Proctor1,2,
Samuel Hansen1, Ángel Chandro-Gómez1,2, Julian Carrivick1
1International Centre for Radio Astronomy Research (ICRAR), M468, University
of Western Australia, 35 Stirling Hwy, Crawley,
WA 6009, Australia.
2ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D).
3Cosmic Dawn Center (DAWN), Denmark.
4Department of Physics & Astronomy, McMaster University, 1280 Main Street W,
Hamilton, ON, L8S 4M1, Canada E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We introduce version 2.0 of the Shark semi-analytic model of galaxy formation
after many improvements to the physics included. The most significant being:
(i) a model describing the exchange of angular momentum (AM) between the
interstellar medium and stars; (ii) a new active galactic nuclei feedback
model which has two modes, a quasar and a radio mode, with the radio mode tied
to the jet energy production; (iii) a model tracking the development of black
hole (BH) spins; (iv) more sophisticated modelling of environmental effects on
satellite galaxies; and (v) automatic parameter exploration using Particle
Swarm Optimisation. We focus on two timely research topics: the structural
properties of galaxies and the quenching of massive galaxies. For the former,
Shark v2.0 is capable of producing a more realistic stellar size-mass relation
with a plateau marking the transition from disk- to bulge-dominated galaxies,
and scaling relations between specific AM and mass that agree well with
observations. For the quenching of massive galaxies, Shark v2.0 produces
massive galaxies that are more quenched than the previous version, reproducing
well the observed relations between star formation rate (SFR) and stellar
mass, and specific SFR and BH mass at $z=0$. Shark v2.0 produces a number
density of massive-quiescent galaxies $>1$ dex higher than the previous
version, in good agreement with JWST observations at $z\leq 5$; predicts a
stellar mass function of passive galaxies in reasonably good agreement with
observations at $0.5<z<5$; and environmental quenching to already be effective
at $z=5$.

###### keywords:

galaxies: formation - galaxies: evolution

††pubyear: 2023††pagerange: Quenching massive galaxies across cosmic time with
the semi-analytic model Shark v2.0–D

##  1 Introduction

Our current theory of galaxy formation and evolution is intimately linked to
the growth of structures in the universe, which are thought to form
hierarchically. The current preferred cosmological model is the $\Lambda$ cold
dark matter ($\Lambda$CDM), in which the cosmic web evolves, for the most
part, by the effect of gravity. Cosmological simulations of galaxy formation
attempt to follow the formation of galaxies as the cosmic web forms from the
very early to the local universe, providing a wide range of predictions across
cosmic time, environment and galaxy populations (see reviews of Somerville et
al. 2015; Vogelsberger et al. 2020).

Among the most popular tools to simulate the formation of galaxies in a
cosmological context are hydrodynamical simulations and semi-analytic models
(SAMs). Both have pros and cons. Hydrodynamical simulations have the main
advantage of solving for the evolution of baryons and DM simultaneously,
avoiding some key simplifications made in SAMs regarding the symmetry of
galaxies and halos and the relevant baryonic components of each. The most
important advantages of SAMs are their speed, ability to thoroughly explore
the parameter space to understand potential degeneracies between physical
models and parameters, and as a result, the very large cosmological boxes they
can be run on, still covering a very large dynamic range in both halo and
stellar mass of the produced galaxies (Baugh, 2006; Benson & Bower, 2010).
Compared to the current generation of cosmological hydrodynamical simulations,
SAMs can push to about two to three orders of magnitude below the stellar mass
resolution of hydrodynamical simulations for the same cosmological volume. For
this reason, SAMs remain an essential part of the toolkit in the quest of
understanding galaxy formation and evolution.

In Lagos et al. (2018, hereafter L18), we introduced the Shark (v1.1) SAM, an
open source, flexible and highly modular SAM. Shark has been extensively used
for many applications (some of which are listed in $\S$ 2). One characteristic
of both SAMs and hydrodynamical simulations is the continuous development of
the physical models included to either introduce new physics into the models,
make more physical assumptions, and/or improve the agreement with observations
once tensions have been identified. Shark is no exception and since L18, the
model has seen continuous improvements on many of the physical models
included, where we recognised too simplistic assumptions were made, and to
address areas of tension that have been identified between the model and
observations. In this paper, we introduce a new version of Shark (v2.0), after
significant development of the physical models, which include, but are not
limited to, more physical models tracking the properties of supermassive black
holes (BHs) and active galactic nuclei (AGN) feedback, the angular momentum
evolution of galaxy components, and environmental effects that affect the
evolution of satellite galaxies. In this paper, we focus on two areas of
tension that have been identified between Shark v1.1 and observations: a
stellar mass-size relation that approximates a single power-law, while
observations display a clear plateau in the size-mass relation associated with
the transition from disk- to bulge-dominated galaxies (Lange et al., 2015);
and an overall scarcity of massive-quiescent galaxies at $z\gtrsim 2$ compared
with observations. More details about these two areas of tension are presented
in $\S$ 2.

The second tension above is particularly interesting as it has been reported
across many cosmological simulations of galaxy formation (Gould et al., 2023;
Valentino et al., 2023). This problem has worsened thanks to the exquisite
observations of the James Webb Space Telescope (JWST). These have revealed
that massive-quiescent galaxies are relatively common at $z>3$ (e.g. Carnall
et al. 2023; Valentino et al. 2023; Nanayakkara et al. 2022; Long et al. 2023)
and more so than previous observational inferences had indicated (Carnall et
al., 2020; Weaver et al., 2022; Gould et al., 2023). These galaxies typically
have stellar masses in excess of $10^{10}\,\rm M_{\odot}$ and number densities
$\gtrsim 10^{-5}\,\rm Mpc^{-3}$. Some of these galaxies have signs of having
ceased their star formation (a.k.a. quenched) recently ($\approx 100$ Myr),
while others are consistent with older stellar population ages ($\sim 1$ Gyr;
Glazebrook et al. 2023). The latter may imply these massive galaxies to have
formed at $z>10$, potentially posing a problem to structure formation in
$\Lambda$CDM (Boylan-Kolchin, 2023).

Although there are still many effects that could lead to systematic errors in
the inferred stellar masses, redshifts and star formation rates (SFRs) of
massive-quiescent galaxies in the observations, we focus here on how the
predictions around this population changes from Shark v1.1 and v2.0 after the
significant revisions of the model. We generalise the problem to understanding
the quenching of massive galaxies across cosmic times. We leave for future
work the critical assessment of systematic biases in the inferred properties
of massive-quiescent galaxies.

This paper is organised as follows. $\S$ 2 briefly describes the Shark model
and the $N$-body DM only simulations we use. $\S$ 3 describes in detail the
significant modifications made to Shark in this new version 2.0 relative to
v1.1, and the parameters we adopt for the default Shark v2.0 model. Minor
modifications are presented in Appendix A. $\S$ 4 presents key results on the
abundance of galaxies across cosmic time and scaling relations, including
structural relations in $\S$ 4.3.2, in the local universe. We also present a
supplementary material with additional comparisons with observations to show
that previous areas of agreement between Shark and observations remain so.
$\S$ 5 focuses on the quenching of galaxies in Shark v2.0 and compares with
v1.1 to reveal how much more efficient quenching is across cosmic time. The
reader interested solely in the problem of massive-quiescent galaxies at
$z\gtrsim 2$ can skip to $\S$ 5.2, where we cover this in detail. $\S$ 6
presents our main conclusions.

##  2 The semi-analytic model Shark

Shark, hosted on GitHub111https://github.com/ICRAR/shark, takes into account
physical processes that we think are critical in shaping the formation and
evolution of galaxies. These are (i) the collapse and merging of DM halos;
(ii) the accretion of gas onto halos, which is modulated by the DM accretion
rate; (iii) the shock heating and radiative cooling of gas inside DM halos,
leading to the formation of galactic disks via conservation of specific
angular momentum of the cooling gas; (iv) star formation (SF) in galaxy disks;
(v) stellar feedback from the evolving stellar populations; (vi) chemical
enrichment of stars and gas; (vii) the growth via gas accretion and merging of
BHs; (viii) heating by AGN; (ix) photoionization of the intergalactic medium;
(x) galaxy mergers driven by dynamical friction within common DM halos which
can trigger starbursts (SBs) and the formation and/or growth of spheroids; (x)
collapse of globally unstable disks that also lead to SBs and the formation
and/or growth of bulges. L18 included several different models for gas
cooling, AGN feedback, stellar and photo-ionisation feedback, and star
formation.

The model presented in L18, has been tested thoroughly across several
publications and shown to reproduce several observed relations, including: the
optical colour bimodality and how this depends on stellar mass (Bravo et al.,
2020); the atomic hydrogen (HI)-halo mass relation and HI clustering (Chauhan
et al., 2020, 2021); the panchromatic emission of galaxies from the far-
ultraviolet (FUV) to the far-infrared (FIR) across cosmic time (Lagos et al.,
2019, 2020; Chen et al., 2023); the redshift distribution of bright FIR
galaxies and the redshift evolution of their number density (Lagos et al.,
2019; Casey et al., 2021; Long et al., 2022), among many other successes. In
addition, Shark has been used to make predictions for the gravitational wave
signal from binary stellar (Rauf et al., 2023) and supermassive (Curyło &
Bulik, 2022) black holes.

Some areas of tension with observations have also been identified, such as a
weak downsizing signal and quenching timescales that are independent of
stellar mass (Bravo et al., 2023); a number density of passive galaxies that
appears to be too low at $z\geq 3$ (Long et al., 2022; Gould et al., 2023); a
stellar size-mass relation that is very close to a single power-law and
massive-end of the stellar mass function (SMF) that is too shallow compared to
observations (Proctor et al. in preparation). Proctor et al. (in preparation)
crucially show that some of these limitations were not easily solved by
modifying the parameters of the model and instead were inherent to the
physical models included in Shark.

The identification of these tensions in addition to the desire to continue to
improve the model have led us to the continuing improvement and inclusion of
new models within Shark. These new models are described in detail in $\S$ 3.
We have solidified these changes into a new version of Shark v2.0 available on
GitHub. Before describing the new models, we introduce the suite of $N$-body
simulations on top of which Shark runs in $\S$ 2.1.

###  2.1 The surfs simulations: halos catalogues and merger trees

Table 1: Simulation parameters of the surfs run used in this paper. Name | L210N1536
---|---
Box size [$\rm cMpc/h$] | 210
Number of particles | $1536^{3}$
Particle Mass [$\rm M_{\odot}/h$] | $2.21\times 10^{8}$
Softening Length [$\rm ckpc/h$] | 4.5

We run Shark over the surfs suite of $N$-body, DM only simulations (Elahi et
al., 2018). Most of the surfs runs have cubic volumes of $210\,\rm cMpc/h$ on
a side, and span a range in particle number, currently up to $8.5$ billion
particles, and adopt a $\Lambda$CDM Planck Collaboration et al. (2016)
cosmology. The cosmological parameters correspond to a total matter, baryon
and $\Lambda$ densities of $\Omega_{\rm m}=0.3121$, $\Omega_{\rm b}=0.0491$
and $\Omega_{\rm L}=0.6879$, respectively, with a Hubble parameter of $H_{\rm
0}=h\,100\,\rm Mpc\,km\,s^{-1}$ with $h=0.6751$, scalar spectral index of
$n_{\rm s}=0.9653$ and a power spectrum normalization of $\sigma_{\rm
8}=0.8150$. surfs was produced using a memory lean version of the gadget2 code
on the Magnus supercomputer at the Pawsey Supercomputing Centre. In this
paper, we use the L210N1536 simulation, with the specifications of Table 1.
surfs produces $200$ snapshots for each simulation, typically having a time
span between snapshots in the range of $\approx 6-80$ Myr.

Merger trees and halo catalogs were constructed using the phase-space finder
VELOCIraptor (Elahi et al., 2019a; Cañas et al., 2019) and the halo merger
tree code TreeFrog, developed to work on VELOCIraptor (Elahi et al., 2019b).
We refer to L18 for more details on how the merger trees and halo catalogs are
constructed for Shark, and to Poulton et al. (2018a); Elahi et al. (2019a, b);
Cañas et al. (2019) for details of VELOCIraptor and TreeFrog.

From the VELOCIraptor catalogues, we take the halo and subhalo masses, and
halo virial radii. Other properties of halos are calculated as described in
$\S$ 4.2 in L18. The existence of central and satellite subhalos in the
VELOCIraptor catalogues leads to Shark galaxies existing in 3 different types:
type = 0 is the central galaxy of the central subhalo, while every other
central galaxy of satellite subhalos are type = 1. If a subhalo merges onto
another one and it is not the main progenitor, it is treated as defunct. All
the galaxies of defunct subhalos are made type = 2 and transferred to the list
of galaxies of the central subhalo of their descendant host halo (see $\S$ 4.1
in L18 for more details).

##  3 New baryon physics models and technical updates in Shark

The modifications presented in this paper have been released in the public
repository of Shark in the version v2.0. Below we introduce the physical
models that signify large improvements from our previous installment of Shark.
In addition, we have implemented small changes to other existing physical
models as described in Appendix A.

###  3.1 ISM-stars angular momentum exchange

In this new version of Shark we include different levels of complexity for the
exchange of specific angular momentum between the cooling gas, interstellar
medium (ISM) and stellar disk. Here, the specific angular momentum will be
referred to as $j\equiv J/M$. The default in L18 assumes the cooling gas to
have the same $j$ of the DM halo,

$j_{\rm cool}=\frac{J_{\rm h}}{M_{\rm halo}},$ (1)

where $J_{\rm h}$ is the halo’s angular momentum, which is calculated from the
mass and halo’s spin parameter, following Mo et al. (1998),

$J_{\rm h}=\frac{\sqrt{2}\,G^{2/3}}{(10\,H(z))^{1/3}}\lambda_{\rm
DM}\,M^{5/3}_{\rm halo}.$ (2)

$j_{\rm cool}$ is then input in the set of ordinary differential equations
(ODEs) that control the exchange of angular momentum (Eqs. (5)-(9)). Section
4.2 of L18 describes how $\lambda_{\rm DM}$ is obtained and Appendix A
presents a minor modification to the calculation of $\lambda_{\rm DM}$
introduced in Shark v2.0.

The gaseous and stellar disks also exchange angular momentum at a rate
$\dot{J}_{\rm g,s}$. In its simplest form,

$\dot{J}_{\rm g,s}=\psi\,j_{\rm cold},$ (3)

where $\psi$ is the instantaneous SFR and $j_{\rm cold}$ (not to be confused
with $j_{\rm cool}$) is the specific angular momentum of the gaseous disk.
This is the default assumption made in L18 and in most SAMs, except for those
that follow disks as a set of annuli that can evolve independently (e.g.
Stringer & Benson 2007; Stevens et al. 2016a).

Mitchell et al. (2018) showed that in the cosmological hydrodynamical
simulation eagle (Schaye et al., 2015; Crain et al., 2015; McAlpine et al.,
2015), the stellar $j$ was systematically lower than in the GALFORM SAM at
fixed stellar mass. Mitchell et al. (2018) suggested that a key modification
required by GALFORM to evolve the specific angular momentum of galaxies more
realistically was to consider the fact that stars form from molecular gas
only. The latter tends to be more concentrated in the centres of galaxies
relative to the total ISM (Lagos et al., 2011). This implies that stars will
be systematically forming from low specific angular momentum in comparison to
the total ISM. Mitchell et al. (2018) showed that this physical effect is
behind the overly large sizes of galaxies in GALFORM. Thus, a more
sophisticated model should include the fact that stars form from molecular
gas, which resides in the high column density regions of disks, and is more
concentrated than the total gas disk.

Here, we implement a more realistic exchange between the specific angular
momentum of the ISM and the stars following the lessons from Mitchell et al.
(2018). This is done by calculating $\dot{J}_{\rm g,s}$ as

$\dot{J}_{\rm g,s}=2\,\pi\,\int_{0}^{\infty}v(r)\,r\,\Sigma_{\rm SFR}\,{\rm
d}r,$ (4)

where $\Sigma_{\rm SFR}$ depends on the assumed SF law and $v(r)$ is the
circular velocity radial profile. We include a boolean parameter in Shark,
angular-momentum-transfer, which set to true means the stellar disk angular
momentum is calculated as Eq. (4), and to false uses Eq. (3) instead. Thus,
the results of L18 can be fully recovered with the latter option.

The half-mass gas and stellar disk sizes are then calculated as $r_{\rm
gas}=f_{\rm norm}\,j_{\rm cold}/V_{\rm circ}$ and $r_{\star}=f_{\rm
norm}\,j_{\star}/V_{\rm circ}$. Here, we set $f_{\rm norm}=0.839$, which is
the value in an idealized exponential disk (Guo et al., 2011).

Before presenting the way we solve for the evolution of the galaxy angular
momentum split by components, we look at the assumptions we make. From the
definition $J=j\,M$, with $J\equiv|\vec{J}|$ and $j\equiv|\vec{j}|$, it
follows that
${\rm\partial}J/{\rm\partial}t=M\,{\rm\partial}j/{\rm\partial}t+j\,{\rm\partial}M/{\rm\partial}t$.
The two terms represent angular momentum transfer by torques and advection,
respectively. We assume that internal to galaxies and from the gas cooling
down to galaxies, ${\rm\partial}j/{\rm\partial}t=0$ (i.e. no torques), and
hence ${\rm\partial}J/{\rm\partial}t=j\,{\rm\partial}M/{\rm\partial}t$. This
is a simplification, as cosmological hydrodynamical simulations have shown
torques can be important between the gas cooling and the galaxy (e.g. Stevens
et al. 2016b). The second important assumption is that all the components in
the galaxy and gas cooling have angular momentum vectors that are aligned.
This is not necessarily the case. Contreras et al. (2017) showed that halos
have an angular momentum vector that can evolve significantly in direction
under the presence of high mass mergers. Quieter assembly histories lead to
smaller changes of the direction of the vector. This in principle could lead
to misalignments between the halo and the galaxy, as shown by Lagos et al.
(2015). The two assumptions above are made to simplify the problem, but we
will consider in the future models that relax these assumptions and study the
effect that would have on the scaling relations analysed in $\S$ 4.3.2.

Simultaneously to the mass and metal exchange (see Eqs. ($49$) to ($58$) in
L18), in the case of SF in disks, we solve for the angular momentum exchange
between these components following the assumptions above, as follows:

$\displaystyle\dot{J}_{\star}$ $\displaystyle=$
$\displaystyle(1-R)\,\dot{J}_{\rm g,s}$ (5) $\displaystyle\dot{J}_{\rm cold}$
$\displaystyle=$ $\displaystyle\dot{M}_{\rm cool}\,j_{\rm
cool}-(1-R+\beta_{\star})\,\dot{J}_{\rm g,s}$ (6) $\displaystyle\dot{J}_{\rm
cold,halo}$ $\displaystyle=$ $\displaystyle-\dot{M}_{\rm cool}\,j_{\rm cool}$
(7) $\displaystyle\dot{J}_{\rm hot,halo}$ $\displaystyle=$
$\displaystyle\dot{m}_{\rm out,\star}\,j_{\rm out}-\dot{m}_{\rm ejec}\,j_{\rm
out}$ (8) $\displaystyle\dot{J}_{\rm ejec}$ $\displaystyle=$
$\displaystyle\dot{m}_{\rm ejec}\,j_{\rm out}.$ (9)

Here, $\dot{J}_{\rm g,s}$ is as described in Eq. (4), $\dot{M}_{\rm cool}$,
$\dot{m}_{\rm out,\star}$ and $\dot{m}_{\rm ejec}$ are the gas cooling rate
(Eq. (5) in L18), outflow rate due to star formation and ejection rate from
the halo, respectively; $\beta_{\star}=\dot{m}_{\rm out,\star}\,\psi^{-1}$,
where $\psi$ is the instantaneous SFR, is the wind mass loading, and $R$ is
the fraction of mass recycled to the ISM (from stellar winds and supernovae;
see § 4.4.6 in L18 for details). In the case of the hot halo and ejected gas
mass components, the angular momentum growth depends on the specific angular
momentum of the outflowing gas. This in principle allows for outflows to
affect the angular momentum of the disk in a differential fashion, which would
be the case if the outflow rate was an explicit function of radius (as it has
been proposed by detailed stellar feedback models, e.g. Creasey et al. 2013;
Hopkins et al. 2012; Lagos et al. 2013). However, in Shark v2.0 we assume
$j_{\rm out}=j_{\rm cold}$, and leave for future work the exploration of how
outflow rates that change radially can impact the angular momentum structure
of galaxies.

Note that implicitly in Eqs. (5)-(9) there is a third assumption, which is
that the relaxation time of the hot halo gas is such that after a SF episode,
it regains the same specific angular momentum as the DM halo. This assumes
that any change in the hot halo gas $j$ due to outflows is negligible.

###  3.2 Dynamical friction leading to galaxy mergers

Inherent limitations of halo and subhalo finders generally imply that subhalos
cease to be tracked when they reach too few particles, or when they become
indistinguishable from the 3D or 6D density background, depending on whether a
3D or 6D structure finder is used. This does not imply that galaxies in those
subhalos should be immediately merge onto the central galaxy, as in some cases
subhalos stop being tracked at large distances from the halo’s centre
(typically 0.5–1 times the virial radius; see Poulton et al., 2018b). For
SAMs, this means that we have to resort to analytic calculations of the
dynamical friction timescale to have a more realistic proxy for when a
satellite galaxy should merge onto the central.

Recent work by Poulton et al. (2020) tracks (sub)halos that are identified in
dense environments in $N$-body simulations, and provides estimates for their
orbital properties. By comparing the merger timescale of these simulated
subhalos with various analytical prescriptions, they are able to identify
regimes where previously published analytical prescriptions do not provide a
realistic estimate of merger timescales. Specifically, they find that in
systems with relatively large host-to-subhalo mass ratios, these analytic
approximations all systematically underpredict the merger timescale. This has
important implications for both the satellite population and the high-mass
galaxy population reported in SAMs, given the importance of mergers in the
mass assembly of large galaxies that is expected from a hierarchical model of
structure formation (Robotham et al., 2014).

Follow-up work presented in Poulton et al. (2021) provided a new analytic
model to compute the dynamical friction timescale ($\tau_{\rm{merge}}$) of
satellite galaxies that much more accurately captures the merger timescales
seen over a wide dynamic range in halo masses. In this work we implement this
new dynamical friction timescale in Shark. The default model presented in L18
adopted the dynamical friction timescale of Lacey & Cole (1993), which
compared to Poulton et al. (2021) merges small galaxies too quickly.

Poulton et al. (2021) formulate the merger timescale in two regimes, based on
their result that $\tau_{\rm{merge}}$ is more strongly dependent on position
for subhalos outside of the virial radius of the host than for subhalos inside
the virial radius. The merger timescale is calculated as

$\displaystyle\tau_{\rm{merge}}=\begin{cases}5.62\sqrt{\frac{R_{\rm
vir,host}}{GM_{\rm encl,host}(r)}}r^{0.8}R_{\rm peri}^{0.2}&\text{for
}r<R_{\rm{vir,host}},\\\ 5.62\frac{R_{\rm vir,host}}{\sqrt{GM_{\rm
vir,host}}}r^{0.3}R_{\rm peri}^{0.2}&\text{for }r\geq R_{\rm{vir,host}},\\\
\end{cases}$ (10)

where $r$ is the position of the subhalo relative to the halo’s centre,
$R_{\rm peri}$ is the pericentric distance from the host centre, $M_{\rm
encl,host}(r)$ is the mass of the host enclosed within a radius $r$ assuming
an Navarro et al. (1997) profile, and $M_{\rm vir,host}$ is the virial mass of
the host halo. These properties are all determined from the underlying
$N-$body simulation. In Shark this model can be adopted by setting
merger_model = poulton20. Note that the default in L18 was merger_model =
lacey93.

###  3.3 AGN feedback: QSO and radio mode feedback

The default model presented in L18 adopted the Croton et al. (2016) AGN
feedback model. That model in itself was inspired by Croton et al. (2006),
which included only a radio mode of AGN feedback. The heating power of AGN
feedback in the radio mode of Croton et al. (2006) was calculated using the
bolometric luminosity of AGN. For the latter, only the BH accretion rate
coming from the hot-halo mode was considered (Eq. $10$ in Croton et al. 2006).
However, BHs at any one time can be growing by multiple channels; in
particular, the accretion rate onto the BHs can be much higher than that
coming from the hot-halo mode only if the galaxy is experiencing a SB (which
in Shark can be either triggered by galaxy mergers or disk instabilities).
There is no physical reason for why in these cases only the fractional
contribution to the accretion rate from the hot-halo mode should be included
in the AGN heating power calculation. In Shark v2.0, we implement a new model
for AGN feedback that includes two modes: a radio and a QSO mode, and uses the
physical properties of the BH to calculate a jet power and a QSO-driven
outflow rate and velocity. Before describing the new AGN feedback models, we
describe how the three fundamental properties of BHs, their mass, accretion
rate and spin, are calculated.

####  3.3.1 BH mass, accretion rate and spin

In addition to BH-BH mergers, in which we instantaneously add the masses of
the BHs that are merging, BHs can also grow via gas accretion in two different
modes. The BH accretion rate due to hot-halo mode is calculated as in L18,

$\dot{m}_{\rm BH,hh}=\kappa\,\frac{15}{16}\pi\,G\,\mu\,m_{\rm
p}\,\frac{\kappa_{\rm B}\,T_{\rm vir}}{\Lambda(T_{\rm vir},Z_{\rm
hot})}\,m_{\rm BH},$ (11)

where $m_{\rm BH}$ is the BH mass, $G$ is the gravitational constant, $\mu$ is
the atomic weight, $m_{\rm P}$ the proton mass, $\kappa_{\rm B}$ is
Boltzmann’s constant, $T_{\rm vir}$ the virial temperature of the halo,
$Z_{\rm hot}$ the halo gas metallicity, $\Lambda(T_{\rm vir},Z_{\rm hot})$ the
cooling function, and $\kappa$ a free parameter. For the case of BH growth
during SBs, we follow L18,

$\delta\,m_{\rm BH,sb}=f_{\rm smbh}\,\frac{m_{\rm gas}}{1+(v_{\rm smbh}/V_{\rm
vir})^{2}},$ (12)

where $m_{\rm gas}$ and $V_{\rm vir}$ are the cold gas mass reservoir of the
SB and the virial velocity of the halo, respectively. $f_{\rm smbh}$ and
$v_{\rm smbh}$ are free parameters. As in L18, we estimate the BH accretion
rate in this mode assuming that the accretion timescale scales with the bulge
dynamical timescale, $\tau_{\rm acc,sb}=e_{\rm sb}\,r_{\rm bulge}/v_{\rm
bulge}$, where $e_{\rm sb}$ is an e-folding parameter. The accretion rate
during SBs is thus,

$\dot{m}_{\rm BH,sb}=\frac{\delta\,m_{\rm BH,sb}}{\tau_{\rm acc,sb}}.$ (13)

The total accretion rate onto the BH at any one time is $\dot{m}_{\rm
BH}=\dot{m}_{\rm BH,hh}+\dot{m}_{\rm BH,sb}$. The BH accretion disk structure
is expected to be a strong function of the accretion rate. We thus define a
normalised accretion rate based on the Eddington luminosity, $L_{\rm Edd}$,

$\dot{m}=\frac{\dot{m}_{\rm BH}}{\dot{m}_{\rm Edd}},$ (14)

where $\dot{m}_{\rm Edd}=L_{\rm Edd}/(0.1\,c^{2})$, and $c$ is the speed of
light.

When $\dot{m}>\dot{m}_{\rm ADAF}$, the accretion disk is expected to be thin
and efficiently cool; this regime is commonly referred to as thin-disk (TD;
Shakura & Sunyaev 1973). If instead, $\dot{m}<\dot{m}_{\rm ADAF}$, the
accretion disk is expected to be unable to cool efficiently by radiation due
to the energy generated by viscosity, commonly referred to as the Advection
Dominated Accretion Flow (ADAF; Rees 1982) regime. It is broadly assumed the
transition between TD and ADAF happens at $\dot{m}_{\rm ADAF}=0.01$. The ADAF
regime according to Mahadevan (1997), can be further subdivided into two
regimes: the lower accretion rate ADAF regime ($\dot{m}<\dot{m}_{\rm
crit,visc}$), in which heating of the electrons is dominated by viscous
heating, and a higher accretion rate ADAF regime ($\dot{m}_{\rm
crit,visc}<\dot{m}<\dot{m}_{\rm ADAF}$), in which the ion-electron heating
dominates the heating of the electrons. On the other extreme, we classify AGN
as super-Eddington (SE) if $\dot{m}>\eta$, with $\eta\sim 1$. In summary, four
BH accretion regimes are defined:

  * •

SE: $\dot{m}>\eta$,

  * •

TD: $\dot{m}\geqslant\dot{m}_{\rm ADAF}$ and $\dot{m}\leq\eta$,

  * •

ADAFhigh: $\dot{m}_{\rm crit,visc}<\dot{m}<\dot{m}_{\rm ADAF}$,

  * •

ADAFlow: $\dot{m}\leq\dot{m}_{\rm crit,visc}$.

The dimensionless BH spin vector, ${\bf a}$, is defined as ${\bf a}\equiv{\bf
J}_{\rm BH}/J_{\rm max}=c\,{\bf J}_{\rm BH}/G\,m^{2}_{\rm BH}$, where ${\bf
J}_{\rm BH}$ is the angular momentum vector of the BH. To calculate ${\bf a}$,
we implement three different models, which in Shark v2.0 can be selected by
setting the variable spin-model. Below, we refer to the norm of the spin
vector as $a=|{\bf a}|$. The list of models are presented below:

  * •

constant. This corresponds to the simplest assumption of a constant spin. The
radiation efficiency of a BH depends on the radius of the last stable orbit,
which in itself depends on the BH spin. The value adopted for the constant
spin in Shark v2.0 is $0.67$, which is equivalent to assuming a constant
radiation efficiency of $0.1$ (Bardeen et al., 1972). (Note that $a=0$ gives
$\epsilon=0.057$ and $a=1$ gives $\epsilon=0.42$; see Eqs. (16)-(19)).

  * •

volonteri07. This is a simple scaling relation inspired by the spin-BH mass
relation found by Volonteri et al. (2007), who presented a model of BH
accretion from a warped disk. The authors found that on average more massive
BHs have a higher spin. We use the average relation they found and scale the
spin directly with the BH mass as

$a=0.305\,{\rm log_{10}}\left(\frac{m_{\rm BH}}{\rm M_{\odot}}\right)-1.7475.$
(15)

We limit $a$ to be in the range $[0,1]$ in this model.

  * •

griffin19. This model is the full implementation in Shark v2.0 of the
sophisticated BH spin development model of Griffin et al. (2019). They
presented a model that follows the changes in BH spin produced by BH-BH merger
and gas accretion. For the latter, Griffin et al. (2019) presented three
different models, which we also implemented in Shark. These are the prolonged
accretion, self-gravitating accretion disk and warped accretion disk. The user
can choose between these three models by setting the variable accretion-disk-
model. Details of these models are presented in Appendix B. Note that in this
model $a$ can take values in the range $[-1,1]$. This model includes the full
implementation of the warped accretion disk model described in Volonteri et
al. (2007), rather than a fit to their results as adopted by spin-
model=volonteri07.

With the BH accretion rate and spin defined, we can calculate the AGN
bolometric luminosity. We first define the radiative accretion efficiency for
a thin accretion disk, $\epsilon_{\rm TD}$,

$\epsilon_{\rm TD}=1-\sqrt{1-\frac{2}{3\,\hat{r}_{\rm lso}}},$ (16)

where $\hat{r}_{\rm lso}$ is the radius of the last stable circular orbit in
units of the gravitational radius of the BH, $r_{\rm G}\equiv G\,m_{\rm
BH}/c^{2}$. We calculate $\hat{r}_{\rm lso}$ following Bardeen et al. (1972)

$\hat{r}_{\rm lso}=3+Z_{2}\pm\sqrt{(3-Z_{1})(3+Z_{1}+2\,Z_{2})},$ (17)

with the negative (positive) sign corresponding to the case when the angle
between the BH spin and the accretion disk is less (larger) than $90$ degrees.
We refer to these two cases as co- and counter-rotation, respectively. $Z_{1}$
and $Z_{2}$ are define as,

$Z_{1}=1+(1-a^{2})^{1/3}\left[(1+a)^{1/3}+(1-a)^{1/3}\right],$ (18)

and

$Z_{2}=\sqrt{3a^{2}+Z^{2}_{1}}.$ (19)

Note that BHs can be counter-rotating only in the griffin19 spin model.

We compute the AGN bolometric luminosity following Griffin et al. (2019) for
the thin disk and ADAF regimes, and following Griffin et al. (2020) for the
super-Eddington regime:

$L_{\rm bol}=\begin{cases}\eta\,\left(1+{\rm
ln}\left(\frac{\dot{m}}{\eta}\,\frac{\epsilon_{\rm
TD}}{0.1}\right)\right)\,L_{\rm Edd},\,\text{in SE}\\\ \epsilon_{\rm
TD}\dot{m}_{\rm BH}\,c^{2},\,\text{in TD}\\\ 0.2\,\epsilon_{\rm
TD}\dot{m}_{\rm BH}\,c^{2}\left(\frac{\dot{m}}{\alpha^{2}_{\rm
ADAF}}\right)\,\left(\frac{\beta}{0.5}\right)\,\left(\frac{6}{\hat{r}_{\rm
lso}}\right),\,\text{in ADAF${}_{\mathrm{high}}$}\\\ 0.0002\,\epsilon_{\rm
TD}\dot{m}_{\rm BH}\,c^{2}\left(\frac{\delta_{\rm
ADAF}}{0.0005}\right)\,\left(\frac{1-\beta}{0.5}\right)\,\left(\frac{6}{\hat{r}_{\rm
lso}}\right),\,\text{in ADAF${}_{\mathrm{low}}$}\end{cases}$ (20)

Here, $\alpha_{\rm ADAF}$ is the viscosity parameter in the ADAF regime,
$\delta_{\rm ADAF}$ is the fraction of viscous energy transferred to the
electrons, which has a value between $0.1$ and $0.5$ (see review of Yuan &
Narayan 2014), and $\beta$ is the ratio of gas pressure to total pressure
(i.e. the sum of gas pressure and magnetic pressure). Following Griffin et al.
(2019), $\beta=1-\alpha_{\rm ADAF}/0.55$, $\alpha_{\rm ADAF}=0.1$ and
$\delta_{\rm ADAF}=0.2$. With these parameters defined, we introduce the
boundary between the two ADAF regimes (Griffin et al., 2019),

$\dot{m}_{\rm crit,visc}=0.001\,\left(\frac{\delta_{\rm
ADAF}}{0.0005}\right)\,\left(\frac{1-\beta}{\beta}\right)\,\alpha^{2}_{\rm
ADAF}.$ (21)

####  3.3.2 AGN radio-mode feedback

In Shark v2.0, the radio mode of AGN feedback is assumed to only occur when
halos have reached a quasi-hydrostatic equilibrium (i.e. halos have developed
a hot gaseous halo). To evaluate this, we use the Correa et al. (2018)
criterion, in which a cooling and heating terms are calculated and compared.
These terms are computed per halo as follows,

$\displaystyle\Gamma_{\rm cool}(M_{\rm halo},r)$ $\displaystyle=$
$\displaystyle M_{\rm hot}\,\frac{n_{\rm H}(r)\,\Lambda(T_{\rm vir},Z_{\rm
hot})}{\mu\,m_{\rm p}},$ (22) $\displaystyle\Gamma_{\rm heat}(M_{\rm halo})$
$\displaystyle=$ $\displaystyle\frac{3\,k_{\rm B}\,T_{\rm vir}}{2\,\mu\,m_{\rm
p}}\,\frac{\Omega_{\rm b}}{\Omega_{\rm m}}\dot{M}_{\rm
halo}\,\left(\frac{2}{3}\,f_{\rm hot}+f^{\rm halo}_{\rm
acc,hot}\right).\,\,\,\,\,\,\,\,\,$ (23)

Here, $M_{\rm halo}$ is the virial mass, $M_{\rm hot}$ is the halo gas mass,
$n_{\rm H}(r)$ is the gas volume density of hydrogen atoms at $r$,
$\Omega_{\rm b}$ and $\Omega_{\rm m}$ are cosmological parameters,
$\dot{M}_{\rm halo}$ is the matter accretion rate onto a halo, $f_{\rm hot}$
is the fraction of gas in the halo relative to the universal baryon fraction
($\equiv M_{\rm hot}/(\Omega_{\rm b}/\Omega_{\rm m}\,M_{\rm halo})$) and
$f^{\rm halo}_{\rm acc,hot}$ is the fraction of the accreting gas that is hot.
The latter fractions were parametrised as a function of the halo mass by
Correa et al. (2018) as follows,

$\displaystyle f_{\rm hot}$ $\displaystyle=$ $\displaystyle
10^{-0.8+0.5\,x-0.05\,x^{2}},$ (24) $\displaystyle f^{\rm halo}_{\rm acc,hot}$
$\displaystyle=$ $\displaystyle\frac{1}{e^{-4.3\,(x+0.15)}+1},$ (25)

with $x$ defined as,

$x={\rm log}_{10}\left(\frac{M_{\rm halo}}{10^{12}\,\rm M_{\odot}}\right).$
(27)

The matter accretion rate onto halos is calculated using Dekel et al. (2009),

$\frac{\dot{M}_{\rm halo}}{\rm M_{\odot}\,Gyr^{-1}}=0.47\,\left(\frac{M_{\rm
halo}}{10^{12}\,\rm
M_{\odot}}\right)^{0.15}\,\left(\frac{1+z}{3}\right)^{2.25}\,\frac{M_{\rm
halo}}{\rm M_{\odot}}.$ (28)

In principle, we can measure $\dot{M}_{\rm halo}$ directly from the
VELOCIraptor and TreeFrog catalogues. However, there are cases of mass
swapping between merger branches which can lead to sudden large changes in the
halo mass. To avoid such big discontinuities in mass, we opt to use the
fitting function above. Chandro-Gómez et al. (in preparation) present a
detailed analysis of the frequency of mass-swapping events for different
simulations, subhalo finders and tree builders, and the effect they can have
on galaxies in Shark.

A halo is considered to be capable of forming a hot halo when $\Gamma_{\rm
heat}>\Gamma_{\rm cool}$. Under this condition, the accumulated shock-heated
gas at the virial radius gains the necessary pressure through external shock-
heating to overcome the energy loss from radiative cooling. In Shark, we
assume the relevant density in Eq. (22) to be the gas density of the halo at
$R_{\rm vir}$ where the shocks for the accreting gas are expected to happen,
which we approximate as $n_{\rm H}(R_{\rm vir})=200\,\rho_{\rm
crit}/\mu\,m_{\rm P}$, where $\rho_{\rm crit}$ is the critical density of the
universe. To allow flexibility in Shark, we include the parameter $\Gamma_{\rm
thresh}$, so that halos with $\Gamma_{\rm cool}/\Gamma_{\rm hot}<\Gamma_{\rm
thresh}$ are considered to have formed a hot halo. Appendix C shows at which
halo mass, most of the halos comply with $\Gamma_{\rm cool}/\Gamma_{\rm
hot}<\Gamma_{\rm thresh}$.

In the default model presented in L18, the heating rate of the AGN in the
radio mode was computed from the AGN bolometric luminosity produced by the
hot-halo mode accretion rate only (Croton et al., 2006),

$\displaystyle\dot{m}_{\rm heat}$ $\displaystyle=$ $\displaystyle\frac{L_{\rm
hh}}{0.5\,V^{2}_{\rm vir}},$ (29) $\displaystyle L_{\rm hh}$ $\displaystyle=$
$\displaystyle 0.1\,\dot{m}_{\rm BH,hh}\,c^{2}.$ (30)

In the new AGN feedback model, we compute the jet power (summed over both
jets, assuming jets to be symmetrical) following Meier (2002),

$\frac{Q_{\mathrm{mech}}}{\rm erg\,s^{-1}}=\begin{cases}2.5\cdot
10^{43}\left(\frac{M_{\mathrm{BH}}}{10^{9}M_{\odot}}\right)^{1.1}\left(\frac{\dot{m}}{0.01}\right)^{1.2}a^{2}\text{
if $\dot{m}\geq\dot{m}_{\rm ADAF}$,}\\\ 2\cdot
10^{45}\left(\frac{M_{\mathrm{BH}}}{10^{9}M_{\odot}}\right)\left(\frac{\dot{m}}{0.01}\right)a^{2}\text{
if $\dot{m}<\dot{m}_{\rm ADAF}$.}\end{cases}$ (31)

In our new model, a fraction $\kappa_{\rm radio}$ of $Q_{\rm mech}$ is used to
offset $\dot{M}_{\rm cool}$, so that the new cooling luminosity is reduced by
$\kappa_{\rm radio}\,Q_{\rm mech}$. If $L_{\rm cool}<\kappa_{\rm
radio}\,Q_{\rm mech}$ the cooling flow is completely shut off. Otherwise, the
cooling rate is defined as

$\dot{M}^{\prime}_{\rm cool}=\dot{M}_{\rm cool}\,\left(1-\frac{\kappa_{\rm
radio}\,Q_{\rm mech}}{L_{\rm cool}}\right).$ (32)

Here, $\kappa_{\rm radio}$ is a free parameter controlling how efficient the
jet power is in heating the halo gas. Note that in principle, $\kappa_{\rm
radio}$ can be $>1$ if the jets are efficient in producing buoyant bubbles
that lift up gas, as recently suggested by Huško & Lacey (2023).

Note that in this model, both BH accretion modes, hot-halo and SB, contribute
to the mechanical power of jets, and hence radio-mode feedback can happen
regardless of the source of gas accretion.

In Shark v2.0, satellite galaxies can continue to accrete gas that is cooling
from the hot halo they have retained (as it is gradually stripped; see $\S$
3.4). This means that satellite galaxies can also undergo radio mode AGN
feedback. The way this operates is the same as above. However, in situations
where $L_{\rm cool}<\kappa_{\rm radio}\,Q_{\rm mech}$ in a satellite galaxy,
we assume that the excess power, $P_{\rm excess}=\kappa_{\rm radio}\,Q_{\rm
mech}-L_{\rm cool}$, can be used to heat up the hot gas reservoir of the
central subhalo. If this is the case, then the effective heating power
affecting the gas cooling onto the central galaxy is larger than $\kappa_{\rm
radio}\,Q_{\rm mech}$ of the central galaxy by $\sum_{i}P_{\rm excess,\it i}$,
where $i$ are all the satellite galaxies with an excess power $>0$. In the
practice $P_{\rm excess}$ makes little difference to the heating rate
offsetting gas cooling of central subhalos.

####  3.3.3 AGN QSO-mode feedback

We consider outflows in QSOs to be driven by direct radiation pressure on dust
grains, as originally introduced by Fabian (1999). For this, we follow the
description of the production of outflows of Murray et al. (2005) and
complement that with the calculations of Ishibashi & Fabian (2015) for the
outflow velocities produced in the case of radiation pressure on dust grains.
In this model we use gas metallicity and gas content as proxies of dust mass
(as described below).

Murray et al. (2005) discussed the role of radiation pressure on dust in
driving galactic-scale winds through momentum deposition. Murray et al. (2005)
assumed an isothermal sphere profile, and argue that in the optically-thick
limit and from the momentum equation of the gas, one can derive a critical
luminosity at which the effective gravity is reduced by the momentum
deposition of the radiation. This critical luminosity for an isothermal sphere
can be written as,

$\frac{L_{\mathrm{crit}}}{3\cdot 10^{46}\ \mathrm{erg}\
\mathrm{s^{-1}}}\approx 10\,f_{\rm gas}\left(\frac{\sigma}{200\ \mathrm{km}\
\mathrm{s^{-1}}}\right)^{4},$ (33)

where $f_{\rm gas}$ is the gas fraction of the bulge and $\sigma$ its velocity
dispersion (which in this case we equate to the stellar velocity dispersion).
If the BH’s bolometric luminosity is in excess of this critical value,
$L_{\mathrm{bol}}>L_{\mathrm{crit}}$, the net motion of the gas in the bulge
would be outwards (i.e., outflow).

If an outflow is produced, and following Murray et al. (2005), it should be
capable of sweeping up the ISM gas that is outside the sublimation radius.
However, the latter tends to be very small ($<100$ pc) and hence, we consider
the total ISM gas mass in the bulge, $m_{\mathrm{gas,b}}$, to participate in
the outflow rate. We take as the relevant timescale for the outflow to be the
Salpeter time (i.e. the time to double the mass of the BH), defined as

$\tau_{\mathrm{Salp}}=43\Gamma^{-1}\ \mathrm{Myr}.$ (34)

where $\Gamma\equiv L_{\mathrm{bol}}/L_{\mathrm{Edd}}$. With this timescale,
we calculate the outflow rate to be

$\dot{m}_{\mathrm{out,QSO}}=\frac{m_{\mathrm{gas,b}}}{\tau_{\mathrm{Salp}}}.$
(35)

The reason why we consider the Salpeter time to be the relevant timescale here
is because this timescale assures us that the BH will grow to values that are
comparable to local BHs as measured in McConnell & Ma (2013), and because the
Salpeter timescale is of a similar magnitude as the duration of SBs. It also
ensures that enough stars are produced in the galaxies in Shark to lie in the
Faber-Jackson (Murray et al. 2005; Power et al. 2011) and BH-bulge mass
relations. Following the arguments presented in Nayakshin et al. (2009), using
$\tau_{\rm Salp}$ in Eq. (35) can potentially underestimate the impact of QSO
feedback in galaxies with small bulges (those with a stellar velocity
dispersion $<150\,\rm km\,s^{-1}$), but it is a fair representation of the
timescale at which BHs grow in more massive bulges. Because QSO feedback is
expected to be significant in massive galaxies only, we consider this
assumption to be reasonable.

We estimate the terminal velocity of the outflow following Ishibashi & Fabian
(2015) in the limit of radiation pressure on dust grains and assuming that the
whole ISM content is contained in an expanding shell:

$\displaystyle v_{\mathrm{out,QSO}}$ $\displaystyle\approx$ $\displaystyle
320\ \mathrm{km}\
\mathrm{s^{-1}}\left(\frac{L_{\mathrm{bol}}}{10^{7}L_{\odot}}\right)^{1/2}\left(\frac{\kappa_{\mathrm{UV}}}{10^{3}\
\mathrm{cm}^{2}\ \mathrm{g}^{-1}}\right)^{1/4}$ (36)
$\displaystyle\cdot\,\left(\frac{m_{\mathrm{gas,b}}}{M_{\odot}}\right)^{-1/4}.$

Here, $\kappa_{\mathrm{UV}}$ is the UV opacity. We use the Ishibashi & Fabian
(2015) approximation $\kappa_{\mathrm{UV}}=10^{3}\,f_{\mathrm{dg,MW}}\
\mathrm{cm}^{2}\ \mathrm{g^{-1}}$, with $f_{\mathrm{dg,MW}}$ being the dust-
to-gas mass ratio relative to the Milky-Way value. Assuming a constant
metallicity-dependent dust-to-gas mass ratio, we can then write
$\kappa_{\mathrm{UV}}=10^{3}\,(Z_{\mathrm{gas,bulge}}/Z_{\odot})\,\mathrm{cm}^{2}\
\mathrm{g}^{-1}$, where $Z_{\mathrm{gas,bulge}}$ is the fraction of mass in
metals relative to the total gas mass in the bulge.

Similar to the description of SB-driven outflows in L18, we define an excess
energy of the QSO-driven outflow that can be used to eject gas from the halo
as,

$E_{\mathrm{excess,QSO}}=\epsilon_{\mathrm{QSO}}\frac{v^{2}_{\mathrm{out,QSO}}}{2}\dot{m}_{\mathrm{out,QSO}}.$
(37)

Here, $\epsilon_{\mathrm{QSO}}$ is a free parameter in which we enclose
variations of geometry and other simplifications of our modelling. We then
define the gas ejection rate due to these QSO outflows as:

$\dot{m}_{\mathrm{ejec,QSO}}=\frac{E_{\mathrm{excess,QSO}}}{0.5\,V^{2}_{\mathrm{circ}}}-\dot{m}_{\mathrm{out,QSO}}$
(38)

where $V_{\mathrm{circ}}$ is the circular velocity of the halo. This can be
reduced to

$\dot{m}_{\mathrm{ejec,QSO}}=\left(\epsilon_{\mathrm{QSO}}\,\frac{v^{2}_{\mathrm{out,QSO}}}{v^{2}_{\mathrm{circ}}}-1\right)\,\dot{m}_{\mathrm{out,QSO}}.$
(39)

To include QSO feedback, we modify the differential equations (Eqs. $49-58$ in
L18) that control the evolution of the stellar ($M_{\star}$), cold gas
($M_{\mathrm{cold}}$), hot halo gas ($M_{\mathrm{hot}}$) and ejected gas
($M_{\mathrm{ejec}}$) masses, and their respective metals ($M^{Z}_{\star}$,
$M^{Z}_{\mathrm{cold}}$, $M^{Z}_{\mathrm{hot}}$ and $M^{Z}_{\mathrm{ejec}}$),
as follows:

$\displaystyle\dot{M}_{\star}$ $\displaystyle=$ $\displaystyle(1-R)\psi$ (40)
$\displaystyle\dot{M}_{\rm cold}$ $\displaystyle=$ $\displaystyle\dot{M}_{\rm
cool}-(1-R+\beta_{\star}+\beta_{\rm QSO})\psi$ (41)
$\displaystyle\dot{M}_{\mathrm{cold,halo}}$ $\displaystyle=$
$\displaystyle-\dot{M}_{\rm cool}$ (42)
$\displaystyle\dot{M}_{\mathrm{hot,halo}}$ $\displaystyle=$
$\displaystyle(\dot{m}_{\mathrm{out,\star}}+\dot{m}_{\mathrm{out,QSO}})$ (43)
$\displaystyle-(\dot{m}_{\mathrm{ejec,\star}}+\dot{m}_{\mathrm{ejec,QSO}})$
$\displaystyle\dot{M}_{\mathrm{ejec}}$ $\displaystyle=$
$\displaystyle\dot{m}_{\mathrm{ejec,\star}}$ (44)
$\displaystyle\dot{M}_{\mathrm{lost}}$ $\displaystyle=$
$\displaystyle\dot{m}_{\mathrm{ejec,QSO}}$ (45)
$\displaystyle\dot{M}_{\star}^{Z}$ $\displaystyle=$
$\displaystyle(1-R)Z_{\mathrm{cold}}\psi$ (46)
$\displaystyle\dot{M}_{\mathrm{cold}}^{Z}$ $\displaystyle=$
$\displaystyle\dot{M}_{\rm cool}Z_{\mathrm{cold,halo}}$ (47)
$\displaystyle+(p-(1+\beta_{\star}+\beta_{\mathrm{QSO}}-R)Z_{\mathrm{cold}})\psi$
$\displaystyle\dot{M}_{\mathrm{cold,halo}}^{Z}$ $\displaystyle=$
$\displaystyle-\dot{M}_{\rm cool}Z_{\mathrm{cold,halo}}$ (48)
$\displaystyle\dot{M}_{\mathrm{hot,halo}}^{Z}$ $\displaystyle=$
$\displaystyle\dot{M}_{\mathrm{hot,halo}}Z_{\mathrm{cold}}$ (49)
$\displaystyle\dot{M}_{\mathrm{ejec}}^{Z}$ $\displaystyle=$ $\displaystyle
Z_{\mathrm{cold}}\dot{m}_{\mathrm{ejec,\star}}$ (50)
$\displaystyle\dot{M}_{\mathrm{lost}}^{Z}$ $\displaystyle=$ $\displaystyle
Z_{\mathrm{cold}}\dot{m}_{\mathrm{ejec,QSO}}$ (51)

where $\beta_{\mathrm{QSO}}\equiv\dot{m}_{\mathrm{out,QSO}}\psi^{-1}$ is the
mass loading due to QSO feedback, $Z_{\mathrm{cold}}\equiv
M_{\mathrm{cold}}^{Z}M_{\mathrm{cold}}^{-1}$ and $Z_{\mathrm{cold,halo}}\equiv
M^{Z}_{\mathrm{cold,halo}}M_{\mathrm{cold,halo}}^{-1}$ are the metallicities
of cold gas in the ISM and the cold gas in the halo (the part actively
cooling), respectively, and $p$ the metal yield. $R$, $\psi$ and
$\beta_{\star}$ were defined when introducing Eqs. (5)-(9). Note that
$\beta_{\rm QSO}$ is only positive in the case of SBs and where
$L_{\mathrm{bol}}>\kappa_{\mathrm{QSO}}L_{\mathrm{crit}}$, and hence during
star formation in the disk, there is no QSO feedback. Note that when gas is
ejected from the halo due to QSO feedback, we assume this gas is lost and
never reincorporated into the halo.

###  3.4 Ram pressure stripping of the halo and ISM gas

Shark v1.1 included a treatment of instantaneous gas stripping of the hot halo
of satellite galaxies. As soon as galaxies became satellites, their hot halo
gas was instantaneously removed and transferred to the hot halo gas of the
central galaxy of the host halo. Here, we include a treatment of ram pressure
stripping (RPS) of both the hot halo and ISM gas of satellite galaxies as
described below.

For the RPS of the hot halo, we follow Font et al. (2008), which follows the
RPS criterion found by McCarthy et al. (2008) in their hydrodynamical
simulations, for a spherical distribution of gas. The halo gas beyond a radius
$r_{\rm sat}$ (measured from the centre of the satellite galaxy) is removed if
the ram pressure at that position exceeds the binding energy the gas feels due
to the satellite’s gravitational potential,

$\rho^{\rm cen}_{\rm halo,gas}(R)\,v^{2}_{\rm sat}>\alpha_{\rm
RPS}\frac{G\,M_{\rm sat}(r_{\rm sat})\,M^{\rm sat}_{\rm halo,gas}}{8\,r^{\rm
sat}_{\rm vir}\,\,r^{3}_{\rm sat}},$ (52)

where $M_{\rm sat}(r_{\rm sat})$ is the total mass enclosed in $r_{\rm sat}$,
$\rho^{\rm cen}_{\rm halo,gas}(R)$ is the central galaxy’s halo gas density at
the position of the satellite galaxy relative to the halo centre, $R$, and
$v_{\rm sat}$ is the velocity of the satellite in the frame of the host halo.
Both $R$ and $v_{\rm sat}$ are the subhalo’s position and velocity for
satellites type 1 (see below for treatment of satellites type 2). The
parameter $\alpha_{\rm ram}=2$ in McCarthy et al. (2008). With the aim of
allowing for flexibility in the code, the latter is left as a free parameter,
which in the default Shark v2.0 is set to $\alpha_{\rm RPS}=1$. Note that
$M^{\rm sat}_{\rm halo,gas}$ and $r^{\rm sat}_{\rm vir}$ are the halo gas mass
and virial radius the satellite galaxy had right before it became a satellite.
The latter assumes that the gas is stripped from the satellite’s halo outside-
in without affecting the density internal to the stripping radius. We find
$r_{\rm sat}$ by assuming equality in Eq. (52) and strip away all the gas that
is at radii $>r_{\rm sat}$ that has not yet been stripped.

For the satellite’s ISM, we assume a similar model for RPS, and remove all the
gas outside $r$, with $r$ being the radius at which

$\rho^{\rm cen}_{\rm halo,gas}(R)\,v^{2}_{\rm sat}=2\pi\,G\Sigma_{\rm
gas}(r)\,\left[\Sigma_{\rm gas}(r)+\Sigma_{\star}(r)\right]$ (53)

where $\Sigma_{\rm gas}$ and $\Sigma_{\star}$ are the ISM’s and stellar
surface densities at $r$. Note that the latter has contributions from both the
disk and bulge components. The ISM in the disk and bulge follows an
exponential profile, with half-mass radii sizes $r_{\rm gas,disk}$ and $r_{\rm
gas,bulge}$. The stellar component of the disk, also follows an exponential
profile with half-mass radius, $r_{\star,\rm disk}$, while the bulge stars
follow a Plummer profile, with half-mass radius $r_{\star,\rm bulge}$. In
Shark v2.0, $r_{\rm gas,disk}$ and $r_{\star,\rm disk}$ are calculated self-
consistently, following the description in $\S$ 3.1, while $r_{\rm
gas,bulge}=r_{\star,\rm bulge}$ and are calculated as described in
$\S\,4.4.12$ of L18. The same profiles are used to calculate $M_{\rm
sat}(r_{\rm sat})$ in Eq. (52).

For the RPS of the ISM, we make the same assumption as for the RPS of the
halo, and assume the gas and stellar profiles internal to $r$ in Eq. (53) are
unchanged under RPS, and hence we keep track of the stripped ISM, so that we
can incorporate it when computing the gas surface density in Eq. (53). We also
assume that the galaxy radii do not change due to RPS. The latter is a
sensible assumption if RPS acts on short timescales (shorter than a dynamical
relaxation timescale).

When a type 1 satellite becomes type 2 (its host subhalo disappears from the
subhalo catalogue), we assume that any left over hot gas is instantaneously
transferred to the central subhalo’s hot gas, but the ISM continues to
experience RPS following the equations above.

###  3.5 Tidal stripping of gas and stars

In addition to the presence of dynamical friction, which ultimately leads to
galaxy mergers, the tidal stripping of these infalling satellite galaxies is
also an important environmental process. Tidally stripped stellar material on
the outskirts of halos typically does not contribute to observational
estimates of the central galaxy stellar mass and instead is associated to an
independent component, which we referred to as “intra-halo stellar mass”. This
can have an important impact on comparing observations of the SMF with model
predictions. Recent work using cosmological hydrodynamical simulations
estimates that up to half of the stellar mass bound to a central subhalo can
be attributed to stellar halos, rather than the central galaxy (e.g. Cañas et
al. 2020; Proctor et al. 2023).

We include in Shark v2.0 a new physical model to describe the tidal effects
felt by satellite galaxies from their central component. This implementation
is based on the work of Errani et al. (2015), who use high resolution $N$-body
simulations of a Milky-Way mass halo to track the evolution of tidal material
that can be associated with dwarf spheroidal galaxies. Although Errani et al.
(2015) did their analysis specifically for a Milky-Way halo, the problem is
highly self-similar and it mostly depends on the fractional mass loss rate
experience by an object orbiting in an (Navarro et al., 1997) profile. Errani
et al. (2015) parameterise the evolutionary tracks of the stripped material
following the work of Peñarrubia et al. (2008),

$\frac{M_{\star,\rm
TS}}{M_{\star,0}}=\frac{2^{\alpha}x^{\beta}}{(1+x)^{\alpha}},$ (54)

where $x=M_{\rm sh}/M_{\rm h,0}$, $M_{\rm h,0}$ is the total halo mass the now
satellite subhalo had at infall, $M_{\rm sh}$ is the current subhalo mass,
$M_{\star,\rm TS}$ is the amount of stellar mass that has been tidally
stripped from the satellite galaxy that is hosted by the subhalo in question,
and $M_{\star,0}$ is the stellar mass the galaxy had at infall. In principle
the halo masses above should not be the total current and at infall ones, but
the mass enclosed within the half-light radii of the galaxy currently or at
infall, respectively. However, Errani et al. (2015) found that the half-light
radius of the galaxy is barely affected by tidal stripping, and hence one can
assume a constant half-light radius to halo’s scale radius ratio. Here we
adopt the parameters $\alpha=3.57$ and $\beta=2.06$, corresponding to a
“cuspy” Navarro et al. (1997) profile and half-light radius to halo’s scale
radius ratio $=0.2$. These values correspond to the largest effect tidal
stripping can have on the stellar content of a satellite galaxy according to
Errani et al. (2015). This allows us to test the maximal effect tidal
stripping can have on our galaxy population. We also implement a minimum
remaining halo mass fraction of 1% (where the fraction gives the difference
between the subhalo mass at the current snapshot and the subhalo mass at
infall; this is a parameter in Shark, named minimum_halo_mass_fraction) to
ensure that a halo can only lose at most 99% of its mass. This model can be
used by setting tidal_stripping = true in the input parameter file.

We implement this in a way that the disk, which has the larger radius, is
stripped first, followed by the bulge. If the galaxy has a cold gas reservoir,
we strip the amount of gas that is outside a radius $r_{\rm str}$ beyond which
a fraction of stellar mass $M_{\star,\rm TS}/M_{\star,0}$ has been stripped.

Table 2: List of new models and parameters included in Shark v2.0. The parameters presented in purple are those that are varied in optim when searching for a best fit to the $z<0.1$ SMF. We present a suggested range for the parameters and in parentheses the value adopted in our best-fitting model. The relevant equation/section for each model or parameter is presented in the third column. Parameters that were introduced in L18 refer to the equations in that paper. Parameter | Suggested value range (adopted) | Equation/Section
---|---|---
Physical Process | AGN Feedback and BH growth |
model | bower06, croton15, lagos22 (lagos22) |  $\S$ 3.3
spin-model | constant, volonteri07, griffin19 (griffin19) |  $\S$ 3.3.1
accretion-disk-model | prolonged, selfgravitydisk, warpeddisk (warpeddisk) | Appendix B (only for spin-model=griffin19)
qso-feedback | true or false (true) |  $\S$ 3.3.3
$\kappa$ |  $10^{-5}-10^{2}$ (10.31) | Eq. (11)
$f_{\rm smbh}$ |  $10^{-5}-10^{-2}$ ($3\times 10^{-3}$) | Eq. (12)
$e_{\rm sb}$ |  $0.5-50$ (15) |  $\S$ 3.3.1
$\dot{m}_{\rm ADAF}$ | 0.01 |  $\S$ 3.3.1
$\delta_{\rm ADAF}$ | 0.1-0.5 (0.2) | Eq. (20)
$\alpha_{\rm ADAF}$ |  $0.05-0.5$ (0.1) | Eq. (20)
$\eta$ | 1-10 (4) |  $\S$ 3.3.1
$\alpha_{\rm TD}$ |  $0.05-0.5$ (0.1) | Appendix B
$\kappa_{\rm radio}$ |  $10^{-1}-3$ ($0.023$) | Eq. (32)
$\Gamma_{\rm thresh}$ |  $0.01-100$ (10) |  $\S$ 3.3.2
$\epsilon_{\rm QSO}$ |  $0.1-100$ (10) | Eq. (37)
Physical Process | Tidal and Ram-pressure Stripping |
tidal-stripping | true or false (true) |  $\S$ 3.5
gradual-stripping-ism | true or false (true) |  $\S$ 3.4
gradual-stripping-halo | true or false (true) |  $\S$ 3.4
$\alpha_{\rm RPS}$ | 0.1-10 (1) | Eq. (52)
minimum-halo-mass-fraction |  $10^{-5}-10^{-1}$ ($10^{-2}$) |  $\S$ 3.5
Physical Process | Dynamical Friction |
merger-timescale-model | lacey93 or poulton20 (poulton20) |  $\S$ 3.2
Physical Process | Chemical Enrichment |
evolving-yield | true or false (true) | Appendix A.1
Physical Process | Star Formation |
$\nu_{\rm SF}$ |  $0.2-1.7\,\rm Gyr^{-1}$ ($1.49\,\rm Gyr^{-1}$) | Eq. (7) in L18
angular-momentum-transfer | true or false (true) |  $\S$ 3.1
Physical Process | Stellar Feedback |
$\beta$ |  $0.5-5$ (3.79) | Eqs. (25)-(28) in L18
$\beta_{\rm min}$ |  $0.01-1$ (0.104) | Appendix A.2
$v_{\rm hot}$ |  $50-500\,\rm km\,s^{-1}$ ($120\,\rm km\,s^{-1}$) | Eqs. (25)-(28) L18
Physical Process | Reincorporation |
$\tau_{\rm reinc}$ |  $1-30$ Gyr (21.53 Gyr) | Eq. (30) in L18
$\gamma$ |  $-3\,\rm to\,0$ (-2.339) | Eq. (30) in L18
$M_{\rm norm}$ |  $10^{9}-10^{12}\,\rm M_{\odot}$ ($1.383\times 10^{11}\,\rm M_{\odot}$) | Eq. (30) in L18

###  3.6 Automatic parameter exploration

We perform an automatic search for a suite of best-fitting parameters by using
a Particle Swarm Optimisation (PSO) python package, optim, introduced in
Proctor et al. (in preparation). We here use the $z<0.1$ SMF of Li & White
(2009) based on SDSS as our primary constraint. In addition to this, optim can
also use as constraints the SMF at $z=0.5,\,1$ and $2$ (from Weaver et al.
2022), the cosmic star formation rate density (CSFRD) of Driver et al. (2018),
the HI mass function of Jones et al. (2018) and the total stellar-size mass
relation of the Galaxy and Mass Assembly (GAMA) survey (see Appendix D).
Proctor et al. (in preparation) introduce in detail optim and these
constraints. Note that optim can be easily used as a standalone tool, and
hence be adapted to work with other SAMs.

Following Proctor et al. (in preparation), we only vary the following
parameters:

  * •

AGN feedback parameters: $\kappa$ (Eq. 11), $\kappa_{\rm radio}$ (Eq. 32), and
$\Gamma_{\rm thresh}$ ($\S$ 3.3.2). These parameters control the efficiency of
gas accretion onto the BH in the hot-halo mode, the efficiency of coupling
between the mechanical power of the AGN and the hot halo of the galaxy, and
the cooling-to-heating specific energy ratio threshold above which hot halos
form, respectively.

  * •

Star formation parameters: $\nu_{\rm SF}$ (Eq. (7) in L18), which controls the
conversion efficiency between the surface density of molecular gas and SFR.

  * •

Stellar feedback parameters: $\beta$ (Eqs. (25)-(28) in L18) and $\beta_{\rm
min}$ (Appendix A). These parameters control the dependence of the outflow
mass loading on the circular velocity of the galaxy, and the minimum value the
mass loading can take, respectively.

  * •

Disk instability parameters: $\epsilon_{\rm disc}$ (Eq. (35) in L18), which
controls the threshold of the stability parameter, below which disks are
considered globally unstable to collapse.

  * •

Gas reincorporation parameters: $\tau_{\rm reinc}$, $\gamma$ and $M_{\rm
norm}$ (Eq. (30) in L18), control how quickly gas that has been removed from
halos due to stellar feedback can be reincorporated into the halo. Once gas is
reincorporated is available for cooling.

Other parameters are left fixed, following the values in Table 2. The ranges
in which we vary the parameters above and the best-fitting values are
presented in Table 2.

Even though we only use the $z=0$ SMF as a numerical constraint for optim, we
visually inspected other results of the model to ensure we were not seriously
compromising their agreement with observations by using the $z=0$ SMF only.
For example, we inspected gas scaling relations, the mass-metallicity relation
and the $z<2$ CSFRD and made sure they look sensible. Some of these results
are presented in the supplementary material.

Note that the above does not mean that we change the parameters from the best-
fitting ones, but that we used the visual inspection of other results to draw
reasonable prior ranges in some of the parameters that we input to optim. An
example of that, is that we limited the stable parameter (see Table 2) to a
maximum of $1.2$. Higher values lead to unphysically large numbers of bulge-
dominated galaxies.

###  3.7 Other technical updates

Many small bug fixes and improvements have been made in the Shark codebase
since version 1.1.0 was released. Below, we highlight the most important ones:

  * •

The memory footprint of Shark has been considerably reduced by reorganising
some of its internal structures, with improvements of about $20$%, depending
on the input data. Shark is usually memory bound, and thus this reduction
allows more Shark instances to run in the same set of resources.

  * •

Executions are now fully reproducible. By default a random seed is used, but
one can be given. The seed used by each run is logged and stored in the output
files.

  * •

Performance has been improved in two fronts: first, we now perform better load
balancing when running in multi-threaded mode, reducing walltime; secondly, an
overall reduction of small, temporary memory allocations allow Shark to run
more streamlined.

  * •

Improved infrastructure: dropped requirement for the HDF5 C++ libraries, moved
our automatic per-commit checks from Travis CI to GitHub actions, added
support for more systems and compilers (e.g. MSVC in Windows), more runtime
information is written (timing, memory usage, etc).

##  4 Characteristics of the galaxy population: abundances and scaling
relations

Figure 1: Galaxy SMF in Shark v2.0 (solid red line) and the default model in
L18 (dashed black line), from $z=0$ to $z=8$, as labelled. We show
observations from Li & White (2009); Bernardi et al. (2013); Driver et al.
(2022) at $z=0$, Muzzin et al. (2013); Thorne et al. (2021); Weaver et al.
(2022) at $0.5\leq z\leq 4$, and from Weaver et al. (2022) only at $z\geq 5$,
as labelled. Li & White (2009), shown in green in the top-left panel, is the
constraint used in optim. Shark v2.0, we also show the SMF using the stellar
mass in galaxies enclosed in an aperture of $30$ kpc (red dotted line, shown
only at $z=0$), and when we apply a random error in the stellar mass of $0.3$
dex at $z\geq 1$ (red dot-dashed line). Figure 2: Galaxy baryon mass function
in Shark v2.0 (red lines) and the default model in L18 (thin black lines), as
labelled, at $z=0$ for galaxies in $4$ bins of host halo mass, as labelled in
each panel. We separate the contribution of central (solid lines) and
satellite (dashed lines) galaxies in each panel. Observations from the galaxy
surveys RESOLVE and ECO from Eckert et al. (2016) are shown as symbols for
central galaxies (blue symbols) and satellites (yellow symbols). The baryon
mass here is defined in the same way as in Eckert et al. (2016), $\rm
M{{}^{\prime}}_{\rm bar}=M_{\star}+1.4\,M_{\rm HI}$, with $M_{\rm HI}$ being
the atomic hydrogen mass of the galaxy.

In this section we discuss fundamental observed properties of galaxies,
comparing the results of the new Shark v2.0 model using the best-fitting
parameters of Table 2 with the default model presented in L18.

###  4.1 Galaxy properties: abundances

####  4.1.1 The stellar mass function

Fig. 1 shows the SMF from $z=0$ to $z=7$ comparing with observations. For
$z=0$ we show the observational estimates of Li & White (2009) which were used
for the parameter tuning, and those of Bernardi et al. (2013); Driver et al.
(2022) to show the large systematic uncertainties that permeate the SMF,
especially at the high-mass end. By construction, Shark v2.0 provides a good
fit to the observations of Li & White (2009), better than was achieved in L18.
This is not only due to the improved automatic parametrisation of the model,
but also to the improved physics. Proctor et al (in preparation) show that
when using Shark v1.1 in tandem with the PSO optim package, we fail to provide
a fit as good at the one found here for Shark v2.0. This boils down to the
inability of the default AGN feedback model adopted in L18 to produce a steep-
enough high-mass end of the SMF.

At $z\geq 0.5$, Shark v2.0 tends to predict a steeper low-mass end of the SMF
compared to L18, due to the adopted parameters leading to a weaker dependence
of the mass loading on the circular velocity of the galaxy (see Fig. 16). At
the high-mass end, Shark v2.0 produces a steeper break of the SMF than L18,
and lower-mass galaxies at fixed number density. At $z\geq 2$, the high-mass
end in Shark v2.0 produces too low number densities of massive galaxies,
$>10^{11}\,\rm M_{\odot}$, but this is easily alleviated if we assume stellar
masses have uncertainties which are Gaussian-distributed with a width of
$\approx 0.3$ dex (see dot-dashed red lines). Uncertainties of that magnitude
are quite common in deriving stellar masses due to the required assumptions
(e.g. SFH, metallicity history, initial mass function, among others) - see
Marchesini et al. (2009) for a thorough list of uncertainties and Robotham et
al. (2020) for a recent example of the uncertainties associated to $z\approx
0$ stellar masses.

At $z\geq 5$, the SMF starts to behave more like a power-law rather than a
Schechter function, lacking a clear stellar mass break. This is accentuated
when we include random errors to the Shark stellar masses. Observations appear
to behave similarly, with a clear stellar mass break present only up to
$z\approx 3-4$. Note also that the low-mass end of the SMF becomes
increasingly steeper with increasing redshifts in Shark (both v2.0 and L18).
The observations appear to follow a similar trend. Better observational
constraints at the low-mass end, however, are needed to confirm this trend.

####  4.1.2 The baryon mass function of satellite and central galaxies

Figure 3: Cosmic SFR (top panel) and stellar mass (bottom panel) density as a
function of redshift. We show the total contribution of all galaxies in black,
galaxy disks in blue and galaxy bulges in red, for the Shark v2.0 default
model (solid lines) and the default model in L18 (dashed lines). In the top
panel we show the observations from Driver et al. (2018); D’Silva et al.
(2023) and the compilation of Adams et al. (2023), as labelled. The latter
includes results from Oesch et al. (2018); Harikane et al. (2023); Bouwens et
al. (2023b, a); Donnan et al. (2023); Pérez-González et al. (2023) and their
own measurements, most of them using the JWST. In the bottom panel we show
observational constraints from Madau & Dickinson (2014); Driver et al. (2018);
Weaver et al. (2022); Santini et al. (2023). Data from Madau & Dickinson
(2014) has been re-scaled to a Chabrier IMF.

Fig. 2 shows the baryon mass function in bins of halo mass, and separating
centrals and satellite galaxies at $z=0$. To allow for a fair comparison with
the observations of Eckert et al. (2016), we define the baryon mass as the sum
of the stellar mass and the HI mass times $1.4$. The latter factor was
introduced by Eckert et al. (2016) to account for the Helium contribution.
Note that Eckert et al. (2016) ignored the contribution of molecular and
ionised gas to the baryon mass. This is a reasonable assumption at $z=0$ given
the small molecular gas mass fractions in gas-rich galaxies (which tend to be
low-mass galaxies) and the overall low gas fractions in the regime where
molecular and atomic gas make a comparable contribution to the gas content of
galaxies (i.e. massive galaxies) (Catinella et al., 2018). The observations of
Eckert et al. (2016) were presented for two complementary surveys RESOLVE and
ECO; the former is more sensitive but of smaller area than the latter. We show
the results from these two surveys using different symbols.

In low-mass groups, $11<{\rm log}_{10}(M_{\rm halo}/\rm M_{\odot})<11.4$ and
$11.4<{\rm log}_{10}(M_{\rm halo}/\rm M_{\odot})<12$, we see that Shark v2.0
is able to produce a much higher number density (by an order of magnitude) of
satellite galaxies across the dynamic range probed by the observations
compared to L18. This is in large part driven by the inclusion of the new
dynamical friction timescale of Poulton et al. (2021), which leads to much
longer dynamical friction timescales in low-mass galaxies than what is
obtained by using Lacey & Cole (1993), allowing them to survive for longer and
exist in bigger quantities. In addition, Shark v2.0 recovers a baryon mass
function for central galaxies that matches the high-mass end of the
observations much better than L18. This has to do with the modified parameters
of stellar feedback, which make it overall less efficient in the new model.

The overall higher number density of satellite galaxies in Shark v2.0 compared
to L18 is still present in halos of masses $12<{\rm log}_{10}(M_{\rm halo}/\rm
M_{\odot})<13.5$ but the difference is much smaller than that seen at lower
halo masses. This is because the Poulton et al. (2021) dynamical friction
timescales get closer to those of Lacey & Cole (1993) for more massive
satellites. In the highest halo mass bin, ${\rm log}_{10}(M_{\rm halo}/\rm
M_{\odot})>13.5$, satellite galaxies in Shark v2.0 display a sharper high-end
cut off than L18, which agrees better with observations. For central galaxies,
the main difference is the slightly higher baryon masses produced by Shark
v2.0 in the $12<{\rm log}_{10}(M_{\rm halo}/\rm M_{\odot})<13.5$ halo mass
bins. Observational errors are quite large in this regime, so within the
uncertainties, both Shark v2.0 and L18 agree with observations reasonably
well. Note that here we have not included any potential confusion between the
tagging of central and satellite galaxies, which tends to be rather large,
especially in low-mass halos (see discussion in Bravo et al. 2020).

###  4.2 Star formation rate and stellar mass cosmic density evolution

Fig. 3 shows the cosmic SFR and stellar mass density evolution from $z=0$ to
$z=15$. Observational constraints from Madau & Dickinson (2014); Driver et al.
(2018); D’Silva et al. (2023); Adams et al. (2023); Weaver et al. (2022);
Santini et al. (2023) are shown. For the CSFRD, Shark v2.0 predicts higher
densities at $z\gtrsim 3$ than L18, driven by both the contribution from SBs
and star formation in disks being higher in the new model than L18. For the
CSFRD in disks, we see that at $z\gtrsim 8$, Shark v2.0 and L18 predict very
similar levels, but the differences remain for the CSFRD in SBs from $z\approx
3.5$ even up to $z=15$. The overall qualitative trend of SBs dominating the
CSFRD at $z\gtrsim 3$ that was presented in L18 is maintained in Shark v2.0.
The new model compares more favourably with observations in general, giving a
higher CSFRD at $z\gtrsim 4$, closer to to current JWST estimates beyond $z=8$
within the uncertainties. At $z\lesssim 4$ Shark v2.0 and v1.1 have similar
CSFRDs.

The overall higher CSFRD at high redshift in Shark v2.0 compared to L18
naturally leads to a cosmic stellar mass density (CSMD) in the early universe
that is higher in Shark v2.0 than L18. It is interesting that just the SB
contribution to the CSFRD and CSMD in Shark v2.0 is similar or even higher
than the total in L18 at $z\gtrsim 7$. This shows that these predictions are
still subject to large variations and that better constraints from the high
redshift universe are needed to narrow down the parameter space that is
plausible in the early universe.

###  4.3 Galaxy scaling relations

####  4.3.1 The stellar-halo mass relation

Figure 4: Top panel: The stellar-halo mass relation for central galaxies in
Shark v2.0 (solid black line) and L18 (dashed black line) at $z=0$. Lines show
the median relation for bins with $\geq 10$ galaxies. For Shark v2.0 we also
show the $16^{\rm th}-84^{\rm th}$ percentile ranges as the shaded region. For
reference, we show the empirical constraints from Moster et al. (2013);
Behroozi et al. (2013), and the more direct observational constraints from
Kravtsov et al. (2018); Taylor et al. (2020); Romeo et al. (2020), as
labelled. Bottom panel: As in the top panel but splitting galaxies between
bulge- (red) and disk-dominated (blue) galaxies, using a disk-to-total stellar
mass ratio (D/T) of $0.5$. Galaxies above (below) the threshold are considered
disk-(bulge-)dominated. Observational constraints from Kravtsov et al. (2018)
and Correa & Schaye (2020) for late- and early-type galaxies, LTGs and ETGs,
respectively, are also shown, as labelled.

The top panel of Fig. 4 shows the stellar-halo mass relation for central
galaxies at $z=0$ for Shark v2.0 and L18. At the low-halo mass end, $\rm
log_{10}(M_{\rm halo}/M_{\odot})<12$, Shark v2.0 produces slightly higher
stellar masses at fixed halo mass than L18, in better agreement with empirical
constraints from abundance matching. This is the direct result of the stellar
feedback parameters adopted in Shark v2.0 compared to L18 (see Table 2). At
$\rm 12<log_{10}(M_{\rm halo}/M_{\odot})<14$, Shark v2.0 produces smaller
stellar masses at fixed halo mass than L18, while at higher halo masses both
models converge to similar values. The difference seen is the result of the
new AGN feedback model leading to a more efficient quenching of central
galaxies in massive halos in Shark v2.0 (discussed in detail in § 5). Note
that this is not just due to the choice of parameters, but the capability of
the new AGN feedback model to more efficiently quench massive galaxies
overall. Proctor et al. (in preparation) show that for the AGN feedback, the
implementation introduced in L18 was not capable of getting a very sharp SMF
break at the high-mass end, while the new model allows for that.

Both Shark v2.0 and L18 produce a stellar-halo mass relation at $\rm
log_{10}(M_{\rm halo}/M_{\odot})>12$ that is steeper than the empirical
relations of Moster et al. (2013); Behroozi et al. (2013). However, compared
to the more direct estimates of Kravtsov et al. (2018), Shark compares
favourably, especially at the galaxy cluster regime. Kravtsov et al. (2018)
argued that the difference with the empirical relations of Moster et al.
(2013); Behroozi et al. (2013) is due to those studies using SMFs that are
severely surface brightness-limited, while the Kravtsov et al. (2018) data
reaches lower surface brightness values, impacting the recovered stellar mass.
We also show the observational measurements of Taylor et al. (2020) combining
weak lensing with stellar mass measurements from the GAMA survey. These again
are significantly different than the empirical relations of Moster et al.
(2013); Behroozi et al. (2013), showing that there are still many systematic
uncertainties in the stellar-halo mass relation which are usually not captured
in the presented errorbars. Individual local Universe galaxies (primarily of
late-type morphology) compiled by Romeo et al. (2020) are also shown,
depicting the significant galaxy-to-galaxy scatter of the relation.

The bottom panel of Fig. 4 shows the stellar-halo mass relation of central
galaxies in Shark v2.0 and v1.1 for disk- and bulge-dominated galaxies,
separately. We define these populations using a disk-to-total stellar mass
threshold of $0.5$. In both Shark versions, and at fixed halo mass, bulge-
dominated galaxies have higher stellar masses than disk-dominated ones at
$M_{\rm halo}<10^{12.5}\,\rm M_{\odot}$. This is due to bulge-dominated
galaxies forming earlier than disk-dominated galaxies and having had more
starbursts during their lifetimes. At $M_{\rm halo}>10^{12.5}\,\rm M_{\odot}$,
the trend reverses in Shark v2.0 but not in v1.1, with disk-dominated galaxies
in v2.0 being more massive at fixed halo mass than bulge-dominated galaxies.
This happens because AGN feedback in v2.0 is not efficient enough in the disk-
dominated galaxies of these high halo masses, so they form stars very
efficiently and end up more massive than the bulge-dominated ones. In Shark
v1.1, the efficiency of AGN feedback is more directly tied to the halo mass
than in the new model, leading to the differences seen at $M_{\rm
halo}>10^{12.5}\,\rm M_{\odot}$.

Interestingly, Correa & Schaye (2020) reported that in the eagle cosmological
hydrodynamical simulations, disk-dominated galaxies have a higher stellar
masses than bulge-dominated galaxies at all halo masses. We find that the
difference with the trend in Shark v2.0 is driven by disk instabilities in the
model. If we completely switch off disk instabilities, we find that the
resulting disk-dominated galaxies have on average higher stellar masses than
bulge-dominated ones at fixed halo mass (not shown here), trend that
qualitatively agrees with eagle. This shows that the exact treatment of disks
during these episodes of global instabilities has an important impact on the
resulting stellar-halo mass relation. This agrees with the conclusions of
Romeo et al. (2020) who argue that disk gravitational instabilities regulate
the stellar-to-halo mass relation in galaxies in the halo mass range of their
sample (see orange stars in Fig. 4).

Observations still have uncertainties that are too large to be able to
confidently verify the conflicting predictions of the stellar-halo mass
relation by galaxy morphology. In fact, Correa & Schaye (2020) show, using
data from the Sloan Digital Sky Survey (SDSS), that depending on how stellar
masses and morphologies are measured, one can get different offsets between
bulge- and disk-dominated galaxies in the stellar-halo mass relation. Future
surveys, such as the 4MOST WAVES (Driver et al., 2019), will provide exquisite
measurements of halo masses down to $\rm log_{10}(M_{\rm
halo}/M_{\odot})\approx 12$ in the local Universe, which will allow a robust
diagnosis of the predictions discussed here.

Figure 5: Size-stellar mass relation for the disk (top panel) and bulge
(middle panel) components of galaxies, and for the entire galaxy in the bottom
panel, at $z=0$, for Shark v2.0 (solid lines) and L18 (dashed lines). In the
top two panels, we plot disk (bulge) half-stellar mass radii vs. disk (bulge)
stellar mass, as labelled, while the bottom panel shows the galaxy half-
stellar mass radius vs total stellar mass. We show bins with $\geq 10$
objects. Lines with the shaded regions show the medians and $16^{\rm
th}-84^{\rm th}$ percentile ranges, respectively. Dotted lines in the top two
panels show the $50^{\rm th}$, $68^{\rm th}$ and $90^{\rm th}$ percentile
regions of the GAMA observations of Lange et al. (2016), while the symbols
show the more recent GAMA results of Robotham et al. (2022); Bellstedt et al.
(2023) (see Appendix D for details). For the latter we show the medians as
symbols, and the $16^{\rm th}-84^{\rm th}$ percentile ranges as errorbars. The
bottom panel shows the mass-weighted stellar size-mass relation of Robotham et
al. (2022); Bellstedt et al. (2023) (see Appendix D for how this was
computed). Figure 6: sAM of the stars as a function of stellar mass (left
panel), sAM of the total ISM and atomic gas a function of total ISM mass
(middle-left panel), sAM of all baryons in the galaxy vs baryon mass (middle-
right panel), and sAM of atomic, molecular gas and stars vs stellar mass
(right panel), for galaxies in Shark v2.0 at $z=0$. In the first three panels
we show the median relations for galaxies with varying disk-to-total stellar
mass ratios, as labelled in the third panel. The shaded region shows the
$16^{\rm th}-84^{\rm th}$ percentile ranges for galaxies with $\rm D/T>0.5$
only. The fourth panel shows the median sAM-stellar mass relations for
galaxies with $\rm D/T>0.25$. In this panel we also show the median relations
of Shark v1.1 (L18) for $j_{\star}$ and $j_{\rm HI}$ (in this model $j_{\rm
H_{2}}\equiv j_{\rm HI}$), as dotted lines. We show observations from
Obreschkow & Glazebrook (2014); Chowdhury & Chengalur (2017); Posti et al.
(2018); Murugeshan et al. (2019); Mancera Piña et al. (2021); Di Teodoro et
al. (2021); Hardwick et al. (2022b, a), as labelled. We highlight the
measurements of Obreschkow & Glazebrook (2014) (OG14) in all panels using
coloured symbols - this is the only sample where measurements of the sAM of
stars, atomic and molecular gas are presented for the same galaxies (shown
with different colours in the right-hand panel, as labelled).

####  4.3.2 Structural properties of galaxies at $z=0$

One of the major modifications in Shark v2.0 is the more sophisticated
modelling of galaxy angular momentum, and the exchange of angular momentum
between the ISM and stars ($\S$ 3.1). Hence, it is natural to explore what the
effect of that is on structural scaling relations of galaxies. Here, we focus
on the size-mass and the specific angular momentum-mass relations at $z=0$.

We start with the size-mass relation presented in Fig. 5. The latter is shown
for the disk and bulge components of galaxies at $z=0$ for Shark v2.0 and L18,
compared with observations from Lange et al. (2016); Robotham et al. (2022);
Bellstedt et al. (2023), in the top two panels. The bottom panel shows the
half-stellar mass radius, $r_{\star}$, vs galaxy stellar mass compared with
the stellar-mass weighted half-mass radius vs stellar mass of Robotham et al.
(2022); Bellstedt et al. (2023). Note that Lange et al. (2016); Robotham et
al. (2022); Bellstedt et al. (2023) use the GAMA survey. In Appendix D we
explain why these measurements are different and the process of measuring the
size-mass relation from the Bellstedt et al. (2023) catalogue. Note that for
both Shark and the observations we use the same Eq. (69) to calculate
$r_{\star}$.

The disk sizes in Shark v2.0 are smaller than L18, especially at the low-mass
end. This is the regime where we expect our new angular momentum treatment to
make a difference. As these galaxies have molecular gas that is more
concentrated than the atomic gas, the new model would produce new stars that
form at overall smaller radii than the old model, which assume all disk
components to have the same angular momentum distribution. The lower disk
sizes in Shark v2.0 agree better with the observations than v1.1 (L18).

The bulge size-mass relation (middle panel in Fig. 5) in Shark v2.0 is steeper
than in L18 with a plateau around a bulge mass of $10^{10}\,\rm M_{\odot}$,
which is roughly when galaxies transition from being disk- to bulge-dominated.
The steeper size-mass relation agrees better with observations than the
prediction of L18, especially on the regime of massive spheroids. Note that
there is a clear systematic difference between the bulge-size mass relations
in the observations shown. This is because defining a bulge and its size is
more difficult than for disks in observations (see discussion in Robotham et
al. 2022).

The galaxy size-mass relation in the bottom panel of Fig. 5 shows that Shark
v2.0 has a clear transition at a stellar mass of $10^{9.7}-10^{10}\,\rm
M_{\odot}$. At lower stellar masses, the disk size dominates $r_{\star}$,
while at higher masses it is the bulge size that determines $r_{\star}$. The
transition zone was not clearly present in L18, in which the size-mass
relation was better described by a single power law. Compared with
observations, Shark v2.0 is in better agreement, although the transition
region in observations starts at lower masses. Note, however, that this is not
part of the diagnostics used to tune the parameters of the model, so the
improvement in the agreement with observations is a success of the new angular
momentum model.

Fig. 6 shows the stellar specific angular momentum (sAM) vs stellar mass,
atomic gas and all ISM sAM vs gas mass, baryon sAM vs baryon mass, and the sAM
of stars, atomic and molecular gas as a function of stellar mass. For the
first three cases, we show the relation for galaxies of varying disk-to-total
stellar mass ratio, D/T, as labelled in the third panel, while for the fourth
panel we show the sAM-stellar mass relations for galaxies with a $\rm
D/T>0.5$. We compare with a large compilation of observations as labelled.

The stellar sAM-mass relation (left panel of Fig. 6) shows a strong dependence
of the zero-point of the relation on D/T, so that galaxies with lower D/T have
a lower sAM at fixed stellar mass. This is the well known morphological
dependence of the sAM-mass relation, first discussed in Fall (1983). Overall,
Shark v2.0 agrees very well with the observations, and even the scatter for
disk-dominated galaxies ($\rm D/T>0.5$) is similar to the scatter reported in
observations, which for the most part contain disk galaxies only. The baryon
sAM-mass relation (middle-right panel of Fig. 6), shows a similar behavior,
but the scatter increases with increasing baryon mass in a way that resembles
what observations report. The difference in baryon sAM between galaxies with
$\rm D/T>0.75$ and $\rm D/T>0$ is larger than the differences obtained in
stellar sAM. The atomic (or total ISM) sAM-mass relation (middle-left panel of
Fig. 6) is the tightest of all, in agreement with what has been reported in
observations. However, we see that Shark v2.0 produces a slightly too shallow
relation. This slope directly depends on what we assume for the halo spin
parameter. L18 assumed a halo spin distribution with a mode that was
independent of halo mass. In Shark v2.0, we instead assume a very weak
dependence of the mode of the spin distribution on halo mass as reported in
Kim et al. (2015). A slightly steeper relation would help the model reproduce
the observed slope, however, we decide not to change the spin-halo mass
dependence arbitrarily. Ideally, we would want to use the spin parameter
directly inferred from the $N$-body simulation. However, such measurements are
very noisy and only reliable for halos resolved with several hundred
particles, which prohibits the use of that in SAMs, in which we use halos
sampled with a number of particles as low as $20$.

The right panel of Fig. 6 illustrates the power of the new angular momentum
exchange model ($\S$ 3.1) by showing the difference between the sAM of the
atomic and molecular gas, and the stars at fixed stellar mass. The offset
between the sAM of atomic and molecular gas increases slightly as the stellar
mass decreases. We also show observations from Obreschkow & Glazebrook (2014),
which include $j_{\star}$, $j_{\rm HI}$, and $j_{\rm H_{2}}$ for the same
sample of galaxies (THINGS, The HI Nearby Galaxy Survey). The difference
between $j_{\star}$ and $j_{\rm HI}$ in THINGS is $\approx 0.3$ dex, while in
Shark v2.0 is closer to $0.5$ dex. In Shark v2.0, $j_{\rm H_{2}}$ is $\approx
0.15$ dex higher than $j_{\star}$ at fixed stellar mass, while in THINGS this
difference is $\approx 0.05$ dex. Nevertheless, within the scatter of the
observations, our predictions agree well. We also show the relations for Shark
v1.1 (L18) in dotted lines, and find clear disagreements with the
observations, which are most evident for $j_{\rm HI}$, which in Shark v1.1 is
$\approx 0.2-0.6$ dex lower than Obreschkow & Glazebrook (2014). Note that
$j_{\star}$ and $j_{\rm HI}$ in Shark v1.1 differ only by $0-0.15$ dex, while
observations consistently prefer $\approx 0.3$ dex across the whole stellar
mass range probed. By definition, $j_{\rm H_{2}}=j_{\rm HI}$ in Shark v1.1,
while in Shark v2.0 $j_{\rm HI}$ is higher than $j_{\rm H_{2}}$ by $0.2-0.5$
dex (with the difference increasing at lower stellar masses), which is similar
to the differences seen in THINGS. The improvements we see in the sAM-mass
scaling relations are a direct result of the new angular momentum exchange
model ($\S$ 3.1).

Overall we see that Shark v2.0 with the default parameters adopted (Table 2)
produces structural scaling relations that are in broad agreement with
observations in the local Universe. We can attribute that in large part to the
new angular momentum exchange model. In the future we also plan to explore
other structural scaling relations, such as the size-mass relation across
cosmic time, the Tully-Fisher relation, and the relation between sizes-mass
and galaxy age reported in Robotham et al. (2022).

####  4.3.3 The BH-stellar mass relation and morphological dependence

Figure 7: BH-stellar mass relation of all galaxies in Shark v2.0 at $z=0$ and
the subsample of galaxies classified as ETGs (those with a D/T$<0.5$) and LTGs
(those with D/T$\geq 0.5$), as labelled. We also show the relation for ETGs
and LTGs in Shark v1.1 (L18) as dashed lines. Observations of ETGs and LTGs
from Sahu et al. (2019) and Davis et al. (2019), respectively, and from
passive (P) and star-forming galaxies (SF) from Terrazas et al. (2017) are
shown with symbols, as labelled. The BH-stellar mass relation of P and SF
galaxies in Terrazas et al. (2017) agree well with those of ETGs and LTGs,
respectively, of Sahu et al. (2019) and Davis et al. (2019).

Fig. 7 shows the BH-stellar mass relation at $z=0$ for Shark v2.0 and v1.1
(L18), split by galaxy morphology. We do the latter classification based on
D/T, with D/T$<0.5$ corresponding to early-type galaxies (ETGs) and D/T$\geq
0.5$ to late-type galaxies (LTGs). Both Shark v2.0 and v1.1 predict LTGs to
have systematically lower BH masses than ETGs at fixed stellar mass. The
offset in the zero-point of the BH-stellar mass relation between the LTG and
ETG populations is larger in Shark v2.0 ($\approx 0.5-1.5$ dex) than in v1.1
($\approx 0.3-0.6$ dex), especially clear at $M_{\star}\lesssim 10^{10.5}\,\rm
M_{\odot}$. Another interesting difference is that Shark v2.0 has LTGs as
massive as $10^{12.2}\,\rm M_{\odot}$ while in v1.1 they at most reach stellar
masses of $\approx 10^{11.5}\,\rm M_{\odot}$.

Observations show that the typical BH mass difference between LTGs and ETGs is
$\approx 1$ dex at fixed stellar mass. Shark v2.0 predicts a difference of
$\approx 0.8-1$ dex in agreement with the observations, while v1.1 prefers
smaller differences of $\approx 0.3-0.6$ dex. Both Shark versions predict the
BH-stellar mass relation of LTGs to have a larger scatter than that of ETGs,
however, this difference in scatter between morphological types is larger in
Shark v2.0 than v1.1. Observations suggest a qualitatively similar difference,
with the scatter of the BH-stellar mass relation of LTGs or star-forming
galaxies being larger than for ETGs or passive galaxies. We also find that
overall the BH-stellar mass relation becomes tighter for both galaxy
populations with increasing stellar mass. More observations are needed to
establish whether that is also the case in the real Universe.

##  5 Quenching of massive galaxies

One of the major developments in Shark v2.0 compared with L18 is the new AGN
feedback model ($\S$ 3.3). As such, we focus here on analysing the effect this
has on the quenching of galaxies, especially massive galaxies. This is also an
area that has attracted a lot of attention due to recent results from the
JWST, which point to massive-quiescent galaxies being more common than
previously thought of at $z>3$ (Carnall et al., 2023; Nanayakkara et al.,
2022; Valentino et al., 2023; Long et al., 2023).

This section is organised as follows: $\S$ 5.1 focuses on the quenching of
galaxies in the local universe, and $\S$ 5.2 analyses how quenching develops
across cosmic time.

###  5.1 Quenching in local universe galaxies

Figure 8: SFR as a function of stellar mass for galaxies at $z=0$. The black
solid and dashed lines show the medians for the galaxies in Shark v2.0 and
v1.1 (L18), respectively, for bins with $\geq 10$ objects. The filled contours
show percentile ranges ranging from $99^{\rm th}$ to $10^{\rm th}$, from the
outer to the inner regions, for Shark v2.0. The brown, dashed line shows the
median SFR-stellar mass relation of SDSS galaxies as reported by Brinchmann et
al. (2004); the blue dot-dashed line shows the main sequence as measured from
GAMA by Davies et al. (2016); the pink dashed line and thin solid lines show
the median and the $99^{\rm th}$, $95^{\rm th}$, $68^{\rm th}$ and $50^{\rm
th}$ percentile ranges of GAMA galaxies from Bellstedt et al. (2023).
Observations of massive galaxies from Terrazas et al. (2017) are also shown as
red symbols as a reference of typical SFRs of massive galaxies that are above
the dynamical range sampled by GAMA. Some of the latter measurements are upper
limits, and are indicated by down-pointing arrows. Figure 9: Specific SFR as a
function of BH mass for galaxies at $z=0$ with stellar masses $>10^{10}\,\rm
M_{\odot}$. The results for Shark v1.1 (L18) and v2.0 are shown as a black
line with the grey shaded region, and red line with the faded red shaded
region, respectively. Lines show medians and the shaded regions show the
$16^{\rm th}-84^{\rm th}$ percentile ranges for bins with $\geq 10$ objects.
Observational estimates from Terrazas et al. (2017) are shown as symbols.
Figure 10: The median SFR as a function of lookback time (LBT) of galaxies at
$z=0$ in different bins of stellar mass (of width $0.25$ dex) for Shark v2.0
(solid lines) and v1.1 (dotted lines). The mean stellar mass of each bin are
labelled next to the corresponding lines.

Fig. 8 shows the SFR-stellar mass plane at $z=0$ in Shark v2.0 and L18.
Observational measurements of the median SFR-stellar mass relation of SDSS
galaxies from Brinchmann et al. (2004), the main sequence in GAMA from Davies
et al. (2016) and the distribution of all GAMA galaxies at $z\leq 0.06$ and
their median SFR-stellar mass plane distribution from the catalogue of
Robotham et al. (2020); Bellstedt et al. (2020b), are also shown. For the
latter, galaxies both on and off the main sequence (MS) are included.

The median SFR at fixed stellar mass is very similar between Shark v2.0 and
L18 in galaxies with $M_{\star}<10^{10}\,\rm M_{\odot}$. Above that stellar
mass, the models diverge in a way that galaxies in Shark v2.0 appear more
quenched than those in L18 by roughly an order of magnitude, at least for
galaxies with stellar masses $10^{10.5}-10^{11.5}\,\rm M_{\odot}$. The most
massive galaxies, $M_{\star}>10^{12}\,\rm M_{\odot}$, in both versions of
Shark have similar SFRs, about an order of magnitude below the MS.

Compared with observations, we see that the MS in both Shark versions agree
well with the measurements of Davies et al. (2016). There is evidence though
that galaxies in GAMA as analysed in Bellstedt et al. (2023) start to, on
average, deviate from the MS at $M_{\star}\gtrsim 10^{9.5}\,\rm M_{\odot}$.
This in Shark v2.0 and L18 happens at higher stellar masses, $M_{\star}\approx
10^{10}\,\rm M_{\odot}$. The deviations from MS in the early SDSS analysis of
Brinchmann et al. (2004) happens at higher stellar masses, closer to what we
get in Shark. Because GAMA’s survey area is too small to have a representative
sample of very massive galaxies, we also show individual measurements of the
SFR and stellar mass of a sample of massive galaxies from Terrazas et al.
(2017). Note that some of these measurements are upper limits.

The quenching of massive galaxies is linked to AGN feedback. It is thus
natural to explore the connection between quenching and BH mass. Fig. 9 shows
the specific SFR ($\equiv\rm SFR/M_{\star}$; sSFR) as a function of the
central BH mass of galaxies with $M_{\star}>10^{10}\rm\,M_{\odot}$ at $z=0$
for Shark v2.0 and L18, compared with the observations of Terrazas et al.
(2017). This plane clearly shows a BH mass scale at which quenching happens in
Shark. This transition mass for Shark v2.0 is $M_{\rm BH}\approx 10^{7.5}\,\rm
M_{\odot}$ and for L18 is $M_{\rm BH}\approx 10^{7.25}\,\rm M_{\odot}$. What
is interesting is that above that mass threshold, Shark v2.0 galaxies reach
much lower sSFRs (i.e. are more quenched) than L18 galaxies, but also the
scatter of sSFR at fixed BH mass is much larger in Shark v2.0 than L18. The
relation obtained by Shark v2.0 resembles much better what the observations of
Terrazas et al. (2017) suggest. Terrazas et al. (2020) used this plane to
diagnose the way AGN feedback acts in the Illustris-TNG simulations, and found
that the transition region in Illustris-TNG, which happens at $M_{\rm
BH}\approx 10^{8}\,\rm M_{\odot}$, was much sharper than observations suggest,
and even bimodal, with galaxies below the transition mass being star-forming,
and right above, being quenched. Shark v2.0 displays a much smoother
transition that agrees better with observations. We came back to the
sSFR-$M_{\rm BH}$ plane in $\S$ 5.2.

Lastly, we analyse the SFR histories (SFHs) of galaxies at $z=0$ in different
bins of stellar mass for Shark v2.0 and L18 in Fig. 10. Galaxies in both Shark
v2.0 and L18 have on average raising SFHs in the stellar mass range
$10^{9}-10^{9.8}\,\rm M_{\odot}$. Galaxies in the next stellar mass bin,
$10^{10}-10^{10.25}\,\rm M_{\odot}$, in both models have a median SFH that is
close to constant over the last $6$ Gyrs. At higher stellar masses there is a
clear peak at lookback times $>9$ Gyrs, followed by a decline in SFR. In L18,
galaxies with $M_{\star}\gtrsim 10^{10.5}\,\rm M_{\odot}$ have relatively
self-similar SFHs, scaled up or down in the overall normalisation. In Shark
v2.0, the SFHs of massive galaxies present a much larger variations at fixed
stellar mass, and on average end up more quenched towards $z=0$. In Shark
v1.1, once quenching kicks in for galaxies with $M_{\star}>10^{10.5}\,\rm
M_{\odot}$, we see that the SFR continues to increase with increasing stellar
mass, while in Shark v2.0, galaxies with $10^{10.5}\,\rm
M_{\odot}<M_{\star}<10^{10.75}\,\rm M_{\odot}$ have higher SFRs at lookback
time $<4$ Gyr than galaxies with $10^{10.75}\,\rm
M_{\odot}<M_{\star}<10^{11.5}\,\rm M_{\odot}$, on average. This agrees
qualitatively with the “downsizing” signal observed in local Universe galaxies
(e.g. Thomas et al. 2010; Bellstedt et al. 2020b), and shows that this new
version of Shark overcomes one of the shortcomings of L18 identified in Bravo
et al. (2022, 2023).

###  5.2 The onset of galaxy quenching

Figure 11: The specific SFR vs BH mass for galaxies with stellar masses $\geq
10^{9}\,\rm M_{\odot}$ from $z=0$ to $z=6$, as labelled, for Shark v1.0
(right) and v2.0 (left). The applied stellar mass limit is to avoid the
emergence of quenched satellite galaxies producing a visible decrease in
specific SFR as the focus here is on AGN feedback. Lines show the medians, and
the shaded regions show the $16^{\rm th}-84^{\rm th}$ percentile ranges. The
latter are only shown for $z=0$, $z=3$ and $z=6$ to avoid crowding the figure.
The vertical dotted line shows for reference a BH mass of
$10^{7.5}\rm\,M_{\odot}$.

To study the onset of quenching in galaxies, we start by exploring the sSFR-BH
mass plane across cosmic time in Fig. 11 for Shark v2.0 and L18. We show how
the median relation for galaxies with $M_{\star}>10^{9}\,\rm M_{\odot}$
evolves from $z=0$ to $z=6$ in both models.

In L18, the BH mass at which galaxies transition from being primarily star-
forming to displaying signs of quenching (i.e. a decrease in sSFR) increases
with increasing redshift. By $z=3$, this happens at $M_{\rm BH}\approx
10^{8.3}\,\rm M_{\odot}$, $\approx 1.2$ dex higher than the transition mass at
$z=0$ for that model. At $z=4$ there are only weak signs of quenching
happening in galaxies with $M_{\rm BH}\gtrsim 10^{9}\,\rm M_{\odot}$. There is
a stark contrast with Shark v2.0, which predicts that the BH transition mass
is relatively redshift independent at $M_{\rm BH}\approx 10^{7.5}\,\rm
M_{\odot}$. The reason for this difference comes down to how radio-mode
feedback was implemented in L18 vs Shark v2.0. As discussed in $\S$ 3.3, the
heating power of AGN feedback in the radio mode of Croton et al. (2006) was
calculated using only the accretion rate coming from the hot-halo mode instead
of the total accretion rate of the BH. The hot-halo mode contributes very
little to the total BH accretion rate at high-z, hence we need to move to very
high BH masses to start to see a larger contribution from that mode, which can
then lead to quenching. In Shark v2.0 we instead are agnostic to the accretion
channel which makes more physical sense, and simply calculate the total jet
power that can be used to offset gas cooling. This leads to an approximately
constant BH mass above which enough jet power is produced to lead to quenching
(or at least to significant deviations from the MS). Combining robust BH mass
and sSFR measurements in observations across cosmic time would greatly help to
constrain the predictions here.

One of the key motivations to study the quenching of massive galaxies comes
from the recent discoveries of a sizeable population of massive-quiescent
galaxies in the early Universe (e.g. Schreiber et al. 2018; Carnall et al.
2020; Weaver et al. 2022; Gould et al. 2023; Nanayakkara et al. 2022; Carnall
et al. 2023; Valentino et al. 2023; Long et al. 2023). These observations have
revealed that these galaxies exist in number densities $\gtrsim 10^{-5}\,\rm
Mpc^{-3}$. We hence investigate what L18 and Shark v2.0 predict in Fig. 12 for
the number density of massive-quiescent galaxies at $0\leq z\leq 5$. We do
this by using two definitions of sSFR, $<10^{-10}\,\rm yr^{-1}$ (as adopted by
Long et al. 2022) and $<10^{-11}\,\rm yr^{-1}$ (a more conservative cut).
These two thresholds represent deviations from the MS at $z=3$ of $\approx
1.3-1.5$ and $2.3-2.5$ dex, respectively. We also investigate two stellar mass
selections, $\gtrsim 10^{10}\,\rm M_{\odot}$ (top panel in Fig. 12) and
$\gtrsim 10^{10.5}\,\rm M_{\odot}$ (bottom panel in Fig. 12), as they
represent typical mass thresholds being adopted in observations.

Figure 12: Number density of passive massive galaxies selected based on two
stellar mass thresholds, $M_{\star}>10^{10}\,\rm M_{\odot}$ (top) and
$M_{\star}>10^{10.5}\,\rm M_{\odot}$ (bottom), as labelled. To select passive
galaxies we use two different thresholds in specific SFR as labelled, which
approximate typical values adopted in the literature. For Shark v2.0 we show
the intrinsic predictions (i.e. taking stellar masses and SFRs without
assuming any errors; thick solid lines), and the prediction after convolving
stellar masses and SFRs with errors that are Gaussian-distributed with a width
$0.3$ dex but no bias (i.e. Gaussian centred at $0$; dashed lines). Shark v1.1
intrinsic predictions are also shown in both panels (thin solid lines, as
labelled in the bottom panel). Observational estimates are also shown: in grey
symbols (pre-JWST results) are from Straatman et al. (2014); Schreiber et al.
(2018); Merlin et al. (2019); Girelli et al. (2019); Carnall et al. (2020);
Weaver et al. (2022); Gould et al. (2023), and green symbols (JWST results)
are from Nanayakkara et al. (2022); Carnall et al. (2023); Valentino et al.
(2023), as labelled at the top panel.  Figure 13: Galaxy SMF in Shark v2.0
(solid red line) and the default model in L18 (solid black line), for galaxies
classified at passive with a simple threshold of $\rm sSFR<10^{-10.75}\,\rm
yr^{-1}$ from $z=0.5$ to $z=5$, as labelled. We also show for Shark v2.0 the
SMF of passive galaxies after applying a random error to the stellar masses
and SFRs of $0.3$ dex (dot-dashed line) and then splitting those galaxies
between centrals (orange) and satellites (green). We show the same for central
and satellite galaxies in L18 as dashed and dotted black lines, respectively.
Observations from Weaver et al. (2022) from the COSMOS survey are also shown.
Vertical dotted lines mark the approximate stellar mass above which
observational estimates of the SMF are reliable according to Weaver et al.
(2022). Figure 14: As Fig. 13 but for centrals and satellites in Shark v2.0
selected to have $\rm sSFR<10^{-10.75}\,\rm yr^{-1}$ (dashed lines), or to
comply with the colour-colour selection of Eq. (55) (solid lines). The latter
is what was used by Weaver et al. (2022) to select passive galaxies and
construct the SMF.

Fig. 12 shows two important differences between L18 and Shark v2.0: (1)
overall the number density of massive-quiescent galaxies in Shark v2.0 is
$\gtrsim 1$ dex higher than L18, with differences being largest for galaxies
with $\rm sSFR<10^{-11}\,yr^{-1}$, $\approx 2$ dex; (2) there is little
difference in number density between selecting galaxies with $\rm
sSFR<10^{-11}\,yr^{-1}$ or $\rm sSFR<10^{-10}\,yr^{-1}$ in Shark v2.0
($\lesssim 0.2$ dex), compared to the large differences seen in L18 ($\gtrsim
0.5$ dex). This tells us that quenched galaxies in Shark v2.0 have on average
sSFR $\gtrsim 1$ dex smaller than L18. We find that the large difference
between Shark v2.0 and L18 is due to the new radio-mode AGN feedback model
rather than the inclusion of QSO-mode AGN feedback. When we turn off the
latter we only see a slightly lower number density of massive-quiescent
galaxies at $z\geq 4$ of $0.2$ dex, and almost no differences at $z<4$.

Overall the number density of massive-quiescent galaxies predicted by Shark
v2.0 are in good agreement with current observational constraints within the
uncertainties. Most of the scatter between current observational estimates
comes from the different employed methods to define “quenched” (e.g. colour-
colour selection vs post-starburst features in the spectrum; see Valentino et
al. 2023 for a discussion). The comparison between observations with
simulations is further complicated by the unconstrained systematic and random
uncertainties in the estimation of stellar mass and SFRs in these galaxies. To
get a sense of what the impact of this could be, we convolve the stellar
masses and SFRs in Shark galaxies with a random Gaussian distribution centred
around $0$ and with a width of $0.3$ dex. We select then massive-quiescent
galaxies based on sSFRs and stellar masses after including errors. This effect
is shown as dashed lines in Fig. 12. The consequence of adding errors can be
quite large, especially for the rarest objects (bottom panel in Fig. 12 shows
$\approx 0.8$ dex difference between having errors or not in the population
with $\rm sSFR<10^{-10}\,yr^{-1}$). We should note that errors of the order of
$0.3$ dex for stellar masses and SFRs are likely lower limits, as even at
$z=0$ with better data quality and multi-wavelength coverage, stellar mass
errors are typically $0.2$ dex (Robotham et al., 2020). A full understanding
of the level of agreement/tension between observations and simulations likely
requires full forward modelling to select galaxies in the same colour-colour
space or with the same spectral features as in observations. We leave that for
future work.

We now investigate the evolution of the SMF of passive galaxies in Fig. 13
from $0.5\leq z\leq 5$. We select passive galaxies as those with a $\rm
sSFR<10^{-10.75}\,\rm yr^{-1}$ following Shuntov et al. (in preparation). We
discuss the effect of the passive selection criteria when discussing Fig. 14.
Comparing Shark v2.0 with L18 first (solid red and dashed black lines,
respectively), we find that at all redshifts Shark v2.0 produces more passive
galaxies across the whole mass range than L18. The differences become larger
with increasing redshift, from $\approx 0.3$ dex at $z=0.5$ to $1.5$ dex at
$z=5$. The difference in number density between the two models is largest for
central galaxies (orange solid vs dotted black lines), with L18 producing a
number density of central, passive galaxies $\lesssim 10^{-5.5}\,\rm Mpc^{-3}$
at $M_{\star}\approx 10^{10.5}\,\rm M_{\odot}$, while Shark v2.0 predicts
number densities of $10^{-4}\rm\,M_{\odot}$. At $z\geq 3$, L18 predicts
virtually no quenched, central galaxies, with number densities $<10^{-6}\,\rm
M_{\odot}$, in contrast with Shark v2.0 which predicts peak number densities
of $\approx 10^{-4.5}\,\rm Mpc^{-3}$ and $\approx 10^{-5.5}\,\rm Mpc^{-3}$ at
$z=3$ and $z=4$, respectively. At $z\leq 3$ a population of passive, low-mass
centrals emerges in both Shark v2.0 and L18 with masses $M_{\star}\lesssim
10^{9}\,\rm M_{\odot}$. These galaxies corresponds to those inhabiting very
low-mass halos that suffer from photo-ionisation from the UV background, whose
efficiency depends on a halo’s circular velocity and redshift (Sobacchi &
Mesinger, 2013), as described in $\S$ 4.4.9 in L18.

The difference between the SMF of passive, satellite galaxies between Shark
v2.0 and L18 comes in part from the new dynamical friction timescale model we
adopted ($\S$ 3.2) in v2.0 producing longer dynamical friction timescales
leading to a longer survival of satellite galaxies than the Lacey & Cole
(1993) model adopted in L18. Satellite galaxies make up most of the population
of passive galaxies at $M_{\star}\lesssim 10^{10.2}\,\rm M_{\odot}$ at $z=0.5$
and at increasingly higher stellar masses with increasing redshift. For
example, at $z\geq 4$, passive, satellite galaxies dominate the number density
across the whole mass range, showing that environment (i.e. RPS and tidal
stripping) can effectively quench galaxies at very early cosmic times.

Fig. 13 shows the observational results of Weaver et al. (2022), in which a
colour-colour selection was used to classify galaxies as passive or star-
forming using the near-UV (NUV), r-band (r), and J-band (J) magnitudes, with
passive galaxies being those with

$\rm(NUV-r)>3\,(r-J)\,+\,1;\,(NUV-r)>3.1,$ (55)

and star-forming galaxies are those that do not comply with the selection
above. Compared with the observations of Weaver et al. (2022), Shark v2.0
performs overall better than L18, reproducing well the population of passive,
central galaxies which likely dominate the high-mass end in the observations.
At $z\geq 3$, Shark v2.0 struggles to reproduce the high-mass end of the
passive SMF, but this is in part due to systematic uncertainties in the
observations. For example, Shuntov et al. (in preparation) using COSMOS-Web
(Casey et al., 2022) show that the high-mass end they recover at those high-
redshift for passive galaxies is about $0.5$ dex lower than that reported by
Weaver et al. (2022) and closer to the Shark v2.0 predictions. Another
important limitation to mention in observations is that the number densities
are unconstrained for galaxies that are below the vertical lines shown in Fig.
13, and hence the onset of environmental quenching in the very early Universe
as predicted by Shark cannot be clearly studied with current observations of
the SMF.

We assess the effect the selection of passive galaxies has on the resulting
SMF, we compare the SMF of galaxies with a sSFR$<10^{-10.75}\,\rm yr^{-1}$ and
that comply with the colour-colour selection of Eq. (55). We compute galaxy
spectral energy distributions following the method described in Lagos et al.
(2019) (referred to as “EAGLE-$\tau$ RR14” in that paper - see Section 6 in
the supplementary material for a short description). The supplementary
material shows that Shark v2.0 produces galaxy luminosities that reproduce
reasonably well the observed luminosity functions of the local Universe from
the NUV to the FIR, the K-band luminosity function evolution up to $z=3$, and
the far-UV luminosity function up to $z=10$.

The SMF of the passive galaxy populations selected by the methods above are
shown in Fig. 14. Focusing first on central galaxies, we see that at $z\leq 1$
Eq. (55) is effective in selecting galaxies of low sSFR and stellar masses
$\gtrsim 10^{10}\,\rm M_{\odot}$, however, it tends to overestimate the number
of low sSFR central galaxies at lower stellar masses by up to $0.8$ dex. The
problem gets worse at increasing redshift, and by $z\geq 3$ most of the
centrals classified as passive by their colour have sSFR$>10^{-10.75}\,\rm
yr^{-1}$. By $z=5$ there are $>1000$ times more galaxies with colours
consistent with being passive but with sSFR$>10^{-10.75}\,\rm yr^{-1}$.
Interestingly, for satellite galaxies we see the opposite. At $z\leq 1$
galaxies with sSFR$<10^{-10.75}\,\rm yr^{-1}$ comply with the colour selection
of Eq. (55), but at higher redshifts, the colour selection underestimates the
number of satellites with low sSFRs. By $z=5$ this underestimation is of
$\approx 1$ dex. To understand where these disparate effects on centrals and
satellite galaxies come from, we studied the median properties of
centrals/satellites with $M_{\star}\geq 10^{10}\,\rm M_{\odot}$. We found that
centrals primarily correspond to star-bursting, dusty galaxies (i.e. most of
their SFR is associated to a burst with super-solar metallicities) at $z>2$.
Satellites at $z\geq 2$ on the other hand, have metallicities that are $0.7-2$
dex below centrals and have their star formation taking place in the galaxy
disk. Hence, Eq. (55) appears to be effective in selecting passive galaxies
when they have gas metallicities close to solar metallicity and they are not
dusty star-forming galaxies (see Lagos et al. 2020 for a detailed analysis of
the contamination of dusty star-forming galaxies in the selection of passive
galaxies based on colour-colour diagrams).

We overall see that the SMF of galaxies selected to be “passive” following the
colour-colour selection of Eq. (55) is in much better agreement with the
observational estimates of Weaver et al. (2022), coming from the same colour-
colour selection (see for example high-mass end at $z=2$ and $z=3$). In the
future we plan to investigate other more modern methods to select passive
galaxies based on colours using independent magnitudes and Bayesian methods to
assign probabilities of being passive (Gould et al., 2023; Long et al., 2023).

Finally, we present predictions for the stellar-to-halo mass relation of
passive and star-forming central galaxies in Shark v2.0 in Fig. 15. We again
test two definitions of passive galaxies as described above (i.e. based on
sSFR and colour), while star-forming galaxies are defined based on their sSFR
distance to the MS or by not complying with the colour-colour selection of Eq.
(55). These two selections are referred to as “sSFR” and “colour” in Fig. 15.
We measure the MS in Shark v2.0 by doing a linear fit to the relation between
$\rm log_{10}(sSFR)$ and $\rm log_{10}(M_{\star})$ for central galaxies with
stellar masses in the range $7\times 10^{8}-10^{10}\,\rm M_{\odot}$ at each of
the shown redshifts. These mass limits are chosen conservatively to avoid
resolution effects at the lower end, and the regime where AGN feedback becomes
efficient at the higher end. Star-forming galaxies, in the case of the “sSFR”
selection, are those with a distance to the main sequence $\geq-0.3$ dex.

Figure 15: Stellar vs halo mass relation for central galaxies in Shark v2.0
selected to be passive by either assuming $\rm sSFR<10^{-10.75}\,\rm yr^{-1}$
(salmon dashed line), or by selecting them based on their NUV-optical colour
(see text for details; red solid line); and star-forming, by either selecting
them to have sSFRs that are $>-0.3$ dex from the MS (light blue dashed line)
or based on their NUV-optical colour (blue solid line). Lines with shaded
regions show the medians and $16^{\rm th}-84^{\rm th}$ percentile ranges,
respectively. Only bins with $\geq 10$ galaxies are shown. This is presented
for $z=1,\,2,\,3$, as labelled. Overall, passive galaxies display a flatter
stellar-halo mass relation than star-forming galaxies, with the difference
being larger when galaxies are selected based on their sSFR rather than
colour.

Fig. 15 shows that passive galaxies display a weaker dependence of their
stellar mass on halo mass at $M_{\rm halo}\gtrsim 10^{11.6}\rm\,M_{\odot}$
than star-forming galaxies. The difference is larger when we use the sSFR to
select galaxies than with the colour-colour selection. Interestingly, the
latter misses most of the central galaxies with low sSFRs that are hosted by
halos with $M_{\rm halo}<10^{11.6}\,\rm M_{\odot}$ at $z=1$ and $z=2$. At
$z=3$ the colour-colour selection picks centrals in those low-mass halos that
are not passive according to their sSFR. The overall different stellar-halo
mass relation of passive and star-forming galaxies should translate into
different clustering signals for these two populations. The flatness of the
stellar-halo mass relations of passive central galaxies leads to their stellar
masses being lower than that of star-forming galaxies at fixed halo mass in
halos with masses $\gtrsim 10^{12.5}\,\rm M_{\odot}$, on average. At halo
masses $10^{11.6}\lesssim M_{\rm halo}/\rm M_{\odot}\lesssim 10^{12.5}$, the
flatness of the stellar-halo mass relation of passive galaxies leads to them
having higher stellar masses than star-forming galaxies. This is very clear at
$z=1$ and apparent at $z=2$. At $z=3$ we have too few high-mass halos to
establish the continuation of such a difference. We will revisit this using
larger cosmological volume simulations in the future.

At $M_{\rm halo}\lesssim 10^{11.6}\,\rm M_{\odot}$, the stellar-halo mass
relation of passive galaxies is very different between the two methods
employed to select them, with the colour-colour selection leading to a
relation that is similar to that of star-forming galaxies (especially clear at
$z=1$ and $z=2$).

One of the only observational measurements we can compare Fig. 15 with were
presented by Cowley et al. (2019), who derived a stellar-to-halo mass relation
of passive and star-forming galaxies from the clustering of these populations
at $z\approx 2-3$. Note that Cowley et al. (2019) use colour selections
similar to Eq. (55) to isolate passive galaxies. Cowley et al. (2019) found
that at $M_{\rm halo}\gtrsim 10^{12.7}\,\rm M_{\odot}$ passive and star-
forming galaxies followed similar relations, but below that mass, passive
galaxies were more massive than star-forming galaxies at fixed halo mass by
$\approx 0.5$ dex. This difference is larger than what we find at $z=2,\,3$
and $M_{\rm halo}\lesssim 10^{12.5}\,\rm M_{\odot}$ in Shark v2.0 (which is
closer to $0.3$ dex), but overall consistent within the scatter. At higher
halo masses, the colour selection indeed leads to stellar-halo mass relations
between star-forming and passive galaxies in Shark v2.0 that are similar (and
even indistinguishable at $z=3$), in agreement with the conclusions of Cowley
et al. (2019). However, if galaxies were selected by sSFR we predict that
different stellar-halo mass relations of passive and star-forming galaxies
should be seen at these high halo masses at least up to $z=2$. Magliocchetti
et al. (2023), also using clustering measurements, presented constraints on
the host halo mass of massive-quiescent galaxies (those with $M_{\star}\approx
10^{10}-10^{11}\,\rm M_{\odot}$) and show that the host halo masses of these
galaxies are likely $\gtrsim 10^{12}\,\rm M_{\odot}$ across a wide redshift
range, $0\leq z\leq 5$. Although we broadly see something similar in Shark
v2.0, a careful comparison requires selecting galaxies using the same colour-
colour selection employed in Magliocchetti et al. (2023) and adding
photometric redshift uncertainties, which we leave for future work.

The predictions of Fig. 15 will be testable in the near future thanks to JWST
programs focused on obtaining large statistical samples of massive-quiescent
galaxies across cosmic time.

##  6 Discussion and conclusions

We have introduced a new version of the Shark SAM (v2.0) after a number of
improvements to the physics included. These changes comprise: (i) a model
describing the exchange of angular momentum between the interstellar medium
and stars that results in different angular momenta for the atomic and
molecular ISM and stars in galaxies; (ii) an updated dynamical friction
timescale of satellite galaxies; (iii) a new AGN feedback model which includes
two modes, a QSO and a radio mode, with the radio mode directly tied to the
jet power production and the QSO mode consisting of a radiation pressure-
driven outflow; (iv) a model that tracks the development of BH spins, which
together with the mass and accretion rate, are used to define different BH
accretion disk states; (v) a model for the gradual ram-pressure stripping of
the hot and cold gas of satellite galaxies; (vi) a model for the tidal
stripping of gas and stars of satellite galaxies; (vii) a method for automatic
parameter exploration of the model using particle swarm optimisation. The
model parameters were chosen to fit the $z<0.1$ SMF of (Li & White, 2009)
only. No high redshift constraints were used.

We showed that Shark v2.0 provides predictions that agree better with
observations than Shark v1.1 (L18). Those include: the evolution of the SMF up
to $z=7$ (Fig. 1), the halo-mass conditional baryon mass function at $z=0$ and
the contribution from satellite and central galaxies (Fig. 2); the stellar
size-mass relation (Fig. 5) and specific angular momentum-mass relation of
different baryon components of galaxies (Fig 6) at $z=0$; the BH-stellar mass
relation of late- and early-type galaxies at $z=0$ (Fig. 7). We show that
these improvements in large part relate to the new physics included in Shark.
Specifically the conditional baryon mass function of satellite galaxies agrees
better with observations in part due to the new dynamical friction timescale;
and the stellar size-mass and specific angular momentum-mass relations
improved significantly thanks to the new angular momentum treatment of galaxy
components.

We presented a detailed analysis of galaxy quenching in Shark v2.0 and
improvements over v1.1 (L18). These included:

  * •

Massive galaxies at $z=0$, $M_{\star}\gtrsim 10^{10.5}\,\rm M_{\odot}$ have
$\approx 1$ dex lower SFRs in Shark v2.0 than v1.1. This allows the model to
better reproduce the SFRs of massive galaxies observed (Fig. 8) and the sSFR-
BH mass relation of galaxies with $M_{\star}>10^{10}\,\rm M_{\odot}$ (Fig. 9)
at $z=0$. The latter shows that Shark v2.0 also produces a much larger sSFR
scatter at fixed BH mass for BH masses $>10^{8}\,\rm M_{\odot}$ than v1.1, in
much better agreement with the observations of Terrazas et al. (2017).

  * •

Shark v2.0 predicts that the transition of galaxies from being star-forming
(i.e. main sequence) to displaying a clear decrease in sSFR relative to the
main sequence happens at roughly the same BH mass at all redshifts ($\approx
10^{7.5}\,\rm M_{\odot}$), while v1.1 requires an increasingly massive BH to
transition to quenched with increasing redshift (Fig. 11). By $z=4$, in Shark
v1.1 only galaxies with a BH mass $\gtrsim 10^{9}\,\rm M_{\odot}$ show signs
of a decreased sSFR.

  * •

Shark v2.0 produces $\approx 1$ dex higher number density of massive-quiescent
galaxies at $z\gtrsim 2$ than v1.1, and those that are classified as being
quiescent display lower sSFRs in Shark v2.0 than v1.1. Our new results agree
well with current observational constraints on the number density of massive-
quiescent galaxies coming from the JWST (Carnall et al., 2023; Nanayakkara et
al., 2022; Valentino et al., 2023; Long et al., 2023) (see Fig. 12). We
highlight that the overall abundance of massive galaxies in Shark v2.0 is
similar to that of v1.1 (L18), but the key difference is in the fraction of
those that are quenched.

  * •

We analyse the SMF of passive galaxies from $z=0.5$ to $z=5$ and show that at
$z\geq 2$ Shark v2.0 produces $\gtrsim 100$ times more central-passive
galaxies than v1.1, and this difference increases with increasing redshift.
Similarly, Shark v2.0 predicts $\approx 10-80$ times more passive satellite
galaxies at $z\geq 3$ than v1.1. These differences tend to disappear towards
$z=0$ with both models predicting similar SMF of passive galaxies (Fig. 13).
Our new model is in better agreement with current observational constraints of
the passive SMF (especially when we select galaxies using the same colour-
colour criterion employed in observations), though there are likely still too
few passive galaxies with masses $>10^{11}\,\rm M_{\odot}$ at $z\geq 4$ in
Shark v2.0.

  * •

We present predictions for the stellar-halo mass relation of star-forming and
passive galaxies at $z\geq 1$ and find clear differences between the two
populations. Passive galaxies tend to display a flat stellar-halo mass
relation at $M_{\rm halo}\gtrsim 10^{11}\,\rm M_{\odot}$, so that at $M_{\rm
halo}\approx 10^{11}-10^{12.5}\,\rm M_{\odot}$ they are more massive than
their star-forming counterparts, while the opposite happens at $M_{\rm
halo}\gtrsim 10^{12.5}\,\rm M_{\odot}$ (Fig. 15). The exact magnitude of this
effect depends on the criteria used to select passive and star-forming
galaxies.

The differences listed above between Shark v1.1 and v2.0 can almost
exclusively be associated with the implementation of the new AGN model
presented in $\S$ 3.3, although other processes, such as the new dynamical
friction timescale, also play a (secondary) role. This work thus demonstrates
the power of SAMs in quickly assessing the possible physical mechanisms behind
current tensions with observations.

In accompanying papers we present a detailed analysis of the BH and AGN
population across cosmic time (Bravo et al. in preparation); and a thorough
analysis of the parameter space of Shark v1.1 and v2.0 (Proctor et al. in
preparation). In the future, we plan to explore different definitions of what
makes a galaxy passive and the impact that has on the number density, SMF
evolution and clustering of passive galaxies; and assess the derived
properties of massive-quiescent galaxies to understand the effect of
systematic uncertainties.

## Data Availability

The surfs halo and subhalo catalogue and corresponding merger trees used in
this work can be accessed from https://tinyurl.com/y6ql46d4. Shark is a public