MBH binary intruders: triple systems from cosmological simulations

ABSTRACT

Massive black hole (MBH) binaries can form following a galaxy merger, but this may not always lead to a MBH binary merger within a Hubble time. The merger time-scale depends on how efficiently the MBHs lose orbital energy to the gas and stellar background, and to gravitational waves (GWs). In systems where these mechanisms are inefficient, the binary inspiral time can be long enough for a subsequent galaxy merger to bring a third MBH into the system. In this work, we identify and characterize the population of triple MBH systems in the Illustris cosmological hydrodynamic simulation. We find a substantial occurrence rate of triple MBH systems: in our fiducial model, 22 per cent of all binary systems form triples, and |$\gt 70{{\ \rm per\ cent}}$| of these involve binaries that would not otherwise merge by z = 0. Furthermore, a significant subset of triples (6 per cent of all binaries, or more than a quarter of all triples) form a triple system at parsec scales, where the three BHs are most likely to undergo a strong three-body interaction. Crucially, we find that the rate of triple occurrence has only a weak dependence on key parameters of the binary inspiral model (binary eccentricity and stellar loss-cone refilling rate). We also do not observe strong trends in the host galaxy properties for binary versus triple MBH populations. Our results demonstrate the potential for triple systems to increase MBH merger rates, thereby enhancing the low-frequency GW signals detectable with pulsar timing arrays and with LISA.

gravitational waves, quasars: supermassive black holes

1 INTRODUCTION

Many theoretical and observational studies have shown that massive black holes (MBH) reside at the centre of most galaxies (e.g. Richstone et al. 1998). There is compelling evidence that the central MBH masses correlate with the luminosity, mass, and velocity dispersion of the galactic stellar bulge (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Ferrarese et al. 2006). This indicates a coordinated evolutionary path for MBHs and their host galaxies.

Mergers constitute an important evolutionary step for the galaxies and their MBHs, not least because they lead to the formation of MBH binaries. These binaries are the loudest GW sources in the Universe, with chirp frequencies ranging from ∼ mHz for ∼10⁶|$\rm M_{\odot }$| MBHs to ∼ nHz for ∼10⁹|$\rm M_{\odot }$| MBHs (Haehnelt 1994; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2004, 2005). Recently, pulsar timing array (PTA) experiments around the globe presented strong evidence for a nHz GW background that is consistent with a population of MBH binaries (Agazie et al. 2023; Antoniadis et al. 2023; Reardon et al. 2023; Xu et al. 2023). In the coming years, the Laser Interferometer Space Antenna (LISA) will be able to detect mHz GWs from MBH mergers in the ≲ 10⁶|$\rm M_{\odot }$| range out to z ∼ 20 (Amaro-Seoane et al. 2017). However, the detection of GWs is only possible if the binary system can reach the GW-dominated regime of inspiral (∼ mpc scales) within a Hubble time.

The time-scales for MBH binaries to inspiral and merge are highly uncertain, which is a major limitation on our ability to predict the GW and EM signatures of the MBH binary population. MBH binary inspiral, also called ‘binary hardening’, is driven by several processes at different spatial scales (Begelman, Blandford & Rees 1980). During the galaxy merger, the MBHs fall toward the galactic centre via dynamical friction (DF, Hernquist 1992; Antonini & Merritt 2012). When the mass enclosed in the binary orbit is comparable to the MBH masses, the binary forms a gravitationally bound system. This typically happens at scales ≲ tens of pc (e.g. Begelman et al. 1980; Quinlan 1996; Yu 2002).

Further hardening of the binary happens via individual stellar scattering events. The region of stellar orbital phase space that can interact with the binary is called the loss cone (LC, Begelman et al. 1980; Quinlan 1996; Quinlan & Hernquist 1997; Merritt 2013). If the LC is depleted more quickly than it can be replenished by two-body stellar relaxation, there is no guarantee that the binary can merge within a Hubble time t_H. This so-called final parsec problem (Begelman et al. 1980; Milosavljević & Merritt 2003) could be ameliorated by efficient LC refilling in galaxies with triaxial, merging, or otherwise asymmetric potentials (Yu 2002; Holley-Bockelmann & Sigurdsson 2006; Berczik et al. 2006; Holley-Bockelmann et al. 2010; Preto et al. 2011; Khan, Just & Merritt 2011; Khan et al. 2016). Binary hardening can also happen through the interactions with a circumbinary disc (CD) in gas-rich mergers (Dotti et al. 2007; Cuadra et al. 2009; Nixon et al. 2011; Goicovic et al. 2017), although it is unclear whether this mechanism can efficiently drive the binary into the gravitational wave (GW) dominated regime (Lodato et al. 2009; Muñoz, Miranda & Lai 2019; Moody, Shi & Stone 2019; Muñoz et al. 2020) and prevent stalling. None the less, many systems may still experience a significant delay between galaxy coalescence and MBH merger; for example, Kelley, Blecha & Hernquist (2017a) find hardening time-scales of many Gyr for some systems even with efficient LC refilling.

Triple MBH interactions provide a potential solution to this problem (Kozai 1962; Ryu et al. 2018; Bonetti et al. 2018b). If the binary inspiral time is longer than the time until the next galaxy merger occurs, a third MBH can enter the system. In cases where this third MBH reaches the host nucleus of the binary before the binary merges, the MBHs can undergo a three-body scattering interaction. Triple interactions can drive stalled binaries to merger on a shorter time-scale through Kozai–Lidov (K-L) oscillations in hierarchical triple systems (Kozai 1962; Lidov 1962; Naoz 2016) or through strong three-body encounters. This can potentially increase the merger rates of MBHs, and if binary stalling is common, triple interactions may be essential for driving MBHs to the GW regime where they can be observed by PTAs and LISA.

According to the K-L mechanism, if an intruder MBH forms a hierarchical triple system with the inner binary at a sufficiently inclined orbit, the system will undergo large oscillations in mutual orbital inclination and in the eccentricity of the inner binary. At very high eccentricities, GW emission can drive the MBHs to rapid merger. One caveat to this picture is that such a hierarchical triple will not evolve in isolation, but rather in the stellar and gaseous background of a galactic nucleus. Particularly in gas-rich galaxies, dynamical interactions with this background may dampen the eccentricity oscillations, thereby weakening the effects of the K-L mechanism. However, as noted above, binary hardening via gas dynamical interactions may also prevent MBH binaries from stalling in the first place. In dry systems, where the stalling scenario is more likely, the oscillations induced by K-L mechanism on the inner binary eccentricity could provide an effective means for rapidly driving the inner binary into the GW-dominated regime.

If the intruder MBH reaches the galactic nucleus at a distance comparable to the inner binary separation, a strong, chaotic, three-body interaction between the MBHs will occur. Typically, such an interaction will result in the slingshot ejection of the lightest MBH from the system, while the more massive pair forms a more tightly bound system that can merge on a much shorter time-scale (Saslaw, Valtonen & Aarseth 1974; Hills 1975; Blaes, Lee & Socrates 2002; Iwasawa, Funato & Makino 2006).

Statistics of close triple MBH encounters were studied in more detail by Hoffman & Loeb (2007). The authors found that the triple interactions increase the coalescence rate of MBHs and cause spikes of GW emission during the close encounters. A more systematic study of the statistical outcomes of triple MBH interactions was done by Bonetti et al. (2016) and Bonetti et al. (2018a), where they included all relativistic corrections up to 2.5PN order, which are very important due to the general relativistic (GR) effects inhibiting the K-L mechanism (Holman, Touma & Tremaine 1997).

Owing to the strong possibility that a triple interaction will result in a slingshot recoil kick to the lightest MBH, these systems are also relevant for their ability to produce offset, wandering MBHs, which may under some circumstances be observable as offset active galactic nuclei (AGNs, e.g. Barrows et al. 2016). A candidate slingshot recoil was recently presented by van Dokkum et al. (2023), though further observations are needed to confirm its nature. Wandering MBHs can also be produced by GW recoil kicks, which result from the asymmetric emission of GWs during a merger between MBHs with unequal masses or spins (Peres 1962; Bekenstein 1973; Loeb 2007; Volonteri & Madau 2008; Blecha & Loeb 2008; Blecha et al. 2011). Constraining the formation of triple MBH systems is therefore important for understanding the relative role of slingshot recoils versus GW recoils in producing a population of wandering BHs, and for determining the impact of such recoil events on the subsequent GW event rate. Sayeb et al. (2021, hereafter SB21) recently implemented a model for MBH binary spin evolution along with a binary inspiral model based on that of Kelley et al. (2017b). One of their key findings was that a non-negligible population of misaligned MBH binaries should exist even if gas dynamical interactions are effective at aligning MBH spins in gas-rich galaxies. Upon merger, these misaligned MBH binaries produce much larger GW recoil kicks, such that the GW recoil distribution always has a high-velocity tail (v_kick > 1000 km s⁻¹).

In order to determine the importance of triple interactions in hastening binary inspiral, producing GW sources, and ejecting MBHs from galactic nuclei, we need to know where and under what circumstances such triple systems actually form. The answer is subject to numerous large uncertainties, not least of which is the uncertainty in the underlying MBH binary inspiral time-scales. Bonetti et al. (2018b) recently combined their triple interaction modelling results with a semi-analytic model for galaxy evolution in order to predict the role of triple encounters in MBH evolution. They concluded that even if all MBH binaries stall, triple encounters should produce a stochastic GW background that is observable with PTAs. A follow-up study concluded that triple interactions could also contribute to a significant fraction of LISA events (Bonetti et al. 2019).

Models of MBH binary inspiral applied to cosmological hydrodynamics simulations also suggest an important role for triple interactions. Using the Illustris cosmological simulations, Kelley et al. (2017a) find that in their fiducial model, roughly half of all MBH binaries have not merged by z = 0, and binary inspiral time-scales of many Gyr are common. Subsequent studies using variations of this model have found similar results (Kelley et al. 2017b; Katz et al. 2020; Sayeb et al. 2021). Kelley et al. (2017a) estimate that about 30 per cent of all MBH binaries experience a subsequent merger with another MBH that overtakes the first binary it merges, and in nearly all of these the first binary is overtaken before z = 0. (Note that these post-processing models evolve the subresolution inspiral of all MBH binaries in isolation.)

The primary goal of this work is to examine in much greater detail the properties and characteristics of the triple MBHs and their host galaxies in Illustris. This constitutes the first dedicated study of triple MBH formation based on a cosmological hydrodynamics simulation. We use the same hardening prescription as in SB21 for a population of merging binaries (based on Kelley et al. 2017b) from Illustris. We identify cases where an intruder takes over a binary and forms a triple system. We look at characteristics of these systems, their environments and their host properties. We also study the dependence of these triples on variations of model parameters. Our results provide a new type of cosmological framework for understanding the role of triple interactions in MBH evolution and GW source populations, including their environments and host galaxy properties.

In Section 2, we describe the MBH population, our hardening prescription, and the method that we use for the identification of the triples. In Section 3, we describe our findings, and finally in Section 4, we discuss the caveats of our finding, and draw conclusions.

2 METHODS

Our study focuses on the population of merging MBHs from the Illustris cosmological hydrodynamic simulation suite¹ (Vogelsberger et al. 2014a; Genel et al. 2014; Nelson et al. 2015). Due to the large scale of these simulations and the computational limitations, a gravitational softening scale is defined for each particle or cell, below which the simulation is not able to resolve the physics. The MBH mergers in Illustris happen when two MBH particles are within a softening length of each other, which is typically of the order of ∼kpc. In reality, however, the MBHs are far from merger at this point. To model binary inspirals at subresolution scales, we implement a prescription that extrapolates the central density profiles of host galaxies from Illustris to calculate hardening rates due to various mechanisms: DF, stellar scattering, circumbinary gas disc-driven hardening, and GW emission (Kelley et al. 2017a, b; Sayeb et al. 2021). To identify triples, we find MBHs that experience more than one merger and track each merger in redshift and MBH separation space (z, a(z)) to find if a merger is overtaken by a subsequent merger (i.e. if the a₁(z) and a₂(z) curves of the two mergers cross; Fig. 1).

Figure 1.

A schematic of MBH merger tree (left) and separation versus time evolution (right) are shown here. The dark circles represent MBHs in the left-hand side. The numbers in each circle represent the MBHs unique ID in the Illustris simulation and the direction of time flow is from top to bottom. The highlighted areas show potential independent triples. Note that the colour of the highlighted area on the merger tree and the schematic curves plotted on the right-hand side are arbitrary; there is no one-to-one correspondence between the merger tree on the left and the schematic curve on the right. However, the solid and dashed linking lines in the diagram on the left do indicate the first and second binaries, respectively. On the right-hand side, the curves show the binary separation as a function of time. The solid curves show the first merger and the dashed lines show the subsequent or second merger. In order to identify the systems that can potentially evolve into a triple MBH, we look for repeated MBH IDs in subsequent mergers (left). If a repeated ID is identified we then look at the separation versus time curve of the current binary (i.e. first binary or inner binary) and that of the subsequent binary (i.e. second binary or outer binary). If the curves cross at some point (denoted with stars) before z = 0 then the system becomes a successful triple. Three different classes of outcomes are shown in the right figure. (a) The blue curves (leftmost solid and dashed line pair) do not intersect with each other, meaning that the first binary is being chased by the intruding MBH (outer binary) and no triple forms. (b) The red curves (center solid and dashed line pair) show a case where the first binary is overtaken by the intruder (outer binary) at small separations (<100 pc); we refer to this as an ‘ST’. (c) The green curves (rightmost solid and dashed line pair) show a binary overtaken by an intruder at large separations, indicating a ‘WT’.

The binary and triple parameter definitions are as follows: m₁ and m₂ stand for the masses of the primary MBH and the secondary MBH in the inner binary. The primary is more massive than the secondary (m₁ > m₂). The mass ratio of the inner binary is defined as q_inner = m₂/m₁, which makes it by definition smaller than unity. The combined mass of the inner binary is defined as M_bin = m₁ + m₂. As for the outer binary, the mass of the intruder is indicated by m₃. For the outer binary mass ratio q_outer, however, we want to distinguish between the ‘intruder’ MBH and the MBH that is a member of the inner binary. Thus, the outer binary mass ratio is determined by:

$$\begin{eqnarray} q_{\rm outer} = \frac{m_3}{M_{\rm bin}} . \end{eqnarray}$$

(1)

The combined mass of the binary and the intruder are defined as M_tot = m₁ + m₂ + m₃.

2.1 Illustris

The Illustris simulations are a suite of cosmological hydrodynamic simulations run using the {arepo} hydrodynamics code (Springel 2010). The {arepo} code is based on an unstructured moving mesh that combines the advantages of smooth particle hydrodynamics (SPH, e.g. Gingold & Monaghan 1977; Lucy 1977) and an Eulerian mesh-based approach (e.g. Berger & Colella 1989). The mesh is formed from Voronoi tessellations based on a set of discrete mesh-generating seeds that can freely move and create a dynamic topology (Springel 2010). Dark matter (DM), star, and MBH particles are superposed on the mesh, where MBHs of mass 1.4 × 10⁵|$\rm M_{\odot }$| are seeded in the centre of haloes with total mass >7.1 × 10¹⁰ |$\rm M_{\odot }$| that do not already have a MBH (Sijacki et al. 2015). The MBHs can merge, accrete gas, and impart AGN feedback energy and momentum to their surroundings. Vogelsberger et al. (2013) describe the details of the subgrid MBH accretion and feedback prescriptions, along with the subgrid prescriptions for other physics including star formation, stellar feedback, chemical enrichment, and metal-line cooling. The Illustris simulation reproduces reasonably well many observed properties of galaxies and their MBHs, including the galaxy stellar mass function, stellar luminosity functions, cosmic star formation history, baryon conversion efficiency, the MBH mass function, and AGN luminosity functions (e.g. Vogelsberger et al. 2014a; Genel et al. 2014; Sijacki et al. 2015).

Galaxies in Illustris are identified via the subfind algorithm (Dolag et al. 2009), which finds gravitationally bound ‘subhaloes’ within friends-of-friends haloes. Throughout this paper, any mention of MBH host galaxies in Illustris refers to the subfind subhaloes that host the MBHs.

The Illustris simulations are run in a cosmological box of side |$L_{\rm box}=75\, h^{-1}$| Mpc, from z = 137 to 0. Throughout this paper, we use the highest resolution run, ‘Illustris-1’ (hereafter ‘Illustris’), which has a DM resolution of 6.3 × 10⁶ |$\rm M_{\odot }$| and a baryonic mass resolution of 1.2 × 10⁶ |$\rm M_{\odot }$|⁠. The simulations assume a WMAP9 cosmology with parameters Ω_m = 0.2865, Ω_Λ = 0.7135, σ₈ = 0.820, and H₀ = 70.4 km s⁻¹ Mpc⁻¹ (Hinshaw et al. 2013).

As in SB21, we use the Illustris simulation for our analysis rather than the more recent IllustrisTNG simulations² (Nelson et al. 2018, 2019; Springel et al. 2018; Marinacci et al. 2018; Pillepich et al. 2018b; Naiman et al. 2018 ; Pillepich et al. 2019) because this enables us to use the existing MBH binary inspiral data from Blecha et al. (2016) and Kelley et al. (2017a). These inspiral models were obtained by extracting stellar, gas, and DM density profiles for the hosts of all merging MBHs and integrating the hardening rates from the initial to final separation for each binary. By using Illustris for this study, we are also able to compare directly with results of previous binary inspiral studies based on this simulation (Blecha et al. 2016; Kelley et al. 2017a, 2018; Katz et al. 2020, SB21). Although IllustrisTNG has been shown to better reproduce numerous stellar and gas properties of galaxies (Pillepich et al. 2018a, b), both Illustris and IllustrisTNG produce MBH populations that broadly agree with empirical constraints (Sijacki et al. 2015; Pillepich et al. 2018a, b; Habouzit et al. 2021). Aside from this the supermassive blackhole (SMBH) seeding in IllustrisTNG is higher than the Illustris by a factor of ≈8 (Vogelsberger et al. 2014b; Pillepich et al. 2018a). This, with our prescription for SMBH population delineated in the next section, would exclude more SMBHs within LISA sensitivity range. The IllustrisTNG simulation suite does, however, include both higher resolution (50 Mpc)³ and lower resolution (300 Mpc)³ volumes in addition to the fiducial (100 Mpc)³ volume, which would enable more detailed studies of the inner structure of galaxies and the highest mass MBH populations, respectively. In future work, we plan to generalize our analysis to other data sets.

2.2 Massive black hole binary population

The time at which Illustris records an MBH merger is the starting point for our post-processing binary evolution model. Thus, we refer to the Illustris merger time as the binary formation time. The MBH binary population studied here is identical to that of SB21. Both in this paper and in SB21, we include a subset of all 23 708 MBH mergers in the Illustris for our study. Illustris uses an MBH re-positioning scheme that places MBHs at the potential minimum of their host haloes. This is done in order to prevent spurious numerical kicks to the MBHs by nearby particles or cells, which have comparable or even greater mass. However, in some unequal mass mergers, a satellite galaxy can lose its MBH to the more massive galaxy on an artificially short time-scale, and in some cases the satellite can then be re-seeded with a new MBH. This process can happen multiple times and create spurious mergers, particularly for MBHs near the seed mass. In order to mitigate the impact of these issues on our results, we follow previous work and exclude all mergers that have MBHs below 10⁶ |$\rm M_{\odot }$| (Blecha et al. 2016; Kelley et al. 2017a, SB21). We note that the excluded range of MBH masses (10⁵–10⁶|$\rm M_{\odot }$|⁠) aligns with the peak sensitivity of LISA to MBH merger events. However, the simple, massive MBH seeding prescriptions in Illustris and other large-volume cosmological simulations limits their predictive ability for LISA in any case (e.g. Salcido et al. 2016; Katz et al. 2020, SB21).

The binary hardening prescription in Kelley et al. (2017a) that we are using requires the binary to have an associated galaxy in the snapshots before and after the merger. This is required to identify MBH hardening environments to calculate DF, stellar scattering, and disc hardening rates. There are further constraints on the number of each particle type in a galaxy for it to be resolvable. We require a minimum of 80 gas cells, 80 star particles, and 300 DM particles for a galaxy to be resolved. A combination of all these constraints reduces the total binary population that we are studying to 9234 systems.

2.3 Hardening prescription

We follow the same hardening prescription as in SB21. This prescription is outlined in Kelley et al. (2017a, b) and implements the four different hardening mechanisms (Begelman et al. 1980) of DF, stellar scattering also knows as LC scattering, CD hardening, and hardening due to emission of GWs. Hardening rates for each are calculated using inward extrapolations of spherically averaged stellar, DM, and gas density profiles. We will briefly go over them and how they are calculated here; we refer the reader to Kelley et al. (2017a, b), and SB21 for more details.

At large binary separations (∼kpc), DF is the only effective binary hardening mechanism (Fig. 2, shaded regions, e.g. Antonini & Merritt 2012; Kelley et al. 2017a). At separations of ∼100 to ∼0.1 pc, LC stellar scattering drives the hardening of the binary. In this paper, we consider two implementations of hardening in the LC regime:

• Variable LC refilling rate, zero eccentricity.

  This model allows for systematic variation of the LC refilling efficiency. For variations of the LC refill fraction in our model, we choose circular orbits (e = 0) for all the binaries. The LC refilling efficiency is based on dynamical models of stellar scattering rates from Magorrian & Tremaine (1999). In a ‘full’ LC, low-angular-momentum stars are replenished as fast as they are scattered out of the LC. If instead, two-body relaxation is the only mechanism by which the LC is replenished (on time-scales ≫t_H) then we have a ‘steady-state’ LC. Our models consider these two scenarios for circular orbits, following Magorrian & Tremaine (1999), as well as intermediate scenarios (a partially full LC). We use the logarithmic ‘refilling fraction’ |$\mathcal {F_{\rm refill}}$|⁠, to interpolate between the stellar flux for a full LC (⁠|$F^{\rm full}_{\rm LC}$|⁠) and that of a steady-state LC (⁠|$F^{\rm eq}_{\rm LC}$|⁠):

  $$\begin{eqnarray} F_{\rm LC}\equiv F^{\rm eq}_{\rm LC}\times \Bigg (\frac{F^{\rm full}_{\rm LC}}{F^{\rm eq}_{\rm LC}}\Bigg)^{\mathcal {F_{\rm refill}}}. \end{eqnarray}$$

  (2)

  Following SB21, we adopt |$\mathcal {F_{\rm refill}}=0.6$| as the fiducial value, corresponding to a partially full LC and we compare to results for |$\mathcal {F_{\rm refill}}=0$| and 1 in Section 3.3.

• Variable binary eccentricity, full LC.

  This model allows for the eccentricity of the binary to evolve over time and is based on the scattering experiments by Sesana, Haardt & Madau (2006). The LC hardening rate da/dt and eccentricity evolution de/dt are:

  $$\begin{eqnarray} \left(\frac{\mathrm{ d}a}{\mathrm{ d}t} \right)_{\rm LC} \equiv -\frac{G\rho }{\sigma }a^2 H , \end{eqnarray}$$

  (3)

  $$\begin{eqnarray} \left(\frac{\mathrm{ d}e}{\mathrm{ d}t} \right)_{\rm LC} \equiv \frac{G\rho }{\sigma }a HK . \end{eqnarray}$$

  (4)

  Here, a is the binary separation e is the binary eccentricity. ρ and σ are stellar density and stellar velocity dispersion profiles calculated from the corresponding host in the Illustris simulation. H and K are dimensionless constants determined by numerical scattering experiments (we use the values from Sesana et al. 2006). In all realizations of the binary population using this model, the initial eccentricity is set to the same value for all binaries; we consider models with initial eccentrities of 0.1–0.99.

  The GW hardening rate and eccentricity evolution are (Peters 1964):

  $$\begin{eqnarray} \left(\frac{\mathrm{ d} a}{\mathrm{ d}t}\right)_{\rm GW}=-\frac{64G^3}{5c^5} \frac{m_1 m_2 \left(m_1 +m_2\right)}{a^3} \frac{(1+ 73 e^2/24+ 37e^4/96)}{(1-e^2)^{7/2}},\\ \end{eqnarray}$$

  (5)

  $$\begin{eqnarray} \left(\frac{\mathrm{ d} e}{\mathrm{ d}t}\right)_{\rm GW}=-\frac{304G^3}{15c^5} \frac{m_1 m_2 \left(m_1 +m_2\right)}{a^4} \frac{(e+ \frac{121}{304}e^3)}{(1-e^2)^{5/2}} . \end{eqnarray}$$

  (6)

  Note that in the LC phase, the eccentricity increases as the binary hardens, while in the GW phase binary rapidly circularizes. Note also that eccentricity does not evolve in the disc phase (CD) in our model.

Figure 2.

The shaded regions in this plot show the inspiral time-scale versus separation for the MBH binaries in our fiducial model, with each colour corresponding to a different hardening mechanism as shown in the legend (from SB21, fig. 2). The scatter plot shows a_triple and the corresponding total inspiral time for the subset of binaries that become triples (specifically, these data correspond to the outer binaries in each triple system). The scatter points are colour-scaled to the logarithm of the mass ratio q_outer of the outer binary and the size of each point correspond to the total mass of the outer binary. The top histogram shows the distribution of a_triple for the triple population for five different subpopulations: all triples, triples with |$\rm q_{\rm outer}\lt 0.1$|⁠, M_tot > 3 × 10⁸, > 3 × 10⁷, and > 3 × 10⁷ |$\rm M_{\odot }$| and  q_outer > 0.1, and triples with M_tot > 3 × 10⁸ M_⊙. Notice the bimodality in the a_triple distribution. This bimodality constitutes the basis of our definition of ST and WT interactions.

2.4 Identifying triples

Because our post-processing binary inspiral model evolves each MBH binary in isolation, we must develop a method for identifying which binary systems involving a common MBH co-exist in time, and which of these co-existing systems are likely to become true MBH triple systems. The left panel in Fig. 1 shows a schematic of a merger tree from Illustris and a schematic separation versus time diagram for the potential triple systems. The highlighted regions indicate potential triples. In order to identify and classify triple systems, we first create a time-ordered array of IDs from the Illustris merger tree. When two MBHs merge in Illustris, one of the MBH IDs survives the merger and is assigned to the remnant MBH. We can see this depicted in Fig. 1. When for instance MBHs with IDs ‘00’ and ‘01’ merge, the remnant is assigned the ID ‘01’. Afterwards, we look at repeated MBH IDs in the merger array. Consecutive mergers involving the same MBH ID with a non-zero overlap in their binary inspiral times are referred to as ‘co-existing’ binaries (recall that the time of the MBH merger in Illustris is the binary formation time in our model). This sample of co-existing binary systems forms the superset of binaries that have the potential to become triple systems.

In some co-existing binary systems, the second binary never overtakes the first. Curve (a) in Fig. 1 shows an example of this case, in which no intersection point (a_triple) is defined. We refer to these systems as ‘failed triples’ (FT). Finally, because we are interested in the population of triple MBH systems in the observable Universe, we further restrict the definition of ‘triples’ to include only those in which the first binary is overtaken by the second before z = 0. Thus, binaries that are overtaken at later times z < 0 are also included in our definition of ‘FT’.

To determine which of these co-existing binaries are most likely to undergo triple interactions, we compare the evolution of binary separation versus time (a, t) for each set of co-existing binaries. For a given set of two co-existing binaries involving a common MBH, we identify the ‘first’ or ‘inner’ binary as the one with the earlier formation time, and the ‘second’ or ‘outer’ binary as the one with the later formation time. The point in time at which the separation of the second binary becomes smaller than the separation of the first binary is defined as the time of formation of the triple (t_triple). The corresponding binary separation is called the triple formation radius a_triple.

We define ‘triples’ any such co-existing binaries in which the first is overtaken by the second before z = 0, such that a_triple is a defined quantity. The schematic curves (b) and (c) in Fig. 1 illustrate two examples of triple formation at small and large binary separations, respectively. We refer to case (b) as a ‘strong triple’ (ST) and case (c) as a ‘weak triple’ (WT), where 100 pc is the critical separation between the two scenarios. STs are of particular interest, as these are the systems in which all three MBHs are most likely to be in close proximity and undergo a dynamical interaction. The motivation for choosing 100 pc for this definition is discussed further in Section 3.1.

We note that these binaries are evolved in isolation in a static background density profile in our models, which is a significant simplifying assumption. For example, we cannot account for situations where a binary is not technically overtaken by a second binary, but where the third MBH may be close enough to have some dynamical influence on the first binary. We also cannot account for merger-induced disturbances to the host galaxy potential that may alter the inspiral time-scales of co-existing binaries. Our primary aim in this study is to characterize the binary systems most likely to become triple systems, and to identify the subset most likely to have ST interactions; a detailed investigation of triple dynamics in time-varying galaxy merger potentials is beyond the scope of this work. However, as discussed below, our findings suggest that these simplifications have a fairly minimal impact on the proportion of binaries that become strongly interacting triple systems.

3 RESULTS

3.1 Characteristics of triple MBH systems

We find that triple MBH systems represent a substantial portion of the binary population. 35 per cent of all binaries in our fiducial model co-exist with a second binary at some point during their inspiral (the superset of all co-existing systems defined above). Most of these co-existing systems (22 per cent of the total binary population) go on to become true triple systems, in which the first binary is overtaken by the second binary before z = 0. As expected, these triple systems disproportionately occur when the first binary is stalled. In SB21, we found that 53 per cent of the binaries in our fiducial model are ‘unmerged’ or ‘stalled’ – that is, they do not merge before z = 0. Based on this fiducial model, 71 per cent of triple interactions involve binaries that would not otherwise merge, demonstrating a significant potential for triple systems to drive these unmerged binaries to merger. The remaining 29 per cent of triple interactions involve binaries that would merge before z = 0 according to our fiducial model, but may instead have a different outcome when triple dynamics are considered (e.g. a hastened merger, or a merger between the intruder MBH and a member of the first binary). The proportion of unmerged and merged binaries from the fiducial model that end up being triple systems is 30 per cent and 13 per cent, respectively. These triple systems, if successful at accelerating the inner binary hardening through either K-L oscillations or three-body scattering, could drive some of the unmerged binaries to merger before z = 0 and increase the total merger rate.

Fig. 2 shows a scatter plot of total inspiral time-scale for triples at the triple formation radius a = a_triple (i.e. the radius at which the second binary overtakes the first, according to our binary inspiral models). These data are overplotted on binary inspiral time-scales versus separation for different hardening mechanisms, for all binaries (co-existing or not). Looking at the scatter plot, we can see that the triple formation radii exhibit a strong bimodality: there are distinct populations of triples that occur at either larger separations (∼ kpc) or smaller separations (∼ pc). As noted in Section 2.4, we adopt a cut-off of 100 pc between these two populations, and we designate the former (overtaken at >100 pc separations, with a median a_triple = 992 pc) as ‘WT’ systems. The latter (overtaken at <100 pc) are designated ‘ST’ systems – these are the triple systems most likely to undergo strong three-body interactions. Most triples (∼74 per cent) form at large separations of a_triple ≳ 100 pc. The ST systems, in contrast, have a median formation radius of a_triple = 0.9 pc.

The choice of the 100 pc cut-off between WT and ST formation radii is validated by its simplicity and its statistical consistency. We performed K-means clustering (Macqueen 1967) on the data in Fig. 2 with features (q, M_total, a) and achieved very similar results to a simple 100 pc cut-off. The largest drop in the K-means inertia, which is an indicator of the compactness of the cluster (smaller inertia indicates a more compact cluster and vice versa), happens for two clusters, and the boundary between the clusters corresponds to ≈ 100 pc. The median a_triple for WT and STs from the clustering analysis is 980 and 0.7 pc, respectively; these values are essentially identical to those obtained above. Therefore, we use the 100 pc threshold for ST and WT classifications henceforth.

In our fiducial model, STs and WTs constitute 6 per cent and 16 per cent of the total binary population and 16 per cent and 47 per cent of the co-existing binary population, respectively; an additional 13 per cent of all binaries are classified as FTs (see Table 1) and are discussed further in Section 3.2. Because the intruder MBH in STs must evolve to much smaller separations before overtaking the inner binary, STs have a much longer delay between the formation of the second binary and the time when the intruder MBH overtakes the first binary: the median time for the intruder MBH in STs to overtake the inner binary is 2.8 × 10⁹ yr. In contrast, WTs form with a vastly shorter median delay time of 6.3 × 10⁷ yr. The corresponding median formation redshifts for STs and WTs are z = 0.7 and 1.1, respectively.

Table 1.

Summary of the categorization of all co-existing binaries, as defined in the text. Co-existing binaries are classified as either triples or FTs. Triples are further categorized into STs and WTs. For FTs, we also show the subsets that nearly miss forming an ST or a WT. For each of these two categories and four subcategories, we report the percentage of all binaries and the percentage of all co-existing binaries that fall in that (sub-)category. We also report the percentage of triples or FTs in each subcategory, and the percentage of each triples (sub-)category for which the inner or outer binaries have a mass ratio q ≥ 0.1.

(Sub-)Category	Per cent of	Per cent of	Per cent of	Per cent with	Per cent with
	all binaries	co-existing binaries	triples or FTs	qin ≥ 0.1	qout ≥ 0.1
Triples	22	63	100	58	64
STs	6	16	74	21	17
WTs	16	47	26	37	47
FTs	13	37	100	–	–
FTs: near-miss STs	1.3	3.7	10	–	–
FTs: near-miss WTs	3.8	11	30	–	–

At the time of formation of the triple, the median total mass of inner binaries is 6.5 × 10⁷ M_⊙, while the intruders have a median mass of 1.0 × 10⁷ M_⊙. The inner binaries have time to accrete gas between the time of binary formation and the time of triple formation and the quoted inner binary mass takes that into account. The median binary masses for the inner binary are similar for WTs and STs: |$M_{\rm bin}=7.0\times 10^7\,\mathrm{ and}\,\,5.8\times 10^7$||$\rm M_{\odot }$|⁠, respectively. Comparatively the median binary mass for the total binary population, without taking into account the accretion from the time of the binary formation until the time of triple formation for triple MBHs, is M_bin = 3.0 × 10⁷|$\rm M_{\odot }$|⁠. In addition, a significant number of the binaries (58 per cent of inner binaries and 64 per cent of outer binaries) are major mergers with q ≥ 0.10 (see Table 1), especially those that form STs. The median inner binary mass ratios are q_inner = 0.10 for WTs and q_inner = 0.33 for STs. For the outer binary, the mass ratios (q_outer) are 0.19 and 0.21 for WTs and STs, respectively.

The size of the points in Fig. 2 corresponds to the total mass of the outer binary, while the colour scale of the points corresponds to the logarithm of the outer binary mass ratio. Additionally, the top histograms in Fig. 2 show the distribution of a_triple for all triples, as well as the subsets in which the outer binary is a major merger (with mass ratio q_outer > 0.1) or has a large total mass of M_tot > 3 × 10⁸|$\rm M_{\odot }$|⁠. We see that most of the triple systems in which the first binary is a minor merger q < 0.1 with a high total binary mass tend to form at large separations. Both the WT and ST subsets span a wide range of inner binary mass ratios and total masses; however, nearly all of the extremely massive systems are in the WT regime and have very low mass ratios. This increases the DF time-scales and leads to a larger triple formation radius. In such cases, if the intruder MBH is also much less massive than the central MBH, the WT may consist of two satellite MBHs orbiting in the massive host halo for more than a Hubble time. Conversely, if the intruder MBH is comparable to or more massive than the central MBH (q_outer ≳ 1), its evolution to form a binary with the central MBH will likely be unhindered by the presence of the much less massive, wide-separation member of the first binary.

Fig. 3 shows a 2D histogram of the mass ratio of the outer binary (q_outer) versus the inner binary (q_inner) for all triple MBHs on the left-hand side and the corresponding total subhalo mass ratios of their hosts on the right-hand side. The q_inner and q_outer mass ratio distributions are fairly similar, but with a large scatter in q_inner versus q_outer for individual systems. On average, though, q_outer in each triple system tends to be slightly higher than q_inner. This may partly reflect the fact that the intruder MBH interacts with the system at a later time and thus has more time to accrete before it enters the system. Another contributing factor can be the fact that inner binaries with low-mass ratios (i.e. small q_inner) in general have longer merger time-scales, and thus are more likely to be overtaken than inner binaries with equal mass ratios (Kelley et al. 2017a, SB21).

Figure 3.

Mass ratio of the outer binary MBH (y-axis) and the inner binary MBH (x-axis) is shown in the left here. On the right, we have the mass ratio of the corresponding subhaloes. The white dashed lines indicates the one-to-one ratio. We see that the q_outer are slightly larger compared to q_outer for the MBHs. For the subhalo mass ratios, we see a similar trend, except that a subset of intruder MBH hosts are very massive relative to the inner binary hosts.

We also examine the mass ratios of the host galaxies. We define q_{in, host} and q_{out, host} in the same manner as q_inner and q_outer, but using the total subhalo masses of the hosts instead of the MBH masses. These host mass ratios largely show a similar trend to the MBH mass ratios, in which the outer binaries in general have higher mass ratios compared to inner binaries. The host mass ratio plot also shows a bimodality in the distribution of outer host mass ratios q_{out, host} compared to q_{in, host}. The very high mass ratio systems with q_{out, host} > 10 constitute |$\sim 10{{\ \rm per\ cent}}$| of the triple population. Most (⁠|$\sim 70{{\ \rm per\ cent}}$|⁠) of these very high mass ratio systems are WT systems.

3.2 Failed triple systems

We now examine in more detail the population of FTs – the ∼1/3 of co-existing MBH binaries that fail to become triples according to the above definitions. The FT MBHs (1170 systems, or 13 per cent of all binaries) consist of two subpopulations. Most of them (79 per cent) are systems in which the inner binary merges before it can be overtaken by the outer binary. In the rest of the FT systems, the outer binary does eventually overtake the inner binary, but not before z = 0. Focusing on the former category, we consider whether some of these FTs may have come close to forming a successful triple. Given that our models include only binary evolution in isolation, whereas triple dynamics will be more complicated, it is possible that some such ‘near misses’ could in reality become triple systems.

We compare the minimum ratio of the outer and the inner binary separations, (a₂/a₁)_min, in order to see how close the outer binary came to overtaking the inner. Of greatest interest are the systems in this category that have a small binary separation ratio (1 < (a₂/a₁)_min ≤ 10) that also have a close inner binary at the point of closest approach (a₁ < 100 pc). These are the FTs that came closest to forming STs without the intruder actually overtaking the inner binary. Only 120 binaries are in this category (10 per cent of all FTs, or 1.3 per cent of all binaries). Given that 6 per cent of all binaries do form STs in our fiducial model, these near misses would at best provide a small but not insignificant enhancement to the ST population, if most of them were actually able to form successful triples. We note that an additional 1.1 per cent of all binaries have a₁ < 100 pc and 10 < (a₂/a₁)_min ≤ 100. Systems that nearly missed forming WTs (1 < (a₂/a₁)_min ≤ 10, while a₁ > 100 pc) are more common, representing 30 per cent of all FTs, or 3.8 per cent of all binaries.

3.3 Dependence on binary inspiral model parameters

Here, we explore the dependence of the triple MBH population characteristics on two key binary inspiral parameters: the stellar LC refilling efficiency |$\mathcal {F_{\rm refill}}$| (as defined in equation 2) and binary eccentricity. Recall that we start the binary hardening process by assigning the same initial eccentricity to all binaries, after which the eccentricity evolves through the LC and GW stages of evolution according to equations (4) and (5). Note that the DF and CD phases of evolution are unaffected by variation in these two binary parameters.

Fig. 4 shows how the various triple subsets depend on the LC refill rate and the initial binary eccentricity. These plots illustrate an important finding of our analysis: the number of triple systems does not depend very strongly on either |$\mathcal {F_{\rm refill}}$| or eccentricity. For WT systems, this is essentially true by definition: at separations >100 pc, binary hardening is not dominated by LC scattering or GW emission, so neither |$\mathcal {F_{\rm refill}}$| or eccentricity has an impact on the occurrence of WTs. Less intuitive is the similar number of ST systems across a wide range in |$\mathcal {F_{\rm refill}}$| and eccentricity values.

Figure 4.

The dependence of the triple population on the LC refilling efficiency (⁠|$\mathcal {F_{\rm refill}}$|⁠, left panel) and binary eccentricity (right panel) is shown here. In both panels, all co-existing binaries are indicated in black (circles for each data point); they are the superset of all other triple populations. STs (dashed magenta with a cross for each data point) and WTs (dashed cyan with a star for each datapoint) constitute the total triple (solid blue with triangles for each datapoint) population. FT (red with squares for each datapoint) are the systems in which the second co-existing binary either never overtakes the first (i.e. the intruder MBH ‘chases’ the inner binary), or it fails to do so before z = 0. Roughly 1/3 of all MBH binaries co-exist at some point with another binary (i.e. their inspirals overlap in time), and most of these become triple systems. A few percent of all MBH binaries form STs at small separations. These results are remarkably insensitive to variation in binary inspiral model parameters.

We can understand the lack of dependence on LC refilling efficiency as a balance between two competing effects: moderate LC refilling can slightly enhance the triple population by bringing in third MBHs more quickly, but in the full LC regime, this is more than offset by the increased likelihood that the inner binary will merge before being overtaken by a third BH. None the less, we find that the population of STs is remarkably robust – at a few percent of the total binary population – whether the LC is full, empty, or somewhere in between.

When |$\mathcal {F_{\rm refill}}$| is increased from zero (a steady-state, empty LC) to 0.6 (a partially full LC), the total number of co-existing binaries decreases slightly, because shorter inspiral time-scales increase the likelihood that the first binary will merge before the second is formed. However, the triple fraction goes up very slightly (by half a percentage point, from 21.6 per cent to 22.1 per cent) over the same range in |$\mathcal {F_{\rm refill}}$| values. This arises from systems in which the increased LC refilling efficiency has an outsized effect on the outer binary relative to the inner binary. In other words, some systems that are FT in an empty LC model become successful triples when moderate LC refilling is assumed.

In contrast, when the LC is always full (⁠|$\mathcal {F_{\rm refill}}=1$|⁠), the faster binary hardening significantly reduces the chance of co-existence. While 38.5 per cent of all binaries co-exist with another binary during their inspiral in the |$\mathcal {F_{\rm refill}}=0$| model, only 27.5 per cent of binaries co-exist in the full LC model. This reduces the fraction of binaries that form ST systems as well: 5.2 per cent and 5.8 per cent of all binaries form STs for |$\mathcal {F_{\rm refill}}=0$| and |$0.6$|⁠, respectively, versus 1.9 per cent of all binaries that form STs for |$\mathcal {F_{\rm refill}}=1$|⁠.

Fig. 4 reveals a similarly weak dependence of triple occurrence on binary eccentricity. Again, the binary eccentricity is not expected to have any effect on WTs in our model, but even for STs, the ST fraction varies by less than 1 per cent when the initial binary eccentricities are varied from e = 0.1 to 0.9. This reflects the fact that eccentricity influences binary hardening rates predominantly in the GW phase, which represents a small portion of the total inspiral time. Only at truly extreme initial eccentricities of e = 0.99 do we see a noticeable drop in the occurrence of triple systems, owing to the strong dependence of GW-driven hardening on eccentricity (equation 5). At very high eccentricities, strong GW emission at pericentre can greatly reduce the merger time-scale and drive the inner binary to merger before it is overtaken by an intruder MBH.

The redshift distribution of binaries, triples and their strong and weak subpopulation are shown in Fig. 5 with the total numbers indicated on the left panel and the fractions on the right panel. The WTs fraction increases over time at the expense of the ST population. The rich environment of early galaxies along with higher merger rates could contribute to the higher proportion of STs.

Figure 5.

The evolution of all binaries and various triple subpopulations are shown here as a function of redshift for the fiducial model. The left plot shows the total numbers and the right plot shows the triple populations as a fraction of the total number of binaries. From the left panel, we can see that the triple population increases as the total binary population increases, as expected, and turns over at z ≲ 1 as binary mergers begin to outweigh new binary formation. The fraction of triples remains fairly flat over time, ∼0.2–0.3. Surprisingly, STs initially outnumber WTs, but this trend reverses at z ∼ 3.7. At redshift close to zero, we see that the majority of triple systems are WTs.

Overall, we see that the triple population is not strongly dependent on the specific hardening model we choose. This stands in stark contrast to the z = 0 population of unmerged binaries, which varies widely depending on the binary inspiral model parameters assumed (most notably the LC refilling efficiency; Kelley et al. 2017a, b, SB21). The fraction of binaries that evolve to become triple systems, however, remains essentially unchanged across a wide range of model variations, except for very high eccentricities and full LCs. Our models therefore predict that a robust ∼ few percent of all MBH binaries will form ST systems. Notably, these results suggest that the fraction of binaries that evolve to become triples is largely determined by the galaxy merger history, rather than the details of the astrophysical environment driving binary inspiral. This underscores the importance of using full cosmological hydrodynamics simulations instead of semi-analytic models to make these predictions, given that these methods have been shown to produce significantly different trends in the galaxy merger rates with redshift and stellar mass (Rodriguez-Gomez et al. 2015).

3.4 Host galaxy properties

Clearly, we can expect that galaxies that host triple MBH systems will generally live in richer cosmic environments than those that do not undergo multiple mergers. Here we examine whether, in addition, the triple MBH host galaxies themselves may have distinctive features. In Fig. 6, we look at key galaxy parameters for all BH binaries in Illustris, along with the subpopulations of first and second binaries in co-existing systems, as well as the isolated binaries that never co-exist with another binary, which comprise 65 per cent of the binary population in our fiducial model.

Figure 6.

Distributions of several properties of MBH binary host galaxies are shown: binary formation redshift (top left panel), sSFR (top right), total galaxy mass M_subhalo (middle left), stellar mass M_* (middle right), and stellar velocity dispersion (bottom left). The black dashed lines show the histograms for the hosts of all 9234 MBH binaries (dashed line). The solid-line histograms show the subpopulations of isolated binaries that never co-exist with another binary (grey), the first binary in each co-existing system (magenta), and the second binary in each system (cyan). Only modest differences are seen between these subpopulations, which reflect the fact that the hosts of first binaries, forming by definition at earlier times, have characteristics more typical of higher redshift galaxies. Note that the discontinuity in the redshift corresponds to two corrupted Illustris snapshots at z ∼ 4 that are excluded from all analysis and do not affect our conclusions.

Fig. 6 shows the host properties for each of these subpopulations in the simulation snapshot immediately following the formation of the MBH binary. Overall, we see that the distributions of host galaxy properties look very similar for isolated binaries versus all binaries, as expected given their dominance of the population. There is, however, a noticeable difference between the first and second binaries in co-existing systems. By definition, the inner binary (first binary) in a co-existing system forms at slightly higher redshift than the outer binary (second binary). This is reflected in the minor differences between their host galaxy properties, such that on average, the hosts of first binaries have lower stellar masses, lower stellar velocity dispersions, and higher specific star formation rates (sSFR, defined as the star formation rate divided by the stellar mass).

Fig. 7 shows the host properties for the first binaries (left) and the second binaries (right) in co-existing and triple systems. We see that the host properties are broadly similar across the subcategories of STs, WTs, FTs, and all co-existing systems. There is a slight shift in the redshift distribution, as the second binaries form later in time, and the stellar mass of the second binary peaks at a higher value. Aside from that all the distributions for the subpopulations of WTs, STs, FTs, and all co-existing systems are qualitatively similar.

Figure 7.

Distributions of host galaxy properties are shown for the first binary (left) and second binary (right) in each co-existing binary system, in the simulation snapshot immediately after the first binary forms. Host properties are plotted in the same manner as in Fig. 6. Here, the blue dashed lines show the host distributions for the each of the first and second binary in each co-existing system, while the solid-line histograms show the host distributions for each of the triple subpopulations: STs (orange), WTs (green), and FT (red). While the distributions are fairly similar overall, the progenitors of STs form at slightly higher redshift in galaxies with slightly lower masses and higher sSFRs. Note that as mentioned before, the discontinuity in the redshift is due to two corrupted Illustris snapshots at z ∼ 4.

Binaries that later become STs after being overtaken by an intruder do tend to form at slightly higher redshift. These ST progenitors also have somewhat lower masses and velocity dispersions and slightly higher sSFRs. Therefore, while none of these differences are large, the host properties that distinguish the first binaries in triple systems from other binaries are even more pronounced for the hosts of first binaries that go on to form STs. This can be understood in part by recalling that WTs have somewhat lower median inner mass ratios than STs (q_inner ∼ 0.1 versus q_inner ∼ 0.3), and that nearly all of the triples with very low mass ratios (q_inner < 0.01) are WTs. Galaxies at z ≳ 1–2 will on average have lower masses (and thus higher minimum mass ratios), higher gas content (yielding higher SFRs and more efficient binary inspiral), and higher galaxy merger rates (such that an intruder is more likely to overtake the first binary before it merges). In addition, because STs typically have a long delay between the first binary formation and the triple formation, some lower redshift binaries that might have become STs are likely to end up instead as FT, because the intruder is unable to overtake the first binary before z = 0.

4 DISCUSSION

In this paper, we have identified and characterized a population of candidate triple MBH systems in the Illustris cosmological simulation, using a new classification scheme that we have developed. We first assign subresolution binary inspiral time-scales to each merging MBH pair from the simulation, using the binary inspiral models of Kelley et al. (2017a). When a given MBH is involved in two successive binary inspirals that co-exist in time, this system is identified as a possible triple MBH. If at any point in time the separation of the second binary becomes smaller than the separation of the first binary, this system is identified as a successful triple MBH. Because this analysis is based on MBH binary models that assume each binary evolves in isolation, these inspiral time-scales do not incorporate the dynamical effects of a third, intruder MBH. We also do not account for any MBH mass growth that may occur during the binary inspiral (cf. Siwek, Kelley & Hernquist 2020). None the less, our findings predict that a non-negligible fraction of MBHs are involved in triple systems at some point in their evolution and that this population is relatively insensitive to assumed binary inspiral model parameters. Our work also highlights some key characteristics of this putative triple MBH population.

We find that more than a third (35 per cent) of all MBH binaries in Illustris do co-exist with another binary, according to our fiducial inspiral model. Most of these (and a quarter of all Illustris MBH binaries) are in fact overtaken by another MBH at some point, and in nearly all such cases (22 per cent of all binaries), the MBH binary is overtaken before z = 0. Our findings are consistent with Kelley et al. (2017a), where they used a similar method to estimate that ∼30 per cent of all Illustris MBH binaries could form triple systems. Similar results were also obtained from the detailed semi-analytic models of Bonetti et al. (2018b, 2019), which indicated that triple interactions could be an important channel for low-frequency GW sources observable with PTAs and LISA.

Another key finding of this work is that the triple MBH population is bimodal. Triple systems can be clearly divided into two main subcategories, which we refer to as ‘STs’ and ‘WTs,’ depending on the binary separation at which the first binary is overtaken by the intruder MBH (defined as the radius of triple formation, or a_triple). The WTs are defined as those that form at large separations, a_triple > 100 pc, with a median of a_triple ≈ 990 pc, and these constitute 16 per cent of the total binary population. STs (a_triple < 100 pc) form at a median separation of a_triple ≈ 0.9 pc, and they constitute 6 per cent of the total binary population. This bimodality can be understood as the result of two different mechanisms that can lead to binary stalling: inefficient DF for inspiraling binaries at large separations, or inefficient stellar LC scattering at small separations. If indeed gaseous CDs do at times increase binary separations rather than drive efficient inspiral, this process may contribute to binary stalling at small scales as well.

We find that STs typically require a lengthy delay of a few Gyr to evolve from the formation of the second binary to the point where the intruder MBH overtakes the first binary. In contrast, WTs (which form at much larger binary separations) have a much shorter typical delay of a few × 10⁷ yr between second binary formation and triple formation. As a result, STs form at lower median redshifts than WTs (z = 0.7 versus 1.1). STs typically have a smaller total binary mass (⁠|$M_{bin}$|⁠) and higher inner mass ratios as compared to the WTs. The difference between the mass ratios for the outer binary is less significant.

Interestingly, we also find that the occurrence rate of triple systems is relatively independent of the binary inspiral parameters in the LC and GW phases. Specifically, we vary the initial binary eccentricity and the LC refilling rate, which by definition in our model do not affect binary evolution in the DF phase. This stage of binary evolution corresponds to the separations at which WTs form. Thus, variations of these two parameters (binary eccentricity and LC refilling rate) do not affect WTs at all. Somewhat surprisingly, though, we find that the incidence of ST is also relatively insensitive to these binary model parameters. In essence, this means that more efficient binary inspiral (due to higher eccentricity or faster LC refilling) shortens merger time-scales and triple formation time-scales by roughly similar amounts. We do, however, see a decrease in triple formation when stellar LCs are assumed to be always full: |$\sim 1.9{{\ \rm per\ cent}}$| of all binaries form ST in this case, as opposed to |$\sim 5\,\mathrm{ per\,cent} - 6{{\ \rm per\, cent}}$| for empty LCs and partially refilled LCs. This indicates that when stellar hardening in the LC phase is very efficient, binaries do tend to merge faster than triples can form. But given the wide range in binary inspiral time-scales between full and empty LCs (Kelley et al. 2017a), the influence on triple occurrence is relatively modest.

Variations in binary eccentricity (which impact evolution in the LC and GW phases) have even less influence on the triple population; only for the extreme model with e = 0.99 do we see an appreciable drop in the formation rate of ST. Note that we do not consider any eccentricity evolution in the DF or gas-driven binary inspiral phases. The CD phase in particular depends on the MBH accretion rates (SB21), so a variation in the accretion model could also change these results. These questions will be explored further in future work.

Regardless, because we find that a consistent ∼ few percent of all MBH binaries are involved in an ST interactions, we can conclude that the occurrence of triple MBH systems is closely tied to the rate of MBH binary formation and by extension to the galaxy merger rate. This highlights the importance of accurate measurements of the galaxy merger rate as an essential foundation for constraining MBH dynamics and GW source populations.

One of the many complexities of triple MBH dynamics is that the galaxy nucleus containing the first binary may be significantly altered by the subsequent galaxy merger that introduces the intruder MBH. Merger-driven gravitational instabilities can perturb the existing stellar distribution and can also drive cold gas inflows, which may trigger new nuclear star formation. Thus, if the first binary’s inspiral had stalled due to an empty stellar LC, a subsequent galaxy merger would provide a possible means for partially refilling the LC on roughly a dynamical time-scale. This could in theory shorten the merger time-scale of the first MBH binary, before the intruder MBH reaches the nucleus to undergo an ST interaction. In some cases, this might even mean that the first MBH binary could merge before a triple system was able to form. As we have noted throughout, a full treatment of triple MBH dynamics is beyond the scope of this work. However, our finding that LC refilling has only a modest impact on the ST population can be used to constrain the influence of galaxy merger dynamics on triple formation.

The most extreme example of the above scenario would be one in which a stalled binary in an initially empty LC is quickly refilled by galaxy merger dynamics (i.e. perturbations to the stellar nucleus, or formation of a new stellar cusp). We can estimate the impact of such a scenario on the triple population by comparing the two extremes of our LC refilling model variation: |$\mathcal {F}_{\rm refill} = 0$| and |$1$|⁠. Fig. 2 shows that the proportion of all binaries that form ST drops from |$\sim 5{{\ \rm per\ cent}}$| to |$\sim 2{{\ \rm per\ cent}}$| when LCs are assumed to be completely full instead of completely empty. Therefore, the impact of galaxy merger dynamics on LC refilling (and thus MBH merger time-scales) would have to be smaller than this. In other words, our results indicate that merger-driven LC refilling could reduce but not decimate the triple MBH population.

We additionally note that, because this study models the evolution of each binary in isolation, some systems that form an FT can eventually develop into either WTs or STs. And some systems that initially form as WTs could potentially develop into STs at a later time. This can happen if the (t, a(t)) components of one binary cross that of the other binary multiple times. For example, if the second binary initially has a higher rate of hardening and the intruder MBH overtakes the first binary at a(t) ∼ 1 kpc, a WT system is formed. If this second binary later stalls at ∼ pc separations, the first binary could in theory overtake it again, leading to an ST interaction. Around |$\sim 19{{\ \rm per\ cent}}$| of the triple population experiences more than one crossing. Of course, in reality the evolution of the inner and outer binaries is not isolated, and their dynamics are more complicated.

Analysis of host galaxy properties reveals some mild distinctions between binary and triple MBH hosts. The first binary in each triple system forms at higher redshift than the second, by definition. Accordingly, the median total and stellar host masses are ∼2–3 times lower when the first binary forms than when the second binary forms. The hosts of the first binaries also have ∼2 times higher sSFR. Among the first binaries in each triple, those that eventually undergo strong interactions with the intruder MBH form at slightly higher redshifts than those that become WTs. The higher major merger rates of high-redshift galaxies are favourable for the formation of triple systems (note that a significant tail of STs form at z < 1). These hosts of STs therefore tend to have lower masses and higher sSFRs than the hosts of WTs. None of these trends in host galaxy properties are dramatic, however, and their distributions all have substantial overlap.

While the WT systems are too widely separated to allow interactions between MBHs, they have an advantage in terms of observability, because their constituent MBHs could be spatially resolved. If more than one of these MBHs is actively accreting, they may appear as a dual or triple AGN system. Such objects are of great interest, in part because they provide information about the early stages of MBH binary formation and inspiral. Moreover, multiple-AGN systems can provide clear examples of AGN fuelling induced by galaxy mergers (Pfeifle et al. 2019). The role of mergers in AGN fuelling is actively debated, with some studies finding a clear merger-AGN connection, and others finding that mergers are not an important channel for triggering AGN (e.g. Ellison et al. 2011; Cisternas et al. 2011; Kocevski et al. 2015; Villforth et al. 2017; Koss et al. 2018). Triple AGN systems, by definition, form in hosts with an active recent merger history; thus, further observations of such systems would offer a unique window into the merger-AGN connection.

Some of the WT systems we identify in Illustris consist of a massive primary MBH, a much less massive secondary MBH, and an intruder MBH that may or may not be massive. In such cases, the low-mass MBH(s) are likely to have long inspiral time-scales and may be difficult to detect unless they are accreting at a substantial fraction of the Eddington rate. Future observational constraints on the prevalence of dual and triple AGN systems, including those with small MBH mass ratios, will enable further studies of the WT population. In addition to revealing key information about the nature of AGN fuelling in galaxy mergers, discoveries of these systems will also provide implicit constraints on the ST population.

The STs can lead to the rapid merger of any of the binary members of the triple, which may greatly reduce the merger time if the binary inspiral had previously been stalled. In a study of the parameter space of triple MBH interactions, Bonetti et al. (2018b) found that ∼ 20 per cent–30 per cent of binaries that would not otherwise have merged within a Hubble time were driven to merger by triple interactions. In such cases where the coalescence happens before z = 0, these triple interactions can increase the MBH merger rates detectable with LISA. Likewise, MBH binaries stalled at ≫ mpc separations will not produce GWs detectable with PTAs, but triple MBH interactions can rapidly harden these binaries and bring them into the GW regime, where they could contribute to the stochastic GW background or even be detectable as continuous-wave GW sources at ∼ nHz frequencies. Meanwhile, in a hierarchical three-body system, the intruder MBH can trigger K-L oscillations, which can rapidly bring the inner MBHs to the GW-dominated regime. This can result in an MBH merger detectable with LISA, possibly preceded by bursts of GW emission at frequencies within the LISA band.

Strong, chaotic three-body interactions can also eject the lightest MBH from the system while further hardening the more massive pair. Thus, in addition to their impact on merger rates, triple interactions can create a population of offset MBHs. Note that a galaxy nucleus still contains one or more MBHs after a slingshot recoil, in contrast to a GW recoil that ejects the merged MBH and leaves behind an empty galactic nucleus. If the remaining MBH pair merges and experiences a large GW recoil kick, however, this could lead to multiple offset AGN and an empty nucleus. Further theoretical studies are needed to understand the relative importance of slingshot versus GW recoil for offset MBH populations and for MBH-galaxy co-evolution.

Our main results can be summarized as follows:

• Significant triple population. 35 per cent of all binary inspirals co-exist when evolved in isolation, and 22 per cent of all binaries form triples, defined as the second binary overtaking the first.

• ST and WT subpopulations. The true triple population can be categorized into two distinct subpopulations of strongly interacting (⁠|$a_{\rm triple}\lt 100\rm pc$|⁠) and weakly interacting (⁠|$a_{\rm triple}\gt 100\rm pc$|⁠) triples based on the radius at which the triple forms. In our fiducial model, the STs and WTs comprise 6 per cent and 16 per cent of the total binary population, respectively.

• Lack of dependence on LC refilling or binary eccentricity. The binary LC refilling parameter has a fairly minor effect on both WT and ST populations. Only for very high values of (F_refill ∼ 1) do we see a noticeable drop in the triple population, which is mostly caused by first binaries merging before they are overtaken by an intruder MBH. Similarly, the triple population and all of its subcategories are remarkably invariant to the initial binary eccentricity. We notice a drop in the triple population only for highly eccentric binaries (e ∼ 1), where very rapid binary evolution drives MBHs to merger before the intruder MBH can intervene.

• Binary and triple MBH host properties. The first binaries and the second binaries in triple systems have slightly different host properties. Hosts of the first binary have lower stellar masses and stellar velocity dispersion and slightly higher sSFR, consistent with their slightly higher formation redshifts. WT and FT have nearly identical distribution of host properties for both the initial binary host and the host for the second merger (which creates the triple system). For STs, on the other hand, the first binary host has slightly smaller subhalo masses and stellar masses. This can be understood as a combination of somewhat larger mass MBH ratios and longer delay times between first binary and ST formations, which favours first binaries that form at slightly higher redshifts.

ACKNOWLEDGEMENTS

This work made use of the python (Rossum & De Boer 1991) programming language along with jupyter notebooks (Kluyver et al. 2016). Visualizations are made with matplotlib library (Hunter 2007). numpy was used for the analysis (Harris et al. 2020). LB acknowledges support from National Science Foundation (NSF) grant AST-1909933, NASA grant 80NSSC20K0502, and from the Research Corporation for Science Advancement via Cottrell Scholar Award 27553.

DATA AVAILABILITY

The underlying data for the Illustris cosmological simulations are publicly available (Nelson et al. 2015). The data underlying the subresolution analysis and MBH binary evolution will be shared on reasonable request to the corresponding author.

Footnotes

1

http://www.illustris-project.org/

2

https://www.tng-project.org