ABSTRACT

With projects such as Laser Interferometer Space Antenna (LISA) and Pulsar Timing Arrays (PTAs) expected to detect gravitational waves from supermassive black hole mergers in the near future, it is key that we understand what we expect those detections to be, and maximize what we can learn from them. To address this, we study the mergers of supermassive black holes in the Illustris simulation, the overall rate of mergers, and the correlation between merging black holes and their host galaxies. We find these mergers occur in typical galaxies along the MBHM* relation, and that between LISA and PTAs we expect to probe the full range of galaxy masses. As galaxy mergers can trigger star formation, we find that galaxies hosting low-mass black hole mergers tend to show a slight increase in star formation rates compared to a mass-matched sample. However, high-mass merger hosts have typical star formation rates, due to a combination of low gas fractions and powerful active galactic nucleus feedback. Although minor black hole mergers do not correlate with disturbed morphologies, major mergers (especially at high-masses) tend to show morphological evidence of recent galaxy mergers which survive for ∼500 Myr. This is on the same scale as the infall/hardening time of merging black holes, suggesting that electromagnetic follow-ups to gravitational wave signals may not be able to observe this correlation. We further find that incorporating a realistic time-scale delay for the black hole mergers could shift the merger distribution towards higher masses, decreasing the rate of LISA detections while increasing the rate of PTA detections.

1 INTRODUCTION

It is well established that supermassive black holes (SMBHs) are found at the centre of massive galaxies (Kormendy & Richstone 1995), and that the black hole (BH) masses correlate strongly with host galaxy properties, suggesting an evolutionary link between galaxy formation and BH growth (e.g. Magorrian et al. 1998; Gebhardt et al. 2000; Graham et al. 2001; Ferrarese 2002; Tremaine et al. 2002; Häring & Rix 2004; Gültekin et al. 2009; Kormendy & Ho 2013; McConnell & Ma 2013; Reines & Volonteri 2015; Greene et al. 2016; Schutte, Reines & Greene 2019). As such, a galaxy merger can provide the opportunity for SMBHs from the respective galaxies to migrate to the new galactic centre, become gravitationally bound, form a binary, and eventually merge (for a seminal discussion, see Begelman, Blandford & Rees 1980).

Over the past several years, gravitational wave (GW) signals from BH mergers have been detected by the LIGO–Virgo collaboration (Abbott et al. 2016), but so far these GWs have been limited to those produced by mergers between stellar mass BHs. The expected mergers between SMBHs at the centres of galaxies would produce much longer wavelength GWs, which the current ground-based interferometers are not sensitive to. However, the planned Laser Interferometer Space Antenna (LISA) space mission will be focused on lower frequency GWs, and is aimed at detecting mergers between SMBHs, with sensitivity peaking at ∼104−107 M (Amaro-Seoane et al. 2017). In addition, Pulsar Timing Array (PTA) observations should be capable of detecting mergers at even higher masses, reaching BHs above 108 M (Desvignes et al. 2016; Reardon et al. 2016; Verbiest et al. 2016; Arzoumanian et al. 2018). The expected detections of GWs from these observations should provide a completely novel and powerful framework to measure BHs and their co-evolution with their host galaxies.

GW detections of SMBH mergers can provide estimates for the rate at which SMBHs merge throughout the Universe (e.g. Klein et al. 2016; Salcido et al. 2016; Kelley et al. 2017b; Ricarte & Natarajan 2018; Katz et al. 2020), the role mergers play in setting the BH–galaxy scaling relations (e.g. Volonteri & Natarajan 2009; Shankar et al. 2016; Simon & Burke-Spolaor 2016), gas environment and accretion efficiencies of BHs (e.g. Kocsis, Yunes & Loeb 2011; Barausse, Cardoso & Pani 2014; Derdzinski et al. 2019), and even the mechanism by which BH seeds form (see e.g. Sesana, Volonteri & Haardt 2007; Ricarte & Natarajan 2018; DeGraf & Sijacki 2020). Beyond information about the BHs themselves, multimessenger studies that combine GW and electromagnetic observations have the potential to analyse the galaxies in which these mergers occur, directly linking supermassive BH mergers with galaxy properties of the host.

Cosmological hydrodynamical simulations that incorporate both galaxy formation and SMBHs provide an excellent tool to investigate the connection between BH mergers detectable by GWs and the properties and histories of the host galaxies in which they are found. Current cosmological simulations (e.g. Dubois et al. 2014; Vogelsberger et al. 2014a; Schaller et al. 2015; Feng et al. 2016; Pillepich et al. 2018a) are able to probe a wide range of spatial scales, with BHs ranging from ∼104−1010 M, combined with galaxies with resolved morphological structure (e.g. Snyder et al. 2015, 2019). Since these simulations self-consistently model the growth of BHs and the co-evolution of the host galaxies in large numbers, they provide an ideal resource to make predictions for what we can expect from upcoming GW detections, associated electromagnetic follow-ups, and how to interpret those observations from a theoretical framework.

In this paper, we use the Illustris simulation (Nelson et al. 2015) to investigate the connection between BH mergers and the galaxies in which they take place, with a particular emphasis on the galaxy morphology of BH merger hosts. Using the detailed BH data in the Illustris simulation, we are able to connect each BH merger to the galaxy in which it occurs, and thus generate a list of every BH merger each galaxy has hosted throughout its history. This allows us to analyse how recent BH mergers correlate with galaxy properties [such as morphology and star formation rate (SFR)], the time-scale over which these correlations survive, and the implications this has for electromagnetic follow-up observations to GW detections.

The outline of this paper is as follows. In Section 2, we discuss the Illustris simulations used for this project and the BH models therein. In Section 3.1, we show the galaxies that are found to host recent BH mergers. In Section 3.2, we investigate how these galaxies evolve before and after the BH merger event. In Section 3.3, we characterize the morphologies of host galaxies, with an emphasis on observationally identifying recent galaxy mergers, and in Section 3.4 we investigate the connection between recent BH mergers and the host galaxy SFRs. In Section 4, we estimate the impact of incorporating inspiral and binary hardening time-scales into the BH merger model, particularly regarding the rate and masses of merging BHs detectable via GWs. Finally, in Section 5 we summarize our conclusions.

2 METHOD

In this work, we primarily use the Illustris1 suite of simulations (Nelson et al. 2015), run using the moving-mesh code arepo (Springel 2010), focusing on the highest resolution simulation run with a periodic box |$106.5\, {\rm Mpc}$| on a side. This model has target gas cell mass |$m_{\rm {gas}}=1.26 \times 10^6\, {\mathrm{M}_\odot }$| and dark matter particle mass |$m_{\rm {DM}}=6.26 \times 10^6\, {\mathrm{M}_\odot }$|⁠, with a standard ΛCDM cosmology: Ωm, 0 = 0.2726, ΩΛ, 0 = 0.7274, Ωb, 0 = 0.0456, σ8 = 0.809, ns = 0.963, and |$H_0=70.4 \, \rm {km}\, \rm {s}^{-1} \, \rm {Mpc}^{-1}$| (consistent with Hinshaw et al. 2013).

The Illustris simulations include detailed models of the physics involved in galaxy formation and evolution, including primordial and metal-line cooling with a time-dependent UV background (Faucher-Giguère et al. 2009) including self-shielding (Rahmati et al. 2013); star formation (SF) with associated supernova feedback (Springel & Hernquist 2003; Springel et al. 2005a); stellar evolution, gas recycling, and metal enrichment (see Wiersma et al. 2009) with mass and metal-loaded outflows (see Oppenheimer & Davé 2008; Okamoto et al. 2010; Puchwein & Springel 2013). For a more complete description of the physics incorporated into these simulations, see Vogelsberger et al. (2014a), Genel et al. (2014), and Sijacki et al. (2015).

Of particular importance to this project is the BH model, which we briefly summarize here (for more complete details, see Sijacki et al. 2015). BHs are treated as collision-less sink particles, which are seeded at mass Mseed = 105h−1 M into any halo with |$M_{\rm {halo}} \gt 5 \times 10^{10}\, h^{-1} {\mathrm{M}_\odot }$| that does not already contain a BH, loosely motivated by the direct collapse model (see e.g. Haehnelt & Rees 1993; Loeb & Rasio 1994; Bromm & Loeb 2003; Begelman, Volonteri & Rees 2006; Regan & Haehnelt 2009), but intended to remain broadly consistent with lighter seed formation models followed by relatively efficient mass growth. After seeding, Illustris grows BHs by gas accretion modelled by a Bondi–Hoyle-like rate (Bondi & Hoyle 1944; Bondi 1952). During gas accretion, three modes of BH feedback are included that deposit energy on to the surrounding material or alter its cooling rate (‘quasar’, ‘radio’, and ‘radiative’, depending on the accretion efficiency). BHs that come within each other’s smoothing lengths merge together into a single BH. The simulation saves the exact time and masses of the merger, giving us precise information for each merger that occurs within the simulation. In addition to the merger events themselves, this provides us with a complete merger tree for every BH, which we use to extract the most recent merger events that occur in each galaxy. We note that Illustris does not incorporate a hardening time for the binary, and instead assumes instantaneous coalescence once the BHs are within smoothing lengths of each other, the implications of which we investigate in Section 4.

In addition to the BH model, Illustris identifies dark matter haloes with a friends-of-friends (FOF) halo finder (Davis et al. 1985), and assigns baryonic matter to the same group as the nearest dark matter particle. After the FOF finder, the subfind algorithm (Springel et al. 2001; Nelson et al. 2015) produces a catalogue of gravitationally self-bound substructures, and the mergers between these subhalos are tracked using the sublink catalogue (Rodriguez-Gomez et al. 2015). The sublink algorithm sums a merit function for every particle found in common between subhalos in different snapshots, and identifies the subhalo descendant as the subhalo with the highest summed merit score. Thus, we have a complete list of galaxy properties (from subfind) and a complete merger tree for all resolved subhaloes (sublink) from which we can connect BHs and their host galaxies, as well as the complete merger history for each such host.

3 GALAXY HOSTS OF SUPERMASSIVE BLACK HOLE MERGERS

3.1 Host galaxies

We start our analysis by looking at the galaxies that host SMBH mergers, in particular focusing on the scaling relation between BH mass (MBH) and galaxy stellar mass (M*). In Fig. 1, we show the relation at z = 0.5, which approximately matches the observed local scaling relation (for more analysis of the scaling relation in Illustris, see Sijacki et al. 2015). Here, we highlight the galaxies that have, in the past 300 Myr, hosted mergers at scales expected to be detectable by PTA or LISA. In blue, we show LISA-detectable mergers, defined as BH mergers with chirp mass |$M_{\rm {chirp}} := (M_1 M_2)^{3/5}/(M_1+M_2)^{1/5} \lt 10^7\, {\mathrm{M}_\odot }$| (with no minimum since the LISA sensitivity extends below the seed mass of BHs in the simulation), and in red the PTA-detectable mergers defined as those with |$M_{\rm {chirp}} \gt 10^8\, {\mathrm{M}_\odot }$|⁠. We also subdivide the mergers based on the ratio between the larger (M1) and smaller (M2) masses involved in the merger, with minor BH mergers defined as BH mass ratios of 0.1 < M2/M1 < 0.25, and major BH mergers defined by M2/M1 > 0.25. As expected, the high-mass PTA-detectable mergers are found exclusively in high-mass galaxies, while lower mass LISA-detectable mergers tend to be found in lower mass galaxies.

MBH−M* relation at z = 0.5. The contours show the overall distribution of galaxies in the Illustris simulation (the outermost contour encompasses $\sim 99.8 {{\ \rm per\ cent}}$ of galaxies, while the innermost encompasses $\sim 16 {{\ \rm per\ cent}}$), with a comparison to the local relation (Kormendy & Ho 2013, the grey-dashed line). The coloured points are the subsample of galaxies that have hosted a BH merger in the past 300 Myr, divided into major (M2/M1 > 0.25, the filled circles) versus minor (0.1 < M2/M1 < 0.25, the open circles) BH mergers, and LISA-type ($M_{\rm {chirp}} \lt 10^7\, {\mathrm{M}_\odot }$, blue) versus PTA-type (Mchirp > 108 M⊙) mergers. MBH here is defined as the mass of the central (i.e. most massive) BH in the galaxy. The solid and open arrows show the typical difference between central BH mass and chirp mass of the merger. The three larger symbols (the star, the upward triangle, and the downward triangle) correspond to specific targets investigated in more detail in Figs 5–7.
Figure 1.

MBHM* relation at z = 0.5. The contours show the overall distribution of galaxies in the Illustris simulation (the outermost contour encompasses |$\sim 99.8 {{\ \rm per\ cent}}$| of galaxies, while the innermost encompasses |$\sim 16 {{\ \rm per\ cent}}$|⁠), with a comparison to the local relation (Kormendy & Ho 2013, the grey-dashed line). The coloured points are the subsample of galaxies that have hosted a BH merger in the past 300 Myr, divided into major (M2/M1 > 0.25, the filled circles) versus minor (0.1 < M2/M1 < 0.25, the open circles) BH mergers, and LISA-type (⁠|$M_{\rm {chirp}} \lt 10^7\, {\mathrm{M}_\odot }$|⁠, blue) versus PTA-type (Mchirp > 108 M) mergers. MBH here is defined as the mass of the central (i.e. most massive) BH in the galaxy. The solid and open arrows show the typical difference between central BH mass and chirp mass of the merger. The three larger symbols (the star, the upward triangle, and the downward triangle) correspond to specific targets investigated in more detail in Figs 57.

We note there are three outliers: LISA-type mergers which occur in high-mass galaxies, with |$M_{\rm {BH}} \gt \gt 10^8\, {\mathrm{M}_\odot }$|⁠, well above our mass scale of |$M_{\rm {chirp}} \lt 10^7\, {\mathrm{M}_\odot }$|⁠. These represent massive galaxies with high-mass central BHs (i.e. high MBH) that host a merger between two low-mass satellite BHs (i.e. low BH merger masses, and thus low Mchirp). This is a rare occurrence, consisting of only |$\sim 2{{\ \rm per\ cent}}$| of LISA-type mergers within Illustris.2 For the remainder of cases, however, we find that Mchirp and MBH are similar, with a mean log(MBH/Mchirp) = 0.5 ± 0.04 (0.6 ± 0.04) for major (minor) BH merger hosts (characterized by filled and open arrows in Fig. 1, respectively). However, Illustris only seeds each galaxy with a single BH, and the minimum masses involved are limited by the seed mass (Mseed = 105h−1 M); thus the frequency of outlier mergers (of BHs undermassive relative to their host galaxy) should be expected to be higher than predicted here. Overall, this suggests that galaxies hosting BH mergers have typical MBH/M* ratios and generally have Mchirp within ∼0.5 dex of MBH, cen. Furthermore, while both LISA and PTAs are limited in the mass ranges they can probe, between the two types of detectors we can expect to probe typical host galaxies across a wide range of masses (covering the entire range resolved by the Illustris simulation).

However, we find that galaxies hosting SMBH mergers are not morphologically typical, although their masses are. To demonstrate this qualitatively, in Fig. 2 we show mock galaxy images in observed frame James Webb Space Telescope bands (neglecting the effect of dust) at z = 0.5 for hosts of recent major, massive BH mergers (ordered by galaxy stellar mass, from least massive in the upper left to most massive in the lower right). Here, we define a recent merger as having occurred in the previous 300 Myr, a major BH merger as one with a mass ratio M2/M1 > 0.25, and a massive BH merger as one in which |$M_2 \gt 10^7\, {\mathrm{M}_\odot }$|⁠. We note that this definition of a massive BH merger is less stringent than the PTA constraint used in Fig. 1, to provide us with a more reasonable sample size for the following statistical analysis, but this change has no qualitative impact on Fig. 1.

Montage (images from Illustris public release, see Nelson et al. 2015) of z = 0.5 galaxies that hosted a BH merger with $M_2 \gt 10^7\, {\mathrm{M}_\odot }$ (thus roughly PTA-type mergers) and M2/M1 > 0.25 in the prior 300 Myr, showing a mix of disturbed and relatively relaxed morphologies. The galaxies are ordered by stellar mass, from least massive in the upper left to most massive in the lower right, with log(M*/M⊙) of the galaxy and a 20 ckpc scale bar provided in each panel. Here, we see the majority of hosts show clearly disturbed morphologies, including asymmetries, tidal tails, light shells, and dual nuclei, indicative of recent merger activity.
Figure 2.

Montage (images from Illustris public release, see Nelson et al. 2015) of z = 0.5 galaxies that hosted a BH merger with |$M_2 \gt 10^7\, {\mathrm{M}_\odot }$| (thus roughly PTA-type mergers) and M2/M1 > 0.25 in the prior 300 Myr, showing a mix of disturbed and relatively relaxed morphologies. The galaxies are ordered by stellar mass, from least massive in the upper left to most massive in the lower right, with log(M*/M) of the galaxy and a 20 ckpc scale bar provided in each panel. Here, we see the majority of hosts show clearly disturbed morphologies, including asymmetries, tidal tails, light shells, and dual nuclei, indicative of recent merger activity.

Fig. 2 shows that from a simple visual analysis, many of the galaxies hosting massive BH major mergers exhibit disturbed morphologies, providing evidence for galaxy mergers in the form of dual nuclei, secondary satellite galaxies, shells, or tidal features. In Figs 3 and 4, we show the galaxies hosting massive major BH mergers at z = 1 and z = 2, respectively. As in the case of the general galaxy population, we see that high-redshift merger hosts tend to be much bluer than those at low-z, and the low-mass galaxies (upper left-hand panels) tend to be bluer than the higher mass galaxies (lower right-hand panels). The higher redshift images show disturbed morphologies, but these disturbances are less ubiquitous than at low z (in particular regarding dual nuclei and distinct satellite galaxies). Thus, we find that galaxies that host massive major BH mergers can show morphological evidence for a galaxy merger, with the strongest evidence occurring at lower redshifts (which we investigate quantitatively in Section 3.3).

As in Fig. 2, but at z = 1. As in Fig. 2, we see clear evidence of disturbed morphologies, and also note that, as expected, the galaxies are noticeably bluer than those at z = 0.5.
Figure 3.

As in Fig. 2, but at z = 1. As in Fig. 2, we see clear evidence of disturbed morphologies, and also note that, as expected, the galaxies are noticeably bluer than those at z = 0.5.

As in Fig. 2, but at z = 2. As in Figs 2 and 3 we see disturbed morphologies, but less universally so across the full sample, with some galaxies displaying unperturbed disc-like structure.
Figure 4.

As in Fig. 2, but at z = 2. As in Figs 2 and 3 we see disturbed morphologies, but less universally so across the full sample, with some galaxies displaying unperturbed disc-like structure.

3.2 Host evolution

To better understand how galaxies hosting SMBH mergers will evolve with time, in Figs 58 we track individual galaxies, over ∼ 1 Gyr, spanning both before and after the merger event. In each figure we show the evolution of host galaxy properties (MDM, M*, Mgas, SFR, and U – V colour) in the leftmost panel, and the galaxy morphology at each snapshot in the right-hand panels. The morphology images are in sequential order starting at the upper left-hand panel, with each panel labelled with the time relative to the BH merger.

The evolutionary history of a galaxy hosting a high-mass major BH merger. Left-hand panel: MDM, M*, Mgas, SFR, and U – V colour for the galaxy as a function of time since the BH merger. Right-hand panels: Stellar composite images from sequential snapshots, corresponding to the circles marking snapshots in the left-hand panel. The time since the BH merger is marked in each panel, with the BH merging at the galaxy centre at time t = 0 Myr. Initially, we see a well-defined nucleus, followed by a disturbed nucleus as a secondary galaxy has its first passage, followed by a defined nucleus with light shells showing evidence of the recent merger. Note that this galaxy corresponds to the star symbol used in Figs 1, 9, 10, and 12.
Figure 5.

The evolutionary history of a galaxy hosting a high-mass major BH merger. Left-hand panel: MDM, M*, Mgas, SFR, and U – V colour for the galaxy as a function of time since the BH merger. Right-hand panels: Stellar composite images from sequential snapshots, corresponding to the circles marking snapshots in the left-hand panel. The time since the BH merger is marked in each panel, with the BH merging at the galaxy centre at time t = 0 Myr. Initially, we see a well-defined nucleus, followed by a disturbed nucleus as a secondary galaxy has its first passage, followed by a defined nucleus with light shells showing evidence of the recent merger. Note that this galaxy corresponds to the star symbol used in Figs 1, 9, 10, and 12.

In Fig. 5, we show the evolutionary history of a galaxy that hosts a merger between BHs with mass |$M_1=6.3 \times 10^8\, {\mathrm{M}_\odot }$| and |$M_2=3.4 \times 10^8\, {\mathrm{M}_\odot }$|⁠. Pre-merger we see a well-defined stellar nucleus. Near the BH merger event we see a disturbed nucleus, produced as the satellite galaxy has a first passage and when the post-merger galaxy settles down. At later times, we see clear light shells as remnants of the galaxy merger, surviving for an extended period (∼500 Myr). We see an increase in SF, with a corresponding change to bluer colour, near the merger time, as the galaxy merger triggers a burst of SF (Springel, Di Matteo & Hernquist 2005b). This is followed by a significant drop in SF and gas mass (Mgas decreases by ∼1.5 orders of magnitude in the |$500\, {\rm Myr}$| following the merger), and a corresponding reddening of the galaxy as the SF falls off.

In Fig. 6, we show the history of a galaxy hosting a more massive BH merger (⁠|$M_1=2.3\times 10^9\, {\mathrm{M}_\odot }$| and |$M_2=1.2\times 10^9\, {\mathrm{M}_\odot }$|⁠). Here, we again see a decreasing SFR and mildly increasing U – V colour, but with an overall SFR that is much lower, and the gas supply remains relatively high throughout the merger period. Morphologically, at the time nearest the merger we see evidence for a dual nucleus, and again there are light shells that survive long past the merger event (though fainter than in Fig. 5).

As in Fig. 5, but for another galaxy. Here, we see less disturbance in the nucleus, but none the less long-lived, well-defined light shells. This galaxy corresponds to the upward pointing triangle symbol in Figs 1, 9, 10, and 12.
Figure 6.

As in Fig. 5, but for another galaxy. Here, we see less disturbance in the nucleus, but none the less long-lived, well-defined light shells. This galaxy corresponds to the upward pointing triangle symbol in Figs 1, 9, 10, and 12.

In Fig. 7, we consider the history of a galaxy hosting a less massive, minor BH merger (⁠|$M_1=1.1 \times 10^8\, {\mathrm{M}_\odot }$| and |$M_2=2.3 \times 10^7\, {\mathrm{M}_\odot }$|⁠, note that this is near the threshold for a major BH merger), which occurs in a very gas-rich galaxy. Here, we again see a reddening of the galaxy; also the morphology is highly disturbed immediately following the BH merger, with tidal features visible for several hundred Myr.

As in Fig. 5, but for a particularly gas-rich galaxy. Here we do not see any light shells, but do observe a highly disturbed morphology caused by the galaxy merger. This galaxy corresponds to the downward pointing triangle symbol in Figs 1, 9, 10, and 12.
Figure 7.

As in Fig. 5, but for a particularly gas-rich galaxy. Here we do not see any light shells, but do observe a highly disturbed morphology caused by the galaxy merger. This galaxy corresponds to the downward pointing triangle symbol in Figs 1, 9, 10, and 12.

We also consider the impact of the BH model used in the simulation. In particular, the Illustris simulation has previously been found to have active galactic nucleus (AGN) feedback that could too easily heat and dilute the central gas in galaxy groups, whereas IllustrisTNG uses a new feedback model to resolve this issue (Weinberger et al. 2017). However, the effect is strongest among highest mass galaxies, so we do not expect significant qualitative impact. None the less, to test the sensitivity to the assumptions made regarding BH modelling in the simulation, we also show an evolutionary history from the TNG100 model of the IllustrisTNG3 suite of simulations (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Springel et al. 2018; Pillepich et al. 2018b), at z = 0 in Fig. 8. Prior to the merger in this example, we see a well-defined spiral galaxy, which is noticeably disturbed near the time of the BH merger. As in Figs 5 and 6, we find a strong decrease in SF post-merger, with a correspondingly significant increase in U – V colour. The similarity between these histories from Illustris and TNG100 suggests that the qualitative results found here do not exhibit a strong sensitivity to the BH model of the simulation. Overall, this subsample suggests that galaxies hosting a massive BH merger generally show morphological evidence of a merger, and tend to redden and decrease in SF during the post-merger phase, which we quantify for the full population in Sections 3.3 and 3.4.

As in Fig. 5, but for a galaxy from the TNG100 simulation at z = 0. As in Fig. 7 we see a clearly disturbed morphology corresponding to the BH merger, and we also witness significant reddening which occurs following the merger.
Figure 8.

As in Fig. 5, but for a galaxy from the TNG100 simulation at z = 0. As in Fig. 7 we see a clearly disturbed morphology corresponding to the BH merger, and we also witness significant reddening which occurs following the merger.

3.3 Host morphologies

Having seen the morphologies of galaxies hosting major massive BH mergers in Figs 24 and individual evolutions of several sample galaxies in Figs 58, we now consider the morphology statistics for the full host sample relative to the full galaxy population.

In order to quantify each galaxy, we use Gini (a measure of how evenly the galaxy flux is distributed, from 1 if all flux is in a single pixel to 0 for a uniform distribution across all pixels; see Lotz, Primack & Madau 2004, for more details) and M20 (defined as the second-order moment of the brightest 20 per cent of pixels) parameters as simple measures of morphology. In Fig. 9, we show the Gini–M20 relation for galaxies at z = 0.5. The coloured contours show the distribution of all galaxies with |$M_* \gt 10^{9.7}\, {\mathrm{M}_\odot }$| in Illustris. ‘Normal’ non-merging galaxies will tend to lie towards the bottom right of the panel, while merging galaxies will be more likely to lie in the upper left portion (Lotz et al. 2004); merging galaxy are more likely to have a high Gini coefficient (i.e. stellar light concentrated in a smaller fraction of total pixels) and a high M20 (i.e. the brightest 20 per cent of pixels are less centrally concentrated) compared to a non-merging galaxy. Here, we plot a dashed line roughly differentiating between morphologies suggesting a merger and those suggesting a quiescent galaxy (Snyder et al. 2019). The data points show the location of individual galaxies hosting major massive BH mergers (red) and major low-mass BH mergers (blue) in the past 300 Myr, and the M20-binned mean Gini values (the coloured error bars) compared to the mean of the full population (the black error bars). Here, we find that recent major BH mergers tend to lie ∼1σ above the general population, and >2σ above for high-mass mergers in high-M20 galaxies, consistent with the earlier qualitative results.

The Gini–M20 relation at z = 0.5 for galaxies with M* ⪆ 109.7 M⊙ (contours, from the outermost encompassing $\sim 99.7 {{\ \rm per\ cent}}$ of galaxies to the innermost encompassing $\sim 30 {{\ \rm per\ cent}}$), for galaxies that have hosted a BH merger in the past 300 Myr with BH mass ratio M2/M1 > 0.1 (the red points) and $M_2/M_1 \gt 0.1 ; M_2 \gt 10^7\, {\mathrm{M}_\odot }$ (the blue points), and error bars that show the mean and standard deviation for each population. Galaxies with a recent massive major BH merger tend to have higher Gini than the general galaxy sample, corresponding to what we would expect a galaxy merger to exhibit. Low-mass major BH mergers have a slightly higher Gini than the full galaxy sample, but the trend is much weaker and only at the order of ∼1σ. The star, the upward triangle, and the downward triangle correspond to the individual histories plotted in Figs 5, 6, and 7, respectively.
Figure 9.

The Gini–M20 relation at z = 0.5 for galaxies with M* ⪆ 109.7 M (contours, from the outermost encompassing |$\sim 99.7 {{\ \rm per\ cent}}$| of galaxies to the innermost encompassing |$\sim 30 {{\ \rm per\ cent}}$|⁠), for galaxies that have hosted a BH merger in the past 300 Myr with BH mass ratio M2/M1 > 0.1 (the red points) and |$M_2/M_1 \gt 0.1 ; M_2 \gt 10^7\, {\mathrm{M}_\odot }$| (the blue points), and error bars that show the mean and standard deviation for each population. Galaxies with a recent massive major BH merger tend to have higher Gini than the general galaxy sample, corresponding to what we would expect a galaxy merger to exhibit. Low-mass major BH mergers have a slightly higher Gini than the full galaxy sample, but the trend is much weaker and only at the order of ∼1σ. The star, the upward triangle, and the downward triangle correspond to the individual histories plotted in Figs 5, 6, and 7, respectively.

To better characterize this difference, we use the merger statistic S (see Snyder et al. 2019), defined as
(1)
(where G is the Gini coefficient) that provides a single parameter to quantify how merger-like the galaxy morphology is, with a positive S(G, M20) corresponding to a merger-like galaxy and negative S(G, M20) being non merger-like.

In Fig. 10, we characterize the sensitivity of the galaxy morphology on time since the BH merger, with the contours showing the time since the most recent BH merger (Δt) against S(G, M20). The error bars show the mean and standard deviation of S(G, M20) binned by time since the most recent non-major BH merger (black), low-mass major BH merger (blue), and high-mass major BH merger (red). For the binned data, we consider the time of the most recent merger, but do not include galaxies if the most recent merger is more ‘significant’ than specified for each population (i.e. the red points calculate the time since the most recent massive major BH merger; the blue points calculate the time since the most recent low-mass major BH merger, but ignore galaxies where there was a more recent high-mass major BH merger).

Merger statistic S (see equation 1) for z = 0.5 galaxies as a function of time since the most recent BH merger within the galaxy (contours, ranging from $\sim 80 {{\ \rm per\ cent}}$ of galaxies to $\sim 11 {{\ \rm per\ cent}}$). Also shown are the mean and standard deviation of S statistics binned by time since the most recent major BH merger (M2/M1 > 0.25, the blue points), since the most recent massive major BH merger ($M_2/M_1 \gt 0.25; M_2 \gt 10^7\, {\mathrm{M}_\odot }$, the red points), and since the most recent non-major BH merger (M2/M1 < 0.25, the black points).
Figure 10.

Merger statistic S (see equation 1) for z = 0.5 galaxies as a function of time since the most recent BH merger within the galaxy (contours, ranging from |$\sim 80 {{\ \rm per\ cent}}$| of galaxies to |$\sim 11 {{\ \rm per\ cent}}$|⁠). Also shown are the mean and standard deviation of S statistics binned by time since the most recent major BH merger (M2/M1 > 0.25, the blue points), since the most recent massive major BH merger (⁠|$M_2/M_1 \gt 0.25; M_2 \gt 10^7\, {\mathrm{M}_\odot }$|⁠, the red points), and since the most recent non-major BH merger (M2/M1 < 0.25, the black points).

For Δt > 1 Gyr, we find no difference between the merger populations, and for Δt > 500 Myr, there are only very minor differences (<<1σ). However, for Δt < 300 Myr, we find that massive major BH merger hosts have a significantly higher S(G, M20), while low-mass major BH merger hosts show a moderately higher S(G, M20), though at <1σ significance.

To further demonstrate this correlation between recent mergers and host morphology, in Fig. 11 we show the z = 0.5 probability (top) and cumulative (bottom) distribution functions of S(G, M20) for all galaxies (black) and for galaxies that, in the past 300 Myr, have hosted a low-mass major BH merger (blue) or high-mass major BH merger (red). As in Figs 9 and 10, we see that low-mass major BH mergers lie ∼1σ above the general galaxy sample, while the high-mass major BH mergers are ∼1.5σ higher, confirming that hosts of major BH mergers are significantly morphologically disturbed.

Probability distribution function (top panels) and cumulative distribution function (bottom panels) of merger statistic S at z = 0.5 (left), z = 1 (middle), and z = 2 (right) for all galaxies (black), galaxies that have hosted high-mass (red) or low-mass (blue) major BH merger in the past 300 Myr. The error bars show the mean and standard deviation of each sample. For most BH mergers, the distribution of the S-statistic is relatively unchanged compared to the general distribution for all galaxies. For high-mass major BH mergers occurring within the past ∼500 Myr (and especially within 300 Myr), the galaxies tend to show morphological evidence of having recently merged. The red-dashed and blue-dashed curves show the weighted distribution of the full galaxy population to provide a mass-matched sample to the merger population (i.e. match of the Mtot distribution of the red and blue populations, respectively).
Figure 11.

Probability distribution function (top panels) and cumulative distribution function (bottom panels) of merger statistic S at z = 0.5 (left), z = 1 (middle), and z = 2 (right) for all galaxies (black), galaxies that have hosted high-mass (red) or low-mass (blue) major BH merger in the past 300 Myr. The error bars show the mean and standard deviation of each sample. For most BH mergers, the distribution of the S-statistic is relatively unchanged compared to the general distribution for all galaxies. For high-mass major BH mergers occurring within the past ∼500 Myr (and especially within 300 Myr), the galaxies tend to show morphological evidence of having recently merged. The red-dashed and blue-dashed curves show the weighted distribution of the full galaxy population to provide a mass-matched sample to the merger population (i.e. match of the Mtot distribution of the red and blue populations, respectively).

In addition, we show the equivalent distributions at z = 1 (middle panels) and z = 2 (right-hand panels). At z = 1, we find that low-mass major BH mergers are no longer significantly different from the general population, while the high-mass major BH mergers are still ∼1σ above the general population. By z = 2, we see that the galaxies typically have the same S(G, M20) distribution regardless of whether they have hosted a recent BH merger or not. Thus, we conclude that recent BH mergers correlate with merger morphologies at low redshift, though above z ∼ 1 there is no longer any significant correlation. We confirm the statistical significance of these claims using a Kolmogorov–Smirnov test, finding that both major BH merger populations at z = 0.5, and the high-mass major BH merger population at z = 1, each have a |$\lt 0.2{{\ \rm per\ cent}}$| likelihood of being drawn from the same distribution as the full galaxy sample (Table 1). We also checked other Δt, and the only longer time-scale for which a merger signal survives is for massive major BH merger hosts at z = 0.5, in which the signal survives for ∼500 Myr (the other cases differ from the general population by <0.5σ). We note that this survival time of ∼ 500 Myr for a morphological signal is on the order of the galaxies’ dynamical time, consistent with a picture in which the galaxies gradually return to a relaxed state.

Table 1.

The KS-test likelihood that the target sample (high- or low-mass major BH merger hosts) comes from the same S(G, M20) distribution as the full galaxy sample (see Fig. 11). At z = 0.5, low mass major BH merger hosts have a signal that survives ∼300 Myr, while the signal in high mass major BH merger hosts survives for ∼500 Myr. At higher redshifts, the signal is weaker; only the high-mass hosts have a signal at z = 1, which survives for ∼300 Myr.

Time since BH merger (Myr)Low-mass major BH mergersHigh-mass major BH mergers
z = 2.00-3000.700.64
z = 1.00-3000.27|$\boxed{1.8 \times 10^{-3}}$|
z = 0.50-300|$\boxed{3.6 \times 10^{-4}}$||$\boxed{5.0 \times 10^{-6}}$|
300-5000.45|$\boxed{3.4 \times 10^{-4}}$|
500-10000.220.16
Time since BH merger (Myr)Low-mass major BH mergersHigh-mass major BH mergers
z = 2.00-3000.700.64
z = 1.00-3000.27|$\boxed{1.8 \times 10^{-3}}$|
z = 0.50-300|$\boxed{3.6 \times 10^{-4}}$||$\boxed{5.0 \times 10^{-6}}$|
300-5000.45|$\boxed{3.4 \times 10^{-4}}$|
500-10000.220.16
Table 1.

The KS-test likelihood that the target sample (high- or low-mass major BH merger hosts) comes from the same S(G, M20) distribution as the full galaxy sample (see Fig. 11). At z = 0.5, low mass major BH merger hosts have a signal that survives ∼300 Myr, while the signal in high mass major BH merger hosts survives for ∼500 Myr. At higher redshifts, the signal is weaker; only the high-mass hosts have a signal at z = 1, which survives for ∼300 Myr.

Time since BH merger (Myr)Low-mass major BH mergersHigh-mass major BH mergers
z = 2.00-3000.700.64
z = 1.00-3000.27|$\boxed{1.8 \times 10^{-3}}$|
z = 0.50-300|$\boxed{3.6 \times 10^{-4}}$||$\boxed{5.0 \times 10^{-6}}$|
300-5000.45|$\boxed{3.4 \times 10^{-4}}$|
500-10000.220.16
Time since BH merger (Myr)Low-mass major BH mergersHigh-mass major BH mergers
z = 2.00-3000.700.64
z = 1.00-3000.27|$\boxed{1.8 \times 10^{-3}}$|
z = 0.50-300|$\boxed{3.6 \times 10^{-4}}$||$\boxed{5.0 \times 10^{-6}}$|
300-5000.45|$\boxed{3.4 \times 10^{-4}}$|
500-10000.220.16

We note that massive BH mergers tend to occur in high-mass galaxies (as shown in Fig. 1), which could potentially contribute to the difference between the general galaxy population and the massive major BH merger hosts. As such, we compute a distribution of the general galaxy population that is a mass-matched equivalent to the merger host sample. The mass-matched equivalent is found by weighting each galaxy in the full population sample by the ratio between a Gaussian best fit to the target sample’s Mtot distribution (neglecting the 5 per cent furthest outliers in each direction) and the full galaxy population’s Mtot distribution. The result is a weighted PDF of the full galaxy sample but with an equivalent galaxy–mass distribution to the target sample. We plot the weighted distributions for these mass-matched samples as dashed lines in Fig. 11. At all redshifts, the mass-matched samples (dashed lines) match the full galaxy population (solid black), which confirms that the S(G, M20) distributions are not dependent on galaxy mass, and the significant increase among BH merger hosts is genuinely due to merging.

3.4 Star formation rates

Given the potential for galaxy mergers to influence typical host SF, and the substantial post-merger decrease in SFR seen in the individual histories in Figs 5, 6, and 8, we also look for correlations between BH mergers and host SFRs. In Fig. 12, we plot the specific SFR (sSFR≔SFR/M*) versus total galaxy mass at z = 0.5. As expected, we see that the sSFR tends to drop among high-mass galaxies (the contours and the black points) since high-mass galaxies tend to be less star forming. We also plot the populations for the high- (red) and low- (blue) mass major BH merger hosts, together with the Mtot-binned mean and standard deviation (error bars). High-mass major BH mergers tend to have similar sSFR to other equivalent mass galaxies (the red versus black error bars). However, low-mass major BH mergers tend to have slightly higher sSFR for their masses (the blue versus black error bars), suggesting that those mergers tend to trigger some additional SF.

Specific star formation rate versus total galaxy mass for the full galaxy sample (contours, with the outermost encompassing $\sim 99.7 {{\ \rm per\ cent}}$ of galaxies, and the innermost encompassing $\sim 46 {{\ \rm per\ cent}}$ of galaxies) and high-mass (red) and low-mass (blue) major BH merger hosts. We find that high-mass major BH mergers tend towards high-mass hosts and thus lower sSFR, but are fully consistent with a mass-matched sample from the full galaxy population.
Figure 12.

Specific star formation rate versus total galaxy mass for the full galaxy sample (contours, with the outermost encompassing |$\sim 99.7 {{\ \rm per\ cent}}$| of galaxies, and the innermost encompassing |$\sim 46 {{\ \rm per\ cent}}$| of galaxies) and high-mass (red) and low-mass (blue) major BH merger hosts. We find that high-mass major BH mergers tend towards high-mass hosts and thus lower sSFR, but are fully consistent with a mass-matched sample from the full galaxy population.

We confirm this in Fig. 13, where we show the cumulative distribution function for sSFR (left column) of the full population (black) compared to the major BH merger populations (red and blue). Given the strong correlation between sSFR and galaxy mass (Fig. 12), it is necessary to compare mass-matched samples rather than comparing with the full galaxy population. As in Fig. 11, for each target sample (line colour) we show the distribution for a mass-matched equivalent from the full galaxy sample (the dashed lines). Here, we find that the hosts of high-mass major BH mergers (solid red) are equivalent to the mass-matched sample from the full population (dashed red). This suggests that hosting a major massive BH merger either tends to have minimal impact on the sSFR, or the effect is short compared to the 300 Myr used for selecting the host sample. Conversely, we find that hosts of low-mass major BH mergers tend to have slightly higher sSFR than the mass-matched sample, suggesting that the galaxy merger triggers an increase in SFR (though only on the order of ∼0.5σ). We find qualitatively similar results across redshifts: higher redshifts have higher sSFRs overall, but high-mass major BH merger hosts match an equivalent-mass sample, whereas low-mass major BH merger hosts tend to have slightly higher sSFR than their equivalent-mass sample.

Cumulative distribution function of sSFR (left) and Mgas/Mbaryon (right). The black lines show the distribution for the full galaxy population. The solid coloured lines show a sample of galaxies hosting BH mergers. The dashed coloured lines show the distribution of the full galaxy population mass weighted to match the distribution of the targeted BH merger sample. Hosts of low-mass major BH mergers tend to have slightly increased sSFR (i.e. SF boost from the galaxy merger). Hosts of high-mass major BH mergers do not exhibit such a boost. At low z, the low gas fraction could explain the lack of SFR boost; at high z, however, the galaxies are gas rich, suggesting another factor must play a role, such as feedback from the high-mass BHs.
Figure 13.

Cumulative distribution function of sSFR (left) and Mgas/Mbaryon (right). The black lines show the distribution for the full galaxy population. The solid coloured lines show a sample of galaxies hosting BH mergers. The dashed coloured lines show the distribution of the full galaxy population mass weighted to match the distribution of the targeted BH merger sample. Hosts of low-mass major BH mergers tend to have slightly increased sSFR (i.e. SF boost from the galaxy merger). Hosts of high-mass major BH mergers do not exhibit such a boost. At low z, the low gas fraction could explain the lack of SFR boost; at high z, however, the galaxies are gas rich, suggesting another factor must play a role, such as feedback from the high-mass BHs.

One possible explanation for why low-mass galaxy mergers may boost the sSFR, while high-mass galaxy mergers do not, is simply due to the gas richness of the merging galaxies. High-mass galaxies, where the high-mass major BH mergers tend to occur, generally have low gas fractions, while the lower mass galaxies that host lower mass BH mergers would be expected to be more gas rich, and thus have more potential to drive a burst of SF (e.g. Genel et al. 2014; Vogelsberger et al. 2014b). Indeed, we find this to be the case at z = 0.5: galaxies hosting high-mass major BH mergers tend to have Mgas/Mb ∼ 0.43 ± 0.22, compared to low-mass major BH merger hosts, which have Mgas/Mb ∼ 0.80 ± 0.16 (see right-hand panels of Fig. 13).

However, we note that this cannot be the sole explanation. In the middle and lower panels of Fig. 13 we show the sSFR and gas fractions at z = 2 and z = 4. Similar to low z, the high-mass major BH merger hosts have the same sSFR as a mass-matched sample, while the low-mass major BH merger hosts tend to have slightly higher sSFR than their mass-matched sample. However, the gas fractions are much higher, with the high-mass hosts at z = 2 having higher gas fractions than the low-mass hosts at z = 0.5. Thus, at high z we have gas-rich galaxy mergers that still do not appear to affect the SF of galaxies hosting high-mass major BH mergers, while there is a weak correlation for low-mass major BH mergers (sSFR is typically ∼0.5σ higher among the low-mass major BH merger hosts). This suggests that the BHs themselves may play a role in preventing the SF from increasing. In particular, we know that feedback from massive BHs can self-regulate (e.g. King 2003; Springel et al. 2005a; Costa et al. 2014; Sijacki et al. 2015), and can suppress SF, while the low-mass BH may not be capable of providing sufficient feedback to do so. Indeed, we find a correlation between the energy emitted by the BH and the position of the galaxy on the star forming main sequence several hundred Myr later, with higher AGN feedback during the merger corresponding to lower sSFR post-merger. This provides additional evidence that feedback from sufficiently massive BHs can help regulate SF following a major galaxy merger, suppressing the boost in SFR one might otherwise expect from a gas-rich galaxy, as shown in Fig. 5.

4 MERGER TIME-SCALES

Having found evidence that galaxies show morphological signatures of a recent galaxy merger over time-scales on the order of 300–500 Myr, we now consider the typical BH merger time-scales and the impact that may have on our results. Within the Illustris simulation (and indeed many similar simulations), a pair of BHs merge as soon as their separation is less than the particle’s smoothing length, rather than incorporating a coalescence time for the binary. Recently, several works have attempted to estimate the expected coalescence time for binary BHs by post-processing cosmological simulations, finding time-scales on the order of 100s of Myr to Gyr for the binary coalescence time (e.g. Blecha et al. 2016; Rantala et al. 2017; Kelley, Blecha & Hernquist 2017a; Mannerkoski et al. 2019; Sayeb et al. 2021). Additionally, the time for a satellite BH to reach the centre of a galaxy and form a binary with the central BH can also be on the order of 100s of Myr to Gyr, based on the dynamical friction time-scale for infall to the galactic centre (e.g. Volonteri et al. 2020). The galaxy structure, e.g. the existence or lack of a dense stellar core, can affect the time satellite BHs spend at large radii (Tremmel et al. 2018,a,b; Barausse et al. 2020), and directly incorporating dynamical friction into cosmological volumes suggests that the orbital decay time-scale may be both substantial and redshift dependent (Bartlett et al. 2021). The Illustris simulation uses a re-centring scheme whereby BHs are re-positioned towards the local potential minimum; this prevents numerical wandering of BHs, but also means the the full infall time to reach the galaxy centre may be notably underestimated. Furthermore, both simulations (e.g. Bellovary et al. 2019) and observations (e.g. Reines et al. 2020) suggest that BHs in dwarf galaxies may frequently be located offset from the galaxy centre, which could further delay any mergers involving low-mass BHs seeded into Illustris (which are initially placed at the galaxy centre). Overall, this suggests that the Illustris simulation likely overestimates the speed with which BHs merge following the merging of their host galaxies, and properly accounting for this has the potential to impact the expected GW detection rate, shift the peak detection time to lower redshift, and prevent GW hosts from being visibly disturbed. Although a complete investigation into accurately estimating the time delays remains beyond the scope of this paper, we address this, by imposing a delay between when the BH particles merge in the simulation (which is closer to when the BH binary may form) and when the final coalescence and GW emission occur.

In Fig. 14, we plot the predicted rate of SMBH mergers throughout the observed Universe as a function of merger redshift (the solid black line), in terms of the rate at which the signals would reach the Earth:
(2)
The first BH merger in the simulation occurs at z ∼ 7.7 (though we note this is limited by the boxsize, resolution and seed criteria used), and the rate then grows with redshift to a peak at z ∼ 2, followed by a decline at later times, broadly consistent with similar analyses (e.g. Salcido et al. 2016; Katz et al. 2020). We also note that the BH merger rate in Illustris (black) is consistent with that from the TNG300 simulation (grey), except that the larger volume TNG300 has slightly more BH mergers at very high redshift. In addition to the Illustris predictions, we also use the larger volume Illustris TNG300 simulation (Springel et al. 2018) to get predicted BH merger rates (grey line). The larger boxsize provides a less noisy redshift distribution, and very slightly higher rates at the highest redshift. Otherwise, the two simulations are largely equivalent, as we would expect given the similar seeding and merging conditions between the two simulations.
The solid lines: The expected BH merger rate signal as a function of redshift. In black, we show the rates from the original simulation (with a comparison to TNG300 in grey), while in red and blue we show the rates if we were to assume a fixed delay time of 500–1000 Myr, respectively. Incorporating a delay time has a significant impact on the earliest BH mergers, but the low-redshift merger rate is largely unaffected by a delay time. The dashed/dotted lines: The BH merger rate signal for mergers corresponding to PTA-type events (dotted; defined as $M_{\rm {chirp}} \gt 10^8\, {\mathrm{M}_\odot }$) and LISA-type events (dashed; defined as $M_{\rm {chirp}} \lt 10^7\, {\mathrm{M}_\odot }$). Note that the rate of LISA-type mergers is limited to those resolved by the Illustris simulation, and that the rate of LISA-detectable mergers is expected to be higher (see main text for further details). Thus, incorporating a BH merger delay can significantly change the merger rates, with a merger delay not only decreasing the rate of high-z mergers (especially at low mass), but also significantly increasing the rate of high-mass mergers at z ∼ 2 (near the merger rate peak).
Figure 14.

The solid lines: The expected BH merger rate signal as a function of redshift. In black, we show the rates from the original simulation (with a comparison to TNG300 in grey), while in red and blue we show the rates if we were to assume a fixed delay time of 500–1000 Myr, respectively. Incorporating a delay time has a significant impact on the earliest BH mergers, but the low-redshift merger rate is largely unaffected by a delay time. The dashed/dotted lines: The BH merger rate signal for mergers corresponding to PTA-type events (dotted; defined as |$M_{\rm {chirp}} \gt 10^8\, {\mathrm{M}_\odot }$|⁠) and LISA-type events (dashed; defined as |$M_{\rm {chirp}} \lt 10^7\, {\mathrm{M}_\odot }$|⁠). Note that the rate of LISA-type mergers is limited to those resolved by the Illustris simulation, and that the rate of LISA-detectable mergers is expected to be higher (see main text for further details). Thus, incorporating a BH merger delay can significantly change the merger rates, with a merger delay not only decreasing the rate of high-z mergers (especially at low mass), but also significantly increasing the rate of high-mass mergers at z ∼ 2 (near the merger rate peak).

In addition to the overall rate, we also show the BH merger rates if the expected inspiral times for the mergers were 500 Myr (the red lines) or 1 Gyr (the blue lines) following when in Illustris a given BH pair merged (which did not incorporate any expected inspiral/hardening time). Here, we use fixed time delays for all BH mergers as a means to show the typical impact that a true merger time-scale may have; a more accurate investigation using variable time delays based on the BH and galaxy properties is beyond the scope of this work, but will be tested in future work. The incorporation of a delay time necessarily has a strong impact on the very early BH mergers. However, we find that at lower redshifts the total BH merger rate evolves slowly enough that adding a delay on the order of 1 Gyr has essentially no impact on the total rate at which we can expect mergers to take place. This is consistent with Weinberger et al. (2018), who showed that in the Illustris TNG300 simulation, BH merger rates binned by mass and redshift have a relatively minor evolution with redshift, and thus concluded that implementing a more realistic inspiral time would not significantly affect the global merger rate.

However, the inspiral time is more than merely a delay on the time at which the BH merger occurs; it also gives the progenitor BHs more time to grow prior to the merger, so we would expect the merger masses to increase with longer delay times (especially as galaxy mergers can often be followed by periods of efficient growth as gas-rich galaxy mergers drive additional gas into the central region where the BHs are accreting). We show the expected effect this will have with the dashed and dotted lines in Fig. 14, where we consider the rate of merging BHs with |$M_{\rm {chirp}} \lt 10^7\, {\mathrm{M}_\odot }$| (dashed, roughly corresponding to LISA-detectable mergers) and |$M_{\rm {chirp}} \gt 10^8\, {\mathrm{M}_\odot }$| (dotted; roughly corresponding to PTA-detectable mergers). We note here that the rate of LISA-type mergers is limited to those BH mergers resolved by the Illustris simulation, which is limited by the seed mass of Mseed = 105h−1 M, whereas the actual detection range of LISA should peak at lower mass scales (see e.g. Dayal et al. 2018; Ricarte & Natarajan 2018; DeGraf & Sijacki 2020); thus the LISA-type merger curve shows how incorporating a delay can impact the rate of the higher mass LISA detections, but should none the less represent a lower limit for the full detection rate.

In the case of the Δt > 0 curves, we assume that the growth of the inspiraling BHs during that delay time is at the equivalent Eddington fraction as the post-merged BH in the original simulation (i.e. we check the fractional accretion in Illustris of the merged BH over the course of Δt, and let each progenitor BH grow by the same fraction). Here, we see that in addition to decreasing the number of high-z BH mergers (especially for low-mass, LISA-type mergers), longer delay times tend to result in more high-mass BH mergers at lower z as the inspiralling BHs have more time to grow prior to merging. In particular, we find that longer merger time-scales that include the associated BH growth lead to an increase in the rate of high-mass (PTA-type) mergers at z ≲ 3.5, especially near the peak at z ∼ 2.

We further consider the effect this may have in Fig. 15, where we show the BH merger rate as a function of chirp mass (Mchirp = (M1M2)3/5/(M1 + M2)1/5), both from the original simulation (black) and for inspiral time delays of Δt = 500 Myr (red) and Δt = 1 Gyr (blue). For these time-delay cases, we consider two growth efficiencies, where the inspiral growth matches either the Eddington fraction or the Bondi fraction of the post-merged BH from the original simulation.

The predicted BH merger rate as a function of chirp mass for the original simulation and after adding 500 and 1000 Myr delay times. For a given Δt, the shaded region spans the range of possible values depending on the accretion efficiency of the binary BHs (between Bondi- and Eddington-like growth; see text for details).
Figure 15.

The predicted BH merger rate as a function of chirp mass for the original simulation and after adding 500 and 1000 Myr delay times. For a given Δt, the shaded region spans the range of possible values depending on the accretion efficiency of the binary BHs (between Bondi- and Eddington-like growth; see text for details).

We expect that the assumption of a fixed Eddington fraction to represent a likely upper bound, as the higher mass from the original simulation (i.e. the combined post-coalescence mass) will tend towards a higher efficiency than the lower mass from the pre-merged BH pair. Conversely, although the simulation uses a Bondi formalism to model accretion, the fixed Bondi rate likely represents an underestimate for the gas accretion, as this assumes each BH is isolated, rather than the combined binary we have here. As such, using both fixed-Eddington and fixed-Bondi fractions provides a range of accretion efficiencies which span a good fraction of the uncertainty within our assumptions.

Here, we again see that although the overall BH merger rate is largely unchanged, the merger masses are affected by the inspiral time, with the longer delays resulting in more massive BH mergers: the mergers peak at |$M_{\rm {chirp}} \sim 2 \times 10^6\, {\mathrm{M}_\odot }$| in the original simulation, |$\sim 7 \times 10^6\, {\mathrm{M}_\odot }$| for a 500 Myr delay, and |$\sim 10^7\, {\mathrm{M}_\odot }$| for a 1 Gyr delay. Furthermore, the delay has the potential to increase the number of massive BH mergers by up to an order of magnitude, depending on the efficiency of growth during the binary phase. However, we again note that the low-mass end is limited by the seed mass of BHs in the Illustris simulation, and would thus expect the actual detection rate for LISA to be significantly higher (see e.g. Dayal et al. 2018; Ricarte & Natarajan 2018; DeGraf & Sijacki 2020).

We show the distribution of BH merger masses in more detail in Fig. 16: in the left-hand panel we plot the distribution of the primary (M1) and secondary (M2) merger masses in Illustris, finding that mergers peak at |$M_1 \sim 5 \times 10^6\, {\mathrm{M}_\odot }$| and |$M_2 \sim 10^6\, {\mathrm{M}_\odot }$|⁠. We also show the fractional change in merger masses when including a delay (⁠|$n(\Delta t=1\, {\rm Gyr}) / n(\Delta t = 0)$|⁠), assuming growth at a fixed Eddington fraction (middle panel) or fixed Bondi fraction (right-hand panel). In both growth models, we find that incorporating a delay time decreases the number of low-mass mergers while increasing high-mass mergers, but with several key differences. The more efficient growth in the fixed Eddington model results in a larger change in the BH merger rate, by as much as 1 dex more high-mass mergers and 1 dex fewer low-mass mergers. However, we note that by construction the fixed Eddington growth model does not affect the mass ratio of the BH mergers, since fixed Eddington fraction growth means both BHs grow proportionally during the binary phase. In contrast, assuming growth based on a fixed Bondi fraction has a smaller overall impact, but does affect the mass ratios. In particular, Bondi growth rates scale with |$M_{\rm {BH}}^2$|⁠, so incorporating a delay assuming Bondi efficiency means that high mass-ratio BH mergers actually feature increasing mass ratios with longer delays (since the less massive progenitor BH is less able to grow during the delay time), while the mass ratio of major mergers are relatively unaffected (since both progenitors grow equally efficiently).

Left-hand panel: The distribution of primary (M1) and secondary (M2) BH merger masses in the Illustris simulation, coloured by the number density of mergers within the simulation. Middle and right-hand panels: The ratio between BH mergers in the original simulation (nIllustris) and after imposing a 1 Gyr delay on all mergers [n(Δt)], assuming growth during the specified Δt continues at the equivalent Eddington rate (middle panel) or equivalent Bondi rate (right-hand panel; see text for details). If growth is Eddington-dependent, both pre-merged black holes grow equivalently, decreasing the number of low-mass mergers with a corresponding increase in high-mass mergers (as each merger moves along a line with slope of 1 in the M2−M1 plane). If growth is Bondi-dependent, then the number of low-mass mergers decreases, but most of the growth occurs in the more massive black hole, such that the main shift is in the M1-direction (i.e. slope less than one). The colour scale shows the number density of BH mergers after imposing a 1 Gyr delay relative to the number density of mergers in the original simulation; red represents an increase in the number of mergers at the given mass scales, while blue shows a decrease.
Figure 16.

Left-hand panel: The distribution of primary (M1) and secondary (M2) BH merger masses in the Illustris simulation, coloured by the number density of mergers within the simulation. Middle and right-hand panels: The ratio between BH mergers in the original simulation (nIllustris) and after imposing a 1 Gyr delay on all mergers [nt)], assuming growth during the specified Δt continues at the equivalent Eddington rate (middle panel) or equivalent Bondi rate (right-hand panel; see text for details). If growth is Eddington-dependent, both pre-merged black holes grow equivalently, decreasing the number of low-mass mergers with a corresponding increase in high-mass mergers (as each merger moves along a line with slope of 1 in the M2M1 plane). If growth is Bondi-dependent, then the number of low-mass mergers decreases, but most of the growth occurs in the more massive black hole, such that the main shift is in the M1-direction (i.e. slope less than one). The colour scale shows the number density of BH mergers after imposing a 1 Gyr delay relative to the number density of mergers in the original simulation; red represents an increase in the number of mergers at the given mass scales, while blue shows a decrease.

The net result of Figs 1416 is that a non-zero BH merger time-scale has the potential to significantly impact the detection rate (both at LISA and PTA sensitivities), and the chirp mass and mass ratios of detectable mergers. However, the magnitude of the changes remains sensitive to the precise accretion prescription during the binary phase, which cannot be fully addressed in a post-processing analysis. Rather, this demonstrates the need for high-resolution simulations that directly incorporate a hardening time-scale and the associated accretion on to each BH of the binary.

In Fig. 17, we plot the typical chirp masses as a function of redshift for the original simulation (red contours) and after incorporating a delay (blue contours). As in Fig. 14, we find the largest effect at the earliest times. High redshift merging BH masses are most strongly impacted by the inspiral time: combining the efficient growth of high-z BHs with a long inspiral time relative to the current BH age leads to a dramatic increase in typical chirp masses. In contrast, at low redshift we see only a relatively minor increase in the typical masses. This impact on typical chirp masses lies within LISA sensitivity (the dashed lines), and so is necessarily when considering expected LISA detections. However, we note that the low-mass end of the contours here are again limited by the simulation resolution, and that LISA should detect significant numbers of BH mergers below the minimum chirp mas probed by Illustris. We also show the evolution of chirp mass for BHs seeded assuming direct collapse BH seed formation, using the results of DeGraf & Sijacki (2020), which show an increasing Mchirp with time. This time evolution is due to the change in seed criteria: the direct collapse conditions (especially the metallicity requirement) are generally more easily satisfied at high redshift. Thus, direct collapse seeds in our simulations tend to form at high redshift and grow relatively efficiently, resulting in few low-mass BHs at low redshift (though other formation pathways may produce additional low-mass low redshift BHs). Both these additions to the BH model (incorporating a coalescence time, and modelling seed formation through a more physically motivated direct collapse model) result in increasing the expected merging BH masses. However, as we see in Fig. 17, the redshifts at which the typical masses increase is reversed: modelling formation via direct collapse leads to an increase of Mchirp with time, whereas incorporating a coalescence time tends to result in a decreasing Mchirp with time. Thus, we note that measuring the redshift evolution of the chirp mass measurement will provide an important means of differentiating between impact of seed formation versus BH merger time-scales.

The redshift–chirp mass distribution of BH mergers in the original Illustris simulation (the red contours) and after applying a 1 Gyr delay (the blue contours), with contour levels encompassing $\sim 82 {{\ \rm per\ cent}}, \sim 55 {{\ \rm per\ cent}}$, and $\sim 15 {{\ \rm per\ cent}}$ of mergers in each case. Typical chirp masses are largely unaffected by delay times at low redshift, but significantly increased in high-z mergers. This is in contrast to the expected impact from changing the seed formation model to only consider direct collapse black holes, which predicts chirp mass to increase with time (green contours; see DeGraf & Sijacki 2020).
Figure 17.

The redshift–chirp mass distribution of BH mergers in the original Illustris simulation (the red contours) and after applying a 1 Gyr delay (the blue contours), with contour levels encompassing |$\sim 82 {{\ \rm per\ cent}}, \sim 55 {{\ \rm per\ cent}}$|⁠, and |$\sim 15 {{\ \rm per\ cent}}$| of mergers in each case. Typical chirp masses are largely unaffected by delay times at low redshift, but significantly increased in high-z mergers. This is in contrast to the expected impact from changing the seed formation model to only consider direct collapse black holes, which predicts chirp mass to increase with time (green contours; see DeGraf & Sijacki 2020).

5 CONCLUSIONS

In this work, we have taken advantage of the cosmological galaxy formation simulation Illustris to investigate the mergers of SMBHs and the connection between BH mergers and the morphologies of the galaxies in which they are found. Our main results are as follows:

  • GW detections of SMBH mergers should probe a large range of BH masses, with PTAs probing the massive end of the MBHM* relation and LISA’s sensitivity probing the low-mass end of the MBHM* relation (beyond that resolved by the Illustris simulation). Most BH mergers involve the central BH of a given galaxy, with Mchirp being within ∼0.5 dex of MBH, cen.

  • By tracking several illustrative BH mergers in time, we find that for z ≲ 1, galaxies hosting massive major BH mergers frequently show morphological signatures of a recent galaxy merger. Prior to BH mergers, galaxies are generally undisturbed. After the merger, however, we find galaxies showing multiple distinct cores (i.e. in-process galaxy merger), disturbed stellar components, tidal tails, and shells of stellar light, indicative of a recent galaxy merger, though some galaxies maintain morphologies lacking direct evidence for a recent galaxy merger.

  • We then statistically characterize galaxies as having morphological evidence of merging using the Gini–M20 relation. Minor BH mergers hosts do not tend to show galaxy merger morphologies, and low-mass major BH mergers only exhibit short-lived merger morphologies at z ∼ 0.5. However, high-mass major BH mergers tend to be found in galaxies that display merger morphologies. The correlation between BH mergers and host morphology is strongest at z ∼ 0.5, and the morphological signal tends to survive for of the order of ∼500 Myr.

  • Although the majority of BH merger hosts show disturbed morphologies at z ∼ 0.5, BH merger hosts represent a minority of disturbed galaxies. Hence, we would expect a BH merger to occur in a galaxy with a disturbed host, but we would not expect most disturbed galaxies to host a recent BH merger with a detectable GW signal.

  • Other than morphology, BH merger hosts tend to have similar properties to non-merger hosts, with the exception of SF. Galaxies hosting low-mass BH mergers tend to have a slightly higher sSFR compared to a mass-matched sample, suggesting that the galaxy merger triggers a burst of SF. In contrast, galaxies hosting high-mass BH mergers have equivalent sSFR to a mass-matched sample. At low z, these high-mass galaxies tend to have low gas fractions, and thus we would not expect to see an increase in SFR as a result of a dry galaxy merger. At high z, the gas fractions are quite high, suggesting that the galaxy merger should trigger an increase in SF. However, this would also trigger increased BH accretion and associated AGN feedback, which can suppress SF, counteracting the boost that would otherwise be found.

  • The ∼500 Myr survival time for morphological evidence in a galaxy following a high-mass major BH merger should include the inspiral/hardening time of the BH binary that is, however, not modelled in our simulations. Incorporating a more realistic time for a satellite BH to fall to the galactic centre and then for the BH binary to merge would suggest that the host galaxy will have time to relax prior to the emission of GWs. As such, it may not be possible for electromagnetic follow-ups to a GW detection to be able to morphologically identify the host galaxy of the BH merger. However, this is sensitive to the time it takes before the GW signal is emitted, and thus confirmation will require new simulations that directly incorporate physically realistic BH merger delay times.

  • A simple post-processing addition of a delay time to all BH mergers has minimal impact on the global BH merger rate at low redshift. However, the delay can affect the masses of the merging BHs, resulting in higher mass BH mergers overall. In particular, the number of low mass BH mergers decreases, while the number of high mass major BH mergers increases, which results in a typical BH merger chirp mass which increases at high z (compared to the z-independent chirp mass in the original Illustris simulation).

  • Accurately incorporating the BH merger time-scale may impact the rate at which we expect to detect GWs. In particular, our analysis suggests that a longer delay time may decrease the number of LISA detections, while increasing the number of detections from PTAs. However, this is highly sensitive to the accretion efficiency of the BHs during the infall and binary phase, further emphasizing the importance of performing simulations which self-consistently incorporate BH merger times and the BH growth that occurs during those times.

In summary, we have found that SMBH mergers occur in galaxies with typical MBHM* ratios, though the morphologies are often disturbed. In particular, we find that the galaxies hosting high-mass major BH mergers show morphological evidence of a recent galaxy merger, which generally survives for ∼500 Myr. However, this survival time-scale should include the inspiral/hardening times which the Illustris simulation does not incorporate. The fact that the BH merger delay time should be on the same order suggests that by the time a GW from a BH merger is emitted, the morphological evidence in the host galaxy may have been lost. However, this is sensitive to the BH merger delay time-scale that is highly dependent on the local medium properties around the binary, emphasizing the need to run further simulations that directly incorporate a realistic model for merger time-scales, to characterize the fraction of galaxies whose morphology may survive until black hole coalescence.

In addition to the correlation to host morphology, we also showed that although a reasonable delay for merging black holes does not dramatically affect the total number of BH mergers, it has the potential to significantly impact the masses when they merge. This is particularly important for upcoming GW surveys, as the additional growth during the binary black hole phase may decrease the rate of detections from LISA, increase the detections from PTAs, and shift the distribution of merging BH masses from both towards higher masses which increase with redshift. The magnitude of these shifts depends on the accretion efficiency of both the binary black hole system as a whole and the relative efficiency of each component black hole. Hence, next generation cosmological simulations able to realistically track SMBH binary inspirals as well as accretion flows on to the binary will play a key role in making quantitatively accurate predictions for upcoming GW surveys.

ACKNOWLEDGEMENTS

We would like to thank Marta Volonteri and Alberto Sesana for useful and interesting discussions that benefited this work, and the anonymous referee whose comments improved this paper. DS would like to thank Prof. D’Anchise for all the help and support in finalizing this manuscript. CD and DS acknowledge the support from the ERC starting grant 638707 ‘Black holes and their host galaxies: co-evolution across cosmic time’ and the STFC. This research used: The Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. Simulations were run on the Harvard Odyssey and CfA/ITC clusters, the Ranger and Stampede supercomputers at the Texas Advanced Computing Center as part of XSEDE, the Kraken supercomputer at Oak Ridge National Laboratory as part of XSEDE, the CURIE supercomputer at CEA/France as part of PRACE project RA0844, and the SuperMUC computer at the Leibniz Computing Centre, Germany, as part of project pr85je. TDM acknowledges funding from NSF ACI-1614853, NSF AST-1616168, NASA ATP 19-ATP19-0084, NASA ATP 80NSSC18K101, and NASA ATP NNX17AK56G. We also acknowledge NASA ATP 80NSSC20K0519.

DATA AVAILABILITY

The data underlying this article were derived from sources in the public domain: The Illustris Project (https://www.illustris-project.org/) and The IllustrisTNG Project (https://www.tng-project.org/)

Footnotes

2

There is a possibility of having numerical issues with halo finders near galaxy mergers, which could potentially result in spurious generation of BH seeds. However, this is rare in simulations that seed based on haloes (rather than subhalos) such as Illustris, and we have explicitly confirmed that the outlying LISA-type mergers are indeed satellite BHs with low masses compared to the host galaxy and not due to spurious seeding.

REFERENCES

Abbott
B. P.
et al. ,
2016
,
Phys. Rev. Lett.
,
116
,
061102

Amaro-Seoane
P.
et al. ,
2017
,
preprint (arXiv:1702.00786)

Arzoumanian
Z.
et al. ,
2018
,
ApJS
,
235
,
37

Barausse
E.
,
Cardoso
V.
,
Pani
P.
,
2014
,
Phys. Rev. D
,
89
,
104059

Barausse
E.
,
Dvorkin
I.
,
Tremmel
M.
,
Volonteri
M.
,
Bonetti
M.
,
2020
,
ApJ
,
904
,
16

Bartlett
D. J.
,
Desmond
H.
,
Devriendt
J.
,
Ferreira
P. G.
,
Slyz
A.
,
2021
,
MNRAS
,
500
,
4639

Begelman
M. C.
,
Blandford
R. D.
,
Rees
M. J.
,
1980
,
Nature
,
287
,
307

Begelman
M. C.
,
Volonteri
M.
,
Rees
M. J.
,
2006
,
MNRAS
,
370
,
289

Bellovary
J. M.
,
Cleary
C. E.
,
Munshi
F.
,
Tremmel
M.
,
Christensen
C. R.
,
Brooks
A.
,
Quinn
T. R.
,
2019
,
MNRAS
,
482
,
2913

Blecha
L.
et al. ,
2016
,
MNRAS
,
456
,
961

Bondi
H.
,
1952
,
MNRAS
,
112
,
195

Bondi
H.
,
Hoyle
F.
,
1944
,
MNRAS
,
104
,
273

Bromm
V.
,
Loeb
A.
,
2003
,
ApJ
,
596
,
34

Costa
T.
,
Sijacki
D.
,
Trenti
M.
,
Haehnelt
M. G.
,
2014
,
MNRAS
,
439
,
2146

Davis
M.
,
Efstathiou
G.
,
Frenk
C. S.
,
White
S. D. M.
,
1985
,
ApJ
,
292
,
371

Dayal
P.
,
Rossi
E. M.
,
Shiralilou
B.
,
Piana
O.
,
Choudhury
T. R.
,
Volonteri
M.
,
2018
,
MNRAs
,
486
,
2336

DeGraf
C.
,
Sijacki
D.
,
2020
,
MNRAS
,
491
,
4973

Derdzinski
A. M.
,
D’Orazio
D.
,
Duffell
P.
,
Haiman
Z.
,
MacFadyen
A.
,
2019
,
MNRAS
,
486
,
2754

Desvignes
G.
et al. ,
2016
,
MNRAS
,
458
,
3341

Dubois
Y.
et al. ,
2014
,
MNRAS
,
444
,
1453

Faucher-Giguère
C.-A.
,
Lidz
A.
,
Zaldarriaga
M.
,
Hernquist
L.
,
2009
,
ApJ
,
703
,
1416

Feng
Y.
,
Di-Matteo
T.
,
Croft
R. A.
,
Bird
S.
,
Battaglia
N.
,
Wilkins
S.
,
2016
,
MNRAS
,
455
,
2778

Ferrarese
L.
,
2002
,
ApJ
,
578
,
90

Gebhardt
K.
et al. ,
2000
,
ApJ
,
539
,
L13

Genel
S.
et al. ,
2014
,
MNRAS
,
445
,
175

Graham
A. W.
,
Erwin
P.
,
Caon
N.
,
Trujillo
I.
,
2001
,
ApJ
,
563
,
L11

Greene
J. E.
et al. ,
2016
,
ApJ
,
826
,
L32

Gültekin
K.
et al. ,
2009
,
ApJ
,
698
,
198

Haehnelt
M. G.
,
Rees
M. J.
,
1993
,
MNRAS
,
263
,
168

Häring
N.
,
Rix
H.-W.
,
2004
,
ApJ
,
604
,
L89

Hinshaw
G.
et al. ,
2013
,
ApJS
,
208
,
19

Katz
M. L.
,
Kelley
L. Z.
,
Dosopoulou
F.
,
Berry
S.
,
Blecha
L.
,
Larson
S. L.
,
2020
,
MNRAS
,
491
,
2301

Kelley
L. Z.
,
Blecha
L.
,
Hernquist
L.
,
2017a
,
MNRAS
,
464
,
3131

Kelley
L. Z.
,
Blecha
L.
,
Hernquist
L.
,
Sesana
A.
,
Taylor
S. R.
,
2017b
,
MNRAS
,
471
,
4508

King
A.
,
2003
,
ApJ
,
596
,
L27

Klein
A.
et al. ,
2016
,
Phys. Rev. D
,
93
,
024003

Kocsis
B.
,
Yunes
N.
,
Loeb
A.
,
2011
,
Phys. Rev. D
,
84
,
024032

Kormendy
J.
,
Ho
L. C.
,
2013
,
ARA&A
,
51
,
511

Kormendy
J.
,
Richstone
D.
,
1995
,
ARA&A
,
33
,
581

Loeb
A.
,
Rasio
F. A.
,
1994
,
ApJ
,
432
,
52

Lotz
J. M.
,
Primack
J.
,
Madau
P.
,
2004
,
AJ
,
128
,
163

Magorrian
J.
et al. ,
1998
,
AJ
,
115
,
2285

Mannerkoski
M.
,
Johansson
P. H.
,
Pihajoki
P.
,
Rantala
A.
,
Naab
T.
,
2019
,
ApJ
,
887
,
35

Marinacci
F.
et al. ,
2018
,
MNRAS
,
480
,
5113

McConnell
N. J.
,
Ma
C.-P.
,
2013
,
ApJ
,
764
,
184

Naiman
J. P.
et al. ,
2018
,
MNRAS
,
477
,
1206

Nelson
D.
et al. ,
2015
,
Astron. Comput.
,
13
,
12

Nelson
D.
et al. ,
2018
,
MNRAS
,
475
,
624

Okamoto
T.
,
Frenk
C. S.
,
Jenkins
A.
,
Theuns
T.
,
2010
,
MNRAS
,
406
,
208

Oppenheimer
B. D.
,
Davé
R.
,
2008
,
MNRAS
,
387
,
577

Pillepich
A.
et al. ,
2018a
,
MNRAS
,
473
,
4077

Pillepich
A.
et al. ,
2018b
,
MNRAS
,
475
,
648

Puchwein
E.
,
Springel
V.
,
2013
,
MNRAS
,
428
,
2966

Rahmati
A.
,
Pawlik
A. H.
,
Raičević
M.
,
Schaye
J.
,
2013
,
MNRAS
,
430
,
2427

Rantala
A.
,
Pihajoki
P.
,
Johansson
P. H.
,
Naab
T.
,
Lahén
N.
,
Sawala
T.
,
2017
,
ApJ
,
840
,
53

Reardon
D. J.
et al. ,
2016
,
MNRAS
,
455
,
1751

Regan
J. A.
,
Haehnelt
M. G.
,
2009
,
MNRAS
,
396
,
343

Reines
A. E.
,
Volonteri
M.
,
2015
,
ApJ
,
813
,
82

Reines
A. E.
,
Condon
J. J.
,
Darling
J.
,
Greene
J. E.
,
2020
,
ApJ
,
888
,
36

Ricarte
A.
,
Natarajan
P.
,
2018
,
MNRAS
,
481
,
3278

Rodriguez-Gomez
V.
et al. ,
2015
,
MNRAS
,
449
,
49

Salcido
J.
,
Bower
R. G.
,
Theuns
T.
,
McAlpine
S.
,
Schaller
M.
,
Crain
R. A.
,
Schaye
J.
,
Regan
J.
,
2016
,
MNRAS
,
463
,
870

Sayeb
M.
,
Blecha
L.
,
Kelley
L. Z.
,
Gerosa
D.
,
Kesden
M.
,
Thomas
J.
,
2021
,
MNRAS
,
501
,
2531

Schaller
M.
,
Dalla Vecchia
C.
,
Schaye
J.
,
Bower
R. G.
,
Theuns
T.
,
Crain
R. A.
,
Furlong
M.
,
McCarthy
I. G.
,
2015
,
MNRAS
,
454
,
2277

Schutte
Z.
,
Reines
A. E.
,
Greene
J. E.
,
2019
,
ApJ
,
887
,
245

Sesana
A.
,
Volonteri
M.
,
Haardt
F.
,
2007
,
MNRAS
,
377
,
1711

Shankar
F.
et al. ,
2016
,
MNRAS
,
460
,
3119

Sijacki
D.
,
Vogelsberger
M.
,
Genel
S.
,
Springel
V.
,
Torrey
P.
,
Snyder
G. F.
,
Nelson
D.
,
Hernquist
L.
,
2015
,
MNRAS
,
452
,
575

Simon
J.
,
Burke-Spolaor
S.
,
2016
,
ApJ
,
826
,
11

Snyder
G. F.
,
Lotz
J.
,
Moody
C.
,
Peth
M.
,
Freeman
P.
,
Ceverino
D.
,
Primack
J.
,
Dekel
A.
,
2015
,
MNRAS
,
451
,
4290

Snyder
G. F.
,
Rodriguez-Gomez
V.
,
Lotz
J. M.
,
Torrey
P.
,
Quirk
A. C. N.
,
Hernquist
L.
,
Vogelsberger
M.
,
Freeman
P. E.
,
2019
,
MNRAS
,
486
,
3702

Springel
V.
,
2010
,
MNRAS
,
401
,
791

Springel
V.
,
Hernquist
L.
,
2003
,
MNRAS
,
339
,
289

Springel
V.
,
White
S. D. M.
,
Tormen
G.
,
Kauffmann
G.
,
2001
,
MNRAS
,
328
,
726

Springel
V.
et al. ,
2005a
,
Nature
,
435
,
629

Springel
V.
,
Di Matteo
T.
,
Hernquist
L.
,
2005b
,
ApJ
,
620
,
L79

Springel
V.
et al. ,
2018
,
MNRAS
,
475
,
676

Tremaine
S.
et al. ,
2002
,
ApJ
,
574
,
740

Tremmel
M.
,
Governato
F.
,
Volonteri
M.
,
Quinn
T. R.
,
Pontzen
A.
,
2018a
,
MNRAS
,
475
,
4967

Tremmel
M.
,
Governato
F.
,
Volonteri
M.
,
Pontzen
A.
,
Quinn
T. R.
,
2018b
,
ApJ
,
857
,
L22

Verbiest
J. P. W.
et al. ,
2016
,
MNRAS
,
458
,
1267

Vogelsberger
M.
et al. ,
2014a
,
MNRAS
,
444
,
1518

Vogelsberger
M.
et al. ,
2014b
,
Nature
,
509
,
177

Volonteri
M.
,
Natarajan
P.
,
2009
,
MNRAS
,
400
,
1911

Volonteri
M.
et al. ,
2020
,
MNRAS
,
498
,
2219

Weinberger
R.
et al. ,
2017
,
MNRAS
,
465
,
3291

Weinberger
R.
et al. ,
2018
,
MNRAS
,
479
,
4056

Wiersma
R. P. C.
,
Schaye
J.
,
Theuns
T.
,
Dalla Vecchia
C.
,
Tornatore
L.
,
2009
,
MNRAS
,
399
,
574

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.