NIHAO XXVI: Nature versus nurture, the Star Formation Main Sequence and the origin of its scatter

It is now well established that most star-forming galaxies lie on the star formation main sequence (SFMS), which is a linear relation (in log-space) between the galaxy’s star formation rate (SFR) and stellar mass (Brinchmann et al., 2004; Noeske et al., 2007). The existence of the SFMS has been observed at $z=0$ (Renzini & Peng, 2015; Elbaz et al., 2011), but also at larger redshifts (Daddi et al., 2007; Elbaz et al., 2007; Schreiber et al., 2015; Tasca et al., 2015). The slope of the SFMS usually has values of 0.5-1, and shows no significant evolution with redshift, while its normalization is increasing with redshift. The scatter of the SFMS is uniform with 0.2-0.3 dex for all observed redshifts. Speagle et al. (2014) compile data from 25 different publications and study the evolution of the SFMS, they find that 0.1 dex of the scatter is due to the usage of different techniques among the different studies (‘inter publication scatter’), after correcting for this effect the intrinsic scatter is 0.2 dex.

Sparre et al. (2015) investigate the SFMS with the Illustris simulations and reproduce the SFMS at $z=0$ and 4, but for intermediate redshifts of $z=1$ and 2 the normalization is slightly lower than observed values. The scatter is about 0.2-0.3 dex and thus consistent with observations. The IllustrisTNG simulations (Donnari et al., 2019) reproduce the SFMS at $z=0$ with 0.3 dex scatter, the slope and normalization are consistent with observations. Their scatter is constant with stellar mass and decreasing with redshift, their normalization is systematically lower than observations at $z\sim 0.75-2$ by 0.2-0.5 dex. Also hydrodynamic simulation from Kannan et al. (2014) and Torrey et al. (2014) produce a tight relation between SFR and stellar mass.

Matthee & Schaye (2019) investigate the scatter of the SFMS and find that galaxies that are located above the SFMS tend to stay there for timescales of about 10 Gyr, thus it seems that “a galaxy’s SFR remembers its past SFR.” They also report that later forming haloes tend to host galaxies with a higher SFR, and that the scatter of the SFMS correlates with the dark matter halo formation time. Dutton et al. (2010) use a semi analytical model to investigate the scatter of the SFMS. They find that the scatter is mostly caused by variations in the galaxy’s gas accretion history, and is moreover dependent on the halo concentration.

In this paper we investigate the SFMS with the NIHAO galaxies. These provide a higher spatial resolution and extend to lower stellar masses than earlier studies. We show that the slope of the SFMS at higher stellar masses ($\gtrsim 10^{9}\,\mathrm{M}_{\odot}$) is already affected by AGN feedback, leading earlier works to predict a slope lower than unity. We show that by incorporating smaller stellar masses and considering galaxies unaffected by AGN feedback, the slope of the SFMS is approximately unity, and thus larger than predicted by earlier works. We also investigate the origin of the scatter of the SFMS. We confirm findings by Matthee & Schaye (2019) that the scatter of the SFMS is correlated with the halo formation time: haloes that form later have a higher than average SFR and vice versa. We provide a new interpretation for this correlation: Galaxies with a higher-than-average SFR (located above the SFMS) can be re-interpreted as galaxies with a lower-than-average stellar mass (located to the left of the SFMS). Thus later (earlier) forming galaxies have a lower (higher) than-average stellar mass, which is simply because they have had less (more) time to form stars than earlier (later) forming galaxies.

The outline of this work is as follows: In Section 2 we introduce the NIHAO simulations and describe how we measure the quantities used in this paper. In Section 3 we calculate the SFMS of the NIHAO galaxies, its slope, normalization and scatter, and compare our results with several observations. In Section 4 we investigate how and why galaxies deviate from the SFMS, i.e. what is causing the scatter in the relation between SFR and stellar mass. In Section 5 we summarize our findings.

We use the NIHAO suite of galaxy simulations, which consists of more than 150 zoom-in simulations of galaxies introduced by Wang et al. (2015) and extended by Blank et al. (2019). The initial conditions originate from cosmological simulations with box sizes of 60, 20 and 15 $\mathrm{Mpc}\,\mathrm{h}^{-1}$ and with $400^{3}$ particles that are evolved until redshift zero, then haloes from these boxes are selected and resimulated individually with a higher resolution and with gas particles. The resolutions of the zoom-in simulations are chosen such that each galaxy is resolved with about $10^{6}$ particles and we resolve the mass profile at $\leq$ 1 per cent of the virial radius. Thus we reach dark matter particle masses of $3.4\times 10^{3}$ to $1.7\times 10^{6}\,\mathrm{M}_{\sun}$ and dark matter softening lengths of 116 to 931 pc. The ratio of dark and gas particle mass equals the cosmological mass ratio of dark matter and baryons of $\Omega_{\mathrm{DM}}/\Omega_{\mathrm{b}}=5.48$, the ratio of dark and gas particle softening length equals $(\Omega_{\mathrm{DM}}/\Omega_{\mathrm{b}})^{1/2}=2.34$. We use a flat LCDM cosmology with parameters from the Planck Collaboration et al. (2014).

NIHAO is simulated with an updated version (Keller et al., 2014) of the TreeSPH code Gasoline2 (Wadsley et al., 2017). For gas cooling we consider hydrogen, helium, and various metal-lines in a uniform ultraviolet ionizing background (Shen et al., 2010), including photoionization, UV background heating (Haardt & Madau, 2012), and Compton cooling. Star formation occurs for gas particles that surpass a density and temperature threshold ($T<15000\,\mathrm{K}$, $n>10.3\,\mathrm{cm}^{-3}$) with a rate of $\dot{M_{\star}}=c_{\star}M_{\mathrm{gas}}t_{\mathrm{dyn}}^{-1}$, where $t_{\mathrm{dyn}}=(4\uppi G\rho)^{-1/2}$ is the gas particle’s dynamical time, $\rho$ its density, $M_{\mathrm{gas}}$ its mass and $c_{\star}=0.1$ the star formation efficiency. To model supernova feedback we use the blastwave formalism of Stinson et al. (2006), where star particles with $8<M_{\star}/M_{\sun}<40$ inject metals and thermal energy to surrounding gas particles 4 Myr after they have formed, subsequently the cooling of these gas particles is delayed for $\sim 30\,\mathrm{Myr}$. Following Stinson et al. (2013) star particles provide ‘early stellar feedback’  before they produce a supernova, where 13 per cent of the total stellar flux of $2\times 10^{50}\,\mathrm{erg}\,\mathrm{M}_{\sun}^{-1}$ is injected into the surrounding gas as thermal energy. The free parameters of the stellar and supernova feedback model have been adjusted to match the $M_{\star}$-$M_{{200}}$ relation for one Milky Way-like galaxy at $z=0$. For more details on the NIHAO project see Wang et al. (2015).

NIHAO has been very successful in reproducing galaxy properties for halo masses of $M_{200}\leq 2\times 10^{12}\,\mathrm{M}_{\sun}$, e.g., the stellar mass versus halo mass relation (Wang et al., 2015), the galaxy velocity function (Macciò et al., 2016), the Tully-Fisher relation (Dutton et al., 2017), the rotation curves of dwarf galaxies (Santos-Santos et al., 2018) and the stellar mass versus black hole mass relation (Blank et al., 2019). Strictly speaking the results presented in this paper pertain to the NIHAO galaxies that we use for our analysis. However, in the previous 25 publications of the NIHAO series we have demonstrated that our simulations produce realistic galaxies, thus we are confident that the conclusions of this paper might also be applied to real galaxies.

We calculate the SFR at a specific time as the mass of all star particles within 20 per cent of the galaxy’s virial radius that have formed within the last 200 Myr, divided by 200 Myr. Some galaxies have simulation outputs with a SFR of zero. In that case we estimate a lower limit of the SFR as the gas particle mass divided by 200 Myr, and show them in our figures as arrows. Our results are robust amongst different time intervals for calculating the SFR. Measuring with a shorter time interval merely increases the incidence of SFRs with value zero.

For Section 3 we use 78 NIHAO galaxies from Wang et al. (2015) and for Section 4 a subset thereof containing 50 galaxies for which we have a dark matter counterpart. Black holes were introduced into NIHAO by Blank et al. (2019), but only for Fig. 2 in Section 3 we use 52 of these galaxies to demonstrate that AGN feedback significantly influences the slope of the SFMS. Otherwise we only use the original NIHAO galaxies from Wang et al. (2015) that do not include black holes, for the following reasons: only galaxies that are actually star forming are part of the SFMS, whereas quenched galaxies form a distinct population in the stellar mass versus SFR plane. This usually poses great difficulties for observations and most simulations alike, as schemes have to be developed to filter quenched galaxies from the sample in order to be able to investigate the SFMS. Furthermore even galaxies that are still within the SFMS can already be significantly affected by AGN feedback and thus distort the view on the SFMS, which we will show in Section 3. The NIHAO sample of galaxies without black holes provide an opportunity to investigate the SFMS unaffected by AGN feedback and unbiased by selection effects and uncertainties due to the exclusion criteria of quenched galaxies. How black holes influence the SFR of galaxies and how AGN feedback quenches star formation will be investigated in a separate publication.

In this section we investigate the SFMS of the NIHAO galaxies, its slope, normalization, scatter and evolution, and compare our findings to several observations. Fig. 1 shows the relation between SFR and stellar mass of the NIHAO galaxies for four different redshifts. At all redshifts (i.e. for all 64 simulation outputs) we fit a linear relation to $\log$ SFR and $\log M_{\star}$, which provides a good match to the simulation data. The arrows denote a lower-limit estimate for the SFR for cases where the SFR is zero, we do not use these values for calculating the linear fit.

We also compare our results to several observations. However, most observations are only available for large stellar masses of $\gtrsim 10^{9}\,\mathrm{M}_{\odot}$, whereas the NIHAO simulations cover stellar masses down to $\sim 10^{4}\,\mathrm{M}_{\odot}$. Behroozi et al. (2013a, b) provide a compilation of various observations from different publications. At $z=0,1,2$ this data fits the NIHAO simulations quite well for stellar masses below $\sim 10^{11}\,\mathrm{M}_{\odot}$. At higher stellar masses the observed values drop to lower SFRs as star formation in these galaxies gets quenched by AGN feedback. This does not happen for the NIHAO galaxies as these do not include black holes. Only the normalization of our fit for stellar masses below $\sim 10^{9}\,\mathrm{M}_{\odot}$ is slightly lower than the Behroozi et al. (2013a, b) data. Also other observations (Speagle et al., 2014; Schreiber et al., 2015; Elbaz et al., 2007, 2011; Renzini & Peng, 2015; Daddi et al., 2007) coincide with the NIHAO simulations.

Fig. 2 shows the parameters of the linear fit, i.e. the slope, normalization (the SFR at stellar mass $10^{10}\,\mathrm{M}_{\odot}$) and the scatter, compared to several observed values: Speagle et al. (2014); Elbaz et al. (2007, 2011); Renzini & Peng (2015); Daddi et al. (2007). Note that Behroozi et al. (2013a, b) and Schreiber et al. (2015) do not fit a linear relation to their data, thus we omit them here.

Our slope is not evolving with time and has values of about unity, and thus is consistently higher than observed values, which range from about 0.5 to 1. This can be explained as follows: most observations cover quite high stellar masses of $\gtrsim 10^{9}\,\mathrm{M}_{\odot}$, where AGN feedback already starts to reduce the SFR of these galaxies, thus their relations are ‘bending down,’  yielding a smaller slope. Thus the ‘natural,’‘zero-order’  slope of the SFMS is around unity, which is then reduced to lower values at high stellar masses due to the influence of AGN feedback. Abramson et al. (2014) also report a declining slope with increasing stellar mass, which they attribute to mass quenching due to the build-up of a stellar bulge.

The normalization in the middle panel of Fig. 2 shows a good agreement of the NIHAO simulations with several observed values. The lower panel of Fig. 2 shows the scatter of the SFMS. Observations usually report a scatter of about 0.2-0.3 dex, that does not depend on redshift or stellar mass. However, the scatter of the NIHAO galaxies is around 0.35 for up to redshift one, and thus higher than observed values of 0.2-0.3. This is because in the NIHAO simulations lower stellar masses ($\lesssim 10^{9}\,\mathrm{M}_{\odot}$), which are usually not covered in most observations, yield a higher scatter than higher stellar masses, which then contributes to the overall scatter shown in Fig. 2. To demonstrate this effect we recalculate the scatter by using only galaxies with stellar masses $>10^{9}\,\mathrm{M}_{\odot}$, shown in the lower panel of Fig. 2 as dotted line. We again omit redshifts with less than 10 galaxies. Indeed the scatter is now much lower with about 0.25 dex and thus in agreement with observed values. Additionally we show in Fig. 3 the scatter of the SFMS as a function of stellar mass for three different redshifts. Each bin of stellar mass contains at least 10 galaxies. As the number of galaxies per bin is low, we estimate the variance of the scatter with the jackknife resampling method. For stellar masses $\gtrsim 10^{8}\,\mathrm{M}_{\odot}$ the scatter is increasing with decreasing stellar mass, and decreasing with redshift. We note that for very low stellar masses ($\lesssim 10^{7}\,\mathrm{M}_{\odot}$) the scatter might be artificially increased due to the stochasticity of star formation, as some galaxies only form few star particles.

To illustrate how AGN feedback reduces the slope of the SFMS we take 52 NIHAO galaxies from Blank et al. (2019) which have been simulated with black holes, and calculate the slope, intercept and scatter of the SFMS for these galaxies. We exclude completely quenched galaxies by only selecting those whose SFR is less than three times the scatter away from the fitted SFMS of Fig. 1. We also exclude those with a stellar mass $<10^{9}\,\mathrm{M}_{\odot}$ which are usually not covered in most observations. We furthermore only calculate these quantities when there are at least 10 galaxies fulfilling these criteria, therefore there is no data for low redshifts where most galaxies are quenched. This slope is shown in the upper panel Fig. 2, it is indeed significantly lower than the slope of around unity for the NIHAO galaxies without black holes, and fits well with the observed slope of Speagle et al. (2014). The intercept does not change significantly, only its standard deviation grows towards smaller redshifts. The galaxies deviate from the SFMS due to AGN feedback, leading to a less tight relation. Also the scatter is at or below the scatter for NIHAO galaxies without black holes, showing that the scatter is not dominated by the effects of AGN feedback.

In this section we investigate why and how galaxies deviate from the SFMS, and what is influencing its scatter. Our fits to the SFMS (Fig. 1) give us an average star formation rate SFR_avg as a function of redshift and stellar mass. Therefore we define the deviation from the average SFR as

Thus an ‘average’  galaxy would have a value of $\Delta\log\mathrm{SFR}=0$ for all times, or fluctuate slightly around zero. However, many galaxies deviate from this expected behavior as shown in Fig. 4, which shows $\Delta\log\mathrm{SFR}$ versus time for three galaxies.

The galaxy g3.44e10¹¹1The name of each galaxy indicates approximately its halo mass at $z=0$. (blue line) indeed fluctuates around zero throughout its lifetime, whereas galaxy g3.21e11 (orange line) is a chronic overproducer of stars, and galaxy g8.26e11 (green line) is a chronic underproducer. What is driving the different behavior of these galaxies? According to Dutton et al. (2010) the scatter of the SFMS is dependent on the halo concentration, therefore we investigate the correlation of $\Delta\log\mathrm{SFR}$ with the galaxies’ $z=0$ NFW (Navarro et al., 1996) halo concentration. We use the halo concentration of the corresponding dark-matter-only (DMO) simulations, as baryonic effects can lead to non-NFW dark matter profiles (e.g., Macciò et al., 2020). From the 78 galaxies used in the previous section 50 have a DMO counterpart.

The histogram in Fig. 5 shows the number of galaxies in each bin of $\Delta\log\mathrm{SFR}$, color coded with the average $z=0$ halo concentration of each bin. Here we use all 50 galaxies at all 64 snapshots, except those that have a SFR of zero.

Chronic overproducers, i.e. galaxies with $\Delta\log\mathrm{SFR}>0$, tend to have lower $z=0$ halo concentrations, whereas chronic underproducers ($\Delta\log\mathrm{SFR}<0$) tend to have higher $z=0$ halo concentrations. Fig. 6 shows a similar histogram, but only for six galaxies with halo mass range of 1-2$\times 10^{11}\,\mathrm{M}_{\sun}$. This shows that the effect shown in Fig. 5 is not caused by variations in galaxy mass. Fig. 7 furthermore shows the $z=0$ halo concentration versus $\Delta\log\mathrm{SFR}$ for all times. A Spearman correlation gives a p-value of 0.014, clearly showing the existence of a correlation between these two quantities.

We next investigate why the $z=0$ halo concentration is correlated with the SFR across cosmic time. The concentration is dependent on the cosmic matter density at which the halo forms, thus the concentration is also correlated with the halo formation time. We follow Wechsler et al. (2002) to calculate the halo formation time (i.e. the scale factor at which the halo forms). The mass growth of a halo as a function of the scale factor $a$ can be described as

where $M_{0}$ is the $z=0$ halo mass, and $a_{\mathrm{c}}$ can be interpreted as the scale factor at which the halo forms. We fit equation (2) to the halo growth histories of the DMO simulations. In Fig. 8 we show the halo growth histories and their fits for five galaxies, showing that equation (2) provides a reasonable description for the halo growth.

The fitting procedure gives us the formation time, i.e. the formation scale factor $a_{\mathrm{c}}$ for each galaxy. We then show the correlation between $z=0$ halo concentration and formation scale factor $a_{\mathrm{c}}$ in Fig. 9, and the correlation between halo formation scale factor $a_{\mathrm{c}}$ and $\Delta\log\mathrm{SFR}$ in Fig. 10. With a p-value of 0.002 this correlation is stronger than the correlation between $\Delta\log\mathrm{SFR}$ and the $z=0$ halo concentration.

Thus haloes that form later, i.e. have a lower halo concentration, have a higher SFR and haloes that form earlier, i.e. have a higher halo concentration, have a lower SFR than an average galaxy. This means that if a galaxy lives above or below the SFMS is set at birth, which can be explained as follows: Galaxies that are overproducers of stars seem to have a higher-than-average SFR, they are shifted above the SFMS. However, a perhaps better interpretation is that they are shifted to the left from the SFMS, thus they do not have a higher-than-average SFR, but a lower-than-average stellar mass (with respect to the SFMS of the same SFR).²²2Although galaxies above the SFMS grow faster (in stellar mass) than galaxies below it, their tracks in the SFR versus $M_{\star}$ plane do not necessarily cross, as the ratio of SFR and $M_{\star}$, i.e. their slope, is always around unity. However, their tracks can cross in the $M_{\star}$ versus time plane. This corresponds to a later formation time, because these galaxies did not have enough time yet to form sufficient stars. Likewise underproducers of stars do not have a lower-than-average SFR, but a higher-than-average stellar mass, because they have formed stars for a longer time. It is of course to some extent arbitrary if the deviation of galaxies from the SFMS is measured in $M_{\star}$ or SFR direction. However, choosing $M_{\star}$ offers a physically more plausible interpretation than choosing the SFR. Finally, according to Wang & Lilly (2020a, b) fluctuations on time-scales < 200 Myr might also contribute to the the scatter of the SFMS; therefore, the effects we describe in this paper might not be the sole cause of the scatter.

Thus the halo formation time influences the scatter of the SFMS. Following Matthee & Schaye (2019) we fit a linear relation to $\Delta\log\mathrm{SFR}$ as a function of halo formation scale factor, thus we can calculate a ‘corrected’  or ‘residual’  SFR as

where the SFRs are measured in $\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$. The upper panel of Fig. 11 shows the scatter of this corrected SFR compared to the original scatter from Fig. 2. The overall scatter is reduced by about 0.05 dex, in accordance with Matthee & Schaye (2019). However, the reduction in the scatter only reaches as far back as 5 Gyr. For earlier times both, the corrected SFR and the original SFR, yield about the same scatter. As a function of stellar mass (lower panel of Fig. 11) the scatter of the corrected SFR is lower only for $z=0$ down to $\sim 10^{8}\,\mathrm{M}_{\odot}$. For $z=1$ and $z=2$ no reduction in the scatter occurs, in agreement with the upper panel of Fig. 11.

In this paper we investigate the SFMS and the origin of its scatter by using the NIHAO suite of galaxy simulations. In the first part we compare the SFMS of the NIHAO galaxies with several observations. They generally agree with each other, but the slope of the SFMS is around unity for all redshifts and thus significantly higher than values derived from observations. This is because galaxies at high stellar masses ($\gtrsim 10^{9}\,\mathrm{M}_{\odot}$) are already influenced by AGN feedback, thus their tracks in the stellar mass versus SFR plane are not linear but ‘bending down,’  leading to lower observed slopes. The NIHAO galaxies without AGN feedback show that the ‘zero-order’  slope of the SFMS, unaffected by AGN feedback, is around unity. We confirm this by recalculating the SFMS with 52 NIHAO galaxies with black holes from Blank et al. (2019), whose slope is significantly smaller than unity and in agreement with observed slopes. Most observations report a scatter of 0.2-0.3 dex (see e.g. Donnari et al., 2019) independent of stellar mass, but the NIHAO simulations show that the scatter is only that small at large stellar masses. As we find that small stellar masses ($\lesssim 10^{9}\,\mathrm{M}_{\odot}$) have larger scatter, we report an overall total scatter of 0.35 dex. Recalculating the scatter by using only galaxies with $M_{\star}>10^{9}\,\mathrm{M}_{\odot}$ gives a scatter of 0.25 which is in agreement with observed values.

In the second part of the paper we investigate the origin of the scatter of the SFMS. We calculate the deviation from the fitted SFMS for each galaxy. This deviation does not fluctuate around zero as would be expected, but is either above or below zero for most of the evolution, i.e. these galaxies are either chronic over- or underproducers of stars. We find that the SFMS deviation correlates with the galaxy’s dark matter halo concentration at $z=0$ (see also Dutton et al., 2010), with overproducers having a low halo concentration and vice versa. We also confirm that halo concentration is anticorrelated with halo formation time (see also Wechsler et al., 2002). Thus a galaxy’s deviation from the SFMS is correlated with the halo formation time (see also Matthee & Schaye, 2019): overproducers of stars form later and underproducers form earlier. What is the origin of this behavior? Does “a galaxy’s SFR remembers its past SFR” (Matthee & Schaye, 2019)? This is not the case. We provide the following simple explanation for this phenomenon: instead of interpreting these galaxies as high-SFR or low-SFR (as they lie above or below the SFMS) they can be interpreted as low-stellar-mass or high-stellar-mass galaxies (as they lie to the left or to the right of the SFMS). That means that later forming galaxies have a lower stellar mass than average, and vice versa. This is easy to comprehend: Later forming galaxies have less time to form stars, thus they have a lower-than-average stellar mass, and earlier forming galaxies have more time to form stars, thus they have a higher-than-average stellar mass. It is therefore their nature (how and when these galaxies form), and not how they were nurtured (their local environment) that sets the journey of a galaxy in the SFR versus stellar mass plane.

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). A part of this research was carried out on the High Performance Computing resources at New York University Abu Dhabi. We used the software package pynbody Pontzen et al. (2013) for our analyses.

The data underlying this article will be shared on reasonable request to the corresponding author.