Galaxies with monstrous black holes in galaxy cluster environments

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.


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 (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 E-mail: vson@strw.leidenuniv.nl been observed, the most interesting being those with extremely high black hole mass (MBH) for their stellar mass. Well-studied examples are NGC 1271 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,  and Graham et al. (2016) both argue that when spheroid/disc decompositions are modelled under different assumptions, the MBH values of e.g. NGC 1277 and NGC 1271 can be brought to within 2-σ of the observed MBH(M ) relation, and are no longer anomalous. It is therefore interesting to test whether such MBH(M ) 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 MBH(M ) 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 MBH > 10 8.5 M although, as mentioned above, such BH outliers have been observed. Plausibly, the lack of highmass 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 ∼ 10 −7 galaxies per 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 MBH(M ) outliers with such high MBH. 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 MBH(M ) outliers with MBH of > 10 8.5 M ) 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 MBH(M ) 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 MBH > 10 8.5 M are found in Hydrangea, it is interesting to see whether the origin of their anomalous MBH(M ) 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 ,1/M ,2 = µ > 0.25, while intermediate mergers are 0.1 < µ < 0.25, where M ,1 and M ,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 MBH > 10 9 M ) 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 MBH(M ) 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 MBH(M ) 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.

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 from 1 M200 = 10 14 to ≈ 2 × 10 15 M . 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 10r200 from each cluster's centre. Particles in the re-simulated higher-resolution regions initially have masses identical to those in the intermediateresolution EAGLE simulations: 1.81 × 10 6 M for gas particles and 9.7 × 10 6 M 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 EA-GLE simulation. Furthermore, between each of the 28 main snapshots, three 'snipshots' were stored, reducing the time resolution to t = 125Myr. 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 = 25Myr 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 metallicitydependent 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 M 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 M h −1 . MBH 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, Cvisc, that influences the MBH(M ) relation (Crain et al. 2015), and was calibrated to reproduce the galaxy stellar mass function (GSMF). To model (unresolved) dynamical friction effects, BHs with MBH < 100 mgas 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 ∆TAGN (= 10 9 K instead of 10 8.5 in EAGLE) and Cvisc(= 2π × 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 r fṁacc,BH c 2 , where r = 0.1 is the radiative efficiency of the accretion disc and 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 MBH(M ) relation at z ≈ 0.

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 par-ticle (if this DM particle belongs to any halo). The SUB-FIND 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 , as the sum over all stellar particles bound to the subhalo. 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.

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 selfbound, 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 MBH 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 ≈ 0.5 − 1% of the BH particles per snapshot, but had no effect on the final sample of BMGgalaxies.

BMG selection
In order to define (positive) MBH(M ) outliers as galaxies with a central SMBH that is much more massive (or equivalently have a stellar mass that is much lower) than expected 2 When using M ,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 BH outliers studied in this work, M ,30 and M are equal at z = 0. The maximum difference encountered was log 10 M /M -log 10 M ,30 /M 10 −2 .
relative to the median MBH − M relation at z = 0, we select galaxies with MBH/M > 0.01 and MBH > 10 7 M . 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 ) BHs. This definition differs from that of B16, who defined outliers as having MBH more than 1.5 dex above the median MBH(M ) relation in the simulation. When using the definition from B16, the fraction of all galaxies with MBH > 10 7 M 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 MBH > 10 7 M 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 MBH (> 10 9 M ).  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.

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 , MBH and MDM 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 A1 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 par-ticles 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 ∼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 ∼ 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. The median relation, shown as a solid red line, is lower for satellite galaxies than the median relation for centrals at M ≈ 10 12 M . This is caused by a small subset (∼ 20) of 'undermassive' BH galaxies (with MBH ≈ 10 8.5 M ) 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 MBH(M ) for both centrals and satellites.
The 5 manually selected BMGs are labelled in Figure  1 with IDs corresponding to those in table A1. They were selected with a preference for a high MBH because of the similarity to the extremely high MBH of several observed MBH(M ) outliers. A bias towards higher MBH at fixed M also leads to a bias towards older galaxies; as the gradient in stellar age in Figure 1 shows, at fixed M , galaxies with higher MBH tend to be older. There is a larger variation in M at fixed MBH for satellite galaxies, with respect to that for centrals at M ∼ 10 10 M . 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 MBH and 3 Videos 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 Figure 1. Galaxy black hole mass as a function of stellar mass for Hydrangea galaxies at z = 0. Circular markers denote satellite galaxies (lower middle panel), whereas triangles correspond to central galaxies (lower left panel). Colours (ranging from blue to yellow) indicate the mean mass-weighted age of the stars in Gyr. The limit above which galaxies are defined as M BH (M ) outliers (M BH /M > 0.01 and M BH > 10 7 M ) is shown as a dashed red line. The median M BH in bins of M is shown as a dash-dotted orange line. We also indicate the 5 'Black hole Monster Galaxies' (BMGs) labelled with their BMG-ID (see Table A1). Galaxies from the (100 Mpc M 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 , and therefore also of supermassive BHs. Because of this larger number of massive galaxies, they are also more likely to contain MBH(M ) outliers at this high-mass end. Outliers with the highest MBH (outliers D and E) have masses comparable to the most recent estimates for NGC 1277 and NGC 1271 (Walsh et al. 2016. For comparison, the outliers identified by B16 in the (100 Mpc) 3 EAGLE simulations that satisfy our requirements of being a MBH(M ) 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 MBH(M ) relation in the simulation at z = 0, which differs from our threshold (see Section 2.2.2).
The distribution of galaxies in the MBH − M plane can be divided into two general regions. As described by Bower et al. (2017), above a halo mass of 10 12 M , stellar feedback is no longer sufficient to prevent gas infall towards the central BH once M reaches ∼ 10 10.5 M . The BH subsequently begins a phase of runaway growth dur-ing 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 ≈ 10 10.5 M . 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 and MBH is increasingly due to mergers.
As an extension to the work done by B16, we find from the MBH−M distribution in Figure 1 that strong MBH(M ) outliers with MBH greater than 10 9 M are expected to exist in cluster environments. These high values are comparable to the most massive observed MBH(M ) outliers. In the next Section we will study the origin of these BMGs, and compare this to results from B16.

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 MBH − M plane (see Figure 1). These 5 outliers, or BMGs, were selected based on their extreme BH masses and similarity to observed MBH 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 MBH − M init plane, where M ,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.

Tidal stripping and early formation
We would like to know whether the physical mechanisms by which these galaxies become such extreme outliers of the MBH(M ) relation are consistent with those found by B16. We therefore investigate the relative contributions of tidal stripping and early formation time to their MBH/M 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 . 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 > 10 9 M , as shown by a KS-test with p 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 ,init,max) that From this it appears that stripping is a significant contributing factor in the formation of BMGs. Even so, galaxies C, D and E remain strong outliers when stripping effects are removed.
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 ,init,max and the M ,init (at z = 0) a good measure of stellar stripping. Following B16, we quantify the impact of stripping by the logarithmic ratio of M ,init,max and M ,init at z = 0: The values of fstrip are shown in Table A1 for all of our BMGs and are also visualized by the coloured lines in Figure 3. The BMGs with the largest fstrip (IDs B and C in Table A1) have lost an order of magnitude in M . 4 The horizontal movement along the M 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 ,init).
B16 found that a secondary factor causing galaxies to become outliers of the MBH(M ) relation is their stellar formation time. At earlier times, the normalization of the MBH(M ) relation was higher. This difference in normalization originates from the difference in the density of the universe which causes galaxies at fixed M 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ṀBH ∝ MBH. However, in practice gas densities are low enough that BH accretion proceeds according to the Bondi-Hoyle rate, which scales aṡ MBH ∝ M 2 BH , which causes a highly non-linear accretion phase. The final MBH 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 ≈ 6 for MBH 10 9 M ). 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 MBH/M ,init ratio against the mean stellar ages (weighted by their initial stellar mass) of central galaxies with MBH > 10 7 M . 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 MBH > 10 7 M , 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 ,init to exclude stripping effects. The median of the MBH/M ,init − age relation (black dashed line) is a measure of the contribution that the age of a galaxy has on the scatter in the MBH − M relation. Here we see a clear positive trend of MBH/M ratio with mean stellar age. This trend represents the evolution of the normalization of the MBH(M ) relation over time.
As in B16, we quantify the contribution of age in a similar manner to that of the tidal stripping contribution.
For a given BMG, fage is the logarithmic difference between the value of the median MBH/M ,init,max -age relation measured at the stellar age of the galaxy (the black dashed line in Figure 4) and α med = -3.15, the median MBH/M ,init value measured for the whole population of central galaxies with MBH > 10 7 M (the dashed red line in Figure 4). The values of fage are shown in Table A1 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 MBH(M ) outliers with MBH > 10 9 M . However, these BMGs still lie on the outer edges of the distribution (at > 1σ 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. ratio as a function of initial stellar mass-weighted ages for central galaxies with M BH > 10 7 M at z = 0. The dash-dotted orange line shows the median relation, while the solid black line represents the median M BH /M ,init value of this entire selection of galaxies. The BMGs are plotted at the maximum M ,init they ever had, in order to remove the effects of tidal stripping. The tops of the coloured, vertical lines correspond to the value of the median M BH /M ,init at the stellar ages of the BMGs. The light and dark gray shaded areas represent the 1 and 2-σ standard deviations from the median relation, respectively. It can be seen that, in particular, outliers C, D and E were formed at relatively early times.

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  A1. 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 MBH/M > 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 fstrip values (0.76 and 1.49 respectively) bot confirm that BMGs A and B are remnants of tidal stripping.
The BMGs with highest MBH are of particular interest as they have similar MBH to observed outliers of interest, since the extraordinarily high MBH/M ratios of these observed galaxies are still debated 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 MBH/M 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 fstrip = 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 MBH > 10 7 M by z = 4. It enters a second rapid growth phase after another major merger event at z ≈ 2.5, causing MBH to increase to ≈ 2 × 10 9 M , 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 MBH/M > 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 . By z = 3 it has reached MBH ≈ 10 8.9 M . It continues to grow in MBH via gas accretion and a second major merger event at z = 2, reaching its maximum MBH of 10 9.0 M by z ≈ 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 MBH 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 MBH = 10 9.45 M . The corresponding AGN feedback blast has cleared all the gas from the galaxy by z = 2.0, and the galaxy passes the MBH/M > 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 MBH ≈ 10 7.4 M and MBH = 10 7.9 M respectively. The latter two mergers increase the initial stellar mass to M = 10 11.5 M and M = 10 11.7 M respectively. These merging events are visible in Figure 2 as the 'staircase' shape that BMG E (solid red line) follows at high MBH.
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 ≈ 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 A1). 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 A1). 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 2; see Table A1). 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 MBH 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 (LAGN > 10 45 erg/s) at z ≈ 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.

Comparison to observations
It is interesting to know whether the origins of the simulated BMGs could be used to explain the origins of observed MBH(M ) outliers. Compact positive outliers from the MBH − M relation are observationally characterized by following the MBH(σ) relation, and by being offset from the Faber-Jackson relation (Walsh et al. , 2016. In order to facilitate a fair comparison to observational outliers, we first evaluate how well the simulations reproduce the empirical MBH(σ) 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 σ = 2E kin, /3M where E kin, is the kinetic energy of stars, and M is the total mass of the stellar particles. Both E kin, and M 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 σ 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 ∼ 50%.
In Figure 5, it can be seen that the median relation from the Hydrangea sample agrees with the empirical MBH(σ) relation from Saglia et al. (2016) (the dash-dotted orange and solid red lines respectively). While NGC 1271 and Mrk 1216 lie within the error of the empirical relation from Saglia et al. (2016), NGC 1277, NGC1600 and M60 UCD1 are outliers above the MBH(σ) relation. For M60 UCD1, the distance above the MBH(σ) 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 Saglia et al. (2016), thereby ruling out definitive comparisons. Furthermore, the median was computer for the complete sample of satellites and centrals with log10(MBH/M ) > 7 at z = 0, while McGee (2013) found that central galaxies have a significantly steeper slope. Splitting the median relation for centrals and satellites would reduce the vertical offset between the median and the BMGs. NGC 1277 and NGC 1600 respectively lie 0.66 and 0.98 dex above the MBH(σ) relation from Saglia et al. (2016), 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 . In Figure 6, we therefore compare the BMGs to the observed σ(M ) relation for early type galaxies (from Cappellari et al. 2013). The σ(M ) relation measured by Cappellari et al. (2013) agrees well with the median value of σ(M ) for the whole observed M 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 σ values are correct, this indicates that they either have a DM halo which is much more massive than typical for their M , 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 . 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 Re = 3.93 kpc and Re = 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 MBH(σ) and σ(M ) relations from Hydrangea are comparable to empirical results within the appropriate σ and M ranges. Hence, BMGs D and E could in particular be relevant for the understanding of the origins of observed compact galaxies with high MBH/M . In terms of MBH and M , BMGs D and E are comparable to the observed outliers NGC 1271, ), NGC 1277(van den Bosch et al. 2012Walsh et al. 2016) and Mrk 1216 (Yıldırım et al. 2015;Walsh et al. 2017). There are no outliers with MBH > 10 10 M , which is similar to the inferred MBH of NGC 1600 (MBH = 1.7 ± 0.15 × 10 10 M , 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 MBH and M .
The BMGs D and E have one-dimensional stellar velocity dispersions σ = 163 and 213 km s −1 , respectively, while the observed outliers NGC 1271, NGC1277 and Mrk 1216 have σ = 276 +73 −4 , 333 ± 12 and 368 ± 3 km s −1 , respectively (see Table A1). 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 σ 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 σ values of BMGs D and E lie within the spread of the velocity dispersion at fixed M ≈ 10 11 M from the ATLAS 3D sample (Cappellari et al. 2013). According to Cappellari et al. (2013), M ∝ σ 2 ×Re, where Re is the projected half-light radius. This implies that at a constant M , the velocity dispersion depends mainly on the size of the galaxy. The observed values of Re for NGC 1271 (Re = 2.2±0.9 kpc) and NGC 1277 (Re = 1.2±0.1 kpc), are two to three times lower than the half-mass radii of BMGs D (Re = 3.93 kpc) and E (Re = 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 . 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. 2015Walsh et al. , 2016Ferré-Mateu et al. 2017). From the analysis of Lagos et al. (2018), BMG E is a slow rotator (λ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. 2016Walsh et al. , 2015Trujillo 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 MBH(σ) and σ(M ) 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 andLagos 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 fstrip and fage. 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. Although the high central density of NGC 1277 argues for a high formation redhift, we stress that this does not exclude the presence of tidal stripping. This attenuates the results from Beasley et al. (2018), who argue NGC 1277 is 'a genuine relic' galaxy that has not been significantly influenced by tidal stripping. It is therefore interesting to study the environments of outliers out to scales of > 1 Mpc, and we encourage observational efforts to do so.

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 α-elements (such as the easily observable Mg), while Type Ia SN, which occur mostly after ∼ 1 Gyr, release primarily iron (Fe) and no α-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 α-enhancement of massive galaxies. We thus expect a positive correlation between [Mg/Fe] and MBH at fixed M . Recently, Martín-Navarro et al. (2016) observed that positive galaxy outliers of the MBH − σ relation are older and more α-enriched. Because of the similarity between the observed MBH − σ and the MBH − M relations, we have examined these relations in the Hydrangea simulations.
To examine the secondary trends in [Mg/Fe] and MBH, 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 [α/Fe], galaxy stellar mass and age are in good agreement with observations of early-type galaxies for galaxies with log10(M /M )> 10.5. Furthermore, when we examine the distribution of simulated galaxies with σ > 100 km s −1 and sSFR > 0.1 M /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 σ, which is in qualitative agreement to the observations. This, in combination with the overlap between the MBH(σ) relation in Hydrangea, and the empirical relation of Saglia et al. (2016) for log10(σ/kms −1 ) > 2.2 (see Figure 5), argue for a meaningful comparison of secondary trends in [Mg/Fe] and MBH between the Hydrangea simulation and observations. In Figure 7 we plot the residuals relative to the median [Mg/Fe](σ) relation as a function of the residuals of the median MBH(σ) relation for galaxies in the Hydrangea simulations at z = 0 with σ > 100 km s −1 , sSFR > 0.1 M /Gyr and MBH > 10 7 M , 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 MBH and the SF history at fixed σ. 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 ∼ 3400 galaxies. The vast majority of the galaxy sample from Martín-Navarro et al. (2016) resides in the domain -0.5 < δlog10(MBH(σ)) < 0.5, which is also the domain in which the median relations are in best agreement.
Although only 2 out of 5 BMGs have σ > 100 km s −1 at z = 0, they appear to be in line with the positive correlation between the scatter in MBH(σ)

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 MBH − M 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 MBH(M ) outliers. We define MBH(M ) outliers within the simulation as those with MBH/M > 0.01 and MBH > 10 7 M .
From this initial set of MBH(M ) 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 MBH −M plane ( Figure  1). They were chosen such that they have masses similar to observed galaxies that are extreme positive outliers of the observed MBH − M 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 MBH(M ) 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 , suggesting that tidal stripping is an important contributor toward galaxies becoming outliers. Outliers typically have M = [10 9 − 10 11 ]M , where our most massive BMG (ID E, see Table A1), reaches a stellar mass of 1.31 × 10 11 M and MBH = 4.05 × 10 9 M (see Figure 1). The simulations therefore reproduce the observed existence of MBH(M ) outliers with MBH > 10 9 M .
• 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. 2008Volonteri et al. , 2016 and their formation at high redshift ( Figure 4). For the BMGs with MBH > 10 9 M , galaxy-galaxy merger events also plausibly play a role in the formation of their anomalously high-mass central BHs.
• Galaxies with higher MBH 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 MBH and M as observed MBH(M ) outlier galaxies NGC 1271 , 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 A1 (see the discussion in Section 4.2). BMG E is a 'survivor' of tidal stripping events (see Figure A1). Its outer stellar halo was completely disrupted by tidal interactions with a proto-cluster at z ≈ 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 (∼1 Mpc) distances, which corresponds to the distance of outlier E to its nearest moremassive 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 highz 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 MBH(M ) outliers with MBH > 10 9 M . Observed outliers with such large MBH are plausibly survivors of tidal stripping and it would be interesting to study their environment out to ∼1 Mpc. The positive correlation between the scatter in MBH(σ) and in [Mg/Fe](σ) indicates a strong relation between the MBH 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. Figure A1. Stellar surface density distribution of BMG-ID E at three different snapshots. These images are created with the publicly available Py-SPHViewer, which uses the Smoothed Particle Hydrodynamics (SPH) interpolation scheme. The colours represent the projected density of stellar particles bound to BMG-ID E, with blue/pink corresponding to low density and yellow corresponding to high-mass density. Black hole particles are indicated with blue dots. The most relevant parameters of the galaxy at each snapshot are written on the top left of each image. A green scale bar shows the scale of the images in comoving kpc. The left and middle panels correspond to z = 0.68 and 0.57 respectively, just before and during the stellar stripping event. The right panel depicts the galaxy at z = 0. This BMG is located in a dense environment right after being stripped. However at z = 0 it has 'escaped' to a low density region, causing it to be the central galaxy of its FOF group. This portrays an example of an BMG surviving the tidal stripping process. Table A1. 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 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. The last two columns (10 and 11) show f strip and fage 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.  >10 0.23 NGC 1277 9.69 ± 0.14 11.08 ± 0.14 1.2 ± 0.1 333 ±12 >10 0.32 Mrk 1216 9.69 ± 0.15 11.30 ± 0.17 2.3 ± 0.1 368 ± 3 12.8 ± 1.5 0.26 ± 0.05