Issue 
A&A
Volume 630, October 2019



Article Number  A144  
Number of page(s)  16  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201935383  
Published online  08 October 2019 
Black hole mass of central galaxies and cluster mass correlation in cosmological hydrodynamical simulations
^{1}
Dipartimento di Fisica dell’Università di Trieste, Sezione di Astronomia, Via Tiepolo 11, 34131 Trieste, Italy
email: gigibassini@gmail.com
^{2}
INAF, Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34131 Trieste, Italy
^{3}
IFPU – Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy
^{4}
INFN, Instituto Nazionale di Fisica Nucleare, Trieste, Italy
^{5}
Universidad Nacional de Córdoba, Observatorio Astronómico de Córdoba, Córdoba, Argentina
^{6}
CONICET, Instituto de Astronomía Teórica y Experimental, Córdoba, Argentina
^{7}
HarvardSmithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA
^{8}
University Observatory Munich, Scheinerstraße 1, 81679 Munich, Germany
^{9}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStraße 1, 85748 Garching bei München, Germany
^{10}
Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 085441001, USA
Received:
1
March
2019
Accepted:
17
August
2019
Context. The correlations between the properties of the brightest cluster galaxy (BCG) and the mass of its central supermassive black hole (SMBH) have been extensively studied from a theoretical and observational angle. More recently, relations connecting the SMBH mass and global properties of the hosting cluster, such as temperature and mass, were observed.
Aims. We investigate the correlation between SMBH mass and cluster mass and temperature, their establishment and evolution. We compare their scatter to that of the classical M_{BH} − M_{BCG} relation. Moreover, we study how gas accretion and BHBH mergers contribute to SMBH growth across cosmic time.
Methods. We employed 135 groups and clusters with a mass range 1.4 × 10^{13} M_{⊙} − 2.5 × 10^{15} M_{⊙} extracted from a set of 29 zoomin cosmological hydrodynamical simulations where the baryonic physics is treated with various subgrid models, including feedback by active galactic nuclei.
Results. In our simulations we find that M_{BH} correlates well with M_{500} and T_{500}, with the scatter around these relations compatible within 2σ with the scatter around M_{BH} − M_{BCG} at z = 0. The M_{BH} − M_{500} relation evolves with time, becoming shallower at lower redshift as a direct consequence of hierarchical structure formation. On average, in our simulations the contribution of gas accretion to the total SMBH mass dominates for the majority of the cosmic time (z > 0.4), while in the last 2 Gyr the BHBH mergers become a larger contributor. During this last process, substructures hosting SMBHs are disrupted in the merger process with the BCG and the unbound stars enrich the diffuse stellar component rather than increase BCG mass.
Conclusions. From the results obtained in our simulations with simple subgrid models we conclude that the scatter around the M_{BH} − T_{500} relation is comparable to the scatter around the M_{BH} − M_{BCG} relation and that, given the observational difficulties related to the estimation of the BCG mass, clusters temperature and mass can be a useful proxy for the SMBHs mass, especially at high redshift.
Key words: galaxies: clusters: intracluster medium / Galaxy: nucleus / methods: numerical / Xrays: galaxies: clusters
© ESO 2019
1. Introduction
It is well known that galaxies of every morphology host super massive black holes (SMBHs) at their center. Interestingly the mass of these SMBHs correlate well with a number of bulge properties of the host galaxy such as bulge stellar mass (e.g., Magorrian et al. 1998; Marconi & Hunt 2003; Häring & Rix 2004; Hu 2009; Sani et al. 2011; Kormendy & Ho 2013 for a review), bulge luminosity (e.g., Kormendy & Richstone 1995; McLure & Dunlop 2002; Hu 2009; Sani et al. 2011), and especially the bulge stellar velocity dispersion (e.g., Ferrarese & Merritt 2000; Gebhardt et al. 2000; Merritt & Ferrarese 2001; Tremaine et al. 2002; Wyithe 2006; Hu 2008; Gültekin et al. 2009; Beifiori et al. 2012; McConnell & Ma 2013). These correlations are as tight as to be often used to estimate the SMBH mass when dynamical measurements are not available and suggest a coevolution between SMBH and the hosting galaxy, although the main physical processes involved are still debated. In the last 20 years many possibilities have been proposed. The most commonly suggested and widely accepted mechanism is the activegalacticnuclei (AGN) feedback. In this scenario gas settles around the SMBH radiating energy at a rate of ∼ηṀc^{2}, with η ∼ 0.1. If the feedback is strong enough to overcome the binding energy, cold gas is expelled from the galaxy halting both star formation and accretion around the central SMBH (e.g., Fabian 1999; Granato et al. 2004; Di Matteo et al. 2005; Hopkins et al. 2006). However, other authors have shown that the M_{BH}–M_{⋆} relation can arise or be contributed by noncausal processes. In this scenario the observed M_{BH}–M_{⋆} relation follows from hierarchical galaxy mergers starting from an uncorrelated distribution of M_{BH} and M_{⋆} (e.g., Peng 2007; Jahnke & Macciò 2011).
For massive elliptical galaxies, it has also been suggested that the main correlation is not with the bulge properties but with the host dark matter halo (e.g., Ferrarese 2002; Bandara et al. 2009; Booth & Schaye 2010; Volonteri et al. 2011). In this respect, between all galaxies a special position is occupied by the brightest cluster galaxies (BCGs), which sit at the center of galaxy clusters and host the most massive SMBHs in the universe. Since the luminosity or mass of these galaxies have been found to relate to the hosting halo both in observational (e.g., Lin & Mohr 2004; Brough et al. 2008) and numerical (e.g., RagoneFigueroa et al. 2018) studies, a correlation between M_{BH} and cluster mass is expected (e.g., Mittal et al. 2009), although the SMBH and cluster global properties ought not be necessarily connected. The argument is not completely settled as observations of another class of objects, bulgeless spiral galaxies, suggest that SMBHs do not correlate directly with dark matter, the latter being estimated from the circular rotation velocities of gas in the outskirts (e.g., Kormendy & Bender 2011; Sabra et al. 2015).
Churazov et al. (2005) proposed a toy model for the evolution of massive elliptical galaxies and their central SMBHs. In this simple model, the equivalence between gas cooling and AGN heating leads to a correlation between SMBH mass and the stellar velocity dispersion (assuming a rough , we can further convert it into a temperature correlation). Gaspari & Sdowski (2017) proposed that the feeding via chaotic cold accretion (CCA) boosts the AGN feedback in a tight selfregulated feedback loop that prevents catastrophic cooling flows for several gigayears, driving a direct physical connection between the SMBH mass and intra cluster medium (ICM) properties, such as plasma temperature and luminosity. With a sample of groups and clusters of galaxies, Bogdán et al. (2018) corroborated this prediction, further finding a scatter in the M_{BH} − T_{500} relation that is lower than that of the M_{BH} − M_{⋆} relation. Phipps et al. (2019) enlarged the sample of Bogdán et al. (2018), focusing on the cluster scales and retrieving the SMBH mass from the AGN fundamental plane (Merloni et al. 2003). Gaspari et al. (2019) based their sample on SMBHs with direct dynamical mass measurements, considering different type of galaxies (BCGs, BGGs, ETGs, S0s, and massive LTGs) and exploring correlations with thermodynamic properties in addition to the Xray luminosity and temperature (such as the gas pressure, gas mass, and the parameter Y_{X} = M_{gas} × T) over different extraction radii (galactic, core, R_{500}). Finally, Lakhchaura et al. (2019) considered a sample of 41 ETGs and confirmed the tight correlation between the SMBH mass and the Xray temperature measured within the effective radius.
All these recent works suggest that the growth of SMBHs in BCG is regulated by physical processes that also influence the thermodynamical properties of the ICM. Various phenomena dominate at different scales, from few kiloparsec to megaparsec. In this work we focus on the largest scale. Specifically, in this paper, we investigate the correlation between the mass of SMBHs in BCGs and the global properties of the hosting cluster, such as temperature and mass measured within R_{500} (defined as the radius of the sphere of mass M_{500} whose enclosed density is 500 times the critical density of the Universe at a specific redshift, 3H(z)^{2}/8πG) by employing a set of cosmological hydrodynamical simulations centered on massive clusters. These zoomin simulations include a number of subgrid models for radiative cooling, star formation and associated feedback, metal enrichment, and chemical evolution and they implement recipes for SMBH accretion and consequent AGN feedback (RagoneFigueroa et al. 2013). In past works with a similar set of simulations, we showed that the AGN feedback at the center of galaxy clusters leads to an appropriate description of the observed ICM thermodynamic quantities (such as entropy, gas density, temperature, and thermal pressure) and metallicity (Rasia et al. 2015; Planelles et al. 2017; Biffi et al. 2017). Given these previous results, we investigate further these simulated regions to study the physical processes that tie SMBH to the hosting cluster. In particular in this work we aim at answering the following questions: (1) Do numerical simulations reproduce the observed T_{500} − M_{BH} and M_{500}–M_{BH} relations? (2) Which are the processes that lead to the observed relations? (3) Do the relations evolve with redshift? (4) Through which channels (e.g., gas accretion or BHBH mergers) do SMBHs grow in time? (5) Is M_{500} as appropriate as M_{BCG} to probe M_{BH}?
The paper is structured as follows: in Sect. 2 we describe the numerical simulations employed. In Sect. 3 we detail how the quantities of interest are computed from simulations and the method employed for linear fitting. In Sect. 4 we present our results, that we discuss in Sect. 5 before concluding in Sect. 6.
2. Simulations
Our analysis is based on a set of 29 cosmological and hydrodynamical zoomin simulations centered on massive galaxy clusters evolved in a Λ cold dark matter (ΛCDM) model with parameters: Ω_{m} = 0.24, Ω_{b} = 0.04, n_{s} = 0.96, σ_{8} = 0.8, and H_{0} = 100 h km s^{−1} Mpc^{−1} = 72 km s^{−1} Mpc^{−1}. These regions were selected from a parent Nbody cosmological volume of 1 h^{−3} Gpc^{3} and resimulated at higher resolution with the inclusion of baryons (for a detailed description of the initial conditions see Bonafede et al. 2011).
The resimulated regions are centered around the 24 most massive clusters of the parental box with mass M_{200} ≥ 8 × 10^{14} h^{−1} M_{⊙} and 5 isolated groups with M_{200} within [1 − 4]×10^{14} h^{−1} M_{⊙}. In the highresolution regions the mass of DM particles is m_{DM} = 8.47 × 10^{8}h^{−1} M_{⊙} and the initial mass of the gas particle is m_{gas} = 1.53 × 10^{8}h^{−1} M_{⊙}. The Plummer equivalent gravitational softening for DM particles is set to ϵ = 5.6 h^{−1} kpc in comoving units at redshift higher than z = 2 and in physical units afterward. The gravitational softening lengths of gas, stars, and black hole particles are fixed in comoving coordinates to 5.6 h^{−1} kpc, 3 h^{−1} kpc, and 3 h^{−1} kpc, respectively.
The simulations were carried out with the code GADGET3, a modified version of the TreePM SmoothedParticleHydrodynamics (SPH) public code GADGET2 (Springel 2005). Our simulations are performed with an improved version that accounts for modifications of the hydrodynamic scheme to better capture hydrodynamical instabilities (see Beck et al. 2016). These changes include a higher order kernel function, a time dependent artificial viscosity scheme, and a time dependent artificial conduction scheme.
The set of zoomin simulations treats unresolved baryonic physics through various subgrid models. A detailed description can be found in Planelles et al. (2017) or Biffi et al. (2017); we briefly summarize here the main aspects. The prescription of metaldependent radiative cooling follows Wiersma et al. (2009). The model of star formation and associated feedback prescriptions are implemented according to the original model by Springel & Hernquist (2003), and metal enrichment and chemical evolution following the formulation by Tornatore et al. (2007). The yields used in our simulations are specified in Biffi et al. (2018). The AGN feedback model is implemented as described in Appendix A of RagoneFigueroa et al. (2013) with one important modification (see, RagoneFigueroa et al. 2018): the distinction between cold mode and hot mode gas accretion (see also Rasia et al. 2015). In practice, the gas accretion is the minimum between the Eddington limit and the αmodified Bondi accretion rate:
with α equal to 10 and 100 for hot (T > 5 × 10^{5} K) and cold (T < 5 × 10^{5} K) gas respectively. The 100× boost during the cold gas accretion mimics the impact of CCA (e.g., Gaspari & Sdowski 2017). M_{BH} is the SMBH mass. All other quantities relate to the gas and are smoothed over 200 gas particles with a kernel centered at the position of the black hole: ρ is the gas density, c_{s} is the sound speed of the gas surrounding the SMBH, and ν_{BH} is the relative velocity between the SMBH and bulk velocity of the gas. We note that for the SMBHs considered in this paper (e.g., SMBHs at the center of BCGs) the radius which enclose 200 gas particles is ∼50 kpc at z = 0. This radius decreases at higher redshift, having a value of ∼30 kpc at z = 2. In order to avoid wandering black holes, they are repositioned at each time step on the position of the most bound particle within the SMBH softening length. This calculation is restricted to particles with relative velocity with respect to the SMBH below 300 s^{−1} km. This condition avoids that the SMBH particle “jumps” into a close flyby structure that would displace it from the cluster center.
A fraction ϵ_{r} of the energy associated to the gas directly fueling the SMBH through accretion is radiated away and a fraction ϵ_{f} of this energy is thermally coupled to the surrounding gas particles. The value of ϵ_{r} is fixed to 0.07, while that of ϵ_{f} depends on the mode of the AGN: during the quasar mode, meaning for Ṁ_{BH}/Ṁ_{Edd} > 0.01, ϵ_{f} = 0.1, while during the radio mode ϵ_{f} is increased to 0.7 (see, RagoneFigueroa et al. 2018). The exact values of both parameters are chosen to reproduce the observed correlation between stellar mass and SMBH mass in galaxies (see Fig. 1). Recently some studies (e.g., Shankar et al. 2016, 2019) pointed out that the currently observed relation between stellar mass and SMBH mass may be biased high. We study the effect of artificially increasing feedback parameters on our results in the appendix.
Fig. 1.
Correlation between stellar mass and SMBH mass in observations and simulations. Small lightblue points represent nonBCG simulated galaxies, large black dots represent simulated BCGs. Yellow, green, red, brown, and orange crosses represent the observational data with their error bars taken from McConnell & Ma (2013), Main et al. (2017), Savorgnan et al. (2016), Bogdán et al. (2018), and Gaspari et al. (2019) respectively. The dashed black line is a linear bestfit of the sample of different type of galaxies by McConnell & Ma (2013). See text for details about M_{BH} and M_{⋆} definition and measurement. 
SMBHs of mass 4.4 × 10^{5} M_{⊙} are seeded at the position of the most bound particle of the structures that, identified by the FriendsofFriends algorithm, simultaneously satisfy all the following conditions: the stellar mass of the structure is greater than 2.2 × 10^{10} M_{⊙} and it is higher than 5% of the darkmatter halo mass; the ratio between the gas mass and the stellar mass is higher than 0.1; no other central SMBH is already present. The mass of the seeding is consistent with the expectation of the direct collapse. Under these conditions, the seeding of the SMBH happens in galaxies that have enough gas to promptly feed it. Two SMBH particles merge whenever their relative velocity v_{rev} is smaller than 0.5 × c_{s} and their distance r is less than twice the SMBH softening length. When a BHBH merger happens, the SMBH particle of the most massive SMBH gains the mass of the merged SMBH. The strategy to position the SMBH, the recipe to implement the AGN feedback, and slightly larger softening lengths are the only differences with respect to the simulation setup of the runs presented in Rasia et al. (2015), Biffi et al. (2017, 2018), Planelles et al. (2017), and Truong et al. (2018).
In Fig. 1 we show the calibration of the AGN feedback model used in the simulations. This is based on the correlation between the stellar mass of galaxies, M_{⋆}, and the mass of their central SMBHs, M_{BH}. In the figure the small lightblue points represent noncentral simulated galaxies identified by Subfind (Dolag et al. 2009), while darkblue dots represent simulated BCGs. The stellar masses of the BCGs are defined as the mass enclosed in a sphere of radius 0.1 × R_{500} around the position of the central SMBH, while total stellar masses of nonBCGs are given as an output by Subfind.
To calibrate the parameters for the AGNfeedback model we aimed at reproducing the entire M_{BH} − M_{⋆} relation including the majority of simulated nonBCG galaxies. These are, in particular, compared with the observational data from McConnell & Ma (2013) represented in the figure by the dashed line.
In the plot we also include other observational data, namely the BCGs from McConnell & Ma (2013) and the samples from Savorgnan et al. (2016), Main et al. (2017), Bogdán et al. (2018), and Gaspari et al. (2019). In Main et al. (2017) SMBH masses are computed from Kband luminosities using the relation log(M_{BH}/M_{⊙}) = − 0.38(±0.06)(M_{K} + 24)+8.26(±0.11) suggested by Graham (2007) and extracted from a sample of elliptical but not BCGs. In all the other works the mass of the SMBHs are derived from dynamical measurements. The BCG masses in McConnell & Ma (2013) and Bogdán et al. (2018) are part of a compilation from previous literature and we refer to the original papers and references therein for further information on the methods employed to infer the stellar masses. In Savorgnan et al. (2016) the stellar masses are computed from bulge luminosities assuming a constant masstolight ratio, while in Gaspari et al. (2019) it is assumed a variable M_{⋆}/L_{K} scaling as a function of the stellar velocity dispersion (see Gaspari et al. 2019 for details). It should be noted that from Savorgnan et al. (2016) we used only ellipticals, that are not necessarily BCGs. In Main et al. (2017) the stellar masses are computed from Kband luminosity using the relation log(M/L_{K}) = − 0.206 + 0.135(B − V) given by Bell et al. (2003).
Our main condition for deciding on the AGN parameters is that the observed correlation between SMBHs mass and nonBCG galaxies (dashed line) passes through the bulk of the simulated galaxies (small blue points). We also care for an overall agreement at the BCG scales but with less emphasis because of the scatter of the observed sample (McConnell & Ma 2013) is high at the high mass end. Regarding the BCGs, we find that the simulated BCGs are in a good agreement with observational data at both ends of the mass range, but that the simulated points tend to stay above observational data around M_{⋆} = 10^{12} M_{⊙}, although still inside their error bars. This discrepancy does not necessarily highlight a poor description of the simulations since several factors need to be considered for a proper comparison. First, the SMBH masses are computed by adopting different methods. For example, those extracted by Main et al. (2017) are calibrated using a relation that does not include BCGs and, indeed, they are more aligned with the nonBCG sample. Second, BCG stellar masses are computed using different apertures in simulations and in the various observational samples. Furthermore, measurements of stellar mass from different works can disagree due to the different assumptions made during the data analysis, such as the assumed initial mass function, the adopted stellar masstolight ratio, distances, and beam aperture. Dynamical measurements can significantly disagree in particular when considering different tracers, such as stars versus circumnuclear gas, including different gas phases as warm and ionized versus cold and molecular gas. The resulting differences among catalogs can be comparable to the separation between simulated and observed data points. An example is clearly represented in the figure by NGC 4889, the galaxy with the most massive SMBH. This object is present in the Savorgnan et al. (2016), McConnell & Ma (2013), and Bogdán et al. (2018) samples and, while M_{BH} is identical because taken from the same source in the literature, the estimated BCG mass can be different even by a factor of ∼2. This emphasizes the intrinsic difficulty in defining the BCG stellar masses and, at the same time, it quantifies a possible level of stellar mass difference among different works.
3. Method and samples
To investigate possible correlations between the central SMBHs and the global cluster properties we need to extract the cluster masses and temperatures from the simulated regions. In addition, we need to calculate the SMBH mass and the two contributions to its growth: the accretion of the surrounding gas and the merger with other SMBHs. In the following, after specifying the definition of the cluster center, we describe how all these quantities are computed.
Cluster center. As mentioned in the previous section, the SMBHs are always positioned at the location of the most bound particle that should identify the center of the host halo. Therefore, we followed the SMBHs back in time to identify the position of the hosting halo center. For this goal, we saved at z = 0 the unique identification number of the most massive SMBH particle which is within 10 kpc from the cluster minimum of the potential well, as identified by Subfind. Subsequently, we tracked it back in time to the epoch of its seeding. At each time, we checked that the SMBH is, as expected, at the minimum of the potential well of the hosting halo and not in a local minimum generated by merging substructures. With this approach, we built a merger tree of the central SMBHs rather than a merger tree of the clusters. We might expect that the two trees differ especially at early epochs (similarly to the small differences between the BCG and the cluster merger trees pointed out in RagoneFigueroa et al. 2018). However, we verified that for 80% of our systems the main progenitor of the SMBH is at the center of the main progenitor of the cluster up to z = 2 and for half of these the two trees coincide till the time when the SMBH is seeded.
Cluster masses. Once the center is defined as above, we considered the total gravitational mass of the cluster within an aperture radius R ≤ R_{500} computed by summing over all the species of particles: dark matter, cold and hot gas, stars, and black holes. At any redshift we considered only clusters with M_{cluster} = M_{500} ≥ 1.4 × 10^{13} M_{⊙} or, equivalently, with at least 20 thousands particles within R_{500}. The properties of the massselected sample are summarized in top part of Table 1.
Number of clusters, range of mass or temperature covered and their mean values for the mass sample (in the first half) and temperature subsample (in the second half).
BCG stellar mass. We defined the mass of the BCG, M_{BCG}, as the stellar mass enclosed in a sphere of radius 0.1 × R_{500} around the cluster center.
SMBH mass. Given the identification number of a SMBH particle, the mass of the SMBH, M_{BH}, at every redshift is quite easily retrieved from the simulation as it is the mass associated to that particle. The total mass of SMBH particles grows in time because of two separate phenomena: through the accretion of the diffuse gas or via BHBH mergers. In our simulations, these are the only possible channels for the SMBH to increase its mass. The accretion mass () is obtained by integrating the accretion rate, information that we saved at each time step. The merged mass () is simply calculated as a difference between the total mass and the accretion mass. As discussed later in the paper, the contribution to the SMBH mass by mergers is negligible at z ≥ 1.5. Therefore the analysis of this component is restricted to lower redshifts.
Temperature. In order to compare our results to those from XMMNewton observations we considered the spectroscopiclike temperature (Mazzotta et al. 2004):
where ρ_{i}, m_{i}, and T_{i} are the density, mass, and temperature of the ith gas particle within R_{500} emitting in the Xray band, that is with T_{i} > 0.3 keV and a cold fraction lower than 10%. We note, in this respect, that according to the effective model by Springel & Hernquist (2003) adopted in our simulations, the gas particles can be multiphase, carrying information on both hot and cold gas. The cold phase provides a reservoir for stellar formation. In order to have a reliable estimation of the temperature inside R_{500} we imposed two conditions: a minimum of 10^{4} hot gas particles and a maximum fraction of 5% of gas particles discarded because too cold. All clusters satisfying these requirements have also M_{500} > 1.4 × 10^{13} M_{⊙}, thus whenever we consider measurements of temperature we refer to a subsample of the mass selectedsample. The properties of the temperatureselected subsample are summarized in the bottom part of Table 1 and the analysis of this subsample is restricted to z ≤ 1 because at the highest redshift bins, z = 1.5 and z = 2, we do not have enough statistics to apply a meaningful analysis.
Bestfitting procedure. For all the considered relations, we looked for a bestfit line in the logarithmic plane:
where log always indicates the decimal logarithm. The temperature, cluster mass, BCG stellar mass, and SMBH mass are always normalized by the same factors, expressed above as X_{0} or Y_{0} and respectively equal to 2 keV, 10^{14} M_{⊙}, 10^{11} M_{⊙}, and 10^{9} M_{⊙}.
To find the bestfit curve, we employed an IDL routine that is resistant with respect to outliers: ROBUST_LINEFIT^{1}. For the simulated data, we always considered the BISECT option, recommended when the errors on X and Y are comparable so there is no true independent variable. This is particularly appropriate in the case of numerical simulations where no errors are linked to measurements. To estimate the error associated with the parameters of the bestfitting relation, we generated ten thousand bootstrap samples by randomly replacing the data. From the resulting distributions we derived the mean values and the standard deviations to be associated, respectively, with the parameters and their errors. All relevant bestfitting coefficients of the linear regressions are reported in Table 2 and will be discussed in the next two sections.
4. Results
4.1. Comparison with observational data
In this section we compare the numerical results to the observations presented in Bogdán et al. (2018), where the correlation between the mass of the SMBHs in BCG and the global temperature of clusters and groups of galaxies is presented. We also complement this analysis with the more recent data of Gaspari et al. (2019)^{2}. Bogdán et al. (2018) derive the temperature from XMMNewton observations of the hot gas; Gaspari et al. (2019) use published data of Chandra and widefield ROSAT/XMM Newton. The extraction region for both is on average ∼0.2R_{500} (core included; such temperature is a good proxy for T_{500}).
In Fig. 2 we show the correlation between M_{BH} and T_{500} for our simulations and for their observational data set. In the figure we also show as a dashed red line the results of the simple toy model by Churazov et al. (2005)^{3}, which even under simplified assumptions reproduces the observed correlation. We find a good agreement between observations and simulations. Nevertheless, a more quantitative comparison between the two samples is difficult for an underrepresentation of clusters with T_{500} > 2 keV in the observational sample that, however, has a good number of systems below 1 keV.
Fig. 2.
Correlation between SMBH mass, M_{BH}, and the clusters temperature, T_{500}. Red and orange crosses refer to observational data from Bogdán et al. (2018) and Gaspari et al. (2019) respectively. Darkblue dots represent simulated clusters in the temperature subsample while cyan stars show the remaining objects of the mass sample. Dashed red line is the prediction of the toy model by Churazov et al. (2005). 
In order to better populate the colder and less massive tail of the simulated cluster distribution we compare the correlation between M_{BH} and M_{500} by using the mass sample rather than the lessnumerous temperature subsample (Table 1). In Bogdán et al. (2018) the total mass was not measured directly from their data but was derived from the temperature via the scaling relation by Lovisari et al. (2015):
For this reason, before analyzing the M_{BH}–M_{500} relation, it is helpful to compare the observed and simulated T_{500}–M_{500} relations. The observed and simulated data sets are shown in Fig. 3 together with the bestfitting linear relations. In case of the observed sample we verified that our fitting procedure, without the BISECT option, returned the same parameters of Eq. (4). In particular we confirm the value of the slope reported in Lovisari et al. (2015): b = 1.65 ± 0.07.
Fig. 3.
Correlation between cluster mass, M_{500}, and cluster temperature, T_{500}. Symbols as in Fig. 2, where the observational data are taken from Lovisari et al. (2015). Dashed lines are the bestfitting lines for both simulated and observed data. 
By looking at Fig. 3, we see a good agreement between simulated and observed clusters in the temperature range covered by Lovisari et al. (2015). However, the extrapolation of their bestfit line suggests a possible difference in the hottestclusters regime. The two slopes agree within 1σ, but the observed clusters have on average slightly higher temperature with respect to simulated clusters at fixed mass. For example, the temperature of observed clusters is 9% higher at M_{500} = 10^{14} M_{⊙}. This feature is not new and has been already noted in Truong et al. (2018), where a similar set of numerical simulations is employed, and, more interestingly, in other numerical analysis, such as Henden et al. (2018). In particular, in their work the authors show that numerical results are in agreement with observations if are considered only cluster masses derived via weak lensing. This suggests that the observed Xray hydrostatic masses are biased low.
Finally, in Fig. 4 we compare the correlation between M_{BH} and the M_{500} as measured in our simulations and as derived in Bogdán et al. (2018) and Gaspari et al. (2019). The results of the comparison are expected from the previous two figures: the simulated data points are in line with observations, especially at high (M_{500} > 3 × 10^{14} M_{⊙}) and low masses (M_{500} < 3 × 10^{13} M_{⊙}). In the intermediate mass range, the few observed data points tend to have slightly higher SMBH masses than the simulated objects. This apparent mismatch is presumably a consequence of the poor statistics of 10^{14} M_{⊙} objects in the observational sample. More unlikely, this feature could suggest a broken power law to describe the M_{BH} − M_{500} relation, but such a drastic change in the slopes is difficult to justify.
4.2. The theoretical M_{BH}–M_{500} relation
Since the simulated sample is in overall good agreement with the observed correlation between the mass of the central SMBHs and the mass of the clusters, we investigate here how single simulated systems evolve throughout time to form, by z = 0, the M_{BH} − M_{500} relation shown in Fig. 5. For this goal, we overplotted three evolutionary tracks of representative SMBHs. These objects are chosen accordingly to their mass at z = 0; specifically, they refer to a small, a mediummass, and a massive SMBH. To have some temporal reference we also indicated four specific times along each line: z = 3, z = 2, z = 1, and z = 0.
Fig. 5.
Evolutionary tracks of three different systems on the M_{BH}–M_{500} plane. Triangles over each line indicate the position of the systems on the plane at z = 3, z = 2, z = 1 and z = 0. Black circles represent our numerical sample at z = 0; the filled ones are systems for which M_{500} is increased by more than 40% in the last Gyr. 
Despite the different final masses, the evolutionary tracks of the three systems have strong similarities which are common also in all the other objects analyzed (not shown for sake of clarity). Three phases are clearly distinguishable. At the highest redshifts, the mass of the SMBHs grows rapidly at almost constant M_{500}. This track begins instantaneously as the SMBHs are seeded in a gas rich region. The SMBHs immediately gain mass at the Eddington limit by the accretion from the abundant surrounding gas, which is mostly cold and thus efficient at increasing the SMBH mass. This phase typically lasts half Gyr and can lead to the formation of SMBHs with a mass already of the order of M_{BH} ≈ 10^{8} − 10^{9} M_{⊙}, in line with other hydrodynamical cosmological simulations (e.g., McAlpine et al. 2018). The fast SMBH growth ends when the M_{BH} is high enough to cause an intense feedback that leads to the ejection of part of the gas outside the shallow potential wells of the hosting galaxies. By then, all SMBHs in our sample are close to the M_{BH} − M_{500} relation. This happens before z = 2 and in some cases even at z > 4.
After this initial phase, the cluster and its central SMBH coevolve, but not with a linear evolutionary track: the increase of the SMBH and cluster masses is not simultaneous. The shape of the tracks, instead, can be described as a stairway: the systems evolve in this plane either at almost constant M_{BH} or at almost constant M_{500}. The former situation occurs during cluster mergers. It starts when merging structures reach and cross R_{500} leading to a quick increment of the total cluster mass and finishes when the secondary objects are fully incorporated. These horizontal shifts in the M_{BH}–M_{500} plane typically last 1 Gyr or less and only in the rare cases when multiple mergers are subsequently taking place they can last up to 2 Gyr. In the following period, spanning from 1 to 3 Gyr, the substructures move towards the center of the cluster and neither the SMBH mass nor the cluster mass change. Eventually, the merging objects reach the core and either feed the central SMBH with gas or induce a BHBH merger or both. The event is captured by the vertical movement in the plot.
All these timeframes are clearly indicative as they depend on several parameters that characterize the merger events such as the mass ratio and the impact parameter. Nonetheless, it is always the case that the mass of the SMBH and of the cluster are for the largest majority of time at the connection between the horizontal and vertical steps rather than along their tracks. This behavior indicates that the scatter of the relation might differ when the sample is selected according to the dynamical status of the SMBH hosts. We expect that the points related to relaxed BCG in relaxed clusters will always be above the points linked to systems where either the BCG or the cluster are experiencing, or just experienced, a merger event. Indeed, we expect that relaxed and perturbed systems are respectively located in our plot on the top and the bottom of the vertical segments. In favor of this picture, it is noticeable that the 20 clusters whose M_{500} grows by more than 40% in the last Gyr (shown as black points in Fig. 5) have SMBHs that on median are 50% smaller than those expected to follow the M_{500}–M_{BH} relation.
4.3. Evolution of the M_{BH}–M_{500} relation
After the inspection on the trajectory of individual objects we study here the evolution of the entire M_{500}–M_{BH} relation. We start with the evaluation of the ratio between the mass gained by clusters and by central SMBHs between z = 2 and z = 0: Δ_{M} = M_{z = 2}/M_{z = 0}, where M refers to either the total mass, M_{500}, or the SMBH mass, M_{BH}. If these two ratios are constant, the slope of the relation will not change. The resulting ratios are shown in Fig. 6 as a function of the cluster total mass reached at z = 0 for the mass sample identified at z = 2. From the plot, we can clearly infer that the variation in total mass between the two epochs is strongly mass dependent. The absolute value of the slope of the bestfitting Δ_{M500} − M_{500} relation is, indeed, greater than 0.75. Clusters with a final mass lower than 10^{14} M_{⊙} increase their total masses by a factor between 3 and 6. Instead, massive clusters with final M_{500} > 10^{15} M_{⊙} increase their mass on average by a factor of about 30 with individual objects that can grow by more than two orders of magnitude. This feature is completely in line with the expectations of hierarchical clustering.
Fig. 6.
Ratios between cluster (SMBH) masses at z = 2 and cluster (SMBH) masses at z = 0 as a function of the cluster masses at z = 0. Clusters are shown as green squares and SMBHs as yellow triangles). The sample used is the massselected sample identified at z = 2. 
The high mass regime is particularly characterized by a large spread that is representative of what we might expect from a volumelimited sample because the most massive objects, M_{500} > 10^{15} M_{⊙} correspond to the most massive systems of the parent volumelimited cosmological box (see Sect. 2). Vice versa, the scatter for the smallest systems is likely underestimated. Indeed, the linear trend is expected to flatten for the lowest masses to a constant growth rate value. On the other side, Fig. 6 also shows that the variation on the SMBHs mass is independent of the cluster mass and that SMBHs grow on average by a factor of about 5−6. As a consequence we expect a marked evolution of the slope of the M_{BH}–M_{500} relation between z = 2 and z = 0.
This is confirmed in Fig. 7 where we plot the bestfitting lines for our mass samples at all redshifts considered. The relations are steeper at higher redshifts: the value at z = 2, b = 1.348, is almost twice the value found at z = 0, b = 0.753. From Figs. 6 and 7 we conclude that the change in the slope is mainly driven by the different evolution rate of the most massive clusters with respect to the smallest objects, trend that is in line with the expectations from ΛCDM cosmology.
Fig. 7.
Correlation between M_{BH} and M_{500} at different redshifts. At every redshift we show the bestfitting relation in the mass range of the respective samples. Namely, we show in green, orange, red, magenta, and blue the mass sample related to z = 2, 1.5, 1, 0.5 and z = 0, respectively. 
4.4. Evolution of SMBH mass
To better understand how the black hole mass evolves with time we separately study the two mechanisms that contribute to the growth of the SMBH mass: the accretion by the diffuse gas, , and the merger with other SMBHs, . First, we analyzed how the mass of single SMBH evolves with time via the two separate channels. As an example in Fig. 8 we plot the evolution of the three SMBHs shown in Fig. 5. As already noted before, the evolution is characterized by an initial phase of intense gas accretion that for these three systems is approximately between z = 5 and z = 3. After that, SMBHs still grow by gas accretion but at a much smaller rate. On the contrary the increase of the SMBH mass due to BHBH mergers becomes more important and it is the main channel of the SMBH mass growth at lower redshifts, that is z ≤ 1 for the two most massive SMBHs and z ≤ 0.5 for the smallest one.
Fig. 8.
Evolution of the mass of black holes showed in Fig. 5. The solid lines represent the mass of the black holes due to gas accretion. Dashed lines represent mass gained via BHBH mergers. 
In Fig. 9 we show the evolution of our complete sample. We plot the median behaviors with a solid line and the 68% of the total sample distributions (from the 16th to the 84th percentiles) with the shaded regions. The total SMBH mass and the masses gained from the two channels are normalized with respect to the total SMBH mass at z = 0. Finally, the dashed vertical lines help to identify three significant times: z = 0.5, 1, and 2.
Fig. 9.
Evolution of SMBH mass divided into two growth channels (gas accretion in red and BHBH mergers in blue) considering the complete sample. In black we also plotted the total SMBH mass. Solid lines represent median values of our sample and shadowed regions represent 16 and 84 percentiles. 
From the SMBH seeding up to z ≈ 2 the total mass of the SMBH grows almost entirely by gas accretion. Half of the final mass gained through gas accretion is, indeed, accumulated before z = 2. Then from z ∼ 2 to z ∼ 1 the mass growth due to BHBH merger becomes more relevant increasing at a rate comparable to the growth rate of . By z = 1 makes up on average 25% of the total mass at that redshift. At lower redshift, z < 0.5, BHBH mergers provide the main channel for SMBH mass growth, and eventually represents the main component of mass gained by z = 0, in line with some previous results (e.g., Volonteri & Ciotti 2013; Dubois et al. 2014; Weinberger et al. 2018). Indeed, the mass accumulated by gas accretion from z = 1 to z = 0 accounts only for 10% of the total final mass, while during the same period the SMBH gains 50% of its final mass via mergers.
The relative importance of the two channels shown in the figure is, however, characterized by a large scatter. We, therefore, explored whether this is due to the broad mass range investigated and, thus, whether the described behavior depends on the mass of the systems. We divided the sample in three mass bins: at z = 0 the least massive objects have M_{500} below 10^{14} M_{⊙}, the most massive above 10^{15} M_{⊙}, and the intermediate in between these two thresholds. We found that the trends of the relative ratio of the two SMBH growth channels are extremely similar to Fig. 9 for the samples of the smallest and intermediate objects. This result is expected for the first mass bin since it is the most numerous, containing 84 objects, but it was not guaranteed for the intermediate sample with only 31 systems. The most massive sample, however, is on average characterized by a continuous and equivalent growth of both channels after z = 2 (see Fig. 10). As a result, at z = 0 accounts for ≈60% of M_{BH}. That said, we also notice that the scatter remains very large and the distributions related to the two channels show a large intersecting area. We, therefore, can conclude that the scatter shown in Fig. 9 is not related to the total mass of the systems.
The stronger relative influence by the SMBH accretion with respect to the BHBH merger in most massive clusters can be due to two different factors: on the one hand AGN feedback is not able to completely balance gas cooling, on the other hand BHBH mergers could be less frequent. In the following we show some evidence that both phenomena are actually in place.
To demonstrate that the AGN are less efficient in regulating the gas cooling in the cores of massive clusters we computed the total energy released by AGN feedback for z < 1 and related it to the gas mass within 0.1 × R_{500}. We find that the ratio of the two quantities is a strongly decreasing function of M_{500} and that it changes by more than a factor of 10 from the least massive to the most massive systems. This suggests that the heating provided by the AGN feedback is relatively smaller for large objects where, therefore, the gas cooling is less contrasted. The central SMBH has therefore more cold gas at its disposal.
To test the reduced frequency of BHBH merger, we computed at z = 0 the number and mass of SMBHs which are inside 0.1 × R_{500} and are not bound to any substructure. We find that 80% of the most massive systems (16 objects over 20) have several SMBHs in that central regions with a total mass greater than 10% of the mass of the central SMBH. Analysing the first mass bin (M_{500} < 10^{14} M_{⊙}), instead, only 10% of clusters have enough SMBHs able to account, all together, for at least 10% of the mass of the central SMBH. This is mostly due to the fact that more massive clusters host more massive and extended BCGs. When the substructures interact with these wellestablished BCGs they are more easily disrupted (see Sect. 4.5 for further details) at a larger radii, preventing or delaying mergers between their SMBHs and the central SMBH.
4.5. Recent growth of SMBH and stellar component
In our simulations the SMBH mass increases by a factor of ∼2.5 between 1 > z > 0. RagoneFigueroa et al. (2018) found instead, selecting the most massive cluster in each Lagrangian region described in Sect. 2, that the central stellar component measured within 30 kpc features a significant smaller growth. We checked that the growth factor of central SMBHs remain unchanged also considering this reduced subsample. This difference is due to the fact that in our simulations many substructures colliding with the BCG at z < 1 are largely disrupted and their stars quickly become part of the intra cluster light (ICL) or settle in the outermost radii of the BCG itself. This feature confirms some previous results. For example, Murante et al. (2007) shows that the bulk of the star component of the ICL originates during the assembly of the most massive galaxies in a cluster (and, in particular, of the BCG) after z ∼ 1.
To visualize this effect, whose detailed study will require a dedicated analysis, we simply compared the evolution of the SMBH mass and the inner stellar component, defined as M_{⋆,in} = M_{⋆}(r/kpc < 30). Furthermore, we added a measure of the outer stellar component, defined as M_{⋆,out} = M_{⋆,not−bound}(30 < r/kpc < 100). The “notbound” identification specifies that we excluded all stars which are gravitationally bound to substructures identified by Subfind and different from the BCG. M_{⋆,out}, therefore, comprises both the ICL and the outermost stellar mass of the largest BCG in the sample. The median values of ΔM_{⋆,in}, ΔM_{BH} and ΔM_{⋆,out} are computed after the normalization to their respective mass values at z = 0.
We considered 30 kpc for M_{⋆, in} because we wanted to evaluate the changes on the stellar component in the immediate surrounding of the SMBH; furthermore, this was used in RagoneFigueroa et al. (2018) as one of the possible definition of the BCG mass. For the outer component, instead, we also considered the mass of the unbound stars measured in other two spherical shells: between 100 kpc and 200 kpc and between 50 kpc and 350 kpc. As clear from Fig. 12, the choice of this region does not substantially change our conclusions: while M_{⋆,in} slowly increases from z = 1 to z = 0, the SMBH mass and rapidly increase. At z = 1, indeed, the quantities are about 80%, 45%, and 35% of their final values, respectively. The remarkable agreement between the two extra definitions of the outer stellar component – M_{⋆, not − bound}(100 < r/kpc < 200) and M_{⋆, not − bound}(50 < r/kpc < 350), thin blue solid and dashed line in Fig. 12 – implies that the growth rate of the ICL is independent on the specific radius used.
The final emerging picture is that many small substructures actually reach the cluster core and merge with the BCG. However, few of their stars remain in the innermost region. The interaction with the BCG causes that most stars of the structures are tidal shocked and stripped. Subsequently, they become gravitationally unbound thereby taking part of the ICL. During the disruption of the substructures, their most massive SMBHs feel the gravitational attraction of the underlying potential and sink towards the minimum of the cluster potential contributing to the growth of the SMBH at the center of the BCG. We note, however, that the modeling of BHBH mergers is very simplistic and could overestimate the efficiency of this physical process which take place at a scale well below the gravitational softening of the simulations.
4.6. The M_{BH}–M_{500} relation for the two SMBH growth channels
Given that at z = 1 and z = 0 the SMBHs have grown from both channels (through gas accretion and BHBH mergers), it is relevant to check whether the SMBH mass of the two channels are separately both related to the total mass or whether only one exhibits a tight correlation while the other mostly contributes to increase the scatter. This possibility is investigated in Fig. 11 where we considered both and as a function of M_{500} at z = 1 (left panel) and at z = 0 (right panel). As we have seen, at z = 1 the gas accretion is the dominant channel. In Sect. 4.4, we saw that this channel shows only a slight increment from z = 1 to z = 0. For this reason, the red points, referring to are substantially unchanged in the two panels. Vice versa, becomes the dominant component at z = 0. The results of the linear fits of the relations are reported in Table 2 also for the other redshifts.
Fig. 11.
Correlation between SMBH mass and M_{500} at z = 1.0 (left panel) and z = 0 (right panel). We plotted in blue the mass gained by mergers and in red the mass gained by gas accretion. Dashed lines represent best fitting lines to the data. 
Fig. 12.
Evolution of M_{⋆,in}, defined as the stellar mass inside a spherical region of radius 30 kpc, SMBH mass, and M_{⋆,out}, defined as the total mass of the stars enclosed in a spherical shell with radii 30 kpc and 100 kpc. All the quantities are normalized to their respective values at z = 0. Solid lines represent median values and shadowed regions represent 16 and 84 percentiles. The three different blue lines represent three definition of M_{⋆,out}. In particular the thin blue solid line is the stellar mass in a spherical shell with radii 100 kpc and 200 kpc while the dashed blue line is the stellar mass in a spherical shell with radii 50 kpc and 350 kpc. 
From the figure, it is evident that both masses correlate well with M_{500} at both redshifts independent of which one of the two is the dominant channel from the SMBH mass growth. From the table, we notice that the slopes of the two relations are consistent within 1 σ being the z = 1 slightly steeper as expected from the Sect. 4.3. Most importantly, the two scatters are similar and both slightly higher than the scatter of the relation of the total SMBH mass (see Table 2).
Finally, we would like to emphasize that Fig. 7 suggests that the M_{BH} − M_{500} relation is already in place at z = 2 when the SMBH mass was almost entirely gained only by gas accretion. These results enlighten that, at least in our simulations, mergers are not essential to establish the relation at first.
5. Discussion
As previously remarked the correlation between the BCG mass and the SMBH mass has been diffusely studied and often used to extract the mass of the SMBH knowing the mass of the hosting galaxy (Sect. 1). In Bogdán et al. (2018), the authors found in their observed sample that the scatter between the SMBH mass and the global cluster properties is tighter by almost 40% than the scatter of the M_{BH} − M_{BCG} relation. Indeed, in their Table 4 they report σ_{MBHMBCG} = 0.61 and σ_{MBHT} = 0.38. The values of the scatters increase in the later analysis by Phipps et al. (2019) because of the method used to derive M_{BH}, the authors suggest. Despite, also in that case the two intrinsic scatters in M_{BH} are comparable between each other. As a consequence, the global properties of cluster within R_{500} are also suitable to estimate the SMBH mass. We show in Sect. 4 that the distribution of our simulated data is in reasonable agreement with the observations by Bogdán et al. (2018) and Gaspari et al. (2019), who used dynamical measurements of M_{BH}. We demonstrate that for our simulated objects there is a clear correlation between the mass of the SMBHs, located at the minimum of the potential well, and the temperature or total mass of the clusters within R_{500}. In this section we discuss the properties of all the relations described in the paper and listed in Table 2.
To this end we refer to Fig. 13 where we show the covariance matrix between the deviations of all the quantities of interest from their bestfitting relations, that are their residual at fixed mass. These are defined as logarithmic differences between the actual value, generically referred as X, and the expected value, X_{FIT}, from the X − M_{500} relations^{4} provided in Table 2:
Fig. 13.
Covariance matrix. Each axis represents the logarithmic difference between the actual value of a quantity X and the expected value from the linear relation (XM_{500}) at its fixed mass, as defined by Eq. (5). Panels above the diagonal refer to z = 1 while panels below the diagonal refer to z = 0. The diagonal panels show the distribution of δ(X) at z = 0. Red points define clusters which have experienced a mass growth of at least 40% during the last Gyr. 
The deviations are computed for each quantity X = T_{500}, M_{BCG}, M_{BH}, , . The panels above the diagonal refer to z = 1 while those below to z = 0. The diagonal panels show the distribution of δ(X) at z = 0. We used the temperature subsample whenever δ(T_{500}) is considered. For each pair of deviations we listed the Spearman correlation coefficients in Table 3 computed at z = 0, 0.5, 1 and we marked in bold the correlations whose probability to be consistent with zero is less than 2 per thousand and its module is greater than 0.4. A strong correlation between δ(M_{BH}) and δ(X) is converted into a small scatter in the relation M_{BH} − X.
Correlation matrix.
As previously commented, the scatter σ_{MBHT500} is comparable to σ_{MBHM500} (see Table 2). As a consequence, the temperature and the total cluster mass are equally good proxy for the SMBH mass. The similarity between these two proxies can be explained by looking at the correlation between δ(T_{500}) and δ(M_{BH}) in Fig. 13. The panels for the two considered redshifts, z = 0 and z = 1, highlight how the variations of cluster temperature are not directly reflected into variations of the SMBH mass. At first sight, this can be surprising as one could expect that both quantities are strongly dependent on the dynamical activities of the cluster core and, especially, sensitive to merger events that impact the innermost region of the clusters. However, the response of the temperature and the SMBH mass to merger is not simultaneous: we saw in Sect. 4.2 that the typical delay between the increase of the M_{BH} after a merger can be around 1−3 Gyr, while the temperature increase typically occurs in less than 1 Gyr from the merging episode. In addition, the rapid increase of the ICM temperature is followed by a small drop caused by the expansion of the shock towards the more external regions. This temperature oscillation along with the mismatch between the two timedelays explain why the correlation coefficients between δ(T_{500}) and all M_{BH} components are always low and with a high probability to be consistent with zero. These characteristics are also in place when the temperature variations are compared with the variation on the SMBH mass due to BHBH mergers.
When comparing, instead, the scatter of the relation between the SMBH mass and the cluster properties with the scatter of the M_{BH} − M_{BCG} relation, we find that σ_{MBHMBCG} is smaller at all times. At redshift z = 0, for example, σ_{MBHM500} is ≈25% larger, although the two scatters are in agreement between 2σ. The difference between the scatters is more significant at z ≥ 1 where typically σ_{MBHM500} is ≈1.4 × σ_{MBHMBCG}. In the covariance matrix formalism this translates into a correlation between δ(M_{BCG}) and δ(M_{BH}). Indeed, at z ≤ 1 the correlation is ≥0.55. Moreover, the correlation is stronger if computed with respect to when it reaches values around 0.7. As we can see from Fig. 13, the correlation at z = 0 is significant for the presence of situations of either premergers or mergers with a large impact parameter when both M_{BCG} and M_{BH} (and in particular ) are smaller than the average of the sample at fixed total mass (e.g., their δ is negative). At z = 1, we also notice a correlation between the variations of the two quantities in the other direction (both δ(M_{BCG}) and greater than zero) implying that the SMBHs that gained mass through z ≥ 1 mergers are hosted in BCG with higher stellar mass with respect to the average of the sample.
In Fig. 5, we show that the scatter of the M_{BH} − M_{500} relation is influenced by the presence of the systems that recently experienced a major merger. Indeed, all objects that increased their total mass by at least 40% in the last Gyr are in the bottom part of the overall distribution (see Sect. 4.2 and Fig. 5). In Fig. 13 we identify these objects as red points. In the majority of the cases their variations with respect to the mean are negative. When we exclude these objects and reapply our fitting procedure^{5}, we find that the scatter of the M_{BH} − M_{500} relation reduces by almost 20% (σ_{MBHM500} = 0.14 ± 0.01), it remains comparable to the recomputed scatter of the M_{BH} − T_{500} relation (σ_{MBHT500} = 0.14 ± 0.02) and consistent with the new σ_{MBHMBCG} = 0.13 ± 0.01. In other words, even if the removal of the dynamically active objects induces a decrease in the scatters of all of the three relations, the most important reduction impacts the scatter at fixed total mass, reducing the gap between σ_{MBHMBCG} and the scatter of the relations involving global cluster properties.
In the interpretation of these results from simulations, we need to recall that the simulated data do not reproduce the observed scatter of the M_{BH} − M_{BCG} relation (Fig. 1 and also RagoneFigueroa et al. 2013; Bogdán et al. 2018). On one hand the growth of simulated SMBHs is regulated by simplistic subgrid models that do not capture all physical processes in place and might lead to a reduced scatter. On the other hand, as explained when discussing Fig. 1, a large portion of the observed scatter around the M_{BH} − M_{BCG} relation can be ascribed to observational uncertainties associated either with the quantity definition (e.g., treatment of intracluster light, BCG boundary definition) or with the measurement procedures (e.g., not fixed aperture mass for the BCG or application of scaling relation to infer the SMBH mass). In simulations, instead, the BCG and SMBH masses are always known and precisely defined (see Sect. 2 on the discussion on the BCG mass determination). These arguments not only provide a possible explanation for the difference between the simulated and observed scatters but also underline that the errors on the measures of M_{BCG} derived from observations are not easily reducible. The estimate of the SMBH mass from the BCG mass can always be subject to these uncertainties. The global cluster properties are also subject to systematics, which however can be treated as follows. A systematic bias on the global temperature can be dealt with precise instrument calibration or with multitemperature fitting. The uncertainties on the total mass are reduced when measurements coming from various wavelengths are combined, such as mass reconstruction from gravitational strong and weak lensing, galaxy dynamics, SZ, and Xray. These considerations, along with the limited difference in the relation scatters, emphasize how the global cluster properties can be powerful proxies for the SMBH mass. This conclusion is even stronger at high redshift, such as z = 1.
6. Conclusion
In this paper we studied the correlation between the mass of SMBH at the center of BCG and global properties of the hosting cluster of galaxies, namely its total mass, M_{500}, and global temperature, T_{500}. Our work is based on 29 zoomin cosmological hydrodynamical simulations carried out with GADGET3, a modified version of the public code GADGET2. This code treats unresolved baryonic physics including AGN feedback through various subgrid models. The parameters used to model the AGN feedback are tuned to appropriately reproduce the scaling relation between the masses of the SMBHs and their hosting galaxy (Fig. 1). For this study, we considered all systems with M_{500} > 1.4 × 10^{13} M_{⊙} identified in the highresolution regions of the resimulated volumes. At z = 0, there are 135 objects.
After showing the agreement between our numerical results and the observational data, we explored how the relation between the SMBH mass and the cluster mass establishes by looking at the coevolution of these quantities in individual systems. We then looked at the evolution of the entire sample considering four different times (z = 0.5, 1, 1.5, and 2). Finally, we characterized the role played by the two channels for the SMBH growth: accretion of gas and BHBH merging. Our main results can be summarized as follows.

The simulated relations between SMBH mass and the global cluster properties (M_{BH}–M_{500} and M_{BH} − T_{500}) are in agreement with the observations by Bogdán et al. (2018) and Gaspari et al. (2019) (Figs. 2 and 4).

The M_{BH} − M_{500} relation at z = 0 originates from a nonsimultaneous growth of the SMBH and of the cluster. In particular, objects evolve on the M_{BH} − M_{500} plane either at almost constant M_{BH} or at almost constant M_{500} (Fig. 5). The rapid increase of the cluster total mass occurs during cluster mergers and typically last 1 Gyr. Subsequently, the substructures move towards the cluster center and, eventually, reach the core feeding the central SMBH with gas and/or inducing a BHBH merger. Clusters that recently experienced major merger events are in general below the mean relation.

The M_{BH} − M_{500} relation of the entire sample shows a degree of evolution: the slope is about 45% smaller at z = 0 with respect to z = 2 (Fig. 7). This evolution naturally arises from hierarchical structure growth as the most massive clusters increase their mass in the period between z = 2 and 0 at a rate which is 10−20 times higher with respect to the smallest objects. Vice versa, the SMBH mass increases by an approximately constant factor (≈5), independent of the mass of the hosting cluster (Fig. 6).

In our simulations, SMBHs grow by two different channels. Gas accretion is the most relevant channel at redshift z > 2 and the only player at the earliest times. The accretion is slowed down only when the SMBHs are massive enough to balance gas cooling via AGN feedback. At lower redshift (z = 1) one quarter of the SMBH mass is ascribed to mergers. From that time to z = 0 the BHBH merger contribution becomes progressively more important. Indeed, mergers contribute by about 60% of the total z = 0 SMBH mass on average. When restricting the analysis to the most massive systems, we find that the accretion onto the SMBH is dominant up to z = 0, both for a reduced power of the AGN heating over the gas cooling and for a less frequent BHBH merger rate. and similarly relate to M_{500} at both z = 0 and z = 1, independently of which SMBH mass growth channel is dominant.

The rapid increase of the SMBH mass at z < 1 due to mergers is much faster than that of the hosting BCG, defined as the stellar mass within the fixed physical aperture of either 30 or 50 kpc and studied in details by RagoneFigueroa et al. (2018). Indeed, in our simulations, we found that at recent times (z < 1) galaxies are frequently merging with the BCG but they get easily disrupted. This process rather than transferring the merging stellar mass to the BCG originates and feeds the diffuse stellar component (the intracluster light), that we simply define here as unbound stars located outside the BCG.

The M_{BH} − M_{500} and M_{BH} − T_{500} relations present a similar scatter (or, equivalently, δ(T_{500}) does not show any significant correlation with δ(M_{BH})), meaning that they are equally valid SMBH mass proxy. On the other hand, δ(M_{BH}) is highly correlated with δ(M_{BCG}). As a consequence the scatter of the M_{BH} − M_{BCG} relation, at z = 0, is ≈25% lower in our simulations, although the two values are in agreement within 2σ. However, it is important to stress that the observed scatter of the M_{BH} − M_{BCG} relation is larger then the simulated one mainly for two reasons: the numerical limitation of a simplistic description of SMBHs growth and the large uncertainties affecting the observational measurements of both SMBH and BCG mass. That said, when the most dynamically active objects are discarded from the sample, the scatters of all relations, M_{BH} − M_{BCG}, M_{BH} − M_{500} and M_{BH} − T_{500}, become similar, thus strengthening the predicting power of the cluster global quantities.
Even though our simplified subgrid models are effective in reproducing the observed relations, before concluding, it is important to list some limitations of our current configuration, which we aim at improving in future work. As first point our resolution is sofar limited to about 5 kpc to obtain a large statistical sample in a cosmological environment. This means that while the largescale inflows can be correctly described, satellite galaxies are resolved with few smoothing lengths implying that the actual level and timing of gas accretion and SMBH mergers might differ from our tracked evolution. Resolving satellite galaxies with less than ∼250 particles results in reducing their concentration and ability to resist to disruption events (van den Bosch & Ogiya 2018). This might be one of the causes of a slightly reduced number of galaxies with M_{⋆} ∼ 10^{11} M_{⊙} in our simulations with respect to the field stellar mass function as measured by Bernardi et al. (2013) (see Fig. B.1). Similar problems affect other simulations (Steinborn et al. 2015) and have been discussed in connection to the mismatch between the substructure distribution measured in observations and simulations (e.g., Grillo et al. 2015; Munari et al. 2016). Explanations for the discrepancy between the observed and simulated stellar mass functions can also be of physical nature rather than exclusively numerical. For example, a more physically motivated parametrization of the AGN feedback could play a relevant role.
Several studies (Gaspari 2016 for a brief review) show that the SMBH accretion in massive galaxies proceeds via chaotic cold accretion (CCA), from the kiloparsec scale down to tens gravitational radii. While our subgrid feeding module tends to mimic CCA in power, the rapid frequency of chaotic events is not fully captured here. CCA leads to the rapid funneling of cold clouds also at low redshift, as the hot halo develops nonlinear thermal instability and a consequent rain of clouds toward the inner SMBH. Such frequent accretion, which is key to quench cooling flows, would likely extend the dominance of gas feeding over mergers even at low z. Our injection of AGN feedback is also simplistic, since we only inject thermal energy within a fixed aperture with a fixed efficiency. The physics of feedback is quite more complex: massive outflows and radio jets inflate bubbles and drive shocks in the ICM, while CCA seems to occur in perpendicular cones (e.g., Gaspari et al. 2018; Yang et al. 2019). This leads to highly variable duty cycle and variations in the core ICM properties (temperature, entropy, etc.). Finally, we can improve the predictive power of future simulations by using physically motivated parameters, such as the horizon mechanical efficiency given in GRMHD BH accretion simulations (e.g., Sdowski & Gaspari 2017).
Equation (7) of Churazov et al. (2005) with δ_{E} = 10^{−4}, t_{9} = 1, Λ(T) as defined in Tozzi & Norman (2001) with Z = 0.3 Z_{⊙}, and T = 3 × 10^{6} σ_{200} K.
The bestfitting relations of the samples derived by excluding the clusters with recent fast accretion are: M_{BH}/(10^{9} M_{⊙}) = 10^{0.49 ± 0.02} × (M_{500}/10^{14} M_{⊙})^{0.77 ± 0.03}; M_{BH}/(10^{9} M_{⊙}) = 10^{0.56 ± 0.02} × (T/2 keV)^{1.32 ± 0.06}; M_{BH}/(10^{9} M_{⊙}) = 10^{−0.41 ± 0.03} × (M_{BCG}/10^{11} M_{⊙})^{1.17 ± 0.04}.
Acknowledgments
We thank the anonymous referee for the careful and constructive reading of the paper and for his/her useful suggestions. We would like to thank Marta Volonteri and Lorenzo Lovisari for helpful feedback during the draft of the paper and Volker Springel for making the GADGET3 code available to us. This project has received funding from: ExaNeSt and Euro Exa projects, funded by the European Union’s Horizon 2020 research and innovation program under grant agreement No 671553 and No 754337, the agreement ASIINAF n.201714H.0, the Agencia Nacional de Promoción Científica y Tecnológica de la República Argentina under the PICT24172013 grant; 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’s Horizon 2020 Research and Innovation Programme under the Marie SklodowskaCurie grant agreement No 734374, PRINMIUR 2015W7KAWC, the INFN INDARK grant, NASA Chandra GO718121X and GO819104X. Simulations have been carried out using MENDIETA Cluster from CCADUNC, which is part of SNCADMinCyT (Argentina) and MARCONI at CINECA (Italy), with CPU time assigned through grants ISCRA C, and through INAFCINECA and University of Trieste – CINECA agreements. The postprocessing has been performed using the PICO HPC cluster at CINECA through our expression of interest.
References
 Bandara, K., Crampton, D., & Simard, L. 2009, ApJ, 704, 1135 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110 [NASA ADS] [CrossRef] [Google Scholar]
 Beifiori, A., Courteau, S., Corsini, E. M., & Zhu, Y. 2012, MNRAS, 419, 2497 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697 [NASA ADS] [CrossRef] [Google Scholar]
 Biffi, V., Planelles, S., Borgani, S., et al. 2017, MNRAS, 468, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Biffi, V., Planelles, S., Borgani, S., et al. 2018, MNRAS, 476, 2689 [NASA ADS] [CrossRef] [Google Scholar]
 Bogdán, Á., Lovisari, L., Volonteri, M., & Dubois, Y. 2018, ApJ, 852, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Bonafede, A., Dolag, K., Stasyszyn, F., Murante, G., & Borgani, S. 2011, MNRAS, 418, 2234 [NASA ADS] [CrossRef] [Google Scholar]
 Booth, C. M., & Schaye, J. 2009, in AIP Conf. Ser., eds. S. Heinz, & E. Wilcots, 1201, 21 [NASA ADS] [Google Scholar]
 Booth, C. M., & Schaye, J. 2010, MNRAS, 405, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Brough, S., Couch, W. J., Collins, C. A., et al. 2008, MNRAS, 385, L103 [NASA ADS] [CrossRef] [Google Scholar]
 Churazov, E., Sazonov, S., Sunyaev, R., et al. 2005, MNRAS, 363, L91 [NASA ADS] [Google Scholar]
 Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497 [NASA ADS] [CrossRef] [Google Scholar]
 Dubois, Y., Volonteri, M., & Silk, J. 2014, MNRAS, 440, 1590 [NASA ADS] [CrossRef] [Google Scholar]
 Fabian, A. C. 1999, MNRAS, 308, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrarese, L. 2002, ApJ, 578, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M. 2016, in Galaxies at High Redshift and Their Evolution Over Cosmic Time, ed. S. Kaviraj, IAU Symp., 319, 17 [NASA ADS] [Google Scholar]
 Gaspari, M., & Sdowski, A. 2017, ApJ, 837, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M., Eckert, D., Ettori, S., et al. 2019, ArXiv eprints [arXiv:1904.10972] [Google Scholar]
 Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Graham, A. W. 2007, MNRAS, 379, 711 [NASA ADS] [CrossRef] [Google Scholar]
 Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Grillo, C., Suyu, S. H., Rosati, P., et al. 2015, ApJ, 800, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198 [NASA ADS] [CrossRef] [Google Scholar]
 Häring, N., & Rix, H.W. 2004, ApJ, 604, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Henden, N. A., Puchwein, E., Shen, S., & Sijacki, D. 2018, MNRAS, 479, 5385 [NASA ADS] [CrossRef] [Google Scholar]
 Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, J. 2008, MNRAS, 386, 2242 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, J. 2009, MNRAS, submitted [arXiv:0908.2028] [Google Scholar]
 Jahnke, K., & Macciò, A. V. 2011, ApJ, 734, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J., & Bender, R. 2011, Nature, 469, 377 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Lakhchaura, K., Truong, N., & Werner, N. 2019, MNRAS, 488, L117 [CrossRef] [Google Scholar]
 Lin, Y.T., & Mohr, J. J. 2004, ApJ, 617, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Lovisari, L., Reiprich, T. H., & Schellenberger, G. 2015, A&A, 573, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285 [NASA ADS] [CrossRef] [Google Scholar]
 Main, R. A., McNamara, B. R., Nulsen, P. E. J., Russell, H. R., & Vantyghem, A. N. 2017, MNRAS, 464, 4360 [NASA ADS] [CrossRef] [Google Scholar]
 Marconi, A., & Hunt, L. K. 2003, ApJ, 589, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10 [NASA ADS] [CrossRef] [Google Scholar]
 McAlpine, S., Bower, R. G., Rosario, D. J., et al. 2018, MNRAS, 481, 3118 [NASA ADS] [CrossRef] [Google Scholar]
 McConnell, N. J., & Ma, C.P. 2013, ApJ, 764, 184 [NASA ADS] [CrossRef] [Google Scholar]
 McLure, R. J., & Dunlop, J. S. 2002, MNRAS, 331, 795 [NASA ADS] [CrossRef] [Google Scholar]
 Merloni, A., Heinz, S., & di Matteo, T. 2003, MNRAS, 345, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D., & Ferrarese, L. 2001, ApJ, 547, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Mittal, R., Hudson, D. S., Reiprich, T. H., & Clarke, T. 2009, A&A, 501, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Munari, E., Grillo, C., De Lucia, G., et al. 2016, ApJ, 827, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Murante, G., Giovalli, M., Gerhard, O., et al. 2007, MNRAS, 377, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Peng, C. Y. 2007, ApJ, 671, 1098 [NASA ADS] [CrossRef] [Google Scholar]
 Phipps, F., Bogdán, Á., Lovisari, L., et al. 2019, ApJ, 875, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Planelles, S., Fabjan, D., Borgani, S., et al. 2017, MNRAS, 467, 3827 [NASA ADS] [CrossRef] [Google Scholar]
 RagoneFigueroa, C., Granato, G. L., Murante, G., Borgani, S., & Cui, W. 2013, MNRAS, 436, 1750 [NASA ADS] [CrossRef] [Google Scholar]
 RagoneFigueroa, C., Granato, G. L., Ferraro, M. E., et al. 2018, MNRAS, 479, 1125 [NASA ADS] [Google Scholar]
 Rasia, E., Borgani, S., Murante, G., et al. 2015, ApJ, 813, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Sabra, B. M., Saliba, C., Abi Akl, M., & Chahine, G. 2015, ApJ, 803, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Sani, E., Marconi, A., Hunt, L. K., & Risaliti, G. 2011, MNRAS, 413, 1479 [NASA ADS] [CrossRef] [Google Scholar]
 Savorgnan, G. A. D., Graham, A. W., Marconi, A., & Sani, E. 2016, ApJ, 817, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Sdowski, A., & Gaspari, M. 2017, MNRAS, 468, 1398 [NASA ADS] [CrossRef] [Google Scholar]
 Shankar, F., Bernardi, M., Sheth, R. K., et al. 2016, MNRAS, 460, 3119 [NASA ADS] [CrossRef] [Google Scholar]
 Shankar, F., Bernardi, M., Richardson, K., et al. 2019, MNRAS, 485, 1278 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Steinborn, L. K., Dolag, K., Hirschmann, M., Prieto, M. A., & Remus, R.S. 2015, MNRAS, 448, 1504 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050 [NASA ADS] [CrossRef] [Google Scholar]
 Tozzi, P., & Norman, C. 2001, ApJ, 546, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, ApJ, 574, 740 [NASA ADS] [CrossRef] [Google Scholar]
 Truong, N., Rasia, E., Mazzotta, P., et al. 2018, MNRAS, 474, 4089 [NASA ADS] [CrossRef] [Google Scholar]
 van den Bosch, F. C., & Ogiya, G. 2018, MNRAS, 475, 4066 [NASA ADS] [CrossRef] [Google Scholar]
 Volonteri, M., & Ciotti, L. 2013, ApJ, 768, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Volonteri, M., Natarajan, P., & Gültekin, K. 2011, ApJ, 737, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Vulcani, B., De Lucia, G., Poggianti, B. M., et al. 2014, ApJ, 788, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberger, R., Springel, V., Pakmor, R., et al. 2018, MNRAS, 479, 4056 [NASA ADS] [CrossRef] [Google Scholar]
 Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009, MNRAS, 399, 574 [NASA ADS] [CrossRef] [Google Scholar]
 Wyithe, J. S. B. 2006, MNRAS, 365, 1082 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, H.Y. K., Gaspari, M., & Marlow, C. 2019, ApJ, 871, 6 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Effects of feedback parameters on SMBH growth
Cosmological hydrodynamical simulations rely on a number of subgrid models employed to simulate physical processes not resolved by the simulation itself (e.g., star formation and stellar feedback, chemical evolution, and AGN feedback). These models depend upon a set of parameters, tuned in order to reproduce some observational constraints. In particular in our simulations the parameters regulating the AGN feedback are chosen in order to reproduce the observed Magorrian relation (Magorrian et al. 1998). However, Shankar et al. (2016, 2019) have suggested that the observed correlation between stellar mass and SMBH mass might be biased high, because the SMBHs sphere of influence needs to be resolved in order to measure SMBHs mass with spatially resolved kinematics.
We, therefore, verified whether our results are influenced by a change of the normalization on the Magorrian relation and in which direction. In our subgrid model the normalization is regulated only by the radiative efficiency value, ϵ_{f}, that is the fraction of the energy radiated by AGN feedback which is thermally coupled to the surrounding gas particles. To perform this test, we used two simulations of the same Lagrangian region centered around a cluster with M_{500} = 3.38 × 10^{14} M_{⊙} at a mass resolution 10 times higher than the one used for the simulations presented in the paper and therefore also more spatially resolved. The two simulations differ only for the feedback efficiency ϵ_{f} of the quasarmode set equal to 0.1 in one simulation and to 0.15 in the other, while in the radio mode the parameter is kept equal to 0.7.
In Fig. A.1 we show the effect of the feedback efficiency on the M_{BH} − M_{⋆} relation. As we commented before the only effect is a change in the normalization, while the slope and the scatter are similar between the two runs. This result has been already studied both from a numerical (e.g., Booth & Schaye 2009) and theoretical (e.g., Churazov et al. 2005) point of view, therefore is not totally unexpected.
Fig. A.1.
Correlation between SMBH mass and stellar mass for the two different values of the radiative efficiency ϵ_{f}. Blue and red solid lines represent best fit to the data. 
As second check, we studied the effect on the M_{BH} − M_{500} relation presented in Sect. 4.1. From Fig. A.2 we see that also in this case a higher feedback efficiency leads to less massive black holes, and thus to a lower normalization of the M_{BH} − M_{500} relation. There is no effect on the total mass of the clusters. Indeed, M_{500} varies by little but not in a systematic direction: it might slightly increase as decrease.
Fig. A.2.
Correlation between SMBH mass and cluster mass. Black points represents simulated data of Fig. 4. Red and blue points represent results of the two simulations of the same region at higher resolution. Dashed lines connect the same systems in the two simulations. 
Finally, we study what is the effect of changing ϵ_{f} in the relative importance of the SMBH growth channels presented in Sect. 4.4: and . In Fig. A.3 we show the equivalent of Fig. 8 for the SMBH at the center of the main cluster in the two simulations. We see that a lower value of ϵ_{f} (red lines) leads to a fraction of mass due to gas accretion which is higher by about 20%. We repeated the same analysis on six SMBHs that we could identify in the Lagrangian region and match as part of corresponding clusters. We find that the median of the ratios between the fraction of SMBH mass due to gas accretion in the two simulations is 0.83, confirming that the plot is representative of the difference in the BH growth.
Fig. A.3.
Evolution of the mass of the same SMBH in the two simulations at high resolution. The solid lines represent the growth of the mass of the SMBHs due to gas accretion. Dashed lines represent the mass gained via BHBH mergers. 
In conclusion, if we chose a larger value of ϵ_{f} to match a lower normalization of the M_{BH} − M_{⋆} relation as suggested by Shankar et al. (2016), we expect that the M_{BH} will decrease but both the stellar mass and the cluster mass remains unchanged. As a consequence the scatter and slope of the M_{BH} − M_{⋆} and the M_{BH} − M_{500} relations will maintain their values. Finally, the limited growth of the SMBH is mostly caused by a much reduced level of gas accretion which is not fully replenished by the BHBH mergers. In fact, even thought the merger events play a more significant role with respect to the accretion, all the merged SMBHs are on their own smaller for the less effective accretion experienced at high redshift.
Appendix B: Galaxy stellar mass function (GSMF)
We compare the high mass end of the GSMF in our simulations with respect to the results presented in Bernardi et al. (2013). In our simulations the GSMF is computed as follows: we considered the region of space within R_{200} of all the 135 groups and clusters used in our paper. In this region we considered all the substructures but the BGGs and BCGs identified by Subfind and computed the stellar mass considering all the star particles which Subfind assign as bound to the substructure. The number of galaxies in each bin is than normalized following the prescription given in the appendix of Vulcani et al. (2014; Eq. (B2)). This factor takes into account the fact that we are in an overdense region with respect to the field.
We note that simulations predict a GSMF lower by a factor of ∼2 when comparing with observations at masses lower than ∼6 × 10^{11} M_{⊙}. This effect has been already discussed when studying the mismatch between the number of substructures in cluster environment as measured in observations and simulations (e.g., Grillo et al. 2015; Munari et al. 2016).
Fig. B.1.
High mass end of Galaxy stellar mass function in simulations (blue points) compared to GSMF of Bernardi et al. (2013). 
All Tables
Number of clusters, range of mass or temperature covered and their mean values for the mass sample (in the first half) and temperature subsample (in the second half).
All Figures
Fig. 1.
Correlation between stellar mass and SMBH mass in observations and simulations. Small lightblue points represent nonBCG simulated galaxies, large black dots represent simulated BCGs. Yellow, green, red, brown, and orange crosses represent the observational data with their error bars taken from McConnell & Ma (2013), Main et al. (2017), Savorgnan et al. (2016), Bogdán et al. (2018), and Gaspari et al. (2019) respectively. The dashed black line is a linear bestfit of the sample of different type of galaxies by McConnell & Ma (2013). See text for details about M_{BH} and M_{⋆} definition and measurement. 

In the text 
Fig. 2.
Correlation between SMBH mass, M_{BH}, and the clusters temperature, T_{500}. Red and orange crosses refer to observational data from Bogdán et al. (2018) and Gaspari et al. (2019) respectively. Darkblue dots represent simulated clusters in the temperature subsample while cyan stars show the remaining objects of the mass sample. Dashed red line is the prediction of the toy model by Churazov et al. (2005). 

In the text 
Fig. 3.
Correlation between cluster mass, M_{500}, and cluster temperature, T_{500}. Symbols as in Fig. 2, where the observational data are taken from Lovisari et al. (2015). Dashed lines are the bestfitting lines for both simulated and observed data. 

In the text 
Fig. 4.
Correlation between SMBH mass, M_{BH}, and cluster mass, M_{500}. Symbols as in Fig. 2. 

In the text 
Fig. 5.
Evolutionary tracks of three different systems on the M_{BH}–M_{500} plane. Triangles over each line indicate the position of the systems on the plane at z = 3, z = 2, z = 1 and z = 0. Black circles represent our numerical sample at z = 0; the filled ones are systems for which M_{500} is increased by more than 40% in the last Gyr. 

In the text 
Fig. 6.
Ratios between cluster (SMBH) masses at z = 2 and cluster (SMBH) masses at z = 0 as a function of the cluster masses at z = 0. Clusters are shown as green squares and SMBHs as yellow triangles). The sample used is the massselected sample identified at z = 2. 

In the text 
Fig. 7.
Correlation between M_{BH} and M_{500} at different redshifts. At every redshift we show the bestfitting relation in the mass range of the respective samples. Namely, we show in green, orange, red, magenta, and blue the mass sample related to z = 2, 1.5, 1, 0.5 and z = 0, respectively. 

In the text 
Fig. 8.
Evolution of the mass of black holes showed in Fig. 5. The solid lines represent the mass of the black holes due to gas accretion. Dashed lines represent mass gained via BHBH mergers. 

In the text 
Fig. 9.
Evolution of SMBH mass divided into two growth channels (gas accretion in red and BHBH mergers in blue) considering the complete sample. In black we also plotted the total SMBH mass. Solid lines represent median values of our sample and shadowed regions represent 16 and 84 percentiles. 

In the text 
Fig. 10.
Same as Fig. 9 but considering only clusters with M_{500} > 10^{15} M_{⊙} at z = 0. 

In the text 
Fig. 11.
Correlation between SMBH mass and M_{500} at z = 1.0 (left panel) and z = 0 (right panel). We plotted in blue the mass gained by mergers and in red the mass gained by gas accretion. Dashed lines represent best fitting lines to the data. 

In the text 
Fig. 12.
Evolution of M_{⋆,in}, defined as the stellar mass inside a spherical region of radius 30 kpc, SMBH mass, and M_{⋆,out}, defined as the total mass of the stars enclosed in a spherical shell with radii 30 kpc and 100 kpc. All the quantities are normalized to their respective values at z = 0. Solid lines represent median values and shadowed regions represent 16 and 84 percentiles. The three different blue lines represent three definition of M_{⋆,out}. In particular the thin blue solid line is the stellar mass in a spherical shell with radii 100 kpc and 200 kpc while the dashed blue line is the stellar mass in a spherical shell with radii 50 kpc and 350 kpc. 

In the text 
Fig. 13.
Covariance matrix. Each axis represents the logarithmic difference between the actual value of a quantity X and the expected value from the linear relation (XM_{500}) at its fixed mass, as defined by Eq. (5). Panels above the diagonal refer to z = 1 while panels below the diagonal refer to z = 0. The diagonal panels show the distribution of δ(X) at z = 0. Red points define clusters which have experienced a mass growth of at least 40% during the last Gyr. 

In the text 
Fig. A.1.
Correlation between SMBH mass and stellar mass for the two different values of the radiative efficiency ϵ_{f}. Blue and red solid lines represent best fit to the data. 

In the text 
Fig. A.2.
Correlation between SMBH mass and cluster mass. Black points represents simulated data of Fig. 4. Red and blue points represent results of the two simulations of the same region at higher resolution. Dashed lines connect the same systems in the two simulations. 

In the text 
Fig. A.3.
Evolution of the mass of the same SMBH in the two simulations at high resolution. The solid lines represent the growth of the mass of the SMBHs due to gas accretion. Dashed lines represent the mass gained via BHBH mergers. 

In the text 
Fig. B.1.
High mass end of Galaxy stellar mass function in simulations (blue points) compared to GSMF of Bernardi et al. (2013). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.