Black hole mass and cluster mass correlation in cosmological hydro-dynamical simulations
Key Words.:galaxy: clusters: intracluster medium – galaxy: nuclei – method: numerical – X-ray: galaxies: cluster
Context:The correlations between the properties of the brightest cluster galaxy (BCG) and the mass of its central super-massive 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 relation. Moreover, we study how gas accretion and BH-BH mergers contribute to SMBH growth across cosmic time.
Methods:We employ 135 groups and clusters with a mass range extracted from a set of 29 zoom-in cosmological hydro-dynamical simulations where the baryonic physics is treated with various sub-grid models, including feedback by active galactic nuclei.
Results:In our simulations we find that well correlates with and , with the scatter around these relations compatible within with the scatter around at . The relation evolves with time, becoming shallower at lower redshift as a direct consequence of hierarchical structure formation. In our simulations, SMBHs mainly grow by gas accretion at redshift . At redshift the main growth channel is instead the BH-BH merging. During this last process, substructures hosting BHs are disrupted in the merger process with the BCG and the unbound stars enrich the diffuse stellar component rather than contribute to increase BCG mass.
Conclusions:From the results obtained in our simulations with simple sub-grid models we conclude that the scatter around the relation is comparable to the scatter around the 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.
It is well known that galaxies of every morphology host super massive black holes (SMBHs) at their center. Interestingly the mass of these black holes111Here and in the following sections with BHs we always refer to super-massive black holes mostly hosted at the center of BCG. (BHs) well correlate with a number of bulge properties of the host galaxy such as bulge stellar mass (eg., MAGORRIAN1998, 2003MARCONI, 2004HARING, 2009HU, 2011SANI and 2013KORMENDY for a review), bulge stellar velocity dispersion (e.g., 2000FERRARESE, 2000GEBHARDT, 2001MERRITT, 2002TREMAINE, 2006WYITHE, 2008HU, 2009GULTEKIN, 2012BEIFIORI, 2013mcconnel) and bulge luminosity (e.g., 1995KORMENDY, 2002MCLURE, 2009HU, 2011SANI). These correlations are as tight as to be often used to estimate the BHs mass when dynamical measurements are not available and suggest a co-evolution between BH 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-galaxy-nuclei (AGN) feedback. In this scenario gas settle around the BH radiating energy at a rate of , with . 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 BH (1999FABIAN, 2004GRANATO, 2005DIMATTEO, 2006HOPKINS). However, other authors have shown that the M-M relation can arise or be contributed by non-causal processes. In this scenario the observed M-M relation follows from hierarchical galaxy mergers starting from an uncorrelated distribution of M and M (eg., 2007PENG, 2011KUND). Moreover, it has also been suggested that, at least for massive elliptical galaxies, the main correlation is not with the bulge properties but with the host dark matter halo (2002FERRARESE, 2009BANDARA, 2010BOOTH, 2011MARTA). However, a handful of observational studies show that there is no strong correlation between BH mass and circular rotation velocities of gas in the outer parts of galaxies without bulge, which the authors use as a proxy of halo mass, leading to the conclusion that BH do not correlate directly with dark matter (2011NATURE, 2015SABRA).
Between all galaxies a special position is occupied by the brightest cluster galaxies (BCGs). These galaxies, set at the center of galaxy clusters, host the most massive BHs, that seem over massive at fixed galaxy stellar velocity dispersion and stellar mass (2000GEBHARDT, 2012LARRONDO, 2015FERREMATEU) with respect to other no-BCGs galaxies. GS17 initially proposed that BHs in BCGs correlate with the core properties () of hot halos (such as the plasma X-ray temperature), where the feeding via chaotic cold accretion (CCA) boosts the AGN feedback in a tight self-regulated feedback loop that prevent catastrophic cooling flows for several gigayears. Albeit in a small sample, 2018BOGDAN went further linking the mass of BHs to the large scale properties of the hosting cluster, such as and , which are tightly related to the gravitational potential of the systems. The authors remarkably found a scatter in lower than that of the relation. This suggests that given their preferential location, at the bottom of the potential well of massive structures, the growth of BHs is influenced by the physical processes that are also regulating the thermo-dynamical properties of the intra-cluster-medium (ICM).
In this work we employ a set of cosmological hydro-dynamical simulations centred on massive clusters to investigate the correlation between the mass of BHs in BCGs and the global properties of the hosting cluster, such as temperature and mass. 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 BH accretion and consequent AGN feedback (2013RAGONE). 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 (2015RASIA, 2017PLANELLES, 2017BIFFI). Given these previous results, we investigate further these simulated regions to study the physical processes that tie BH to the hosting cluster. In particular in this work we aim at answering the following questions: 1) do numerical simulations reproduce the observed and - relations? 2) Which are the processes that lead to the observed relations? 3) Do the relations evolve with redshift? 4) Through which channels (i.e. gas accretion vs merger) do SMBHs grow in time? 5) Is as appropriate as to probe ?
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.
Our analysis is based on a set of 29 cosmological and hydro-dynamical zoom-in simulations centered on massive galaxy clusters evolved in a CDM model with parameters: , , , and km s Mpc km s Mpc. These regions are selected from a parent -Body cosmological volume of 1 Gpc and re-simulated at higher resolution with the inclusion of baryons (for a detailed description of the initial conditions see Bonafede).
The re-simulated regions are centered around the 24 most massive clusters of the parental box with mass222 indicates the mass measured within the radius at which the enclosed density is times the critical density of the Universe at the specific redshift, /. M and 5 isolated groups with within M. In the high-resolution regions the mass of DM particles is M and the initial mass of the gas particle is M. The Plummer equivalent gravitational softening for DM particles is set to in comoving units at redshift higher than and in physical units afterward. The gravitational softening lengths of gas, stars and black hole particles are fixed in comoving coordinates to kpc, kpc, and kpc, respectively.
The simulations are carried out with the code GADGET-3, a modified version of the Tree-PM Smoothed-Particle-Hydrodynamics (SPH) public code GADGET2 (GADGET2005). Our simulations are performed with an improved version that accounts for modifications of the hydrodynamic scheme to better capture hydro-dynamical instabilities (see BECK2016). 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 2017PLANELLES or 2017BIFFI; we briefly summarize here the main aspects. The prescription of metal-dependent radiative cooling follows WIERSMA2009. The model of star formation and associated feedback prescriptions are implemented according to the original model by SPRINGEL2003 and metal enrichment and chemical evolution following the formulation by LUCA2007. The yields used in our simulations are specified in 2018BIFFI. The AGN feedback model is implemented as described in Appendix A of 2013RAGONE with one important modification (see, 2018CINTHIA): the distinction between cold mode and hot mode gas accretion (see also 2015RASIA). 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 ( K) and cold ( K) gas respectively. is the BH 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, is the sound speed of the gas surrounding the BH, and is the relative velocity between the BH and bulk velocity of the gas. 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 BH softening length. This calculation is restricted to particles with relative velocity with respect to the BH below s km. This condition avoids that the BH particle ”jumps” into a close flyby structure that would displace it from the cluster center.
A fraction of the energy associated to the gas directly fueling the BH through accretion is radiated away and a fraction of this energy is thermally coupled to the surrounding gas particles. The value of is fixed to , while that of depends on the mode of the AGN: during the quasar mode, i.e., for , , while during the radio mode is increased to (see, 2018CINTHIA). The exact values of both parameters are chosen to reproduce the observed correlation between stellar mass and BH mass in galaxies (see Fig. 1).
BHs of mass 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 and it is higher than 5 percent of the dark-matter halo mass; the ratio between the gas mass and the stellar mass is higher than , no other central BHs are already present. These conditions allow the seeding of the BH in galaxies that have enough gas to promptly feed the BH.
Two BH particles merge whenever their relative velocity is smaller than and their distance is less than twice the BH softening length. When a BH-BH merger happens, the BH particle of the most massive BH gains the mass of the merged BH.
The strategy to position the BH and the recipe to implement the AGN feedback are the only difference with respect to the simulation set-up of the runs presented in 2015RASIA, 2017BIFFI, 2017PLANELLES, 2018BIFFI, 2018TRUONG.
2.1 Calibration of AGN feedback model
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, , and the mass of their central BHs, . In the figure the small light-blue points represent non-central simulated galaxies identified by Subfind (Dolag2009), while dark-blue dots represent simulated BCGs. The stellar masses of the BCGs are defined as the mass enclosed in a sphere of radius around the position of the central BH, while total stellar masses of non-BCGs are given as an output by Subfind.
To calibrate the parameters for the AGN-feedback model we aim at reproducing the entire relation including the majority of simulated non-BCG galaxies. These are, in particular, compared with the observational data from 2013mcconnel represented in the figure by the dashed line.
In the plot we also include other observational data, namely the BCGs from 2013mcconnel and the samples from Savorgnan2016, 2017MAIN and 2018BOGDAN. In 2017MAIN BH masses are computed from -band luminosities using the relation suggested by 2007GRAHAM and extracted from a sample of elliptical but not BCGs. In all the other works the mass of the BHs are derived from dynamical measurements. The BCG masses in 2013mcconnel and 2018BOGDAN 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 Savorgnan2016 the stellar masses are computed from bulge luminosities assuming a constant mass-to-light ratio. Note that from Savorgnan2016 we used only ellipticals, that are not necessarily BCGs. In 2017MAIN the stellar masses are computed from K-band luminosity using the relation log given by 2003KBAND.
Our main condition for deciding on the AGN parameters is that the observed correlation between BH 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 (2013mcconnel) is high in the high mass end. Regarding the BCGs, we found 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 , 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 BH masses are computed by adopting different methods. For example, those extracted by 2017MAIN 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 and the adopted stellar mass-to-light ratio. The resulting differences between different 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 BH. This object is present in the Savorgnan2016, 2013mcconnel and 2018BOGDAN samples and, while M is identical because taken from the same source in the literature, the estimated BCG mass can be different even by a factor of . 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 BH mass and the two contributions to its growth: the accretion of the surrounding gas and the merger with other BHs. 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 BHs are always positioned at the location of the most bound particle that should identify the center of the host halo. Therefore, we follow the BHs back in time to identify the position of the hosting halo center. For this goal, we save at z=0 the unique
identification number of the most massive BH particle which is within 10 kpc from the cluster minimum of the potential well, as identified by Subfind. Subsequently, we track it back in time to the epoch of its seeding. At each time, we check that the BH 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 build a merger tree of the central BHs 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 2018CINTHIA). However, we verify that for 80% of our systems the main progenitor of the BH is at the center of the main progenitor of the cluster up to and for half of these the two trees coincide till the time when the BH is seeded.
Cluster masses. Once the center is defined as above, we consider the total gravitational mass of the cluster within an aperture radius 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 or, equivalently, with at least 20 thousands particles within . The properties of the mass-selected sample are summarized in top part of Table 1.
BCG stellar mass. We define the mass of the BCG as the stellar mass enclosed in a sphere of radius around the cluster center.
BH mass. Given the identification number of a BH particle, the mass of the BH, , at every redshift is quite easily retrieved from the simulation as it is the mass associated to that particle. The total mass of BH 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 BH to increase its mass. The accretion mass () is obtained by integrating the accretion rate, information that we save 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 BH mass by mergers is negligible at . 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 consider the spectroscopic-like temperature (2004Mazzotta):
where , , and are the density, mass, and temperature of the gas particle within emitting in the X-ray band, i.e. with keV and a cold fraction333According to the effective model by SPRINGEL2003 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. lower than 10 per cent. In order to have a reliable estimation of the temperature inside we impose two conditions: a minimum of hot gas particles and a maximum fraction of 5 per cent of gas particles discarded because too cold. All clusters satisfying these requirements have also , thus whenever we will consider measurements of temperature we will 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 because at the highest redshift bins, and , we do not have enough statistics to apply a meaningful analysis.
Best-fitting procedure. For all the considered relations, we look for a best-fit line in the logarithmic plane:
where always indicates the decimal logarithm. The temperature, cluster mass, BCG stellar mass and BH mass are always normalized by the same factors, expressed above as or and respectively equal to 2 keV, , , and .
To find the best-fit curve, we employ an IDL routine that is resistant with respect to outliers: ROBUST_LINEFIT444https://idlastro.gsfc.nasa.gov/ftp/pro/robust/robust_linefit.pro. For the simulated data, we always consider 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 generate ten thousand bootstrap samples by randomly replacing the data. From the resulting distributions we derive 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 LABEL:tab:all and will be discussed in the next two sections.
4.1 Comparison with observational data
In this section we compare the numerical results to the observations presented in 2018BOGDAN, where the correlation between the mass of the SMBHs in BCG and the global temperature of clusters and groups of galaxies is presented. The temperature is derived from XMM-Newton observations of the hot gas inside without any core exclusion.
In Fig. 2 we show the correlation between and for our simulations and for their observational data set. 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 keV in the observational sample that, however, has a good number of systems below 1 keV.
In order to better populate the colder and less massive tail of the simulated cluster distribution we compare the correlation between and by using the mass sample rather than the less-numerous temperature sub-sample (Table 1). In 2018BOGDAN the total mass was not measured directly from their data but was derived from the temperature via the scaling relation by 2015LOVISARI:
For this reason, before analyzing the - relation, it is helpful to compare the observed and simulated - 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 verify that our fitting procedure555For this fit we did not consider the BISECT option because in the observed sample the temperature is the true independent variable. returned the same parameters of Eq. 4. In particular we confirm the value of the slope reported in 2015LOVISARI: . The results of the best-fitting line of the simulated data are reported in Table LABEL:tab:all.
By looking at Fig. 3, we see a good agreement between simulated and observed clusters in the temperature range covered by 2015LOVISARI. 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 per cent higher at . This feature is not new and has been already noted in 2018TRUONG, where a similar set of numerical simulations is employed, and, more interestingly, in other numerical analysis, such as FABLE. In particular, in their work the authors showed that numerical results are in agreement with observations if are considered only cluster masses derived via weak lensing. This suggest that X-ray hydrostatic masses are biased low.
Finally, in Fig. 4 we compare the correlation between and the as measured in our simulations and as derived in 2018BOGDAN. The results of the comparison are as expected from the previous two figures: the simulated data points are in line with observations, especially at high () and low masses (). In the intermediate mass range, the few observed data points tend to have slightly higher BH masses than the simulated objects. This apparent mis-match is presumably a consequence of the poor statistics of objects in the observational sample. More unlikely, this feature could suggest a broken power law to describe the relation, but such a drastic change in the slopes is difficult to justify.
4.2 The theoretical - relation
Since the simulated sample is in overall good agreement with the observed correlation between the mass of the central BH and the mass of the clusters, we investigate here how single simulated systems evolve throughout time to form, by , the relation shown in Fig.5. For this goal, we over-plot three evolutionary tracks of representative SMBHs. These objects are chosen accordingly to their mass at ; specifically, they refer to a small, a medium-mass, and a massive BH. To have some temporal reference we also indicate four specific times along each line: , , and .
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 BHs grows rapidly at almost constant . This track begins instantaneously as the BHs are seeded in a gas rich region. The BHs immediately gain mass at the Eddington limit by the accretion from the abundant surrounding gas. This phase typically lasts half Gyr and can lead to the formation of BHs with a mass already of the order of . The fast BH growth ends when the 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 BHs in our sample are close to the relation. This happens before and in some cases even at .
After this initial phase, the cluster and its central SMBH co-evolve, but not with a linear evolutionary track: the increase of the BH 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 or at almost constant . The former situation occurs during cluster mergers. It starts when merging structures reach and cross leading to a quick increment of the total cluster mass and finishes when the secondary objects are fully incorporated. These horizontal shifts in the - 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 substructure moves towards the center of the cluster and neither the BH mass nor the cluster mass change. Eventually, the merging object reaches the core and either feeds the central BH with gas or induces 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 masses 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 BH 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 grows by more than 40 per cent in the last Gyr (shown as black points in Fig. 5) have BHs that on median are 50 per cent smaller than those expected to follow the M-M relation.
4.3 Evolution of the M-M relation
After the inspection on the trajectory of individual objects we study here the evolution of the entire M-M relation. We start with the evaluation of the ratio between the mass gained by clusters and by central BHs between and : , where refers to either the total mass or the BH mass . If these two ratios are constant, the slope of the relation will not change. The resulting ratios are shown in Fig. 6 as function of the cluster total mass reached at for the mass sample identified at .
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 relation is, indeed, greater than 0.75. Clusters with a final mass lower than increase their total masses by a factor between 3 and 6. Instead, massive clusters with final increase their mass on average by a factor of about 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.
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 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 BH mass is independent of the cluster mass and that BHs grow on average by a factor of about 5-6. As a consequence we expect a marked evolution of the slope of the - relation between and .
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 , , is almost twice the value found at , . From Fig. 6 and Fig. 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.
4.4 Evolution of BHs 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 BH mass: the accretion by the diffuse gas, , and the merger with other BHs, .
First, we analyze how the mass of single BHs evolves with time via the two separate channels. As an example in Fig. 8 we plot the evolution of the three BHs 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 and . After that, BHs still grow by gas accretion but at a much smaller rate. On the contrary the increase of the BH mass due to BH-BH mergers becomes more important and it is the main channel of the BH mass growth at lower redshifts, i.e. for the two most massive BHs and for the smallest one.
In Fig. 9 we show the evolution of our complete sample. We plot the median behaviors with a solid line and the 68 per cent of the total sample distributions (from the 16 to the 84 percentiles) with the shaded regions. The total BH mass and the masses gained from the two channels are normalized with respect to the total BH mass at . Finally, the vertical dashed lines help to identify three significant times: , 1, and 2.
From the BH seeding up to the total mass of the BH grows almost entirely by gas accretion. Half of the final mass gained through gas accretion is, indeed, accumulated before . Then from to the mass growth due to BH-BH merger becomes more relevant increasing at a rate comparable to the growth rate of . By makes up on average 25 per cent of the total mass at that redshift. At lower redshift, , BH-BH mergers constitute the main channel for BH mass growth, and eventually represents the main component of mass gained by , in line with some previous results (see 2013MARTA, 2014DUBOIS). Indeed, the mass accumulated by gas accretion from to accounts only for per cent of the total final mass, while during the same period the BH gains per cent 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, explore whether this is due to the broad mass range investigated and, thus, whether the described behaviour depends on the mass of the systems. We divided the sample in three mass bins: at the least massive objects have below , the most massive above , and the intermediate in between these two thresholds. We found that the trends of the relative ratio of the two BH 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 (see Fig. 10). As a result, at accounts for per cent of . 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 BH accretion with respect to the BH-BH merger can be due to two different factors: on the one hand in massive clusters 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 compute the total energy released by AGN feedback for and related it to the gas mass within . We find that the ratio of the two quantities is a strongly decreasing function of 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 BH has therefore more cold gas at its disposal.
To test the reduced frequency of BH-BH merger, we compute at the number and mass of BHs which are inside and are not bound to any substructure. We find that 80 per cent of the most massive systems (16 objects over 20) have several BHs in that central regions with a total mass greater than 10 percent of the mass of the central BH. Analysing the first mass bin (), instead, only 10 per cent of clusters have enough BHs able to account, all together, for at least 10 per cent of the mass of the central BH. 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 BHs and the central BH.
4.5 Recent growth of BH and stellar component
In our simulations the BH mass increases by a factor of between . 2018CINTHIA found instead, on the same sample of simulated objects666Note that in their work the authors selected only the most massive cluster in each Lagrangian region described in Sect. 2. However, we have checked that the growth factor of central BHs remain unchanged also considering this reduced subsample., that the central stellar component measured within 30 kpc features a significant smaller growth. This difference is due to the fact that in our simulations many substructures colliding with the BCG at 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, 2007MURANTE show 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 .
To visualize this effect, whose detailed study will require a dedicated analysis, we simply compare the evolution of the BH mass and the inner stellar component, defined as . Furthermore, we add a measure of the outer stellar component, defined as . The ‘not-bound’ identification specifies that we exclude all stars which are gravitationally bound to all substructures identified by Subfind and different from the BCG. , therefore, comprises both the ICL and the outermost stellar mass of the largest BCG in the sample. The median values of , and are computed after the normalization to their respective mass values at .
We consider 30 kpc for because we wanted to evaluate the changes on the stellar component in the immediate surrounding of the BH; furthermore, this was used in 2018CINTHIA 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 slowly increases from to , the BH mass and rapidly increase. At , indeed, the quantities are about 80 percent, 45 and 35 per cent of their final values, respectively. The remarkable agreement between the two extra definitions of , (), thin blue solid line in Fig. 12 and ), thin blue 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 BHs 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 relation for the two BH-growth channels
Given that at and the BH has grown from both channels (through gas accretion and BH-BH mergers), it is relevant to check whether the BH 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 consider both and as a function of M at (left panel) and at (right panel). As we have seen, at the gas accretion is the dominant channel. In Sect. 4.4, we saw that this channel shows only a slight increment from to . For this reason, the red points, referring to are substantially unchanged in the two panels. Viceversa, becomes the dominant component at . The results of the linear fits of the relations are reported in Table LABEL:tab:all also for the other redshifts.
From the figure, it is evident that both masses correlate well with M at both independent of which one of the two is the dominant channel from the BH mass growth. From the table, we notice that the slopes of the two relations are consistent within 1 being the 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 BH mass (see Table LABEL:tab:all).
Finally, we would like to emphasize that Fig. 7 suggests that the relation is already in place at when the BH 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.
As previously remarked the correlation between the BCG mass and the BH mass has been diffusely studied and often used to extract the mass of the BH knowing the mass of the hosting galaxy (Sect. 1). In 2018BOGDAN, the authors found in their observed sample that the scatter between the BH mass and the global cluster properties is tighter by almost 40 per cent than the scatter of the relation. Indeed, in their Table 4 they report and . As a consequence, the global properties of cluster within are also suitable to estimate the BH mass. We show in Sect. 4 that the distribution of our simulated data is in reasonable agreement with 2018BOGDAN observations and that for our simulated objects there is a clear correlation between the mass of the BHs, located at the minimum of the potential well, and the temperature or total mass of the clusters within . 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. These are defined as logarithmic differences between the actual value, generically referred as , and the expected value from the relations777The relation has previously been inverted into . () provided in Table 2:
The deviations are computed for each quantity , , , , . The panels above the diagonal refer to while those below to . The diagonal panels show the distribution of at . We use the temperature subsample whenever is considered. For each pair of deviations we list the Spearman correlation coefficients in Table 3 computed at , , 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 and is converted into a small scatter in the relation .
|0.498— 0.072 — 0.255|
|0.322 — 0.210 — 0.355||0.557 — 0.693 — 0.694|
|0.310 — 0.243 — 0.267||0.288 — 0.405 — 0.485||0.684 — 0.804 — 0.865|
|0.230 — 0.087 — 0.277||0.612 — 0.751 — 0.679||0.777 — 0.753 — 0.737||0.190 — 0.304 — 0.384|
5.1 Scatters of relations
As previously commented, the scatter is comparable to (see Table 2). As a consequence, the temperature and the total cluster mass are equally good proxy for the BH mass. The similarity between these two proxies can be explained by looking at the correlation between and in Fig. 13. The panels at both redshifts highlight how the variations of cluster temperature are not directly reflected into variations of the BH 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, we saw in Sect 4.2 that the typical delay between the increase of the 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 and all 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 BH mass due to BH-BH mergers.
When comparing, instead, the scatter of the relation between the BH mass and the cluster properties with the scatter of the relation, we find that is smaller at all times. At redshift , for example, is per cent larger, although the two scatters are in agreement between . The difference between the scatters is more significant at where typically is . In the covariance matrix formalism this translates into a correlation between and . Indeed, at 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 is significant for the presence of situations of either pre-mergers or mergers with a large impact parameter when both and (and in particular ) are smaller than the average of the sample at fixed total mass (i.e. their is negative). At , we also notice a correlation between the variations of the two quantities in the other direction (both and greater than zero) implying that the BHs that gained mass through 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 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 per cent 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 procedure888The best-fitting relations of the samples derived by excluding the clusters with recent fast accretion are:
., we find that the scatter of the relation reduces by almost 20 per cent (), it remains comparable to the re-computed scatter of the relation () and consistent with the new . 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 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 relation (Fig. 1 and also 2013RAGONE; 2018BOGDAN). On one hand the growth of simulated BHs 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 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 BH mass). In simulations, instead, the BCG and BH masses are always known and precisely defined (see Sect. 2.1 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 derived from observations are not easily reducible. The estimate of the BH 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 BH mass. This conclusion is even stronger at high redshift, such as .
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, , and global temperature, . 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 identified in the high-resolution regions of the re-simulated volumes. At , 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 (, and ). Finally, we characterized the role played by the two channels for the BH growth: accretion of gas and BH-BH merging. Our main results can be summarised as follows.
The simulated relations between SMBH mass and the global cluster properties (- and ) are in agreement with the observations by 2018BOGDAN (Fig. 2 and Fig. 4).
The relation at originates from a non-simultaneous growth of the BH and of the cluster. In particular, objects evolve on the plane either at almost constant or at almost constant (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 BH with gas and/or inducing a BH-BH merger. Clusters that recently experienced major merger events are in general below the mean relation.
The relation of the entire sample shows a degree of evolution: the slope is about 45 per cent smaller at with respect to (Fig. 7). This evolution naturally arises from hierarchical structure growth as the most massive clusters increase their mass in the period between and at a rate which is 10-20 times higher with respect to the smallest objects. Viceversa, the BH mass increases by an approximately constant factor (), independent of the mass of the hosting cluster (Fig. 6).
In our simulations, BHs grow by two different channels. Gas accretion is the most relevant channel at redshift and the only player at the earliest times. The accretion is slowed down only when the BHs are massive enough to balance gas cooling via AGN feedback. At lower redshift () one quarter of the BH mass is ascribed to mergers. From that time to the BH-BH merger contribution becomes progressively more important. Indeed, mergers contribute by about 60 percent of the total BH mass on average. When restricting the analysis to the most massive systems, we find that the accretion onto the SMBH is dominant up to , 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 at both and , independently of which BH-mass growth channel is dominant.
The rapid increase of the BH mass at 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 2018CINTHIA. Indeed, in our simulations, we found that at recent times () 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 and relations present a similar scatter (or, equivalently, does not show any significant correlation with ), meaning that they are equally valid BH mass proxy. On the other hand, is highly correlated with . As a consequence the scatter of the relation, at , is lower in our simulations, although the two values are in agreemen within . However, it is important to stress that the observed scatter of the relation is larger then the simulated one mainly for two reasons: the numerical limitation of a simplistic description of BHs growth and the large uncertainties affecting the observational measurements of both BH and BCG mass. That said, when the most dynamically active objects are discarded from the sample, the scatters of all relations, , and , 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 obtaining a large statistical sample in a cosmological environment. This means that while the large-scale inflows are resolved satellite galaxies are resolved with few smoothing lengths implying that the actual level and timing of gas accretion and BH mergers might differ from our tracked evolution. Furthermore, the modelling of the gas accretion is quite simple. Several studies (2016GASPARI for a brief review) show that the SMBH accretion in massive galaxies proceeds via chaotic cold accretion (CCA), from the kpc scale down to tens gravitational radii. While our modified Bondi accretion tends to mimic CCA in power, the frequency and rapid boosting of chaotic events are not 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 and boosted accretion can more effectively quench cooling flows. 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., 2018GASPARI, 2019TREMBLAY). 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., 2017GASPARI).