Morphological evolution of supermassive black hole merger hosts and multimessenger signatures

With projects such as Laser Interferometer Space Antenna (LISA) and Pulsar Timing Arrays 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 that these mergers occur in typical galaxies along the $M_{\rm{BH}}-M_*$ relation, and that between LISA and PTAs we expect to probe the full range of galaxy masses. As galaxy mergers can trigger increased 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 AGN 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 survives for ~500 Myr. This is on the same scale as the infall/hardening time of the merging black holes, suggesting that electromagnetic followups to gravitational wave signals may not be able to observe this correlation. We further find that incorporating a realistic timescale delay for the black hole mergers could shift the distribution of merger masses toward higher-masses, decreasing the rate of LISA detections while increasing the rate of PTA detections.


INTRODUCTION
It is well established that supermassive black holes (SMBHs) are found at the center 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;McConnell & Ma 2013;Kormendy & Ho 2013;Reines & Volonteri 2015;Greene et al. 2016;Schutte et al. 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 et al. 1980).
Over the past several years, gravitational wave (GW) signals from BH mergers have been detected by the LIGO-Virgo collab-oration (Abbott et al. 2016), but so far these gravitational waves have been limited to those produced by mergers between stellar mass BHs. The expected mergers between SMBHs at the centers of galaxies would produce much longer wavelength gravitational waves, 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 ∼ 10 4 − 10 7 M (Amaro- Seoane et al. 2017). In addition, Pulsar Timing Array observations should be capable of detecting mergers at even higher masses, reaching BHs above 10 8 M (Verbiest et al. 2016;Desvignes et al. 2016;Reardon 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.
Gravitational wave 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;Simon & Burke-Spolaor 2016;Shankar et al. 2016), gas environment and accretion efficiencies of BHs (e.g. Kocsis et al. 2011; Barausse et al. 2014;Derdzinski et al. 2019), and even the mechanism by which BH seeds form (see, e.g. Sesana et al. 2007;Ricarte & Natarajan 2018;DeGraf & Sijacki 2020). Beyond information about the BHs themselves, multimessenger studies which 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 which incorporate both galaxy formation and SMBHs provide an excellent tool to investigate the connection between BH mergers detectable by gravitational waves and the properties and histories of the host galaxies in which they are found. Current cosmological simulations (e.g. Vogelsberger et al. 2014a;Dubois et al. 2014;Schaller et al. 2015;Feng et al. 2016;Pillepich et al. 2018b) are able to probe a wide range of spatial scales, with BHs ranging from ∼ 10 4 − 10 10 M , combined with galaxies with resolved morphological structure (e.g. Snyder et al. 2015Snyder et al. , 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 gravitational wave detections, associated electromagnetic followups, and how to interpret those observations from a theoretical framework. In this paper, we use the Illustris simulation  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 analyze how recent BH mergers correlate with galaxy properties (such as morphology and star formation rate), the timescale over which these correlations survive, and the implications this has for electromagnetic followup 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 which 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 star formation rates. In Section 4 we estimate the impact of incorporating inspiral and binary hardening timescales into the BH merger model, particularly with regards to the rate and masses of merging BHs detectable via gravitational waves. Finally, in Section 5 we summarize our conclusions.
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 collisionless sink particles, which are seeded at mass M seed = 10 5 h −1 M into any halo with M halo > 5 × 10 10 h −1 M which 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 et al. 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 which deposit energy onto the surrounding material or alter its cooling rate ("quasar", "radio", and "radiative", depending on the accretion efficiency). BHs which 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 which 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 which 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 catalog of gravitationally self-bound substructures, and the mergers between these subhalos are tracked using the SubLink catalog (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. The contours show the overall distribution of galaxies in the Illustris simulation (the outermost contour encompasses ∼ 99.8% of galaxies, while the innermost encompasses ∼ 16%), with a comparison to the local relation (Kormendy & Ho 2013, dashed grey line). The coloured points are the subsample of galaxies which have hosted a BH merger in the past 300 Myr, divided into major (M 2 /M 1 > 0.25, filled circles) vs. minor (0.1 < M 2 /M 1 < 0.25, open circles) BH mergers, and LISA-type (M chirp < 10 7 M , blue) vs PTA-type (M chirp > 10 8 M ) mergers. M BH 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 (star, upward triangle, downward triangle) correspond to specific targets investigated in more detail in Figs. 5-7.

Host galaxies
We start our analysis by looking at the galaxies which host SMBH mergers, in particular focusing on the scaling relation between BH mass (M BH ) 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 which have, in the past 300 Myr, hosted mergers at scales expected to be detectable by Pulsar Timing Arrays (PTA) or LISA. In blue we show LISA-detectable mergers, defined as BH mergers with chirp mass M chirp (M 1 M 2 ) 3/5 /(M 1 + M 2 ) 1/5 < 10 7 M (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 chirp > 10 8 M . We also subdivide the mergers based on the ratio between the larger (M 1 ) and smaller (M 2 ) masses involved in the merger, with minor BH mergers defined as BH mass ratios of 0.1 < M 2 /M 1 < 0.25, and major BH mergers defined by M 2 /M 1 > 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.
We note there are three outliers: LISA-type mergers which occur in high-mass galaxies, with M BH >> 10 8 M , well above our mass scale of M chirp < 10 7 M . These represent massive galaxies with high-mass central BHs (i.e. high M BH ) which host a merger between two low-mass satellite BHs (i.e. low BH merger masses, and thus low M chirp ). This is a rare occurrence, consisting of only ∼ 2% of LISA-type mergers within Illustris 2 . For the remainder of cases, however, we find that M chirp and M BH are similar, with a mean log(M BH /M chirp ) = 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 (M seed = 10 5 h −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 M BH /M * ratios and generally have M chirp within ∼ 0.5 dex of M BH,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 supermassive black hole 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 M 2 /M 1 > 0.25, and a massive BH merger as one in which M 2 > 10 7 M . 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. 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 panels) tend to be bluer than the higher mass galaxies (lower right 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 which 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).

Host evolution
To better understand how galaxies hosting supermassive black hole mergers will evolve with time, in Figs. 5-8 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 (M DM , M * , M gas , star formation rate, and U-V colour) in the leftmost panel, and the galaxy morphology at each snapshot in the  Nelson et al. 2015) of z = 0.5 galaxies that hosted a BH merger with M 2 > 10 7 M (thus roughly PTA-type mergers) and M 2 /M 1 > 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 scalebar 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. right hand panels. The morphology images are in sequential order starting at the upper left panel, with each panel labelled with the time relative to the BH merger.
In Fig. 5 we show the evolutionary history of a galaxy which hosts a merger between black holes with mass M 1 = 6.3 × 10 8 M and M 2 = 3.4 × 10 8 M . 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 star formation, with a corresponding change to bluer colour, near the merger time, as the galaxy merger triggers a burst of star formation (Springel et al. 2005a). This is followed by a significant drop in star formation (SF) and gas mass (M gas decreases by ∼1.5 orders of magnitude in the 500 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 × 10 9 M and M 2 = 1.2 × 10 9 M ). Here we again see a decreasing star formation rate (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 ev-  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. idence for a dual nucleus, and again there are light shells which survive long past the merger event (though fainter than in Fig. 5).
In Fig. 7 we consider the history of a galaxy hosting a less massive, minor BH merger (M 1 = 1.1 × 10 8 M and M 2 = 2.3 × 10 7 M , 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.
We also consider the impact of the black hole model used in the simulation. In particular, the Illustris simulation has previously been found to have Active Galactic Nucleus (AGN) feedback which 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. Nonetheless, to test the sensitivity to the assumptions made regarding black hole modeling in the simulation, we also show an evolutionary history from the TNG100 model of the IllustrisTNG 3 suite of simulations (Marinacci et al. 2018;Springel et al. 2018;Nelson et al. 2018;Pillepich et al. 2018a;Naiman et al. 2018), 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  Fig. 2, but at z = 2. As in Figs. 2-3 we see disturbed morphologies, but less universally so across the full sample, with some galaxies displaying unperturbed disc-like structure.
here do not exhibit a strong sensitivity to the black hole 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.

Host morphologies
Having seen the morphologies of galaxies hosting major massive BH mergers in Figs. 2-4 and individual evolutions of several sample galaxies in Figs. 5-8, 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 et al. 2004, for more details) and M 20 (defined as the secondorder moment of the brightest 20% of pixels) parameters as simple measures of morphology. In Fig. 9 we show the Gini-M 20 relation for galaxies at z = 0.5. The coloured contours show the distribution of all galaxies with M * > 10 9.7 M in Illustris. "Normal" nonmerging galaxies will tend to lie toward 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 M 20 (i.e. the brightest 20% of pixels are less centrally concentrated) compared to a nonmerging 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 datapoints 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 M 20 -binned mean Gini values (coloured errorbars) compared to the mean of the full population (black errorbars). 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-M 20 galaxies, consistent with the earlier qualitative results.
To better characterize this difference, we use the merger statistic S (see Snyder et al. 2019), defined as: (where G is the Gini coefficient) which provides a single parameter to quantify how merger-like the galaxy morphology is, with a positive S (G, M 20 ) corresponding to a merger-like galaxy and negative S (G, M 20 ) 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, M 20 ). The error bars show the mean and standard deviation of S (G, M 20 ) 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. red points calculate the time since the most recent massive major BH merger; blue points calculate the time since the most recent lowmass major BH merger, but ignore galaxies where there was a more recent high-mass major BH merger).
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, M 20 ), while low-mass major BH merger hosts show a moderately higher S (G, M 20 ), 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, M 20 ) for all galaxies (black) and for galaxies which, 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.
In addition, we show the equivalent distributions at z = 1 (middle panels) and z = 2 (right panels). At z = 1 we find that lowmass 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, M 20 ) 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 < 0.2% likelihood of being drawn from the same distribution as the full galaxy sample (Table 1). We also checked other ∆t, and the only longer timescale 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 Right panels: Stellar composite images from sequential snapshots, corresponding to the circles marking snapshots in the left 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. Here we see less disturbance in the nucleus, but nonetheless long-lived, well defined light shells. This galaxy corresponds to the upward 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.  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. 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.
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 which 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 M tot distribution (neglecting the 5% furthest outliers in each direction) and the full galaxy population's M tot 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, M 20 ) distributions are not dependent on galaxy mass, and the significant increase among BH merger hosts is genuinely due to merging.

Star formation rates
Given the potential for galaxy mergers to influence typical host star formation, 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 star formation rates. In Fig. 12, we plot the specific star formation rate (sSFR SFR/M * ) vs. total galaxy mass at z = 0.5. As expected, we see that the sSFR tends to drop among high-mass galaxies (contours and 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 M tot -binned mean and standard deviation (error bars). High-mass major BH mergers tend to have similar sSFR to other equivalent mass galaxies (red vs. black error bars). However low-mass major BH mergers tend to have slightly higher sSFR for their masses (blue vs. black error bars), suggesting that those mergers tend to trigger some additional star formation.
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 (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 equivalentmass sample whereas low-mass major BH merger hosts tend to have slightly higher sSFR than their equivalent-mass sample.
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 which host lower-mass BH mergers would be expected to be more gas-rich, and thus have more potential to drive a burst of star formation (e.g. Vogelsberger et al. 2014b;Genel et al. 2014). Indeed, we find this to be the case at z = 0.5: galaxies hosting high-mass major BH mergers tend to have M gas /M b ∼ 0.43 ± 0.22, compared to lowmass major BH merger hosts, which have M gas /M b ∼ 0.80 ± 0.16 (see right 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 which still do not appear to affect the star formation 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 ma-  Table 1. The KS-test likelihood that the target sample (high-or low-mass major BH merger hosts) comes from the same S (G, M 20 ) distribution as the full galaxy sample (see Fig. 11). At z = 0.5, low mass major BH merger hosts have a signal which 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. Figure 9. The Gini -M 20 relation at z = 0.5 for galaxies with M * 10 9.7 M (contours, from the outermost encompassing ∼ 99.7% of galaxies to the innermost encompassing ∼ 30%), for galaxies which have hosted a BH merger in the past 300 Myr with BH mass ratio M 2 /M 1 > 0.1 (red points) and M 2 /M 1 > 0.1; M 2 > 10 7 M (blue points), and error bars which 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, upward triangle, and downward triangle correspond to the individual histories plotted in Fig. 5, 6, and 7, respectively.
jor BH merger hosts). This suggests that the black holes themselves may play a role in preventing the star formation from increasing. In particular, we know that feedback from massive black holes can self-regulate (e.g. King 2003;Springel et al. 2005b;Costa et al. 2014;Sijacki et al. 2015), and can suppress SF, while the low mass black hole may not be capable of providing sufficient feedback to do so. Indeed, we find a correlation between the energy emitted by the black hole 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 black holes can help regulate star formation following a major galaxy merger, suppressing the boost in SFR one might otherwise expect from a gas-rich galaxy, as shown in Fig. 5. ) for z = 0.5 galaxies as a function of time since the most recent BH merger within the galaxy (contours, ranging from ∼ 80% of galaxies to ∼ 11%). Also shown are the mean and standard deviation of S statistics binned by time since the most recent major BH merger (M 2 /M 1 > 0.25, blue points), since the most recent massive major BH merger (M 2 /M 1 > 0.25; M 2 > 10 7 M , red points), and since the most recent non-major BH merger (M 2 /M 1 < 0.25, black points).

MERGER TIMESCALES
Having found evidence that galaxies show morphological signatures of a recent galaxy merger over timescales on the order of 300-500 Myr, we now consider the typical BH merger timescales and the impact that may have on our results. Within the Illustris simulation (and indeed many similar simulations), a pair of black holes merge as soon as their separation is less than the particle's smoothing length, rather than incorporating a coalescence time for the binary.  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 which 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 dashed red and blue 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 M tot distribution of the red and blue populations, respectively). . Specific star formation rate vs. total galaxy mass for the full galaxy sample (contours, with the outermost encompassing ∼ 99.7% of galaxies, and the innermost encompassing ∼ 46% of galaxies) and highmass (red) and low-mass (blue) major BH merger hosts. We find that highmass major BH mergers tend toward high-mass hosts and thus lower sSFR, but are fully consistent with a mass-matched sample from the full galaxy population.
both substantial and redshift-dependent (Bartlett et al. 2021). The Illustris simulation uses a re-centering scheme whereby black holes are re-positioned toward the local potential minimum; this prevents numerical wandering of black holes, but also means the the full infall time to reach the galaxy center may be notably under-estimated. Furthermore, both simulations (e.g. Bellovary et al. 2019) and ob-servations (e.g. Reines et al. 2020) suggest that black holes in dwarf galaxies may frequently be located offset from the galaxy center, which could further delay any mergers involving low-mass black holes seeded into Illustris (which are initially placed at the galaxy center). 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 supermassive black hole mergers throughout the observed universe as a function of merger redshift (solid black line), in terms of the rate at which the signals would reach Earth: 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 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.
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. In addition to the overall rate, we also show the BH merger rates if the expected inspiral times for the mergers were 500 Myr (red lines) or 1 Gyr (blue lines) following when in Illustris a given black hole 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 which a true merger timescale 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 black holes more time to grow prior to the merger, so Figure 14. 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 or 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.
Dashed/dotted lines: The BH merger rate signal for mergers corresponding to PTA-type events (dotted; defined as M chirp > 10 8 M ) and LISA-type events (dashed; defined as M chirp < 10 7 M ). Note that the rate of LISAtype 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).
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 black holes are accreting). We show the expected effect this will have with dashed and dotted lines in Fig. 14, where we consider the rate of merging black holes with M chirp < 10 7 M (dashed, roughly corresponding to LISA-detectable mergers) and M chirp > 10 8 M (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 M seed = 10 5 h −1 M , whereas the actual detection range of LISA should peak at lower-mass scales (see, e.g. Ricarte & Natarajan 2018;Dayal et al. 2018;DeGraf & Sijacki 2020); thus the LISAtype merger curve shows how incorporating a delay can impact the rate of the higher-mass LISA detections, but should nonetheless 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 black holes during that delay time is at the equivalent Eddington fraction as the post-merged black hole in the original simulation (i.e. we check the fractional accretion in Illustris of the merged black hole over the course of ∆t, and let each progenitor black hole grow by the same fraction). Here we see that in addition to decreasing the number of high-z BH mergers (especially for lowmass, LISA-type mergers), longer delay times tend to result in more high-mass BH mergers at lower-z as the inspiralling black holes have more time to grow prior to merging. In particular, we find that longer merger timescales which include the associated black hole 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 (M chirp = (M 1 M 2 ) 3/5 /(M 1 + M 2 ) 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 black hole from the original simulation.
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 toward a higher efficiency than the lower mass from the pre-merged black hole 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 black hole 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 chirp ∼ 2 × 10 6 M in the original simulation, ∼ 7 × 10 6 M for a 500 Myr delay, and ∼ 10 7 M 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 black holes in the Illustris simulation, and would thus expect the actual detection rate for LISA to be significantly higher (see, e.g. Ricarte & Natarajan 2018;Dayal et al. 2018;DeGraf & Sijacki 2020).
We show the distribution of BH merger masses in more detail in Fig. 16: in the left panel we plot the distribution of the primary (M 1 ) and secondary (M 2 ) merger masses in Illustris, finding that mergers peak at M 1 ∼ 5 × 10 6 M and M 2 ∼ 10 6 M . We also show the fractional change in merger masses when including a delay (n(∆t = 1 Gyr)/n(∆t = 0)), assuming growth at a fixed Eddington fraction (middle panel) or fixed Bondi fraction (right 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 black holes 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 2 BH , 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 black hole 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).
The net result of Figs. 14-16 is that a non-zero BH merger timescale 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 which directly incorporate a hardening timescale and the associated accretion onto each black hole 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 black holes with a long inspiral time relative to the current black hole 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 (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 black holes seeded assuming direct collapse black hole (DCBH) seed formation, using the results of DeGraf & Sijacki (2020), which show an increasing M chirp 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 black holes 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 modeling 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: modeling formation via direct collapse leads to an increase of M chirp with time, whereas incorporating a coalescence time tends to result in a decreasing M chirp with time. Thus, we note that measuring the redshift evolution of The ratio between BH mergers in the original simulation (n Illustris ) 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 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 M 2 − M 1 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 M 1 -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 17. The redshift -chirp mass distribution of BH mergers in the original Illustris simulation (red contours) and after applying a 1 Gyr delay (blue contours), with contour levels encompassing ∼ 82%, ∼ 55%, and ∼ 15% 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). the chirp mass measurement will provide an important means of differentiating between impact of seed formation vs. BH merger timescales.

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.
• Gravitational wave detections of supermassive black hole mergers should probe a large range of BH masses, with PTAs probing the massive end of the M BH − M * relation and LISA's sensitivity probing the low mass end of the M BH − M * relation (beyond that resolved by the Illustris simulation). Most BH mergers involve the central black hole of a given galaxy, with M chirp being within ∼0.5 dex of M BH,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-M 20 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 which 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 gravitational wave signal.
• Other than morphology, BH merger hosts tend to have similar properties to non-merger hosts, with the exception of star formation. 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 star formation. 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 which is however not modelled in our simulations. Incorporating a more realistic time for a satellite black hole to fall to the galactic center and then for the BH binary to merge would suggest that the host galaxy will have time to relax prior to the emission of gravitational waves. As such, it may not be possible for electromagnetic followups 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 which 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 black holes, 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 timescale may impact the rate at which we expect to detect gravitational waves. 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 black holes during the infall and binary phase, further emphasizing the importance of performing simulations which self-consistently incorporate BH merger times and the BH growth which occurs during those times.
In summary, we have found that supermassive black hole mergers occur in galaxies with typical M BH − M * 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 timescale 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 gravitational wave 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 timescale which is highly dependent on the local medium properties around the binary, emphasizing the need to run further simulations which directly incorporate a realistic model for merger timescales, 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 toward 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 supermassive black hole binary inspirals as well as accretion flows onto the binary will play a key role in making quantitatively accurate predictions for upcoming gravitational wave surveys.

ACKNOWLEDGMENTS
We would like to thank Marta Volonteri and Alberto Sesana for useful and interesting discussions which 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 Su-perMUC 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.illustrisproject.org/) and The IllustrisTNG Project (https://www.tngproject.org/)