Massive Black Hole Mergers with Orbital Information: Predictions from the ASTRID Simulation

The Astrid simulation is a large-scale cosmological hydrodynamic simulation in a $250\,{\rm Mpc}/h$ box with $2\times 5500^{3}$ particles. Astrid contains a statistical sample of halos which can be compared to future survey data from JWST, while resolving galactic halos down to $10^{9}\,M_{\odot}$ (corresponding to 200 dark matter particles). Astrid has been run from $z=99$ to $z=3$. It contains models for inhomogeneous hydrogen and helium reionization, baryon relative velocities and massive neutrinos, as well as ’full-physics’ star formation model, BH accretion and associated supernova and AGN feedback respectively. The BH model includes mergers driven by dynamic friction rather than repositioning. Our treatment of MBHs largely follows the BlueTides  simulation in terms of the BH accretion and feedback, which is based on the earlier work by Di Matteo et al. (2005); Springel et al. (2005). Compared with BlueTides , we slightly changed the seeding scheme of MBHs by drawing the seed mass from a power-law distribution instead of using a universal seed mass. Furthermore, we use a dynamical friction model (tested and validated in Chen et al., 2021) to evolve the binary black holes and include the sinking and merger of MBHs in the simulation in a more physical way. Here we briefly summarize the black hole seeding and dynamics treatment in Astrid, and refer to Bird et al. (2021) and Ni et al. (2021) for detailed presentations of physical models for star formation and black holes.

There are a total of 445635 BH mergers in the simulation for $z>3$. For each merger event we extract the relevant environmental variables (the density profiles of gas, dark matter and stars, and the stellar velocity dispersion) from the nearest snapshot before and after the merger. The snapshots used are separated by $\sim 20$ Myrs. In a small fraction of cases, the mergers take place within $20\,{\rm Myrs}$ after one of the MBHs are born, and so we cannot find the corresponding MBH in the previous snapshot. We remove these mergers from the catalog, after which 440999 mergers remain.

From the snapshots immediately before and after the merger, we identify the host halos and subhalos containing the binaries using FOF and SubFind, respectively. Out of the mergers that remain in the catalog, we further remove those not associated with any halo/subhalo, and those whose host galaxy has less than $200$ star particles. The hosts for these binaries are not well resolved in our simulation, so we cannot reliably compute the binary hardening time in post processing. This leaves us with a final catalog of 430938 black hole merger candidates. For each host halo identified, we define the halo center as the position of the particle with the minimum potential, and calculate galaxy properties such as the density profiles and half-mass radius with respect to this point.

In Figure 1, we show the last few orbits of a few selected BH pairs in our merger catalog plotted on their host galaxies’ stellar distributions. The distance from left to right of each image is $8\,{\rm ckpc}/h$. The brightness corresponds to the stellar density, and the colours show the stellar age with older stars being redder and younger stars being whiter. The red curves are the BH pairs’ positions relative to their center of mass.

Given the images in Figure 1, a natural way of measuring the orbital eccentricity is to use the shape of the orbits just before the numerical merger, and this is the first approach we take.

On $\sim{\rm kpc}$ scales, since most orbits are not Keplerian except those of the most massive BHs and the orbits are constantly shrinking, the BH orbits do not fit an ellipse. Instead, they exhibit the feature of a Rosetta orbit (the feature is most prominent in e.g. second row, second column of Figure 1, although standard Rosetta orbits do not shrink over time). For orbits resulted from the spherically symmetric potential, we can characterize the eccentricity by the size of the inner radius and the size of the outer radius. More specifically, for each orbit, we define $\Delta r_{2}$ to be position of the secondary BH with respect to the center of mass, and we take the local minimum of $\Delta r_{2}$ as the (generalized) periapsis of the orbit, and the local maximum of $\Delta r_{2}$ as the apoapsis. Then, we represent the orbital eccentricity of the binary by the generalized eccentricity, defined for a spherical potential as: (Binney & Tremaine, 2008):

where $r_{\rm apo}$ and $r_{\rm per}$ are the apoapsis and the periapsis of the orbit, respectively. To distinguish between the measurement of the two methods, we will use the subscript ”sh” to refer to the measurements from this shape-based method. We average the eccentricity measurements over the last three orbits. We note, however, that the distribution in eccentricity does not change significantly when we take the average of the last one, two or three orbits.

In addition to the shape-based measurement in §3.1, we also calculate the generalized orbital eccentricity by simply solving the orbital equation. Using these two independent methods we will then be able to compare the robustness of the BH orbit eccentricity distribution measurement from the simulations.

When the BH merger occurs in the simulation, the separation between the black hole pair is $\sim 3$ ckpc/$h$. At this distance, the gravitational potential is dominated by the surrounding stars and dark matter particles instead of the BHs themselves. Under such circumstances, the orbit of the satellite BH is non-Keplerian, as we have shown in Section 3.1. In the case of a spherical potential, the (generalized) apoapsis and periapsis can be obtained by solving the generalized orbital equation (Binney & Tremaine, 2008):

where $\phi(r)$ is the gravitational potential computed from the density profile of surrounding particles, $E$ is the total energy per unit mass and $J$ is the angular momentum per unit mass of the secondary black hole with respect to the host galaxy center. The larger root of the equation corresponds to $r_{\rm apo}$ and the smaller root is $r_{\rm peri}$.

When solving Equation 5, we take $E$ and $J$ to be the average energy and angular momentum over the last half-orbit (i.e. from the last local maximum to the last local minimum of the BH separation) of the BH. We did not take the average over a more extended period of time because the black hole pair is constantly losing energy. After getting the two apses, we again use Equation 4 to calculate the orbital eccentricity. We refer to this method as the energy method, and use subscript ”en” when showing results.

Figure 2 shows a comparison between the (generalized) eccentricity measurements from the shape method and the energy method. The left panel shows the distribution of eccentricities for all the mergers in the simulation. The measurements from both methods show that the BH binary population dominated by highly eccentric orbits, with a peak at $\epsilon\sim 0.85$ for the shape-based method and $\sim 0.75$ for the energy-based method. Comparing the two distributions, we see that the shape measurement generally produces a distribution with higher eccentricities than the energy method. In the middle panel we show a scatter plot of the eccentricities from the two measurements. There is a positive correlation between the two eccentricities, with the majority of the measurements close to the diagonal line. This means that the two measurements are not only close in distribution, but also yield correlated results for each individual orbit. Similar to what shown by the 1D distributions, the shape method predicts higher values of eccentricity for most pairs than the energy method (typically $\sim$ 10% lower).

In addition to the eccentricity, in the right panel we further compare the apoapses and periapases from the two measurements. Overall, we can see that the apoapsis peaks around $1\sim 3$ kpc, while the periapsis peaks around $0.1\sim 0.7$ kpc. Again there is a good alignment between the two measurements, with the peaks distributed close to the diagonal line. In the majority of cases the shape measurement gives an larger apoapsis value.

To estimate the binary hardening time, we will use the measured binary eccentricities as an input to the model. By doing so, we do not consider any time-evolution of the binary eccentricity due to dynamical friction beyond the point of numerical merger (Colpi et al., 1999; Hashimoto et al., 2003). In particular, we will only show results using the values from the energy-based method ($\epsilon_{\rm en}$) in later sections, and we have tested that the effect of using the shape-based values is minor compared with the uncertainties from other sources (e.g. density in the central region of the galaxy).

In Astrid, we have already accounted for the dynamical friction timescale above the resolution limit, leading to significant delays of in-simulation mergers compared to the traditional MBH repositioning methods. However, dynamical friction will continue to dominate over other delay processes on scales of $10\sim 100$ pc (e.g. Kelley et al., 2017a), which is beyond our current resolution. In this section, we will compute the unresolved DF timescales for the MBH mergers, and compare with the in-simulation DF timescale.

For the in-simulation DF time, we measure it in the following way: for each black hole pair that merge within the simulation at $z_{\rm merge}$, we track their trajectories before the merger event, and find the redshift $z_{\rm encounter}$ at which they are first within $2\epsilon_{g}$ of each other. $z_{\rm encounter}$ is the approximate time at which the BHs would merge if we did not account for the dynamical friction time at all (note that under the reposition model, BHs usually merge even before $z_{\rm encounter}$). We consider the time difference between $z_{\rm encounter}$ and $z_{\rm merge}$ as the in-simulation DF time, $T_{\rm DF,sim}$. Among all BH mergers in the simulation, 5713 mergers ($\sim 1.4\%$ of the whole merger population) happen at the first encounter.

For the post-processed DF time $T_{\rm DF,post}$, we adopt the treatment from Merritt (2013) and Dosopoulou & Antonini (2017), who modifies the Chandrasekhar formalism (e.g. Chandrasekhar, 1943; Binney & Tremaine, 2008) to include the effect of the secondary BH embedded in a tight core of stars brought in from the secondary galaxy. The increased dynamical friction allows the secondary to sink faster towards the primary galaxy’s center, and thus the resulting dynamical friction time is less than the prediction from the canonical Binney & Tremaine (2008) treatment assuming a bare BH. In Dosopoulou & Antonini (2017) the assumption is that the mass of stars bound to the secondary BH is 1000 times the mass of the BH itself, and the resulting dynamical friction timescale is:

where $\text{log}(\Lambda)$ is the Coulomb logarithm, $M_{2}$ is the mass of the secondary black hole. For the initial separation $r$, we use the radius of the circular orbit that has the same energy as the last orbit of the binary before the (numerical) merger. Note that the model in Equation 6 does not account for the effect of non-circular orbits on the DF time. Taking the eccentricity into consideration can further reduce the estimated DF time (e.g. Taffoni et al., 2003). Following the method in Chen et al. (2021),we compute the Coulomb logarithm by:

where $M_{\rm 2}$ is the mass of the secondary black hole and $v_{\rm BH}$ is the velocity of the secondary black hole with respect to the host galaxy center.

Figure 3 shows the comparison between the in-simulation dynamical friction time and the post-processed dynamical friction time from above. The top panel shows the overall distributions of the DF times. The two distributions are on the same order of magnitude at around $10^{2}$ Myrs, with a range from 10 Myrs to 1 Gyrs. For most BH pairs, $T_{\rm DF,sim}$ is longer than $T_{\rm DF,post}$. This means that accounting for dynamical friction to the simulation, we have already included about half of the total dynamical friction delay effects. Note that both DF timescales are shorter than 1 Gyr. In the case of the resolved DF time, this is mainly due to the fact that most of the black holes have not existed for more than 1 Gyr at $z=3$.

In the bottom panel of Figure 3, we show the relation between the two DF times and the binary mass ratio $q$. Here we only show the times involving at least one non-seed MBH, defined as mergers with $M_{1}>2M_{1,\rm seed}$. This is because our merger population is dominated by seed MBHs which have not grown out of their dynamical mass and thus the in-simulation DF time estimation is not exact. From Equation 6, we can see that the DF time is correlated with the mass of the primary galaxy (and MBH) through $\sigma$, and that it is inversely proportional to the secondary black hole mass. Hence, we expect that minor mergers will have longer decay timescales, and in the plot we do see a negative correlation between $T_{\rm DF,post}$ and $q$. For the in-simulation DF, although this relation is not imposed explicitly, we still observe a negative correlation between $q$ and the DF time. This indicates that the negative correlation is still captured by the in-simulation dynamical friction modeling.

Once the two MBHs become gravitationally bound, the dynamical friction formalism is no longer a valid approximation, and individual interactions between singular stars and the binary must be considered. These interactions extract angular momentum from the binary, driving them closer to each other (e.g. Merritt, 2013; Vasiliev & Merritt, 2013). This regime is the loss-cone scattering (LC) regime, which refers to the specific cone in parameter space where stars have to exist in order to extract angular momentum from the binary (e.g. Frank & Rees, 1976; Lightman & Shapiro, 1977). On even smaller scales, the binary will enter the gravitational wave regime where it will evolve until coalescence. Once the binary enters the gravitational wave regime, its dynamics follow the formalism of Peters (1964) at small separations of $10-1000$ Schwarzschild radii.

For the loss-cone scattering and gravitational-wave hardening phase, we adopt the analytical prescription in Vasiliev et al. (2015) (V15 hereafter). However, the time estimation in Equation (25) of V15 assumes a single family of $M_{\rm tot}-\sigma_{\rm inf}-r_{\rm inf}$ relation, and thus it may over-simplify the properties of the galaxies hosting the merger events. Hence we will adopt the V15 formalism but with some slight changes, so that we can use the host galaxy properties measured from the simulation. In this Section, we first explain how we measure the relevant galaxy properties, then we give our binary hardening time estimation by combining analytic modeling with the measured properties.

We calculate the rate by integrating the number of mergers in the simulation over redshifts, incorporating the cosmic volume at the given redshift:

where ${\rm d}V_{c}$ is the comoving volume element of the universe at a given redshift and $n(z)$ is the number of mergers at that redshift. The $1/(1+z)$ term redshifts the infinitesimal time element in ${\rm d}z/{\rm d}t$ to the observer frame time interval.

To calculate this rate from our simulation, we take the finite-interval approximation:

where $\Delta z$ is the width of the redshift bin, $N(z)$ is the total number of mergers within that redshift bin, and $V_{\rm sim}=(250\,{\rm Mpc}/h)^{3}$ is the volume of our simulation in comoving units.

To clearly see the effect of each stage of the delay, we calculate three different rates. We first compute the ”Sim” rate which uses the numerical merging time as the redshift of the merger (also see DeGraf et al., prep). Then we add the post-processed DF time to the numerical merger time to compute the ”DF-only” rate. Finally, we further account for the binary hardening timescales and calculate the ”DF+Hard” rate.

In the left panel of Figure 8, we plot the merger rates with different levels of post-processed delays, for the whole merger population in Astrid. First, we notice that the number of mergers keeps increasing with decreasing redshift for all three models. This is because we keep seeding BHs as structures form and grow, and that we have not reached the peak in seeding rate at $z=3$. Without considering any post-processed delays (”Sim”), we expect a total of $\sim 1.8$ mergers per year of observation down to $z=3$. The post-processed DF time does not significantly impact the total observed merger rate (”DF-only”), with a $\sim 50\%$ decrease at the highest redshift ($z\sim 8$). The binary hardening time has the most significant impact on the merger rate at all redshifts (”DF+Hard”). We see that the merger rate is reduced by a factor of $3\sim 7$ after adding the delay from binary hardening. The resulting merger rate is $0.3\sim 0.7$ at $z>3$. Here the upper limit is given by assuming $r_{\rm ext}=2\epsilon_{g}$ and the lower limit is given by $r_{\rm ext}=1.5\epsilon_{g}$. On the bottom panel, we show the ratio between the delayed merger rate and the simulation merger rate as a function of redshift. For both DF-only and DF+Hard delays, the fractional rates get higher at lower redshifts. This is a result of the high-redshift mergers being pushed down to low redshifts.

On the right panel, we show the mass distribution of the two MBHs involved in each merger. The dashed lines corresponds to the simulation merger without any delays, and the solid lines shows the distribution of the merger population after the DF+hardening delays. First, we can see that both before and after the delay, the merger population is dominated by seed-mass mergers (the ones enclosed by the vertical dashed lines), with $M_{1}$ evenly distributed across the seed masses and $M_{2}$ concentrated on the lower-mass end of the seeds. It is also this seed-mass merger population that gets suppressed the most by the delay. From the ratio between the mass functions shown in the bottom panel, we see that for the seed-mass mergers, only $\sim 15\%$ still merge at $z>3$ after the delays, whereas at the high-mass end this fraction increases to $50\%$.

In order to disentangle different merger populations, in Figure 9 we further split the rate by how many seed MBHs are involved in the merger. The left panel shows the merger rates for the seed-mass population, where the masses of both MBHs are below two times their seed masses. This population makes up $\sim 60\%$ of the mergers. At $z>5$, the seed-seed mergers are strongly suppressed by the binary hardening delays because the stellar density is relatively low. The middle panel shows the mergers with the only more massive MBH grown beyond two times its seed mass. At $z=3$, the rate from this group is comparable to the rate from the seed-seed mergers. However, the number decreases more steeply as we go to higher redshifts. Compared to the seed-seed mergers, this group has a higher mass ratio and thus a longer DF time. The effect of the binary hardening delay, however, is smaller because of the higher-density in the remnant galaxy. Finally, on the right panel, we show the more massive and major mergers. Compared to the previous two groups where at least one seed-mass MBH is involved in the merger, the mergers from this group is $\sim 6$ times lower. The effect of delay is also the smallest. In particular, we noticed that the DF-only rate is very similar to the simulation rate. Even for this group where the effect of delays is the smallest, the merger rate is still suppressed by $>50\%$ at each redshift compared to the simulation merger rate.

Figure 10 shows the distribution of MBH mergers on the $M_{\rm tot}-z_{\rm merge}$ plane for both the simulation and delayed mergers, color coded by the number of mergers per Myr. Without any delay, the majority of the merger events are seed-seed mergers around $z=3-4$. As we would expect from the black hole mass growth over time, we see more massive mergers at lower redshifts. The middle panel shows the same merger population with the post-merger DF time added. As was discussed in the previous paragraph and in Section 4.1, the post-processed DF peaks around $200$ Myrs and does not significantly delay the mergers. Here, we see a slight shift of the merger population towards lower redshift.

On the right panel, we show the distribution after considering the DF delay and hardening phase. Note that since the final simulation output is at $z=3$, all the data points at $z<3$ are the results of delayed $z>3$ numerical mergers, and are not representative of all merger events at $z<3$. Compared with the other two panels, we see a significant shift of the mergers towards lower redshifts. The population that is most significantly shifted are the low-mass mergers with $M_{\rm tot}<10^{6.5}M_{\odot}$, while the most massive binaries are still able to merge at relatively high redshifts. This is a consequence of the large hardening time scale of smaller BHs associated with lower $\rho_{\rm inf}$.

From the previous section, we have seen that while some low-mass mergers are significantly delayed and do not merge at $z>3$, $\sim 15\%$ of them still do. For the non-seed mergers, although the delay is generally less significant, we still see a $50\%$ decrease in merger rate when accounting for the delays. Now we will investigate which part of the merger population gets significantly delayed, and which still manages to merge at high redshifts.

In Figure 11 we show the properties of MBHs involved in both the simulation mergers and the delayed mergers. The top row shows the properties of the non-seed mergers, and the bottom row shows the properties of the seed-seed mergers. We start by looking at the mass distribution of galaxies hosting the mergers (shown in the first column). For the simulation mergers consisting of two non-seed MBHs, the masses of the host galaxies peak at $4\times 10^{9}M_{\odot}$. For systems that still merge after the delays, we see a clear shift towards the higher-end in stellar masses with a peak at $\sim 10^{10}M_{\odot}$. This is because for more massive galaxies, the high stellar density enables more efficient hardening through loss-cone scattering, and thus the delay time is shorter (also see Figure 7). For mergers involving two MBH seeds shown on the bottom, we observe a similar trend. Overall, seed mergers reside in less massive galaxies with stellar masses below $4\times 10^{8}\,M_{\odot}$. The delayed merger events also pick up the more massive galaxy population out of the simulation mergers with galaxy masses distributed around $10^{9}\,M_{\odot}$.

In addition to the stellar environment which plays an important role in the delay time estimation, the seeding redshift of the MBHs can also affect whether the two MBHs still merge at a high redshift after the delay. This is shown in the second column of Figure 11. While the seeding redshift of the simulation merger MBHs is $z\sim 7$, the MBHs involved in delayed mergers are seeded as early as $z=10$. For the seed-seed mergers shown on the bottom, the overall seeding redshift is lower, but we also see a shift towards higher redshift when comparing the delayed mergers to the simulation mergers. The bias towards early MBH seeding for delayed mergers is also correlated with the higher host galaxy mass we have seen earlier: because the delayed mergers favor earlier seeds, they also tend to reside in galaxies that are massive enough at high redshifts to host an MBH seed.

On the right two columns, we examine the properties of other MBHs embedded in the host galaxy of the mergers. The third column shows the total number of MBHs embedded in the host galaxy of the merger, in the snapshot immediately following the numerical merger (so the merging MBHs will be counted as one object). The fourth column is the mass ratio between all MBHs in the host galaxy and the merging system. For both the seed and non-seed merger populations, the merging system is the sole MBH in the host galaxy in a majority of mergers. For the non-seed population, there is still a $>50\%$ fraction of mergers happening next to a third MBH (or more). Interestingly, the delayed merger systems favor galaxies with more MBHs besides the merging ones (also correlated with larger galaxy masses). Nonetheless, the merging system is still the most massive MBH in its host galaxy in the most cases when we look at the $M_{\rm BH,gal}/M_{\rm tot}$ ratio.

When constrained to seed-seed mergers, we see that $\sim 70\%$ of the mergers are the single MBH in the host galaxy. The delayed mergers also tend to pick out the galaxies with more MBHs compared to the simulation mergers. However, contrary to the non-seed case where the merging MBH is more massive than the other MBHs in the same galaxy, for seed-seed mergers that do occur near a third MBH, the mass of the third MBH is more likely to be larger. This can be seen from the fact that the $N_{\rm BH,gal}$ distribution is more peaked at $N_{\rm BH,gal}>1$ compared to the $M_{\rm BH,gal}/M_{{\rm tot}}>1$ distribution (it means that if there is a third MBH, its mass can be larger than $M_{\rm tot}$ in some cases, resulting in the longer tail of $M_{\rm BH,gal}/M_{{\rm tot}}$).

From the investigations above, we conclude that the $z>3$ mergers after the DF and hardening delay make up a small and biased sample of the simulation mergers. In particular, they are systems with MBHs seeded earlier and embedded in more massive galaxies compared to the overall simulation merger population. Moreover, the majority of the merger remnant is the only MBH in its host galaxy, especially for the seed-mass mergers. However, the delayed mergers tend to pick out more systems that have other nearby MBHs in the remnant galaxy compared to the overall simulation merger population.

MBH binaries provide a variety of signals measurable by LISA since their chirp evolution in the frequency domain occurs near the low-frequency band edge of the LISA sensitivity curve. Binaries with $10^{5}-10^{7}M_{\odot}$ total mass will provide a measurable inspiral, merger, and ringdown leading to signals out to the cosmic horizon (Amaro-Seoane et al., 2017). The binary inspiral is the initial stage of binary black hole coalescence when the two MBHs orbit one another at separations greater than the innermost stable circular orbit (ISCO; $R=6GM_{\rm BH}/c^{2}$). At these separations, the orbit is usually treated with a post-Newtonian formalism. The merger stage follows the binary inspiral with a highly non-linear relativistic process. This process continues until the binary components form a single event horizon, leading to ringdown.

We use the characteristic strain, $h_{s}$, to model the binary signal which accounts for the time the binary spends in each frequency bin (Finn & Thorne, 2000). The characteristic strain is given by (e.g. Moore et al., 2015):

where $\tilde{h}(f)$ represents the Fourier transform of a time domain signal.

To generate the waveforms, we use the phenomenological waveform PhenomD (Husa et al., 2016; Khan et al., 2016) implemented within the gwsnrcalc Python package (Katz & Larson, 2019). The input parameters are the binary masses, merging redshift, and the dimensionless spins of the binary. For the MBH masses, we do not account for mass growth after the numerical merger. However, we note that the MBH can potentially gain a significant fraction of its mass during the $>1\,{\rm Gyr}$ of time in the dynamical friction (e.g. Banks et al., 2021) or loss-cone scattering phase. The dimensionless spin $a$ characterizes the alignment of the spin angular momentum with the orbital angular momentum, and the value of $a$ ranges from $-1$ to $1$. However, we do not have any information on the spin of the SMBHs in our simulation. Therefore, following the argument in Katz et al. (2020), we assume a constant dimensionless spin of $a_{1}=a_{2}=0.8$ for all binaries (e.g. Miller, 2007; Reynolds, 2013).

In Figure 12, we show the distribution of the merging frequency $f_{\rm merge}$ and the strain at this frequency for all binaries in the simulation, before and after applying the DF+Hardening delay models. To evaluate the detectability of the population with LISA, we also plot the proposed LISA sensitivity curves. We use the LISA sensitivity configuration from the LISA Mission Proposal (Amaro-Seoane et al., 2017), and we use $h_{N}=\sqrt{S_{N}}$ (Moore et al., 2015) to convert from the proposed power spectral density $S_{N}$ to strain $h_{N}$.

In the left panel, we show example waveforms for binaries of different masses but similar numerical merging time. The thick curve shows the waveform assuming $\epsilon=0$, with the dot representing the merging frequency $f_{\rm merge}$. We will discuss the thin lines with non-zero eccentricities in later sections. From the example waveforms, we see that at a fixed source redshift, the more massive binary has a higher strain amplitude. However, this does not necessarily lead to a more significant detection, because the lower frequency at which the wave is emitted falls into the region where the LISA sensitivity is worse. Out of these three binaries, the two least massive binaries are detectable by LISA while the most massive one is not. After the DF and hardening delays, all curves have higher strain amplitudes, as the strength of the signal is negatively correlated with redshift.

After looking at individual cases, we turn to the whole binary population. On the right panel, we show the distribution of $f_{\rm merge}$ and $h_{s}(f_{\rm merge})$ for Astrid mergers, after the post-processed delays. We have masked the signals from $z<3$ mergers in light grey, as they are purely due to the post-processed delays, and are not part of our simulation predictions. The majority of merger events within the simulation lie above the LISA sensitivity curve. From example waveforms, we see that once any given GW signal crosses the detector sensitivity curve, the ratio of the signal to the sensitivity curve rapidly increases by a few orders of magnitude. Since the merger population is dominated by seed-seed mergers, we see a peak around $f_{\rm merge}\sim 10^{-2}$Hz, corresponding to the example green curve. Finally, we demonstrate the shift of the signal due to the delay model by the colored arrows. The tail of the arrows indicates the location of the frequency/strain before the delays. The head of the arrows are the signals after the delays. We see that in the example cases, the signal shifts to the high-strain, high-frequency region of the plane. This is mainly because of the delay of the mergers from $z>3$ to $z<3$.

In the previous section we have shown a single $h_{s}-f$ relation by assuming circular orbits for the binaries. In this section, we will utilize the eccentricity measured from the simulation when calculating the strain and signal-to-noise ratio (SNR) for each binary.

The GW strain from an individual, eccentric source can be related to that of a circular source as (e.g. Amaro-Seoane et al., 2010; Kelley et al., 2017b):

where $h_{\rm r,circ}$ is the characteristic strain of a circular source given by Equation 18, $g(n,\epsilon)$ is the GW frequency distribution function given by Equation 20 in Peters & Mathews (1963) with $\sum_{n=1}^{\infty}g(n,\epsilon)=F(\epsilon)$, where $F(\epsilon)$ is defined by Equation 12.

During the GW-driven inspiral, the orbital eccentricity also evolves according to Peters (1964) Equation (5.7), such that it decays towards zero as the binary inspirals toward merger. This will affect the orbital frequency by:

where $\epsilon_{0}$ is the initial eccentricity at the reference frequency $f_{0}$.

In Figure 12, the multiple thin lines are the waveforms from higher-order harmonics assuming eccentric orbits. For circular orbits, the GW is emitted at a single frequency at a fixed separation, while the eccentric binaries emit GW at higher-order harmonics at a given time. One consequence of this is that the energy dissipated in higher-order harmonics is below the detection sensitivity, and thus the signal will be smaller compared with the circular orbits.

We note that when using the simulation measurement of the orbital eccentricity as the initial eccentricity in the inspiral phase, we did not account for the possible increase in $\epsilon$ during the loss-cone scattering phase (see, e.g. Sesana, 2010; Kelley et al., 2017b).

Although the strain in Figure 12 is a good estimation of the detectibility of a circular binary, for the eccentric case a more careful prediction comes from the signal-to-noise ratio (SNR). The SNR is estimated by integrating the ratio of the signal to the noise in the frequency domain. The sky, orientation, and polarization averaged SNR is given by :

where $f_{\rm start}=f(t_{\rm start})$ and $f_{\rm end}=f(t_{\rm end})$, with $t_{\rm start}$ and $t_{\rm end}$ representing the starting and ending time of when the signal is observed. Note that here we are assuming eccentric waveforms for the binaries, and thus $h_{s}$ is given by the sum over different modes following Equation 19. As it is not computationally feasible to sum an infinite number of modes, we truncate the sum in Equation 19 at $n=50$ and we have checked that the difference between the first $50$ and the first $100$ modes is less than $5\%$.

For the current configuration, we assume that the LISA observation lasts for 4 years. We further assume a most optimistic SNR for all mergers by taking $t_{\rm end}=t_{\rm peak}$ and $t_{\rm start}=t_{\rm peak}-4{\rm yrs}$. Under this assumption, we are always integrating the part of the waveform where the strain is maximized. However, as was discussed in Salcido et al. (2016) and Katz et al. (2020), the actual SNR may be smaller if there is an offset between the LISA observation window and the merger time of the binary.

Figure 13 shows the distribution of the SNR computed for all mergers in the simulation. The left column shows the joint distribution of SNR and the merging redshift. The top row is the SNR computed with merging redshifts before the DF and hardening delays, and the bottom row is the SNR after the delay time is applied. As was expected from the simpler calculation shown in Figure 12, the majority of the binary population in the simulation has a SNR larger than the LISA detection threshold of 8 (plotted as dashed gray lines). The ones that falls below the SNR cut are mainly massive mergers with $M_{\rm tot}>10^{7}\,M_{\odot}$. When we account for the delays, the mergers are pushed towards lower redshifts, and the resulting SNR is higher for each event. However, when restricted to $z>3$, the

The middle panel shows the effect of delays and the SNR cut on mergers with different masses. The SNR cut removes all mergers with $M_{\rm tot}>10^{8}\,M_{\odot}$ from the LISA-detectable population. On the low-mass end of MBH mergers, the reduction results from the DF and binary hardening delays. Combining both the delays and SNR cut, we see that the overall detectable mergers at $z>3$ are $\sim 15\%$ of the original Astrid merger population across all masses. The seed-mass mergers still dominate over other events even though they are most strongly suppressed by the delays. Finally, on the right panel we show the mass distribution of the two MBHs involved in each detectable events. The majority of these events are expected to be mergers from two seed-mass MBHs. On the high-mass end, the detectable events has a mass ratio of $q\sim 1$ (close to the diagonal line). Based on these results, the likelihood that a LISA detection comes from mergers of MBH seeds is high, but the detectable MBH seed mergers is only a small sample of the seed MBH pairs and the associated galaxy mergers.

Here, accounting for the delay time to merger affects the resulting SNR more than the eccentricity. The eccentricity itself, however, may affect the prospects for multi-messenger follow-up. For example, eccentric binaries may spend a shorter amount of time in the LISA band compared to circular binaries. Spin-orbit interactions in eccentric binaries may change the orbital inclination with respect to the line-of-sight, which may also play a role in detectability and sky localization. We will explore such effects and their implications for multi-messenger follow-up in a companion paper.