An evolving and mass dependent ${\rm \sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation for galaxies

TL;DR

This study uses EAGLE simulations to analyze how the scatter in specific star formation rates varies with galaxy mass and redshift, revealing an evolving U-shaped relation influenced by feedback processes.

Contribution

It introduces an evolving, mass-dependent ${ m \sigma_{sSFR}}$-${ m M_{ ext{star}}}$ relation, highlighting the roles of AGN and supernova feedback in galaxy star formation diversity.

Findings

The ${ m \sigma_{sSFR}}$-${ m M_{ ext{star}}}$ relation is U-shaped with a minimum at ${ m 10^{9.5} M_{ ext{sun}}}$.

Feedback mechanisms significantly influence the scatter at different mass scales.

Methodologies for measuring SFR/M* greatly affect the observed relation.

Abstract

The scatter (${\rm\sigma_{\text{sSFR}}}$) of the specific star formation rates (sSFRs) of galaxies is a measure of the diversity in their star formation histories (SFHs) at a given mass. In this paper we employ the EAGLE simulations to study the dependence of the ${\rm \sigma_{\text{sSFR}}}$ of galaxies on stellar mass (${\rm M_{\star}}$) through the ${\rm \sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation in $ {\rm z \sim 0-4}$. We find that the relation evolves with time, with the dispersion depending on both stellar mass and redshift. The models point to an evolving U-shape form for the ${\rm \sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation with the scatter being minimal at a characteristic mass $M^{\star}$ of ${\rm 10^{9.5}}$ ${\rm M_{\odot}}$ and increasing both at lower and higher masses. This implication is that the diversity of SFHs increases towards both at the low- and high-mass ends. We find that active galactic nuclei feedback is important for increasing the ${\rm \sigma_{\text{sSFR}}}$ for high mass objects. On the other hand, we suggest that SNe feedback increases the ${\rm \sigma_{\text{sSFR}}}$ of galaxies at the low-mass end. We also find that excluding galaxies that have experienced recent mergers does not significantly affect the ${\rm \sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation. Furthermore, we employ the combination of the EAGLE simulations with the radiative transfer code SKIRT to evaluate the effect of SFR/stellar mass diagnostics in the ${\rm \sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation and find that the ${\rm SFR/M_{\star}}$ methodologies (e.g. SED fitting, UV+IR, UV+IRX-$\beta$) widely used in the literature to obtain intrinsic properties of galaxies have a large effect on the derived shape and normalization of the ${\rm \sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation.

Tables

Table 1. Table 1: The EAGLE cosmological simulations used for this work

Run	L [Mpc]	NTOT	Feedback
L100N1504-Ref	100	2 $\times$ ${1504}^3$	AGN $+$ SN
L100N1504-Ref+SKIRT	100	2 $\times$ ${1504}^3$	AGN $+$ SN
L50N752-Ref	50	2 $\times$ ${752}^3$	AGN $+$ SN
L50N752-NoAGN	50	2 $\times$ ${752}^3$	No AGN $+$ SN
L50N752-OnlyAGN	50	2 $\times$ ${752}^3$	AGN $+$ No SN

Table 2. Table 2: Fraction of passive galaxies excluded in order to define a main sequence

z	0	0.350	0.615	0.865	1.400	2.000	3.000	4.000
${\sigma }_{\text{text}sSFR,MS}$-$M_{\star }$ sSFR Cut ($y{r}^{-1}$)	${10}^{-11.0}$	${10}^{-10.8}$	${10}^{-10.3}$	${10}^{-10.2}$	${10}^{-9.9}$	${10}^{-9.6}$	${10}^{-9.4}$	${10}^{-9.1}$
$F_{Passive}$, $\log \left(M_{\star }/M_{\odot }\right)={10}^{8.0-9.5}$	0.08	0.06	0.10	0.07	0.04	0.08	0.08	0.05
$F_{Passive}$, $\log \left(M_{\star }/M_{\odot }\right)={10}^{9.5-10.5}$	0.12	0.10	0.09	0.07	0.04	0.08	0.12	0.09
$F_{Passive}$, $\log \left(M_{\star }/M_{\odot }\right)={10}^{10.5-11.0}$	0.32	0.24	0.37	0.32	0.27	0.33	0.35	0.29
${\sigma }_{\text{text}sSFR,MS,Moderate}$-$M_{\star }$ sSFR Cut ($y{r}^{-1}$)	Guo et al. (2015)	${10}^{-11.0}$	${10}^{-11.0}$	${10}^{-11.0}$	${10}^{-11.0}$	${10}^{-11.0}$	${10}^{-11.0}$	${10}^{-11.0}$
$F_{Passive}$, $\log \left(M_{\star }/M_{\odot }\right)={10}^{8.0-9.5}$	0.02	0.03	0.02	0.03	0.02	0.02	0.01	0.01
$F_{Passive}$, $\log \left(M_{\star }/M_{\odot }\right)={10}^{9.5-10.5}$	0.03	0.05	0.03	0.02	0.01	0.01	0.01	0.01
$F_{Passive}$, $\log \left(M_{\star }/M_{\odot }\right)={10}^{10.5-11.0}$	0.17	0.16	0.17	0.10	0.04	0.02	0.01	0.01

Table 3. Table 3: Summary of the different observations used for this work.

Publication	Redshift range	Technique to obtain	$⟨{\sigma }_{\text{text}sSFR}⟩$ $\pm$ Uncertainty, ${\sigma }_{\text{text}sSFR}$ with $M_{\star }$ $\uparrow$ [dex]
	Stellar mass range	sSFRs and SFRs	Intrinsic or Observed ${\sigma }_{\text{text}sSFR}$, Shape of ${\sigma }_{\text{text}sSFR}$-$M_{\star }$
Noeske et al. (2007)	z = 0.32, 0.59, 1.0	EL+UV+$I{R}_{24\mu m}$	0.3 $\pm 0.05$ [0.3 $\to$ 0.3]
	$\log \left(M_{\star }/M_{\odot }\right)=9.5-11.45$		Observed, Constant
Rodighiero et al. (2010)	z = 1.47, 2.2	UV+$I{R}_{24\mu m}$	0.3 $\pm 0.05$ [0.3 $\to$ 0.3]
	$\log \left(M_{\star }/M_{\odot }\right)=9.5-11.45$		Observed, Constant
Guo et al. (2013)	z = 0.7	UV+$I{R}_{24\mu m}$	0.24 $\pm 0.04$ [0.182 $\nearrow$ 0.307]
	$\log \left(M_{\star }/M_{\odot }\right)=9.75-11.25$		Observed, increases with $M_{\star }$
Schreiber et al. (2015)	z = 0.5, 1.0, 1.5, 2.2, 3.0	$UV+I{R}_{SED}$	0.35 $\pm 0.03$ [0.29 $\nearrow$ 0.37]
	$\log \left(M_{\star }/M_{\odot }\right)=9.45-10.95$		Observed, increases with $M_{\star }$
Ilbert et al. (2015)	z = 0.3, 0.7, 0.9, 1.3,	$UV+I{R}_{SED}$	0.33 $\pm 0.03$ [0.22 $\nearrow$ 0.481]
	$\log \left(M_{\star }/M_{\odot }\right)=9.75-11.25$		Observed, increases with $M_{\star }$
Willett et al. (2015)	z	$U{V}_{SED}+H\alpha$	0.33 $\pm 0.03$ [0.52 $\searrow$ 0.37 $\nearrow$ 0.48]
	$\log \left(M_{\star }/M_{\odot }\right)=8.35-11.5$		Observed, U-Shape
Guo et al. (2015)		$H\alpha +I{R}_{22\mu m}$	0.44 $\pm 0.012$ [0.366 $\nearrow$ 0.557]
	$\log \left(M_{\star }/M_{\odot }\right)=8.85-10.75$		Observed, increases with $M_{\star }$
Kurczynski et al. (2016)	z = 0.75, 1.25, 1.75, 2.25, 2.75	SED fitting	0.40 $\pm 0.02$ [0.404 $\searrow$ 0.315 $\nearrow$ 0.435]
	$\log \left(M_{\star }/M_{\odot }\right)=6.85-10.25$		Intrinsic, redshift/mass dependent
Santini et al. (2017)	z = 1.65, 2.5, 3.5	UV + $\beta$ slope	0.42 $\pm 0.05$ [0.54 $\searrow$ 0.31]
	$\log \left(M_{\star }/M_{\odot }\right)=6.85-10.25$		Observed, decreases with $M_{\star }$
Boogaard et al. (2018)	z = 0.1 - 0.9	H$\alpha$+H$\beta$	0.44 ${\pm }_{0.04}^{0.05}$ [0.44 $\to$ 0.44]
	$\log \left(M_{\star }/M_{\odot }\right)=8.0-10.5$		Intrinsic, constant
Davies et al. (2019)	z 0.1	$SF{R}_{H\alpha }$	0.66 $\pm 0.02$ [0.74308562 $\searrow$ 0.53393775 $\nearrow$ 0.70720613]
	$\log \left(M_{\star }/M_{\odot }\right)=7.5-11.0$		Observed, U-shape
Davies et al. (2019)	z 0.1	$SF{R}_W$	0.44 $\pm 0.02$ [0.42797872 $\searrow$ 0.3490565 $\nearrow$ 0.53074792]
	$\log \left(M_{\star }/M_{\odot }\right)=7.5-11.0$		Observed, U-Shape

Full text

Cosmological simulations — Star formation — Galaxies — Surveys

††journal: ApJ

1 Introduction

The scatter ( ${\rm\sigma_{\text{sSFR}}}$) of the specific star formation rate (the ratio of star formation rate and stellar mass)-stellar mass (${\rm M_{\star}}$) relation provides a measurement for the variation of star formation across galaxies of similar masses with physical mechanisms important for galaxy evolution make their imprint to it. These processes include gas accretion, minor mergers, disc dynamics, halo heating, stellar feedback and AGN feedback. The above are typically dependent on galaxy stellar mass and cosmic epoch (Cano-Díaz et al., 2016; Abbott et al., 2017; Wang et al., 2017; Chiosi et al., 2017; García et al., 2017; Eales et al., 2018; Qin et al., 2018; Sánchez et al., 2018; Cañas et al., 2018; Blanc et al., 2019) and possibly affect the shape of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ differently. However, it is difficult to study and especially quantify the effect of the different prescriptions important for galaxy evolution to the scatter solely by using insights from observations.

In addition, the shape of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation and the value of the dispersion is a matter of debate in the literature. The scatter is usually reported to be constant (${\rm\sim 0.3}$ dex) with stellar mass in most studies, especially those which address the high and intermediate redshift Universe (${\rm z>1}$). For example, Rodighiero et al. (2010) and Schreiber et al. (2015), using mostly UV-derived SFRs, suggest that the dispersion is independent on galaxy mass and constant (${\rm\sim 0.3}$ dex) over a wide $M_{\star}$ range for ${\rm z\sim 2}$ star forming galaxies (${\rm 10^{9}-10^{11}}$ ${\rm M_{\odot}}$). Whitaker et al. (2012) reported a variation of 0.34 dex from Spitzer MIPS observations. Similarly, Noeske et al. (2007) and Elbaz et al. (2007) reported a 1${\rm\sigma}$ dispersion in ${\rm\log(SFR)}$ of around ${\rm 0.3}$ dex at ${\rm z\sim 1}$ for their flux-limited sample. However, other studies suggest that the dispersion tends to be larger for more massive objects and in the lower-redshift Universe (Guo et al., 2013; Ilbert et al., 2015) implying that mechanisms important for galaxy evolution are prominent and contribute to a variety of star formation histories for massive galaxies. On the other hand, Santini et al. (2017) suggested that the scatter decreases with increasing mass and this implies that mechanisms important for galaxy formation are giving a larger diversity of SFHs to low mass objects. In addition, Boogaard et al. (2018) using the MUSE Hubble Ultra Deep Field Survey suggest that the intrinsic scatter of the relation at the low mass end is $\sim 0.44$ dex and larger than what is typically found at higher masses. In disagreement with all the above, Willett et al. (2015) using the Galaxy Zoo survey for ${\rm z<0.085}$ find a dispersion that actually decreases with mass at ${\rm 10^{8}-10^{10}}$ ${\rm M_{\odot}}$ from ${\rm\sigma=0.45}$ dex to ${\rm 0.35}$ dex and increases again at ${\rm 10^{10}-10^{11.5}}$ ${\rm M_{\odot}}$ to reach a scatter of ${\rm\sim 0.5}$ dex. All the above observational studies have conflicting results and this is possibly because they are affected by selection effects, uncertainties originating from star formation rate diagnostics (Katsianis et al., 2016; Davies et al., 2017), different separation criteria for passive/star forming galaxies (Renzini & Peng, 2015) and usually focus on different redshifts and masses. In addition, the observed scatter can be different than the intrinsic value. More specifically, at the high-mass end an increased scatter can be inferred due to the uncertainties in removing passive objects, while at the low-mass end an increased scatter can be due to poor signal-to-noise ratio 111According to Kurczynski et al. (2016) the intrinsic scatter is $\sim 0.10$/$\sim 0.15$ dex lower than the observed value at $z\sim 0.5-1.0$/$\sim 2.5-3.0$.. Due to the above conflicting results and limitations it is almost impossible to decipher if there is an evolution in the scatter of the relation and if it is mass/redshift dependent or not, solely by relying on observations.

Cosmological simulations are able to reproduce realistic star formation rates and stellar masses of galaxies (Tescari et al., 2014; Katsianis et al., 2015; Furlong et al., 2015; Katsianis et al., 2017b; Pillepich et al., 2018) and thus are a valuable tool to address the questions related to the shape of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$, the value of the dispersion, and the way mechanisms important to galaxy evolution affect it. Simulations have limitations in resolution and box size. Thus they suffer from small number statistics of galaxies at a given mass, especially at the high-stellar mass end and cosmic variance due to finite box-size. However, despite their limitations the retrieved properties of galaxies do not suffer from poor S/N at the low mass end or uncertainties brought by different methodologies employed in observational studies (Katsianis et al., 2016, 2017a) and thus can provide a useful guide to future surveys or address controversies in galaxy formation Physics. Dekel et al. (2009) pointed out that the scatter of the specific star formation rate - stellar mass relation in cosmological simulations is about 0.3 dex and driven mostly by the galaxies’ gas accretion rates. Hopkins et al. (2014) using the FIRE zoom-in cosmological simulations studied the dispersion in the SFR smoothed over various time intervals and pointed out that the star formation main sequence and distribution of specific SFRs emerge naturally from the shape of the galaxies’ star formation histories, from $M_{\star}\sim 10^{8}-10^{11}$ $M_{\odot}$ at $z\sim 0-6$. The authors suggested that the scatter is larger on small timescales and masses, while dwarf galaxies ($<10^{8}$) exhibit much more bursty SFHs (and therefore larger scatter) due to stochastic processes like star cluster formation, and their associated feedback. Matthee & Schaye (2018) argued that the scatter of the main sequence ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation, defined by a sSFR cut in galaxies, is mass dependent and decreasing with mass at $z\sim 0$, while presenting a comparison between the EAGLE reference model and SDSS data. The authors suggested that the scatter of the relation originates from a combination of fluctuations on short time-scales (ranging from 0.2-2 Gyr) that are presumably associated with self-regulation from cooling, star formation and outflows (which nature is stochastic), but is dominated by long time-scale variations (Hopkins et al., 2014; Torrey et al., 2017), which is governed by the SFHs of galaxies, especially at high masses ($\log(M_{\star}/M_{\odot})>10.0)$222The authors pointed out that the total scatter of $\sim 0.4$ is driven by a combination of short and long time-scale variations, while for massive galaxies ($\log(M_{\star}/M_{\odot})=10.0-11.0$) the contribution of stochastic fluctuations (Kelson, 2014) is not significant ($<$ 0.1 dex). For objects with $\log(M_{\star}/M_{\odot})<10.0$ the contribution of fluctuations on short time scales, which nature is more stochastic, becomes relatively more important ($\sim$ 0.2 dex).. Dutton et al. (2010) using a semi-analytic model suggested that the scatter of the SFR sequence appears to be invariant with redshift and with a small value of ${\rm<0.2}$ dex.

The Virgo project Evolution and Assembly of GaLaxies and their Environments simulations (EAGLE, Schaye et al., 2015; Crain et al., 2015) is a suite of cosmological hydrodynamical simulations in cubic, periodic volumes ranging from 25 to 100 comoving Mpc per side. The reference model reproduces the observed star formation rates of ${\rm z\sim 0-8}$ galaxies (Katsianis et al., 2017b) and the evolution of the stellar mass function (Furlong et al., 2015). In addition, EAGLE allows us to investigate this problem with superb statistics (several thousands of galaxies at each redshift) and investigate different configurations that include different subgrid Physics. All the above provide a powerful resource for understanding the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation, address the shortcomings of observations, study its evolution across cosmic time and decipher its shape.

In this paper, we examine the dependence of the sSFR dispersion on ${\rm M_{\star}}$ using the EAGLE simulations (Schaye et al., 2015; Crain et al., 2015; Katsianis et al., 2017b). In section 2 we present the simulations used for this work. In section 3 we discuss the evolution of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation (subsection 3.1 presents the reference model) and how different feedback prescriptions (subsection 3.2) and ongoing mergers (subsection 3.3) affect its shape. In addition, we employ the EAGLE+SKIRT data (Camps et al., 2018), which represent a post-process of the simulations with the radiative transfer code SKIRT, in order to decipher how star formation rate and stellar mass diagnostics affect the relation in section 4. Finally in section 5 we draw our conclusions. Studies of the dispersion that rely solely on 2d scatter plots (i.e. displays of the location of the individual sources in the plane) are not able to provide a quantitative information of how galaxies are distributed around the mean sSFR and cannot account for galaxies that could be under-sampled or missed by selection effects so we extend our analysis of the dispersion of the sSFRs at different mass intervals on their distribution/histogram, namely the specific Star Formation Rate Function (sSFRF) by comparing the results of EAGLE with the observations present in Ilbert et al. (2015). In appendix A we present the evolution of the simulated specific star formation rate function in order to present how the sSFRs are distributed.

2 The EAGLE simulations used for this work

The EAGLE simulations track the evolution of baryonic gas, stars, massive black holes and non-baryonic dark matter particles from a starting redshift of ${\rm z=127}$ down to ${\rm z=0}$. The different runs were performed to investigate the effects of resolution, box size and various physical prescriptions (e.g. feedback and metal cooling). For this work we employ the reference model (L100N1504-Ref), a configuration with smaller boxsize (50 Mpc) but same resolution and physical prescriptions (L50N752-Ref), a run without AGN feedback (L50N752-NoAGN) and a simulation without SN feedback but AGN included (L50N752-OnlyAGN). We outline a summary of the different configurations in Table 1.

\startlongtable

The EAGLE reference simulation has $2\,\times\,1504^{3}$ particles (an equal number of gas and dark matter elements) in an L = 100 comoving Mpc volume box, initial gas particle mass of ${\rm m_{g}}=1.81\times 10^{6}\,{\rm M_{\odot}}$, and mass of dark matter particles of ${\rm m_{g}}=9.70\times 10^{6}\,{\rm M_{\odot}}$. The simulations were run using an improved and updated version of the N-body TreePM smoothed particle hydrodynamics code GADGET-3 (Springel, 2005) and employ the star formation recipe of Schaye & Dalla Vecchia (2008). In this scheme, gas with densities exceeding the critical density for the onset of the thermo-gravitational instability (${\rm n_{H}\sim 10^{-2}-10^{-1}\,cm^{-3}}$) is treated as a multi-phase mixture of cold molecular clouds, warm atomic gas and hot ionized bubbles which are all approximately in pressure equilibrium (Schaye, 2004). The above mixture is modeled using a polytropic equation of state ${\rm P=k\rho^{\gamma_{eos}}}$, where P is the gas pressure, ${\rm\rho}$ is the gas density and ${\rm k}$ is a constant which is normalized to ${\rm P/k=10^{3}\,cm^{-3}\,K}$ at the density threshold ${\rm n_{H}^{\star}}$ which marks the onset of star formation. The simulations adopt the stochastic thermal feedback scheme described in Dalla Vecchia & Schaye (2012). In addition to the effect of re-heating interstellar gas from star formation, which is already accounted for by the equation of state, galactic winds produced by Type II Supernovae are also considered. EAGLE models AGN feedback by seeding galaxies with BHs as described by Springel (2005), where seed BHs are placed at the center of every halo more massive than $10^{10}$ M${}_{\rm\odot}/{\rm h}$ that does not already contain a BH. When a seed is needed to be implemented at a halo, its highest density gas particle is converted into a collisionless BH particle inheriting the particle mass. These BHs grow by accretion of nearby gas particles or through mergers. A radiative efficiency of ${\rm\epsilon_{r}=0.1}$ is assumed for the AGN feedback. Other prescriptions such as inflow-induced starbursts, stripping of gas due to different interactions between galaxies, stochastic IMF sampling or variations to the AGN feedback prescription such as torque-driven accretion models (Anglés-Alcázar et al., 2017) or kinetic feedback (Weinberger et al., 2017) are not currently modeled in EAGLE. The EAGLE reference model and its feedback prescriptions have been calibrated to reproduce key observational constraints, into the present-day stellar mass function of galaxies (Li & White, 2009; Baldry et al., 2012), the correlation between the black hole and bulge masses (McConnell & Ma, 2013) and the dependence of galaxy sizes on mass (Baldry et al., 2012) at $z\sim 0$. Alongside with these observables the simulation was able to match many other key properties of galaxies in different eras, like molecular hydrogen abundances (Lagos et al., 2015), colors and luminosities at ${\rm z\sim 0.1}$ (Trayford et al., 2015), supermassive black hole mass function (Rosas-Guevara et al., 2016), angular momentum evolution (Lagos et al., 2017), atomic hydrogen abundances (Crain et al., 2017), sizes (Furlong et al., 2017), SFRs (Katsianis et al., 2017b), Large-scale outflows (Tescari et al., 2018) and ring galaxies (Elagali et al., 2018). In addition, Schaller et al. (2015) pointed out that there is a good agreement between the normalization and slope of the main sequence present in Chang et al. (2015) and the EAGLE reference model. Katsianis et al. (2016) demonstrated that cosmological hydrodynamic simulations like EAGLE, Illustris (Sparre et al., 2015) and ANGUS (Tescari et al., 2014) produce very similar results for the SFR-$M_{\star}$ relation with a normalization being in agreement with that found in observations at $z\sim 0-4$ (Kajisawa et al., 2010; De Los Reyes et al., 2014; Bauer et al., 2013; Salmon et al., 2015) and a slope close to unity. In this work, galaxies and their host halos are identified by a friends-of-friends (FoF) algorithm (Davis et al., 1985) followed by the SUBFIND algorithm (Springel et al., 2001; Dolag & Stasyszyn, 2009) which is used to identify substructures or subhalos across the simulation. The star formation rate of each galaxy is defined to be the sum of the star formation rate of all gas particles that belong to the corresponding subhalo and that are within a 3D aperture with radius 30 kpc (Schaye et al., 2015; Crain et al., 2015; Katsianis et al., 2017b).

3 The evolution of the intrinsic $\sigma_{\text{sSFR}}$-$M_{\star}$ relation

In this section we present the evolution of the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation in order to quantify and decipher its evolution and its dependence (or not) upon stellar mass and redshift. In subsection 3.1 we present the results of the EAGLE reference model and the compilation of observations used in this work, while in subsections 3.2 and 3.3 we focus on the effect of feedback and mergers on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation, respectively. For the simulations, we split the sample of galaxies at each redshift in 30 stellar mass bins from $\log(M_{\star}/M_{\odot})\sim 6.0$ to $\log(M_{\star}/M_{\odot})\sim 11.5$ (stellar mass bins of 0.18 dex at $z\sim 0$) and measure the 1 $\sigma$ standard deviation $\sigma(\log_{10}sSFR)$ in each bin.

We compare our simulated results with a range of observational studies in which usually different authors employ different techniques to exclude quiescent objects in their samples. In order to select only Star Forming Galaxies (SFGs) the authors may select only blue cloud galaxies (Peng et al., 2010), or use the BzK two-color selection (Daddi et al., 2007a), or the standard BPT (Baldwin et al., 1981) criterion, or employ the rest-frame UVJ selection (Whitaker et al., 2012; Schreiber et al., 2015) or an empirical color selection (Rodighiero et al., 2010; Guo et al., 2013; Ilbert et al., 2015) or specify a sSFR separation criterion (Guo et al., 2015). All these criteria should ideally cut out galaxies with low sSFR, but the thresholds differ significantly in value from one study to another, with some being redshift dependent and others not (Renzini & Peng, 2015). In order to surpass this complication and uncertainty of the effectiveness of excluding “passive objects” in observational studies, in the following subsection we present:

• •

the $\sigma_{\text{sSFR}}$-$M_{\star}$ of the full (Star forming $+$ Passive) EAGLE sample,

• •

the $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ relation of a “main sequence” defined by excluding passive objects with sSFRs $<10^{-11}$yr*-1*, $<10^{-10.8}$yr*-1*, $<10^{-10.3}$ yr*-1*, $<10^{-10.2}$ yr*-1*, $<10^{-9.9}$ yr*-1*, $<10^{-9.6}$yr*-1*, $<10^{-9.4}$yr*-1*, $<10^{-9.1}$yr*-1* for redshifts $z=0$, $z=0.35$, $z=0.615$, $z=0.865$, $z=1.485$, $z=2.0$, $z=3.0$, $z=4.0$, respectively (Furlong et al., 2015; Matthee & Schaye, 2018) and

• •

the $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ relation of a “main sequence” defined by the more conservative sSFR cuts (that ensures a more complete SF sample in the expense of some possible passive galaxy contamination) of $<10^{-11.0}yr^{-1}$ for $z>0$ and $\log 10(sSFR)<0.18\times\log 10(M_{\star})-4.5\,Gyr^{-1}$ for $z=0$ (Ilbert et al., 2015; Guo et al., 2015).

The differences between observations used in this work in terms of assumed methodology to exclude (or not) quiescent objects is described in the previous paragraph and table 2. The different data sets and methodologies to obtain SFRs for the same studies are described below and table 3. Noeske et al. (2007) used 2095 field galaxies from the Wavelength Extended Groth Strip International Survey (AEGIS) and derived SFRs from emission lines, GALEX, and Spitzer MIPS and 24$\mu$m photometry. Guo et al. (2013) used 12,614 objects from the multi-wavelength data set of COSMOS, while SFRs are obtained combining 24 $\mu$m and UV luminosities. Schreiber et al. (2015) used GOODS-North, GOODS-South, UDS, and COSMOS extragalactic fields and derived SFRs using UV+FIR luminosities. Ilbert et al. (2015) based their analysis on a 24 $\mu$m selected catalogue combining the COSMOS and GOODS surveys. The authors estimated SFRs by combining mid- and far-infrared data for 20,500 galaxies. Willett et al. (2015) used optical observations in the SDSS DR7 survey, while stellar masses and star formation rates are computed from optical diagnostics and taken from the MPA-JHU catalogue (Salim et al., 2007). Guo et al. (2015) used the SDSS data release 7, while in their studies SFRs are estimated from H$\alpha$ in combination with 22 $\mu$m observation from the Wide-field Infrared Survey Explorer. Santini et al. (2017) used the Hubble Space Telescope Frontier Fields, while SFRs were estimated from observed UV rest-frame photometry (Meurer et al., 1999; Kennicutt & Evans, 2012). Davies et al. (2019) used 9,005 galaxies from the GAMA survey (Driver et al., 2011, 2016a). The SFR indicators used are described in length in Driver et al. (2016b) and involve a) the SED fitting code magphys, b) a combination of Ultraviolet and Total Infrared (UV+TIR), c) H$\alpha$ emission line d) the Wide-field Infrared Survey Explorer (WISE) W3-band (Cluver et al., 2017), and e) extinction-corrected u-band luminosities derived using the GAMA rest-frame u -band luminosity and u-g colors.

3.1 The evolution of the $\sigma_{\text{sSFR}}$-$M_{\star}$, $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ relations in the EAGLE model

In Fig. 1 we present the evolution of the $\sigma_{\text{sSFR}}$-$M_{\star}$, which includes both passive and star forming objects (represented by the red solid line), and the “main sequence” $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ (represented by the black dotted line) and $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ (represented by the blue dashed line) relations in the EAGLE L100N1504-Ref at $\log(M_{\star}/M_{\odot})\sim 7.0$ to $\log(M_{\star}/M_{\odot})\sim 11.5$ and compare them with observations. The two vertical dotted lines represent the mass resolution limit of 100 baryonic particles ($\log(M_{\star}/M_{\odot})\sim 8.25$) and the statistic limit where there are fewer than 10 galaxies at the low- and high-mass ends (Furlong et al., 2015; Katsianis et al., 2017b). The shaded regions represent the 95 % bootstrap confidence interval for 5000 re-samples for the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation.

Starting from redshift $z=4.0$ (top left panel of Fig. 1) we see that the $\sigma_{\text{sSFR}}$-$M_{\star}$ of the reference model has a U-shape form. The dispersion decreases with mass at the $\log(M_{\star}/M_{\odot})\sim 8.5-9.5$ interval from $\sigma_{\text{sSFR}}=0.4$ to $0.2$ dex while it increases with mass at the $\log(M_{\star}/M_{\odot})\sim 9.5-10.5$ interval from $\sigma_{\text{sSFR}}\sim 0.2$ to $0.6$ dex. For the “main sequence”, $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ relation, defined by the exclusion of objects with $<10^{-9.1}yr^{-1}$ (Furlong et al., 2015; Matthee & Schaye, 2018), the scatter increases more moderately at the high mass end (from 0.2 dex to 0.3 dex), since passive objects which would increase the dispersion are excluded using a sSFR cut. We note that at this era the fraction of quiescent galaxies is expected to be small, thus the exclusion of quiescent objects should not affect significantly the relation (especially at the low mass end), and it is very possible that the above selection criterion is too strict. However, a more moderate cut of $<10^{-11.0}yr^{-1}$ (Ilbert et al., 2015) results in a relationship that is closer to that of the full EAGLE sample since the exclusion of quenched objects is less severe. We note that the observations of Santini et al. (2017) are broadly consistent with the $\sigma_{\text{sSFR,MS,moderate}}$-$M_{\star}$ (represented by the black dotted line) and $\sigma_{\text{sSFR}}$-$M_{\star}$ (represented by the red solid line) relations (green filled circles representing the observations of Santini et al. (2017) within 0.1 dex with respect the simulated results). A similar behavior is found for lower redshifts up to $z\sim 2.0$. This is possibly due to the fact that the moderate sSFR cut of $<10^{-11.0}yr^{-1}$ (Ilbert et al., 2015) resembles more closely the selection performed by Schreiber et al. (2015) and Santini et al. (2017).

At redshift $z\sim 1.4$ (left medium panel of Fig. 1) we find that there is an increment of scatter with mass for the EAGLE $\sigma_{\text{sSFR}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ relations at the high mass end ($\log(M_{\star}/M_{\odot})\sim 9.5-11.0$) from $\sim 0.2$ dex to $0.45$ and $0.65$ dex, respectively. On the other hand, the $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ relation has an almost constant scatter of $\sim 0.2$ dex with mass at $\log(M_{\star}/M_{\odot})\sim 10.0-11.0$. The observations for this mass interval (Rodighiero et al., 2010; Guo et al., 2013; Ilbert et al., 2015; Schreiber et al., 2015) typically lay between the two “main sequence” relations $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS}}$-$M_{\star}$, something that is possibly related to the uncertainties of removing quiescent objects in the literature (Renzini & Peng, 2015). The observations of Ilbert et al. (2015) and Schreiber et al. (2015) imply an increasing scatter with mass, while those of Noeske et al. (2007) and Rodighiero et al. (2010) a constant ($\sim 0.3$ dex). On the other hand, for the low- mass end ($\log(M_{\star}/M_{\odot})\sim 8-9.5$) there is a decrement with mass according to the EAGLE reference model in agreement with Santini et al. (2017). Similarly with higher redshifts the $\sigma_{\text{sSFR}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ relations have a U-shape form. The same behavior is found for lower redshifts up to $z\sim 0.35$, which reflects the fact that both low and high mass galaxies have a larger scatter/diversity of star formation histories than $M^{\star}$ (characteristic mass) objects. For $z\sim 0$ (bottom right panel of Fig. 1) the scatter is constant with mass for both the $\sigma_{\text{sSFR}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ relations for the $\log(M_{\star}/M_{\odot})\sim 8.5-9.5$ interval at $\sim 0.4$ dex. The scatter increases for the high mass end to $0.9$ dex when both passive and star forming galaxies are included. The increment is more moderate when cuts similar to the ones of Guo et al. (2015) are applied. On the other hand, the scatter decreases with mass for the “ main sequence” $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ relation from 0.4 dex to 0.2 dex.

We note that the EAGLE reference model suggests that both $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ relations are evolving with redshift and are not time independent. In Fig. 2 we present the evolution of the $\sigma_{\text{sSFR,MS,Moderate}}$ at $z\sim 0-4$ for the stellar masses of $\log(M_{\star}/M_{\odot})\sim 8.5$ (blue solid line), $\log(M_{\star}/M_{\odot})\sim 9.5$ (black dotted line) and $\log(M_{\star}/M_{\odot})\sim 10.5$ (red dashed line). The blue circles, black squares and red triangles represent the compilation of observations of (Guo et al., 2013; Ilbert et al., 2015; Schreiber et al., 2015; Willett et al., 2015; Santini et al., 2017) at $\log(M_{\star}/M_{\odot})\sim 8.5$, $\log(M_{\star}/M_{\odot})\sim 9.5$ and $\log(M_{\star}/M_{\odot})\sim 11.0$, respectively. Both in simulations and observations the $\sigma_{\text{sSFR,MS,Moderate}}$ at $\log(M_{\star}/M_{\odot})\sim 8.5$ increases steadily from $\sim$ 0.33 dex at $z\sim 4$ to $\sim$ 0.55 dex. For $\log(M_{\star}/M_{\odot})\sim 9.5$ the scatter remains almost constant at $\sim 0.2$ dex for $z\sim 4.0-1.5$ but increases up to 0.35 at $z\sim 0$ for the redshift interval of $z\sim 1.5-0$. Last, the scatter increases from 0.25 dex to 0.45 dex at $\log(M_{\star}/M_{\odot})\sim 11.0$.

In Fig. 3 we present the evolution of the $\sigma_{\text{sSFR,MS}}$ at $z\sim 0-4$ for the stellar masses of $\log(M_{\star}/M_{\odot})\sim 8.5$ (blue solid line), $\log(M_{\star}/M_{\odot})\sim 9.5$ (black dotted line) and $\log(M_{\star}/M_{\odot})\sim 11.0$ (red dashed line). The scatter increases at $z\sim 0-4$ from $\sim 0.25$ dex to $\sim 0.35$ dex at $\log(M_{\star}/M_{\odot})\sim 8.5$ and from $\sim 0.2$ dex to $\sim 0.3$ dex at $\log(M_{\star}/M_{\odot})\sim 9.5$. However, the $\sigma_{\text{sSFR,MS}}$ remains almost constant at $0.25$ dex at $\log(M_{\star}/M_{\odot})\sim 10.5$.

In conclusion, the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation (represented by the red solid line) has a U-shape form at all redshifts with the scatter typically decreasing at the $\log(M_{\star}/M_{\odot})\sim 8-9.5$ mass interval and increasing at $\log(M_{\star}/M_{\odot})\sim 9.5-11.0$. This implies a diversity of star formation histories for both the low and high mass ends. The above is supported by recent observations, while a complementary work using data from the GAMA survey demonstrates as well the U-shape form of the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation at $z\sim 0$ (Davies et al., 2019). We support that the shape found is not an observational effect and can be found as well in cosmological hydrodynamic simulations like EAGLE or semi-analytic models like SHARK (Lagos et al., 2018a; Davies et al., 2019). Galaxies at the low- and high- mass ends have a larger diversity of SFRs than intermediate mass objects, implying that there are multiple pathways for low and high mass galaxies to evolve. In the following sections we will demonstrate which prescriptions drive the U-Shape form of the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation. We note that the scatter increases with redshift and evolves with time (Fig. 2 and Fig. 3). The above findings are in agreement with the work of Kurczynski et al. (2016), which suggests a moderate increase in scatter with cosmic time from 0.2 to 0.4 dex across the epoch of peak cosmic star formation. When moderate sSFR cuts are employed in order to define a main sequence (e.g. $\sigma_{\text{sSFR,MS}}$-$M_{\star}$, $\sigma_{\text{sSFR}}$-$M_{\star}$, Ilbert et al., 2015; Guo et al., 2015) the U-shape form of the relation persist at $z>0.5$ but is less visible at lower redshifts.

3.2 The effect of the AGN and SN feedback on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation.

In this subsection we investigate the effect of Active Galactic Nuclei (AGN) feedback and Supernovae (SN) feedback on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation. To do so we compare 3 different configurations that have the same resolution and box size but have different subgrid physics recipes. These include:

• •

L50N752-Ref, which is a simulation with the same feedback prescriptions and resolution as the EAGLE reference model L100N1504-Ref (dark solid line of Fig. 4),

• •

L50N752-NoAGN, which has the same Physics and resolution as L50N752-Ref but does not include AGN feedback (dotted red line of Fig. 4),

• •

L50N752-OnlyAGN, which has the same Physics, boxsize and resolution as L50N752-Ref but does not include SN feedback (blue dashed line of Fig. 4).

In Fig. 4 we present the effect of feedback prescriptions on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation in the EAGLE simulations. Guo et al. (2013) and Guo et al. (2015) find an increasing scatter with mass and suggest that halo/stellar-mass dependent processes such as disk instabilities, bar-driven tidal disruption, minor mergers and major mergers/interactions or stellar feedback are important for large objects. However, the comparison between the reference model (L50N752-Ref, black solid line) and the configuration which does not include the AGN feedback prescription (L50N752-NoAGN, red dotted line) suggests that the effect of the AGN feedback mechanism is mostly responsible for increasing the scatter of the relation for high mass objects ($M_{\star}>10^{9.5}$ $M_{\odot}$). The prescription increases the diversity in star formation histories at the $\log(M_{\star}/M_{\odot})\sim 9.5-11.5$ regime (Figs 4 and 5) by creating a large number of quenched objects at the high mass end. We note that the absence of SN feedback would result in a more aggressive AGN feedback mechanism, which would significantly increase the dispersion for high mass galaxies (blue dashed line). This shows the interplay between the two different prescriptions and the finding is in agreement with the work of Bower et al. (2017), which suggests that black hole growth is suppressed by stellar feedback. If there are no galactic winds to decrease the accretion rate for a galaxy with a supermassive black hole, then the later will become bigger and its AGN feedback mechanism will affect more severely the sSFR of the object333In our AGN feedback scheme halos more massive than the $10^{10}$ $M_{\odot}/h$ threshold are seeded with a Supermassive Black Hole (SMBH), thus even relatively small galaxies at the $\log(M_{\star}/M_{\odot})\sim 8.5-9$ range which reside halos with $>10^{10}$ $M_{\odot}/h$ are affected, since the absence of the SNe feedback which would be important at this mass interval allows a fast and significant SMBH growth, which otherwise would not be possible.. The effect of the AGN feedback for the case of the L50N752-OnlyAGN run (in which SNe feedback is absent) is significant at z $\sim 1-4$, an epoch when the SFRs of simulated objects would be influenced significantly from SNe feedback (For more details see Fig. 7 in Katsianis et al. (2017b) which describes the effect of SNe feedback on the star formation rate function). In contrast with $\sigma_{\text{sSFR}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ (Figs 4 and 5), we find that the $\sigma_{\text{sSFR,MS}}$-$M_{\star}$ relation is not affected by feedback mechanisms (Fig. 6). The different fractions of quenched objects between configurations which employ different feedback prescriptions does not affect the comparison between them since the quenched objects, which increase the scatter, would in any simulation be excluded with the strict selection criterion adopted.

Santini et al. (2017) used the Hubble Space Telescope Frontier Fields to study the main sequence and its scatter. In contrast with Guo et al. (2013) and Ilbert et al. (2015) the authors found a decreasing scatter with mass at all redshifts they considered and they suggested that this behavior is a consequence of the smaller number of progenitors of low mass galaxies in a hierarchical scenario and/or of the efficient stellar feedback processes in low mass halos. Comparing the reference model (L50N752-Ref, black solid line) with the configuration which does not include the SN feedback mechanism (L50N752-OnlyAGN, blue dashed line) we see that the effect of this mechanism is indeed to increase the scatter of the relation with decreasing mass for low mass objects ($M_{\star}<10^{8.5}$ $M_{\odot}$). In addition, in low mass galaxies discrete gas accretion events trigger bursts of star formation which inject SNe feedback. Since feedback is very efficient in low mass galaxies this largely suppresses star formation until new gas is accreted (Faucher-Giguère, 2018, e.g.). However, the finding is below the resolution limit of 100 particles and higher resolution simulations are required to investigate this in the future (Figs 4 and 5).

In conclusion, AGN and SN feedback are playing a major role in producing the U-shape to the $\sigma_{\text{sSFR}}$-$M_{\star}$ and $\sigma_{\text{sSFR,MS,Moderate}}$-$M_{\star}$ relations described in subsection 3.1 and drive the evolution of the scatter. Both prescriptions give a range of SFHs both at the low- and high-mass ends with AGN feedback increasing the scatter mostly at the $\log(M_{\star}/M_{\odot})\sim 9.5-11.5$ interval.

3.3 The effect of excluding mergers on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation.

Guo et al. (2015) pointed out that in massive galaxies interactions like minor and major mergers can induce starbursts followed by strong stellar feedback that can contribute significantly to the spread of sSFRs. In contrast, according to the authors, lower-mass galaxies are supposed to be less affected by the above and thus should have a smaller dispersion of SFHs leading to a scatter that increases with mass for the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation. On the other hand, Peng et al. (2010) suggest that interacting/merging low-mass satellite galaxies are sensitive to environmental quenching and this could input a significant dispersion to the sSFRs at the low-mass end of the distribution. Orellana et al. (2017) reported that interactions between galaxies can affect the scatter for a range of masses. As Guo et al. (2015) suggest the effect of the above mechanisms to the sSFR dispersion are difficult to examine.

Qu et al. (2017) and Lagos et al. (2018b) studied the impact of mergers on mass assembly and angular momentum on the EAGLE galaxies. The authors found that the reference model is able to reproduce the observed merger rates and merger fractions of galaxies at various redshifts (Conselice et al., 2003; Kartaltepe et al., 2007; Ryan et al., 2008; Lotz et al., 2008; Conselice et al., 2009; de Ravel et al., 2009; Williams et al., 2011; López-Sanjuan et al., 2012; Bluck et al., 2012; Stott et al., 2013; Robotham et al., 2014; Man et al., 2016). Thus EAGLE can be used to study the effect of mergers on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation. We identify mergers using the merger trees available in the EAGLE database. These merger trees were created using the D-Trees algorithm of Jiang et al. (2014). Qu et al. (2017) described how this algorithm was adapted to work with EAGLE outputs. We consider that a merger (major or minor) occurs when the stellar mass ratio between the two merging systems, $\mu=M_{2}/M_{1}$ is above 0.1 (Crain et al., 2015). The separation criterion, $R_{merge}$, is defined as $R_{merge}=5\times R_{1/2}$, where $R_{1/2}$ is the half-stellar mass radius of the primary galaxy (Qu et al., 2017). The above selection method to identify mergers and to separate into major or minor mergers is widely assumed in the EAGLE literature (Jiang et al., 2014; Crain et al., 2015; Qu et al., 2017; Lagos et al., 2018b). The fraction of mergers in the reference model at $z\sim 0$ increases towards higher masses (Lagos et al., 2018b). Thus, it would be expected that the mechanism affects mostly the SFHs of high mass objects (Guo et al., 2015).

In Fig. 7 we present the effect of mergers on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation in the EAGLE simulations at $z\sim 0-4$. In the bottom panel the black line represents the reference model at $z\sim 0$, in which all galaxies are considered. The red dotted line represents the same but when ongoing mergers and objects that have experienced merging from $z\sim 0.365$ are not included in the analysis. The blue dashed line describes the results when objects that have experienced merging from $z\sim 0.5$ are excluded. The magenta represents the same when galaxies that have experienced merging from $z\sim 0.615$ are excluded. The above analysis allows us to quantify the effect of ongoing (at $z=0$), ongoing + recent ($z\sim 0-0.35$) and ongoing + recent + past ($z\sim 0-0.65$) mergers to the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation and has been done similarly for $z\sim 0,\,0.9,\,2,\,4$. We see that according to the reference model, mergers do not induce a significant dispersion in the star formation histories of galaxies. The above findings can be seen at all redshifts considered. This implies that recent mergers, despite their importance for galaxy formation and evolution, do not impart a significant scatter on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation.

4 The effect of SFR and stellar mass diagnostics on the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation.

To obtain the intrinsic properties of galaxies, observers have to rely on models for the observed light. Stellar masses are typically calculated via the Spectral Energy Density (SED) fitting technique, while for the case of SFRs different authors employ different methods [e.g. Conversion of IR+UV luminosities to SFRs (Arnouts et al., 2013; Whitaker et al., 2014), SED fitting (Bruzual & Charlot, 2003), conversion of UV, H$\alpha$ and IR luminosities (Katsianis et al., 2017b)]. However, there is an increasing number of reports that different techniques give different results, most likely due to systematic effects affecting the derived properties (Bauer et al., 2011; Utomo et al., 2014; Fumagalli et al., 2014; Katsianis et al., 2016; Davies et al., 2016, 2017). Boquien et al. (2014) argued that SFRs obtained from SED modeling, which take into account only FUV and U bands are overestimated. Hayward et al. (2014) noted that the SFRs obtained from IR luminosities (e.g. Noeske et al., 2007; Daddi et al., 2007a) can be artificially high. Ilbert et al. (2015) compared SFRs derived from SED and UV+IR, and find a tension reaching 0.25 dex. Guo et al. (2015) suggested that sSFR based on mid-IR emission may be significantly overestimated (Salim et al., 2009; Chang et al., 2015; Katsianis et al., 2017a, b). All the above uncertainties on the determination of intrinsic properties could possibly affect the observed/derived $\sigma_{\text{sSFR}}$-$M_{\star}$ relation.

Camps et al. (2016), Trayford et al. (2017) and Camps et al. (2018) presented a procedure to post-process the EAGLE galaxies and produce mock observations that describe how galaxies appear in various bands (e.g. GALEX-FUV, MIPS160, SPIRE500). The authors did so by performing a full 3D radiative transfer simulation to the EAGLE galaxies using the SKIRT code (Baes et al., 2003, 2011; Stalevski et al., 2012; Camps & Baes, 2015; Peest et al., 2017; Stalevski et al., 2016, 2017; Behrens et al., 2018). In this section we use the artificial SEDs present in Camps et al. (2018), in order to study how typical SFR/${\rm M_{\star}}$ diagnostics affect the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation. We stress that the EAGLE objects that were post-processed by SKIRT were galaxies with stellar masses $M_{\star}>10^{8.5}M_{\odot}$, above the resolution limit of 100 gas particles, and sufficient dust content. We use the Fitting and Assessment of Synthetic Templates (FAST) code (Kriek et al., 2009) to fit the mock SEDs to identify the SFRs and stellar masses of the EAGLE+SKIRT objects like in a range of observational studies (González et al., 2012, 2014; Botticella et al., 2017; Soto et al., 2017; Aird et al., 2018). Doing so enables us to evaluate the effect of different SFR/stellar mass diagnostics on the derived $\sigma_{\text{sSFR}}$-$M_{\star}$ relation and thus isolating the systematic effect on the $\sigma_{\text{sSFR}}$. We assume an exponentially declining SFH [${\rm SFR\sim exp(-t/tau)}$]444We note that this parameterization despite the fact that is commonly used in the literature could misinterpret old stellar light for an exponentially increasing contribution originating from a younger stellar population. This can underestimate by a factor of two both SFRs and stellar masses. In addition, exponentially SFHs may not be representative in describing star forming galaxies (Ciesla et al., 2017; Iyer & Gawiser, 2017; Carnall et al., 2018; Leja et al., 2018). (Longhetti & Saracco, 2009; Michałowski et al., 2012; Botticella et al., 2012; Fumagalli et al., 2016; Blancato et al., 2017; Abdurro’uf & Akiyama, 2019), a Chabrier IMF (Chabrier, 2003), a Calzetti et al. (2000) dust attenuation law (Mitchell et al., 2013; Sklias et al., 2014; Cullen et al., 2018; McLure et al., 2018b) and a metallicity 0.02 $Z_{\odot}$ (Ono et al., 2010; Greisel et al., 2013; Chan et al., 2016; McLure et al., 2018a). The above data are labelled in this work as L100N1504-Ref+SKIRT.

The black line of Fig. 8 represents the intrinsic scatter of the sSFRs of the EAGLE+SKIRT objects (intrinsic SFRs and stellar masses). The vertical lines define the mass interval at which the above data fully represent the total EAGLE sample. The shaded grey region represents the 95 % bootstrap confidence interval for 5000 re-samples for the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation. For clarity we present the above only for the reference model. The main results are the following:

• •

The magenta dashed line describes the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation when SFRs and stellar masses are both inferred using the SED fitting technique from the mock survey. The method is used in a range of observational studies to obtain the ${\rm SFR-M_{\star}}$ relation and its scatter (de Barros et al., 2014; Salmon et al., 2015). According to the EAGLE+SKIRT data we see that the shape remain relatively unchanged with respect to the intrinsic relation (black solid line), but the scatter is slightly underestimated ($\sim$ 0.5 dex) at the high mass end at $z\sim 4$ but overestimated by $\sim 0.15$ dex at the mass interval of $\log(M_{\star}/M_{\odot})\sim 10.0-11.0$. The above gives the false impression that the scatter increases more sharply with mass.

• •

The dark green dashed line represents the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation when SFRs are obtained using the FUV luminosity (Kennicutt & Evans, 2012) and the IRX-$\beta$ relation (Meurer et al., 1999; Bouwens et al., 2012; Katsianis et al., 2017a) while the stellar masses are calculated through the SED fitting technique. This combination to obtain properties is widely used in the literature to estimate the ${\rm SFR-M_{\star}}$ relation and its scatter (Santini et al., 2017). We see that this method implies an artificially higher $\sigma_{\text{sSFR}}$ (with respect the intrinsic-black solid line) at the low mass end. This gives an artificial mass independent ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation with a scatter of $\sim 0.35$ dex, which is not evolving significantly from $z\sim 2$ to $z\sim 0$.

• •

The red solid line is the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation retrieved when stellar masses are calculated via SED fitting and SFRs by combining the FUV and Infrared luminosity estimated from the $24\mu m$ luminosity (Kennicutt & Evans, 2012; Dale & Helou, 2002). Combining IR and UV luminosities to obtain SFRs in observations is a clasic method in the literature (Daddi et al., 2007b; Santini et al., 2009; Heinis et al., 2014) to obtain the ${\rm SFR-M_{\star}}$ relation. At $z\sim 4$ the scatter is underestimated with respect the reference black line by $0.2$ dex at high masses and imply a dispersion with a constant scatter around 0.15 dex. At lower redshifts there is a good agreement (within 0.05 dex) from the reference intrinsic black line.

• •

The blue dotted line describes the $\sigma_{\text{sSFR}}$-$M_{\star}$ relation when SFRs are derived from UV$+$TIR and stellar masses from the SED fitting technique. According to the EAGLE+SKIRT data this technique agrees well with the intrinsic relation, except redshift 4 where the derived relation is mass independent with a scatter of 0.2 dex. Similarly with the magenta line, which represents the results from SED fitting) the scatter is overestimated with respect the reference black line by $0.15$ dex for high mass objects at the mass interval of $\log(M_{\star}/M_{\odot})\sim 10.0-11.0$ for $z\sim 0.9$.

In conclusion, according to the EAGLE+SKIRT data, the inferred shape and normalization of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation can be affected by the methodology used to derive SFRs and stellar masses in observations. This can affect the conclusions about its shape and it is important for future observations to investigate this further (Davies et al., 2019). However, we note that having access to IR data and deriving SFRs and stellar masses from SED fitting or combined UV+IR luminosities typically give a ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation close to the intrinsic simulated relation and can probe successfully the shape of the relation for $\log(M_{\star}/M_{\odot})\sim 9.0-11.0$.

5 Conclusion and discussion

The ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation reflects the diversity of star formation histories for galaxies at different masses. However, it is difficult to decipher the true shape of the relation, the intrinsic value of the scatter and which mechanisms important for galaxy evolution govern it, solely by relying on observations. In this paper we presented the evolution of the intrinsic ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation employing the EAGLE suite of cosmological simulations and a compilation of multiwavelength observations at various redshifts. We deem the EAGLE suite appropriate for this study as it is able to reproduce the observed star formation rate and stellar mass functions (Furlong et al., 2015; Katsianis et al., 2017b) for a wide range of SFRs, stellar masses and redshifts. The investigation is not limited by the shortcomings encountered by galaxy surveys and address a range of redshifts and mass intervals. Our main conclusions are summarized as follows:

• •

In agreement with recent observational studies (Guo et al., 2013; Ilbert et al., 2015; Willett et al., 2015; Santini et al., 2017) the EAGLE reference model suggests that the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation is evolving with redshift and the dispersion is mass dependent (Section 3.1). This is in contrast with the widely accepted notion that the dispersion is mass/redshift independent with a constant scatter ${\rm\sigma_{\text{sSFR}}\sim 0.2-0.3}$ (Noeske et al., 2007; Elbaz et al., 2007; Rodighiero et al., 2010; Whitaker et al., 2012). We find that the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation has a U-shape form with the scatter increasing both at the high and low mass ends. Any interpretations of an increasing (Guo et al., 2013; Ilbert et al., 2015) or decreasing dispersion (Santini et al., 2017) with mass may be misguided, since they usually focus on limited mass intervals (Subsection 3.2). The finding about the U-shape form of the relation is supported by results relying on the GAMA survey (Davies et al., 2019) at $z\sim 0$.

• •

AGN and SN feedback are driving the shape and evolution of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation in the simulations (Subsection 3.2) . Both mechanisms cause a diversity of star formation histories for low mass (SN feedback) and high mass galaxies (AGN feedback).

• •

Mergers do not play a major role on the shape of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation (Subsection 3.3).

• •

We employ the EAGLE/SKIRT mock data to investigate how different ${\rm SFR/M_{\star}}$ diagnostics affect the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation. The shape of the relation remains relatively unchanged if both the SFRs and stellar masses are inferred through SED fitting, or combined UV+IR data. However, SFRs that rely solely on UV data and the IRX-${\rm\beta}$ relation for dust corrections imply a constant scatter with stellar mass with almost no redshift evolution. Methodology used to derive SFRs and stellar masses can affect the inferred ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation in observations and thus compromising the robustness of conclusions about its shape and normalization.

Acknowledgments

We would like to thank the anonymous referee for their suggestions and comments which improved significantly our manuscript. In addition, we would like to thank Jorryt Matthee for discussions and suggestions. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). A.K has been supported by the Tsung-Dao Lee Institute Fellowship, Shanghai Jiao Tong University and CONICYT/FONDECYT fellowship, project number: 3160049. G.B. is supported by CONICYT/FONDECYT, Programa de Iniciacion, Folio 11150220. V.G. was supported by CONICYT/FONDECYT iniciation grant number 11160832. X.Z.Z. thanks supports from the National Key Research and Development Program of China (2017YFA0402703), NSFC grant (11773076) and the Chinese Academy of Sciences (CAS) through a grant to the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. AK would like to thank his family and especially George Katsianis, Aggeliki Spyropoulou John Katsianis and Nefeli Siouti for emotional support. He would also like to thank Sophia Savvidou for her IT assistance.

Appendix A The evolution of the specific star formation rate function.

In this appendix, we base our analysis of the dispersion of the sSFRs at different mass intervals on their distribution/histogram, namely the specific Star Formation Rate Function (sSFRF), following Ilbert et al. (2015), because studies of the dispersion that rely solely on 2d scatter plots (i.e. displays of the location of the individual sources in the plane) are not able to provide a quantitative information of how galaxies are distributed around the mean sSFR and cannot account for galaxies that could be under-sampled or missed by selection effects. In section 3.1 we present the evolution of the ${\rm\sigma_{\text{sSFR}}}$-${\rm M_{\star}}$ relation at ${\rm z\sim 0-4}$ in order to visualize the scatter across galaxies, its shape and its evolution. We present the distribution of the sSFR of the EAGLE reference model (L100N1504-Ref) and compare it with observations (Ilbert et al., 2015) in Fig. 9 (${\rm z\sim 0.8-1.4}$) and Fig. 10 (${\rm z\sim 0.2-0.6}$). The EAGLE SFRs are reported to be 0.2 dex lower than observations but are able to replicate the observed evolution and shape of the Cosmic Star Formation Rate Density (Furlong et al., 2015; Katsianis et al., 2017b), evolution of star formation rate function (Katsianis et al., 2017b) and SFR-BHAR (McAlpine et al., 2017). Following, McAlpine et al. (2017) we decrease the observed sSFRs by 1.58 in order to focus solely on the scatter of the distribution and its shape. We note that this discrepancy may have its roots to the methodologies used by observers to obtain the intrinsic SFRs and stellar masses (Katsianis et al., 2017a).

In the top panels of Fig. 9 and Fig. 10 we present the comparison between the simulated and observed data. The black solid line corresponds to the results from the EAGLE reference model, while the orange dashed line (double-exponentional) and red dotted lines (log-normal profile) represent the fits for the observed distribution. The green vertical line mark the limits of the observations. The cyan area represents the 95 $\%$ bootstrap confidence interval for 1000 re-samples of the EAGLE SFRs, while the black errorbars represent the 1 sigma poissonian errors. We can see that the agreement between the simulated and observed sSFR functions is typically good at low mass objects but breaks down for galaxies more massive than ${\rm\log(M_{\star}/M_{\odot})>10}$, at ${\rm z>0.8}$ with the L100N1504-Ref distribution being shifted to lower sSFR and having larger peak values in comparison with the observations. Ilbert et al. (2015) were not able to directly constrain the full shape of the sSFR function, despite the fact that they combined both GOODS and COSMOS data. In most of the redshift and mass bins, the sSFR function is incomplete below the peak in sSFR. The authors tried to discriminate between a log-normal and a double-exponential profile but were not able to sample sufficiently low sSFR to see any advantage of using either one or the other parametrization. They argued that the fit with a double-exponential function is more suitable than the log-normal function at ${\rm z\sim 0}$. The EAGLE reference model points to the direction that the sSFRF at the different mass bins and redshifts follows a double-exponential function. However, for higher redshifts the simulated distributions are slightly flatter than the double-exponential fits of the observations. A double-exponential profile, that is not commonly used to describe the sSFR distribution (Ilbert et al., 2015), allows a significant density of star-forming galaxies with a low sSFR and the confirmation of this shape from future observations is important.