© ESO 2020

1. Introduction

Galaxy clusters are the most massive objects in the universe; they are located at the nodes of the cosmic web and are characterised by very high densities. Given their extreme nature, they are very important probes for both cosmology (Kravtsov & Borgani 2012) and galaxy formation models, where the evolution of galaxies and thus their resulting properties depend on a range of environmental effects.

In the local universe, galaxy clusters appear as massive virialised objects. They are characterised by the presence of a hot (∼10⁸ K) diffuse plasma, namely the intracluster medium (ICM), which can be studied through X-ray observations and the Sunyaev-Zeldovich effect (Rosati et al. 2002; Carlstrom et al. 2002). Galaxies, whose velocity dispersion is consistent with the ICM temperature, are mainly bulge dominated and ellipticals with reduced ongoing star formation, especially in the cluster core region. The stellar population is typically very old, with stars nearly as old as the Universe. Indeed observations (Mancone et al. 2010; Wylezalek et al. 2014; Foltz et al. 2015) and theoretical semi-analytical models (e.g., De Lucia & Blaizot 2007) suggest that galaxies are undergoing passive evolution since z ∼ 1 and have formation times of z_f ≳ 2.

As opposed to the local universe, where all galaxy clusters show similar properties, at z ≳ 1.4 clusters behave as a rather diverse population. Some observations report the discovery of already mature clusters, with an enhanced fraction of red and quenched galaxies with respect to the field, at least in the core (Papovich et al. 2010; Strazzullo et al. 2010, 2013; Gobat et al. 2011; Tanaka et al. 2013a; Newman et al. 2014; Andreon et al. 2014; Cooke et al. 2016), while others clearly show a mixed population of quenched and star forming galaxies (Tanaka et al. 2013b; Brodwin et al. 2013; Gobat et al. 2013; Hatch et al. 2017; Strazzullo et al. 2016). A few observations also suggest a reverse star formation rate (SFR)–density relation at z ≳ 1.5, where the specific star formation rate (sSFR) increases toward the core of the cluster (Tran et al. 2010; Santos et al. 2014, 2015; Smith et al. 2019).

At even higher redshift, z ≳ 2, systems lack the presence of a massive virialised halo, and instead show multiple halos spread over large scales (Hayashi et al. 2012; Lemaux et al. 2014; Kubo et al. 2015; Casey et al. 2015; Casey 2016). This is in line with theoretical expectations from numerical simulations, which predict a hierarchical formation of massive clusters formed by the assembly of smaller halos that at z ∼ 2 might occupy a region as large as 20 comoving Mpc (cMpc). (Chiang et al. 2013; Muldrew et al. 2015; Contini et al. 2016).

Since at this redshift (proto)clusters are not virialised, it is difficult to detect them with techniques based on the ICM properties. Therefore, different methods have been used, such as galaxy over-densities. However, this approach can bias the results depending on the galaxy properties used for the selection. In this respect, an important population of galaxies are dusty star forming galaxies (DSFGs; see Casey et al. 2014), which are highly star forming and heavily obscured by dust, emitting in the far-infrared (FIR) and sub-millimetre bands. These galaxies represent the strongest starbursts and are expected to be the progenitors of local massive ellipticals (Cimatti et al. 2008; Ricciardelli et al. 2010; Fu et al. 2013; Ivison et al. 2013; Toft et al. 2014; Gómez-Guijarro et al. 2018). They trace the dusty star-forming phase of protoclusters, and their expected short star-burst phase of a few hundred million years (Granato et al. 2004) makes them relatively rare objects in the sky. Despite their relative scarcity, DSFGs have been successfully used to identify dense and highly star forming environments up to redshift z  ∼  4 (Clements et al. 2014; Oteo et al. 2018; Miller et al. 2018) and have been observed in several previously identified high-redshift protoclusters (Chapman & Casey 2009; Dannerbauer et al. 2014; Umehata et al. 2015; Coogan et al. 2018; Lacaille et al. 2019; Smith et al. 2019).

So far, numerical simulations have been unable to reproduce the high SFRs observed in protoclusters characterised by overdensities of DSFGs (Granato et al. 2015), as simulations fail to predict sufficiently high peaks of star formation activity at early epochs. This result adds to the long-standing difficulty for cosmological simulations to reproduce the star formation properties of galaxies, such as the main sequence (MS) of star forming galaxies, around the peak of the cosmic SFR density (Davé et al. 2016, 2019; McCarthy et al. 2017; Donnari et al. 2019). Indeed, Granato et al. (2015) found that the bulk of star formation within the observed putative progenitors of massive galaxy clusters occurred at higher rates and lasted a shorter time than in simulations. This conclusion was based on the observations available at that time which were characterised by low angular resolution and SFR integrated on the Mpc scale. In the last few years, with instruments like ALMA, it has been possible both to resolve single sources within protoclusters and to obtain information on the galaxy cold gas content (Wang et al. 2018; Gómez-Guijarro et al. 2019; Hill et al. 2020). On the simulations side, progress has been made to increase numerical resolution, needed to resolve higher density peaks and related higher SFRs. Therefore, it is now possible to carry out a deeper inspection of the capability of simulations to reproduce protocluster properties.

In this work we make use of 12 simulations out of the set of 29 hydrodynamical zoom-in simulations of galaxy clusters named Dianoga to investigate the predictive power of state-of-the-art cosmological simulations around the peak of the SFR in the protocluster stage of structure formation. In particular, we aim to compare both the integrated values of SFRs in protocluster regions and the protocluster galaxy properties to recent observations that are now available. The set of simulations used is particularly suited to this aim as it comprises massive objects that can only be found in suitable numbers within large cosmological boxes of ∼1 h⁻¹ Gpc (per side). Moreover, simulations were carried out at a ten-times-higher resolution than before (i.e. Granato et al. 2015).

The paper is structured as follows: in Sect. 2 we describe the simulations setup, with particular focus on the AGN feedback implementation and the subgrid star formation model. In Sect. 3 we show brightest cluster galaxy (BCG) properties and stellar mass function at z = 0. In Sects. 4 and 5 we compare the predicted SFRs in protocluster regions with the available observations at z ∼ 2 and z ∼ 4, respectively, and we analyse the MS of star forming galaxies at both redshifts. In Sect. 6 we show the evolution of the mass-normalised SFR in clusters and protoclusters. In Sect. 7 we study the gas-related properties of our simulated galaxies in comparison with observations. In Sect. 8 we summarise the main results of our analysis and present our main conclusions.

2. Simulations

In this section we describe the set of numerical simulations used in this work, and in particular we detail the observational constraints used to calibrate the subgrid model of AGN feedback. We also briefly review the star formation model of Springel & Hernquist (2003) implemented in our code.

2.1. Setup of simulations

The analysis presented in this paper is based on a set of 12 hydrodynamical zoom-in simulations evolved in a ΛCDM cosmology, with parameters: Ω_m = 0.24, Ω_b = 0.037, n_s = 0.96, σ₈ = 0.8, and H₀ = 100 h km s⁻¹ Mpc⁻¹ = 72 km s⁻¹ Mpc⁻¹. These are part of a sample of 29 simulations, the Dianoga set, described in Rasia et al. (2015), Planelles et al. (2017), Biffi et al. (2017, 2018), Ragone-Figueroa et al. (2018), and Bassini et al. (2019), but have a ten-times-higher mass resolution and small differences in the code (see below) with respect to the cited works. We refer to the previous set of simulations as low-resolution (LR) simulations throughout the paper. The regions are extracted from a parent dark-matter(DM)-only simulation of 1 h⁻¹ Gpc (per side). From this cosmological box the 24 most massive clusters (M₂₀₀ >  8 × 10¹⁴ h⁻¹ M_⊙)¹ were selected together with five randomly chosen, smaller objects (M₂₀₀ ∈ [1 − 4]×10¹⁴ h⁻¹ M_⊙). Their Lagrangian regions, of radius about five times the virial radius of the selected clusters, were re-simulated with the inclusion of baryons and at a greater resolution (see Bonafede et al. 2011 for a full description of the re-simulation procedure). The 12 simulations used for this work include the 5 less massive clusters and 7 massive clusters. The masses of the particles in the high-resolution region are m_DM = 8.44 × 10⁷ h⁻¹ M_⊙ for DM and m_gas = 1.56 × 10⁷ h⁻¹ M_⊙ for the initial gas particles. The Plummer equivalent gravitational softening adopted for DM particles is 4.2 h⁻¹ comoving kpc (ckpc) at z >  2 and 1.4 h⁻¹ physical kpc (pkpc) otherwise. The softening lengths for gas, star, and black hole (BH) particles are 1.4, 0.35, and 0.35 h⁻¹ pkpc, respectively.

The simulations are carried out with the code GADGET-3, a modified version of the Tree-PM smoothed-particle hydrodynamics (SPH) public code GADGET2 (Springel 2005). We employed the hydrodynamical scheme presented in Beck et al. (2016), where a higher order interpolating kernel function is implemented, along with a time-dependent artificial viscosity and a time-dependent artificial conduction, which improve the SPH performance in capturing discontinuities and the development of gas-dynamical instabilities. The unresolved baryonic physics, mostly involving the stellar component and the activity of the BH population, is treated with subgrid models. The prescription of metal-dependent radiative cooling follows Wiersma et al. (2009). The model of star formation and associated feedback is implemented according to the original model by Springel & Hernquist (2003); see below for extra details. For the metal enrichment and chemical evolution we follow the formulation by Tornatore et al. (2007). The stellar yields and a deeper description of metals are specified in Biffi et al. (2018, see also Planelles et al. 2017; Biffi et al. 2017). The treatment of BHs and associated AGNs is detailed below (Sect. 2.3) after a brief summary of the star formation model.

2.2. Star formation

Here we review the main aspects of the subgrid model for star formation which we discuss throughout the paper. For a full explanation we refer to the original paper by Springel & Hernquist (2003). In this model, each SPH particle samples a region of the interstellar medium (ISM), and is subdivided into a cold phase and a hot phase, characterised by densities ρ_c and ρ_h, in pressure equilibrium one with each other. The total density associated to a particle is the sum of the two: ρ = ρ_c + ρ_h. A SPH particle has a non-null fraction of cold gas (i.e. ρ_c >  0), and thus becomes multiphase whenever its density is higher than a given threshold ρ_thr. Given the cold fraction, a numerical instantaneous SFR is associated to each multiphase particle²:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\rho }_{\star }}{\mathrm{d}t}={\dot{\rho }}_{\star }=(1-\beta )\frac{{\rho }_{\mathrm{c}}}{t_{\star }},\end{array}$$(1)

where t_⋆ is the characteristic timescale for star formation, while β is the fraction of massive stars that are expected to instantly explode as supernovae and depends on the chosen IMF. In this work, we employ a Chabrier initial mass function (IMF, Chabrier 2003). The parameter t_⋆ follows the expression:

$$\begin{array}{c}\hfill t_{\star }(\rho )=t_0^{\star }{\left(\frac{\rho }{{\rho }_{\mathrm{thr}}}\right)}^{-1/2},\end{array}$$(2)

and $t_0^{\star }$ is set to 1.5 Gyr in order to match the observed Kennicutt relation (Kennicutt 1998). In practice, varying this parameter directly leads to a variation of the numerical star formation efficiency (SFE). Indeed, using Eqs. (1) and (2):

$$\begin{array}{c}\hfill \mathrm{SFE}=\frac{{\dot{\rho }}_{\star }}{{\rho }_{\mathrm{c}}}=\frac{1-\beta }{t_0^{\star }}{\left(\frac{\rho }{{\rho }_{\mathrm{thr}}}\right)}^{1/2}·\end{array}$$(3)

Given the SFR, part of the cold clouds is evaporated by means of supernovae feedback:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\rho }_{\mathrm{c}}}{\mathrm{d}t}=-A\beta \frac{{\rho }_{\mathrm{c}}}{t_{\star }},\end{array}$$(4)

where A is the efficiency of evaporation that determines the efficiency of thermal supernovae feedback and is taken to be a function of the local gas density, A ∝ ρ^−4/5. These equations give rise to the self-regulated cycle of star formation: high cold cloud density leads to a high SFR, that in turn means more feedback and cloud evaporation. When cloud evaporates, the SFR decreases and material is returned to the hot phase, increasing ρ_h. Finally, a higher density means a higher cooling rate, which causes more gas to condense in cold clouds, and the cycle restarts.

2.3. AGN feedback

The AGN feedback is inspired by the original model developed by Springel et al. (2005) and is implemented according to the scheme described in Ragone-Figueroa et al. (2013) with two main modifications. First, in the feeding process we differentiate between hot and cold accretion (see Sect. 2.3.2, Eq. (5)). Second, we do not impose a temperature threshold to define multiphase gas particles and the energy released by AGN feedback is not used to evaporate the cold phase of gas particles. This second modification is motivated by the fact that this implementation results in better agreement between the simulated galaxy stellar mass function (GSMF) and the observed one (see Sect. 3.1), with one side effect being that overly massive BCGs are produced at z = 0 (see Sect. 3.2.1).

2.3.1. Black hole seeding and positioning

Briefly, during run time we identify groups of particles using the Friends of Friends (FoF) algorithm (Huchra & Geller 1982). In practice, two DM particles are considered to be part of the same group if their distance is less than a fixed parameter, referred to as the linking length, which is commonly defined as a fraction of the mean interparticle distance, $\overline{d}$. We fix this parameter to $0.16\times \overline{d}$. Hence, BHs in our simulations are spawned at the centre of each FoF group (defined as the position of the most bound particle) with a seed mass of 5.5 × 10⁵ M_⊙ whenever all of the following conditions are simultaneously fulfilled: (i) the total stellar mass is higher than 2.8 × 10⁹ M_⊙; (ii) the stellar-to-DM mass ratio is higher than 0.05; (iii) the gas mass is equal to or larger than 10% of the stellar mass; and (iv) no other central BH is already present. As the simulation evolves, we avoid the presence of wandering BHs by adopting a different strategy from that of Ragone-Figueroa et al. (2018). In this latter study, the authors pinned the BHs, meaning that they were re-positioned at each time step at the location of the most bound particle of a halo. Here, instead, we assign a large dynamical mass to the BH and use low values of star and BH particle softening lengths. In other words, the BH dynamical mass is fixed equal to the DM particle mass until it outgrows that value and the softening values are four times smaller than before, once re-scaled to the higher resolution. These numerical prescriptions are sufficient to mimic dynamical friction without the necessity to explicitly include a dynamical friction force (Steinborn et al. 2016). Even though this scheme performs overall well at the current numerical resolution, BH centring remains a major challenge for our, and presumably all, numerical simulations, and it is still possible for a BH to move from the centre of a structure. This is particularly problematic in cluster simulations, where the absence of AGN feedback at the centre of the BCG leads to catastrophic cooling, resulting in high BCG mass and SFR of the order of ∼10³ M_⊙ yr⁻¹ at z ∼ 0. Indeed, in 1 out of the 12 simulations used for this work, the proto-BCG looses its central BH at z ∼ 4, and is characterised by an incredibly high SFR at z = 0. Even though it is not an issue for the conclusions of this work (see Sect. 7), the problem needs to be addressed in future simulations.

2.3.2. AGN accretion and feedback

Once BHs are seeded, they grow by two different channels: accretion of the surrounding gas and BH–BH mergers. The former follows the Eddington-limited alpha-enhanced Bondi accretion rate (Bondi 1952) formula:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Bondi},\alpha }=\alpha \frac{4\pi G^2{M}_{\mathrm{BH}}^2\rho }{{(c_{\mathrm{s}}^2+v_{\mathrm{BH}}^2)}^{3/2}},\end{array}$$(5)

with α equal to 10 and 100 for hot (T >  5 × 10⁵ K) and cold (T <  5  ×  10⁵ K) gas, respectively (Steinborn et al. 2015). In Eq. (5), all gas-related quantities (sound speed, c_s, bulk gas velocity relative to BH velocity, v_BH, and gas density, ρ) are smoothed over 200 gas particles with a kernel function centred at the position of the BH. Given the gas accretion onto the BH particle, AGN energy feedback is given by

$$\begin{array}{c}\hfill \dot{E}={ϵ}_{\mathrm{r}}{ϵ}_{\mathrm{f}}\dot{M}{c}^2,\end{array}$$(6)

where Ṁ is the minimum between Eq. (5) and the Eddington accretion rate, Ṁ = min(Ṁ_Bondi,α, Ṁ_Eddington), and the energy is distributed and thermally coupled to the nearest 200 gas particles. In Eq. (6), ϵ_r is the fraction of mass transformed into radiation energy and ϵ_f is the fraction of radiated energy thermally coupled to the gas particles.

In the previous version of the code the energy was used to evaporate the cold fraction of the multiphase gas particles (see Sect. 2.2 for a brief explanation of multiphase particles; for a comprehensive review we refer to the original paper of Springel & Hernquist 2003), while in the current setup we couple the energy only to the hot phase of each gas particle. The effects of this choice on our results are presented in Sect. 7.3.

Black hole particles can also grow by BH–BH mergers. In our implementation of the subgrid model, two BHs are allowed to merge whenever all the following conditions are fulfilled: (i) v_rel <  0.5 × c_s; (ii) r_rel <  3.5 h⁻¹ pkpc; and (iii) $\left|V_{\text{pot},\text{rel}}\right|+v_{\text{rel}}^2<0.5\times c_{\text{s}}^2$, where v_rel and r_rel are the relative velocity and position between the two BH particles, c_s is the sound speed, and V_pot, rel is the difference between the gravitational potentials computed at the positions of the BH particles.

2.3.3. AGN feedback calibration

The values of ϵ_f and ϵ_r are chosen in order to reproduce the observed correlation between BH mass and galaxy stellar mass (Magorrian et al. 1998). In particular, we aim to reach agreement with the observational results reported by McConnell & Ma (2013) and more recently by Gaspari et al. (2019). The value chosen for ϵ_r is 0.07 independently of Ṁ_Bondi,α, while ϵ_f is lower than 1 only in quasar mode when Ṁ_Bondi,α/Ṁ_Eddington > 0.01 (ϵ_f = 0.15). In Fig. 1 we show numerical results in comparison with observations. In the plot, dark-blue squares represent central galaxies, defined as galaxies at the centre of groups with at least 100 substructures. All other simulated galaxies are represented by light-blue points. Stellar masses of non-central galaxies are computed considering the star particles associated to the substructures identified by the code Subfind (Dolag et al. 2009) and within a sphere of 50 pkpc. Stellar masses of central galaxies are computed considering an aperture of 0.15 × R₅₀₀ to match the aperture used by Gaspari et al. (2019). We note that even though simulations correctly reproduce the normalisation of the observed correlation, the scatter is still under-reproduced, especially at the high mass end. Indeed, the intrinsic scatter around the observed correlation of Gaspari et al. (2019) is σ = 0.40 ± 0.03, while it is a factor of two lower in our simulations (σ = 0.20)³. Even though the observed scatter can be marginally boosted by uncertainties on the assumptions made to obtain these quantities from observational data (e.g. star formation history, IMF, metallicity, etc.), the most probable explanation for this difference is that the subgrid models adopted do not capture the diversity of conditions of BH accretion and AGN feedback at small scales.

Fig. 1. Correlation between the galaxies stellar mass and the central SMBHs mass. Observational data are taken from McConnell & Ma (2013) (dashed black line) and from Gaspari et al. (2019) (red circles). The simulated stellar masses for satellite galaxies (cyan points) are obtained considering the star particles bound to the substructure (accordingly to Subfind) and within 50 pkpc from its centre. The mass of the central galaxies (dark-blue squares) is obtained by summing over all stellar particles within an aperture of 0.15 × R500.

3. Galaxy cluster population at z = 0

In this section we analyse GSMF and BCG properties at z = 0 in our simulations in comparison with observations.

3.1. Galaxy stellar mass function

We start by comparing the observed and simulated GSMFs. First we need to take into account the fact that our simulations are centred on galaxy clusters, while we compare to data obtained for the field. For this reason, we expect to have in simulations a significant higher normalisation of the GSMF. For the purpose of this study, we define galaxy clusters as spherical regions enclosing a mean matter density Δ times the critical density, ρ_crit. Using the relation between ρ_crit and the mean cosmic matter density $\overline{\rho }$:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{crit}}(z)=\overline{\rho }(z)\times [{\mathrm{\Omega }}_{\mathrm{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}]/{\mathrm{\Omega }}_{\mathrm{M}}{(1+z)}^3,\end{array}$$(7)

and assuming that galaxies follow the DM distribution, we have to normalise our GSMF by

$$\begin{array}{c}\hfill N_{\mathrm{norm}}=\mathrm{\Delta }\times [{\mathrm{\Omega }}_{\mathrm{M}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}]/{\mathrm{\Omega }}_{\mathrm{M}}{(1+z)}^3;\end{array}$$(8)

see Vulcani et al. (2014), Appendix B. In practice, we consider the most massive cluster in each of our simulated regions and all the galaxies but the BCG within R₁₀₀ (i.e. Δ = 100), and then we normalise each mass bin by N_norm.

The stellar mass of each galaxy is computed using three different definitions: the sum of all stellar particles that Subfind associates to a substructure and the same sum limited to stellar particles within 30 pkpc and 50 pkpc from the centre. The results are shown in Fig. 2, where we plot the GSMF obtained with the three definitions of stellar mass. Results from simulations agree overall well with observations from Bernardi et al. (2013) starting from the stellar mass of galaxies that we consider well-resolved in our simulated set (M_⋆ >  10¹⁰ M_⊙). The main difference is around 10¹⁰ M_⊙, where simulations have too many galaxies. As noted by Henden et al. (2020), larger values of the softening length numerically decrease the number of galaxies at these masses. However, an investigation of the stellar feedback will be needed to better describe the GSMF at our low-mass end. Regarding the different definitions of stellar mass, there is no statistical difference in the results once a fixed aperture is used, as Subfind likely also associates a fraction of the intracluster light (ICL) to massive galaxies. The good agreement with observational data is an improvement with respect to the GSMF presented in Bassini et al. (2019) and obtained with the set of simulations at lower resolution. In that case, we showed that the GSMF was a factor of approximately two below the results of Bernardi et al. (2013) at M_⋆ ∼ 10¹¹ M_⊙. This difference, as we discuss in the following section, is not due to the increased resolution but to the different implementation of the AGN feedback discussed in Sect. 2.3.

Fig. 2. GSMF at z = 0. Observational data are taken from Bernardi et al. (2013) (black solid line). Simulation data are derived considering as stellar mass the sum of all stellar particles bound to the galaxy by Subfind (red triangles), and the same sum restricted to particles within 50 pkpc (green hexagon) and 30 pkpc (blue squares). Error bars are computed assuming Poissonian errors. The simulated GSMF is normalised following Eq. (8). Filled and empty marks represent the mass bins with respectively more than and less than ten galaxies.

3.2. Properties of the BCG

In this section we study the properties of our BCGs at z = 0. In particular we study their mass and their SFR.

3.2.1. M_⋆,BCG − M₅₀₀ correlation

In Fig. 3 we show the correlation between the stellar mass of the BCG and M₅₀₀. In this plot, M_⋆,BCG is defined as the total stellar mass associated to the halo (according to Subfind) and within a 2D aperture of 50 pkpc. We checked that using different line-of-sight directions brings no more than 40% difference in the estimated stellar mass, with a median difference of 2% for our sample. We note that the total mass represents the BCG stellar mass plus the ICL present along the line of sight. In the plot, we also highlights the properties of the BCG that lost its central BH at z ∼ 4 (see Sect. 2.3.1).

Fig. 3. Correlation between BCG stellar mass and M500 at z = 0. Observations are taken from DeMaio et al. (2018) (blacks quares) and Kravtsov et al. (2018) (black triangles). The simulated values are shown as blue points. The red hexagon refers to the BCG that lost its central BH (see Sect. 2.3.1). The orange line is the fit to LR simulations (Ragone-Figueroa et al. 2018). BCG masses are obtained summing over all stellar particles bound to the main subhalo of a group or cluster by Subfind (BCG+ICL) and within a 2D aperture of 50 pkpc.

From the figure, we see that our simulations tend to have overly massive BCGs with respect to observational data (a factor of two at M₅₀₀ = 3 × 10¹⁴ M_⊙), in line with the findings of other groups (e.g. Bahé et al. 2017; Pillepich et al. 2018; Henden et al. 2020). Moreover, current simulations have more massive BCGs than the LR simulations from Ragone-Figueroa et al. (2018, orange line in the plot), being a factor of 2.3 more massive at M₅₀₀ = 3 × 10¹⁴ M_⊙. A possible explanation could be the different resolution. However, in order to investigate this hypothesis we ran a simulation at a ten-times-lower mass resolution and obtained only a ∼12% lower BCG mass. Hence the BCG mass is very stable against the mass resolution of the simulation (see also Ragone-Figueroa et al. 2018).

To check whether or not the different results with respect to LR simulations (i.e., Ragone-Figueroa et al. 2018) are due to the different implementations of the AGN feedback adopted, we ran two more simulations: one at the current mass resolution and the other with a ten-times-lower resolution, both implementing the same AGN prescription as that in Ragone-Figueroa et al. (2018). In this setup, gas particles need to be colder than a fixed temperature threshold to be considered multiphase and the energy released by the AGN feedback is used to evaporate the cold gas. Even in this case, we do not find any trend with resolution, the results being in agreement within 5%. Moreover, the BCG masses obtained in these last two runs are in agreement within 30% with Ragone-Figueroa et al. (2018) results (orange line in Fig. 3), pointing out that the BCG mass is most sensitive to the prescription adopted for AGN feedback.

Interestingly, with the Ragone-Figueroa et al. (2018) setup, together with a lower BCG mass, we also get a lower normalisation for the GSMF, in line with the results of Bassini et al. (2019). Therefore, we conclude that the different BCG masses are not due to the increased resolution but to the different prescription for the AGN feedback and that, with the current implementation of feedback, it is difficult to simultaneously reproduce both the observed M_BCG − M₅₀₀ relation and the GSMF.

3.2.2. Star formation rate of BGCs

In Fig. 4 we show the SFR of our simulated BCGs in comparison with observational data. Observations are taken from McDonald et al. (2018) and constitute a subsample of the BCGs used in Fraser-McKelvie et al. (2014). The original selection has been made considering all galaxy clusters in a volume-limited sample, z <  0.1, with a measured X-ray luminosity in the ROSAT 0.1 − 2.4 keV band L_X >  10⁴⁴ erg s⁻¹. This luminosity cut ensures a completeness > 80% for the cluster sample (see Fraser-McKelvie et al. 2014). Moreover, selecting the sample on cluster properties rather than BCG properties enables us to correctly account for low values of SFR and non-detections. However, it has been noted that the SFRs published in Fraser-McKelvie et al. (2014) lack important k-corrections which lead to biased results (see Green et al. 2016). McDonald et al. (2018) recomputed the SFRs using 12 μm flux following the procedure of Green et al. (2016) for all the clusters with L_X >  3.3 × 10⁴⁴ erg s⁻¹. The final sample comprises 33 objects and is complete above the cut in X-ray luminosity. In grey we show the results of our simulations. We emphasise that the BCGs are taken from the same simulated regions at different redshifts, and therefore are not independent. In blue we plot the median value with 16th and 84th percentiles. To mimic the selection of McDonald et al. (2018), we considered only the 11 clusters with M₅₀₀ >  2.8  ×  10¹⁴ M_⊙ (M₂₀₀ ≳ 4  ×  10¹⁴ M_⊙), which corresponds to L_X >  3.3  ×  10⁴⁴ erg s⁻¹ following the correlation between L_X and M₅₀₀ showed by Truong et al. (2018). The simulations used in the work of Truong et al. (2018) are not the same as those used for this work, but we checked that using a relation based on our clusters leads to the same final sample. The SFR in simulations is the instantaneous SFR predicted by the effective model for multiphase particles computed considering all the particles bound to the group by Subfind and within a 2D aperture of 30 pkpc. We employed this aperture to directly compare with other numerical simulations, after checking that the aperture choice does not affect our conclusions. Our simulations present a high residual SFR at these low redshifts (∼10 M_⊙ yr⁻¹ against the average observed SFR of ∼0.3 M_⊙ yr⁻¹). Considering a 3D aperture of 30 pkpc brings us to the same conclusions, as the bulk of SFR is located near the centre of the cluster and the median difference between a 3D and a 2D aperture is 25%. Similar results were also found by other groups. Henden et al. (2020) made a similar analysis considering all the clusters with M₂₀₀ >  10¹⁴ M_⊙ and computed the instantaneous SFR within a spherical aperture of 30 kpc at z = 0.2. Their results show that their BCGs form stars at a rate similar to ours when they do not have a null SFR (see their Fig. 8).

Fig. 4. Star formation rate of BCGs in observations and simulations. Grey circles are BCGs of our simulations from different snapshots, while the grey triangle represents the BCG that lost its central BH at z ∼ 4 (see Sect. 2.3.1). BCGs from the same snapshot are shifted only for visualisation purposes. The median values are shown as blue circles and the vertical bars indicate the range between the 16th and 84th percentiles. A 2D aperture of 30 pkpc is used. Red squares are BCGs from the sample of McDonald et al. (2018) (see text for more details).

Even though the disagreement between simulations and observations is quite large (a factor of ∼30 at z ∼ 0 according to our results), it is important to keep in mind that measuring the SFR of galaxies is always a non-trivial task, especially in the case of BCGs due to their low SFR values and to the crowded environment. In the particular case of the sample used in our comparison, the SFR is obtained from the 12 μm luminosity through the relation derived by Cluver et al. (2014). However, this relation is calibrated on star-forming galaxies and flattens at SFR < 5 M_⊙ yr⁻¹. For this reason the values of SFR inferred from observations of BCGs are likely to be underestimated, thus possibly alleviating the tension outlined in Fig. 4. The same caution should be exercised when considering the conclusions drawn in the following paragraph.

In Fig. 5 we plot the sSFR for our simulations and the McDonald et al. (2018) observations. As we did for the SFR, we checked that the choice of the aperture does not affect our results. Indeed, the results obtained using an aperture of 30 and 50 pkpc are in agreement within 30% at z = 0. Also in this case, numerical simulations appear to be an order of magnitude above observations. Similar results were also found by other groups. Davies et al. (2019) showed that Eagle simulation presents a sSFR ∼10⁻¹¹ yr⁻¹ at M₂₀₀ ∼ 10¹⁴ M_⊙ and z = 0, and that IllustrisTNG BCGs have a sSFR of ∼2.5 × 10⁻¹² at M₂₀₀ ∼ 10¹⁴ M_⊙, which is a factor of three higher than the median value of the McDonald et al. (2018) BCGs. Moreover, considering that the BCG sSFR is an increasing function of mass in IllustrisTNG (see Fig. 12 of Davies et al. 2019) and that the sample of McDonald et al. (2018) is of massive clusters (L_X >  3.3 × 10⁴⁴ erg s⁻¹), the factor of approximately three is probably a lower limit.

Fig. 5. Specific SFR of BCGs in observations and simulations. Grey circles are BCGs of our simulations from different snapshots (blue circles are median values with 16th and 84th percentiles), while the grey triangle is used for the BCG that lost its central BH at z ∼ 4 (see Sect. 2.3.1). BCGs from the same snapshot are shifted only for visualisation purposes. A 2D aperture of 30 pkpc is used. Red squares are BCGs from the sample of McDonald et al. (2018) (see text for more details).

For our simulations, a possible solution to this mismatch could be a more effective AGN feedback. However, Ragone-Figueroa et al. (2018) found a very similar result (sSFR ∼ 1.5 × 10⁻¹¹ yr⁻¹ at z = 0) with an implementation of the AGN feedback that is more effective in quenching star formation (see also results in Sects. 3.1 and 3.2.1). Moreover, with the current scheme a more effective feedback could (at least in principle) reduce the gap in the SFR between simulated and observed BCGs, at the cost of an increase in the difference in the GSMF (see Sect. 3.1). Therefore, a better solution would be to modify the prescription of the AGN feedback in order to have more efficient quenching only for massive galaxies.

4. Protoclusters at z ∼ 2

In this section we compare our simulations to observational results of protocluster regions of different sizes and identified at z  ∼  2. As we are interested in comparing values of SFR, we include only protocluster regions with coverage in the FIR and submillimetre bands, because the total star formation budget is mainly contributed by its obscured component. In particular, we compare with the observations of Clements et al. (2014), Dannerbauer et al. (2014), Wang et al. (2016), Kato et al. (2016), Coogan et al. (2018), and Gómez-Guijarro et al. (2019).

4.1. Protocluster SFR within ∼1 pMpc

In Fig. 6 we compare the observed star formation obtained by Clements et al. (2014), Dannerbauer et al. (2014), and Kato et al. (2016) with the same quantity computed in simulations. With the red bands we show the protocluster regions at 1 ≲ z ≲ 3 identified as clumps in Planck 857 GHz band by Clements et al. (2014). This frequency is suitable for identifying DSFGs, which trace star-bursting phases of protoclusters (e.g. Granato et al. 2004). Clements et al. (2014) retrieve the far-infrared luminosity, L_FIR, by fitting the spectral energy distribution (SED) with a modified black body formula, and then they compute the SFR using the relation given by Bell (2003), which assumes a Salpeter IMF (Salpeter 1955). In order to compare with our simulations, which instead adopt a Chabrier IMF (Chabrier 2003), we divide the Clements et al. (2014) values of SFR by 1.74. We also include in our comparison the analysis that these latter authors performed on the fields previously introduced by Stevens et al. (2010), who analysed the density flux of submillimetre galaxies (SMGs) obtained at 850 μm in five fields centred on QSOs in the redshift range 1.7 <  z <  2.8. The SFRs for these fields are computed following the procedure already described, with the only difference being that the FIR luminosity is computed from the F₈₅₀ flux using an Arp 220 spectral template. Also in this case we corrected for the choice of the IMF. Among the four clumps and five fields analysed in Clements et al. (2014) we include in our comparison the six residing at 1.74 <  z <  2.27 (the values of SFR for the four fields are 2240, 2330, 3966, and 4095 M_⊙ yr⁻¹, and thus they overlap in Fig. 6). The physical volume used to compute the SFRs in Clements et al. (2014) clumps is 4.2 Mpc³, equal to a sphere of 1 Mpc radius, while the fields are characterised by a volume of 1.4 Mpc³, equivalent to a sphere of 0.7 Mpc radius.

Fig. 6. Star formation rate of protocluster regions at z ∼ 2 in observations and simulations within an aperture of ∼2 pMpc. Red bands refer to two clumps from Clements et al. (2014), black solid lines refer to four fields from Stevens et al. (2010) and analysed by Clements et al. (2014). The blue square highlights to the Spiderweb structure (Dannerbauer et al. 2014). The green square and green band show the two protoclusters analysed by Kato et al. (2016) (HS1700 and 2QZCluster, respectively). Black circles and triangles refer to numerical simulations, where the SFR is plotted against protocluster mass (see text). We used black circles for groups which end up in the central cluster of the region at z = 0, and black triangles otherwise.

Dannerbauer et al. (2014) studied the FIR properties of the protocluster associated to the radio-galaxy HzRG MRC1138-262 at z  =  2.16 (also known as Spiderweb galaxy, Pentericci et al. 2000; Miley et al. 2006). This structure is characterised by an overdensity of Lyman alpha emitters (LAEs), and has been studied at different submillimetre wavelengths (see Dannerbauer et al. 2014, and references there in). Observations in the FIR (100, 160, 250, 350, 500, and 850 μm) were used to fit the SED of detected sources assuming a grey body formula, which was used to compute the total FIR luminosity and the correlated SFR through the relation given by Kennicutt (1998). The resulting SFR computed within a sphere of 1 pMpc radius, corrected to a Chabrier IMF, is ∼3600 M_⊙ yr⁻¹ and is showed as a blue square in Fig. 6. Numerical simulations suggest that this protocluster is the progenitor of a massive z = 0 cluster (Saro et al. 2009), with a predicted mass of M₂₀₀ ∼ 10¹⁵ M_⊙ yr⁻¹. Therefore, this structure is a candidate progenitor of the massive simulated clusters used in this study. Finally, we also show the observations in the FIR of another two known protoclusters: 2QZCluster (z = 2.23, Matsuda et al. 2011) and HS1700 (z = 2.3, Steidel et al. 2005). Kato et al. (2016) used SPIRE bands (250, 350, and 500 μm) to obtain a colour-selected sample of DSFGs possibly associated to these two protoclusters. In their work, these latter authors found overdensities of DSFGs in both protoclusters regions, even though the redshift of these sources is not yet confirmed. Assuming a grey-body spectrum, these latter authors computed L_IR in the ∼1 pMpc region containing the highest number of DSFGs, obtaining values very similar to the one reported by Clements et al. (2014). We show Kato et al. (2016) results in Fig. 6 as a green square (HS1700) and green band (2QZCluster). In both cases the upper limit on the value of SFR is obtained assuming that all the detected sources are within the protocluster, while the lower limit is obtained subtracting field average values (see Kato et al. 2016 for further details).

Regarding the data from simulations, we considered at z = 2 the five most massive groups identified by Subfind in each of the analysed regions. The mass M₅₀₀ of the group is given by Subfind, while the SFR is the sum of the instantaneous SFR of all gas particles within a sphere of 1 Mpc radius from the centre of the group. This aperture matches the volume adopted for Clements et al. (2014) clumps, and is slightly larger than the volume of Stevens et al. (2010) fields. In addition, we adopt another possible definition of the SFR: the SFR averaged over ∼100 Myr. Although this may be the optimal choice when comparing to DSFGs, it does not quantitatively affect our results as our most highly star forming protocluster region is characterised by a 2% difference in the SFR when the two methods are used. For this reason we only use the instantaneous SFR throughout the paper. In Fig. 6 we also differentiate between the progenitors of the most massive cluster at the centre of each region by z = 0 and the groups that will form other objects. From the plot it is clear that the simulated protoclusters do not reproduce the high SFR observed, as the difference between the highest value of SFR rate reached within our set of simulated protoclusters (∼1300 M_⊙ yr⁻¹) is a factor of approximately five lower than the SFR measured in one of the Clements et al. (2014) clumps (∼7000 M_⊙ yr⁻¹), and a factor of approximately three lower than the SFR measured within the protocluster associated to the Spiderweb galaxy (∼3600 M_⊙ yr⁻¹). This result is in agreement with the conclusions of Granato et al. (2015), who used a set of simulations with the same initial conditions used here but at ten-times-lower mass resolution and a previous version of our code (the main difference being the prescription for AGN feedback and BH repositioning; see Ragone-Figueroa et al. 2013) together with dust reprocessing and radiative transfer post-processing performed with GRASIL-3D (Domínguez-Tenreiro et al. 2014) to directly compare FIR fluxes with Clements et al. (2014). Granato et al. (2015) concluded that simulations fail to reproduce the observed fluxes at z = 2 by a factor ≳3 − 4. Given the results shown in Fig. 6, we conclude that the results of Granato et al. (2015) hold at higher resolution and are not dependent on the particular prescription adopted for AGN feedback.

4.2. Protocluster SFR within ∼100 pkpc

In Fig. 7 we compare simulations with observations by Wang et al. (2016) (blue square), Coogan et al. (2018) (red square), and Gómez-Guijarro et al. (2019) (green bands).

Fig. 7. Star formation rate of protocluster regions at 2 <  z <  2.6 in observations and simulations within an aperture of ∼100 pkpc. Green bands refer to two protoclusters from Gómez-Guijarro et al. (2019), the blue square refers to Wang et al. (2016), and the red square refers to Coogan et al. (2018). Black circles and triangles refer to numerical simulations, where the SFR is plotted against protocluster core mass. We used black circles for groups which end up in the central cluster of the region at z = 0, and black triangles otherwise.

In particular, Wang et al. (2016) recently discovered a cluster (CL J1001+0220) at z = 2.506. This structure, detected as an overdensity of distant red galaxies (DRGs), appears as a massive, virialised halo. Through an analysis of the velocity dispersion, stellar mass content, and also the detected X-ray emission, Wang et al. (2016) estimated the cluster mass to be 10^{13.9 ± 0.2} M_⊙. However, differently from local clusters, CL J1001+0220 is characterised by a high fraction of massive (M_⋆ >  10¹¹ M_⊙) star forming galaxies. The SFR in the 80 pkpc core region, computed from FIR luminosity and corrected to a Chabrier IMF, is estimated to be ∼2000 M_⊙ yr⁻¹. Moreover, the fraction of starbursting galaxies is ∼25%, much higher than the value in the field which is about 3% accordingly to Schreiber et al. (2015). Gómez-Guijarro et al. (2019) spectroscopically confirmed two protoclusters through CO emission lines (a third one was analysed in the same work, but it was associated with the well-known CL J1001+0220 cluster in the COSMOS field; see Wang et al. 2016, 2018). These objects were previously recognised as separate sources by Bussmann et al. (2015) who observed with ALMA 29 bright DSFGs taken from the Hermes Survey (Oliver et al. 2012). The two new protoclusters are composed of four and five gas rich DSFGs over a region of 125 pkpc and 64 pkpc at z = 2.171 and z = 2.602, respectively. The L_IR used to compute the SFRs of single sources following Kennicutt (1998) are computed integrating the SED fitted from the available IR flux measurements at 24, 250, 350, 500 μm, and 3 mm.

Finally, we also include in the comparison the observations by Coogan et al. (2018) of the protocluster Cl J1449+0856, identified by Gobat et al. (2011) as an overdensity of IRAC colour-selected galaxies. The protocluster has also been detected from the X-ray emission, from which a mass has been estimated in the range [4 − 6]×10¹³ M_⊙ (Valentino et al. 2016). Coogan et al. (2018) employed ALMA observations of 870 μm continuum and the CO(4–3) emission line to compute the SFR within the ∼0.08 pMpc² cluster central region (see Coogan et al. 2018 or Strazzullo et al. 2018 for more details on the computation of the SFR). The value of SFR reported in Fig. 7 is the sum of the SFR obtained for obscured (∼400 M_⊙ yr⁻¹) and unobscured (∼60 M_⊙ yr⁻¹) star formation. As a final remark, we note that the values of SFR computed through the 870 μm continuum rely on the assumed SED template. The values reported in Fig. 7 are obtained considering a template for MS galaxies. Considering a template typical of starburst galaxies, the value of SFR would be twice as high (e.g. Strazzullo et al. 2018).

In numerical simulations we computed the instantaneous SFR considering a 2D aperture of 90 pkpc (around the mean value of the four observed protoclusters that we use as comparison) and integrating 1 pMpc along the line of sight. The choice of the projected distance does not affect our results as most of the stars are produced at the centre of the protoclusters. Indeed, we verified that the results are quantitatively the same integrating up to 3 pMpc along the line of sight. Similarly to previous results, the most highly star forming region within our set of simulations (SFR ∼ 500 M_⊙ yr⁻¹) underpredicts the highest observed SFR (SFR ∼ 2000 M_⊙ yr⁻¹) by a factor of approximately four.

4.3. Main sequence of star forming galaxies

To explore a possible origin for the difference between SFRs in simulations and observations, in Fig. 8 we show the observed and simulated correlation between stellar mass and SFR in galaxies. In particular, we compare with the MS of star forming field galaxies as derived by Whitaker et al. (2014), obtained considering star forming galaxies selected in UVJ colours in the redshift range 2 <  z <  2.5, and with the MS at z ∼ 2.3 from Schreiber et al. (2015), who used a similar approach but considered only photometry at rest-frame wavelengths larger than 30 μm to avoid pollution from AGNs. To compare these latter results with our simulations and the results of Whitaker et al. (2014) we corrected the Schreiber et al. (2015) MS to a Chabrier IMF. In the plot we also include single galaxies from the protoclusters of Gómez-Guijarro et al. (2019) and the cluster of Wang et al. (2018). For our simulations, we use the redshift-dependent threshold in sSFR of Pacifici et al. (2016), which mimics the selection in UVJ colours, to distinguish quiescent from star forming galaxies.

Fig. 8. Star formation rate as a function of galaxy stellar mass at z ∼ 2.3. Red solid and dashed lines are observational data from Whitaker et al. (2014) and Schreiber et al. (2015), respectively. Green hexagons and blue squares are galaxies from the protoclusters of Gómez-Guijarro et al. (2019) and the cluster of Wang et al. (2018), respectively. Grey points are galaxies in our simulations. Black dashed line fix the distinction between active and passive galaxies (Pacifici et al. 2016). Black points represent median values of star forming galaxies with 16th and 84th percentiles. Both SFRs and stellar masses are computed considering a 3D aperture of 30 pkpc.

As we can see from the plot, the SFRs of simulated galaxies are below the observed relation by a factor of approximately three. This is a known discrepancy between simulations (and also semi-analytical models) and observations (see, e.g. Davé et al. 2016; Xie et al. 2017). Indeed, around the peak of the cosmic SFR density, simulations show a normalisation for MS of star forming galaxies that is typically a factor of approximately two to three lower than observations. In Fig. 9 we show the results from different numerical simulations and semi-analytical models. All but our simulations refer to cosmological boxes, meaning that a possible bias in our results toward a lower MS normalisation can be expected when comparing with other simulations. However, observational works do not find differences between the MS computed in different environments (e.g. Koyama et al. 2013). This result is also in line with the predictions of the GAEA semi-analytical model (Hirschmann et al. 2016).

Fig. 9. Main sequence of star forming galaxies at z ∼ 2. Red triangles are observational data from Whitaker et al. (2014). The black line shows median values for our simulations. Coloured solid and dashed lines are data from other cosmological simulations and semi-analytical models respectively. In particular: Eagle (orange solid line, Guo et al. 2016), TNG300 (red solid line, Donnari et al. 2019), Simba (yellow solid line, Davé et al. 2019), Galform (green dashed line, Guo et al. 2016), L-galaxies (dark green dashed line, Guo et al. 2016), and GAEA (blue dashed line, Hirschmann et al. 2016). For the GAEA model we also show the results obtained considering only galaxies that at z = 0 are within galaxy clusters with mass > 1014.25 M⊙ (see text for more details).

In particular, we computed the MS of star forming galaxies considering only the main progenitors of the galaxies that are found in a cluster with mass M_vir >  10^14.25 M_⊙ at z = 0 (GAEA clusters in the plot), finding no more than a 30% difference with respect to the MS obtained considering all active galaxies in the simulation.

The discrepancy outlined in Fig. 9 is an interesting feature given the fact that this difference persists also in numerical simulations which reproduce the GSMF at every redshift (e.g. Davé et al. 2019). A systematic factor of approximately three in the galaxy SFR at z  ∼  2, which naturally arises in cases where the observed normalisation of the MS is matched, could alleviate the discrepancy between the SFR in simulated and observed protoclusters (see Figs. 6 and 7). Moreover, looking at individual galaxies in observed protocluster regions in Fig. 8, we see that most of them are above the MS with also a few galaxies classified as starburst. In this respect it is interesting to note that galaxies within the cluster identified by Wang et al. (2016), detected also in X-ray and hence probably in a more mature evolutionary stage with respect to the structures identified by Gómez-Guijarro et al. (2019), are characterised by higher masses and are scattered around the observed MS. On the contrary, galaxies within Gómez-Guijarro et al. (2019) structures have lower masses and very high SFRs, all above the MS. The level of SFR of these galaxies is not reproduced by simulations that, besides underpredicting the normalisation of the MS, do not exhibit strong starbursts (see Sect. 7.1).

5. Protoclusters at z ∼ 4

In this section we compare our simulations to observational results of protocluster regions identified at z ∼ 4. In particular we compare with the observations by Oteo et al. (2018) and Miller et al. (2018). Before discussing our results, it is important to make a few considerations. First, the protoclusters studied at z ∼ 2 in the previous section come from relatively small surveys, the largest being the one analysed by Clements et al. (2014). This survey encompasses 90 deg², that in the redshift range 0.76 <  z <  2.3 corresponds to ∼0.6 h⁻³ cGpc³ in our cosmology, and thus is smaller than the cosmological box from which the simulated clusters are extracted (see also Granato et al. 2015). This is not true for the protoclusters studied by Oteo et al. (2018) and Miller et al. (2018). The first has been identified within the H-ATLAS fields, corresponding to a total sky area of ∼600 deg² (Ivison et al. 2016), while the second comes from a catalogue from ∼770 deg² of the South Pole Telescope Sunyaev-Zel’dovich (SPT-SZ) survey (Mocanu et al. 2013). The comoving volume corresponding to the H-ATLAS fields in the redshift range 2.7 <  z <  6.4, corresponding to the redshift spanned by the ultra-red galaxies selected with Herschel by Ivison et al. (2016), is ∼10 h⁻³ cGpc³, about ten times larger than the box from which the simulated clusters are extracted. Therefore, Oteo et al. (2018) and Miller et al. (2018) structures could be sufficiently rare not to be sampled by our simulations. Moreover, the protoclusters analysed at z ∼ 2 include a few bona fide z = 0 massive clusters (Dannerbauer et al. 2014; Wang et al. 2016; Coogan et al. 2018). However, this may not be the case for the ones observed by Oteo et al. (2018) and Miller et al. (2018), as it is not guaranteed that a halo of M_halo ∼ 10¹³ M_⊙ at z ∼ 4 will eventually evolve to a Coma-like structure at z = 0. Indeed, numerical simulations suggest that the value of the mass of the most massive halo in a protocluster region at z ∼ 4 is not enough to safely predict the cluster mass by z = 0. An analysis of the large-scale structure, such as the galaxy overdensity over a scale of ∼5 pMpc, would be needed to place better constraints on the final cluster mass (see Chiang et al. 2013). It is important to keep this in mind when comparing observations and simulations.

The observations by Oteo et al. (2018) and Miller et al. (2018) of two highly star forming protocluster cores at z ∼ 4 and z ∼ 4.3 are shown in Fig. 10 with a green line and blue squares respectively. The protocluster core presented by Oteo et al. (2018) was firstly detected as part of an overdensity of DSFGs in the wide-field LABOCA (a low-resolution bolometer camera on the APEX telescope) map at 870 μm. Subsequent observations with ALMA at 2 mm and 3 mm revealed that the most luminous source consists of at least 11 separate sources, of which ten were spectroscopically confirmed to be at z = 4.002 and are therefore part of the same structure. The L_IR for the ALMA resolved sources are computed considering the flux density at 2 mm (∼400 μm at rest frame) and assuming an ALESS template for the SED. The resulting SFR (corrected to a Chabrier IMF) in the 260 pkpc × 310 pkpc central region is ∼3700 M_⊙ yr⁻¹. We note that this value is highly uncertain as it depends on the assumed SED template. However, Oteo et al. (2018) showed that within a large range of SED templates, only the one reported by Pearson et al. (2013) yields a lower SFR than the one obtained with ALESS (by a factor of 0.66). The protocluster from Miller et al. (2018), SPT2349-56, was first detected by the South Pole Telescope (SPT). Subsequent follow up observations with LABOCA and ALMA led to the identification of 14 sources in an extremely small area (∼130 kpc diameter) at z = 4.3. Star formation rates were derived from 870 μm flux density (S_870 μm) assuming a SFR-to-S_870 μm ratio of 150 ± 50 M_⊙ yr⁻¹ mJy⁻¹, which is typical for SMGs. The gas mass of all 14 galaxies was computed from CO(4–3) line luminosity (converted to CO(1–0) line luminosity through the ratio r_4, 1 = 0.41 ± 0.07) assuming a CO/H₂ conversion factor of ${\alpha }_{\mathrm{CO}}=0.8\frac{M_{\odot }}{\mathrm{K}\mathrm{km}{\mathrm{s}}^{-2}{\mathrm{pc}}^2}$ and through the relation:

$$\begin{array}{c}\hfill M_{\mathrm{gas}}={\alpha }_{\mathrm{CO}}{L}_{\mathrm{CO}(1-0)}^{\prime }.\end{array}$$(9)

Fig. 10. Star formation rate as a function of M500 at z ∼ 4.3. The blue square and green line are the observed values of Miller et al. (2018) and Oteo et al. (2018) protoclusters respectively. Black symbols refer to the SFR computed in a cylinder 1 pMpc long and within a circular aperture of 130 pkpc in our simulations. The five most massive groups of each region are shown. We used black circles for groups which end up in the central cluster of the region at z = 0, and black triangles otherwise.

When the CO(4–3) line is not detected, [CII] line luminosity is converted to CO(4–3) using the average CO(4–3)/[CII] ratio for their detected sample (Miller et al. 2018).

The SFR in simulations is computed considering a 2D aperture of 130 pkpc and integrating 1 pMpc along the line of sight. We also verified that the choice of the length of the cylinder does not quantitatively affect our results. In particular, integrating over the whole box along the line of sight and using different orientations for the cylinder axis induce differences in the measured SFR not higher than 50% at SFR > 400 M_⊙ yr⁻¹. This suggests that at this redshift, the star formation takes place only in the densest simulated regions. As we did at z ∼ 2, for each of our regions we selected the five most massive groups at z = 4.3. Also in this case, the most star forming group within our simulations differs from observations by a factor of approximately four.

However, it is important to note that both SPT2349-56 and the protocluster observed by Oteo et al. (2018) represent very rare objects. In fact, none of our simulated protoclusters at z ∼ 4 has more than seven star forming galaxies with a mass higher than 10¹⁰ M_⊙, while Oteo et al. (2018) and Miller et al. (2018) spectroscopically confirmed 10 and 14 sources, respectively. Thus, we conclude that among the main progenitors of our 12 clusters (7 of which very massive), we do not have a structure with the same number of star forming galaxies as these observed protoclusters. Even though this could certainly be due to the limited statistics of the simulated volumes, it could also be related to the star formation subgrid model, which does not correctly describe galaxy properties at this redshift, or both. Nevertheless, assuming that doubling the number of star forming galaxies within a protocluster to match the observed number within SPT2349-56 would also double the total SFR, the SFR of the ‘boosted’ simulated protoclusters would still be a factor of approximately two lower than the observed SFR.

In Fig. 11 we show the MS of star forming galaxies at 4 <  z <  4.8. In particular we compare our simulations with the MS for the field as observed by Steinhardt et al. (2014), together with detected galaxies in SPT2349-56 as reported by Hill et al. (2020), who updated the values of SFR of Miller et al. (2018) and estimated the mass of single galaxies by dynamical methods (through their measured line widths). As we can see from the plot, simulations show relatively good agreement with observations with no statistical difference in the normalisation of the MS. Therefore, we cannot explain the difference we observe in terms of SFR as being due to a systematic offset on the SFR–M_⋆ plane. However, if we look at the galaxies in SPT2349-56, we see that they are scattered around the MS with also a few strong starbursts. On the contrary, galaxies in our simulated protoclusters are mainly MS galaxies with a very small scatter. Therefore, we conclude that at z ∼ 4 our simulations fail to reproduce the high SFR observed because they are unable to produce strong starbursts lying well above the MS.

Fig. 11. Star formation rate as a function of galaxy stellar mass at z ∼ 4.3. The red line shows observational data from Steinhardt et al. (2014). Orange dots represent galaxies of SPT2349-56 as analysed in Hill et al. (2020). Grey points are galaxies in our simulations. The black dashed line fixes the distinction between quiescent and star forming galaxies (Pacifici et al. 2016). Black points represent median values with 16th and 84th percentiles for star forming galaxies. Green circles are galaxies from the simulated protoclusters shown in Fig. 10.

6. Redshift evolution of mass-normalised SFR

In this section we show the evolution of the SFR once normalised by the cluster mass, Σ(SFR)/M_cl. The results are shown in Fig. 12.

Fig. 12. Star formation rate normalised by cluster mass as a function of redshift. Black squares represent median values from Dianoga simulations (grey points). See the text for a complete explanation of sample selection. The dashed black line is the best fit to simulations. Coloured points are observational data from Popesso et al. (2012), Ma et al. (2015), Smail et al. (2014), Santos et al. (2015), Wang et al. (2016), Miller et al. (2018), and Smith et al. (2019). The solid black line ∼(1 + z)7 shows an empirical fit to data suggested by Cowie et al. (2004) and Geach et al. (2006).

Σ(SFR)/M_cl is an increasing function of redshift with an observationally driven empirical parametrisation of (1 + z)ⁿ (e.g. Cowie et al. 2004; Geach et al. 2006). Popesso et al. (2012) used a sample of nine groups in the redshift range 0.1 <  z <  1.6 and nine clusters in the redshift range 0.1 <  z <  0.85 and derived the best fit to be (213 ± 44)×z^{1.33 ± 0.34} and (66 ± 23)  ×  z^{1.77 ± 0.36} for groups and clusters, respectively. Popesso et al. (2012) also showed that there is no evidence for a significant Σ(SFR)/M_cl − M_cl or M_cl − z correlation, concluding that Σ(SFR)/M_cl − z is a genuine correlation (not driven by a decreasing mass evolution with redshift within their sample). Recent observations of highly star forming protocluster regions suggest a stronger evolution with redshift, ∝(1 + z)⁷ (e.g. Smail et al. 2014; Ma et al. 2015; Santos et al. 2015; Smith et al. 2019), in line with the trend found by Cowie et al. (2004) for the number of star-forming ultraluminous infrared galaxies (ULIRGs) in the redshift range 0 <  z <  1.5. In Fig. 12 we add the previously cited protocluster regions at high redshift. We note that for these protoclusters the SFR is computed within an aperture of ∼1 pMpc, with two exceptions: Miller et al. (2018) computed the SFR within a 2D aperture of ∼130 pkpc and Wang et al. (2016) computed the SFR in a 2D aperture of ∼80 pkpc. All the SFRs derived assuming the Salpeter IMF are converted to a Chabrier IMF.

To build the comparison we considered clusters and groups in our simulations with a mass threshold varying with redshift. In particular, the minimum mass considered is M₂₀₀ >  10¹⁴ at z = 0, decreasing linearly with redshift up to M₂₀₀ >  10¹³ at z = 4.3. We made this choice to mimic as close as possible the minimum mass of the protoclusters and clusters in Fig. 12 at all redshifts. To mark a few examples, Popesso et al. (2012) clusters are in the mass range [3.9, 27.6]×10¹⁴ M_⊙; the cluster by Smith et al. (2019) at z ∼ 2 has an estimated mass of 0.5 × 10¹⁴ M_⊙; and the protocluster by Miller et al. (2018) at z ∼ 4.3 has an estimated mass of 1.16 ± 0.70 × 10¹³ M_⊙. We use M₂₀₀ (like Popesso et al. 2012) as an estimate of the cluster mass and we compute the instantaneous SFR considering all gas particles within R₂₀₀. This choice confers the advantage that it matches the aperture used by Wang et al. (2016) and Miller et al. (2018), the two observations at the highest redshifts (R₂₀₀ ∼ 300 pkpc and R₂₀₀ ∼ 150 pkpc at z = 3 and z = 4 respectively).

Simulated clusters show a clear evolution with redshift; however, the trend is shallower than in observations and is better described by ∝(1 + z)^{3.84 ± 0.15}. In particular, simulations predict a higher SFR at low redshift, reflecting the results already discussed in Fig. 4. At redshift z >  2, the predicted SFR is much lower, mirroring the outcome of the discussion of the previous sections. A similar mismatch with respect to observations has also been pointed out by Ragone-Figueroa et al. (2018) for the sSFR of BCGs of our lower resolution simulations.

7. Discussion

In recent years, a good number of observational studies have confirmed the detection of protocluster regions characterised by SFRs from several hundreds to several thousands of M_⊙ yr⁻¹ (Clements et al. 2014; Dannerbauer et al. 2014; Umehata et al. 2015; Wang et al. 2016; Oteo et al. 2018; Coogan et al. 2018; Miller et al. 2018; Gómez-Guijarro et al. 2019; Smith et al. 2019; Lacaille et al. 2019). These high values of SFR are often dominated by DSFGs, with typical SFRs from ∼100 M_⊙ yr⁻¹ to ∼1000 M_⊙ yr⁻¹. The physical reason of these high values of SFR is debated. Some observations suggest that starburst galaxies and SMGs are characterised by a high star formation efficiency (e.g. Daddi et al. 2010). Other recent observations suggest that the starbursting phase of galaxies is related to high gas fractions (e.g. Scoville et al. 2016; Gómez-Guijarro et al. 2019). Finally, some observations suggest that starburst galaxies are characterised by both high star formation efficiency and high gas fraction (Genzel et al. 2015; Béthermin et al. 2015).

In this section we focus on starburst galaxies in our simulations and the differences in terms of star formation efficiency and gas fraction with respect to galaxies in the observed protoclusters used as references in Sects. 4 and 5. We also investigate the implications for the subresolution model of star formation adopted in our simulations.

7.1. Starburst galaxies in numerical simulations

Galaxies are usually defined as starburst depending on how much their SFR is above the SFR of MS galaxies with the same mass and at the same redshift. Here, following Schreiber et al. (2015), we consider the threshold SFR/SFR_MS >  4. In the M_⋆−SFR plane, starburst galaxies do not only represent the tail of the MS distribution. Indeed several studies have shown that at fixed stellar mass and redshift, the distribution of galaxies around the MS is better described by a double Gaussian (Sargent et al. 2012; Schreiber et al. 2015), where the second component describes the population of starburst galaxies. This population is estimated to comprise 3% of star forming galaxies without significant redshift dependence (Schreiber et al. 2015).

Following the works mentioned above we study the starburst population in our simulations by plotting the distribution of galaxies around the MS in two mass bins at z = 2; see Fig. 13. Since in principle it is not guaranteed that the most star forming galaxies will be within protocluster regions, we also plot the results from the Magneticum simulations⁴. The Magneticum simulations are a set of hydrodynamical simulations of different cosmological volumes (Hirschmann et al. 2014; Ragagnin et al. 2017) performed with the same GADGET-3 code used in our simulations (see Hirschmann et al. 2014 for the differences in the AGN feedback implementation). From the Magneticum set, we consider the Box2 and Box2b (352 and 640 h⁻¹ Mpc respectively). The mass resolution is m_DM = 6.9 × 10⁸ h⁻¹ M_⊙, a factor of ten lower than the one used for this work.

Fig. 13. Star formation rate distribution of star forming galaxies at fixed stellar mass at z = 2. Blue points refer to Dianoga simulations, green squares to Magneticum Box2b, and red triangles to Magneticum Box2. SFRMS is computed independently for every simulation. NMS is the number of galaxies within the bin corresponding to SFR = SFRMS. Only bins with at least ten galaxies are shown. Coloured solid lines are Gaussian fits to simulations. Vertical black dashed line define the threshold above which data are used to estimate the fit.

In Fig. 13 the value of SFR_MS is computed for each simulation considering only active galaxies (see Sect. 4, Fig. 8). The two mass bins analysed are chosen following Sargent et al. (2012). We analysed only the two lower mass bins as at higher masses the number of galaxies in the Dianoga simulations is too low for a statistical analysis. Blue points are Dianoga simulations, red triangles are results for Box2, and green squares are results for Box2b. Only bins with at least ten galaxies are plotted. Solid lines are Gaussian fits to the data, which were obtained considering only galaxies with SFR >  0.5 × SFR_MS. This cut, marked as a vertical black dashed line in the plot, is needed to avoid considering galaxies on their way to being quenched (but still selected as active by our cut in sSFR) in the Gaussian fit. The three simulations are in very good agreement in both mass bins, despite the different box sizes and environment, and can all be fitted with a single Gaussian with a standard deviation estimated to be 0.19 <  σ <  0.21. This value is in agreement with the results of Sargent et al. (2012) (σ = 0.188), but is lower than the estimate of Schreiber et al. (2015) (σ = 0.31). Finally, the fraction of starburst galaxies (i.e. SFR/SFR_MS >  4) is 0.03%< f_SB <  0.2%, at least one order of magnitude lower than what was estimated by Schreiber et al. (2015).

As a final warning, it is also important to keep in mind that the observational estimates of the values of SFR, in particular for starburst galaxies, are affected by a number of uncertainties. A relevant role is played by the assumption of the IMF, as studies on the chemical abundances and abundance ratios suggest that starburst galaxies are characterised by a top-heavy stellar IMF (e.g. Romano et al. 2017, 2019).

If this proves to be correct, it could have important implications for the estimated values of SFR, affecting as a consequence also the conclusions reached in this section. Indeed, a more top-heavy IMF would lower the SFR values obtained through observations, while affecting numerical predictions to a much lesser extent (see also the discussion by Granato et al. 2015).

7.2. Star formation efficiency and gas fraction

In this section we study the cold gas and SFR properties of our simulated galaxies to determine which variable is more related to the underestimated normalisation of the MS at z ∼ 2 (see Sect. 4) and to the absence of starburst galaxies.

In Figs. 14 and 15 we study the gas fraction and the depletion time in simulations and observations at z ∼ 2 and z ∼ 4 respectively. The two are defined as

$$\begin{array}{c}\hfill f_{\mathrm{gas}}=\frac{M_{\mathrm{gas}}}{M_{\mathrm{gas}}+M_{\star }},\end{array}$$(10)

Fig. 14. Galaxy correlations at z = 2.3. Top panel: gas fraction as a function of stellar mass. Bottom panel: depletion time as a function of stellar mass. Grey circles refer to Dianoga simulations at z = 2.3. Green hexagons and blue squares are data from Gómez-Guijarro et al. (2019) and Wang et al. (2018) respectively. The orange dashed line is the functional form of Liu et al. (2019) for MS galaxies at z = 2.3, while the shaded region encompasses galaxies with an SFR that is four times lower and higher than MS galaxies. Red dashed lines are obtained combining the MS by Whitaker et al. (2014) and the integrated Kennicutt-Schmidt law from Sargent et al. (2014, orange line in the plot).

Fig. 15. Galaxy correlations at z = 4.3. Top panel: gas fraction as a function of stellar mass. Bottom panel: depletion time as a function of stellar mass. Grey circles refer to Dianoga simulations at z = 4.3. Brown circles are data from Hill et al. (2020). The orange dashed line is the functional form of Liu et al. (2019) for MS galaxies at z = 4.3, while the shaded region encompass galaxies with an SFR that is four times lower and higher than MS galaxies.

where M_gas is the cold gas mass that in simulations is computed considering only the cold phase of SPH particles, and

$$\begin{array}{c}\hfill t_{\mathrm{dep}}=\frac{M_{\mathrm{gas}}}{\mathrm{SFR}}·\end{array}$$(11)

We plot at z ∼ 2 and at z ∼ 4 all the galaxies in the observed protoclusters used in Sects. 4 and 5, respectively. Moreover, we add the functional forms of f_gas and t_dep from Liu et al. (2019), which depend on redshift, stellar mass, and relative distance from the MS (SFR/SFR_MS). The red dashed line in Fig. 14 was computed as follows: given a galaxy stellar mass, we assume that the expected SFR for a MS galaxy is given by the MS relation of Whitaker et al. (2014). Given the SFR, we assume that the mass of molecular gas is given by Eq. (4) of Sargent et al. (2014) for normal galaxies. Combining M_⋆, M_gas, and SFR we obtain the expected values of f_gas and t_dep for normal star forming galaxies. This procedure is supported also by recent observations with ALMA, which show that massive (M_⋆ >  10¹⁰ M_⊙) MS galaxies obey the star formation law for star-forming galaxies (Liu et al. 2019). From Figs. 14 and 15 we see that observational data are not in agreement among each other. In particular, data from Hill et al. (2020) at z ∼ 4 suggest that starburst galaxies are characterised by a shorter depletion time (and thus a high star formation efficiency) than normal star forming galaxies. This is supported by other observations (i.e. Daddi et al. 2010, Sargent et al. 2014, Liu et al. 2019). On the other hand, data from Gómez-Guijarro et al. (2019) and Wang et al. (2018) at z ∼ 2.4 are characterised by a high gas fraction but a star formation efficiency (or depletion time) that is consistent with the integrated Kennicutt-Schmidt law for normal star forming galaxies. This is also consistent with numerical simulations, where the position on the MS depends on the gas fraction (see Fig. 8 of Davé et al. 2019). The difference among the observational results is largely due to the different values of the parameter α_CO used to derive the molecular gas mass from the CO line luminosity $L_{\text{CO}(\text{1-0})}^{\prime }$ through Eq. (9). Gómez-Guijarro et al. (2019) used α_CO = 3.5, typical for normal star forming galaxies at solar metallicity. Wang et al. (2018) adopted the mass–metallicity relation by Genzel et al. (2015) to retrieve the galaxy metallicity for the members of their cluster, and then computed the metalicity-dependent value of α_CO following Genzel et al. (2015) and Tacconi et al. (2018). The values reported for the α_CO are in the range [4.06, 4.12], again consistent with normal star forming galaxies at solar metallicity. Hill et al. (2020), on the other hand, used α_CO = 1, typical for high-redshift SMGs. Due to these arguments, it remains unclear whether starburst galaxies in protocluster environments are mainly driven by a high gas fraction or a high star formation efficiency, as the results strongly depend on the assumptions needed to derive gas-related quantities.

If we look at the population of normal star forming galaxies at z ∼ 2 (red and orange dashed lines in Fig. 14) we see that our simulations match the observed star formation efficiency. On the other hand, the simulated gas fractions are consistently lower than observations. In particular, at M_⋆ = 10¹⁰ M_⊙, observations show a M_gas that is higher by a factor of ∼3.5. This factor is very similar to the difference in the observed and simulated MS (see Fig. 8). Thus, the lower normalisation in the simulated MS seems to be driven by an underestimated gas fraction.

It is also interesting to study the correlation between the MS galaxies and gas-related properties in simulations. In Fig. 16 we show the MS colour-coded with respect to f_gas (upper panel) and t_dep (lower panel). At fixed stellar mass, galaxies below the MS are characterised by both a long depletion time and low gas fraction. Vice versa, galaxies above the MS have both short depletion times and high gas fraction. This visual impression is also confirmed by the computation of the Pearson’s correlation coefficient for the two relations: Log(SFR/SFR_MS)−Log(f_gas) and Log(SFR/SFR_MS)−Log(t_dep). The results are r = 0.62 and r = −0.63 for the two relations respectively. Moreover, results of the two linear regressions suggest that in our simulations the position on the MS scales with f_gas and t_dep with similar slopes: SFR/SFR${}_{\mathrm{MS}}\propto f_{\mathrm{gas}}^{0.85}$ and SFR/SFR${}_{\mathrm{MS}}\propto t_{\mathrm{dep}}^{-1}$.

Fig. 16. Top panel: 2D histogram of MS star forming galaxies in simulations at z = 2.3. Each bin is colour-coded with the respective median value of fgas Bottom panel: same as upper panel, colour-coded with respect to tdep.

7.3. Simulation tests

In the previous sections we showed that our simulations underestimate the normalisation of the MS relation at z ∼ 2 by a factor of approximately three. Moreover, we have seen that the star formation model (Springel & Hernquist 2003) implemented in our code, with the current choice for the model parameters set to reproduce quiescent mode of star formation, does not reproduce the observed population of starburst galaxies at z >  2. While the normalisation of the MS seems to be mainly related to the gas fraction, it remains unclear whether we miss starburst galaxies because we do not correctly sample the star formation efficiency, the gas fraction, or both. Therefore, we performed a set of simulations aiming at checking whether the fraction of starburst galaxies and MS normalisation are sensitive to the choice of the parameters of the subgrid model. All the following simulations are performed for only one of our regions, a cluster with M₂₀₀ = 5.4 × 10¹⁴ M_⊙ at z = 0. In Figs. 17–19 we show the results for the MS, SFE, and gas fraction at z ∼ 3, respectively. The choice of the redshift is somewhat arbitrary, as we do not aim to compare simulations with particular observational data but to study the effect of different parameters on our results. Each panel refers to a different simulation, while we plot with grey circles the results for the reference simulation used in the previous sections. In the following we briefly discuss the specific changes for each simulation and the effects on the results.

Fig. 17. Main sequence of star forming galaxies at z = 3 for different simulations. Grey points refer to the results relative to the same region used for the tests with the set up used for this work. Different panels refer to: t0 0.3x: shorter timescale for star formation; SFTh 10x: increased density threshold for star formation; SFTht0: increased density threshold for star formation and shorter star formation timescale; A0 0.1x: reduced supernova thermal feedback; Tthr: AGN feedback implementation as in Ragone-Figueroa et al. (2018); No-AGN: no AGN feedback.

Fig. 18. Correlation between cold gas mass and SFR at z = 3 for different simulations. Grey points refer to the results relative to the same region used for the tests with the set up used for this work. Different panels refer to: t0 0.3x: shorter timescale for star formation; SFTh 10x: increased density threshold for star formation; SFTht0: increased density threshold for star formation and shorter star formation timescale; A0 0.1x: reduced supernova thermal feedback; Tthr: AGN feedback implementation as in Ragone-Figueroa et al. (2018); No-AGN: no AGN feedback.

Fig. 19. Correlation between stellar mass and cold gas mass at z = 3 for different simulations. Grey points refer to the results relative to the same region used for the tests with the set up used for this work. Different panels refer to: t0 0.3x: shorter timescale for star formation; SFTh 10x: increased density threshold for star formation; SFTht0: increased density threshold for star formation and shorter star formation timescale; A0 0.1x: reduced supernova thermal feedback; Tthr: AGN feedback implementation as in Ragone-Figueroa et al. (2018); No-AGN: no AGN feedback.

7.3.1. Increasing the star formation efficiency (t₀ 0.3x)

We recall that in the Springel & Hernquist (2003) model, the characteristic time for star formation, t_⋆, is $t_{\star }\propto t_0^{\star }{t}_{\text{dyn}}$, where $t_0^{\star }$ is a parameter usually tuned to reproduce the Kennicutt relation (see Sect. 2.2). Here we increase the efficiency to match the observed SFE of Hill et al. (2020) (see Fig. 15) by lowering $t_0^{\star }$ by a factor of three. The results of this test are shown in the top-left panel of the figures. From Fig. 18 we see that indeed the SFE is higher, but there is little difference in the MS (see Fig. 17). Moreover, there is no difference in the fraction of starburst galaxies. In fact, the model is so tightly self-regulated that in response to a high SFE we have a lower gas fraction (see Fig. 19), resulting in similar SFRs.

7.3.2. Increasing the star formation threshold (SFTh 10x and SFTh t₀)

In the SFTh 10x simulation we increased the density threshold, ρ_thr, by a factor of ten; this threshold is used to decide whether or not a gas particle becomes multiphase (we reiterate the fact that only multiphase particles can form stars; see Sect. 2.2 and Springel & Hernquist 2003). Increasing this threshold should allow a larger reservoir of gas to accumulate and reach higher densities before starting to produce stars, increasing the gas fraction and the overall SFR. However, from the top-right panel of Figs. 17–19 we see that we do not have major differences in terms of MS normalisation, SFE, and gas fraction. The only appreciable difference is the reduction of the most massive galaxies and the increase of passive galaxies. Indeed, higher densities at the centre of galaxies also means more gas accretion onto the central BH and consequently a stronger AGN feedback.

In the SFTh t₀ run (central-left panel) we increased both the density threshold for multiphase particles by a factor of ten and the SFE by a factor of three. Again, the self-regulation of the star formation model and the AGN feedback prevent any appreciable difference with respect to our fiducial run.

7.3.3. Increasing timescale for cold gas evaporation (A₀ 0.1x)

Following Springel & Hernquist (2003), even if the subgrid model is explicitly constructed to reproduce quiescent star formation, starburst should arise whenever the timescale for star formation is shorter than the timescale for the evaporation of cold gas. In fact, in this regime self-regulation is expected to break down with cold gas transformed into stars before it can be evaporated by stellar feedback. In practice, the relation that should be satisfied is:

$$\begin{array}{c}\hfill \frac{t_{\mathrm{c}}}{t_{\star }}={\left(\frac{\rho }{{\rho }_{\mathrm{thr}}}\right)}^{4/5}\frac{1}{\beta A_0}>1,\end{array}$$(12)

where ρ_thr is the density threshold for a particle to become multiphase, β is the fraction of stars that instantly die as supernovae, and A₀ is a parameter of the model that defines the energy of supernovae used to evaporate cold gas. In this test we reduced the value of A₀ by a factor of ten. From the results shown in the central-right panels we see no improvement in terms of starburst galaxies. Thus, even if we checked that single gas particles satisfy Eq. (12), this is not sufficient to have a high-enough integrated value of SFR.

7.3.4. Varying AGN feedback implementation (Tthr)

To quantify the effect of a specific aspect of the AGN feedback implementation on our results we also ran a simulation with the same AGN feedback prescription as that of Ragone-Figueroa et al. (2018). We recall that in that setup there is an extra condition on the temperature (T <  T_thr) to consider a particle as multiphase and that the energy released by AGN feedback is used to evaporate molecular clouds, while in the current implementation the energy is coupled only to the hot phase of multiphase particles. From the bottom-left panels of Figs. 17–19 we can see that the only difference with respect to our fiducial run is that in this case we have less massive galaxies. This is expected from the results shown in Sects. 3.1 and 3.2.1, where it is clear that the feedback implementation of Ragone-Figueroa et al. (2018) is more effective in quenching star formation.

7.3.5. No AGN feedback (No-AGN)

Finally, we also performed a simulation without AGN feedback (bottom-right panels). This is of course to test an extreme scenario, as the absence of AGN feedback would result in GSMF, BCG masses, and SFR inconsistent with low-redshift observations. From Fig. 17, we see that in this run we have fewer galaxies on their way to becoming passive, and more massive galaxies, as expected. However, the MS retain the same normalisation and there is no sign of an increased fraction of starburst galaxies. Moreover, it is interesting to note that in the No-AGN run the SFE is higher in the low-mass regime (see Fig. 18). This difference is due to the fact that without AGN feedback the gas reaches higher density, especially in the low-mass regime where the feedback is more efficient in expelling gas outside the shallow potential wells of galaxies.

8. Conclusions

Here, we study the SFR of simulated protocluster regions and the gas properties of protocluster galaxies in the redshift range 2 <  z <  4, and we compared them with observations. Our work is based on a subsample of the Dianoga simulations (Bonafede et al. 2011). In particular, we used 12 clusters, 7 of which very massive (M₂₀₀ >  8 × 10¹⁴ h⁻¹ M_⊙). The simulations are carried out with GADGET3, a modified version of the public code GADGET2, which implements a SPH scheme for hydrodynamics and treats the unresolved baryonic physics through various subgrid models. In particular, we use the Springel & Hernquist (2003) model for star formation and a thermal AGN feedback. Sections 3.1 and 3.2.1 present the degrees of freedom of the AGN feedback implementation. With the implementation of Ragone-Figueroa et al. (2018), where a temperature threshold is used to define multiphase gas particles and the energy released by AGN feedback is used to evaporate their cold phase, we match the observed correlation between cluster and BCG mass, but the normalisation of the galaxy stellar mass function is lower by a factor of approximately two with respect to observations (see Bassini et al. 2019, Appendix B). On the other hand, without the temperature threshold and coupling the energy released by AGN feedback only with the hot phase of gas particles, we match the GSMF but we get overly massive BCGs (a factor of ∼2 more massive; see Sect. 3.2.1). In this paper, we use the latter implementation, which maximises the value of the SFR in the high-redshift regime (z ∼ 2 − 4) in which we are interested. Our main results can be summarised as follows:

• At z ∼ 2, simulations under-predict the SFR of highly star forming protocluster regions by a factor of four or even larger, in line with the results we presented in Granato et al. (2015) based on a larger set of lower resolution simulations. This result is indeed stable against numerical resolution and is the combination of two effects: (i) simulations under-predict the normalisation of the MS at 2 <  z <  2.5 by a factor of three; (ii) simulations predict a fraction of starburst galaxies, defined as galaxies with an SFR at least four times higher than MS galaxies, of [0.2%−0.03%], which is at least a factor of ten lower than what was found by recent observations (Schreiber et al. 2015). We verified that this result is independent of the environment by performing the same analysis on the Magneticum cosmological boxes of 352 and 640 h⁻¹ Mpc per side (Hirschmann et al. 2014; Ragagnin et al. 2017).

• At z ∼ 4, simulations correctly reproduce the MS normalisation, but fail to reproduce the starburst population. Indeed, simulations under-predict the SFR of highly star forming protocluster regions by a factor of four.

• In our simulations, the normalisation of the MS strongly depends on the gas fraction. Comparison with observations suggests that simulations under-predict the gas fraction in galaxies at the peak of the cosmic star formation rate density and consequently the normalisation of the MS.

• In numerical simulations, the position on the MS depends on both the gas fraction and the star formation efficiency. However, observations of galaxy properties in dense environments are affected by uncertainties on the assumptions needed to derive gas-related quantities. Therefore, it remains unclear whether simulations under-reproduce starburst galaxy populations because of a low gas fraction or a low star formation efficiency.

• Our results indicate that the adopted model of star formation (i.e. Springel & Hernquist 2003) efficiently reproduces the self-regulated evolution of quiescent low-redshift star formation but is not suitable for capturing violent events like high-redshift starbursts. We verified that our results are robust and the conclusions hold for a wide range of values of the model parameters and do not depend on the implementation of the AGN feedback.

Finally, we note that even though simulations tend to under-reproduce the level of SFR at high redshift, the stellar mass at z = 0 is even higher than what is suggested by observations (see Fig. 3). Therefore, as already pointed out by Granato et al. (2015), the star formation history of protoclusters must be characterised by peaks that are higher and shorter than those of numerical simulations. Given the results obtained in this work, it seems unfeasible to achieve this goal without any substantial modification of the model of star formation, as imposing a self-regulated regime of star formation does not allow us to reach high-enough values of SFR. The large amount of data that is becoming available at high redshift from instruments like ALMA will help to put constraints on high-redshift galaxy properties and to accordingly improve the degree of realism of star formation models implemented in cosmological simulations of galaxy formation.

Acknowledgments

We thank the anonymous referee for the careful and constructive reading of the paper and for his/her useful suggestions. We thank L. Boco, L. Pantoni, and M. Valentini for helpful discussions. We would like to thank Volker Springel for making the GADGET-3 code available to us. We thank Romeel Davé for sharing data from Simba simulations; Quan Guo for sharing data from EAGLE simulation and GALFORM and L-GALAXIES semi-analytical models; Gabriella De Lucia and Fabio Fontanot for sharing data from GAEA semi-analytical model. VB acknowledges support by the DFG project nr. 415510302. This project has received funding from: ExaNeSt and Euro Exa projects, funded by the European Union Horizon 2020 research and innovation program under grant agreement No 671553 and No 754337, the agreement ASI-INAF n.2017-14-H.0; the Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina (CONICET); the Secretaría de Ciencia y Técnica de la Universidad Nacional de Córdoba – Argentina (SeCyT); the European Union Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement No 734374, PRIN-MIUR 2015W7KAWC, the INFN INDARK grant. NRN acknowledges financial support from the “One hundred top talent program of Sun Yat-sen University” grant N. 71000-18841229. Simulations have been carried out using MENDIETA Cluster from CCAD-UNC, which is part of SNCAD-MinCyT (Argentina); MARCONI at CINECA (Italy), with CPU time assigned through grants ISCRA B, and through INAF-CINECA and University of Trieste – CINECA agreements; at the Tianhe-2 platform of the Guangzhou Supercomputer Center by the support from the National Key Program for Science and Technology Research and Development (2017YFB0203300). The post-processing has been performed using the PICO HPC cluster at CINECA through our expression of interest.

References

All Figures

Fig. 1. Correlation between the galaxies stellar mass and the central SMBHs mass. Observational data are taken from McConnell & Ma (2013) (dashed black line) and from Gaspari et al. (2019) (red circles). The simulated stellar masses for satellite galaxies (cyan points) are obtained considering the star particles bound to the substructure (accordingly to Subfind) and within 50 pkpc from its centre. The mass of the central galaxies (dark-blue squares) is obtained by summing over all stellar particles within an aperture of 0.15 × R500.

In the text

Fig. 2. GSMF at z = 0. Observational data are taken from Bernardi et al. (2013) (black solid line). Simulation data are derived considering as stellar mass the sum of all stellar particles bound to the galaxy by Subfind (red triangles), and the same sum restricted to particles within 50 pkpc (green hexagon) and 30 pkpc (blue squares). Error bars are computed assuming Poissonian errors. The simulated GSMF is normalised following Eq. (8). Filled and empty marks represent the mass bins with respectively more than and less than ten galaxies.

In the text

Fig. 3. Correlation between BCG stellar mass and M500 at z = 0. Observations are taken from DeMaio et al. (2018) (blacks quares) and Kravtsov et al. (2018) (black triangles). The simulated values are shown as blue points. The red hexagon refers to the BCG that lost its central BH (see Sect. 2.3.1). The orange line is the fit to LR simulations (Ragone-Figueroa et al. 2018). BCG masses are obtained summing over all stellar particles bound to the main subhalo of a group or cluster by Subfind (BCG+ICL) and within a 2D aperture of 50 pkpc.

In the text

Fig. 4. Star formation rate of BCGs in observations and simulations. Grey circles are BCGs of our simulations from different snapshots, while the grey triangle represents the BCG that lost its central BH at z ∼ 4 (see Sect. 2.3.1). BCGs from the same snapshot are shifted only for visualisation purposes. The median values are shown as blue circles and the vertical bars indicate the range between the 16th and 84th percentiles. A 2D aperture of 30 pkpc is used. Red squares are BCGs from the sample of McDonald et al. (2018) (see text for more details).

In the text

Fig. 5. Specific SFR of BCGs in observations and simulations. Grey circles are BCGs of our simulations from different snapshots (blue circles are median values with 16th and 84th percentiles), while the grey triangle is used for the BCG that lost its central BH at z ∼ 4 (see Sect. 2.3.1). BCGs from the same snapshot are shifted only for visualisation purposes. A 2D aperture of 30 pkpc is used. Red squares are BCGs from the sample of McDonald et al. (2018) (see text for more details).

In the text

Fig. 6. Star formation rate of protocluster regions at z ∼ 2 in observations and simulations within an aperture of ∼2 pMpc. Red bands refer to two clumps from Clements et al. (2014), black solid lines refer to four fields from Stevens et al. (2010) and analysed by Clements et al. (2014). The blue square highlights to the Spiderweb structure (Dannerbauer et al. 2014). The green square and green band show the two protoclusters analysed by Kato et al. (2016) (HS1700 and 2QZCluster, respectively). Black circles and triangles refer to numerical simulations, where the SFR is plotted against protocluster mass (see text). We used black circles for groups which end up in the central cluster of the region at z = 0, and black triangles otherwise.

In the text

Fig. 7. Star formation rate of protocluster regions at 2 <  z <  2.6 in observations and simulations within an aperture of ∼100 pkpc. Green bands refer to two protoclusters from Gómez-Guijarro et al. (2019), the blue square refers to Wang et al. (2016), and the red square refers to Coogan et al. (2018). Black circles and triangles refer to numerical simulations, where the SFR is plotted against protocluster core mass. We used black circles for groups which end up in the central cluster of the region at z = 0, and black triangles otherwise.

In the text

Fig. 8. Star formation rate as a function of galaxy stellar mass at z ∼ 2.3. Red solid and dashed lines are observational data from Whitaker et al. (2014) and Schreiber et al. (2015), respectively. Green hexagons and blue squares are galaxies from the protoclusters of Gómez-Guijarro et al. (2019) and the cluster of Wang et al. (2018), respectively. Grey points are galaxies in our simulations. Black dashed line fix the distinction between active and passive galaxies (Pacifici et al. 2016). Black points represent median values of star forming galaxies with 16th and 84th percentiles. Both SFRs and stellar masses are computed considering a 3D aperture of 30 pkpc.

In the text

Fig. 9. Main sequence of star forming galaxies at z ∼ 2. Red triangles are observational data from Whitaker et al. (2014). The black line shows median values for our simulations. Coloured solid and dashed lines are data from other cosmological simulations and semi-analytical models respectively. In particular: Eagle (orange solid line, Guo et al. 2016), TNG300 (red solid line, Donnari et al. 2019), Simba (yellow solid line, Davé et al. 2019), Galform (green dashed line, Guo et al. 2016), L-galaxies (dark green dashed line, Guo et al. 2016), and GAEA (blue dashed line, Hirschmann et al. 2016). For the GAEA model we also show the results obtained considering only galaxies that at z = 0 are within galaxy clusters with mass > 1014.25 M⊙ (see text for more details).

In the text

Fig. 10. Star formation rate as a function of M500 at z ∼ 4.3. The blue square and green line are the observed values of Miller et al. (2018) and Oteo et al. (2018) protoclusters respectively. Black symbols refer to the SFR computed in a cylinder 1 pMpc long and within a circular aperture of 130 pkpc in our simulations. The five most massive groups of each region are shown. We used black circles for groups which end up in the central cluster of the region at z = 0, and black triangles otherwise.

In the text

Fig. 11. Star formation rate as a function of galaxy stellar mass at z ∼ 4.3. The red line shows observational data from Steinhardt et al. (2014). Orange dots represent galaxies of SPT2349-56 as analysed in Hill et al. (2020). Grey points are galaxies in our simulations. The black dashed line fixes the distinction between quiescent and star forming galaxies (Pacifici et al. 2016). Black points represent median values with 16th and 84th percentiles for star forming galaxies. Green circles are galaxies from the simulated protoclusters shown in Fig. 10.

In the text

Fig. 12. Star formation rate normalised by cluster mass as a function of redshift. Black squares represent median values from Dianoga simulations (grey points). See the text for a complete explanation of sample selection. The dashed black line is the best fit to simulations. Coloured points are observational data from Popesso et al. (2012), Ma et al. (2015), Smail et al. (2014), Santos et al. (2015), Wang et al. (2016), Miller et al. (2018), and Smith et al. (2019). The solid black line ∼(1 + z)7 shows an empirical fit to data suggested by Cowie et al. (2004) and Geach et al. (2006).

In the text

Fig. 13. Star formation rate distribution of star forming galaxies at fixed stellar mass at z = 2. Blue points refer to Dianoga simulations, green squares to Magneticum Box2b, and red triangles to Magneticum Box2. SFRMS is computed independently for every simulation. NMS is the number of galaxies within the bin corresponding to SFR = SFRMS. Only bins with at least ten galaxies are shown. Coloured solid lines are Gaussian fits to simulations. Vertical black dashed line define the threshold above which data are used to estimate the fit.

In the text

Fig. 14. Galaxy correlations at z = 2.3. Top panel: gas fraction as a function of stellar mass. Bottom panel: depletion time as a function of stellar mass. Grey circles refer to Dianoga simulations at z = 2.3. Green hexagons and blue squares are data from Gómez-Guijarro et al. (2019) and Wang et al. (2018) respectively. The orange dashed line is the functional form of Liu et al. (2019) for MS galaxies at z = 2.3, while the shaded region encompasses galaxies with an SFR that is four times lower and higher than MS galaxies. Red dashed lines are obtained combining the MS by Whitaker et al. (2014) and the integrated Kennicutt-Schmidt law from Sargent et al. (2014, orange line in the plot).

In the text

Fig. 15. Galaxy correlations at z = 4.3. Top panel: gas fraction as a function of stellar mass. Bottom panel: depletion time as a function of stellar mass. Grey circles refer to Dianoga simulations at z = 4.3. Brown circles are data from Hill et al. (2020). The orange dashed line is the functional form of Liu et al. (2019) for MS galaxies at z = 4.3, while the shaded region encompass galaxies with an SFR that is four times lower and higher than MS galaxies.

In the text

	Fig. 16. Top panel: 2D histogram of MS star forming galaxies in simulations at z = 2.3. Each bin is colour-coded with the respective median value of fgas Bottom panel: same as upper panel, colour-coded with respect to tdep.
In the text

Fig. 17. Main sequence of star forming galaxies at z = 3 for different simulations. Grey points refer to the results relative to the same region used for the tests with the set up used for this work. Different panels refer to: t0 0.3x: shorter timescale for star formation; SFTh 10x: increased density threshold for star formation; SFTht0: increased density threshold for star formation and shorter star formation timescale; A0 0.1x: reduced supernova thermal feedback; Tthr: AGN feedback implementation as in Ragone-Figueroa et al. (2018); No-AGN: no AGN feedback.

In the text

Fig. 18. Correlation between cold gas mass and SFR at z = 3 for different simulations. Grey points refer to the results relative to the same region used for the tests with the set up used for this work. Different panels refer to: t0 0.3x: shorter timescale for star formation; SFTh 10x: increased density threshold for star formation; SFTht0: increased density threshold for star formation and shorter star formation timescale; A0 0.1x: reduced supernova thermal feedback; Tthr: AGN feedback implementation as in Ragone-Figueroa et al. (2018); No-AGN: no AGN feedback.

In the text

Fig. 19. Correlation between stellar mass and cold gas mass at z = 3 for different simulations. Grey points refer to the results relative to the same region used for the tests with the set up used for this work. Different panels refer to: t0 0.3x: shorter timescale for star formation; SFTh 10x: increased density threshold for star formation; SFTht0: increased density threshold for star formation and shorter star formation timescale; A0 0.1x: reduced supernova thermal feedback; Tthr: AGN feedback implementation as in Ragone-Figueroa et al. (2018); No-AGN: no AGN feedback.

In the text