An Evolving and Mass-dependent σsSFR–M_⋆ Relation for Galaxies

The scatter (σ_sSFR) of the specific star formation rate (sSFR) (the ratio of star formation rate (SFR) and stellar mass)-stellar mass (M_⋆) relation provides a measurement of the variation of star formation across galaxies of similar masses with physical mechanisms important for galaxy evolution making their imprint to it. These processes include gas accretion, minor mergers, disk dynamics, halo heating, stellar feedback, and active galactic nucleus (AGN) feedback. The above are typically dependent on galaxy stellar mass and cosmic epoch (Cano-Díaz et al. 2016; Abbott et al. 2017; Chiosi et al. 2017; García et al. 2017; Wang et al. 2017; Cañas et al. 2019; Eales et al. 2018; Qin et al. 2018; Sánchez et al. 2018; Blanc et al. 2019), and may affect the shape of the σ_sSFR–M_⋆ differently. However, it is difficult to study—and especially so to 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 σ_sSFR–M_⋆ relation and the value of the dispersion is a matter of debate in the literature. The scatter is usually reported to be constant (∼0.3 dex) with stellar mass in most studies, especially those addressing the high- and intermediate-redshift universe (z > 1). For example, Rodighiero et al. (2010) and Schreiber et al. (2015), using mostly UV-derived SFRs, suggest that the dispersion is independent of galaxy mass and constant (∼0.3 dex) over a wide M_⋆ range for z ∼ 2 star-forming galaxies (SFGs) (10⁹–10¹¹ M_⊙). 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σ dispersion in $\mathrm{log}(\mathrm{SFR})$ of around 0.3 dex at z ∼ 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 give a larger diversity of star formation histories (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 ∼0.44 dex, which is larger than what is typically found at higher masses. In disagreement with all of the above, Willett et al. (2015), using the Galaxy Zoo survey for z < 0.085, find a dispersion that actually decreases with mass at 10⁸–10¹⁰ M_⊙ from σ = 0.45 to 0.35 dex and increases again at 10¹⁰–10^11.5 M_⊙ to reach a scatter of ∼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), and different separation criteria for passive/SFGs (Renzini & Peng 2015). Furthermore, they usually focus on different redshifts and masses, and 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.¹⁵ Because of the above conflicting results and limitations, if one relies solely on observations, it is almost impossible to decipher whether there is an evolution in the scatter of the relation and whether or not it is dependent on mass/redshift.

Cosmological simulations are able to reproduce realistic star formation rates and stellar masses of galaxies (Furlong et al. 2015; Katsianis et al. 2015; Tescari et al. 2014; Katsianis et al. 2017a; Pillepich et al. 2018), and thus are a valuable tool to address the questions related to the shape of the σ_sSFR–M_⋆, 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 signal-to-noise at the low-mass end or uncertainties brought by different methodologies employed in observational studies (Katsianis et al. 2016, 2017b), and thus they can provide a useful guide to future surveys or address controversies in galaxy formation physics. Dekel et al. (2009) point 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, have studied the dispersion in the SFR smoothed over various time intervals and point 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_⋆ ∼ 10⁸–10¹¹ M_⊙ at z ∼ 0–6. The authors suggest that the scatter is larger on small timescales and masses, while dwarf galaxies (<10⁸) exhibit much more bursty SFHs (and therefore larger scatter) due to stochastic processes like star cluster formation, and their associated feedback. Matthee & Schaye (2019) argue that the scatter of the main sequence σ_sSFR–M_⋆ relation, defined by a sSFR cut in galaxies, is mass-dependent and decreasing with mass at z ∼ 0, while presenting a comparison between the Evolution and Assembly of GaLaxies and their Environments (EAGLE) reference model and Sloan Digital Sky Survey (SDSS) data. The authors suggest that the scatter of the relation originates from a combination of fluctuations on short timescales (ranging from 0.2 to 2 Gyr) that are presumably associated with self-regulation from cooling, star formation, and outflows (which have a stochastic nature), but is dominated by long-timescale variations (Hopkins et al. 2014; Torrey et al. 2017) that are governed by the SFHs of galaxies, especially at high masses ($\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 10.0)$.¹⁶ Dutton et al. (2010), using a semi-analytic model, suggest that the scatter of the SFR sequence appears to be invariant with redshift and with a small value of <0.2 dex.

The Virgo project EAGLE (Crain et al. 2015; Schaye 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 z ∼ 0–8 galaxies (Katsianis et al. 2017a) 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 σ_sSFR–M_⋆ 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 M_⋆ using the EAGLE simulations (Crain et al. 2015; Schaye et al. 2015; Katsianis et al. 2017a). In Section 2, we present the simulations used for this work. In Section 3, we discuss the evolution of the σ_sSFR–M_⋆ relation (Section 3.1 presents the reference model) and how different feedback prescriptions (Section 3.2) and ongoing mergers (Section 3.3) affect its shape. In addition, we employ the EAGLE+SKIRT data (Camps et al. 2018), which represent a post-processing 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 quantitative information regarding how galaxies are distributed around the mean sSFR, and cannot account for galaxies that could be undersampled 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 the Appendix, we present the evolution of the simulated sSFRF, in order to illustrate how the sSFRs are distributed.

The EAGLE simulations track the evolution of baryonic gas, stars, massive black holes, and non-baryonic dark matter particles, from a starting redshift of z = 127 down to 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 a smaller box size (50 Mpc) but the same resolution and physical prescriptions (L50N752-Ref), a run without AGN feedback (L50N752-NoAGN), and a simulation without SN feedback but with AGN included (L50N752-OnlyAGN). We outline a summary of the different configurations in Table 1.

The EAGLE reference simulation has 2 × 1504³ particles (an equal number of gas and dark matter elements) in an L = 100 comoving Mpc volume box, initial gas particle mass of m_g = 1.81 × 10⁶ M_⊙, and mass of dark matter particles of m_g = 9.70 × 10⁶ M_⊙. 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 (n_H ∼ 10⁻²–10⁻¹ cm⁻³) is treated as a multiphase 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 $P=k{\rho }^{{\gamma }_{\mathrm{eos}}}$, where P is the gas pressure, ρ is the gas density, and k is a constant that is normalized to P/k = 10³ cm⁻³ K at the density threshold ${n}_{{\rm{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 reheating 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 black holes (BHs) as described by Springel (2005), where seed BHs are placed at the center of every halo more massive than 10¹⁰ M_⊙/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 _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, or stochastic initial mass function (IMF) sampling, and 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 ∼ 0. Alongside these observables, the simulation was able to match many other key properties of galaxies in different eras, such as molecular hydrogen abundances (Lagos et al. 2015), colors and luminosities at z ∼ 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. 2017a), large-scale outflows (Tescari et al. 2018), and ring galaxies (Elagali et al. 2018). In addition, Schaller et al. (2015) point 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) have 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_⋆ relation, with a normalization in agreement with that found in observations at z ∼ 0–4 (Kajisawa et al. 2010; Bauer et al. 2013; De Los Reyes et al. 2014; 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 are within a 3D aperture with radius 30 kpc (Crain et al. 2015; Schaye et al. 2015; Katsianis et al. 2017a).

In this section, we present the evolution of the σ_sSFR–M_⋆ relation, in order to quantify and decipher its evolution and its dependence (or not) upon stellar mass and redshift. In Section 3.1, we present the results of the EAGLE reference model and the compilation of observations used in this work, while Sections 3.2 and 3.3 are respectively focused on the effects of feedback and mergers on the σ_sSFR–M_⋆ relation. For the simulations, we split the sample of galaxies at each redshift into 30 stellar mass bins from $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 6.0$ to $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 11.5$ (stellar mass bins of 0.18 dex at z ∼ 0) and measure the 1σ standard deviation $\sigma ({\mathrm{log}}_{10}\mathrm{sSFR})$ in each bin.

We compare our simulated results with a range of observational studies in which different authors usually employ different techniques to exclude quiescent objects in their samples. In order to select only SFGs, the authors may select only blue cloud galaxies (Peng et al. 2010), or use the B–Z versus Z–K two-color selection (Daddi et al. 2007b) the standard Baldwin–Phillips–Terlevich (Baldwin et al. 1981) criterion, the rest-frame U–V versus V–J 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 the uncertainty regarding the effectiveness of excluding "passive objects" in observational studies, in the following subsection we present:

• 1.  
  The σ_sSFR–M_⋆ of the full (Star-forming + Passive) EAGLE sample.
• 2.  
  The σ_sSFR,MS–M_⋆ relation of a "main sequence" defined by excluding passive objects with sSFRs <10⁻¹¹ yr⁻¹, <10^−10.8 yr⁻¹, <10^−10.3 yr⁻¹, <10^−10.2 yr⁻¹, <10^−9.9 yr⁻¹, <10^−9.6 yr⁻¹, <10^−9.4 yr⁻¹, <10^−9.1 yr⁻¹ 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 2019).
• 3.  
  The σ_{sSFR,MS,Moderate}–M_⋆ relation of a "main sequence" defined by more conservative sSFR cuts (that ensure a more complete SF sample at the expense of some possible passive galaxy contamination) of < 10^−11.0 yr⁻¹ for z > 0 and $\mathrm{log}10(\mathrm{sSFR})\lt 0.18\times \mathrm{log}10({M}_{\star })\,-4.5\,{\mathrm{Gyr}}^{-1}$ for z = 0 (Guo et al. 2015; Ilbert et al. 2015).

The differences between observations used in this work in terms of assumed methodology to exclude (or include) 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 in Table 3. Noeske et al. (2007) used 2095 field galaxies from the All-wavelength Extended Growth Strip International Survey (AEGIS) and derived SFRs from emission lines, GALEX, Spitzer MIPS, and 24 μm photometry. Guo et al. (2013) used 12,614 objects from the multiwavelength data set of COSMOS, while SFRs are obtained by combining 24 μ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 μm selected catalog 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 catalog (Salim et al. 2007). Guo et al. (2015) used SDSS Data Release 7, while SFRs are estimated from Hα in combination with 22 μm observation from the Wide-field Infrared Survey Explorer (WISE). Santini et al. (2017) used the Hubble Space Telescope Frontier Fields, while SFRs are estimated from observed UV rest-frame photometry (Meurer et al. 1999; Kennicutt & Evans 2012). Davies et al. (2019) used 9005 galaxies from the Galaxy And Mass Assembly (GAMA) survey (Driver et al. 2011, 2016b). The SFR indicators used are described at length in Driver et al. (2016a). They include: (a) the Spectral Energy Density (SED) fitting code magphys, (b) a combination of Ultraviolet and Total Infrared (UV+TIR), (c) the Hα emission line, (d) the 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.

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
σsSFR,MS–M⋆ sSFR Cut (yr−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
FPassive, $\mathrm{log}({M}_{\star }/{M}_{\odot })={10}^{8.0\mbox{--}9.5}$	0.08	0.06	0.10	0.07	0.04	0.08	0.08	0.05
FPassive, $\mathrm{log}({M}_{\star }/{M}_{\odot })={10}^{9.5\mbox{--}10.5}$	0.12	0.10	0.09	0.07	0.04	0.08	0.12	0.09
FPassive, $\mathrm{log}({M}_{\star }/{M}_{\odot })={10}^{10.5\mbox{--}11.0}$	0.32	0.24	0.37	0.32	0.27	0.33	0.35	0.29
σsSFR,MS,Moderate–M⋆ sSFR Cut (yr−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
FPassive, $\mathrm{log}({M}_{\star }/{M}_{\odot })={10}^{8.0\mbox{--}9.5}$	0.02	0.03	0.02	0.03	0.02	0.02	0.01	<0.01
FPassive, $\mathrm{log}({M}_{\star }/{M}_{\odot })={10}^{9.5\mbox{--}10.5}$	0.03	0.05	0.03	0.02	0.01	0.01	0.01	<0.01
FPassive, $\mathrm{log}({M}_{\star }/{M}_{\odot })={10}^{10.5\mbox{--}11.0}$	0.17	0.16	0.17	0.10	0.04	0.02	0.01	0.01

Note. The fraction of passive galaxies (F_Passive) at each redshift excluded in order to define a main sequence. We adopt the effect of two different sSFR cuts. The first criterion (Furlong et al. 2015; Matthee & Schaye 2019) is used to define the σ_sSFR,MS–M_⋆ relation, while the second, more moderate criterion (Guo et al. 2015; Ilbert et al. 2015) is used to define the σ_{sSFR,MS,Moderate}–M_⋆ relation.

Download table as:  ASCIITypeset image

Table 3.  Summary of the Different Observations Used for This Work

Publication	Redshift Range	Technique to Obtain	$\langle {\sigma }_{\mathrm{sSFR}}\rangle$± Uncertainty, ${\sigma }_{\mathrm{sSFR}}$ with M⋆ $\uparrow$ (dex)
	Stellar Mass Range	sSFRs and SFRs	Intrinsic or Observed σsSFR, Shape of σsSFR–M⋆
Noeske et al. (2007)	z = 0.32, 0.59, 1.0	EL+UV+IR24 μm	0.3 ± 0.05 [0.3 $\to$ 0.3]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=9.5\mbox{--}11.45$		Observed, Constant
Rodighiero et al. (2010)	z = 1.47, 2.2	UV+IR24 μm	0.3 ± 0.05 [0.3 $\to$ 0.3]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=9.5\mbox{--}11.45$		Observed, Constant
Guo et al. (2013)	z = 0.7	UV+IR24 μm	0.24 ± 0.04 [0.182 $\nearrow$ 0.307]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=9.75\mbox{--}11.25$		Observed, increases with M⋆
Schreiber et al. (2015)	z = 0.5, 1.0, 1.5, 2.2, 3.0	$\mathrm{UV}+{\mathrm{IR}}_{\mathrm{SED}}$	0.35 ± 0.03 [0.29 $\nearrow$ 0.37]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=9.45\mbox{--}10.95$		Observed, increases with M⋆
Ilbert et al. (2015)	z = 0.3, 0.7, 0.9, 1.3,	$\mathrm{UV}+{\mathrm{IR}}_{\mathrm{SED}}$	0.33 ± 0.03 [0.22 $\nearrow$ 0.481]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=9.75\mbox{--}11.25$		Observed, increases with M⋆
Willett et al. (2015)	z < 0.085	${\mathrm{UV}}_{\mathrm{SED}}+{\rm{H}}\alpha$	0.33 ± 0.03 [0.52 $\searrow$ 0.37 $\nearrow$ 0.48]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=8.35\mbox{--}11.5$		Observed, U-Shape
Guo et al. (2015)	0.01 < z < 0.03	${\rm{H}}\alpha +{\mathrm{IR}}_{22\mu {\rm{m}}}$	0.44 ± 0.012 [0.366 $\nearrow$ 0.557]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=8.85\mbox{--}10.75$		Observed, increases with M⋆
Kurczynski et al. (2016)	z = 0.75, 1.25, 1.75, 2.25, 2.75	SED fitting	0.40 ± 0.02 [0.404 $\searrow$ 0.315 $\nearrow$ 0.435]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=6.85\mbox{--}10.25$		Intrinsic, redshift/mass-dependent
Santini et al. (2017)	z = 1.65, 2.5, 3.5	UV + β slope	0.42 ± 0.05 [0.54 $\searrow$ 0.31]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=6.85\mbox{--}10.25$		Observed, decreases with M⋆
Boogaard et al. (2018)	z = 0.1–0.9	Hα+Hβ	0.44 $\,\pm {\,}_{0.04}^{0.05}$ [0.44 $\to$ 0.44]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=8.0\mbox{--}10.5$		Intrinsic, constant
Davies et al. (2019)	z < 0.1	${\mathrm{SFR}}_{{\rm{H}}\alpha }$	0.66 ± 0.02 [0.74308562 $\searrow$ 0.53393775 $\nearrow$ 0.70720613]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=7.5\mbox{--}11.0$		Observed, U-shape
Davies et al. (2019)	z < 0.1	SFRW	0.44 ± 0.02 [0.42797872 $\searrow$ 0.3490565 $\nearrow$ 0.53074792]
	$\mathrm{log}({M}_{\star }/{M}_{\odot })=7.5\mbox{--}11.0$		Observed, U-Shape

Note. Summary of the different observations used in this work. Column 1 gives the publication name; Column 2 (top), the redshift range; Column 2 (bottom), the stellar mass range; Column 3, the technique to obtain galaxy sSFRs and SFRs; Column 4 (top), the average σ_sSFR ± uncertainty and the behavior of the scatter with increasing mass at the lowest redshift considered by the authors (Column 2, top); and Column 4 (bottom), the type of scatter used by the authors (intrinsic or observed) and the shape of the σ_sSFR–M_⋆.

Download table as:  ASCIITypeset image

3.1. The Evolution of the σ_sSFR–M_⋆, σ_sSFR,MS–M_⋆, and σ_{sSFR,MS,Moderate}–M_⋆ Relations in the EAGLE Model

In Figure 1, we present the evolution of the σ_sSFR–M_⋆, which includes both passive and star-forming objects (represented by the red solid line), and the "main sequence" σ_sSFR,MS–M_⋆ (represented by the black dotted line) and σ_{sSFR,MS,Moderate}–M_⋆ (represented by the blue dashed line) relations in the EAGLE L100N1504-Ref at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 7.0$ to $\mathrm{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 ($\mathrm{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. 2017a). The shaded regions represent the 95% bootstrap confidence interval for 5000 resamples of the σ_sSFR–M_⋆ relation.

Zoom In Zoom Out Reset image size

Starting from redshift z = 4.0 (top left panel of Figure 1), we see that the σ_sSFR–M_⋆ of the reference model has a U-shaped form. The dispersion decreases with mass at the $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8.5\mbox{--}9.5$ interval from σ_sSFR = 0.4 to 0.2 dex, while it increases with mass at the $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5\mbox{--}10.5$ interval from σ_sSFR ∼ 0.2 to 0.6 dex. For the "main sequence," σ_sSFR,MS–M_⋆ relation, defined by the exclusion of objects with <10^−9.1 yr⁻¹ (Furlong et al. 2015; Matthee & Schaye 2019), the scatter increases more moderately at the high-mass end (from 0.2 to 0.3 dex), because passive objects that would increase the dispersion are excluded using a sSFR cut. We note that the fraction of quiescent galaxies is expected to be small at this era, ergo the exclusion of quiescent objects should not significantly affect 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⁻¹ (Ilbert et al. 2015) results in a relationship that is closer to that of the full EAGLE sample, because the exclusion of quenched objects is less severe. We note that the observations of Santini et al. (2017) are broadly consistent with the σ_{sSFR,MS,moderate}–M_⋆ (represented by the black dotted line) and σ_sSFR–M_⋆ (represented by the red solid line) relations (green filled circles representing the observations of Santini et al. (2017) within 0.1 dex with respect to the simulated results). A similar behavior is found for lower redshifts up to z ∼ 2.0. This is possibly due to the fact that the moderate sSFR cut of < 10^−11.0 yr⁻¹ (Ilbert et al. 2015) more closely resembles the selection performed by Schreiber et al. (2015) and Santini et al. (2017).

At redshift z ∼ 1.4 (middle left panel of Figure 1), we find that there is an increment of scatter with mass for the EAGLE σ_sSFR–M_⋆ and σ_{sSFR,MS,Moderate}–M_⋆ relations at the high-mass end ($\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5\mbox{--}11.0$) from ∼0.2 dex to 0.45 and 0.65 dex, respectively. On the other hand, the σ_sSFR,MS–M_⋆ relation has an almost constant scatter of ∼0.2 dex with mass at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 10.0\mbox{--}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 σ_{sSFR,MS,Moderate}–M_⋆ and σ_sSFR,MS–M_⋆, something that may be 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) imply a constant (∼0.3 dex). On the other hand, for the low-mass end ($\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8\mbox{--}9.5$), the EAGLE reference model indicates a decrement with mass, in agreement with Santini et al. (2017). Similarly, at higher redshifts, the σ_sSFR–M_⋆ and σ_{sSFR,MS,Moderate}–M_⋆ relations have a U-shaped form. The same behavior is found for lower redshifts up to z ∼ 0.35, which reflects the fact that both low- and high-mass galaxies have a larger scatter/diversity of star formation histories than M^⋆ (characteristic mass) objects. For z ∼ 0 (bottom right panel of Figure 1), the scatter is constant with mass for both the σ_sSFR–M_⋆and σ_{sSFR,MS,Moderate}–M_⋆ relations for the $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8.5\mbox{--}9.5$ interval at ∼0.4 dex. The scatter increases to 0.9 dex for the high-mass end when both passive and SFGs 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" σ_sSFR,MS–M_⋆ relation from 0.4 to 0.2 dex.

We note that the EAGLE reference model suggests that both σ_{sSFR,MS,Moderate}–M_⋆ and σ_sSFR,MS–M_⋆ relations are evolving with redshift and are not independent of time. In Figure 2, we present the evolution of the σ_{sSFR,MS,Moderate} at z ∼ 0–4 for the stellar masses of $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8.5$ (blue solid line), $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5$ (black dotted line), and $\mathrm{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), and Santini et al. (2017) at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8.5$, $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5$, and $\mathrm{log}({M}_{\star }/{M}_{\odot })\,\sim 11.0$, respectively. In both simulations and observations, the σ_{sSFR,MS,Moderate} at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8.5$ increases steadily from ∼0.33 dex at z ∼ 4 to ∼0.55 dex. For $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5$, the scatter remains almost constant at ∼0.2 dex for z ∼ 4.0–1.5, but increases up to 0.35 at z ∼ 0 for the redshift interval of z ∼ 1.5–0. Last, the scatter increases from 0.25 to 0.45 dex at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 11.0$.

Zoom In Zoom Out Reset image size

In Figure 3, we present the evolution of the σ_sSFR,MS at z ∼ 0–4 for the stellar masses of $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8.5$ (blue solid line), $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5$ (black dotted line), and $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 11.0$ (red dashed line). The scatter increases at z ∼ 0–4 from ∼0.25 to ∼0.35 dex at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8.5$ and from ∼0.2 to ∼0.3 dex at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5$. However, the σ_sSFR,MS remains almost constant at 0.25 dex at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 10.5$.

Zoom In Zoom Out Reset image size

In conclusion, the σ_sSFR–M_⋆ relation (represented by the red solid line) has a U-shaped form at all redshifts, with the scatter typically decreasing at the $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 8\mbox{--}9.5$ mass interval and increasing at $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5\mbox{--}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 also demonstrates the U-shaped form of the σ_sSFR–M_⋆ relation at z ∼ 0 (Davies et al. 2019). We have shown that the shape found is not an observational effect and can also be found in cosmological hydrodynamic simulations like EAGLE or semi-analytic models like SHARK (Lagos et al. 2018b; 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-shaped form of the σ_sSFR–M_⋆ relation. We note that the scatter increases with redshift and evolves with time (Figures 2 and 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., σ_sSFR,MS–M_⋆, σ_sSFR–M_⋆, Guo et al. 2015; Ilbert et al. 2015), the U-shaped form of the relation persists at z > 0.5 but is less visible at lower redshifts.

3.2. The Effect of AGN and SN Feedback on the σ_sSFR–M_⋆ Relation

In this subsection, we investigate the effect of feedback from AGNs and supernovae (SNe) on the σ_sSFR–M_⋆ relation. To do so, we compare three different configurations that have the same resolution and box size but have different subgrid physics recipes. These include:

• 1.  
  L50N752-Ref, which is a simulation with the same feedback prescriptions and resolution as the EAGLE reference model L100N1504-Ref (dark solid line in Figure 4).
• 2.  
  L50N752-NoAGN, which has the same physics and resolution as L50N752-Ref, but does not include AGN feedback (dotted red line in Figure 4).
• 3.  
  L50N752-OnlyAGN, which has the same physics, box size, and resolution as L50N752-Ref, but does not include SN feedback (blue dashed line in Figure 4).

Zoom In Zoom Out Reset image size

In Figure 4, we present the effect of feedback prescriptions on the σ_sSFR–M_⋆ 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 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 that 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_⋆ > 10^9.5 M_⊙). The prescription increases the diversity in star formation histories at the $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5\mbox{--}11.5$ regime (Figures 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 latter will become larger and its AGN feedback mechanism will more severely affect the sSFR of the object.¹⁷ The effect of the AGN feedback for the case of the L50N752-OnlyAGN run (in which SNe feedback is absent) is significant at z ∼ 1–4, an epoch when the SFRs of simulated objects would be influenced significantly by SNe feedback (for more details, see Figure 7 in Katsianis et al. (2017a), which describes the effect of SNe feedback on the star formation rate function). In contrast with σ_sSFR–M_⋆ and σ_{sSFR,MS,Moderate}–M_⋆ (Figures 4 and 5), we find that the σ_sSFR,MS–M_⋆ relation is not affected by feedback mechanisms (Figure 6). The different fractions of quenched objects between configurations that employ different feedback prescriptions does not affect the comparison between them, because the quenched objects, which increase the scatter, would be excluded in any simulation with the strict selection criterion adopted.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 that 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_⋆ < 10^8.5 M_⊙). In addition, discrete gas accretion events in low-mass galaxies trigger bursts of star formation that inject SNe feedback. Because feedback is very efficient in low-mass galaxies, this largely suppresses star formation until new gas is accreted (e.g., Faucher-Giguère 2018). However, the finding is below the resolution limit of 100 particles, and higher-resolution simulations are required to investigate this in the future (Figures 4 and 5).

In conclusion, AGN and SN feedback play a major role in producing the U-shape of the σ_sSFR–M_⋆ and σ_{sSFR,MS,Moderate}–M_⋆ relations described in Section 3.1, and drive the evolution of the scatter. Both prescriptions give a range of SFHs at both the low- and high-mass ends, with AGN feedback increasing the scatter mostly at the $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.5\mbox{--}11.5$ interval.

3.3. The Effect of Excluding Mergers on the σ_sSFR–M_⋆ 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 those authors, lower-mass galaxies are supposed to be less affected by the above; thus, they should have a smaller dispersion of SFHs, leading to a scatter that increases with mass for the σ_sSFR–M_⋆ 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) report that interactions between galaxies can affect the scatter for a range of masses. As Guo et al. (2015) suggest, the effects of the above mechanisms upon the sSFR dispersion are difficult to examine.

Qu et al. (2017) and Lagos et al. (2018a) 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; Lotz et al. 2008; Ryan et al. 2008; Conselice et al. 2009; de Ravel et al. 2009; Williams et al. 2011; Bluck et al. 2012; López-Sanjuan 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 σ_sSFR–M_⋆ 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, μ = M₂/M₁, is above 0.1 (Crain et al. 2015). The separation criterion, R_merge, is defined as R_merge = 5 × 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 them 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. 2018a). The fraction of mergers in the reference model at z ∼ 0 increases toward higher masses (Lagos et al. 2018a). Thus, it would be expected that the mechanism affects mostly the SFHs of high-mass objects (Guo et al. 2015).

In Figure 7, we present the effect of mergers on the σ_sSFR–M_⋆ relation in the EAGLE simulations at z ∼ 0–4. In the bottom panel, the black line represents the reference model at z ∼ 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 ∼ 0.365 are not included in the analysis. The blue dashed line describes the results when objects that have experienced merging from z ∼ 0.5 are excluded. The magenta represents the same when galaxies that have experienced merging from z ∼ 0.615 are excluded. The above analysis allows us to quantify the effect of ongoing (at z = 0), ongoing + recent (z ∼ 0–0.35), and ongoing + recent + past (z ∼ 0–0.65) mergers to the σ_sSFR–M_⋆ relation, and has been done similarly for z ∼ 0, 0.9, 2, and 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 σ_sSFR–M_⋆ relation.

Zoom In Zoom Out Reset image size

vTo obtain the intrinsic properties of galaxies, observers have to rely on models for the observed light. Stellar masses are typically calculated via the SED fitting technique, while different authors employ different methods to calculate SFRs: e.g., conversion of IR+UV luminosities to SFRs (Arnouts et al. 2013; Whitaker et al. 2014); SED fitting (Bruzual & Charlot 2003); or conversion of UV, Hα, and IR luminosities (Katsianis et al. 2017a)). 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; Fumagalli et al. 2014; Utomo et al. 2014; Davies et al. 2016, 2017; Katsianis et al. 2016). Boquien et al. (2014) argue that SFRs obtained from SED modeling, which take into account only FUV and U bands, are overestimated. Hayward et al. (2014) note that the SFRs obtained from IR luminosities (e.g., Noeske et al. 2007; Daddi et al. 2007b) 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) suggest that sSFR based on mid-IR emission may be significantly overestimated (Salim et al. 2009; Chang et al. 2015; Katsianis et al. 2017a, 2017b). All the above uncertainties on the determination of intrinsic properties could possibly affect the observed/derived σ_sSFR–M_⋆ relation.

Camps et al. (2016), Trayford et al. (2017), and Camps et al. (2018) have 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, MIPS₁₆₀, SPIRE₅₀₀). The authors have done 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, 2016, 2017; Camps & Baes 2015; Peest et al. 2017; Behrens et al. 2018). In this section, we use the artificial SEDs presented in Camps et al. (2018), in order to study how typical SFR/M_⋆ diagnostics affect the σ_sSFR–M_⋆ relation. We stress that the EAGLE objects that were post-processed by SKIRT were galaxies with stellar masses M_⋆ > 10^8.5 M_⊙, above the resolution limit of 100 gas particles, and having sufficient dust content. We use the Fitting and Assessment of Synthetic Templates code (Kriek et al. 2009) to fit the mock SEDs to identify the SFRs and stellar masses of the EAGLE+SKIRT objects, as has been done 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 σ_sSFR–M_⋆ relation and thus isolate the systematic effect on the σ_sSFR. We assume an exponentially declining SFH [SFR ∼ exp(−t/tau)]¹⁸ (Longhetti & Saracco 2009; Botticella et al. 2012; Michałowski 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_⊙ (Ono et al. 2010; Greisel et al. 2013; Chan et al. 2016; McLure et al. 2018a). The above data are labeled in this work as L100N1504-Ref+SKIRT.

The black line of Figure 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 gray region represents the 95% bootstrap confidence interval for 5000 resamples for the σ_sSFR–M_⋆ relation. For clarity, we present the above only for the reference model. The main results are the following:

• 1.  
  The magenta dashed line describes the σ_sSFR–M_⋆ 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 SFR–M_⋆ relation and its scatter (de Barros et al. 2014; Salmon et al. 2015). According to the EAGLE+SKIRT data, we see that the shape remains relatively unchanged with respect to the intrinsic relation (black solid line), whereas the scatter is slightly underestimated (∼0.5 dex) at the high-mass end at z ∼ 4 but overestimated by ∼0.15 dex at the mass interval of $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 10.0\mbox{--}11.0$. The above gives the false impression that the scatter increases more sharply with mass.
• 2.  
  The dark green dashed line represents the σ_sSFR–M_⋆ relation when SFRs are obtained using the FUV luminosity (Kennicutt & Evans 2012) and the IRX–β relation (Meurer et al. 1999; Bouwens et al. 2012; Katsianis et al. 2017b) while the stellar masses are calculated through the SED fitting technique. Using this combination to obtain properties is widely done in the literature, in order to estimate the SFR–M_⋆ relation and its scatter (Santini et al. 2017). We see that this method implies an artificially higher σ_sSFR (with respect the intrinsic black solid line) at the low-mass end. This gives an artificial mass-independent σ_sSFR–M_⋆ relation with a scatter of ∼0.35 dex, which does not evolve significantly from z ∼ 2 to z ∼ 0.
• 3.  
  The red solid line is the σ_sSFR–M_⋆ relation retrieved when stellar masses are calculated via SED fitting and SFRs by combining the FUV and Infrared luminosity estimated from the 24 μm luminosity (Dale & Helou 2002; Kennicutt & Evans 2012). Combining IR and UV luminosities to obtain SFRs in observations is a classic method used in the literature (Daddi et al. 2007a; Santini et al. 2009; Heinis et al. 2014) to obtain the SFR–M_⋆ relation. At z ∼ 4, the scatter is underestimated, with respect the reference black line, by 0.2 dex at high masses. This implies a dispersion with a constant scatter around 0.15 dex. At lower redshifts, there is a good agreement (within 0.05 dex) with the reference intrinsic black line.
• 4.  
  The blue dotted line describes the σ_sSFR–M_⋆ 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 in the case of redshift 4, where the derived relation is mass-independent with a scatter of 0.2 dex. The scatter of the magenta line, which represents the results from SED fitting, is similarly overestimated (by 0.15 dex) with respect to the reference black line, for high-mass objects at the mass interval of $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 10.0\mbox{--}11.0$ at z ∼ 0.9.

Zoom In Zoom Out Reset image size

In conclusion, according to the EAGLE+SKIRT data, the inferred shape and normalization of the σ_sSFR–M_⋆ relation can be affected by the methodology used to derive SFRs and stellar masses in observations. This can affect 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 gives a σ_sSFR–M_⋆ relation close to the intrinsic simulated relation, and can successfully probe the shape of the relation for $\mathrm{log}({M}_{\star }/{M}_{\odot })\sim 9.0\mbox{--}11.0$.

The σ_sSFR–M_⋆ 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, if solely relying on observations. In this paper, we have presented the evolution of the intrinsic σ_sSFR–M_⋆ 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. 2017a) for a wide range of SFRs, stellar masses, and redshifts. The investigation is not limited by the shortcomings encountered by galaxy surveys and addresses a range of redshifts and mass intervals. Our main conclusions are summarized as follows:

• 1.  
  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 σ_sSFR–M_⋆ 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 independent of mass/redshift, with a constant scatter σ_sSFR ∼ 0.2–0.3 (Elbaz et al. 2007; Noeske et al. 2007; Rodighiero et al. 2010; Whitaker et al. 2012). We find that the σ_sSFR–M_⋆ relation has a U-shaped 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, because they usually focus on limited mass intervals (Section 3.2). The finding about the U-shaped form of the relation is supported by results relying on the GAMA survey (Davies et al. 2019) at z ∼ 0.
• 2.  
  AGN and SN feedback drive the shape and evolution of the σ_sSFR–M_⋆ relation in the simulations (Section 3.2). Both mechanisms cause a diversity of star formation histories for low-mass (SN feedback) and high-mass galaxies (AGN feedback).
• 3.  
  Mergers do not play a major role in the shape of the σ_sSFR–M_⋆ relation (Section 3.3).
• 4.  
  We employ the EAGLE/SKIRT mock data to investigate how different SFR/M_⋆ diagnostics affect the σ_sSFR–M_⋆ 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–β relation for dust corrections imply a constant scatter with stellar mass with almost no redshift evolution. The methodology used to derive SFRs and stellar masses can affect the inferred σ_sSFR–M_⋆ relation in observations, and thus compromise the robustness of conclusions about its shape and normalization.

We would like to thank the anonymous referee for suggestions and comments that significantly improved 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. is thankful for support 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. M.S. acknowledges support by the Ministry of Education, Science, and Technological Development of the Republic of Serbia through the projects Astrophysical Spectroscopy of Extragalactic Objects (176001) and Gravitation and the Large Scale Structure of the Universe (176003). A.K. 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.

In this appendix, we base our analysis of the dispersion of the sSFRs at different mass intervals on their distribution/histogram, namely the 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 unable to provide quantitative information regarding how galaxies are distributed around the mean sSFR and cannot account for galaxies that could be undersampled or missed by selection effects. In Section 3.1, we present the evolution of the σ_sSFR–M_⋆ relation at z ∼ 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 Figure 9 (z ∼ 0.8–1.4) and Figure 10 (z ∼ 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. 2017a), the evolution of the star formation rate function (Katsianis et al. 2017a), and the rates of star formation and black hole accretion (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 in the methodologies used by observers to obtain the intrinsic SFRs and stellar masses (Katsianis et al. 2017b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the top panels of Figures 9 and 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-exponential) 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 resamples of the EAGLE SFRs, while the black error bars represent the 1σ Poissonian errors. We can see that the agreement between the simulated and observed sSFR functions is typically good for low-mass objects, but breaks down for galaxies more massive than $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 10$ at 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 parameterization. They argued that the fit with a double-exponential function is more suitable than the log-normal function at z ∼ 0. The EAGLE reference model indicates that the sSFRF of 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, which is not commonly used to describe the sSFR distribution (Ilbert et al. 2015), allows a significant density of SFGs with a low sSFR, and the confirmation of this shape in future observations is important.