Galaxies with monstrous black holes in galaxy cluster environments

TL;DR

This study uses cosmological simulations to investigate the origins of galaxies with exceptionally large black holes, revealing early mergers and tidal stripping as key factors in their development.

Contribution

It identifies the formation pathways of extreme black hole outliers in galaxy clusters, linking early mergers and tidal effects to their growth.

Findings

Early mergers at high redshift contribute to black hole growth.

Tidal stripping influences the outlier status of galaxies.

Alpha enhancement correlates with black hole growth.

Abstract

Massive early-type galaxies follow a tight relation between the mass of their central supermassive black hole ($\rm M_{BH}$) and their stellar mass ($\rm M_{\star}$). The origin of observed positive outliers from this relation with extremely high $\rm M_{BH}$ ($> 10^{9} M_{\odot}$) remains unclear. We present a study of such outliers in the Hydrangea/C-EAGLE cosmological hydrodynamical simulations, designed to enable the study of high-mass galaxy formation and evolution in cluster environments. We find 69 $M_{\rm BH}(M_{\star})$ outliers at $z=0$, defined as those with $ \rm M_{BH} >10^{7} M_{\odot}$ and $\rm M_{BH}/\rm M_{\star}> 0.01$. This paper focusses on a sample of 5 extreme outliers, that have been selected based on their $\rm M_{BH}$ and $\rm M_{\star}$ values, which are comparable to the most recent estimates of observed positive outliers. This sample of 5 outliers, classified as `Black hole monster galaxies' (BMGs), was traced back in time to study their origin and evolution. In agreement with the results of previous simulations for lower-mass $\rm M_{BH}(\rm M_{\star})$ outliers, we find that these galaxies became outliers due to a combination of their early formation times and tidal stripping. For BMGs with $\rm M_{BH} > 10^9 M_{\odot}$, major mergers (with a stellar mass ratio of $\mu > 0.25$) at early times ($z>2$) precede the rapid growth of their supermassive BHs. Furthermore, the scatter in the relation between $\rm M_{BH}$ and stellar velocity dispersion, $\sigma$, correlates positively with the scatter in [Mg/Fe]($\sigma$). This indicates that the alpha enhancement of these galaxies, which is closely related to their star formation history, is related to the growth of their central BHs.

Tables

Table 1. Table 1: The five BMGs and their properties. Column 1 shows the BMG ID (for simulated galaxies) or name (for observations). From left to right: BH mass (column 2), stellar mass (column 3), half-mass radii for the simulated galaxies, or effective radius for observations (column 4), one dimensional stellar velocity dispersion, which is computed from the average of the three dimensional velocity dispersion (column 5), mass-weighted stellar age (column 6), and [Mg/Fe] abundance (column 7). Column 8 shows the redshift of the snapshot at which each galaxy crossed the a M BH / M ⋆ > 0.01 subscript 𝑀 BH subscript 𝑀 ⋆ 0.01 M_{\rm BH}/M_{\star}>0.01 threshold. These redshifts have an uncertainty of 500 Myr. Column 9 is the distance to the nearest more massive neighbour of each BMG at z = 0 𝑧 0 z=0 . The last two columns (10 and 11) show f strip subscript 𝑓 strip f_{\rm strip} and f age subscript 𝑓 age f_{\rm age} from equations 1 and 2 , respectively. The last four rows show the observational values for NGC1271, NGC1277 and Mrk 1216 as determined by Walsh et al. ( 2015 ) , Walsh et al. ( 2016 ) , Ferré-Mateu et al. ( 2015 ) , Ferré-Mateu et al. ( 2017 ) and Martín-Navarro et al. ( 2016 ) .

\pbox3cmName
ID	\pbox3cm log10(M${}_{\mathrm{BH}}{}^{}/M_{\odot }$)	\pbox3cmlog10($M_{\star }$/$M_{\odot }$)	\pbox3cm$R_e$
$\left[\mathrm{kpc}\right]$	\pbox3cm$\sigma$
$\left[{\mathrm{kms}}^{-1}\right]$	\pbox3cmAge
$\left[\mathrm{Gyr}\right]$	\pbox3cm[Mg/Fe]	\pbox3cmOutlier z	\pbox3cm$D_{\mathrm{host}}$
$\left[\mathrm{kpc}\right]$	\pbox3cm$f_{\mathrm{strip}}$	\pbox3cm$f_{\mathrm{age}}$
A	7.60	9.48	1.49	33	10.02	0.02	0.00	25	0.76	0.30
B	8.40	9.27	2.28	34	9.94	-0.01	0.04	30	1.49	0.38
C	9.02	9.49	1.53	55	10.40	0.40	2.35	60	1.22	0.59
D	9.32	11.04	3.93	163	11.99	0.36	0.62	475	0.03	0.59
E	9.61	11.12	2.61	214	10.96	0.31	2.35	1,150	0.43	0.55
NGC 1271	$9.48\pm 0.15$	$11.18\pm 0.20$	$2.2\pm 0.9$	276${}_{-4}{}^{}{}_{}{}^{+73}$	>10	0.23
NGC 1277	$9.69\pm 0.14$	$11.08\pm 0.14$	$1.2\pm 0.1$	333 $\pm$12	>10	0.32
Mrk 1216	9.69 $\pm$ 0.15	11.30 $\pm$ 0.17	$2.3\pm 0.1$	368 $\pm$ 3	12.8 $\pm$ 1.5	0.26 $\pm$ 0.05

Equations

$$\text{centering}{f}_{\mathrm{strip}}=\mathrm{lo}{\mathrm{g}}_{10}\left(\frac{{\mathrm{M}}_{\star ,\mathrm{init},\max }}{{\mathrm{M}}_{\star ,\mathrm{init}}}\right)\text{@add@centering}$$

$$\text{centering}{f}_{\mathrm{age}}=\mathrm{lo}{\mathrm{g}}_{10}\left[\mathrm{med}{\left(\frac{{\mathrm{M}}_{\mathrm{BH}}}{{\mathrm{M}}_{\star ,\mathrm{init}}}\right)}_{\mathrm{age},\mathrm{i}}\right]-{\alpha }_{\mathrm{med}}.\text{@add@centering}$$

Full text

††pubyear: 2017††pagerange: Galaxies with monstrous black holes in galaxy cluster environments –A

1 Introduction

The theory of galaxy formation is one of the most active topics of modern-day cosmology. To aid our understanding of the physical processes involved, many empirical scaling relations between galaxy properties have been deduced. One of the strongest of such relations is that between the mass of a supermassive black hole (BH) and the stellar mass of its host galaxy, $M_{\star}$ (for a review on supermassive BH scaling relations, see Kormendy & Ho, 2013). This strong correlation hints at a co-evolution between the central BH and the galaxy as a whole.

Recently, a handful of outliers to this relation have been observed, the most interesting being those with extremely high black hole mass ($M_{\rm BH}$) for their stellar mass. Well-studied examples are NGC 1271 (Graham et al., 2016; Walsh et al., 2015), NGC 1277 (van den Bosch et al., 2012; Trujillo et al., 2014; Walsh et al., 2016), and Mrk 1216 (Ferré-Mateu et al., 2017). The origin (and even the existence) of such over-massive BHs is highly debated. For example, Savorgnan & Graham (2016) and Graham et al. (2016) both argue that when spheroid/disc decompositions are modelled under different assumptions, the $M_{\rm BH}$ values of e.g. NGC 1277 and NGC 1271 can be brought to within 2-$\sigma$ of the observed $M_{\rm BH}(M_{\star})$ relation, and are no longer anomalous. It is therefore interesting to test whether such $M_{\rm BH}(M_{\star})$ outliers are expected to exist in current galaxy formation theories, which can most accurately be tested via cosmological simulations.

Barber et al. (2016, hereafter B16) studied positive outliers of the $M_{\rm BH}(M_{\star})$ relation in the (100 Mpc)3-volume EAGLE simulation (Schaye et al., 2015; Crain et al., 2015). They confirmed that such outliers do indeed exist, and arise due to their early formation time and/or the tidal stripping of stars. The importance of tidal stripping was also pointed out by Volonteri et al. (2008) and Volonteri et al. (2016). However, Barber et al. (2016) did not find outliers with $M_{\rm BH}>10^{8.5}\rm M_{\sun}$ although, as mentioned above, such BH outliers have been observed. Plausibly, the lack of high-mass outliers in B16, is due to the limited volume of the EAGLE simulation, as Saulder et al. (2015) found that the frequency of massive compact high stellar velocity dispersion galaxies (such as NGC1271 and NGC1271) is $\sim 10^{-7}$ galaxies per $\rm Mpc^{-3}$, implying that $<1$ of such galaxies should be found in the (100 Mpc)3 volume studied by B16. A larger simulation volume would thus be required to find $M_{\rm BH}(M_{\star})$ outliers with such high $M_{\rm BH}$. In addition, simulations of representative volumes, such as EAGLE, are dominated by field galaxies, while we expect to find more high-mass BHs (and thus also $M_{\rm BH}(M_{\star})$ outliers with $M_{\rm BH}$ of > $10^{8.5}\rm M_{\sun}$) in galaxy cluster environments due to the higher densities and galaxy masses. Moreover, tidal stripping is expected to be enhanced in such galaxy cluster environments due to the stronger tidal forces.

The Cluster-EAGLE (C-EAGLE) project (Bahé et al., 2017; Barnes et al., 2017), a suite of hydrodynamical, cosmological zoom-in simulations of galaxy clusters and their surrounding environments, is ideal for the purpose of this research. This work aims to study the origin and properties of the $M_{\rm BH}(M_{\star})$ outliers in the cluster environments of the Hydrangea simulations (Bahé et al., 2017, a subset of the C-EAGLE simulations, for which the high-resolution zoom-in region extends to 10 times the virial radius) in a stellar mass regime closer to that of observed galaxies with overly-massive BHs. We compare our results to those from B16 for the (100 Mpc)3 EAGLE volume and investigate the physical origin of galaxies with extreme BH masses. When outliers with $M_{\rm BH}>10^{8.5}\rm M_{\odot}$ are found in Hydrangea, it is interesting to see whether the origin of their anomalous $M_{\rm BH}(M_{\star})$ could be explained by the same arguments as B16, or whether additional reasoning is required. Recently, McAlpine et al. (2018) and Steinborn et al. (2018) investigated the effect of major mergers on triggering the rapid growth phase of central black holes. Major mergers are defined as $M_{\star,1}/M_{\star,2}=\mu>0.25$, while intermediate mergers are $0.1<\mu<0.25$, where $M_{\star,1}$ and $M_{\star,2}$ are the the stellar masses of the less and more massive progenitor of the galaxy in question. They both argue that major galaxy-galaxy interactions are not statistically dominant mechanisms for fuelling nuclear activity or triggering the rapid growth phase at $z>2$. We therefore investigate whether our most extremely positive outliers (with $M_{\rm BH}>10^{9}\rm M_{\odot}$) have experienced any major mergers.

This paper is organized as follows: we discuss the Hydrangea simulations and our galaxy tracing methods in Section 2. The simulated $M_{\rm BH}(M_{\star})$ distribution at $z=0$ is evaluated in Section 3. With the purpose of comparing the simulated BMGs to recently observed outliers, we continue with a more detailed evaluation of the 5 most interesting BMGs in Section 4, where we also consider the effects of tidal stripping and early formation time on $M_{\rm BH}(M_{\star})$ outliers in Section 4.1. We conclude with a short discussion on the [Mg/Fe] ratio of galaxies in Hydrangea (subsection 4.4). Our main findings are summarized in Section 5.

Throughout this paper, the following cosmological parameters, consistent with Schaye et al. (2015) and Planck Collaboration et al. (2014) are assumed: $\Omega_{\rm m}$ = 0.307, $\Omega_{\rm b}$ = 0.04825, $\Omega_{\Lambda}$ = 0.693, $\sigma_{8}$ = 0.8344, and $h=\rm H_{0}$ /(100 km s*-1* Mpc*-1*) $=0.6777$.

2 Simulations and Analysis

2.1 The Hydrangea simulations

The EAGLE project is a suite of hydrodynamical simulations aimed at understanding how galaxies form and evolve from shortly after the Big Bang to the present day (Schaye et al., 2015; Crain et al., 2015). The public release is described in McAlpine et al. (2016). The C-EAGLE (short for Cluster-EAGLE) project (Barnes et al., 2017; Bahé et al., 2017) is a suite of cosmological hydrodynamical zoom-in simulations that adopts the EAGLE galaxy formation model to study the formation and evolution of 30 galaxy clusters with masses ranging from111$M_{200}$ is defined as the mass within $r_{200}$, which is the the radius at which the mean enclosed mass density $\rho$ is equal to 200 times the critical density ($\rho_{\rm crit}=3H^{2}/8\pi G$) of the Universe. $M_{200}=10^{14}$ to $\approx 2\times 10^{15}\rm M_{\odot}$. The Hydrangea project (Bahé et al., 2017) is a core part of the C-EAGLE project, which consists of 24 of the 30 total galaxy clusters, simulated with high-resolution zoom-in regions extending out to 10$r_{200}$ from each cluster’s centre. Particles in the re-simulated higher-resolution regions initially have masses identical to those in the intermediate-resolution EAGLE simulations: $1.81\times 10^{6}\rm M_{\odot}$ for gas particles and $9.7\times 10^{6}\rm M_{\odot}$ for dark matter particles.

Particle information was saved in regular time step intervals of 500 Myr in 28 snapshots, running from $z=14.0$ to $z=0$. Additionally, snapshots 19 ($z=0.37$) and 26 ($z=0.1$) were created to enable comparison to the EAGLE simulation. Furthermore, between each of the 28 main snapshots, three ‘snipshots’ were stored, reducing the time resolution to $t=125\rm Myr$. These snipshots contain only the most important, and most rapidly time-varying, quantities, such as particle positions and velocities (similar to EAGLE; see Schaye et al., 2015). For the 1 Gyr intervals at lookback times of 0–1, 5–6, and 8–9 Gyr, 19 additional snipshots were stored, boosting the time resolution to $t=25\rm Myr$ for those epochs.

A full description of the model is given by Schaye et al. (2015), while a complete overview of the C-EAGLE project and Hydrangea simulations is given by Barnes et al. (2017) and Bahé et al. (2017), respectively. Here we will summarize the aspects that are most relevant for our study.

Star particles are formed stochastically from gas particles once their density exceeds the metallicity-dependent threshold of Schaye (2004). The star formation rate is obtained by converting the Kennicutt-Schmidt law into a pressure law (Schaye & Dalla Vecchia, 2008), and stellar mass loss from massive stars, AGB stars and type Ia and type II supernovae is simulated by the redistribution of stellar mass and metals among neighbouring gas particles (Wiersma et al., 2009b). The mass of the individual elements lost through stellar feedback processes (i.e. winds from AGB stars, massive stars, and core collapse supernovae) is computed using the nucleosynthetic yields and metallicity-dependent lifetimes of Marigo (2001) and Portinari et al. (1998). Stellar energetic feedback is implemented by injecting thermal energy into the neighbouring gas particles, following Dalla Vecchia & Schaye (2012). During each stochastic feedback event, the temperature increase of the receiving gas is kept fixed, while the number of influenced gas particles depends on the local gas density and metallicity, and was calibrated to reproduce the observed $z=0.1$ galaxy stellar mass function and galaxy sizes (Crain et al., 2015). The elements most important for radiative cooling (H, He, C, N, O, Ne, Mg, Si, and Fe; Wiersma et al., 2009a) are traced individually. Metallicities mentioned in this work refer to stellar metallicity. For SN Type Ia, yields are taken from the W7 model of Thielemann et al. (2003).

Following Springel et al. (2005), a BH is seeded in a Friends of Friends (FOF) halo that does not already contain a BH, once its total mass exceeds $10^{10}\rm M_{\odot}$ $h^{-1}$. The gas particle with the highest density is converted into a BH particle with an initial sub-grid seed mass of $10^{5}\rm M_{\odot}$ $h^{-1}$. $M_{\rm BH}$ can subsequently grow through gas accretion and BH-BH merging (Schaye et al., 2015; Rosas-Guevara et al., 2015). The accretion model is based on a sub-grid accretion disc with an adjustable viscosity parameter, $C_{\rm visc}$, that influences the $M_{\rm BH}(M_{\star})$ relation (Crain et al., 2015), and was calibrated to reproduce the galaxy stellar mass function (GSMF). To model (unresolved) dynamical friction effects, BHs with $M_{\rm BH}$ < 100 $m_{\rm gas}$ are forced to migrate towards the minimum of the gravitational potential of the halo.

BHs merge when they are spatially separated by less than the BH smoothing kernel and three times the gravitational softening length, and their relative velocity is smaller than the circular velocity at the smoothing length of the most massive BH of the pair. The latter condition prevents BHs from merging during flybys.

C-EAGLE uses EAGLE’s ‘AGNdT9’ model to simulate AGN feedback, which differs from the reference model used in EAGLE only by a different value of $\Delta T_{\rm AGN}$ (= $10^{9}$ K instead of $10^{8.5}$ in EAGLE) and $C_{\rm visc}$(= $2\pi\times 10^{2}$), since this produces more realistic hot gas haloes of galaxy groups. This reduces numerical cooling losses but still results in somewhat too large hot gas mass fractions in clusters (Barnes et al., 2017). AGN-feedback is implemented by stochastically heating the surrounding gas particles, where the stochastic probability is chosen to ensure that on average the energy injection rate is $\epsilon_{r}\epsilon_{f}\dot{m}_{\rm acc,BH}c^{2}$, where $\epsilon_{r}=0.1$ is the radiative efficiency of the accretion disc and $\epsilon_{f}=0.15$ is the fraction of the radiated energy that is coupled to the interstellar medium. The efficiency of BH feedback was calibrated to yield a match to the observed $M_{\rm BH}(M_{\star})$ relation at $z\approx 0$.

2.2 Galaxy identification and corrections

Dark matter haloes are identified using the FOF algorithm on dark matter particles with linking length $0.2$ times the mean inter-particle spacing (Davis et al., 1985). Baryon particles are then attached to the halo of the nearest DM particle (if this DM particle belongs to any halo). The SUBFIND algorithm (Springel et al., 2001; Dolag et al., 2009) is used subsequently to identify gravitationally self-bound substructures within these FOF haloes, i.e. objects that can be identified as galaxies. Subhaloes of fewer than 20 particles are not catalogued as they are poorly sampled. The subhalo in the FOF-group that contains the particle with the deepest gravitational potential is defined as the central galaxy of that group, the other subhaloes are labelled as satellites. We define a subhalo’s stellar mass, $M_{\star}$, as the sum over all stellar particles bound to the subhalo. 222When using $M_{\star,30}$, defined as the sum over all stellar particles bound to the subhalo which are within a 30 proper kpc aperture centred on the potential minimum, the simulations are consistent with the observational definition of stellar mass based on Petrosian radii (Schaye et al., 2015; Furlong et al., 2015). For the majority of the extreme $M_{\rm BH}$ outliers studied in this work, $M_{\star,30}$ and $M_{\star}$ are equal at $z=0$. The maximum difference encountered was $\rm log_{10}M_{\star}/\rm M_{\odot}$ - $\rm log_{10}M_{\star,30}/\rm M_{\odot}$ $\lesssim 10^{-2}$. To avoid contamination from low resolution regions of the simulation, we only include galaxies which have always been > 2 cMpc away from a low resolution boundary particle at $z=0$.

2.2.1 Subhalo corrections

Since SUBFIND identifies galaxies as bound structures in the simulation, it can occur that the central BH in a subhalo temporarily becomes its own bound structure with a small fraction of the stars in the central region. This results in subhaloes seeming to spuriously pop in and out of existence over time (B16). While such structures may indeed be self-bound, they would not be identified as individual galaxies in observations, so it is crucial to this work to remove such structures from the data by merging them back into the surrounding galaxy. As in Schaye et al. (2015), subhaloes separated by a distance less than 3 proper kpc and smaller than the half-mass radius of the star particles were merged in post-processing.

The SUBFIND algorithm is known to occasionally assign BHs to the host of a satellite galaxy instead of to the galaxy in whose centre it resides (B16). In this case, the satellite is typically left without any black hole particles, causing its $M_{\rm BH}$ to temporarily drop to zero, typically during pericentric passages while orbiting its central. We correct for this BH misassignment following B16. We select candidate subhaloes for each BH particle at each snapshot, at a distance $d$ smaller than the half-mass radius of the relevant subhalo and $d<3$ kpc (in proper distance). If the BH was not part of the sub-group of one of these candidates, we selected the most massive candidate as the correct subhalo for the BH and reassigned the particle to that subhalo. This reassignment affected $\approx 0.5-1\%$ of the BH particles per snapshot, but had no effect on the final sample of BMG-galaxies.

2.2.2 BMG selection

In order to define (positive) $M_{\rm BH}(M_{\star})$ outliers as galaxies with a central SMBH that is much more massive (or equivalently have a stellar mass that is much lower) than expected relative to the median $M_{\rm BH}-M_{\star}$ relation at $z=0$, we select galaxies with $M_{\rm BH}/M_{\star}>0.01$ and $M_{\rm BH}>10^{7}\rm M_{\odot}$. The former criterion was chosen to select only the strongest outliers in the simulation, while the latter avoids resolution effects caused by seed mass ($10^{5}h^{-1}M_{\sun}$) BHs. This definition differs from that of B16, who defined outliers as having $M_{\rm BH}$ more than 1.5 dex above the median $M_{\rm BH}(M_{\star})$ relation in the simulation. When using the definition from B16, the fraction of all galaxies with $M_{\rm BH}>10^{7}\rm M_{\odot}$ at $z=0$ that would be classified at outliers in EAGLE and Hydrangea respectively amounts to $0.075\%$ and $10\%$. Our new definition of being an outlier corresponds to $0.035\%$ and $1.8\%$ of the galaxies with $M_{\rm BH}>10^{7}\rm M_{\odot}$ at $z=0$ for EAGLE and Hydrangea respectively. This new definition thus prevents a large fraction of galaxies from being classified as outliers, while still including the galaxies with high $M_{\rm BH}$ (> $10^{9}\rm M_{\odot}$). For comparison, the value of the median $M_{\rm BH}/M_{\star}\approx$ 0.001 at $z=0$ for $M_{\star}\approx 10^{11}\rm M_{\odot}$.

Since our new definition still results in quite a large set (69) of $M_{\rm BH}(M_{\star})$ outliers, we manually select the most extreme positive outliers of interest for further analysis, i.e. those with the highest ratio of $M_{\rm BH}/M_{\star}$. We selected 16 of the most extreme outliers of the $M_{\rm BH}(M_{\star})$ relation, based on a visual inspection of their locations in the $M_{\rm BH}-M_{\star}$ plane (see Figure 1), which will hereafter be referred to as “Black hole monster galaxies” (BMGs) (using the term first coined by Kormendy & Ho, 2013).

The stellar mass surface density images of the main progenitors of each of these BMGs were carefully inspected at each snapshot by means of the merger tree algorithm as described below in Section 2.2.3. Additionally, the stellar surface mass density images of the region around each main progenitor were studied at each output, exploiting the high time resolution of the combined snap- and snipshots.

2.2.3 Tracing the main progenitors

The main progenitors of these “BMGs” were primarily traced using the “merger tree” algorithm from Bahé et al. (in prep.), which is similar to that used for EAGLE (McAlpine et al., 2016; Qu et al., 2017). In essence, this algorithm identifies candidate descendants of a subhalo in the subsequent snapshot by following its five per cent most bound collisionless particles (i.e. DM, stars, and BHs) as tracers. The actual descendant subhalo is then chosen according to the fraction of tracer particles it contains. In the case of multiple candidate progenitors for the same descendant, the descendant subhalo is chosen according to the total mass contributed by each of them.

The BMGs typically experience complex interactions with their environment, making it difficult to trace their progenitors between snapshots by means of a single algorithm. Hence we confirmed the main progenitors based on an additional visual analysis of the stellar density profiles. The $M_{\star}$, $M_{\rm BH}$ and $M_{\rm DM}$ surface density profiles of the BMGs and their environment were visualized at each output (including snipshots), so that they could be subjected to careful inspection by eye (see Figure 8 for a smoothed example of such an image). For some complicated events, the tracing algorithm did not yield the correct progenitor masses. For example, this can occur if the BH particle becomes unbound between snapshots, or when a large fraction of the stellar or DM particles was temporarily (for one or two snapshots) bound to a more massive neighbouring galaxy. A manual inspection of the stellar density, gas density and the BH particles of the region (of up to $\sim$150 pkpc) surrounding the subhalo in question was therefore decisive in determining the correct progenitor mass in such cases.

It was found that 9 of the 16 BMG-subhaloes lie very close (at a distance $d<15$ pkpc) to a more massive neighbour at $z=0$. These 9 subhaloes fell just outside of our merger and BH reassignment requirements of $d<3$ pkpc, but were clearly an artefact of the SUBFIND algorithm or in the final stages of a galaxy-galaxy merger as described above. One additional subhalo was classified as being in the final stages of a merging event between two massive ($M_{\star}\sim 10^{11}$ M*☉*) galaxies. These high masses have a larger gravitational influence radius and therefore reside at a larger distance (>50 kpc) from its merging companion. Lastly, one subhalo was found to consist of only 23 baryonic particles at $z=0$, where there was no significant stellar over-density at the location of this subhalo. It was therefore classified as a naked/ejected BH without a real galaxy hosting it.

We concluded that these 11 subhaloes would not be considered an BMG in the observational sense and should therefore not be considered as individual galaxies. We thus excluded them from further analysis, and focus on the remaining 5 BMGs. 333Videos of the stellar density profiles, including an extended description of each of the initial sample of 16 subhaloes can be found at https://sites.google.com/view/vson-umbh-galaxies/homepage

3 Resulting $M_{\rm BH}(M_{\star})$ Outliers at $z=0$

Figure 1 shows the $M_{\rm BH}-M_{\star}$ plane for Hydrangea central and satellite galaxies. A red dashed line marks our threshold for selecting $M_{\rm BH}(M_{\star})$ outliers as described in Section 2.2.2. At $z=0$, 77$\%$ of the $M_{\rm BH}(M_{\star})$ outliers are satellites and typically have $M_{\star}=[10^{9}-10^{11}]\rm M_{\odot}$. The median relation, shown as a solid red line, is lower for satellite galaxies than the median relation for centrals at $M_{\star}\approx 10^{12}\rm M_{\odot}$. This is caused by a small subset ($\sim 20$) of ‘undermassive’ BH galaxies (with $M_{\rm BH}\approx 10^{8.5}\rm M_{\odot}$) that are most likely caused by the ejection of black holes through tidal encounters of similar-mass black holes. Throughout this paper, we compare outliers to the total median $M_{\rm BH}(M_{\star})$ for both centrals and satellites.

The 5 manually selected BMGs are labelled in Figure 1 with IDs corresponding to those in table 1. They were selected with a preference for a high $M_{\rm BH}$ because of the similarity to the extremely high $M_{\rm BH}$ of several observed $M_{\rm BH}(M_{\star})$ outliers. A bias towards higher $M_{\rm BH}$ at fixed $M_{\star}$ also leads to a bias towards older galaxies; as the gradient in stellar age in Figure 1 shows, at fixed $M_{\star}$, galaxies with higher $M_{\rm BH}$ tend to be older. There is a larger variation in $M_{\star}$ at fixed $M_{\rm BH}$ for satellite galaxies, with respect to that for centrals at $M_{\star}\sim 10^{10}\rm M_{\odot}$. This could be an indicator of tidal stripping in such cluster environments. However, it could also plausibly be caused by the quenching of star formation in satellites as described by e.g. Bahé et al. (2017).

In Figure 1, we also show histograms of the $M_{\rm BH}$ and $M_{\star}$ distributions of both the (100 Mpc)3 EAGLE volume and the Hydrangea simulations. Unsurprisingly, the Hydrangea simulations contain more galaxies than EAGLE at high $M_{\star}$, and therefore also of supermassive BHs. Because of this larger number of massive galaxies, they are also more likely to contain $M_{\rm BH}(M_{\star})$ outliers at this high-mass end. Outliers with the highest $M_{\rm BH}$ (outliers D and E) have masses comparable to the most recent estimates for NGC 1277 and NGC 1271 (Walsh et al., 2016, 2015). For comparison, the outliers identified by B16 in the (100 Mpc)3 EAGLE simulations that satisfy our requirements of being a $M_{\rm BH}(M_{\star})$ outlier are shown in Figure 1 by star-shaped symbols and are labelled with the outlier IDs from B16. We remind the reader that B16 identified outliers as galaxies that lie 1.5 dex above the median $M_{\rm BH}(M_{\star})$ relation in the simulation at $z=0$, which differs from our threshold (see Section 2.2.2).

The distribution of galaxies in the $M_{\rm BH}-M_{\star}$ plane can be divided into two general regions. As described by Bower et al. (2017), above a halo mass of $10^{12}\rm M_{\odot}$, stellar feedback is no longer sufficient to prevent gas infall towards the central BH once $M_{\star}$ reaches $\sim$ $10^{10.5}\rm M_{\odot}$. The BH subsequently begins a phase of runaway growth during which it accretes gas until it is massive enough to regulate its own growth (this process is described in more detail in Weinberger et al., 2017; Bower et al., 2017, and in Section 4.1). The rapid growth phase causes a nearly vertical trend at $M_{\star}$ $\approx$ $10^{10.5}\rm M_{\odot}$. Once the BH is massive enough to regulate its own growth via AGN feedback, the galaxy follows the shallower high-mass 1:1 relation (i.e. parallel to the gray dashed lines), where growth in $M_{\star}$ and $M_{\rm BH}$ is increasingly due to mergers.

As an extension to the work done by B16, we find from the $M_{\rm BH}-M_{\star}$ distribution in Figure 1 that strong $M_{\rm BH}(M_{\star})$ outliers with $M_{\rm BH}$ greater than $10^{9}\rm M_{\odot}$ are expected to exist in cluster environments. These high values are comparable to the most massive observed $M_{\rm BH}(M_{\star})$ outliers. In the next Section we will study the origin of these BMGs, and compare this to results from B16.

4 Evolution and origin of black hole monster galaxies

As described in Section 2.2.2, we manually select 5 outliers of interest based on a visual inspection of their locations in the $M_{\rm BH}-M_{\star}$ plane (see Figure 1). These 5 outliers, or BMGs, were selected based on their extreme BH masses and similarity to observed $M_{\rm BH}$ outliers. We have traced their main progenitors back in time to study their origin and evolution.

In Figure 2 we show the evolution of the main progenitors for the five BMGs in the $M_{\rm BH}-M_{\star}$init plane, where $M_{\star,\rm init}$, the initial stellar mass, is the the total stellar mass using the mass of each star particle at birth, thereby excluding mass loss due to stellar evolution. Individual snapshots are indicated with circles. The typical movement of BMGs in the diagram over time is indicated by black arrows and the last snapshot ($z=0$) is indicated with a star. These tracks were constructed by examining the merger trees and the stellar mass density profiles of all BMGs, as described in Section 2.2.3.

4.1 Tidal stripping and early formation

We would like to know whether the physical mechanisms by which these galaxies become such extreme outliers of the $M_{\rm BH}(M_{\star})$ relation are consistent with those found by B16. We therefore investigate the relative contributions of tidal stripping and early formation time to their $M_{\rm BH}/M_{\star}$ ratio.

The outliers are in general relatively compact with respect to the complete sample of galaxies at $z=0$. That is, $78\%$ of the outliers have a stellar half-mass radius below the median at the same $M_{\star}$. Furthermore, we compared the distances of satellite outliers to their nearest more massive neighbour. The outliers lie significantly closer to their more massive neighbours than other satellites with $M_{\star}>10^{9}\rm M_{\odot}$, as shown by a KS-test with $p\ll 0.01$. These results point towards tidal stripping as an important contributor in the process of galaxies becoming outliers.

To investigate the extent to which stripping of stellar material due to tidal forces has influenced our BMGs, we compute the maximum initial stellar mass ($M_{\star,\rm init,max}$) that each BMG reached during its evolution. The unbinding of star particles is the only effect that causes galaxies to lose initial stellar mass, making the difference between $M_{\rm\star,\rm init,max}$ and the $M_{\star,\rm init}$ (at $z=0$) a good measure of stellar stripping. Following B16, we quantify the impact of stripping by the logarithmic ratio of $M_{\star,\rm init,max}$ and $M_{\star,\rm init}$ at $z=0$:

[TABLE]

The values of $f_{\rm strip}$ are shown in Table 1 for all of our BMGs and are also visualized by the coloured lines in Figure 3. The BMGs with the largest $f_{\rm strip}$ (IDs B and C in Table 1) have lost an order of magnitude in $M_{\star}$. 444BMG A would not be considered an outlier when evaluated at its initial stellar mass at $z=0$. Though it was pushed to the edge of our outlier definition by tidal stripping of its $M_{\star}$, it was a additional loss in $M_{\star}$ due to stellar evolution that made it an outlier. The horizontal movement along the $M_{\star}$ plane in Figure 3 shows a consistent result to B16: most of the BMGs are heavily affected by tidal stripping. However galaxies C, D and E are still considered outliers (according to our $z=0$ definition of being an outlier), even after removing the effect of tidal stripping (i.e. considering their maximum $M_{\star,\rm init}$).

B16 found that a secondary factor causing galaxies to become outliers of the $M_{\rm BH}(M_{\star})$ relation is their stellar formation time. At earlier times, the normalization of the $M_{\rm BH}(M_{\star})$ relation was higher. This difference in normalization originates from the difference in the density of the universe which causes galaxies at fixed $M_{\star}$ to have deeper potential wells (Bower et al., 2017).

When star formation no longer provides enough feedback to regulate the infall of circumgalactic gas onto the galaxy, the central BH enters a rapid growth phase. Black holes have a theoretical maximum accretion rate (the Eddington rate) which scales as $\dot{M}_{\rm BH}\propto M_{\rm BH}$. However, in practice gas densities are low enough that BH accretion proceeds according to the Bondi-Hoyle rate, which scales as $\dot{M}_{\rm BH}\propto M_{\rm BH}^{2}$, which causes a highly non-linear accretion phase. The final $M_{\rm BH}$ depends on the binding energy of the halo, thus at higher $z$ the BH has to grow to a higher mass before it becomes self-regulating (Booth & Schaye, 2010, 2011). Isolated galaxies that formed at high $z$ therefore become outliers at $z=0$, while the rest of the population lowers the normalization. McAlpine et al. (2018) also confirm that the more massive the BH is today, the earlier it began its rapid growth phase (typically $z\approx 6$ for $M_{\rm BH}$ $10^{9}\rm M_{\odot}$). Note that our definition of an outlier exclusively refers to $z=0$: our $1\%$ cut-off is not redshift-dependent.

To test the effect of formation time on our 5 BMGs, we plot in Figure 4 the $M_{\rm BH}/M_{\star,\rm init}$ ratio against the mean stellar ages (weighted by their initial stellar mass) of central galaxies with $M_{\rm BH}>10^{7}\rm M_{\odot}$. Here we assume that the majority of central galaxies at $z=0$ have not been affected significantly by stellar stripping and only include galaxies with $M_{\rm BH}>10^{7}\rm M_{\odot}$, where most galaxies have already been through the rapid BH growth phase, as seen in Figure 1. The BMGs are plotted using their maximum $M_{\rm\star,\rm init}$ to exclude stripping effects. The median of the $M_{\rm BH}/M_{\star,\rm init}-{\rm age}$ relation (black dashed line) is a measure of the contribution that the age of a galaxy has on the scatter in the $M_{\rm BH}-M_{\star}$ relation. Here we see a clear positive trend of $M_{\rm BH}/M_{\star}$ ratio with mean stellar age. This trend represents the evolution of the normalization of the $M_{\rm BH}(M_{\star})$ relation over time.

As in B16, we quantify the contribution of age in a similar manner to that of the tidal stripping contribution.

[TABLE]

For a given BMG, $f_{\rm age}$ is the logarithmic difference between the value of the median $M_{\rm BH}/M_{\rm\star,init,max}$ - age relation measured at the stellar age of the galaxy (the black dashed line in Figure 4) and $\alpha_{\rm med}$ = -3.15, the median $M_{\rm BH}/M_{\star,init}$ value measured for the whole population of central galaxies with $M_{\rm BH}>10^{7}\rm M_{\odot}$ (the dashed red line in Figure 4). The values of $f_{\rm age}$ are shown in Table 1 for all five of the BMGs. They are illustrated as vertical coloured lines in Figure 4, which can easily be compared to the stripping effects for each of these BMGs (horizontal coloured lines) in Figure 3. From Figures 3 and 4, it can be seen that those BMGs (C, D and E) which could not be explained by tidal stripping have formed at relatively early times. This shows that the effects found by B16 also apply for to the origin of $M_{\rm BH}(M_{\star})$ outliers with $M_{\rm BH}>10^{9}\rm M_{\odot}$. However, these BMGs still lie on the outer edges of the distribution (at $>1\sigma$ from the median relation), even after the effects of both stripping and high stellar age are accounted for. In the next section we study these objects in more detail to better understand their origin.

4.2 Detailed analysis of selected BMGs

The formation and evolution of the 5 BMGs (selected based on their resemblance to observational outliers, see Section 2.2.2 for the full definition) were investigated in more detail by studying their morphology and environments through visual inspection of the density of stellar, and BH particles throughout the simulation. An example of such a stellar density profile is shown in Appendix A where we depict the outlier ID E at redshifts $z=0.68,0.57$ and $0.00$.

The properties of the BMGs are summarized in Table 1. The majority of these galaxies have lost a significant fraction of their maximum stellar mass due to tidal stripping. BMG IDs A and B lie $<100$ kpc from their nearest more massive neighbour at $z=0$. They have crossed above the $M_{\rm BH}/M_{\star}>0.01$ threshold between $z=0.036$ and $z=0$, which leads us to believe that they are still in the process of being tidally stripped. Their stellar density plots and $f_{\rm strip}$ values (0.76 and 1.49 respectively) bot confirm that BMGs A and B are remnants of tidal stripping.

The BMGs with highest $M_{\rm BH}$ are of particular interest as they have similar $M_{\rm BH}$ to observed outliers of interest, since the extraordinarily high $M_{\rm BH}/M_{\star}$ ratios of these observed galaxies are still debated (Savorgnan & Graham, 2016; Graham et al., 2016). The early formation time dominantly contributes to the outlier origins of BMGs D and E. For BMG C, it is tidal stripping that is the dominant process, although its early formation time also strongly contributes to its outlier origin. However, after removing the effects of stripping and age, the $M_{\rm BH}/M_{\star}$ ratios of BMGs are still on the outskirts of the distribution (see Figure 4). We will now discuss the origin of these three BMGs in detail.

BMG-D, in contrast to BMGs A and B, has barely been influenced by tidal stripping, which is reflected in its value of $f_{\rm strip}=0.03$ (see also Figures 2 and 3). This galaxy hosts a very old BH that entered the rapid growth phase directly following a major merger event at $z=5$, reaching $M_{\rm BH}$ $>10^{7}\rm M_{\odot}$ by $z=4$. It enters a second rapid growth phase after another major merger event at $z\approx 2.5$, causing $M_{\rm BH}$ to increase to $\approx 2\times 10^{9}$ $M_{\sun}$, within less than 125 Myr. This triggers an AGN feedback blast wave that removes all the gas from the galaxy’s environment, thereby quenching star formation (around $z=2$). The galaxy resides on the outskirts of a cluster. It has passed through this cluster around $z=0.6$, at which point it crossed the $M_{\rm BH}/M_{\star}>0.01$ threshold. At $z=0$, the galaxy has relaxed and it is left as an isolated massive galaxy. We therefore classify it as a relic of the high-redshift ($z>2$) Hydrangea Universe.

Similar to BMG D, BMG C hosts a very old central BH. This central BH entered the rapid growth phase between $z=6.7$ and $z=5.5$, immediately following a merger event involving two galaxies with near equal $M_{\star}$. By $z=3$ it has reached $M_{\rm BH}\approx 10^{8.9}$ M*☉*. It continues to grow in $M_{\rm BH}$ via gas accretion and a second major merger event at $z=2$, reaching its maximum $M_{\rm BH}$ of $10^{9.0}\rm M_{\odot}$ by $z\approx 1.9$. At redshift $z=0.6$ it encounters a Brightest Cluster Galaxy (BCG), which starts to strip its stellar mass. This stripping continues up to $z=0$, where we see a clear ‘nugget’ stellar core on the outskirts of the BCG. We therefore classify this galaxy as a severely stripped relic of the high-redshift ($z>2$) Universe.

The most notable of the BMGs is outlier ID E, since it is not only the BMG with the highest $M_{\rm BH}$ but it is also the only BMG that is classified as a central galaxy at $z=0$. Like BMGs C and D, it contains an old central BH. An intermediate merger at $z=2.7$ triggers the rapid growth phase of its central BH, increasing its black hole mass to $M_{\rm BH}=10^{9.45}\rm M_{\odot}$. The corresponding AGN feedback blast has cleared all the gas from the galaxy by $z=2.0$, and the galaxy passes the $M_{\rm BH}/M_{\star}>0.01$ limit around the same redshift. A third and fourth intermediate merger occur at $z=2.3$ and $z=1.3$, where the merging companions host black holes of $M_{\rm BH}\approx 10^{7.4}\rm M_{\odot}$ and $M_{\rm BH}=10^{7.9}\rm M_{\odot}$ respectively. The latter two mergers increase the initial stellar mass to $M_{\star}=10^{11.5}\rm M_{\odot}$ and $M_{\star}=10^{11.7}\rm M_{\odot}$ respectively. These merging events are visible in Figure 2 as the ‘staircase’ shape that BMG E (solid red line) follows at high $M_{\rm BH}$.

From $z=0.7$ onwards, BMG E interacts with a nearby proto-cluster, which causes the outer stellar halo to completely disconnect from the galaxy-core (the central stars and BH) at $z\approx 0.6$. This event decreased its half-mass radius by $83\%$ and causes BMG E to lose 65% of its stellar mass (this event is shown in the left-hand and middle panels of Figure 8). The outer halo, instead of being shredded, continues to exist as a separate, core-less subhalo. The central BH and the stellar core remain as a compact, fairly isolated object that escapes from this proto-cluster and has time to re-adjust towards an equilibrium state such that it does not appear particularly disturbed at $z=0$ (see the right-hand panel of Figure 8). Although this would be considered an unprobable event, it is not necessarily unphysical. We thus conclude that although this galaxy is an isolated (with its nearest more massive neighbouring subhalo at a distance of > 1.15 Mpc) central at $z=0$, it has experienced significant environmental influence as a satellite in the past. Similar to BMG C, this galaxy is therefore classified as a severely stripped relic of the high-redshift ($z>2$) Universe.

BMG C, D and E became outliers immediately after their BH went through the rapid growth phase at very high redshift ($z\gtrsim 2$; see Table 1). All three of these galaxies have furthermore experienced multiple merger events at high redshift ($z>2$). Since these BMGs entered the rapid growth phase within 500 Myrs of a merger event, it is possible that these merger events created initial conditions of the rapid growth phase that drove their $M_{\rm BH}$ to anomalously high values. McAlpine et al. (2018) predict that galaxy interactions, while relevant at low $0<z<2$, are less important for triggering the rapid growth phase at high redshifts ($z>2$). Steinborn et al. (2018) investigated the effect of major mergers on the AGN luminosity in the Magneticum Pathfinder simulations. They found that $>50\%$ of the luminous AGN ($L_{\rm AGN}>10^{45}$ erg/s) at $z\approx 2$ have experienced a merger within the last 0.5 Gyr. However, they argue that, to some extent, this result reflects the intrinsically high merger rates of massive galaxies in which luminous galaxies typically reside. Steinborn et al. (2018) furthermore found that although mergers increase the probability for nuclear activity by a factor of three, they still play only a minor role ($<20\%$) in causing nuclear activity in the overall AGN population. These major merger events immediately preceding the rapid growth phase at $z>2$ would thus be considered rare, or at least uncommon coincidences. It is therefore plausible that if major mergers fuelled the central BHs of BMGs C, D and E, this drove them to anomalously high BH masses, which in turn causes them to still be offset in Figure 4.

4.3 Comparison to observations

It is interesting to know whether the origins of the simulated BMGs could be used to explain the origins of observed $M_{\rm BH}(M_{\star})$ outliers. Compact positive outliers from the $M_{\rm BH}-M_{\star}$ relation are observationally characterized by following the $M_{\rm BH}(\sigma)$ relation, and by being offset from the Faber-Jackson relation (Walsh et al., 2015, 2016). In order to facilitate a fair comparison to observational outliers, we first evaluate how well the simulations reproduce the empirical $M_{\rm BH}(\sigma)$ and Faber-Jackson relations (see Figures 5 and 6). Since these relations could depend sensitively on how either property is measured, observed outliers should only be compared to their respective empirical relation, and likewise, simulated BMGs to their respective relation from the Hydrangea sample.

In the Hydrangea simulation, the stellar velocity dispersion is defined as $\sigma=\sqrt[]{2E_{kin,\star}/3M_{\star}}$ where $E_{kin,\star}$ is the kinetic energy of stars, and $M_{\star}$ is the total mass of the stellar particles. Both $E_{kin,\star}$ and $M_{\star}$ were summed over all the stellar particles that are part of the subhalo, as identified by the SUBFIND algorithm. This one dimensional velocity dispersion implicitly assumes an isotropic stellar velocity dispersion. From Figure 5, it is immediately apparent that the velocity dispersion of BMGs A, B and C are remarkably low. Evaluating their $\sigma$ values over the course of the simulation showed that these galaxies had an order of magnitude higher velocity dispersion in the past. For BMG A and B, this drop in velocity dispersion occurs simultaneously to the tidal stripping events described in subsection 4.2. This may indicate that the high-velocity stars have been stripped from A and B. As discussed in section 4.2, these two galaxies are most likely still in the process of being tidally stripped, which implies that they are not in equilibrium at $z=0$. Furthermore, if the velocity dispersions of the BMGs were to be measured in a more observationally consistent way, their values could increase. This is especially relevant for BMGs A, B and C, since these galaxies reside in a dense environment ($<60$ pkpc from their host), and an observational aperture measurement will include more high-velocity background particles. To test the latter hypothesis, the velocity dispersions were re-computed within a spherical aperture of 5 kpc (including all star particles irrespective of their subhalo membership). With this, the velocity dispersions are indeed higher at $z=0$ (43 kms*-1*, 95 kms*-1* and 100 kms*-1* for A, B and C respectively). As can be seen in Figure 5, when using the re-computed stellar velocity dispersions, BMGs A, B and C lie closer to the median of the Hydrangea sample. For BMGs D and E, recomputing the velocity dispersions within a 5 kpc aperture made no significant difference. A velocity dispersion measurement using a cylindrical aperture was also tested. This increased the velocity dispersions by another $\sim 50\%$. This implies that the velocity dispersions of NGC 1271 and NGC 1277, which both reside in a dense cluster environment, could also be biased high due to the observations being unable to distinguish between host and satellite.

In Figure 5, it can be seen that the median relation from the Hydrangea sample agrees with the empirical $M_{\rm BH}(\sigma)$ relation from McConnell & Ma (2013) (the dash-dotted orange and solid red lines respectively). While NGC 1271 and Mrk 1216 lie along the empirical relation from McConnell & Ma (2013), NGC 1277, NGC1600 and M60 UCD1 are outliers above the $M_{\rm BH}(\sigma)$ relation. For M60 UCD1, the distance above the $M_{\rm BH}(\sigma)$ relation is comparable to the number of dex that BMGs B and C lie above the median relation of Hydrangea. However, we caution that the velocity dispersions of A, B and C lie outside of the observed range of velocity dispersions from McConnell & Ma (2013), thereby ruling out definitive comparisons. NGC 1277 and NGC 1600 respectively lie 0.66 and 0.98 dex above the $M_{\rm BH}(\sigma)$ relation from McConnell & Ma (2013), which is comparable to the offsets of BMGs D and E from the median relation in Hydrangea (1.15 and 1.02 dex, respectively).

The luminosities needed for the computation of the Faber-Jackson relation of the Hydrangea sample require stellar population synthesis modelling of the stellar particles, which is beyond the scope of this work. However, to first order, luminosity is proportional to $M_{\star}$. In Figure 6, we therefore compare the BMGs to the observed $\sigma(M_{\star})$ relation for early type galaxies (from Cappellari et al., 2013). The $\sigma(M_{\star})$ relation measured by Cappellari et al. (2013) agrees well with the median value of $\sigma(M_{\star})$ for the whole observed $M_{\star}$ range (the solid red and dash-dotted orange lines in Figure 6). Simulated BMGs A, B and C lie outside the domain of this relation. The observed outliers NGC 1277 and Mrk 1216 are positive outliers ($0.27$ dex each) to the relation from Cappellari et al. (2013). Assuming that these observed $\sigma$ values are correct, this indicates that they either have a DM halo which is much more massive than typical for their $M_{\star}$, or that they are out of equilibrium. I.e., they have very recently been stripped dramatically, but the (remaining) stars have not yet responded to the shallower potential well. The former explanation seems very unlikely for NGC 1277, since we expect satellites to have less DM at fixed $M_{\star}$. The latter explanation, is in line with what we found for BMG E (see section 4.1). Simulated BMGs D and E have a small offset from the median of the Hydrangea sample (0.02 and 0.12 dex above the median respectively), which is comparable to the offset of NGC 1600 and NGC 1271 to the empirical relation (0.04 and 0.17 dex respectively). BMGs D and E are very compact galaxies (with $R_{e}=3.93$ kpc and $R_{e}=2.61$ kpc). It could be that the stellar velocity dispersions of very compact galaxies are underestimated due to numerical resolution effects, since at lower resolution, the softening artificially smooths out the density in the centre, while at higher resolution, the central potential well should become deeper, and so the velocity dispersion of the galaxy, in the very centre, might increase.

Figures 5 and 6 show that the $M_{\rm BH}(\sigma)$ and $\sigma(M_{\star})$ relations from Hydrangea are comparable to empirical results within the appropriate $\sigma$ and $M_{\star}$ ranges. Hence, BMGs D and E could in particular be relevant for the understanding of the origins of observed compact galaxies with high $M_{\rm BH}/M_{\star}$. In terms of $M_{\rm BH}$ and $M_{\star}$, BMGs D and E are comparable to the observed outliers NGC 1271, (Walsh et al., 2015), NGC 1277 (van den Bosch et al., 2012; Walsh et al., 2016) and Mrk 1216 (Yıldırım et al., 2015; Walsh et al., 2017). There are no outliers with $M_{\rm BH}>10^{10}\rm M_{\odot}$, which is similar to the inferred $M_{\rm BH}$ of NGC 1600 ($M_{\rm BH}=1.7\pm 0.15\times 10^{10}\rm M_{\sun}$, Thomas et al., 2016). This can be explained by the general lack of galaxies in this mass range (see Figure 1). A simulation with even better statistics at the high-mass end would be required to investigate such systems.

We will now compare the BMGs D and E to observed outliers with similar $M_{\rm BH}$ and $M_{\star}$.

The BMGs D and E have one-dimensional stellar velocity dispersions $\sigma$ = 163 and 213 km s*-1*, respectively, while the observed outliers NGC 1271, NGC1277 and Mrk 1216 have $\sigma$ = $276^{+73}_{-4}$, $333\pm 12$ and $368\pm 3$ km s*-1*, respectively (see Table 1). The quoted stellar velocity dispersions of NGC 1271 and 1277 specifically exclude the sphere of influence of the black hole. However, this is not the case for Mrk 1216, where $\sigma$ is defined as the average stellar velocity dispersion within the effective radius.

As discussed above, it could be that the stellar velocity dispersion is underestimated for galaxies as compact as our BMGs, due to numerical resolution effects. However, the simulated $\sigma$ values of BMGs D and E lie within the spread of the velocity dispersion at fixed $M_{\star}\approx 10^{11}\rm M_{\odot}$ from the ATLAS3D sample (Cappellari et al., 2013). According to Cappellari et al. (2013), M $\propto\sigma^{2}\times R_{e}$, where $R_{e}$ is the projected half-light radius. This implies that at a constant $M_{\star}$, the velocity dispersion depends mainly on the size of the galaxy. The observed values of $R_{e}$ for NGC 1271 ($R_{e}$ = $2.2\pm 0.9$ kpc) and NGC 1277 ($R_{e}$ = $1.2\pm 0.1$ kpc), are two to three times lower than the half-mass radii of BMGs D ($R_{e}$ = 3.93 kpc) and E ($R_{e}$ = 2.61 kpc). This difference in size may account for the lower velocity dispersions in our simulated BMGs relative to these observed examples, despite having very similar $M_{\star}$.

The ages of both NGC 1277 and NGC 1271 are estimated to be greater than 10 Gyr (Ferré-Mateu et al., 2015), which is consistent with the ages of all of our simulated BCGs. The age of Mrk 1216 is even estimated to be >12 Gyr (Ferré-Mateu et al., 2017). This is comparable to BMG D, which is primarily an outlier due to its early formation time.

Another observed characteristic of compact outliers is their fast-rotator structure (Yıldırım et al., 2015; Walsh et al., 2015, 2016; Ferré-Mateu et al., 2017). From the analysis of Lagos et al. (2018), BMG E is a slow rotator ($\lambda_{R}$ close to 0), which adds to the potential differences in detail discussed above.

All of the BMGs in the simulation reside in a very rich cluster environment. It is thus not surprising that most outliers have been heavily stripped during their lifetimes. NGC 1271 and NGC 1277 are both embedded in the nearby Perseus galaxy cluster and are relatively compact galaxies (effective radius smaller than 2.5 kpc Walsh et al., 2016, 2015; Trujillo et al., 2014), however they do not show any obvious signs of having been stripped. Mrk 1216 resides in a significantly different environment; it has only two neighbouring galaxies within 1 Mpc (Yıldırım et al., 2015). This galaxy is also relatively compact, though deep field imaging again show no apparent signatures of interaction or stripping (Ferré-Mateu et al., 2017).

As discussed above, BMGs D and E share characteristics with the observed outliers NGC 1271, NGC 1277 and Mrk 1216, but also differ in some details such as their offsets from the $M_{\rm BH}$($\sigma$) and $\sigma$($M_{\star}$) relations or kinematic structure (fast vs. slow rotator). It is unclear whether these differences are related to their BMG nature, or reflect unrelated shortcomings of the simulation model. Further, detailed work would be required to investigate this issue further (see also van de Sande et al., 2018 and Lagos et al., 2018).

The similarities described above, between NGC 1271, NGC 1277, Mrk 1216 and the simulated BMGs D and E could hint towards a similar origin: relics of the $z>2$ Universe. It is important to note that these ‘high-z relics’ may also have been significantly stripped, as is the case for both BMG C and E, which have high values for both $f_{\rm strip}$ and $f_{\rm age}$. The fact that the nearest more massive neighbour of BMG E lies at > 1 Mpc (which is the case for Mrk 1216) does not necessarily imply that they have not undergone tidal stripping at earlier times, as they may have drifted away from these more massive neighbours. NGC 1271, NGC 1277 and Mrk 1216 may also have gone through a similar tidal stripping period but ‘survived’ to eventually end up seemingly undisturbed at $z=0$. It is therefore interesting to study the environments of outliers out to scales of > 1 Mpc, and we encourage observational efforts to do so.

4.4 Alpha enhancement

The [Mg/Fe] ratio is an indicator of how extended the SF period of a galaxy has been, since core-collapse supernovae (SNe), which take place shortly after a star formation episode, release $\alpha$-elements (such as the easily observable Mg), while Type Ia SN, which occur mostly after $\sim 1$ Gyr, release primarily iron (Fe) and no $\alpha$-elements. The more abruptly SF is quenched, the less time there is to form Fe via Type Ia SNe and to incorporate it into stars. Segers et al. (2016b) have shown that in EAGLE, AGN feedback is responsible for its success in reproducing the observed $\alpha$-enhancement of massive galaxies. We thus expect a positive correlation between [Mg/Fe] and $M_{\rm BH}$ at fixed $M_{\star}$. Recently, Martín-Navarro et al. (2016) observed that positive galaxy outliers of the $M_{\rm BH}-\sigma$ relation are older and more $\alpha$-enriched. Because of the similarity between the observed $M_{\rm BH}-\sigma$ and the $M_{\rm BH}-M_{\star}$ relations, we have examined these relations in the Hydrangea simulations.

To examine the secondary trends in [Mg/Fe] and $M_{\rm BH}$, it is important that the simulated Mg/Fe relation matches observations. Segers et al. (2016a) and Segers et al. (2016b) showed that the EAGLE simulated trends between [$\alpha$/Fe], galaxy stellar mass and age are in good agreement with observations of early-type galaxies for galaxies with log10($M_{\star}$/$\rm M_{\odot}$)> 10.5. Furthermore, when we examine the distribution of simulated galaxies with $\sigma$ > 100 km s*-1* and sSFR > 0.1 $\rm M_{\odot}$/Gyr (in order to facilitate comparison with the observational data) it appears that BMGs D and E have relatively high [Mg/Fe] values at fixed $\sigma$, which is in qualitative agreement to the observations. This, in combination with the overlap between the $M_{\rm BH}(\sigma)$ relation in Hydrangea, and the empirical relation of McConnell & Ma (2013) for log10($\sigma$/kms*-1*) > 2.2 (see Figure 5), argue for a meaningful comparison of secondary trends in [Mg/Fe] and $M_{\rm BH}$ between the Hydrangea simulation and observations.

In Figure 7 we plot the residuals relative to the median [Mg/Fe]$(\sigma)$ relation as a function of the residuals of the median $M_{\rm BH}(\sigma)$ relation for galaxies in the Hydrangea simulations at $z=0$ with $\sigma>100$ km s*-1*, sSFR > 0.1 $\rm M_{\odot}$/Gyr and $M_{\rm BH}>10^{7}\rm M_{\odot}$, where the first two constraints facilitate comparison to the observational data and the latter avoids resolution effects. In this way, we can investigate the relation between $M_{\rm BH}$ and the SF history at fixed $\sigma$. The observed relation from Martín-Navarro et al. (2016) is shown in red, with the dashed lines showing the uncertainty in the fit.

It can be seen that the scatter in $M_{\rm BH}(\sigma)$ correlates positively with the scatter in Mg/Fe, in qualitative agreement with the findings of Martín-Navarro et al. (2016). The difference in functional form may be due to selection effects, which are also reflected in the sample sizes; Martín-Navarro et al. (2016) uses a very heterogeneous sample of 57 galaxies, while we include $\sim 3400$ galaxies. The vast majority of the galaxy sample from Martín-Navarro et al. (2016) resides in the domain -0.5 < $\delta\rm log_{10}(M_{\rm BH}(\sigma))$ < 0.5, which is also the domain in which the median relations are in best agreement.

Although only 2 out of 5 BMGs have $\sigma>100$ km s*-1* at $z=0$, they appear to be in line with the positive correlation between the scatter in $M_{\rm BH}(\sigma)$ and the scatter in Mg/Fe. The [Mg/Fe] abundance of outlier D ([Mg/Fe] = 0.36) and E ([Mg/Fe] = 0.31) are similar to the best-fitting value from Martín-Navarro et al. (2016) for NGC 1277 ([Mg/Fe] = 0.32) and Mrk 1216 (0.26 $\pm$ 0.05 Ferré-Mateu et al. 2017). Note that uncertainties in stellar population modelling lead to systematic errors in observed alpha-enhancement values of $\approx 0.1$ dex, while nucleosynthetic yields (used in the simulations) are uncertain by a factor $\approx 2$ (Wiersma et al., 2009b). Thus, BMGs D and E have [Mg/Fe] consistent with the observed outliers. This points towards a similar SF history.

5 Summary

We have used the Hydrangea simulations (Bahé et al., 2017), part of the C-EAGLE project (Barnes et al., 2017), of massive galaxy clusters and their surroundings to study the evolution and origin of galaxies with central SMBHs that are overmassive relative to the general $M_{\rm BH}-M_{\star}$ relation at $z=0$ in galaxy cluster environments. This work builds on B16, who used the (100 Mpc)3 EAGLE simulation to study the origin of overmassive black holes with smaller masses and in less overdense environments than studied here. We probe galaxies in a mass regime comparable to observed $M_{\rm BH}(M_{\star})$ outliers. We define $M_{\rm BH}(M_{\star})$ outliers within the simulation as those with $M_{\rm BH}/M_{\star}>0.01$ and $M_{\rm BH}>10^{7}\rm M_{\odot}$.

From this initial set of $M_{\rm BH}(M_{\star})$ outliers, a selection by eye of the five most extreme outliers of this relation, which we term “Black hole monster galaxies” (BMGs), has been made based on their location in the $M_{\rm BH}-M_{\star}$ plane (Figure 1). They were chosen such that they have masses similar to observed galaxies that are extreme positive outliers of the observed $M_{\rm BH}-M_{\star}$ relation. We have traced the BMGs through time to study their formation and evolution. Our main findings are as follows:

• •

We identified 69 galaxies that fulfil our $M_{\rm BH}(M_{\star})$ outlier criteria at $z=0$, of which 77$\%$ are satellites. The outliers lie closer to their host and are significantly more compact than other satellites of similar $M_{\star}$, suggesting that tidal stripping is an important contributor toward galaxies becoming outliers. Outliers typically have $M_{\star}=[10^{9}-10^{11}]\rm M_{\odot}$, where our most massive BMG (ID E, see Table 1), reaches a stellar mass of $1.31\times 10^{11}\rm M_{\odot}$ and $M_{\rm BH}=4.05\times 10^{9}\rm M_{\odot}$ (see Figure 1). The simulations therefore reproduce the observed existence of $M_{\rm BH}(M_{\star})$ outliers with $M_{\rm BH}>10^{9}\rm M_{\odot}$.

• •

In agreement with B16, our results suggest that the origin of galaxies with such over-massive BHs is primarily due to a combination of tidal stripping (Figure 3, see also Volonteri et al., 2008, 2016) and their formation at high redshift (Figure 4). For the BMGs with $M_{\rm BH}>10^{9}\rm M_{\odot}$, galaxy-galaxy merger events also plausibly play a role in the formation of their anomalously high-mass central BHs.

• •

Galaxies with higher $M_{\rm BH}$ tend to have a higher [Mg/Fe] ratio at fixed stellar velocity dispersion, in agreement with the observational results of Martín-Navarro et al. (2016) (Figure 7). This correlation is thought to underline the strong influence of AGN feedback on the star formation history of high-mass galaxies (Segers et al., 2016b).

• •

BMGs D and E have similar $M_{\rm BH}$ and $M_{\star}$ as observed $M_{\rm BH}(M_{\star})$ outlier galaxies NGC 1271 (Walsh et al., 2015), NGC 1277 (Walsh et al., 2016) and Mrk 1216 Walsh et al. (2017). Both D and E host a BH that entered the rapid growth phase at $z>2$ and both exhibit high [Mg/Fe] values, consistent with being “relics” of the high-$z$ ($z>2$) Hydrangea universe. The properties of these BMGs are summarized in Table 1 (see the discussion in Section 4.2). BMG E is a ‘survivor’ of tidal stripping events (see Figure 8). Its outer stellar halo was completely disrupted by tidal interactions with a proto-cluster at $z\approx 0.5$, and the remaining core does not appear particularly disturbed at $z=0$. We therefore encourage observational efforts to search for outliers in clusters even out to large ($\sim$1 Mpc) distances, which corresponds to the distance of outlier E to its nearest more-massive neighbour at $z=0$. Due to the similarities between the observed outliers and the simulated BMGs, a similar origin is plausible. This means that in addition to being a high-$z$ relic, outliers could very well have been stripped at an earlier time, survived, and appear undisturbed at $z=0$. On top of this, all BMGs were found to have experienced a merger event immediately preceding the rapid growth phase of their central BH. It is therefore plausible that major merger events triggering the rapid growth phase of BHs could contribute to the formation of anomalously massive SMBHs

Based on hydrodynamical, cosmological simulations we thus expect the existence of observed $M_{\rm BH}(M_{\star})$ outliers with $M_{\rm BH}>10^{9}\rm M_{\odot}$. Observed outliers with such large $M_{\rm BH}$ are plausibly survivors of tidal stripping and it would be interesting to study their environment out to $\sim$1 Mpc. The positive correlation between the scatter in $M_{\rm BH}(\sigma)$ and in Mg/Fe indicates a strong relation between the $M_{\rm BH}$ and the time-scale of the quenching process for massive galaxies. It would furthermore be interesting for future studies to investigate the relation between merger-driven rapid growth phases at high redshifts ($z>2$) and anomalously high-mass SMBHs.

Acknowledgements

This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility. DiRAC is part of the National E-Infrastructure. This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. STK and DJB acknowledge support from STFC through grant ST/L000768/1. RAC is a Royal Society University Research Fellow. The Hydrangea simulations were in part performed on the German federal maximum performance computer ‘Hazel Hen’ at the maximum performance computing centre Stuttgart (HLRS), under project GCS-HYDA/ID44067 financed through the large-scale project ‘Hydrangea’ of the Gauss Centre for Supercomputing. Further simulations were performed at the Max Planck Computing and Data Facility in Garching, Germany. This project received funding from the EU Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement 747645 (ClusterGal). This work has furthermore benefited from the use of Py-SPHViewer (Benitez-Llambay 2015).

Appendix A BMG properties

In Table 1 we list all of the relevant properties of the 5 BMGs and observed outliers. Additionally, we list information on some of the most notable observed $M_{\rm BH}(M_{\star})$ outliers. For NGC 1271 and NGC1277, the [Mg/Fe] values are the best-fit values derived by Martín-Navarro et al. (2016). In column 7, the ‘outlier $z$’ represents the lower limits of a 500 Myr time step.

In Figure 8 we have plotted the stellar mass density of outlier ID E and its surroundings at three different snapshots throughout the simulation. The remaining images are available at https://sites.google.com/view/vson-umbh-galaxies/homepage