RABBITS – I. The crucial role of nuclear star formation in driving the coalescence of supermassive black hole binaries

In this study of the ‘Resolving supermAssive Black hole Binaries In galacTic hydrodynamical Simulations’ (RABBITS) series, we focus on the hardening and coalescing process of supermassive black hole (SMBH) binaries in galaxy mergers. For simulations including different galaxy formation processes (i.e. gas cooling, star formation, SMBH accretion, stellar and AGN feedback), we systematically control the effect of stochastic eccentricity by fixing it to similar values during the SMBH hardening phase. We find a strong correlation between the SMBH merger time-scales and the presence of nuclear star formation. Throughout the galaxy merging process, gas condenses at the centre due to cooling and tidal torques, leading to nuclear star formation. These recently formed stars, which inherit low angular momenta from the gas, contribute to the loss cone and assist in the SMBH hardening via three-body interactions. Compared to non-radiative hydrodynamical runs, the SMBH merger time-scales measured from the runs including cooling, stellar and SMBH physical processes tend to be shortened by a factor of ∼ 1.7. After fixing the eccentricity to the range of e ∼ 0 . 6 – 0 . 8 during the hardening phase, the simulations with AGN feedback reveal merger time-scales of ∼ 100 – 500 Myr for disc mergers and ∼ 1 – 2 Gyr for elliptical mergers. With a semi-analytical approach, we find that the torque interaction between the binary and its circumbinary disc has minimal impact on the shrinking of the binary orbit in our retrograde galaxy merger. Our results are useful in improving the modelling of SMBH merger time-scales and gravitational wave event rates.


INTRODUCTION
Supermassive black hole (SMBH, with mass in the range of MBH ∼ 10 5 -10 10 M ⊙ ) binaries within galaxy mergers are the primary target for low-frequency gravitational wave (GW) observatories.These include the ongoing nano-Hertz Pulsar Timing Array (see Burke-Spolaor et al. 2019 for a review, and Agazie et al. 2023a,b;Antoniadis et al. 2023;Reardon et al. 2023;Xu et al. 2023 for recent results) and the upcoming milli-Hertz space-borne detectors such as the Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017, 2023), TianQin (Luo et al. 2016), andTaiji (Ruan et al. 2020).Understanding the coalescence process of SMBH binaries in detail is pivotal for predicting SMBH merger time-scales and event rates (Begelman et al. 1980;Amaro-Seoane et al. 2023).
⋆ Email: shliao@nao.cas.cnDuring galaxy mergers, the central SMBHs initially sink to the remnant nucleus due to dynamical friction from stars, dark matter, and gas (Chandrasekhar 1943;Ostriker 1999), forming a gravitationally bound binary.The separation of the binary SMBHs continues to drop rapidly, initially because of dynamical friction and, later on, through three-body interactions with the surrounding stars, also referred to as gravitational slingshot interactions (e.g.Mikkola & Valtonen 1992;Quinlan 1996;Sesana et al. 2006).When the semimajor axis of the SMBH binary reaches the hard binary separation (usually at ∼pc scales), where the specific binding energy of the binary equals the specific kinetic energy of stars, it forms a hard binary.The hard binary can eject nearby stars at high velocities beyond its influence radius (i.e. the radius within which the enclosed stellar mass is twice the binary mass), effectively excluding them from further interactions.The ejections of stars can lead to the formation of a large stellar density core in massive elliptical galaxies (e.g.Milosavljević & Merritt 2001;Milosavljević et al. 2002;Merritt 2006;Rantala et al. 2018;Frigo et al. 2021;Nasim et al. 2021).Further shrinking of the binary orbit necessitates replenishing the loss cone, i.e. the low angular momentum region of the phase space where the stars are prone to interact with the binary (see Merritt 2013, for a review).Whether the two SMBHs of a hard binary can reach small enough separations (i.e.≲ mpc) to enter the GW emission-driven phase (Peters & Mathews 1963;Peters 1964) has important implications for GW observations.
Previous gas-free N -body studies found that a hard SMBH binary might stall at a separation of roughly one parsec due to the depletion and inefficient refilling of the loss cone.This is referred to as the final parsec problem (Milosavljević & Merritt 2003).This phenomenon was seen in idealized spherical systems, where the refilling process occurs through the two-body relaxation mechanism.This gives rise to a hardening time-scale that strongly depends on the number of star particles within the system, i.e. ∼N (Makino & Funato 2004;Berczik et al. 2005).If we extrapolate the result to the real galaxy systems in the Universe, where N is large as it represents the number of actual stars, this time-scale can surpass the Hubble time, implying that the binary could fail to get through the final parsec and thus could not enter the GW-driven regime.
Several possible solutions have been proposed to overcome the final parsec problem bottleneck.Especially, departures from spherical symmetry have emerged as an effective solution to this problem (Merritt & Poon 2004;Berczik et al. 2006;Khan et al. 2011;Preto et al. 2011).Specifically, in triaxial potentials, which are a natural outcome of galaxy mergers, there is a large population of stars with centrophilic orbits which can travel close to the centre and interact with the binary, facilitating the two SMBHs to enter the GW-driven phase.The Brownian motion, i.e. the random wandering of the binary's centre of mass (CoM) triggered by encounters with stars, might also help the binary to interact with more stars (Quinlan & Hernquist 1997;Merritt 2001;Chatterjee et al. 2003;Bortolas et al. 2016).The interaction with a third SMBH, introduced through a subsequent galaxy merger, can excite the binary eccentricity to high values (e.g. the von Zeipel- Kozai-Lidov mechanism, von Zeipel 1910;Kozai 1962;Lidov 1962), substantially reducing the merger time-scale and thus enhancing the SMBH merger rate (e.g.Blaes et al. 2002;Iwasawa et al. 2006;Hoffman & Loeb 2007;Kulkarni & Loeb 2012;Ryu et al. 2018;Bonetti et al. 2019;Mannerkoski et al. 2021).Moreover, the presences of massive perturbers (e.g.Perets & Alexander 2008;Bortolas et al. 2018;Arca Sedda et al. 2019) and galaxy rotation (Holley-Bockelmann & Khan 2015; Mirza et al. 2017) can also contribute to the shrinking of a hard binary.It is important to note that real galaxy systems may experience a combination of these mechanisms concurrently.
Apart from the aforementioned dissipationless gravitational mechanisms, dissipative gas processes are also expected to affect the shrinking of the binary orbit, especially in gas-rich disc galaxy mergers.The gas dynamical friction and torques from circumnuclear discs have been suggested to play an important role in bringing two SMBHs to form a binary (e.g.Armitage & Natarajan 2002;Escala et al. 2004Escala et al. , 2005;;Dotti et al. 2006;Mayer et al. 2007;Chapon et al. 2013;Souza Lima et al. 2017).During the hard binary phase, the SMBH binary further interacts with the gaseous circumbinary disc (CBD, see Lai & Muñoz 2023, for a recent review).Small-scale hydrodynamical CBD simulations have revealed that the gas accretion from the CBD drives the binary toward equal mass (e.g.Artymowicz & Lubow 1996;Hayasaki et al. 2007;Farris et al. 2014;Duffell et al. 2020).Furthermore, the torques from the CBD can lead to either a shrinkage or an expansion of the binary orbit, depending on the disc and binary properties (e.g.MacFadyen & Milosavljević 2008;Roedig et al. 2012;Miranda et al. 2017;Muñoz et al. 2019Muñoz et al. , 2020;;Duffell et al. 2020;Heath & Nixon 2020;Tiede et al. 2020;Franchini et al. 2021Franchini et al. , 2023;;Siwek et al. 2023;Tiede & D'Orazio 2024).Note that the exact effects from CBDs remain a subject of ongoing debate.
In this work, we introduce the 'Resolving supermAssive Black hole Binaries In galacTic hydrodynamical Simulations' (RAB-BITS) series of studies, aiming to investigate the orbital evolution of SMBH binaries down to the GW emission phase in galacticscale hydrodynamical simulations of both disc and elliptical galaxy mergers.In a companion paper, Liao et al. (2023b) (hereafter Paper II), we demonstrate that the central properties of galaxies and the sinking and formation of SMBH binaries can be significantly influenced by galaxy formation processes, especially the active galactic nuclei (AGNs) feedback process and its specific numerical implementations.
In this study (Paper I of the RABBITS series), we focus on the evolution of SMBH binaries during the hardening and GW emission phases.We examine how various galaxy formation processes impact the binary hardening.Previous studies have shown that the eccentricity of an SMBH binary, when it becomes hard, exhibits stochastic behaviour due to the random encounters with stars (e.g.Quinlan & Hernquist 1997;Berczik et al. 2005;Nasim et al. 2020;Gualandris et al. 2022;Rawlings et al. 2023).This stochasticity in binary eccentricity can affect the estimation of the SMBH merger time-scale and may contaminate the effects of galaxy formation processes.In this study, we systematically mitigate the impact of stochastic eccentricities by fixing them to similar values in the different simulations.
The results in this work reveal a strong correlation between the merger time-scales and the nuclear star formation during the galaxy merger.This suggests another crucial process for accelerating the hardening of SMBH binaries, which is particularly important in gas-rich disc mergers.
This paper is organized as follows.In Section 2, we describe the numerical code and the simulation details.The evolution of orbital parameters and the comparison of SMBH merger time-scales are presented in Section 3. The physical processes driving SMBH coalescence are analysed and discussed in Section 4. We summarize and conclude in Section 5. Appendices A and B offer details on the eccentricity resetting and the resolution study, respectively.

NUMERICAL SIMULATIONS
As we focus on the evolution of SMBH binaries during the hardening and GW emission phases in this study, we initiate the simulations by selecting a snapshot from a parent galaxy merger simulation at the point where the SMBH binary becomes hard (t = t hard ).Subsequently, the binary eccentricity is reset to a desired value, and the simulation is evolved forward until the two SMBHs eventually merge.In this section, we first provide an overview of our numerical code, then detail the information of the parent galaxy merger runs, and finally introduce the simulations that were rerun from t hard .
Following Merritt (2013), in this work, the hard binary separation is defined as where G is the gravitational constant, M red = MBH,1MBH,2/(MBH,1 + MBH,2) is the reduced mass of the SMBH binary, MBH,1 and MBH,2 are the masses of the two SMBHs (assuming MBH,2 ⩽ MBH,1), and σ⋆ is the onedimensional stellar velocity dispersion within the half stellar mass radius.The hard binary separation is roughly the radius where the specific binding energy of the SMBH binary equals the specific kinetic energy of the surrounding stars.An SMBH binary is termed a hard binary when its semimajor axis becomes smaller than R hard .

Numerical code
In this work, the numerical code we used to perform the simulations is the KETJU code (Rantala et al. 2017), which is based on the GADGET-3 code (i.e. the successor of GADGET-2 described by Springel 2005) and utilizes the algorithmically regularized integrator MSTAR (Rantala et al. 2020) to resolve the SMBH dynamics.1To account for the orbital evolution of SMBH binaries due to GW emission, the KETJU code also includes the post-Newtonian (PN) corrections up to the order of 3.5PN for SMBH binaries (Thorne & Hartle 1985;Mora & Will 2004;Blanchet 2014).To compute the gas dynamics, we use the SPHGAL smoothed particle hydrodynamics (SPH) implementation (Hu et al. 2014) with the pressureentropy formulation.We adopt the Wendland C 4 kernel (Dehnen & Aly 2012) and set the number of SPH neighbours to N ngb = 100.
For the galaxy formation subgrid model including radiative cooling, star formation, and stellar feedback, the KETJU code adopts the subgrid model which was originally developed by Scannapieco et al. (2005Scannapieco et al. ( , 2006) ) and later improved in Aumer et al. (2013) and Núñez et al. (2017).This model has been used both in isolated and merger galaxy simulations (Eisenreich et al. 2017;Liao et al. 2023a) and in cosmological zoom-in simulations (Mannerkoski et al. 2021(Mannerkoski et al. , 2022)).The abundances of 11 chemical elements (i.e.H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) are traced in the simulation, and the metal-dependent cooling rate tables are adopted from Wiersma et al. (2009).The star formation conditions are as follows: the gas density ρgas ⩾ 2.2 × 10 −24 g cm −3 (i.e. the hydrogen number density nH ⩾ 1 cm −3 ), the gas temperature Tgas ⩽ 1.2 × 10 4 K, and the gas flow should be convergent (i.e. the velocity divergence ∇ • vgas ⩽ 0).When a gas particle satisfies these criteria, it is stochastically converted into a star particle.For stellar feedback, we consider the effects from both Type Ia and Type II supernova explosions and the winds of asymptotic giant branch (AGB) stars.When stellar feedback is performed, star particles give masses (i.e.metals) and thermal and kinetic energy to their surrounding gas particles.We refer the reader to Núñez et al. (2017) for more details of the stellar feedback model.
For SMBH accretion, we use the SMBH binary accretion subgrid model introduced in Liao et al. (2023a), which extends the widely used Bondi-Hoyle-Lyttleton (BHL) accretion (Hoyle & Lyttleton 1939;Bondi & Hoyle 1944;Bondi 1952) model into the SMBH binary phase and incorporates the preferential accretion subgrid model of CBDs (see Lai & Muñoz 2023, for a recent review).Specifically, when two SMBHs form a gravitationally bound binary, the total accretion rate is determined using the BHL formula with the gas properties computed at the binary's CoM position.Then, the gas accretion on to each SMBH is distributed based on the fitting formula provided by Duffell et al. (2020).As a result, the secondary SMBH (i.e. the one with the lower mass) experiences a higher accretion rate, facilitating the evolution of the binary towards equal masses.This binary accretion model provides more physically motivated SMBH mass evolution in the binary phase, which is important in modelling the AGN feedback and the GW induced recoil velocities.
To study the impact of different AGN feedback implementations, we consider both pure thermal feedback and pure kinetic feedback.For the thermal AGN feedback model, we follow the implementation of Springel et al. (2005).The thermal feedback is continuously performed at each time-step when the SMBH particle is active.Specifically, the amount of feedback energy is first computed according to the SMBH accretion rate, then this energy is added to the internal energy of the surrounding gas particles according to their SPH kernel weights.For the kinetic AGN feedback model, we adopt an implementation similar to the kinetic feedback mode proposed in Weinberger et al. (2017).In contrast to the aforementioned thermal AGN feedback, the kinetic AGN feedback here is implemented using a pulsed approach, i.e. at each time-step when the SMBH is active, the feedback energy -derived from the SMBH accretion rate -is added to a feedback energy reservoir, and the energy in the reservoir is released only after it reaches a predefined threshold.When the release occurs, the surrounding gas particles get kicked and their kinetic energies are boosted.The direction of the kick velocity is equally probable (50 per cent chance) to be either parallel or anti-parallel to either the angular momentum direction of the kicked gas particle with respect to the SMBH (in the single SMBH phase) or the direction of the orbital angular momentum of the SMBH binary (in the binary SMBH phase).More details on the code implementation and the adopted parameter values can be found in Liao et al. (2023a).

Parent runs
The parent galaxy merger simulations are described and analysed in detail in Paper II, consisting of two equal-mass galaxy mergers: the disc-disc merger 'DD-11-G5' and the elliptical-elliptical merger 'EE-11-G5'.These simulation names follow the convention of 'progenitor galaxy types-galaxy mass ratio-orbit geometry' introduced in Liao et al. (2023a).The so-called G5 retrograde orbit, adopted from Naab & Burkert (2003), is specified by the inclination angles between the spin plane and the orbit plane (i) and the arguments of pericentre (ω) for the two galaxies (Toomre & Toomre 1972) orbit geometry is adopted as it exhibits modest starbursts during the merger, which helps to make the simulations computationally efficient.
The initial disc and elliptical galaxies have similar masses in their different components (stars, gas, and dark matter) but differ in their morphology.Specifically, the initial disc galaxy consists of a stellar disc, a stellar bulge, a gas disc, a central SMBH, and a dark matter halo.Conversely, the initial elliptical galaxy is composed of a stellar bulge, a hot gas halo, a central SMBH, and a dark matter halo.For both galaxies, the masses for the dark matter, stellar, and gaseous components, and the SMBH are MDM = 2.5 × 10 12 M ⊙ , M⋆ = 1.2 × 10 11 M ⊙ , Mgas ≈ 2.5 × 10 10 M ⊙ , and MBH = 7.5 × 10 7 M ⊙ , respectively.
The initial particle masses of dark matter, stars, and gas are mDM = 1.6 × 10 6 M ⊙ , m⋆ = 10 5 M ⊙ , and mgas = 10 5 M ⊙ , respectively.The numbers of particles for different galaxy components in the initial conditions are NDM = 3.2 × 10 6 , N⋆ = 2.4 × 10 6 , and Ngas ≈ 5 × 10 5 .The softening lengths adopted Table 1.Simulation information.From left to right, the galaxy merger, the simulation set, the included physical processes, the time when the SMBH binary becomes hard (t hard ), the eccentricity at t hard from the parent simulation, the number of star particles at t hard , the number of gas particles at t hard , and the number of rerun realizations.The abbreviations for the physical processes: Grav -gravity, SPH -smoothed particle hydrodynamics, Cool -gas cooling, Star -star formation and stellar feedback, ThmAGN -SMBH accretion and thermal AGN feedback, KinAGN -SMBH accretion and kinetic AGN feedback.
For both the DD-11-G5 and EE-11-G5 galaxy mergers, we perform simulation sets including different physical processes: (i) NoGas.This is a pure gravity simulation, with all gas particles in the initial condition being converted into star particles.
(ii) NoCool.Compared to the NoGas case, this run includes SPH computations for gas particles with an adiabatic equation of state.
(iii) CoolStarNoAGN.Compared to NoCool, this run further includes gas cooling, star formation, and stellar feedback.However, no SMBH accretion and AGN feedback processes are considered in this run.
(iv) CoolStarThmAGN.Compared to CoolStarNoAGN, this run further includes SMBH accretion and the pure thermal AGN feedback process.
(v) CoolStarKinAGN.In contrast to CoolStarThmAGN, this run incorporates pure kinetic AGN feedback instead.
In total, the parent runs comprises ten simulations.Each simulation starts from the cosmic time of t0 = 10.7 Gyr, and is evolved for a duration of 3 or 4 Gyr (depending on when the two SMBHs finally merge).As we show in Paper II, the properties of our initial galaxies and the galaxy merger remnants from our simulations which include the AGN feedback process (i.e.CoolStarThmAGN and CoolStarKinAGN) are in good agreement with the observations.
For different parent runs, we summarize in Table 1 the simulation time 2 (t hard ), the binary eccentricity (e parent hard ), and the number of star (N⋆) and gas particles (Ngas) from the first snapshot at which the two SMBHs reach the hard binary phase.We refer the interested reader to Paper II for more details and analyses of the parent runs.

2.3
Reruns with e hard = 0.6 The eccentricity of a binary is a key parameter that strongly influences the SMBH merger time-scales.According to Peters (1964), 2 The simulation time is defined as the time span after t 0 .
the rates of decay for the semimajor axis (a) and the eccentricity (e) due to GW emission can be approximated by (i.e. at the 2.5PN level) Here, M bin = MBH,1 + MBH,2 is the total SMBH mass, c is the speed of light.When the eccentricity of a binary is closer to 1, the binary enters the GW emission-dominated phase earlier due to the factor of (1 − e 2 ) −7/2 , resulting in a shorter SMBH merger timescale.In the parent simulations, distinct runs exhibit quite different eccentricities at the time when the two SMBHs become hard (see the column of e parent hard in Table 1).Previous numerical simulations have demonstrated that the eccentricity of a binary at t hard displays inherent stochasticity, arising from random stellar encounters (e.g.Quinlan & Hernquist 1997;Berczik et al. 2005;Nasim et al. 2020;Gualandris et al. 2022;Rawlings et al. 2023).This imposes a limitation on directly comparing SMBH merger time-scales across different parent runs, as well as the study of how factors beyond GW emission influence these time-scales.
To mitigate the impact from the stochastic nature of eccentricity, we manually reset e hard to the same value in different runs, and rerun the simulations starting from t hard and continuing until the point of final SMBH merger.To reset e hard , we take the snapshot of t = t hard from each simulation, and keep the SMBH mass, the binary CoM position and velocity, and the semimajor axis unchanged.The coordinates and velocities of the two SMBHs are adjusted to match the specific relative angular momentum corresponding to the given new eccentricity.The eccentricity of a binary is determined by the magnitude of the specific relative angular momentum, therefore, we can randomly assign different directions to the specific relative angular momentum (or equivalently the orbital plane of the binary) to generate different realizations.The details of how the eccentricity is reset can be found in Appendix A.
In this study, we opt to reset e hard to an intermediate value of 0.6.This choice ensures that the GW emission dominated phase is neither triggered prematurely (making it challenging to disentangle the effects of eccentricity and GW emission from those ascribable to galaxy formation processes) nor delayed excessively (resulting in a longer time for the binary to merge and becoming computationally inefficient), making it well-suited for this pilot investigation.
Once the initial condition for a rerun is set up, we use exactly the same numerical code and parameters as in the original parent run to evolve the simulation until the SMBHs merge.For each parent run, we conduct ten rerun realizations (i.e. they differ in the directions of the binary's specific relative angular momentum) to further minimize the influence of the eccentricity stochasticity on our conclusions.This approach enables us to estimate both the mean and standard deviation of the SMBH merger time-scales, facilitating direct comparisons across different simulations.

Simulation overview
In Fig. 1, we use the EE-11-G5 CoolStarKinAGN run as an illustrative example to provide an overview of our simulations.In row (a), we plot the evolution of the stellar components from the parent run.At around t hard ≈ 2.0 Gyr, the SMBH binary in the centre of the galaxy remnant reaches the hard binary separation, and its eccentricity is then e parent hard = 0.92; see the middle and right panels in row (b).For the rerun, we take the snapshot at t = t hard from the parent run, and then manually reset the binary eccentricity to e hard = 0.6 (the left panel of row b), after which we evolve the simulation until the point of the SMBH merger.
In row (c) of Fig. 1, we compare the evolution of the orbital parameters from the parent run (light blue) and one of the reruns (purple).As we reduce e hard from 0.92 to 0.60, the SMBH merger time-scales increase from ∼1 to ∼2 Gyr.We notice that although we reset the initial eccentricity to a value of 0.6 in the rerun, the eccentricity gradually evolves to a higher value of ∼0.8, before the GW emission starts to dominate and circularizes the binary.This originates from the three-body interactions between the binary and star particles.Similar phenomenon has been observed in previous pure N -body simulations (e.g.Berczik et al. 2005;Merritt et al. 2007), in which a circular binary is initially placed at the centre of an isolated galaxy and the interactions with the passing stars rapidly induce nonzero eccentricities.It has also been shown in the three-body scattering experiments that the eccentricities in the stellar hardening phase tend to grow slightly over time (e.g.Mikkola & Valtonen 1992;Quinlan 1996;Sesana et al. 2006).The rotation of the galaxy remnant can further shape the eccentricity evolution, i.e. the counter-rotation (co-rotation) between the stellar nucleus and the SMBH binary can lead to an increase (decrease) of the binary eccentricity (Sesana et al. 2011).However, we anticipate a mini-mal role of rotation in this run, given that the dimensionless spin parameter (Bullock et al. 2001) of the stellar nucleus at t hard is relatively low, e.g.within 1 kpc (500 pc), it is λ⋆ = 0.0016 (0.0022).3For comparison, the typical spin parameters of random motionsupported dark matter haloes and rotation-supported galactic discs are ∼0.05 and ∼0.4,respectively (Mo et al. 2010).
To further validate this stellar-driven explanation, we conducted a test using the same initial condition as the rerun, but this time incorporating softened gravity between the SMBH and star particles.Note that the interactions between SMBHs remained the same as in the rerun, i.e. non-softened gravity with PN corrections up to the order of 3.5PN.The evolution of orbital parameters in this test run is displayed in row (c) using an olive colour.It is evident that in the absence of strong interactions with stars (i.e.softened gravity), the binary stalls at ∼1 pc.In addition, the binary eccentricity remains a constant around 0.6, providing further support for the stellar-driven eccentricity evolution.

Star formation and central stellar density at t hard
In Fig. 2, we show the stellar distributions of galaxy merger remnants at t hard , which are the initial conditions for the reruns.The spherically averaged stellar density profiles centred on the binary CoM are displayed in row (a), with solid and dotted lines showing the results from all stars and new stars, respectively.Here, the new stars are defined as those formed after the onset of the parent simulations, i.e. stars with a formation time tSF > t0.
It is obvious that the central stellar densities are significantly affected by the formation of new stars, i.e. with more nuclear star formation, the central stellar density is markedly boosted.This is particularly evident in the disc-disc galaxy merger runs, where the absence of AGN feedback in the CoolStarNoAGN run leads to unrealistically high star formation in the centre.Once thermal AGN feedback is included, the central stellar density (r ≲ 0.1 kpc) decreases, but it is still about one order of magnitude higher than in the NoGas and NoCool cases which do not include star formation.The kinetic AGN feedback is more effective in suppressing star formation than the thermal case, resulting in a moderate central stellar density in the CoolStarKinAGN run.
In the elliptical-elliptical galaxy merger, the CoolStarNoAGN run also exhibits a large fraction of nuclear star formation, as there is no AGN feedback to prevent the cooling flows.Once AGN feedback is included, the central densities of the new stars are significantly suppressed.Especially for the case of kinetic AGN feedback, as shown in Paper II, the star formation is completely shut down roughly after the first pericentre passage of the two galaxies.Consequently, all runs except CoolStarNoAGN exhibit very similar stellar density profiles.
The differences in the central stellar distributions are also evident from the projected plots in row (b).Moreover, as shown in Paper II, the central mass distributions at t hard in all runs exhibit clear triaxiality, a natural outcome of galaxy mergers.Compared to spherical potentials, triaxial potentials can assist the binary hardening process by providing more stars with centrophilic orbits, thus avoiding the final parsec problem (Merritt & Poon 2004;Berczik et al. 2006;Khan et al. 2011;Preto et al. 2011).
From the comparisons above, we can conclude that nuclear star formation plays an important role in setting the environment in which the hard SMBH binary evolves.

SMBH merger time-scales from reruns
The time-evolution of the orbital parameters from all the reruns is summarized in Fig. 3.Each simulation set is represented by a distinct colour, and within each set, the results of all ten realizations are plotted.
Although all reruns start with an eccentricity of e hard = 0.6, the interactions with stellar particles result in varying eccentricities among the different realizations, as discussed in Section 3.1.Prior to the GW emission dominated phase, which ultimately circularizes the binary, the eccentricity spread in each simulation set falls approximately within the range of e ∼ 0.6-0.8.This consistency in the eccentricity range allows us to alleviate the impact of eccentricity stochasticity and directly compare the mean SMBH merger time-scales across different simulation sets.In Fig. 3, the mean and standard deviation of the SMBH coalescence time, computed from ten realizations, are indicated by filled circles and error bars for different simulation sets.The detailed SMBH merger time-scales, Tm = t coal − t hard (with t coal being the SMBH coalescence time), are summarized and compared in Table 2 and Fig. 4.
In the case of the disc-disc galaxy merger, the NoGas and No-Cool runs yield similar SMBH merger time-scales of ∼1 Gyr.However, when additional galaxy formation processes are included, the merger time-scales exhibit a clear dependence on the central stellar density, which is affected by the AGN feedback process.In the absence of AGN feedback, the central stellar density is markedly boosted, and the SMBH merger time-scale is notably short, i.e.only ∼20 Myr.When thermal AGN feedback is introduced, the central stellar density is lowered, and this time-scale is extended to ∼140 Myr.With the more effective kinetic AGN feedback, the merger time-scale further increases to ∼450 Myr.
For the elliptical-elliptical galaxy merger, the Cool-StarNoAGN run reveals the shortest merger time-scale of ∼540 Myr, followed by the CoolStarThmAGN run with Tm ∼ 1.4 Gyr.The other runs produce merger time-scales that are similar within the 1σ standard deviation, i.e. ∼2 Gyr.Kinetic AGN feedback proves more effective in maintaining red and quiescent elliptical galaxies (as shown in Paper II).The similarity in SMBH merger time-scales between the CoolStarKinAGN and NoGas runs supports the gas-free assumption, which has been adopted in many previous studies of elliptical galaxies (e.g.Quinlan & Hernquist 1997;Milosavljević & Merritt 2001;Khan et al. 2011;Rantala et al. 2019;Frigo et al. 2021;Nasim et al. 2021).
Across all simulation sets, the SMBH merger time-scales display a consistent pattern, with gas-rich disc galaxy mergers exhibiting shorter time-scales than their gas-poor elliptical counterparts.Specifically, in the simulations with AGN feedback, Tm is ∼100-500 Myr in disc galaxy mergers, while it is ∼1-2 Gyr for elliptical galaxy mergers.It is important to keep in mind that these values are based on the reruns with eccentricities in the range of ∼0.6-0.8 (prior to the GW emission-dominated phase).The quantitative SMBH merger time-scales will vary if we reset e hard to other values in the rerun experiments.
By comparing the SMBH merger time-scales in Fig. 4 to the central stellar densities illustrated in Fig. 2, we notice a pronounced correlation between the SMBH merger time-scales and the nuclear stellar densities, i.e. higher nuclear stellar densities correspond to shorter SMBH merger time-scales.This is present in all simulation sets (including both DD-11-G5 and EE-11-G5 runs).As the variation of the nuclear stellar density here is primarily driven by star formation, this correlation hints that nuclear star formation plays an important role in influencing the binary hardening, especially in gas-rich disc galaxy mergers.We will scrutinize this picture quantitatively in the following section.
Finally, we note that we have performed a resolution study based on the CoolStarKinAGN run.This involved both increasing the mass resolution by a factor of two and decreasing the mass resolutions by factors of two and four (see Appendix B for details).We found that for the DD-11-G5 runs, there exists a mild resolution dependence for the SMBH merger time-scales, where higher mass resolutions correspond to slightly shorter Tm.However, the merger time-scale from the runs presented above agrees well with that from the higher resolution runs.For the EE-11-G5 runs, no resolution dependence is observed for the merger time-scales.

PHYSICAL PROCESSES DRIVING SMBH HARDENING
The relations between the dynamical quantities (i.e. the specific energy ε and specific angular momentum l) and the orbital parameters (i.e. a and e) of a binary system are given by Taking the logarithm and then the time derivative for both sides, we get Ṁbin where the over-dot symbol (˙) denotes the time derivative.If we consider the approximation of ė ≈ 0, then This suggests that the changes in both the total mass of the binary and its specific angular momentum can trigger changes in its orbital evolution.Specifically, SMBH mass accretion leads to a shrinking of the orbit; the decrease in the specific angular momentum due to three-body interactions contributes to the shrinking of the binary orbit; and the change in specific angular momentum due to accretion or torque interactions with surrounding gas (e.g.CBD) will also cause orbital evolution.Furthermore, gas accretion on to SMBHs can introduce additional Brownian motion to the binary and increase the loss cone angular momentum; AGN feedback can alter directly the distribution of surrounding gas and subsequently other galaxy components, thereby impacting the triaxial potential.
All of these effects can potentially influence the loss cone refilling rate, and thus affect the binary hardening rate.At the final stage of the binary evolution, the loss of energy and angular momentum due to GW emission dominates and the two SMBHs coalesce rapidly.Hence, the inclusion of gas and the related galaxy formation processes greatly increase the complexity of the SMBH binary hardening process.In the following subsections, we first analyse the effects of stellar interactions, then study the combined impact from the various galaxy formation processes, and finally assess the influence of the unresolved binary-CBD torque interaction.

Effect of stellar interactions
In Fig. 5, we plot the time-evolution of the inverse semimajor axis, 1/a, and the hardening rate, Different simulation sets are distinguished by different colours, and all ten rerun realizations are plotted for each set.Note that the hardening rate is estimated by performing a linear fit to 1/a(t) in every interval of ∆t = Tm,max/N bin , where Tm,max is the maximum merger time-scale across ten realizations and N bin = 100.We have performed tests with other values for N bin (e.g.50 and 200) and confirmed that the presented results are not sensitive to the adopted bin number.
Prior to the dominance of GW emission, which results in a rapid orbital decay, the inverse semimajor axis tends to be a linear function of time (or equivalently the hardening rate tends to be constant).This is particularly noticeable in simulation sets where the three-body interaction is the only mechanism driving SMBH hardening, such as the NoGas and NoCool runs.We also notice that the simulation sets with shorter SMBH merger time-scales generally exhibit higher hardening rates at t hard .For instance, in the DD-11-G5 CoolStarNoAGN run, which possesses the shortest Tm, the hardening rate is s ∼ 3 pc −1 Myr −1 at t hard .On the other hand, in the EE-11-G5 NoCool run, which has the longest Tm, the hardening rate is much lower, i.e. s ∼ 0.01 pc −1 Myr −1 .
From three-body scattering experiments, it has been shown that the hardening rate is a simple function of the stellar properties (e.g.Quinlan 1996;Sesana et al. 2006), where the dimensionless hardening parameter H is approximately a constant once the binary becomes hard, and ρ⋆ and σ⋆ are the density and velocity dispersion of the stellar background, respectively.
In Fig. 6, we compare the relation between s and Gρ⋆/σ⋆ at t hard from our reruns to the fitting relation from three-body scattering experiments (Sesana et al. 2006).For each simulation set, the mean and standard deviation of the hardening rate is computed using all ten realizations; the mean density, ρ⋆, and the onedimensional velocity dispersion, σ⋆, are computed using the star particles within the influence radius, rm.Here, the influence radius is defined as the distance from the binary CoM, within which the enclosed stellar mass is twice the binary mass, i.e.M⋆(< rm) = 2M bin .
The relatively good agreement between our simulations and the three-body scattering experiments suggests that the SMBH binary-star interactions are the primary mechanism in driving the binary hardening at t hard .Fig. 6 also clearly shows that compared to the case of elliptical galaxy mergers, the more rapid SMBH hardening in disc galaxy mergers is due to the higher central stellar density (or more precisely the higher ρ⋆/σ⋆).
In Fig. 7, we further study the time-evolution of the loss cones during the interval from t hard until the semimajor axis shrinks by a factor of 10 (before entering the GW-driven phase).To achieve better time resolution, we have randomly chosen one realization from each simulation set and increased its snapshot output frequency if necessary. 4It is worth noting that other realizations with lower snapshot cadence show qualitatively similar results.The top panels show the fraction of stellar mass enclosed within the loss cone, denoted as M⋆,LC/M⋆, where M⋆ is the total stellar mass in the simulation and M⋆,LC is the total mass of the stars in the loss cone.Here, the loss cone angular momentum is estimated as LLC = √ 2GM bin a, and a star particle with specific angular momentum 5 L⋆ < LLC is classified as being part of the loss cone Myr.In the reruns of the CoolStarNoAGN, CoolStarThmAGN, and Cool-StarKinAGN simulation sets of the DD-11-G5 merger with high snapshot cadence, the output intervals are reduced to 0.2, 2, and 5 Myr, respectively.In the EE-11-G5 CoolStarNoAGN rerun, the output interval is reduced to 10 Myr.The remaining simulations sets keep the original output interval. 5The specific angular momentum of a star particle is computed with respect to the CoM of the SMBH binary.We have also tested with the CoM of all star particles, and the results are similar.population.Overall, the loss cone stellar mass fractions decrease as the binary shrinks and more star particles are ejected from the vicinity of the binary.In disc galaxy mergers, one would expect that the runs including gas cooling and star formation have in general an enhanced fraction of stars in the loss cone.Notably, the CoolStarK-inAGN run exhibits the highest mass fractions, primarily attributing to its high LLC -an outcome of a larger binary mass (arising from the binary becoming hard later and having a longer time to accrete more gas; see Paper II for the detailed SMBH growth history) and a larger hard binary separation (see Fig. 3).In contrast, the mass fractions are more similar across the different simulation sets in the elliptical galaxy mergers due to the smaller differences in star formation and SMBH mass growth.
The bottom panels of Fig. 7 display the mass fractions of new stars (i.e.stars with tSF > t0) in the loss cone.In disc galaxy mergers, ∼80% of the stars in the loss cone of the CoolStarNoAGN run  (Sesana et al. 2006).Note that Sesana et al. (2006) only provide fitting formulae of H for circular binaries with different mass ratios, and thus here H = 11.4 is computed assuming an eccentricity e = 0 and a mass ratio of q = 1.According to fig. 3 of Sesana et al. (2006), at t hard , the hardening rate parameters for different eccentricities are fairly close.The agreement between our simulation results and the three-body scattering experiments, which only considers stellar scattering, suggests that the SMBH binarystar interactions are the primary mechanism in driving binary hardening at t = t hard in our simulations.are newly formed, followed by the CoolStarThmAGN and Cool-StarKinAGN runs with fractions of ∼50% and ∼15%, respectively.For the elliptical galaxy mergers, the CoolStarNoAGN run exhibits a mass fraction of new stars fluctuating around ∼20% and the fraction reduces as the binary hardens; the CoolStarThmAGN only shows non-zero fractions in a few snapshots.This comparison highlights the crucial role of new stars in replenishing the loss cone and accelerating the process of binary hardening, particularly for gas-rich disc galaxy mergers.The analyses in this subsection point to the following two step picture in relation to SMBH binary hardening following a galaxy merger: firstly, due to radiative cooling and tidal torques, gas condenses to the galaxy centre and forms stars, boosting the central stellar density.Secondly, these recently formed stars -inheriting the low angular momenta from the gas -contribute to the loss cone and assist in the SMBH hardening process via three-body interactions.This process is particularly important in gas-rich disc galaxy mergers, resulting in shorter SMBH merger time-scales, compared to the gas-poor elliptical galaxy mergers.The comparison across different simulation sets highlights the importance of modelling galaxy formation processes, which affect gas cooling and star formation, in studying the SMBH coalescence process.
It is worth noting that Liao et al. (2023a) already demonstrated a tendency for the stellar particles ejected by the SMBH binary during the hardening phase to have young ages, thereby providing additional support for the picture described above.The simulation results from Khan et al. (2016) also strengthen this picture.By extracting a galaxy merger remnant at z ∼ 3 from a cosmological simulation and further evolving an SMBH binary in the remnant centre with the direct N -body code ϕ-GPU, Khan et al. (2016) found that due to the high central stellar density, which is attributed to recent star formation, the SMBH merger time-scale is as short as ∼10 Myr.Intriguingly, in the absence of AGN feedback in both their parent cosmological simulation and zoomed simulations, their SMBH merger time-scale is quite close to that of our DD-11-G5 CoolStarNoAGN run (i.e.∼10 Myr versus ∼20 Myr).However, our results suggest that once AGN feedback is included, the SMBH merger time-scales become much longer at ∼100-500 Myr.

Combined effects of galaxy formation processes
During the hardening phase, as the SMBH binary is dynamically ejecting stars via the gravitational slingshot interaction, the various galaxy formation processes also come into play.Specifically, SMBHs accrete gas and release feedback energy, which increases the SMBH masses, introduces Brownian motions, and alters the surrounding matter distribution, directly affecting the binary hardening.Simultaneously, gas cools and turns into stars stochastically (if the star formation conditions are met) and stars influence surrounding gas via stellar feedback processes.These combined effects can further impact the distribution of matter and the repopulation of the loss cone.Given that all of these processes occur simultaneously and interplay with each other, it is challenging to distinctly isolate each effect.
To study the combined impact of these galaxy formation processes on the SMBH coalescence, we perform the 'NoCool test' reruns -starting from the same initial conditions as the Cool-StarNoAGN, CoolStarThmAGN, and CoolStarKinAGN reruns but running them with the NoCool model, which only considers gravity and SPH with an adiabatic equation of state.These NoCool test reruns (ten realizations for each set) are plotted with a brown colour in Fig. 5.In the corresponding 'Info' panels, we summarize the SMBH merger time-scales for both the original reruns (Tm) and the NoCool test reruns (T ′ m ), including both the mean and the standard deviation.In addition, we calculate the ratio of their mean merger time-scales (T ′ m /Tm).Moreover, we provide the averaged star formation rate (SFR) over Tm and the averaged SMBH mass growth factor (M bin,f /M bin,i with M bin,i being the binary mass at t hard and M bin,f being the binary mass at t hard + Tm) for the original reruns.
By comparing the NoCool test reruns with the original reruns, we can obtain an estimation of the impact from the galaxy formation processes.6Interestingly, all runs except DD-11-G5 CoolStarNoAGN and EE-11-G5 CoolStarKinAGN have an SMBH merger time-scale ratio of T ′ m /Tm ∼ 1.7, indicating that, on average, the inclusion of galaxy formation processes has shortened the SMBH merger time-scales by a factor of ∼1.7.In the case of DD-11-G5 CoolStarNoAGN, the high hardening rate due to the existing stellar background overshadows the impact of gas cooling and star formation on the SMBH merger time-scale.In the case of EE-11-G5 CoolStarKinAGN, the low star formation rate (∼0.01 M ⊙ yr −1 ) and minimal SMBH mass growth (∼1%) due to the effective kinetic AGN feedback result in a smaller reduction in the merger time-scale (only by a factor of ∼1.1).

Effects from the unresolved binary-CBD torque interactions
In our subgrid model for binary SMBH accretion, we incorporate the preferential accretion mechanism motivated by the small-scale hydrodynamical CBD simulations, i.e. the secondary SMBH tends to have a higher accretion rate than the primary SMBH.However, our SMBH subgrid model does not include the binary-CBD torque interaction, which can also affect the orbital evolution of the binary.This omission stems partly from the ongoing debates within this field and the large uncertainty in the existing results.To elaborate, early CBD literature suggested a negative gravitational torque from the CBD acting upon the binary, which caused the binary's orbit to shrink (e.g.Artymowicz et al. 1991;MacFadyen & Milosavljević 2008).On the contrary, some recent CBD simulations have proposed a positive torque: while the torque from the outer CBD is negative, as previously suggested, the torque from the gas within the inner streams and mini-discs can be overwhelmingly positive.This results in a net positive torque and, consequently, the expansion of the binary (e.g.Miranda et al. 2017;Tang et al. 2017;Moody et al. 2019;Muñoz et al. 2019Muñoz et al. , 2020;;Duffell et al. 2020).Subsequent simulation studies have demonstrated that whether an SMBH binary shrinks or expands can depend on various factors, in-cluding the specific disc parameters such as the disc thickness and viscosity (e.g.Heath & Nixon 2020;Tiede et al. 2020;Dittmann & Ryan 2022, 2023), the binary properties like the SMBH mass ratio and binary eccentricity (e.g.D 'Orazio & Duffell 2021;Siwek et al. 2023), and the detailed simulation setups (e.g.2D versus 3D, inclusion of CBD self-gravity or not, and fixed versus live SMBH binaries; see e.g.Franchini et al. 2021Franchini et al. , 2022Franchini et al. , 2023)).
In the following, we estimate the impact of the binary-CBD torque interaction on our binary hardening results.To achieve this, we compare the binary inspiral/outspiral rate, ȧ/a, measured from CBD simulations with the inspiral rate directly measured from our simulations.
From CBD simulations, it has been found that under the binary-CBD interaction, an eccentric and equal-mass binary with e ≳ 0.1 evolves towards an equilibrium eccentricity of eeq ∼ 0.5 while a nearly circular and equal-mass binary (e ≲ 0.1) stays in a circular orbit with eeq ∼ 0 (e.g.D 'Orazio & Duffell 2021;Zrake et al. 2021;Siwek et al. 2023).This implies that our simulated binary, which has MBH,2/MBH,1 ≈ 1 and e hard = 0.6, will settle at e ≈ 0.5 within a few Myr (Siwek et al. 2023) when the binary-CBD interaction is considered.
In the CBD study, the inspiral/outspiral rate of the binary orbit is usually expressed using the specific accretion rate of the binary system.For an equal-mass binary with e = 0.5, Muñoz et al. (2019) Note that in both relations, ȧ/a is positive, implying that the binary expands due to the torque interaction with its CBD.It is important to compare the outspiral rate from CBD simulations, which consider only the binary-CBD interaction, with the inspiral rate from our simulations, which consider the binary-star interaction and the effects from other galaxy formation processes, to see which mechanism dominates and whether our simulated binary shrinks or expands.
This comparison is summarized in Fig. 8 for the reruns with AGN feedback.Note that different rerun realizations of the same simulation set exhibit qualitatively similar results, therefore only one representative realization from each set is displayed here.The upper three rows present the time-evolution of the binary accretion rate, the binary mass, and the Eddington ratio of the binary system.All runs show binary accretion rates below the Eddington limit, with f Edd ≲ 0.3.The fourth row displays the magnitude of the inspiral rate (i.e.− ȧ/a) from our simulations (in red) and the estimated CBD-induced outspiral rates according to Eqs (11) and (12) (in olive and cyan respectively).It is evident that in all runs, the inspiral rates attributed to the binary-star interaction and other galaxy formation processes are higher than the outspiral rates resulting from the binary-CBD torque interaction.According to the bottom panel, the former is at least a few times higher than the latter.For runs with low accretion rates (e.g.EE-11-G5 CoolStarK-inAGN), the ratio is notably within the range of ∼10 2 -10 4 .This quantitative comparison suggests that for our galaxy mergers, prior to the GW emission-dominated phase, the primary driving mechanisms for binary orbital evolution are the binarystar slingshot interaction and other galaxy formation processes.The binary-CBD torque interaction appears to play only a complementary and secondary role.In gas-rich mergers, despite the high binary accretion and CBD-induced outspiral rates, efficient star formation boosts the central stellar density and the inspiral rate.Conversely, gas-poor mergers exhibit low binary accretion rates, resulting in even lower CBD-induced outspiral rates despite the low inspiral rates from binary-star interactions.Our results overall agree with the work of Bortolas et al. (2021), which suggests that CBDs have a minimal impact on binary evolution when accretion rates are significantly lower than the Eddington limit.

CONCLUSIONS
In this work, we utilized the simulations of the RABBITS series to study how various galaxy formation processes (e.g.cooling, star formation, SMBH accretion, stellar and AGN feedback) influence the orbital evolution of SMBH binaries during the hardening and GW emission phases.Especially, we systematically mitigated the impact of the stochastic binary eccentricity at t hard by adjusting it to the same value across the different simulations.
Our major findings are summarized as follows.
(i) Gas-rich disc galaxy mergers exhibit significantly shorter SMBH merger time-scales than gas-poor elliptical mergers.For eccentricities of e ∼ 0.6-0.8during the hardening phase, the simulations with AGN feedback reveal merger time-scales of Tm ∼ 100-500 Myr for disc mergers and Tm ∼ 1-2 Gyr for elliptical mergers.Elliptical merger runs with kinetic AGN feedback yield similar time-scales to gas-free runs, supporting the gas-free assumption for simulating elliptical galaxy mergers adopted in previous studies.
(ii) Nuclear (or centrally concentrated) star formation has a strong influence on SMBH merger time-scales.Throughout the galaxy merging process, the gas condensation at the galaxy centre due to radiative cooling and tidal torques leads to star formation in the nuclear region of the galaxy, boosting the central stellar density (i.e.within the radius of a few hundred parsecs).These recently formed stars, inheriting low angular momenta from the gas, contribute to the loss cone, assisting binary hardening, especially in gas-rich disc mergers.This also explains why the simulations with cooling and star formation but without AGN feedback (i.e. the CoolStarNoAGN runs) exhibit the shortest merger time-scales.
(iii) After an SMBH binary becomes hard, compared to nonradiative hydrodynamical runs (i.e. the NoCool test runs), galaxy formation processes collectively shorten the SMBH merger timescale by a factor of ∼1.7 on average.
(iv) In our simulated retrograde galaxy mergers, the orbital evolution of SMBH binaries is primarily driven by binary-star gravitational slingshot interactions and other galaxy formation processes, with the binary-CBD torque interaction only playing a minor role.The outspiral evolution due to the binary-CBD torque interaction has minimal impact on the shrinking of the orbit and SMBH merger time-scales in our studied mergers.
Our study highlights the crucial role of nuclear star formation in driving the coalescence of SMBHs, especially in gas-rich disc galaxy mergers.It also reinforces the importance of including galaxy formation processes in improving the predictions of SMBH merger time-scales.
In the current work, we have considered only the G5 galaxy merger geometry and fixed e hard to a single value of 0.6.To further explore the parameter space and investigate the relative roles of the binary-star and binary-CBD interactions, in our future RABBITS studies, we will consider other merger geometries (see e.g.Naab & Burkert 2003), which can result in different strengths of tidal response and nuclear starburst, and adjust e hard to other values.These upcoming investigations are expected to further deepen our understanding in the hardening and coalescing process of SMBH binaries, a key target of low-frequency GW observations.Table B1.Details of the CoolStarKinAGN simulations with different numerical resolutions.From the left, the galaxy merger, the simulation name, the number of star particles at t hard (N⋆), the number of gas particles at t hard (Ngas), the number of dark matter particles (N DM ), the averaged star/gas particle mass (m⋆ or mgas), the dark matter particle mass (m DM ), the stellar softening length (ϵ⋆), the gas softening length (ϵgas), the dark matter softening length (ϵ DM ), and the mean and standard deviation of the SMBH merger time-scale (Tm).4, but for the CoolStarKinAGN resolution test runs.For higher resolution runs, the standard deviation of the SMBH merger time-scales tends to be smaller.In the DD-11-G5 runs, the merger timescales show a weak resolution dependence (i.e. with shorter Tm for runs with higher resolutions), but the merger time-scale from the Fiducial-Np runs agrees well with that from the Hi-res-2Np runs within the one-sigma scatter.In the EE-11-G5 simulations, there is no clear resolution dependence for the merger time-scales.The hardening rates here are estimated by performing linear fits to 1/a(t) in the interval ∆t from t hard to the time when a(t) = a(t hard )/2.The DD-11-G5 and EE-11-G5 results are shown with blue and orange colours, respectively.The mean and standard deviation computed from the ten realizations of each resolution set are plotted with the filled circles/diamonds and error bars.For the DD-11-G5 runs, the hardening rates show a dependence on particle number in the low-resolution runs, however, the Fiducial-Np run converges well to the Hi-res-2Np run.

Galaxy
For the EE-11-G5 runs, there is no dependence of hardening rates on particle number.
the interval from t hard to the point when a(t) = a(t hard )/2.In the DD-11-G5 reruns, the hardening rates show a dependence on particle number in the low-resolution runs (i.e. with higher hardening rates for reruns with more particles), yet the Fiducial-Np rerun converges to the Hi-res-2Np well.Note that the particle number dependence observed in the low-resolution runs here is not related to the N -dependence observed in spherical galaxies from pure Nbody simulations (e.g.Makino & Funato 2004;Berczik et al. 2005).
In the latter case, the N -dependence arises from the refilling of the loss cone due to two-body relaxation, exhibiting an opposite Ndependence (i.e.lower hardening rates in simulations with more particles).Here, the galaxy remnant is triaxial, and the resolution dependence in the low-resolution runs originates from the galaxy formation subgrid processes.Specifically, as shown in appendix A of Paper II, the two low-resolution DD-11-G5 parent runs exhibit lower star formation rates, leading to lower central stellar densities and consequently lower hardening rates.
In the EE-11-G5 runs, there is no dependence of the hardening rates on particle number.Note that from pure N -body simulations, the lack of particle number influence on binary hardening rates is expected for triaxial remnants of galaxy mergers (Merritt & Poon 2004;Berczik et al. 2006;Khan et al. 2011;Preto et al. 2011).
To conclude, the results (such as the binary orbital evolution and SMBH merger time-scales) from our fiducial resolution reruns, which are presented in the main text, converge well with those from the reruns conducted with a higher mass resolution.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Overview of the simulations used in this study.(a) Evolution of the stellar components in the CoolStarKinAGN parent run of the EE-11-G5 galaxy merger.(b)The right panel displays the projected stellar density in the 2 × 2 kpc region centered on the SMBH binary at t = t hard .In the middle panel, the two SMBHs are represented by black filled circles, and their first few orbits after the time t hard from the parent run are shown.The red cross indicates the centre of mass of the binary.In this panel, the SMBH orbits are depicted on the orbital plane.The left panel illustrates the SMBHs and their first few orbits from the rerun, where the binary eccentricity was reset to 0.6.(c) Evolution of the PN-corrected orbital parameters in the parent run (light blue), the rerun (purple), and a test rerun with softened gravity between the SMBH and star particles (olive).The comparison between the rerun and the test rerun clearly shows that in the absence of strong interactions with stars (i.e.softened gravity), the binary stalls at ∼1 pc and its eccentricity remains constant.

Figure 2 .
Figure 2. Stellar distributions of galaxy remnants at t = t hard from simulations including different galaxy formation processes.(a) Spherically averaged stellar density profiles at t = t hard from the DD-11-G5 (left) and EE-11-G5 (right) runs.These profiles are centred on the CoM of the binary.The NoGas, NoCool, CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are displayed with blue, orange, green, red, and purple lines, respectively.The solid lines plot the total stellar density profiles, while the dotted lines show the density profiles of newly formed stars.Here, the new stars are defined as the stars formed after the start of the parent simulations, i.e. with star formation time t SF > t 0 .(b) The 2 × 2 kpc projected stellar distributions (on the xy-plane) centred on the SMBH binary from the DD-11-G5 (upper) and EE-11-G5 (lower) runs.From left to right, the NoGas, NoCool, CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are plotted.

Figure 3 .
Figure3.Time-evolution of the orbital parameters in simulations that start from snapshots when the two SMBHs become dynamically hard (with the initial orbital eccentricity reset to 0.6).The NoGas, NoCool, CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are shown with blue, orange, green, red, and purple colours, respectively.For each simulation set, all ten realizations are shown here.In the eccentricity panels, the filled circles and error bars mark the mean and the standard deviation of the SMBH coalescence time from ten realizations, and they are vertically shifted for different simulation sets in order to provide a clear illustration.The horizontal dashed lines mark e = 0.6.

Figure 4 .
Figure 4. SMBH merger time-scales from different runs.The blue circles and the orange diamonds show the mean merger time-scales (averaged over ten realizations) from DD-11-G5 and EE-11-G5 runs, respectively.For each set of runs, the error bar shows the standard deviation of the merger timescales.

TFigure 5 .
Figure 5.Time evolution of the inverse semimajor axis (1/a) and hardening rate (s) from the DD-11-G5 (upper) and EE-11-G5 runs (lower).From left to right, all ten realizations of the NoGas (blue), NoCool (orange), CoolStarNoAGN (green), CoolStarThmAGN (red), and CoolStarKinAGN (purple) simulation sets are plotted.For each simulation set, the top panel summarizes the following information: the mean and standard deviation of the merger time-scales Tm, the averaged SFR over Tm (if star formation is included), and the averaged SMBH mass growth factor M bin,f /M bin,i (if SMBH accretion is included); the middle and bottom panels plot 1/a and s, respectively.In the CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN panels, the brown lines show the results from the NoCool tests, which start from the corresponding initial conditions but evolve with gravity and SPH only.The merger time-scale, T ′ m , of the NoCool test run and the ratio of the mean merger time-scales, T ′ m /Tm, are summarized in the top panel.

Figure 6 .
Figure6.Relation between the hardening rate and the stellar properties within the influence radius rm at t = t hard .The DD-11-G5 and EE-11-G5 runs are plotted with filled circles and diamonds, respectively.Simulations including different physical processes are distinguished by colours following the same convention used in the other figures.The dashed line shows the hardening rate relation obtained from three-body scattering experiments(Sesana et al. 2006).Note thatSesana et al. (2006) only provide fitting formulae of H for circular binaries with different mass ratios, and thus here H = 11.4 is computed assuming an eccentricity e = 0 and a mass ratio of q = 1.According to fig.3ofSesana et al. (2006), at t hard , the hardening rate parameters for different eccentricities are fairly close.The agreement between our simulation results and the three-body scattering experiments, which only considers stellar scattering, suggests that the SMBH binarystar interactions are the primary mechanism in driving binary hardening at t = t hard in our simulations.

Figure 7 .
Figure7.Evolution of loss cones for the DD-11-G5 (left) and EE-11-G5 (right) runs.The top panels display the stellar mass fractions within the loss cone, while the bottom panels plot the mass fractions of new stars (i.e.stars formed after t 0 ) in the loss cone.Different simulations are distinguished by colours following the convention used in previous figures.The x-axis is the ratio between the semimajor axis at t hard and the semimajor axis at a given time t.Only the time period from t hard until the semimajor axis shrinks by a factor of 10 (before the domination of GW emission in driving the binary hardening) is plotted here.
provided the following fitting relation from their simulations

Figure 8 .
Figure8.Comparison between our simulations and the expected effects from the CBD torques.The two left columns show the results from the CoolStarTh-mAGN (red) and CoolStarKinAGN (purple) runs of the DD-11-G5 merger, while the two right columns plot those runs of the EE-11-G5 merger.From top to bottom, we plot the total SMBH binary accretion rate ( Ṁbin ), the SMBH binary mass (M bin ), the Eddington ratio (f Edd ≡ Ṁbin / ṀEdd ), the magnitude of the SMBH binary inspiral/outspiral rate (| ȧ/a|), and the ratio between the inspiral rate from our simulations and the outspiral rate expected from the CBD torques.In the bottom two rows, the olive and cyan curves show the results of the outspiral SMBH binary from CBD simulations ofMuñoz et al. (2019) andSiwek et al. (2023), respectively.Note that the plotted physical quantities have been averaged over 1 Myr, 5 Myr, 10 Myr, and 10 Myr for the runs from left to right.

Figure B1 .Figure B2 .
Figure B1.Similar as Fig. 3, but for the CoolStarKinAGN resolution study runs.The Low-res-Np/4, Low-res-Np/2, Fiducial-Np, and Hi-res-2Np sets are plotted with the blue, orange, green, and red colours, respectively.For each resolution set, the filled circle and the error bar show the mean SMBH coalescence time over ten realizations and their standard deviation.The DD-11-G5 test runs show a weak resolution dependence for the SMBH merger time-scales, while the merger time-scales in the EE-11-G5 runs do not exhibit resolution dependences; see Fig. B2 for details.

Figure B3 .
Figure B3.Binary hardening rates from the CoolStarKinAGN reruns with different particle numbers.The hardening rates here are estimated by performing linear fits to 1/a(t) in the interval ∆t from t hard to the time when a(t) = a(t hard )/2.The DD-11-G5 and EE-11-G5 results are shown with blue and orange colours, respectively.The mean and standard deviation computed from the ten realizations of each resolution set are plotted with the filled circles/diamonds and error bars.For the DD-11-G5 runs, the hardening rates show a dependence on particle number in the low-resolution runs, however, the Fiducial-Np run converges well to the Hi-res-2Np run.For the EE-11-G5 runs, there is no dependence of hardening rates on particle number.

Table 2 .
Merger time-scales (mean ± standard deviation) from the reruns.The values are computed from ten realizations.