© 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 co-evolution 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 active-galactic-nuclei (AGN) feedback. In this scenario gas settles around the SMBH radiating energy at a rate of ∼ηṀc², 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 non-causal 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., Ragone-Figueroa 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 self-regulated 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₅₀₀ 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 X-ray 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₅₀₀). Finally, Lakhchaura et al. (2019) considered a sample of 41 ETGs and confirmed the tight correlation between the SMBH mass and the X-ray 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 thermo-dynamical 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₅₀₀ (defined as the radius of the sphere of mass M₅₀₀ whose enclosed density is 500 times the critical density of the Universe at a specific redshift, 3H(z)²/8πG) by employing a set of cosmological hydro-dynamical simulations centered on massive clusters. These zoom-in simulations include a number of sub-grid 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 (Ragone-Figueroa 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₅₀₀ − M_BH and M₅₀₀–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 BH-BH mergers) do SMBHs grow in time? (5) Is M₅₀₀ 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 hydro-dynamical zoom-in 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, σ₈ = 0.8, and H₀ = 100 h km s⁻¹ Mpc⁻¹ = 72 km s⁻¹ Mpc⁻¹. These regions were selected from a parent N-body cosmological volume of 1 h⁻³ Gpc³ and re-simulated at higher resolution with the inclusion of baryons (for a detailed description of the initial conditions see Bonafede et al. 2011).

The re-simulated regions are centered around the 24 most massive clusters of the parental box with mass M₂₀₀ ≥ 8 × 10¹⁴ h⁻¹ M_⊙ and 5 isolated groups with M₂₀₀ within [1 − 4]×10¹⁴ h⁻¹ M_⊙. In the high-resolution regions the mass of DM particles is m_DM = 8.47 × 10⁸h⁻¹ M_⊙ and the initial mass of the gas particle is m_gas = 1.53 × 10⁸h⁻¹ M_⊙. The Plummer equivalent gravitational softening for DM particles is set to ϵ = 5.6 h⁻¹ 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⁻¹ kpc, 3 h⁻¹ kpc, and 3 h⁻¹ kpc, respectively.

The simulations were carried out with the code GADGET-3, a modified version of the Tree-PM Smoothed-Particle-Hydrodynamics (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 hydro-dynamical 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 zoom-in simulations treats unresolved baryonic physics through various sub-grid 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 metal-dependent 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 Ragone-Figueroa et al. (2013) with one important modification (see, Ragone-Figueroa 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:

(1)

with α equal to 10 and 100 for hot (T >  5 × 10⁵ K) and cold (T <  5 × 10⁵ 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 re-positioned 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⁻¹ 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, Ragone-Figueroa 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 light-blue points represent non-BCG 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 best-fit of the sample of different type of galaxies by McConnell & Ma (2013). See text for details about MBH and M⋆ definition and measurement.

SMBHs of mass 4.4 × 10⁵ M_⊙ are seeded at the position of the most bound particle of the structures that, identified by the Friends-of-Friends algorithm, simultaneously satisfy all the following conditions: the stellar mass of the structure is greater than 2.2 × 10¹⁰ M_⊙ and it is higher than 5% of the dark-matter 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 BH-BH 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 set-up 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 light-blue points represent non-central simulated galaxies identified by Subfind (Dolag et al. 2009), while dark-blue dots represent simulated BCGs. The stellar masses of the BCGs are defined as the mass enclosed in a sphere of radius 0.1 × R₅₀₀ around the position of the central SMBH, while total stellar masses of non-BCGs are given as an output by Subfind.

To calibrate the parameters for the AGN-feedback model we aimed at reproducing the entire M_BH − M_⋆ relation including the majority of simulated non-BCG 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 K-band 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 mass-to-light 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 K-band 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 non-BCG 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¹² 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 non-BCG 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 mass-to-light 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 Ragone-Figueroa 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₅₀₀ 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₅₀₀ ≥ 1.4 × 10¹³ M_⊙ or, equivalently, with at least 20 thousands particles within R₅₀₀. The properties of the mass-selected sample are summarized in top part of Table 1.

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₅₀₀ 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 BH-BH 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 XMM-Newton observations we considered the spectroscopic-like temperature (Mazzotta et al. 2004):

(2)

where ρ_i, m_i, and T_i are the density, mass, and temperature of the ith gas particle within R₅₀₀ emitting in the X-ray 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₅₀₀ we imposed two conditions: a minimum of 10⁴ hot gas particles and a maximum fraction of 5% of gas particles discarded because too cold. All clusters satisfying these requirements have also M₅₀₀ > 1.4 × 10¹³ M_⊙, thus whenever we consider measurements of temperature we refer to a subsample of the mass selected-sample. The properties of the temperature-selected 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.

Best-fitting procedure. For all the considered relations, we looked for a best-fit line in the logarithmic plane:

(3)

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₀ or Y₀ and respectively equal to 2 keV, 10¹⁴ M_⊙, 10¹¹ M_⊙, and 10⁹ M_⊙.

To find the best-fit curve, we employed an IDL routine that is resistant with respect to outliers: ROBUST_LINEFIT¹. 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 best-fitting 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 best-fitting 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)². Bogdán et al. (2018) derive the temperature from XMM-Newton observations of the hot gas; Gaspari et al. (2019) use published data of Chandra and wide-field ROSAT/XMM Newton. The extraction region for both is on average ∼0.2R₅₀₀ (core included; such temperature is a good proxy for T₅₀₀).

In Fig. 2 we show the correlation between M_BH and T₅₀₀ 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)³, 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 under-representation of clusters with T₅₀₀ >  2 keV in the observational sample that, however, has a good number of systems below 1 keV.

Fig. 2. Correlation between SMBH mass, MBH, and the clusters temperature, T500. Red and orange crosses refer to observational data from Bogdán et al. (2018) and Gaspari et al. (2019) respectively. Dark-blue 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₅₀₀ by using the mass sample rather than the less-numerous temperature sub-sample (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):

(4)

For this reason, before analyzing the M_BH–M₅₀₀ relation, it is helpful to compare the observed and simulated T₅₀₀–M₅₀₀ relations. The observed and simulated data sets are shown in Fig. 3 together with the best-fitting 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, M500, and cluster temperature, T500. Symbols as in Fig. 2, where the observational data are taken from Lovisari et al. (2015). Dashed lines are the best-fitting 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 best-fit line suggests a possible difference in the hottest-clusters 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₅₀₀ = 10¹⁴ 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 X-ray hydrostatic masses are biased low.

Finally, in Fig. 4 we compare the correlation between M_BH and the M₅₀₀ 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₅₀₀ > 3 × 10¹⁴ M_⊙) and low masses (M₅₀₀ < 3 × 10¹³ M_⊙). In the intermediate mass range, the few observed data points tend to have slightly higher SMBH masses than the simulated objects. This apparent mis-match is presumably a consequence of the poor statistics of 10¹⁴ M_⊙ objects in the observational sample. More unlikely, this feature could suggest a broken power law to describe the M_BH − M₅₀₀ relation, but such a drastic change in the slopes is difficult to justify.

Fig. 4. Correlation between SMBH mass, MBH, and cluster mass, M500. Symbols as in Fig. 2.

4.2. The theoretical M_BH–M₅₀₀ 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₅₀₀ relation shown in Fig. 5. For this goal, we over-plotted three evolutionary tracks of representative SMBHs. These objects are chosen accordingly to their mass at z = 0; specifically, they refer to a small, a medium-mass, 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 MBH–M500 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 M500 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₅₀₀. 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⁸ − 10⁹ 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₅₀₀ relation. This happens before z = 2 and in some cases even at z >  4.

After this initial phase, the cluster and its central SMBH co-evolve, 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₅₀₀. The former situation occurs during cluster mergers. It starts when merging structures reach and cross R₅₀₀ 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₅₀₀ 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 BH-BH merger or both. The event is captured by the vertical movement in the plot.

All these time-frames 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₅₀₀ 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₅₀₀–M_BH relation.

4.3. Evolution of the M_BH–M₅₀₀ relation

After the inspection on the trajectory of individual objects we study here the evolution of the entire M₅₀₀–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₅₀₀, 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 best-fitting Δ_M500 − M₅₀₀ relation is, indeed, greater than 0.75. Clusters with a final mass lower than 10¹⁴ M_⊙ increase their total masses by a factor between 3 and 6. Instead, massive clusters with final M₅₀₀ > 10¹⁵ 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 mass-selected 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 volume-limited sample because the most massive objects, M₅₀₀ >  10¹⁵ M_⊙ correspond to the most massive systems of the parent volume-limited cosmological box (see Sect. 2). Vice versa, the scatter for the smallest systems is likely under-estimated. 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₅₀₀ relation between z = 2 and z = 0.

This is confirmed in Fig. 7 where we plot the best-fitting 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 MBH and M500 at different redshifts. At every redshift we show the best-fitting 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 BH-BH 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 BH-BH 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 BH-BH 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 BH-BH 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, BH-BH 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₅₀₀ below 10¹⁴ M_⊙, the most massive above 10¹⁵ 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.

Fig. 10. Same as Fig. 9 but considering only clusters with M500 >  1015 M⊙ at z = 0.

The stronger relative influence by the SMBH accretion with respect to the BH-BH 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 BH-BH 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₅₀₀. We find that the ratio of the two quantities is a strongly decreasing function of M₅₀₀ 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 BH-BH merger, we computed at z = 0 the number and mass of SMBHs which are inside 0.1 × R₅₀₀ 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₅₀₀ <  10¹⁴ 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 well-established 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. Ragone-Figueroa 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 “not-bound” 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 Ragone-Figueroa 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 BH-BH 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₅₀₀ 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 BH-BH 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₅₀₀ 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 M500 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₅₀₀ 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₅₀₀ 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 σ_MBH|MBCG = 0.61 and σ_MBH|T = 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₅₀₀ 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₅₀₀. 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 best-fitting 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₅₀₀ relations⁴ provided in Table 2:

(5)

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 (X-M500) 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₅₀₀, 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₅₀₀) 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.

As previously commented, the scatter σ_MBH|T500 is comparable to σ_MBH|M500 (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₅₀₀) 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 time-delays explain why the correlation coefficients between δ(T₅₀₀) 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 BH-BH 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 σ_MBH|MBCG is smaller at all times. At redshift z = 0, for example, σ_MBH|M500 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 σ_MBH|M500 is ≈1.4 × σ_MBH|MBCG. 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 pre-mergers 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₅₀₀ 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⁵, we find that the scatter of the M_BH − M₅₀₀ relation reduces by almost 20% (σ_MBH|M500 = 0.14 ± 0.01), it remains comparable to the re-computed scatter of the M_BH − T₅₀₀ relation (σ_MBH|T500 = 0.14 ± 0.02) and consistent with the new σ_MBH|MBCG = 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 σ_MBH|MBCG 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 Ragone-Figueroa 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 intra-cluster 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 multi-temperature 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 X-ray. 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₅₀₀, and global temperature, T₅₀₀. Our work is based on 29 zoom-in cosmological hydrodynamical simulations carried out with GADGET-3, a modified version of the public code GADGET-2. This code treats unresolved baryonic physics including AGN feedback through various sub-grid 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₅₀₀ >  1.4 × 10¹³ M_⊙ identified in the high-resolution regions of the re-simulated 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 co-evolution 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 BH-BH merging. Our main results can be summarized as follows.

• The simulated relations between SMBH mass and the global cluster properties (M_BH–M₅₀₀ and M_BH − T₅₀₀) 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₅₀₀ relation at z = 0 originates from a non-simultaneous growth of the SMBH and of the cluster. In particular, objects evolve on the M_BH − M₅₀₀ plane either at almost constant M_BH or at almost constant M₅₀₀ (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 BH-BH merger. Clusters that recently experienced major merger events are in general below the mean relation.

• The M_BH − M₅₀₀ 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 BH-BH 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 BH-BH merger rate. and similarly relate to M₅₀₀ 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 Ragone-Figueroa 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 intra-cluster light), that we simply define here as unbound stars located outside the BCG.

• The M_BH − M₅₀₀ and M_BH − T₅₀₀ relations present a similar scatter (or, equivalently, δ(T₅₀₀) 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₅₀₀ and M_BH − T₅₀₀, become similar, thus strengthening the predicting power of the cluster global quantities.

Even though our simplified sub-grid 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 so-far limited to about 5 kpc to obtain a large statistical sample in a cosmological environment. This means that while the large-scale 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¹¹ 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 GR-MHD BH accretion simulations (e.g., Sdowski & Gaspari 2017).

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 GADGET-3 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 ASI-INAF n.2017-14-H.0, the Agencia Nacional de Promoción Científica y Tecnológica de la República Argentina under the PICT2417-2013 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 Sklodowska-Curie grant agreement No 734374, PRIN-MIUR 2015W7KAWC, the INFN INDARK grant, NASA Chandra GO7-18121X and GO8-19104X. Simulations have been carried out using MENDIETA Cluster from CCAD-UNC, which is part of SNCAD-MinCyT (Argentina) and MARCONI at CINECA (Italy), with CPU time assigned through grants ISCRA C, and through INAF-CINECA and University of Trieste – CINECA agreements. The post-processing has been performed using the PICO HPC cluster at CINECA through our expression of interest.

References

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 sub-grid 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₅₀₀ = 3.38 × 10¹⁴ 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 quasar-mode 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₅₀₀ 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₅₀₀ relation. There is no effect on the total mass of the clusters. Indeed, M₅₀₀ 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 BH-BH 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₅₀₀ 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 BH-BH 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₂₀₀ 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¹¹ 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

All Figures

Fig. 1. Correlation between stellar mass and SMBH mass in observations and simulations. Small light-blue points represent non-BCG 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 best-fit of the sample of different type of galaxies by McConnell & Ma (2013). See text for details about MBH and M⋆ definition and measurement.

In the text

Fig. 2. Correlation between SMBH mass, MBH, and the clusters temperature, T500. Red and orange crosses refer to observational data from Bogdán et al. (2018) and Gaspari et al. (2019) respectively. Dark-blue 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, M500, and cluster temperature, T500. Symbols as in Fig. 2, where the observational data are taken from Lovisari et al. (2015). Dashed lines are the best-fitting lines for both simulated and observed data.
In the text

	Fig. 4. Correlation between SMBH mass, MBH, and cluster mass, M500. Symbols as in Fig. 2.
In the text

	Fig. 5. Evolutionary tracks of three different systems on the MBH–M500 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 M500 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 mass-selected sample identified at z = 2.
In the text

	Fig. 7. Correlation between MBH and M500 at different redshifts. At every redshift we show the best-fitting 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 BH-BH mergers.
In the text

	Fig. 9. Evolution of SMBH mass divided into two growth channels (gas accretion in red and BH-BH 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 M500 >  1015 M⊙ at z = 0.
In the text

	Fig. 11. Correlation between SMBH mass and M500 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 (X-M500) 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 BH-BH 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