Black hole mergers from dwarf to massive galaxies with the NewHorizon and Horizon-AGN simulations

Horizon-AGN simulates a large volume, (142 comoving Mpc)³, at a relatively low spatial and mass resolution: refinement is permitted down to $\Delta x=1$ kpc, the dark matter particle mass is $8\times 10^{7}\,{\rm M_{\odot}}$, the stellar particle mass is $2\times 10^{6}\,{\rm M_{\odot}}$, and the MBH seed mass is $10^{5}\,{\rm M_{\odot}}$.

Gas cooling is modelled down to $10^{4}\,\rm K$ with curves from Sutherland & Dopita (1993). Heating from a uniform UV background takes place after redshift $z_{\rm reion}=10$ following Haardt & Madau (1996). The gas follows an equation of state for an ideal monoatomic gas with an adiabatic index of $\gamma_{\rm ad}=5/3$.

Star formation is modelled with a Schmidt relation adopting a constant star formation efficiency $\epsilon_{*}=0.02$ (Kennicutt, 1998; Krumholz & Tan, 2007) in regions which exceed a gas hydrogen number density threshold of $n_{0}=0.1\,\rm H\,cm^{-3}$ following a Poisson random process (Rasera & Teyssier, 2006; Dubois & Teyssier, 2008). Mechanical energy injection from Type Ia SNe, Type II SNe and stellar winds is included assuming a Salpeter (1955) initial mass function with cutoffs at $0.1\,\rm M_{\odot}$ and $100\,\rm M_{\odot}$.

MBHs are created in cells where the gas density is larger than $n_{0}$, and where the gas velocity dispersion is larger than $100\,\rm km\,s^{-1}$; an exclusion radius of 50 comoving kpc is imposed to avoid formation of multiple MBHs in the same galaxy. MBH formation is stopped at $z=1.5$. The accretion rate follows a Bondi-Hoyle-Littleton rate modified by a factor $\alpha=(n/n_{0})^{2}$ when $n>n_{0}$ and $\alpha=1$ otherwise (Booth & Schaye, 2009) in order to account for the inability to capture the multiphase nature of the interstellar gas. The effective accretion rate onto MBHs is capped at the Eddington luminosity with a radiative efficiency of 0.1. For luminosities above 1% of the Eddington luminosity, 15% of the MBH emitted luminosity is isotropically coupled to the gas within $4\Delta x$ as thermal energy.

At lower luminosities feedback takes a mechanical form, with 100% of the power injected into a bipolar outflow with a velocity of $10^{4}\,\rm km\,s^{-1}$, injected in a cylinder with radius $\Delta x$ and height $2\,\Delta x$.

To avoid spurious motions of MBHs due to finite force resolution effects, we adopt an explicit drag force of the gas onto the MBH (Dubois et al., 2012) with the same boost factor $\alpha$ used for accretion. This gas dynamical friction is expressed as $F_{\rm DF}=f_{\rm gas}4\pi\alpha\rho_{\rm gas}(GM_{\rm BH}/\bar{c}_{s})^{2}$, where $\rho_{\rm gas}$ is the mass-weighted mean gas density within a sphere of radius $4\,\Delta x$ and $f_{\rm gas}$ is a factor function of the mach number ${\mathcal{M}}=\bar{u}/\bar{c}_{s}$, which accounts for the extension and shape of the wake (Ostriker, 1999) and takes a value between 0 and 2 for an assumed Coulomb logarithm of 3 (Chapon et al., 2013). Including this drag force allows us to avoid pinning MBHs to galaxy centers, i.e. constantly repositioning them at the minimum of the local potential, which causes unnatural dynamics (Tremmel et al., 2015). See Dubois et al. (2013) for additional details.

NewHorizon is a smaller and higher resolution volume: it resimulates an “average” density sphere with a radius 10 comoving Mpc of Horizon-AGN, focusing on field galaxies. The most massive halo at $z=0.45$ is $5.9\times 10^{12}\,{\rm M_{\odot}}$. The dark matter mass resolution in the zoomed region is $M_{\rm DM,hr}=1.2\times 10^{6}\,\rm M_{\odot}$, stellar mass resolution is $10^{4}\,\rm M_{\odot}$, $\Delta x=40$ pc, and the MBH seed mass $10^{4}\,{\rm M_{\odot}}$. In the high-resolution region a passive refinement scalar is injected, with a value of $1$ within that region, and we only allow for refinement when the value of this scalar is above³³3The initial pure Lagrangian volume changes size and shape over time and the Eulerian volume gets polluted by low-resolution dark matter particles over time. a value of $0.01$. We use this passive scalar also to measure the refined volume, which collapses somewhat from the initial volume of $\sim 4000$ Mpc³ under the effect of gravity.

Gas cooling is modelled with rates tabulated by Sutherland & Dopita (1993) above $10^{4}\,\rm K$ and those from Dalgarno & McCray (1972) below $10^{4}\,\rm K$. The uniform UV background is as in Horizon-AGN, except that we include self-shielding where the gas density is larger than $n_{\rm shield}=0.01\,\rm H\,cm^{-3}$, i.e. UV photo-heating is reduced by a factor $\exp(-n/n_{\rm shield})$.

Star formation follows a Schmidt law, but with a density threshold of $n_{0}=10\,\rm H\,cm^{-3}$, and with a varying star formation efficiency that depends on the local turbulent Mach number and Jeans length (see e.g. Trebitsch et al., 2018). The initial mass function differs from Horizon-AGN, adopting a Chabrier functional form (Chabrier, 2005), with cutoffs at 0.1 and 150 $\,{\rm M_{\odot}}$. We include the mechanical SN feedback model from Kimm et al. (2015), which models the energy and momentum conserving phases of the explosion separately. This SN feedback is more effective than the purely kinetic formulation used in Horizon-AGN (Kimm et al., 2015).

Black holes form in cells where both the gas and stellar densities are above the threshold for star formation and the stellar velocity dispersion is larger than $20\,\rm km\,s^{-1}$, again with an exclusion radius of 50 comoving kpc. An explicit drag force of the gas onto MBHs is implemented, with the same boost as in Horizon-AGN. Since NewHorizon has sufficient spatial resolution to model some of the multiphase nature of the gas, accretion is modelled via a non-boosted Bondi-Hoyle-Littleton rate ($\alpha=1$), capped at the Eddington luminosity.

Spin evolution via gas accretion and MBH-MBH mergers are followed explicitly in NewHorizon, using the implementation described in Dubois et al. (2014a). For MBH-MBH mergers we adopt for the final spin a fit motivated by general-relativistic simulations (Rezzolla et al., 2008). For accretion, the direction of the angular momentum of the accreted gas is used to decide whether the accreted gas feeds an aligned or misaligned disc (King et al., 2005), respectively spinning up or down the MBH (Bardeen, 1970) for accretion above 1% of the Eddington luminosity. At lower luminosities the MBH spin up (down) rate is motivated by jets tapping energy from MBH spins (Blandford & Znajek, 1977; Moderski & Sikora, 1996) and follows the results from McKinney et al. (2012), where we fit a fourth-order polynomial to their sampled values (from their table 7, A$a$N100 runs, where $a$ is the value of the MBH spin). The radiative efficiency is calculated individually for each MBH based on its spin, thus the Eddington mass accretion rate will vary as a MBH spin evolves, and it is reduced linearly with the Eddington ratio for MBHs below 1% of the Eddington luminosity following Benson & Babul (2009). AGN feedback also depends on spin, being 15% of the spin-dependent emitted power for the thermal feedback, injected within a sphere of radius $\Delta x$, and following the results of magnetically chocked accretion discs of McKinney et al. (2012, from their table 5 for A$a$N100 runs). The bi-polar outflows are modelled as in Horizon-AGN. We also enforce the refinement within a region of radius $4\Delta x$ around the MBH at the maximum allowed level of refinement.

In Horizon-AGN we identify dark matter haloes and galaxies with the AdaptaHOP halo finder (Aubert et al., 2004). The density field used in AdaptaHOP is smoothed over 20 particles. In Horizon-AGN we fix the density threshold at 178 times the average total matter density and require 50 dark matter particles for identification for both dark matter haloes and 50 star particles for galaxies. In NewHorizon only haloes with average density larger than 80 times the critical density are considered, and the minimum number of particles per halo is 100. For galaxies in NewHorizon, we select them using the HOP finder, requiring at least 50 star particles for a galaxy. The difference with AdaptaHOP is that substructures are not extracted from the main component. This is necessary as the high resolution achieved in the NewHorizon galaxies allows the formation of dense star forming clumps, which will otherwise be removed with AdaptaHOP.

NewHorizon is a zoom simulation embedded in a larger cosmological volume filled with lower resolution dark matter particles, we need to consider only “pure” haloes, devoid of low resolution dark matter particles, and the galaxies embedded in these haloes. Low-resolution particles are more massive than high-resolution particles, therefore the presence of low- and high-resolution particles in the same halo would induce numerical issues; using only pure galaxies/halos the dynamics of the galaxies under examination is robust. The purity criterion only affects halos at the edge of the refined volume, where lower-resolution particles from outside can get mixed into the high-resolution region. At the final redshift of $z=0.45$, this gives 762 uncontaminated galaxies with stellar mass larger than $5\times 10^{5}\,\rm M_{\odot}$ and 238, 87 and 17 with masses larger than $10^{8}\,{\rm M_{\odot}}$, $10^{9}\,{\rm M_{\odot}}$ and $10^{10}\,{\rm M_{\odot}}$ respectively.

In Horizon-AGN, we obtain the centre of haloes and galaxies using the shrinking sphere approach proposed by Power et al. (2003) to get the correct halo centre. In NewHorizon we use the shrinking sphere approach for haloes and for galaxies with a reduction factor of the radius of respectively 10 and 30 per cent at each shrinking step, and with a stopping search radius of respectively 0.5 and 1 kpc. This rather large value for galaxies (with respect to their typical size) in NewHorizon allows us to avoid converging towards dense and massive but randomly located clumps, and provides a qualitatively better centring over the total stellar distribution. Despite this, defining the galaxy “centre” for low-mass and high-z galaxies remains challenging.

Using the information in the simulation and post-processing techniques we create catalogues of MBH mergers. To select merging MBHs, we use the information on all MBHs at each synchronised (coarse) timestep of the corresponding simulation, every $\sim 0.5$ Myr in NewHorizon and $\sim 0.6-0.7$ Myr in Horizon-AGN. Only MBH information is available at each coarse timestep, as galaxy and halo information is only available at full outputs which are saved much less frequently.

A pair of MBHs are merged in the simulation when a MBH present at a given time disappears in the succeeding time step; the companion MBH is located by searching for the MBH that was closest to the vanished one, keeping in mind that the simulation merges MBHs when they are separated by $\leqslant 4\Delta x$, corresponding to 4 kpc for Horizon-AGN and 160 pc for NewHorizon, and they are energetically bound in vacuum. These separations, especially for Horizon-AGN, are much larger than those where MBHs effectively merge in normal conditions and in section 3.2 we discuss how we include additional timescales to account for this. We refer to these as “numerical mergers” and the initial lists contains 542 MBH mergers for NewHorizon and 85397 for Horizon-AGN. After a numerical merger, the remnant acquires the ID of the most massive MBH in the pair. This is the starting general catalogue for Horizon-AGN (“all”) and we use the subscript “${\rm in}$” for quantities (redshift, cosmic time, MBH masses, accretion rates, galaxy masses) at the initial time of the numerical merger.

For NewHorizon, which is a zoom simulation, this initial list also includes MBHs that are located in polluted haloes, i.e., haloes that contain low-resolution dark matter particles and whose evolution is therefore untrustworthy. The first task is therefore to remove MBHs in polluted haloes from our sample. This is performed by spatially matching MBHs with galaxies in unpolluted haloes. For each MBH pair we request that the primary is located within $\max(10R_{\rm eff},4\Delta x)$ from a galaxy in an unpolluted halo, where $R_{\rm eff}$ is the geometric mean of the projected half-mass radius over the 3 cartesian axes. A galaxy is associated to a halo if the galaxy centre is within 10% of the virial radius of the dark matter halo. The initial list of 542 numerical mergers is thus reduced to 385. For NewHorizon we consider this to be the starting general catalogue (“all”).

In both NewHorizon and Horizon-AGN, MBHs are not pinned in the centre of galaxies and haloes, so an additional concern is removing possible “spurious” numerical mergers. In a realistic situation it is unlikely that MBHs merge far from the centre of the host galaxy (Merritt et al., 2001). This is because MBHs have a small impact parameter for binding to another MBH: they must pass within each other’s sphere of influence. In the simulation, for slowly moving MBHs the impact parameter is instead $4\Delta x$, much larger than the physical impact parameter under realistic conditions. The energy condition ensures that in vacuum the MBHs would bind, but not that they would merge, since much energy has to be extracted from the binary before they are sufficiently close that gravitational waves become effective in driving the MBHs to coalescence. Processes that extract energy from the orbit of a MBH or from a binary, first dynamical friction, then hardening by scattering off single stars and torques from circumbinary disc, are far more effective in higher surrounding stellar and/or gas densities. For off-centre MBHs and binaries, the orbital evolution is therefore slower than for centrally located MBHs. This is particularly important for NewHorizon, where the low initial MBH seed mass can make dynamics erratic (Pfister et al., 2019), while for Horizon-AGN the main worry, as we will see below, is the correction for dynamical delays, even in the case of “central" mergers.

To identify possible spurious numerical mergers we cross-correlate MBHs with galaxies, and we select pairs where the MBHs are within $\max(2R_{\rm eff},4\Delta x)$ (“sel”). In NewHorizon galaxy properties are output every $\sim 15$ Myr, while in Horizon-AGN every $\sim 150$ Myr. In NewHorizon we consider the galaxy output closest in time to the numerical merger. Given the sparsity of Horizon-AGN outputs, galaxies and MBHs can have moved considerably between the two outputs surrounding the time of the MBH merger; for this simulation we apply the criterion $\max(2R_{\rm eff},4\Delta x)$ on the primary MBH at its position in the galaxy outputs before and after the merger.

In NewHorizon, this catalogue contains 314 MBHs, and in Horizon-AGN it contains 64505 MBHs. Numerical noise in MBH dynamics arises when the ratio of MBH mass to stellar particle mass is too small, which results in a lack of dynamical friction from particles (Tremmel et al., 2015; Pfister et al., 2019) that would otherwise offset this noise. This spurious noise can keep MBHs away from galaxy centres, and, hence, we define “centre” very generously and we consider sel as the reference sample. We discuss alternative choices in Appendix C, noting that choices can make a significant difference – up to an order of magnitude – in the final results. The merger rate is still likely to be considered a lower limit. We apply the same selection criterion of $\max(2R_{\rm eff},4\Delta x)$ in matching MBHs to galaxies at all subsequent steps of the post-processing described below.

As discussed, the simulations numerically “merge” MBHs when their separation is much larger than the typical separation, of order of a millipc, where emission of gravitational waves becomes effective and brings the binary to coalescence in less than a Hubble time (Colpi, 2014). The catalogues have to be post-processed to obtain more realistic merging timescales and merger rates. Various approaches for adding merger timescales in semi-analytical models (Volonteri et al., 2003; Barausse, 2012; Bonetti et al., 2019) and in numerical simulations by post-processing have been proposed (Kelley et al., 2017; Krolik et al., 2019; Katz et al., 2020; Sayeb et al., 2020). Any correction will necessarily be more arbitrary for low-resolution simulations, since the available information must be extrapolated over a much larger parameter space.

We construct the history of all galaxies, from their birth to the time the simulation ends ($z=0$ for HAGN and $z=0.45$ for NewHorizon), and associate numerical MBH mergers to the galaxy mergers from which they originate.

The history of each galaxy is contained in the simulation’s merger trees: we use TreeMaker (Tweed et al., 2009) on all galaxies for Horizon-AGN and at the final redshift for NewHorizon. We define the main descendant of a galaxy as the one that shares the most mass in the following output. Two galaxies are defined as merged when their main descendant is the same, and this also defines the time of the galaxy merger.

To associate MBH and galaxy mergers, for each MBH merger we select the first output where each of the MBHs in the binary can be associated to a galaxy using a threshold of $2R_{\rm eff}$ to search for a MBH in a galaxy. We then search in the merger tree to identify a galaxy merger that involves descendants of the initial galaxies, and check that the two original MBHs were hosted in the merging galaxies. The match MBH-galaxy at the time of the galaxy merger uses a threshold of $2R_{\rm eff}$ for Horizon-AGN, but we allow the search out to $10R_{\rm eff}$ of galaxies in unpolluted haloes in NewHorizon, with in practice 90% of matches having distances less than $4.25R_{\rm eff}$. The reason is that, close to a merger, the galaxy centre is not a well-defined quantity because of morphological disturbances, and this is more evident in high-resolution simulations. With this procedure, for most numerical MBH mergers we identify a galaxy merger for which we know the redshift and the properties of the merging galaxies.

Horizon-AGN and NewHorizon are simulations with very different characteristics: Horizon-AGN simulates a large volume at low spatial and mass resolution, while NewHorizon simulates a relatively small average volume at high spatial and mass resolution. Horizon-AGN is therefore unable to correctly resolve the low-mass galaxies hosting the low-mass MBHs that LISA can detect: the sweet spot for LISA’s sensitivity is at about $10^{5}-10^{6}\,{\rm M_{\odot}}$. Conversely, NewHorizon does not contain any massive galaxies hosting MBHs massive enough to contribute to the PTA signal, which is dominated by MBHs with mass $>10^{8}\,{\rm M_{\odot}}$ (Sesana et al., 2008). This is exemplified in Fig. 2 where we show the distribution of the MBH masses in the two simulations. We note that $M_{2}$ is not followed self-consistently after the numerical merger: at any time after the numerical merger, $M_{\rm BH}$ is the sum of $M_{1}$ and $M_{2}$ plus any mass accreted after that time.

In Horizon-AGN the distribution of masses in the full (all) sample, shown as an orange dash-dotted curve, is relatively flat from the seed mass out to $10^{7}\,{\rm M_{\odot}}$, and drops afterwards. When we remove spurious merger events, a large fraction of low-mass MBH mergers disappear (sel, red dashed cuve). When we look at $t_{\rm in}+t_{\rm df}$ (dark red dotted curve), MBHs with mass $<10^{6}\,{\rm M_{\odot}}$ virtually disappear, while at the high mass end there are two competing effects: on the one hand some of largest MBHs are removed from the sample, because the time between $t_{\rm in}$ and $z=0$ is too short, on the other hand MBHs grow during $t_{\rm df}$, extending the distribution to larger masses. Finally, we look at the mass distribution at $t_{\rm in}+t_{\rm df}+t_{\rm bin}$, where we find a similar behaviour. In NewHorizon we find a more limited loss of MBHs at the low-mass end because fewer spurious MBH mergers occur in the outskirts of galaxies, but the overall trends are similar.

Fig. 3 shows the relation between the masses of MBH binaries and the total stellar mass of their host galaxies at different times: the numerical merger, $t_{\rm in}$, after the dynamical friction phase, $t_{\rm in}+t_{\rm df}$, and after the binary evolution, $t_{\rm in}+t_{\rm df}+t_{\rm bin}$. The observational samples are from Reines & Volonteri (2015) and Baron & Ménard (2019). Reines & Volonteri (2015) include both quiescent galaxies with dynamically measured MBH masses, typically hosted in spheroids, and active galaxies with single-epoch MBH masses estimates using broad AGN lines (type I AGN), hosted in both spheroids and disc galaxies. Baron & Ménard (2019) propose a novel way to estimate MBH masses from narrow AGN lines (type II AGN) and apply it to AGN at $z<0.3$.

In Horizon-AGN at given galaxy mass the binary masses are consistent with the masses of single MBHs in the simulation, and show the same trends as single MBHs, i.e., they reproduce well the data at the high-mass end but they overpredict MBH masses at the low mass end, and the distribution has a smaller scatter than in observations (see Volonteri et al. (2016) for a complete discussion). In NewHorizon, MBH growth is delayed by SN feedback (Dubois et al., 2015) and MBH growth picks up only in galaxies more massive than $5\times 10^{9}\,{\rm M_{\odot}}$ (Dekel et al., 2019) when major gas inflows allow MBHs to grow significantly above their seed masses. MBHs in the most massive galaxies by the end of the simulation, $z=0.45$ have masses in agreement with observational samples, but there are no MBHs with mass $>10^{6}$ in galaxies with stellar mass $10^{9}-10^{10}\,{\rm M_{\odot}}$. We note that this delayed MBH growth leads to a better match of the theoretical AGN luminosity function with observations in NewHorizon (Beckmann et al. in prep.) than in Horizon-AGN (Volonteri et al., 2016), where it is overestimated.

The typical mass ratio of merging binaries is an important information in preparation for LISA: at the time of writing, waveforms are available for mass ratios $q\sim 1$ and $q\ll 1$, but the range $q=10^{-1}-10^{-6}$ falls in between those that can be studied with numerical relativity, post-Newtonian techniques and gravitational self-force (Barack et al., 2019). Understanding if these mass ratios are statistically significant is therefore of great interest to the waveform community. Unfortunately, the mass ratio of the merging binaries is well defined only at $t_{\rm in}$: at that point in the simulation $M_{1}$ and $M_{2}$ are merged and therefore subsequently the two masses cannot be followed separately although we artificially consider the two MBHs as still under dynamical evolution. We can however follow the total mass evolution of the binary, by tracking the numerically merged black hole.

For reference we define the mass ratios at $t>t_{\rm in}$ in three different ways: either we consider that the mass ratio has not changed, i.e., that $M_{2}$ has grown at the same pace as $M_{1}$, or we consider that $M_{2}$ has not changed and all mass growth occurred on $M_{1}$, or that there has been a coup (Van Wassenhove et al., 2014) with a swap between $M_{1}$ and $M_{2}$ so that all subsequent growth occurred on $M_{2}$. The truth could be anywhere in between and also outside this range, but we refrain from trying to speculate further, although we note that during both the dynamical friction phase and the evolution in a circumbinary disc $M_{2}$ is expected to grow faster than $M_{1}$ (Callegari et al., 2011; Noble et al., 2012; Capelo et al., 2015; Farris et al., 2015).

The mass ratios are shown in Fig. 4 for the full population at $t_{\rm in}$ (orange dot-dashed) and for the MBHs that form a binary at $t_{\rm in}+t_{\rm df}$. Adding dynamical delays shifts the Horizon-AGN distribution to larger mass ratios, while it has the opposite effect in NewHorizon. In Horizon-AGN this is due to the long dynamical friction time when $M_{2}$ is small, while in NewHorizon, where a large fraction of MBHs can bind, the effect is due to the growth of $M_{\mathrm{BH}}$ with respect to $M_{2}$ (thick lines) and to mergers with initial $q\approx 1$ involving two MBHs with $M_{\mathrm{BH}}\approx 10^{4}\,{\rm M_{\odot}}$, which have long merging timescales. In summary, a large fraction of the binaries should have mass ratio between 0.1 and 1, but a tail at $q<0.1$ cannot be excluded.

The MBH merger rate from the two simulations, defined as the rate measurable from an observer on Earth over the whole sky (Haehnelt, 1994) is shown in Fig. 5. First, we show how removing “spurious” mergers affects the results. We compare the numerical MBH merger rates for catalogues all and sel: in Horizon-AGN almost 80% of mergers occur outside $2R_{\rm eff}$. This is not the case for NewHorizon, although in both simulations most mergers within $2R_{\rm eff}$ actually occur between half and twice the galaxy effective radius.

We stress that NewHorizon provides only a lower limit to the merger rate, since it simulates an average region of the Universe and biased regions provide a significant contribution to the merger rate even taking into account for their rarity (Sesana et al., 2005). Although the simulations should not be combined since they are not self-similar in sub-grid physics, the merger rate of MBHs for LISA should be higher than the sum of the rates of MBHs with mass $<10^{7}\,{\rm M_{\odot}}$ from NewHorizon and Horizon-AGN: NewHorizon does not have biased regions, and Horizon-AGN does not resolve low-mass galaxies. Furthermore, as noted in section 3.1, the lack of a direct implementation of dynamical friction from stars and dark matter in the simulations – only dynamical friction from gas is included on-the-fly – also goes in the direction of reducing the number of MBH mergers, since with additional friction the MBHs would have a smoother dynamics, bind more easily, and numerical mergers would be facilitated.

Including dynamical delays in post-processing can severely reduce the raw merger rate from numerical mergers (Katz et al., 2020). Adding delays shifts the merger rate to lower redshift, with a peak at $z\sim 1-2$. This is favorable for electromagnetic counterpart searches, since the counterparts will be brighter and with better sky localization from LISA (McGee et al., 2020). We postpone a complete study of the counterparts to a future study.

The merger rate is generally lower than that predicted by semi-analytical models, mostly because only few semi-analytical models include dynamical modelling of MBH orbital decay and binary evolution (Volonteri et al., 2003; Barausse, 2012; Bonetti et al., 2019), and even in these cases the early evolution is normally approximated by integrated dynamical friction timescales such as Eq. 1 that do not capture the more erratic and stochastic behaviour seen in simulations (Pfister et al., 2019; Bortolas et al., 2020). Our study includes integrated dynamical friction timescales in post-processing, but before the numerical merger dynamics is calculated directly in the simulation: MBHs respond to inhomogeneous, time-varying conditions and the direct implementation of dynamical friction, although only from gas, allows us to account for evolution of the MBH mass and of the environment on-the-fly. Recently Barausse et al. (2020) include in a semi-analytical additional kpc-scale delays before formation of binaries, and show that the merger rate is then greatly reduced.

Furthermore, the highest merger rate in semi-analytical models is predicted when including the mergers of the remnants of Population III stars in very high redshift galaxies (Sesana et al., 2007; Klein et al., 2016; Ricarte & Natarajan, 2018; Dayal et al., 2019), which we do not treat in our simulations (see Katz et al., 2020, for a careful comparison between simulations and semi-analytical models).

The predicted merger rate from Horizon-AGN is also somewhat lower than in EAGLE or Illustris (Salcido et al., 2016; Katz et al., 2020) because in our simulations MBHs are not constantly repositioned at the center of the potential minimum, which makes MBHs merge immediately after galaxy mergers: even adding delays in postprocessing simulations with MBH repositioning facilitate mergers. In both Horizon-AGN and NewHorizon the large scale dynamics is driven by the explicit use of a drag force, and we include delays in post-processing only below resolution.

The most important point in the comparison between Horizon-AGN and NewHorizon, however, is that a small, high-resolution simulation like NewHorizon predicts higher merger rates than a large, low-resolution simulation like Horizon-AGN. The reason is the ability to track mergers of low-mass galaxies hosting MBHs. Although the occupation fraction of MBHs in galaxies is predicted to decrease with galaxy mass (Volonteri et al., 2008a; van Wassenhove et al., 2010), observationally it seems to be between 0.1 and 1 in galaxies with mass $\sim 10^{9}\,{\rm M_{\odot}}$ at $z=0$ (Miller et al., 2015; She et al., 2017) and between 0.5 and 1 when considering local galaxies within 4 Mpc and mass $\sim 10^{9}-10^{10}\,{\rm M_{\odot}}$ and published MBH dynamical masses or limits (Nguyen et al., 2018). We refer to Greene et al. (2019) for a review. At $z=0.45$ the occupation fraction of MBHs within $2R_{\rm eff}$ in NewHorizon is 0.1 at $M_{\rm gal}=10^{6}\,{\rm M_{\odot}}$ and it reaches unity at $M_{\rm gal}=10^{9}\,{\rm M_{\odot}}$, while it becomes about 2, i.e., there are on average two MBHs within a galaxy, at $M_{\rm gal}=10^{10.5}\,{\rm M_{\odot}}$. Tracing MBHs in dwarf galaxies is crucial.

In the context of LISA’s science, simulations like Horizon-AGN, EAGLE or Illustris, which do not resolve dwarf galaxies, will underestimate the merger rate of LISA’s MBHs. A simulation like Romulus (Tremmel et al., 2017), with volume and resolution intermediate between Horizon-AGN and NewHorizon, and also seeding MBHs in low-mass galaxies like NewHorizon, is a good compromise. Ideally, we would like a simulation with the resolution of NewHorizon and the volume of Horizon-AGN run to $z=0$: unfortunately this is currently computationally unfeasible.

Beyond semi-analytical models and simulations, another technique often used to estimate the MBH merger rate is through the galaxy merger rate, assuming scaling relations between galaxies and MBHs (Sesana, 2013; Simon & Burke-Spolaor, 2016; Chen et al., 2019). We test this approach here, by relating MBH mergers to the galaxy merger they originated from. In our general picture, at the beginning of a galaxy merger, each MBH sits at the centre of the merging galaxy, and for mass ratios of the merging galaxies $>1:4$ we expect that the two MBHs end up in the centre of the merger in a bound binary (Begelman et al., 1980; Callegari et al., 2009; Van Wassenhove et al., 2012; Pfister et al., 2017). In the case of low-mass MBHs, however, this picture may not capture the full behaviour. Low-mass MBHs may not reside in the centre of their host galaxy (Bellovary et al., 2019) and/or their dynamical evolution can be subject to disturbances (Pfister et al., 2019) preventing the formation of the binary even in mergers between similar mass galaxies.

We compare the whole distribution of galaxy stellar masses and mass ratios to that of galaxy mergers that lead to MBH mergers in Fig. 6. Additional parameters play a role, e.g., the compactness of the satellite galaxy (Tremmel et al., 2018), but these quantities are hard to measure in observations, therefore we focus here only on masses and mass ratios, since these are “macroscopic” quantities than can be measured in observations, although even measuring a galaxy stellar mass mass requires some modelling from imaging and spectroscopy.

We consider here “numerical MBH mergers”, while the results when we include dynamical friction are discussed further down. We define $M_{\rm gal,1}$ as the most massive of the two merging galaxies, but in $\sim$10% of cases the least massive MBH is hosted in the most massive galaxy of the pair (violet histograms in Fig. 6). Each histogram is normalised to the total number of objects in the respective catalogue, and in Horizon-AGN the galaxy mergers identified as sources of a MBH merger represent 16% of the total number of mergers for catalogue sel, and 37% of catalogue all; the fractions for NewHorizon are even smaller (3.5% and 4% respectively): not all galaxy mergers lead to a MBH merger. The reason is twofold: the occupation fraction of MBHs is below unity in low-mass galaxies (Volonteri et al., 2008b), so the probability of a merger between two galaxies each hosting a MBH is also much lower than unity (as a first approximation this probability scales as the occupation fraction squared, although biased regions experience an enhanced number of mergers). Furthermore, even if galaxies merge, the MBHs can be stalled at large separations in the galaxy even for major mergers (Tremmel et al., 2018).

As expected, galaxy mergers that generate MBH mergers are a biased sample with respect to the general merging population, specifically they have a higher mass ratio compared to the full distribution: in Horizon-AGN the mean mass ratio is 0.05 for the general population, and 0.17 for sel, for NewHorizon they are respectively 0.004 and 0.1 so that even minor mergers, i.e., with mass ratio $<1:4$ contribute to sourcing MBH mergers, and in fact they represent more than 50% of the sample.

In Horizon-AGN, the masses of galaxies sourcing MBH mergers are somewhat larger than those of the full merging population. The reason is that the occupation fraction of MBHs is larger in more massive galaxies (for Horizon-AGN see Figures 9 and 10 in Volonteri et al., 2016, for NewHorizon the trend is similar at $M_{\rm gal}>10^{8.5}\,{\rm M_{\odot}}$ and follow the same slope at lower galaxy mass, reaching $\sim 0.1-0.3$ at $M_{\rm gal}=10^{6}\,{\rm M_{\odot}}$. Note that both in Horizon-AGN and in NewHorizon some galaxies host multiple MBHs). In the case of NewHorizon, the increase in galaxy mass seems to disappear, but this is because most mergers involving galaxies with mass $>10^{10}\,{\rm M_{\odot}}$ are very minor mergers. If we restrict the mass distribution to mergers with mass ratio $>0.1$ the masses of primary galaxies in mergers leading to MBH mergers are 0.5 dex larger than the mean galaxy masses of the global merging population (the same is true if we perform the same mass ratio cut in Horizon-AGN).

When we require that the dynamical friction timescales added in post-processing bind the MBHs within the Hubble time at $z=0$, the mean mass ratios increase to 0.22 and 0.14 for Horizon-AGN and NewHorizon, while the galaxy masses decrease slightly. Dynamical friction timescales are shorter the larger the mass of the infalling MBH, and therefore the more massive its host galaxy, $M_{\rm gal,2}$, in general, while they are shorter the smaller $M_{\rm gal,1}$ is, via $\sigma$ in Eq. 1. Furthermore, as mergers involving massive galaxies happen at later times, there is less time to complete the dynamical friction phase by $z=0$.

In Fig. 7 we show the time delay between when galaxies merge and when their MBHs are numerically merged, $\Delta t=t_{\mathrm{nmerg,BH}}-t_{\mathrm{merg,gal}}$, for the sel selection. We recall that the time of the merger is defined as when two galaxies have the same main descendant, i.e., only one galaxy is identified by the halo finder. These delays do not include any of the timescales added in post-processing, therefore they are lower limits to the time when MBHs actually merge. Overall, there is a wide range of delays at any given galaxy mass or mass ratio and there is no simple fit/trend to describe the delays.

Typically when mergers occur, the signatures of these events in the remnant morphologies, structures, and kinematics will last for around 1 - 1.5 Gyr, at most (Conselice, 2006). This is however for the central parts of the remnant galaxy which is most readily visible, and often the only parts visible. Outer tidal features may persist for some time, and shells for far longer (Pop et al., 2017) but both of these are very difficult to observe as they only exist in low surface brightness light, which even for nearby galaxies is not easily seen even when using extremely deep exposures and performing careful sky subtraction. The structural peculiarities in the central parts of these remnant mergers will last for a maximum of 2 Gyr or so. If we use the merger time-scale from Snyder et al. (2017) and Duncan et al. (2019) of 2.4 Gyr $(1+z)^{-2}$ we obtain the result shown in the bottom panel of Fig. 7 demonstrating that this merger time has elapsed before most galaxies will merge their black holes. This 2.4 Gyr time, scaled by redshift, is the pair timescale to go from 30 kpc to effectively the start of the merger process. This is therefore the time for two galaxies to merge, but not necessarily the time it would take to coalesce into a single system. After this 2.4 Gyr time, scaled by redshift, there will be some time afterwards when the galaxy is still morphologically and structurally distorted, before it dynamically relaxes. However, this timescale would be about $1.5(1+z)^{-2}\,{\rm Gyr}$ (e.g., Snyder et al., 2017). This time-scale is shorter at higher redshifts, meaning that the merger signatures will disappear even faster at earlier times.

Overall, taking $t_{\rm obs}=2.4(1+z)^{-2}\,{\rm Gyr}$, the values shown in Fig. 7 are upper limits to the time we could identify these systems as mergers through pair selection or a morphological approach such as the CAS (concentration C, asymmetry A, and clumpiness S) method for ongoing mergers, or the use of visual morphologies or machine learning techniques to identify peculiar structures (e.g., Conselice, 2003). For the simulation, we consider the time elapsed from the galaxy merger and the MBH numerical merger, $\Delta t$, which sets a lower limit on MBH merger timescales as we do not include any of the delays discussed in Section 3.2. The choices have been very conservative: for MBHs we do not include all the time delays after the numerical merger, and therefore underestimate the time between galaxy and MBH merger. For galaxies, the $2.4(1+z)^{-2}$ Gyr timescale is the time to start galaxy coalescence from when they are separated by 30 kpc. The $1.5(1+z)^{-2}$ Gyr timescale is the post-merger time over which discernible features are observable. The time we consider for the MBHs is the post-merger phase, therefore it should be compared to the $1.5(1+z)^{-2}$ Gyr timescale, but we considered the $2.4(1+z)^{-2}$ Gyr timescale to be more conservative and have a more robust result.

Using this criterion, all MBHs below the yellow horizontal line will merge after the host galaxy can be observed to be in interaction. The difference between the two simulations can be ascribed to two main reasons: (i) Horizon-AGN includes mergers at $z<1$ for which $t_{\rm obs}$ is longer, and (ii) at a fixed galaxy mass MBHs are lighter in NewHorizon, therefore the orbital decay longer: the simulations include on-the-fly dynamical friction following (Ostriker, 1999), where the deceleration depends linearly on the mass of the MBH.

The result of this is that if the LISA sources are identified in a galaxy it is very likely that the system will be dynamically and morphologically relaxed given the long time-period between the galaxy merger and the merger of the central black holes, as long as no further galaxy merger intervenes during $\Delta t$. Massive black holes in massive galaxies at low redshift, i.e. those relevant for PTA, are more likely to be found in galaxies that show signs of interaction caused by the same galaxy merger that supplied the MBH for the MBH merger. MBHs hosted in relatively small galaxies at high redshift, which are most relevant for LISA, will merge in galaxies where the signs of the interaction that generated the merger have been erased. However, at $z>1$ the galaxy merger rate is very high, and a galaxy can experience one or more galaxy mergers during $\Delta t$, as can be seen in the examples in Section 4.4. It is important, however, to realise that the disturbed morphology of the galaxy is not caused by the same galaxy merger from which the MBH merger originated: without taking this effect into account risks mis-associating galaxy and MBH mergers, which would lead to incorrectly inferring short merger timescales for MBHs.

We conclude this section by comparing the galaxy⁵⁵5For a general discussion on galaxy mergers in simulations, see, e.g., Rodriguez-Gomez et al. (2015). We only note there that the galaxy merger rate of NewHorizon and Horizon-AGN is in good agreement with observations (Duncan et al., 2019). and MBH merger rate in Fig. 8. As already noted, the total number of galaxy mergers is larger than the total number of MBH mergers, because not all galaxies host MBHs, and not all galaxy mergers lead to a MBH-MBH merger. If we consider only numerical mergers, i.e., without including delays in post-processing, the change between the galaxy and MBH merger rate is mostly in the normalization, because at least some of the delays shown in Fig. 7 are short. Adding additional time to the delay between the numerical merger and the formation of the binary or the merger of the binary causes a shift in both normalization and shape, i.e., MBH mergers are shifted to later cosmic times. We show some reference cuts in mass and mass ratio to galaxy mergers inspired by the distributions shown in Fig. 6. Given that delay times are not easily connected to basic observable galaxy properties (mass, mass ratio), converting a galaxy merger rate into a MBH merger rate is non-trivial. Although the full physical picture of linking MBH and galaxy mergers is complex, it is generally possible to infer statistical MBH merger rates from galaxy merger rates, but the cuts in mass and mass ratio should be tuned to the particular MBH properties of interest.

As Fig. 9 shows, the probability for two galaxies to merge and then to contain a MBH numerical merger within 3 Gyr is highest for the higher stellar mass systems, and appears to drop in the Horizon-AGN simulation at lower masses, although this drop is less steep for the NewHorizon simulation, where low-mass galaxies are better resolved and have a higher MBH occupation fraction. This probability takes into account that not all galaxy mergers lead to a MBH merger, explaining why although most MBHs have $\Delta t=t_{\mathrm{nmerg,BH}}-t_{\mathrm{merg,gal}}<3$ Gyr in Fig. 7, the probability is generally $<1$, except at the high-mass end for mass ratios $>0.25$. What this plot shows is that we are more likely to find a MBH merger in more massive galaxies. Observationally, this is important if we want to trace the merger history of MBHs from examining the merger history of galaxies (e.g., Conselice et al. 2020 in prep). Interestingly, this plot also shows that there is a higher probability for MBH mergers to occur for major mergers where the $q_{\rm gal}=M_{gal,2}/M_{gal,1}\leqslant 1$ ratio is larger than for more minor mergers. For shorter time periods, this decreases at lower galaxy masses and lower mass ratio of galaxy mergers, and for longer time periods there is a slight increase, with differences being more pronounced for NewHorizon.

While the observability of electromagnetic counterparts, from the MBHs themselves or from the host galaxies, will be the subject of a future work, let us study some representative cases in NewHorizon, to highlight the importance of dynamical delays in determining what type of galaxies host MBH mergers. In the following figures, galaxies are represented with false-colour maps of their surface brightness through rest-frame u-g-r filters including the absorption by dust.

In the first example of a MBH merger, of MBHs ‘936’ and ‘41’ (Fig. 10), which originated from a merger between two galaxies with mass $2\times 10^{9}\,{\rm M_{\odot}}$ and $10^{10}\,{\rm M_{\odot}}$ at $z=3.18$, the MBHs are numerically merged at $z=2.6$, with an ensuing dynamical friction timescale of about 10.8 Gyr, meaning that it finishes at $z\sim 0.01$, beyond the end time of the simulation. Disturbed features in the host galaxy, with mass $6\times 10^{10}\,{\rm M_{\odot}}$, at $z=1.19$ are not related to the memory of the initial merger at $z=3.18$ but to a merger with mass ratio 0.73 at $z=1.22$.

The second example is another merger involving MBH ‘936’, this time with MBH ‘96’ (Fig. 11). It started with the same galaxy merger as ‘936’ and ‘41’, with MBHs ‘41’ and ‘96’ in the same galaxy at that time, but the MBH numerical merger happens much later, at $z=1.38$, when $2.67$ Gyr have elapsed from the initial galaxy merger. The galaxy, however, has experienced another 0.2 mass ratio merger at $z=1.39$. Dynamical friction ends at $z=0.6$, when the galaxy is almost $7\times 10^{10}\,{\rm M_{\odot}}$.

In the third MBH merger, of MBH ‘796’ and MBH ‘9’ (Fig. 12), the initial galaxy merger took place at $z=3.76$, with a mass ratio of 0.77, between two galaxies with mass $\sim 10^{9}\,{\rm M_{\odot}}$. The two MBHs are numerically merged at $z=1.58$, $\sim 2.53$ Gyr later, when the galaxy is experiencing a new merger, with a galaxy 5 times lighter. The ensuing dynamical friction and binary evolution timescales end at $z=1.19$, corresponding to a cosmic time of $\sim 5.29$ Gyr. During the latter phase of dynamical evolution, the galaxy has grown to $4\times 10^{10}\,{\rm M_{\odot}}$ and experienced four additional mergers with mass ratio $>0.1$.

The fourth case, MBHs ‘166’ and ‘656’ (Fig. 13), starts with the merger of two galaxies with mass ratio 0.32 at $z=2.24$, where the least massive galaxy has a mass of $10^{10}\,{\rm M_{\odot}}$. The numerical merger happens 0.5 Gyr later, when the galaxy has not relaxed yet from the merger, and our estimates for subsquent delays are very short, less than 0.1 Myr each, since the MBHs numerically merge very close to the centre of a galaxy with high stellar and gas density. This is a genuine case where the host galaxy is disturbed because of the galaxy merger from which the MBH merger originated.

The fifth example, of MBH ‘455’ merging with MBH ‘156’ (Fig. 14), originated in a very minor galaxy merger at $z=7.4$ with the most massive galaxy weighing $2.5\times 10^{7}\,{\rm M_{\odot}}$ and the least massive $10^{6}\,{\rm M_{\odot}}$. The two MBHs are numerically merged at $z=6.4$, $\sim 0.18$ Gyr later, and despite the short intervening time the host galaxy has experienced a further merger with a mass ratio of 0.18. The dynamical friction timescale is $\sim 1.87$ Gyr, ending at $z=2.43$. When the MBHs form a binary, the host galaxy has increased its mass to $7\times 10^{9}\,{\rm M_{\odot}}$. By this time the host galaxy has experienced two additional mergers with mass ratio $>0.1$ and is on its way to another major galaxy merger, but none of these are the origin of the MBH merger in question.

We end therefore with a word of caution on looking for the host galaxies of LISA MBH mergers among galaxies with signs of interactions: a galaxy could be involved in a merger at the time of a MBH binary coalescence, but it would generally be a random coincidence, due to the elevated merger rate of galaxies at high redshift.