Star Formation in Simulated Galaxies: Understanding the Transition to Quiescence at $3\times10^{10}$ M$_\odot$

TL;DR

This study uses hydrodynamical cosmological simulations to explore how star formation in galaxies is regulated by gas availability, black hole activity, and environment, revealing a critical black hole mass for quenching star formation.

Contribution

The paper introduces a model linking black hole mass to star formation quenching, matching observed galaxy transition mass and reproducing observed star formation rates.

Findings

Gas mass controls star formation in low-mass galaxies.

Black holes >10^7.5 M_sun quench star formation in high-mass galaxies.

Model aligns with observed star formation rates from GAMA and ALFALFA surveys.

Abstract

Star formation in galaxies relies on the availability of cold, dense gas, which, in turn, relies on factors internal and external to the galaxies. In order to provide a simple model for how star formation is regulated by various physical processes in galaxies, we analyse data at redshift $z=0$ from a hydrodynamical cosmological simulation that includes prescriptions for star formation and stellar evolution, active galactic nuclei (AGN), and their associated feedback processes. This model can determine the star formation rate (SFR) as a function of galaxy stellar mass, gas mass, black hole mass, and environment. We find that gas mass is the most important quantity controlling star formation in low-mass galaxies, and star-forming galaxies in dense environments have higher SFR than their counterparts in the field. In high-mass galaxies, we find that black holes more massive than $\sim10^{7.5}$ M$_\odot$ can be triggered to quench star formation in their host; this mass scale is emergent in our simulations. Furthermore, this black hole mass corresponds to a galaxy bulge mass $\sim2\times10^{10}$ M$_\odot$, consistent with the mass at which galaxies start to become dominated by early types ($\sim3\times10^{10}$ M$_\odot$, as previously shown in observations by Kauffmann et al.). Finally, we demonstrate that our model can reproduce well the SFR measured from observations of galaxies in the GAMA and ALFALFA surveys.

Tables

Table 1. Table 1 : Fitted parameters of the models presented in Section 3 . The parameters denote: C 𝐶 C – normalisation; α 𝛼 \alpha – coefficient of log ⁡ M ∗ subscript 𝑀 \log M_{*} ; β 𝛽 \beta – coefficient of log ⁡ s 5 subscript 𝑠 5 \log s_{5} , γ 𝛾 \gamma – coefficient of log ⁡ M gas subscript 𝑀 gas \log M_{\rm gas} ; k 1 subscript 𝑘 1 k_{1} – coefficient of log ⁡ M BH subscript 𝑀 BH \log M_{\rm BH} (low mass); k 2 subscript 𝑘 2 k_{2} – coefficient of log ⁡ M BH subscript 𝑀 BH \log M_{\rm BH} (high mass); log ⁡ M b subscript 𝑀 b \log M_{\rm b} – break mass separating low- and high-mass BHs; Δ logM subscript Δ logM \Delta_{\rm logM} – width of transition from low- to high-mass BHs. Three values for each parameter are given for each model: the first is the best-fit value from the simulated data; the second is the mean and standard deviation from the full bootstrap parameter distributions; and the third is the modal value of the distributions with asymmetric errors denoting the 16 th superscript 16 th 16^{\rm th} and 84 th superscript 84 th 84^{\rm th} percentiles. The distribution of Δ logM subscript Δ logM \Delta_{\rm logM} in model 3 is bimodal, so asymmetric error bars are unavailable; this value is marked † .

$x$	Best Fit	$\overline{x}\pm {\sigma }_x$	Mode ${}_{-}{}^{}{}_{}{}^{+}$ distribution width
Model 1 (equation (6))
$C$	$-2.50$	$-2.59\pm 0.18$	$-{2.58}_{-0.18}^{+0.18}$
$\alpha$	$0.21$	$0.21\pm 0.03$	${0.21}_{-0.03}^{+0.03}$
$\beta$	$-0.17$	$-0.15\pm 0.04$	$-{0.16}_{-0.03}^{+0.04}$
$\gamma$	$1.03$	$1.04\pm 0.05$	${1.04}_{-0.05}^{+0.05}$
$k_1$	$0.00$	$-0.02\pm 0.03$	$-{0.01}_{-0.04}^{+0.03}$
$k_2$	$-1.11$	$-1.99\pm 3.13$	$-{1.12}_{-0.82}^{+0.52}$
$\log M_b$	$7.31$	$7.79\pm 0.64$	${7.29}_{-0.15}^{+0.46}$
Model 2 (equation (7))
$C$	$-2.72$	$-2.72\pm 0.15$	$-{2.72}_{-0.14}^{+0.15}$
$\alpha$	$0.21$	$0.21\pm 0.03$	${0.21}_{-0.03}^{+0.03}$
$\beta$	$-0.12$	$-0.12\pm 0.04$	$-{0.12}_{-0.03}^{+0.04}$
$\gamma$	$1.05$	$1.05\pm 0.05$	${1.04}_{-0.05}^{+0.06}$
$k_1$	$-0.05$	$-0.05\pm 0.02$	$-{0.04}_{-0.03}^{+0.02}$
Model 3 (equation (8))
$C$	$-3.12$	$-3.08\pm 0.19$	$-{3.08}_{-0.18}^{+0.19}$
$\alpha$	$0.21$	$0.22\pm 0.03$	${0.21}_{-0.03}^{+0.03}$
$\beta$	$-0.16$	$-0.16\pm 0.03$	$-{0.16}_{-0.03}^{+0.03}$
$\gamma$	$1.05$	$1.05\pm 0.05$	${1.06}_{-0.05}^{+0.05}$
$k_1$	$0.04$	$0.06\pm 0.05$	${0.08}_{-0.08}^{+0.03}$
$k_2$	$-1.19$	$-1.25\pm 0.85$	$-{0.82}_{-0.53}^{+0.31}$
$\log M_b$	$7.68$	$7.66\pm 0.15$	${7.58}_{-0.07}^{+0.14}$
${\Delta }_{\mathrm{logM}}$	$1.22$	$1.55\pm 0.76$	${2.06}^{†}$

Equations

$${\dot{M}}_{\mathrm{acc}}=\frac{4\pi G^2{M}_{\mathrm{BH}}^2\rho }{{\left(c_s^2+v^2\right)}^{3/2}},$$

$${\dot{M}}_{\mathrm{Edd}}=\frac{4\pi G{M}_{\mathrm{BH}}{m}_{\mathrm{p}}}{{ϵ}_{\mathrm{r}}{\sigma }_{\mathrm{T}}c}\equiv \frac{M_{\mathrm{BH}}}{t_{\mathrm{Sal}}},$$

$$\log {\mathrm{SFR}}_{\mathrm{fit}}=-9.32+0.92\log M_{\ast },$$

$$\mathrm{SFR}=A\times M_{\ast }^{\alpha }\times s_5^{\beta }\times M_{\mathrm{gas}}^{\gamma }\times f(M_{\mathrm{BH}}),$$

$$\begin{array}{rl}\log \mathrm{SFR}=C& +\alpha \log \left(\frac{M_{\ast }}{10^{10}{\mathrm{M}}_{\odot }}\right)+\beta \log \left(\frac{s_5}{10^2\mathrm{kpc}}\right)\\ & +\gamma \log \left(\frac{M_{\mathrm{gas}}}{10^7{\mathrm{M}}_{\odot }}\right)+\log f\left(\frac{M_{\mathrm{BH}}}{10^8{\mathrm{M}}_{\odot }}\right).\end{array}$$

$$\log f_1(M_{\mathrm{BH}})=\{\begin{array}{ll}{k}_1\log M_{\mathrm{BH}}& M_{\mathrm{BH}}<M_{\mathrm{b}},\\ k_2\log M_{\mathrm{BH}}+(k_1-k_2)\log M_{\mathrm{b}}& M_{\mathrm{BH}}\ge M_{\mathrm{b}},\end{array}$$

$$\log f_2(M_{\mathrm{BH}})=k_1\log M_{\mathrm{BH}},$$

$$\begin{array}{rl}\log f_3& (M_{\mathrm{BH}})=\frac{1}{2}(k_1+k_2)\log M_{\mathrm{BH}}\\ & +\frac{1}{2}(k_2-k_1){\Delta }_{\mathrm{logM}}\ln \left(\cosh \left(\frac{\log M_{\mathrm{BH}}-\log M_{\mathrm{b}}}{{\Delta }_{\mathrm{logM}}}\right)\right).\end{array}$$

$$\begin{array}{rl}\log & \left({\mathrm{SFR}}_{\mathrm{model}}\right)=-2.50+0.21\log \left(\frac{M_{\ast }}{10^{10}{\mathrm{M}}_{\odot }}\right)\\ & -0.17\log \left(\frac{s_5}{100\mathrm{kpc}}\right)+1.03\log \left(\frac{M_{\mathrm{HI}}}{0.75\times 10^7{\mathrm{M}}_{\odot }}\right).\end{array}$$

Full text

††pagerange: Star Formation in Simulated Galaxies: Understanding the Transition to Quiescence at $3\times 10^{10}{\,\rm M}_{\odot}$–Star Formation in Simulated Galaxies: Understanding the Transition to Quiescence at $3\times 10^{10}{\,\rm M}_{\odot}$

1 Introduction

Understanding how galaxies evolve is a central problem in astronomy. It has been recognised for some time that energy feedback from stars and super-massive black holes (BH) in active galactic nuclei (AGN) is responsible for suppressing star formation in low- and high-mass galaxies, respectively. The transition between these two regimes occurs at a galaxy stellar mass $M_{*}\sim 3\times 10^{10}{\,\rm M}_{\odot}$ (e.g., Kauffmann et al., 2003; Baldry et al., 2006). The aim of this study is to improve our understanding of the physical processes that influence star formation in all galaxies and give rise to this transition scale.

Stars are formed according to an initial mass function (IMF), which describes the distribution of stellar masses formed (e.g., Salpeter, 1955; Chabrier, 2003; Kroupa, 2008; Kroupa et al., 2013). Very massive, but short-lived, stars explode as core-collapse supernovae, heating the surrounding interstellar medium (ISM), and enriching it with metals. Lower mass stars can also produce supernovae on longer (Gyr) timescales via white dwarf or neutron star progenitors. The injection of energy into the ISM by both stellar feedback (e.g., Heckman et al., 2000; Pettini et al., 2000; Ohyama et al., 2002) and AGN feedback (Lynds, 1967; Kraft et al., 2009; Feruglio et al., 2010; Cicone et al., 2012; Tombesi et al., 2013; Teng et al., 2014) can suppress further star formation on galactic scales (see also, e.g., Bicknell et al., 2000; Silk, 2013; Zubovas et al., 2013; Shabala et al., 2015; Bieri et al., 2016, who discuss the possibility of AGN-induced star formation on small scales). Stellar and AGN feedback can drive galactic outflows, which have been reproduced in cosmological simulations (e.g., Taylor & Kobayashi, 2015b). By quenching star formation, AGN feedback in massive galaxies causes them to have lower specific star formation rates, enhanced [$\alpha$/Fe], and redder colours than lower mass galaxies, which may lead to the observed downsizing phenomenon (Cowie et al., 1996; Juneau et al., 2005; Bundy et al., 2006; Stringer et al., 2009).

Stars form in the densest regions of giant molecular clouds in the ISM (Elmegreen & Scalo, 2004; Mac Low & Klessen, 2004; McKee & Ostriker, 2007; Federrath & Klessen, 2012; Padoan et al., 2014). AGN activity, too, relies on the accretion of gas onto the central super-massive BH. The availability of gas for both star formation and AGN therefore plays an important role in the subsequent evolution of a galaxy.

Most previous theoretical work has sought to explain observed star formation rates (SFRs) as a fraction of the molecular gas mass per (density-dependent) free-fall time (Krumholz et al., 2012; Federrath, 2013; Salim et al., 2015; Semenov et al., 2016). Such studies can explain local SFRs on the scale of giant molecular clouds, but a full understanding of the processes that affect the global galactic SFR is lacking. Observational data showed that there are two important, empirical relations: the Kennicutt-Schmidt law (Schmidt, 1959; Kennicutt, 1989; Kennicutt & Evans, 2012) and the star formation main sequence (SFMS; e.g., Elbaz et al., 2011; Zahid et al., 2012; Renzini & Peng, 2015). Most previous theoretical work has sought to explain the former, but a full understanding of the latter is lacking. Feldmann et al. (2016) use high-resolution simulations from the Feedback in Realistic Environments (FIRE) project (Hopkins et al., 2014) to study the processes governing the SFR of galaxies at $z\sim 2$, but do not explore the relative importance of these processes.

Galaxies form from massive gas clouds at high redshift, but can both gain and lose gas over the course of their life. Simulations based on the CDM cosmology predict large-scale flows of cold gas along filaments in the cosmic web (e.g., Dekel et al., 2009; Cen, 2014) that can efficiently supply a galaxy with gas without being shock-heated to the virial temperature (e.g., Birnboim et al., 2016). Groups and clusters are found at more massive nodes in the web, and can be fed by a greater number of filaments.

Gas is also transported into galaxies via mergers. So-called ‘wet’ (two gas-rich galaxies) and ‘damp’ (one gas-rich and one gas-poor galaxy) mergers can provide a galaxy with gas for both star formation and AGN activity, as well as potentially altering the morphology. During a merger, gas can be efficiently transported to the centre of the new galaxy, providing fuel for star formation and AGN activity (Toomre & Toomre, 1972; Combes et al., 1990; Mihos & Hernquist, 1994, 1996; Barnes & Hernquist, 1996; Hopkins et al., 2008). Major mergers, in which the galaxies have similar masses, are typically thought to be the main drivers of galaxy evolution (in terms of morphology, SFR, and AGN activity) in the local universe (e.g., Darg et al., 2010), though minor mergers contribute significantly to the cosmic star formation budget at low redshift (e.g., Kaviraj, 2014). In cosmological simulations, both major and minor mergers occur according to the hierarchical clustering of galaxy halos.

Note that, as well as the feedback processes described above, galaxies can lose gas when falling into a cluster. Tidal interactions with cluster members and the dark matter halo itself can have an effect (e.g., Moore et al., 1996), and ram pressure stripping of the interstellar medium by the hot intra-cluster medium can remove much of the halo gas (e.g., Gunn & Gott, 1972; Dressler & Gunn, 1983; Gavazzi et al., 1995; Boselli & Gavazzi, 2006). These effects are not well reproduced in smoothed particle hydrodynamics (SPH) simulations.

It is clear that there are a number of complex and competing processes that alter the amount of gas in a galaxy that can form stars or fuel BHs. What is less apparent is which of these processes is more important, and in what circumstances, in determining the SFR? In this paper, we aim to quantitatively answer this using cosmological hydrodynamical simulations that self-consistently solve the relevant physics. Section 2 gives a brief overview of the simulations analysed. In Section 3, we describe the analysis methodology including model formulae, and determine the model parameters in Section 4. Section 5 presents a discussion of the results, and in Section 6 we test our model against observational data. Finally, we give our conclusions in Section 7.

2 Simulations

The simulations used in this paper were introduced in Taylor & Kobayashi (2015a); they are a pair of cosmological, chemodynamical simulations, one of which includes a model for AGN feedback (Taylor & Kobayashi, 2014), but otherwise have indentical initial conditions and physics. Our simulation code is based on the SPH code gadget-3 (Springel, 2005), updated to include: star formation (Kobayashi et al., 2007), energy feedback and chemical enrichment from supernovae (SNe II, Ibc, and Ia, Kobayashi, 2004; Kobayashi & Nomoto, 2009) and hypernovae (Kobayashi et al., 2006; Kobayashi & Nakasato, 2011), and asymptotic giant branch (AGB) stars (Kobayashi et al., 2011); heating from a uniform, evolving UV background (Haardt & Madau, 1996); metallicity-dependent radiative gas cooling (Sutherland & Dopita, 1993); and a model for BH formation, growth, and feedback (Taylor & Kobayashi, 2014), described in more detail below. We use the IMF of stars from Kroupa (2008) in the range $0.01-120{\,\rm M}_{\odot}$, with an upper mass limit for core-collapse supernovae of $50{\,\rm M}_{\odot}$.

The initial conditions for both simulations consist of $240^{3}$ particles of each of gas and dark matter in a periodic, cubic box $25\,h^{-1}$ Mpc on a side, giving spatial and mass resolutions of $1.125\,h^{-1}$ kpc and $M_{\rm DM}=7.3\times 10^{7}\,h^{-1}{\,\rm M}_{\odot}$, $M_{\rm g}=1.4\times 10^{7}\,h^{-1}{\,\rm M}_{\odot}$, respectively. We employ a WMAP-9 $\Lambda$CDM cosmology (Hinshaw et al., 2013) with $h=0.7$, $\Omega_{\rm m}=0.28$, $\Omega_{\Lambda}=0.72$, $\Omega_{\rm b}=0.046$, and $\sigma_{8}=0.82$.

BHs form from gas particles that are metal-free and denser than a specified critical density, mimicking the most likely formation channels in the early Universe via direct collapse of a massive gas cloud (e.g., Bromm & Loeb, 2003; Koushiappas et al., 2004; Agarwal et al., 2012; Becerra et al., 2015; Regan et al., 2016; Hosokawa et al., 2016) or as the remnant of Population III stars (e.g., Madau & Rees, 2001; Bromm et al., 2002; Schneider et al., 2002). The BHs grow through gas accretion and mergers. The accretion rate is estimated assuming Bondi-Hoyle accretion:

[TABLE]

where $M_{\rm BH}$ is the mass of the BH, $\rho$ the local gas density, $c_{s}$ the local sound speed, and $v$ the speed of the BH relative to the local gas particles. We limit the accretion rate to the Eddington rate, given by

[TABLE]

where $m_{\rm p}$ is the proton mass, $\sigma_{\rm T}$ is the Thompson cross section, $\epsilon_{\rm r}$ denotes the radiative efficiency of the BH, which we set to $0.1$, and $t_{\rm Sal}\approx 45\,{\rm Myr}$ is the Salpeter time. Two BHs merge if their separation is less than the gravitational softening used, and their relative speed is less than the local sound speed. A fraction of the energy liberated by gas accretion is coupled to neighbouring gas particles in a purely thermal form.

3 A new simple model for the SFR

The main aim of this paper is to gain a quantitive understanding of the relative importance of how various physical processes affect the SFR of galaxies at the present day. We start from the SFMS, which relates the SFR and stellar mass of galaxies. Observations show that star-forming, late-type galaxies (LTG) form a narrow sequence with more massive galaxies having higher SFR (e.g., Elbaz et al., 2011; Zahid et al., 2012; Renzini & Peng, 2015), while early-type galaxies (ETG) occupy a region of parameter space below this sequence (e.g., Wuyts et al., 2011; Renzini & Peng, 2015). In our simulations, AGN feedback is responsible for quenching star formation and moving galaxies off the SFMS (Taylor & Kobayashi, 2016). Fig. 1 shows our simulated SFMS both with (black stars) and without (red diamonds) AGN feedback. The solid red line in Fig. 1 shows the linear bisector fit to the data from the simulation without AGN feedback, given by

[TABLE]

where $M_{*}$ is the galaxy stellar mass (see Taylor & Kobayashi (2016) for a detailed comparison between our simulated SFMS and observational data). Using only the data from this simulation to obtain equation (3) ensures that all galaxies are star-forming and lie on the SFMS, giving a relationship with no AGN contamination.

However, as discussed in Section 1, other factors, such as the mass of gas in a galaxy, and the environment the galaxy exists in, may also affect the SFR. Feedback from star formation, supernovae, and AGN activity can remove gas from galaxies (Taylor & Kobayashi, 2015b), depending on the galaxy mass and strength of feedback. Although the feedback energy in our AGN model is calculated from the instantaneous accretion rate, the current mass of gas in the galaxy likely depends more on the BH mass, since this reflects the entire accretion history of the BH. Similarly, the total stellar mass better reflects the integrated influence of stellar feedback on the gas than the current SFR.

Therefore, in order to quantify the relative importance of these processes, we suggest the following simple model to describe SFR in all galaxies:

[TABLE]

where $s_{5}$ is the 3-dimensional 5th-nearest neighbour distance, which we showed in Taylor & Kobayashi (in prep.) to be a good measure of environmental density (though numerous measures of environmental density have been proposed; see e.g., Haas et al., 2012), $M_{\rm gas}$ is the mass of gas within the galaxy, and $A$ is a constant. The function $f(M_{\rm BH})$ will be discussed in detail below.

To perform the fitting of our simulated data to this model, it is convenient to work with the logarithm111Throughout this paper we adopt the standard notation $\log x\equiv\log_{10}x$ and $\ln x\equiv\log_{e}x$. of equation (4):

[TABLE]

The coefficients $C$, $\alpha$, $\beta$, and $\gamma$ are constants that will be determined by fitting (note that $\alpha$, $\beta$ and $\gamma$ are the scaling exponents for $M_{*}$, $s_{5}$, and $M_{\rm gas}$, respectively, as defined in equation (4)). $M_{*}$, $s_{5}$, $M_{\rm gas}$, and $M_{\rm BH}$ are normalised by representative quantities, such that the distribution of $\log\left(M_{*}/10^{10}{\,\rm M}_{\odot}\right)$ (and similarly for $s_{5}$, $M_{\rm gas}$, and $M_{\rm BH}$) has a modal value of about zero; this too is a fitting convenience, and affects only the value of fitted parameter $C$.

The function $f(M_{\rm BH})$ in equations (4) and (5) reflects the fact that BHs must grow sufficiently massive before their feedback energy can not be efficiently radiated away. This is illustrated in Fig. 2, which shows the residuals of SFR in the simulation with AGN compared with the SFR estimated from equation (3), as a function of BH mass. There are two clear regimes: at $M_{\rm BH}\lesssim 10^{7}{\,\rm M}_{\odot}$, the residuals are distributed around [math], while at $M_{\rm BH}\gtrsim 10^{7}{\,\rm M}_{\odot}$ the residuals decrease with increasing $M_{\rm BH}$, and the SFR can be as much as two orders of magnitude below the SFMS.

This suggests the following form for $f(M_{\rm BH})$:

[TABLE]

where $k_{1}$ and $k_{2}$ are constant, and $M_{\rm b}$ is a break mass at which $f(M_{\rm BH})$ turns over222Given the normalisation of $M_{\rm BH}$ to $10^{8}{\,\rm M}_{\odot}$ in equation (5), $M_{\rm b}$ is also measured relative to $10^{8}{\,\rm M}_{\odot}$. In subsequent sections, the absolute value of $M_{\rm b}$ will be given, with the factor $10^{8}$ taken into account.. In Section 4, we also test the following forms of $f(M_{\rm BH})$:

[TABLE]

[TABLE]

Function $f_{2}$ is the special case of $f_{1}$ in which $k_{2}=k_{1}$, i.e. there is no break. $f_{3}$ is a continuous, smooth analogue to $f_{1}$, in which the gradient varies smoothly from $k_{1}$ to $k_{2}$, with the additional parameter $\Delta_{\rm logM}$ controlling the width over which the transition occurs. Equation (8) reproduces $f_{1}$ in the limit $\Delta_{\rm logM}\rightarrow 0$. The constant of integration333$f_{3}$ reproduces $f_{1}$ to within a constant. Matching at $\log M_{\rm BH}\rightarrow\pm\infty$ gives this constant as $\frac{1}{2}(k_{2}-k_{1})(\Delta_{\rm logM}\ln 2-\log M_{\rm b})$. In practice, we absorb this into $C$ to avoid unnecessary correlation between $C$ and the other parameters. that would appear in equation (8) can simply be absorbed by $C$ in equation (5). We show in Fig. 3 an example of each of the three functions described above.

Values of $C$, $\alpha$, $\beta$, $\gamma$, $k_{1}$, $k_{2}$, $\log M_{\rm b}$, and $\Delta_{\rm logM}$ are found using the amoeba routine (Press et al., 1992), minimising the quantity $\sum_{i}\left(\log{\rm SFR}_{i}-\log{\rm SFR}_{i,{\rm fit}}\right)^{2}$. In order to estimate uncertainties on the fitted parameters we employ a bootstrap resampling technique whereby the parameters are re-derived for a random selection of the data, with repeats, the same size as the original dataset. We do this for $5\times 10^{6}$ resamplings, in order to fully sample the parameter space of the models with the most free parameters.

4 Model Fits

Fitting the SFR of our simulated galaxies using the method described in Section 3 to the models described in equations (5)-(8) results in the parameters given in Table 1; the full parameter distributions for $f_{1}$ are shown in Fig. 4. The uncertainties are derived from the bootstrap resampling procedure described in Section 3; the first column of Table 1 gives the best-fit parameters to the simulated data, the second gives the mean and standard deviation from the bootstrap distributions, and third shows the modal value with asymmetric errors denoting the $16^{\rm th}$ and $84^{\rm th}$ percentiles. For all three models, the values of $C$, $\alpha$, $\beta$, and $\gamma$ are consistent with one another, and are well-constrained by the data.

We find that SFR depends most strongly on the amount of gas in a galaxy, scaling linearly with $M_{\rm gas}$ ($\gamma\approx 1$). Additionally, galaxies with large stellar mass and high environmental density (low $s_{5}$) show greater SFR, with the dependence on galaxy mass being the more important (i.e., $\alpha>|\beta|$). However, with $f_{1}$ and $f_{3}$, this is tempered by the strong negative correlation between SFR and $M_{\rm BH}$ in galaxies with BHs more massive than $M_{\rm b}$ (i.e., $k_{2}$ is large and negative). Such galaxies tend to be found in the densest environments, at the centre of clusters, implying that satellite galaxies within the cluster have the highest SFRs, while field galaxies tend to have lower SFRs at a given stellar mass. This is consistent with the observational result of Koyama et al. (2013), who found that the SFR of star-forming galaxies increases with environmental density (see also Peng et al., 2010; Wijesinghe et al., 2012, who find weak or no dependence of SFR on environment for star-forming galaxies). In all models, star formation in galaxies with low-mass BHs is not strongly affected, as indicated by the value $k_{1}\sim 0$, i.e. SFR follows the SFMS for $M_{\rm BH}\lesssim M_{\rm b}$ (see Fig. 2).

The relatively large uncertainties on the values of both $k_{2}$ $\log M_{\rm b}$ in $f_{1}$ are due to their distributions, shown in Fig. 4, being highly extended, and, in the case of $\log M_{\rm b}$, nearly bimodal. In addition to the main peak around $\log M_{\rm b}\approx 7.3$, there is an extended plateau of values near $\log M_{\rm b}\sim 8$. $\log M_{\rm b}\gtrsim 8$ is found when very few of the 38 galaxies with $M_{\rm BH}>10^{7}{\,\rm M}_{\odot}$ are included in a bootstrap realisation of the data set. In such cases, none of the data constrains $M_{\rm b}$ other than being larger than the largest $M_{\rm BH}$, and so $M_{\rm b}$ can take any arbitrarily large number without affecting how well the model fits the data. This is also responsible for the very extended distribution of $k_{2}$ seen in Fig. 4.

Function $f_{3}$ produces fairly consistent values with those of $f_{1}$, favouring a slightly higher value $\log M_{\rm b}\approx 7.6$, as well as predicting a wide range of black hole masses over which AGN feedback quenches star formation with $\Delta_{\rm logM}\sim 1-2$. Similar caveats apply to $k_{2}$, $\log M_{\rm b}$, and $\Delta_{\rm logM}$ in this model as for $k_{2}$ and $\log M_{\rm b}$ in $f_{1}$. Fortunately, however, the conclusions we draw in subsequent sections do not depend on the exact values of these parameters.

Fig. 5 shows the residuals for each of the models $f_{1}$ to $f_{3}$ as a function of $M_{\rm BH}$. The mean and standard deviation are shown for bins with width 0.5 dex, and the red band shows the range of $M_{\rm b}$ from Table 1. At $\log M_{\rm BH}\lesssim 7$, the three models have residuals that are distributed fairly uniformly around 0, though $f_{2}$ may show a very slight increasing trend with $M_{\rm BH}$. At higher masses, $f_{2}$, which assumes that there is no break mass $M_{\rm b}$, shows a significant systematic decrease in residuals with mass that is not seen in the residuals of $f_{1}$. A similar trend is seen for the residuals of $f_{3}$, though the residuals are consistent with 0 over the full range of $M_{\rm BH}$.

5 Discussion

5.1 BH growth and AGN feedback

We focus now on understanding the implications of the value of $\Delta_{\rm logM}$ from $f_{3}$. It is instructive to estimate a timescale over which BHs grow by a factor $10^{\Delta_{\rm logM}}\sim 10-100$. In Fig. 6 we show the mass assembly history of all simulated BHs with final masses greater than $10^{6}{\,\rm M}_{\odot}$, and the horizontal dashed and dot-dashed lines show $M_{\rm BH}=M_{\rm b}\times 10^{\pm\frac{1}{2}\Delta_{\rm logM}}$ assuming the best-fit and modal values, respectively. It is clear that, regardless of the exact values of $M_{\rm b}$ and $\Delta_{\rm logM}$ used, the timescale for BHs to grow by a factor $10^{\Delta_{\rm logM}}$ is at least several Gyr. This should not be interpreted, however, as the timescale over which an individual BH quenches star formation in its host, since we showed in Taylor & Kobayashi (2015a) that star formation in individual galaxies can be quenched abruptly. Rather, this shows that once they grow to $\sim M_{\rm b}$, BHs have the potential to exert influence over the evolution of their galaxy.

The fact that an individual BH can quench star formation in its host much more quickly than it takes to grow by a factor $10^{\Delta_{\rm logM}}$ once it reaches $M_{\rm BH}\approx M_{\rm b}$ implies that attaining a mass $M_{\rm b}$ is a necessary but not sufficient requirement for AGN feedback to become effective. This may suggest that local properties within a galaxy can affect exactly when AGN feedback shuts off star formation, or that a ‘trigger’ for strong AGN feedback may be required once the BH is sufficiently massive, with gas-rich galaxy mergers being a likely candidate (e.g., Comerford et al., 2015; Gatti et al., 2016). The cause of strong AGN activity, such as the AGN-driven outflows described in Taylor & Kobayashi (2015b), will be investigated in detail in a future work.

5.2 Galaxy transition at $3\times 10^{10}{\,\rm M}_{\odot}$

Figure 5 shows that the model that does not include a BH break mass, $f_{2}$, does not describe our simulated data well for $M_{\rm BH}\gtrsim M_{\rm b}$. Therefore we can conclude that the influence of BHs on star formation changes once the BH grows above $M_{\rm b}$. There is an empirical relationship between $M_{\rm BH}$ and the stellar bulge mass of its host galaxy, $M_{\rm bulge}$. Converting the range of fitted values of $M_{\rm b}$ from both $f_{1}$ and $f_{3}$ ($M_{\rm b}\approx 2-5\times 10^{7}{\,\rm M}_{\odot}$) into a stellar bulge mass using the relation given in Graham & Scott (2015) gives $M_{\rm bulge}\approx 1-2\times 10^{10}{\,\rm M}_{\odot}$. At these stellar masses, galaxies typically have bulge-to-disc ($B/D$) ratios of $0.3-0.4$ (Shen et al., 2003), corresponding to bulge-to-total ($B/T$) ratios of $0.23-0.29$ (note that ETGs can have much larger $B/T\sim 0.7$; see, e.g., Morselli et al., 2017). For such values of $B/T$, $M_{\rm b}$ corresponds to a galaxy mass that is consistent with the galaxy mass of $M_{*}=3\times 10^{10}{\,\rm M}_{\odot}$ found by Kauffmann et al. (2003) to separate ETGs and LTGs.

This consistency between our simulations and observational data implies that feedback from super-massive BHs may be responsible for the transition from LTGs to ETGs at $M_{*}=3\times 10^{10}{\,\rm M}_{\odot}$. Such a result is in broad agreement with other studies; Keller et al. (2016) find that even very strong supernova feedback in the form of superbubbles is inefficient in galaxies $M_{*}\gtrsim 4\times 10^{10}{\,\rm M}_{\odot}$, arguing that AGN feedback must dominate in more massive galaxies. Similarly, Bower et al. (2017) find that BH accretion is suppressed by supernova feedback in galaxies less massive than $3\times 10^{10}{\,\rm M}_{\odot}$, with AGN feedback becoming important once supernova feedback can no longer remove gas from the galaxy potential. In this paper, we have not investigated why strong AGN feedback is triggered in massive galaxies, finding only that it is triggered once a BH grows to $M_{\rm b}\approx 2-5\times 10^{7}{\,\rm M}_{\odot}$. The reason for the triggering will be investigated in detail in a future work.

5.3 Connection to the Kennicutt-Schmidt Law

The Kennicutt-Schmidt (KS) law is an empirical relation between SFR surface density $\Sigma_{\rm SFR}$, and gas surface density, $\Sigma_{\rm gas}$ (Schmidt, 1959; Kennicutt, 1989; Kennicutt & Evans, 2012). With $\gamma\approx 1$ (the exponent of $M_{\rm gas}$ in our model), dividing equation (4) by galaxy area gives $\Sigma_{\rm SFR}\sim\Sigma_{\rm gas}$. Observations typically find a super-linear relation when all gas is considered (e.g., Kennicutt, 1989, 1998), as we have done, while a linear or sub-linear relation tends to be found when only dense gas is considered (e.g., Bigiel et al., 2008; Shetty et al., 2013, 2014, but see also Liu et al. 2011; Momose et al. 2013). Shetty et al. (2014) analysed disk galaxies from the STING survey and concluded that galactic properties other than $\Sigma_{\rm gas}$ affect $\Sigma_{\rm SFR}$; such properties were included in our model, which may explain the difference between our linear relationship and observations (though there is significant scatter in the KS relation; see e.g. Krumholz et al., 2012; Federrath, 2013; Salim et al., 2015). We note for completeness that if we only consider SFR and $M_{\rm gas}$, we obtain a slope $1.54\pm 0.03$, in much better agreement with the observed super-linear relation of Kennicutt (1989, 1998).

5.4 The star formation main sequence

Regardless of the form of $f(M_{\rm BH})$ used, Table 1 shows that $\alpha=0.21$, i.e. $\log{\rm SFR}\sim 0.21\log M_{*}$. This is in contrast to observational estimates, which find a present-day slope around 0.75 – 1 (e.g., Elbaz et al., 2007, 2011; Zahid et al., 2012; Renzini & Peng, 2015). This is due to the correlations between stellar mass and other physical properties that are not accounted for in the standard SFMS, but are in this model, as was the case for the KS relation above. Our results suggest that stellar mass is not as important in determining galactic SFR as implied by the standard SFMS, and that there is an environmental dependence as well as strong gas mass dependence on the SFR of star-forming galaxies.

6 Comparison of our new SFR Model to observational data

We test our star formation models by comparing to observations. Ideally, we would wish to repeat the analysis of Section 4 on observational data to critically appraise our star formation model and our simulations. However, there are no overlapping surveys that can provide all of SFR, $M_{*}$, $s_{5}$, $M_{\rm gas}$, and $M_{\rm BH}$ for a large number of galaxies. By searching the literature and various current observational databases, we were able to find a sample of galaxies for which all the required measurements of SFR, $M_{*}$, $s_{5}$, and $M_{\rm gas}$ were available (except $M_{\rm BH}$). We use this set of galaxies to compare their measured SFR with that estimated from our model.

The Galaxy And Mass Assembly survey (GAMA, Driver et al., 2011) is an optical survey using the Anglo-Australian Telescope and AAOmega spectrograph (Sharp et al., 2006). It is $>$98 per cent complete to a depth of $r<19.4$ or $19.8$. The survey is spread over three $12^{\circ}$ by $4^{\circ}$ equatorial regions. Environmental density data are only available for 16,062 galaxies in the region centred at 15h right ascension. In this paper, we use aperture-corrected $M_{*}$ (Taylor et al., 2011), SFR (Gunawardhana et al., 2011), and $s_{5}$ (Brough et al., 2013) data for these galaxies.

For the gas masses, we use data from the Arecibo Legacy Fast ALFA survey (ALFALFA, Giovanelli et al., 2005). ALFALFA is a wide-field survey using the Arecibo radio telescope to map the 21 cm line of atomic hydrogen. Corresponding HI masses, $M_{\rm HI}$ are derived in Haynes et al. (2011); to convert this to total gas mass, we make the simplistic assumptions that 1) $M_{\rm HI}$ is equal to the total mass of hydrogen in the ISM, and 2) the ISM has primordial composition, so that $M_{\rm gas}=M_{\rm HI}/0.75$. We associate HI detections from ALFALFA with an optical counterpart in GAMA if they are separated by less than $0\aas@@fstack{\circ}01$ in RA and Dec. Matching the catalogues in this way gives 14 galaxies with measured $M_{*}$, SFR, $s_{5}$, and $M_{\rm gas}$.

We estimate the SFR of these galaxies assuming the best-fit parameters of Model 1, neglecting any contribution from BHs, i.e.,

[TABLE]

The most massive of these 14 galaxies has $M_{*}=1.7\times 10^{10}{\,\rm M}_{\odot}$, and so we expect to be in the parameter space where BH feedback does not play a crucial role for the SFR, i.e. $\log f\sim k_{1}M_{\rm BH}\sim 0$. Higher-mass galaxies with $M_{*}$, $s_{5}$, $M_{\rm gas}$, and simultaneous $M_{\rm BH}$ measurements are needed in the future to test our predicted dependence of the SFR on BH mass.

In Fig. 7 we show the measured SFR (SFRobs) and the SFR predicted by our model (SFRmodel) for these galaxies (filled circles with error bars). The agreement between the observed SFR and those predicted by our model is extremely encouraging, considering the assumptions used to calculate $M_{\rm gas}$ and the fact that $s_{5}$ from the GAMA catalogue is projected, whereas we used the 3-D $s_{5}$ in deriving our model parameters. Also shown in Fig. 7 are the predicted SFRs if $M_{\rm gas}$ or $M_{\rm gas}$ and $s_{5}$ are excluded (diamond and star symbols, respectively). In both cases, the model SFR is significantly less than is observed, and shows little evidence of correlation with the observed values. This reinforces the fact that $M_{\rm gas}$ is the most important quantity in determining SFR. The inclusion of $s_{5}$ alters the predicted SFR very little, but may be more important in more extreme environments such as galaxy groups and clusters (e.g., Schaefer et al., 2017).

In the coming decade, vast radio surveys such as the Square Kilometre Array (SKA) and its precursors will provide HI masses for tens of thousands of galaxies, allowing for a more rigorous test of our model. However, black hole masses are more difficult to measure, and are available for only a relatively small number of galaxies in the local Universe. Directly validating our full model including the effects of BHs will be a challenging, but important task for the near future.

7 Conclusions

We have analysed data from a cosmological hydrodynamical simulation that includes a detailed prescription for star formation and stellar feedback, as well as BH formation and AGN feedback. Our aim was to understand better the relative importance of factors that affect SFR of galaxies. We suggested that SFR could be influenced by the mass of the host galaxy, its environmental density, its gas content, and the mass of its BH. We proposed a relatively simple form for the relation between SFR and host galaxy properties (equation (5)), with which we were able to make quantitative comparisons.

We find that once a BH has grown sufficiently massive ($\sim 2-5\times 10^{7}{\,\rm M}_{\odot}$), it can be triggered to shut off star formation through AGN feedback; the details of such triggering will be investigated in future work. Such a BH mass corresponds, via the Maggorian relation (Graham & Scott, 2015), to a galaxy bulge mass of $\sim 1-2\times 10^{10}{\,\rm M}_{\odot}$, which is consistent with the observed galaxy mass above which ETGs make up the dominant fraction of the galaxy population ($\sim 3\times 10^{10}{\,\rm M}_{\odot}$, Kauffmann et al., 2003). It is important to note that these masses are emergent in our simulations; none of the input parameters of the baryon physics is chosen to yield such a transition mass, but are constrained by other observations.

The SFR depends strongly on the amount of gas in a galaxy, scaling approximately linearly with gas mass, while the dependence on stellar mass is weaker than predicted from the standard SFMS since we take into account the correlations between SFR and the additional physical quantities $M_{\rm gas}$, $s_{5}$, and $M_{\rm BH}$ by fitting all simultaneously. Our simulation cannot resolve the cold molecular gas that would form stars, and in this analysis we have treated all gas equally, regardless of temperature or density. In light of this, it will be useful in future works to analyse galaxies simulated with sufficient resolution to distinguish between gas phases.

We find that star-forming galaxies in high density regions have larger SFR than star-forming galaxies in the field at given mass, in agreement with observations (e.g., Koyama et al., 2013). This would most likely evolve with redshift, with massive galaxies showing the highest SFR in the past. The existence and strength of any such evolution will be investigated in a future work. Furthermore, the size of our simulation box limits the size of our most massive simulated cluster, and it would be informative to investigate if the trends seen here hold for more extreme environments.

We use observational data from the GAMA and ALFALFA surveys to compare the SFR observed with that predicted by our new model (see Fig. 7). There is excellent agreement between these two quantities, lending credence to the model we have presented. The sample is limited by the available observations, but will increase to tens of thousands in the era of the SKA. Further validation of our model will only be possible with a catalogue of BH masses for large numbers of galaxies.

Acknowledgements

We thank the referee for their useful comments, which improved the quality of this paper. C.F. gratefully acknowledges funding provided by the Australian Research Council’s Discovery Projects (grants DP150104329 and DP170100603). The simulations presented in this work used high performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Partnership for Advanced Computing in Europe (PRACE grant pr89mu), the Australian National Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, in the framework of the National Computational Merit Allocation Scheme and the ANU Allocation Scheme. This work has made use of the University of Hertfordshire Science and Technology Research Institute high-performance computing facility. PT thanks S. Lindsay for helpful discussions. Finally, we thank V. Springel for providing GADGET-3.