RABBITS – II. The impact of AGN feedback on coalescing supermassive black holes in disc and elliptical galaxy mergers

In this study of the ‘Resolving supermAssive Black hole Binaries In galacTic hydrodynamical Simulations’ (RABBITS) series, we investigate the orbital evolution of supermassive black holes (SMBHs) during galaxy mergers. We simulate both disc and elliptical galaxy mergers using the KETJU code, which can simultaneously follow galaxy (hydro-)dynamics and small-scale SMBH dynamics with post-Newtonian corrections. With our SMBH binary sub-grid model, we show how active galactic nuclei (AGNs) feedback affects galaxy properties and SMBH coalescence. We find that simulations without AGN feedback exhibit excessive star formation, resulting in merger remnants that deviate from observed properties. Kinetic AGN feedback proves more effective than thermal AGN feedback in expelling gas from the centre and quenching star formation. The different central galaxy properties, which are a result of distinct AGN feedback models, lead to varying rates of SMBH orbital decay. In the dynamical friction phase, galaxies with higher star formation and higher SMBH masses possess denser centres, become more resistant to tidal stripping, experience greater dynamical friction, and consequently form SMBH binaries earlier. As AGN feedback reduces gas densities in the centres, dynamical friction by stars dominates over gas. In the SMBH hardening phase, compared to elliptical mergers, disc mergers exhibit higher central densities of newly formed stars, resulting in accelerated SMBH hardening and shorter merger time-scales (i.e. ≲ 500 Myr versus ≳ 1 Gyr). Our findings highlight the importance of AGN feedback and its numerical implementation in understanding the SMBH coalescing process, a key focus for low-frequency gravitational wave observatories.

Galaxy mergers naturally occur within the framework of hierarchical structure assembly in the cold dark matter scenario (see e.g.Mo, van den Bosch & White 2010;Frenk & White 2012, for reviews).When two massive galaxies merge, it is predicted that their central SMBHs undergo a three-phase evolutionary process (Begelman, Blandford & Rees 1980): (i) at ∼kpc scales, the two SMBHs sink to the galaxy remnant centre due to dynamical friction from stars and gas (Chandrasekhar 1943;Ostriker 1999), eventually forming a gravitationally bound binary.(ii) At ∼pc scales, the SMBH binary continues to harden (i.e. its orbit shrinks) by losing energy and angular momentum to stars through gravitational slingshot interactions (e.g.Mikkola & Valtonen 1992;Quinlan 1996;Sesana, Haardt & Madau 2006).The SMBH binary will also interact with the circumbinary gas disc (see Lai & Muñoz 2023, for a review).(iii) At ∼mpc scales, the binary evolution is dominated by gravitational wave (GW) emission (Peters & Mathews 1963;Peters 1964) until coalescence.The GW emitted from these coalescing SMBH binaries are key targets for the ongoing and upcoming lowfrequency GW observatories, such as pulsar timing array (Burke-Spolaor et al. 2019;Agazie et al. 2023) operating in the nano-hertz frequency range, and the future space GW observatories working in the milli-hertz frequency range, including the Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017, 2023), TianQin (Luo et al. 2016), andTaiji (Ruan et al. 2020).
AGN feedback plays an important role in shaping the general properties of merging galaxies.Early simulations using the thermal AGN feedback model (Springel et al. 2005) revealed that in gasrich disc galaxy mergers, the intense nuclear gas inflows induced by tidal interactions trigger starbursts and provide fuel to the central SMBHs.Subsequently, the high SMBH accretion rates result in powerful AGN feedback which in turn expels the gas, reduces nuclear star formation, and quenches the galaxy remnant (e.g.Di Matteo et al. 2005;Springel et al. 2005;Hopkins et al. 2006;Johansson et al. 2009).Later simulations demonstrated that different implementations of AGN feedback can significantly impact specific galaxy properties (e.g.Choi et al. 2012Choi et al. , 2014Choi et al. , 2015;;Debuhr et al. 2012;Wurster & Thacker 2013;Barai et al. 2014;Eisenreich et al. 2017).For instance, compared to thermal AGN feedback, the kinetic (or mechanical) AGN feedback leads to highervelocity gas outflows, reduced thermal gas X-ray luminosity, and greater AGN variability (e.g.Choi et al. 2012Choi et al. , 2014;;Debuhr et al. 2012; Barai et al. 2014).Recently, an AGN-driven wind injection model has been introduced to further improve the implementation of AGN feedback in galaxy formation simulations (see e.g.Costa et al. 2020;Torrey et al. 2020;Bollati et al. 2023a).
However, one limitation of the aforementioned hydrodynamical galaxy merger simulations is the adoption of a simplified ad hoc merger criterion.In these simulations, two SMBHs are considered to merge instantaneously when their distance falls below a certain predefined threshold (e.g. the SMBH smoothing length and/or gravitational softening length), typically at the kpc or sub-kpc scales. 1his approach means that the three-phase coalescence process of SMBHs is not fully resolved in these simulations.Consequently, the GW-related predictions (e.g.SMBH merger rates) from cosmological galaxy formation simulations, which often employ similar SMBH merger criteria as galaxy merger simulations, rely on either an ad hoc constant delay time (e.g.Salcido et al. 2016;DeGraf et al. 2021) or post-processing semi-analytical models (e.g.Kelley et al. 2017;Volonteri et al. 2020;Chen et al. 2022;Li et al. 2022).
Efforts have been made to improve the modelling of SMBH pair (i.e. two close SMBHs that have not yet formed a gravitationally bound binary) or binary dynamics in galactic-scale sim-ulations.These improvements can be categorized into two groups: elliptical and disc galaxy mergers.
Elliptical galaxy mergers.For gas-poor elliptical galaxies, the approximation of gas free is typically adopted and pure Nbody simulations are used to study the SMBH coalescence process in such systems (e.g.Quinlan & Hernquist 1997;Milosavljević & Merritt 2001;Berczik et al. 2006;Khan et al. 2011;Rantala et al. 2019;Frigo et al. 2021;Nasim et al. 2021).Usually, the dynamical friction and binary hardening phases are resolved in these simulations and some of them also include post-Newtonian (PN) corrections to model the GW emission phase (e.g.Berentzen et al. 2009; Khan et al. 2012a;Rantala et al. 2017;Mannerkoski et al. 2019).From these collisionless simulations, it was found that the SMBH mergers typically take place on a ∼Gyr time-scale unless the SMBH binaries are highly eccentric (e.g.Khan et al. 2012a;Vasiliev et al. 2015;Rantala et al. 2017).
However, all of these elliptical merger simulations do not consider gas and the associated physical processes (e.g.radiative cooling, SMBH accretion, and AGN feedback).In fact, even when we consider previous elliptical merger simulations that did not resolve small-scale SMBH dynamics (e.g.González-García & van Albada 2003;Boylan-Kolchin et al. 2005;Naab et al. 2006), there is still a scarcity of literature on elliptical mergers that include gas (for a few examples, see Johansson et al. 2009;Sinha & Holley-Bockelmann 2009).How does the presence of hot gas haloes in elliptical galaxies and the related physical processes impact the evolution of merging galaxies and coalescing SMBHs?Addressing this question can shed light on the validity of the gas-free approximation commonly adopted in previous elliptical merger simulations.
However, the range of included physical processes varies across these disc galaxy simulations.Especially, only a handful of simulations (e.g.Callegari et al. 2011;Van Wassenhove et al. 2014;Pfister et al. 2017;Souza Lima et al. 2017;Bollati et al. 2023b;Chen et al. 2023) incorporate AGN feedback process.Recently, it was suggested that AGN feedback can weaken the gas dynamical friction on SMBHs (Sijacki et al. 2011;Souza Lima et al. 2017;Bollati et al. 2023b).Moreover, AGN feedback can affect star formation and gas outflows in galaxies, thus altering the galactic structure and potential, which in turn should impact the mo-tion of SMBH binaries.How exactly does AGN feedback, along with its numerical implementations, influence the SMBH orbital evolution in different phases?This question remains relatively unexplored and necessitates further investigation through additional simulation studies.
In addition, the above mentioned disc galaxy hydrodynamical simulations do not resolve the SMBH binary hardening and GW emission phases, i.e. the SMBH binaries or pairs in these simulations cease to shrink as they approach the adopted softening lengths (usually ≳ pc).Some attempts have been made to address the binary hardening and coalescence processes within gasrich galaxy mergers.For example, Khan et al. (2016) selected a disc-disc galaxy merger from the Argo cosmological hydrodynamical simulation and resimulated the very central region of the galaxy remnant and the SMBH binary evolution using the direct N -body code ϕ-GPU with PN corrections.However, it is worth noting that their parent simulations did not include AGN feedback.Additionally, prior to the direct N -body calculations, the authors manually converted gas particles into stars (see similar treatments in Khan et al. 2012bKhan et al. , 2018)).Koehn et al. (2023) extracted galaxies from the ROMULUS25 cosmological simulations and then resimulated the dynamics of SMBH triplets with direct N -body codes.The question then arises: How to resolve the SMBH binary dynamics in hydrodynamical simulations with a more self-consistent approach, avoiding these oversimplified approximations?This represents yet another challenge to be explored.
To address the questions mentioned above, we introduce the series of studies, 'Resolving supermAssive Black hole Binaries In galacTic hydrodynamical Simulations' (RABBITS).In the RAB-BITS studies, we employ the KETJU code (Rantala et al. 2017) to accurately model the dynamics of SMBH binaries in hydrodynamical simulations of galaxy mergers, including both disc and elliptical galaxies.
Based on the GADGET-3 code (last described by Springel 2005) which includes smoothed particle hydrodynamics (SPH) and galaxy formation subgrid models (Hu et al. 2014;Eisenreich et al. 2017), the KETJU code further incorporates the high-accuracy algorithmically regularized integrator MSTAR (Rantala et al. 2020) andPN corrections (Blanchet 2014).This unique hybrid approach enables us to simultaneously capture galactic-scale galaxy formation processes and resolve the small-scale dynamics of SMBHs.The KETJU code has been successfully used to simulate the complex evolution of SMBHs in galaxy group environments in cosmological zoom-in simulations (Mannerkoski et al. 2021(Mannerkoski et al. , 2022)).Recently, Liao et al. (2023) (hereafter L23) introduced a new accretion model for SMBH binaries to the KETJU code, which helps to improve the modelling of SMBH binaries in gas-rich galaxy mergers.Notably, L23 incorporated both thermal and kinetic AGN feedback, allowing us to explore the impact of different implementations of AGN feedback.
Our focus in this paper (Paper II of the RABBITS series) is on studying how AGN feedback and its numerical implementations affect both galaxy properties and the orbital evolution of SMBH binaries.In a companion paper, Liao et al. (2024) (hereafter Paper I), we focus on the impact of nuclear star formation on the hardening of SMBH binaries.The structure of this paper is as follows.In Section 2, we describe the numerical code, the initial conditions, and the simulation details.The impact of AGN feedback on the properties of merging galaxies is presented in Section 3. In Section 4, we discuss the SMBH orbital decay in both the dynamical frictiondominated phase and the SMBH binary hardening and coalescence phase.We summarize and conclude in Section 5.In Appendices A and B, we study the impact of numerical resolution and perform simulation tests of isolated elliptical galaxies.

Numerical code
All the simulations in this work are performed with the KETJU code (Rantala et al. 2017), which is based on the GADGET-3 code (i.e. the updated version of GADGET-2 described in Springel 2005) and incorporates the high-accuracy algorithmically regularized MSTAR integrator (Rantala et al. 2020) to resolve the SMBH dynamics. 2he code also includes PN corrections up to the order of 3.5PN for SMBH binaries (Thorne & Hartle 1985;Mora & Will 2004;Blanchet 2014), which allows capturing the binary evolution down to the separation regime where the hardening process is predominantly driven by GW emission.The mass loss of the SMBH merger remnant due to GW emission is modelled using the fitting formulae from Zlochower & Lousto (2015).
The KETJU code uses the SPHGAL implementation (Hu et al. 2014) for solving gas hydrodynamics.The adopted SPH kernel is the Wendland C 4 kernel and we set the number of SPH neighbours to N ngb = 100.To model the metal-dependent radiative cooling, star formation, and stellar feedback processes, the 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).The cooling tables are adopted from Wiersma et al. (2009).Stars form from gas particles in a stochastic implementation and the star formation criteria are: gas density ρgas ⩾ 2.2 × 10 −24 g cm −3 (or equivalently hydrogen number density nH ⩾ 1 cm −3 ), gas temperature Tgas ⩽ 1.2 × 10 4 K, and the gas particle must be part of an convergent flow (i.e. the velocity divergence ∇ • vgas ⩽ 0).The stellar feedback model includes the feedback from the Type Ia (SNIa) and Type II (SNII) supernovae explosions and the slow winds of asymptotic giant branch (AGB) stars.
For the SMBH subgrid model, we use the SMBH binary accretion model introduced in L23, 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 circumbinary discs (Duffell et al. 2020).Specifically, when an SMBH is in the single SMBH phase, its accretion rate is computed following the traditional BHL model (Springel et al. 2005). 3When two SMBHs form a bound binary, their total accretion rate is computed using the BHL formula and the gas properties at the centre of mass (CoM) of the SMBH binary, i.e.
where α = 25 is the boost factor, G is the gravitational constant, M bin = MBH,1 + MBH,2 (assuming MBH,1 ⩾ MBH,2) is the total mass of the SMBH binary, MBH,1 and MBH,2 are the masses of the primary SMBH and the secondary SMBH, ρCoM and cs,CoM are the gas density and the sound speed at the CoM position of the binary, and v rel,CoM is the magnitude of the relative velocity between the binary's CoM velocity and the gas velocity at the CoM position of the binary.The total accretion rate is capped by the Eddington limit of the binary system.Then, to distribute the total accretion rate between each SMBH of the binary system, we use the fitting formula from circumbinary disc simulations derived in Duffell et al. ( where q ≡ MBH,2/MBH,1 ∈ (0, 1] is the SMBH mass ratio.As a result, the secondary SMBH gets a higher accretion rate compared to the primary SMBH, which is called preferential accretion and drives the SMBH binary toward equal mass.This binary accretion model provides more physically motivated SMBH mass evolution in the binary phase, which is important for modelling the AGN feedback and the GW induced recoil velocities. To study the impact of different AGN feedback implementations, we consider either pure thermal feedback or pure kinetic feedback in this work.For the thermal AGN feedback model, we follow the implementation in Springel et al. (2005).At each timestep, the total amount of energy given by is added to the internal energy of the surrounding gas particles according to their SPH kernel weights.Here, ϵr = 0.1 is the radiative efficiency, ϵ f,thm = 0.02 is the thermal AGN feedbak efficiency, c is the speed of light, and ∆t is the time-step.
For the kinetic AGN feedback model, we adopt an implementation similar to the kinetic mode in Weinberger et al. (2017).At each time-step, the amount of feedback energy computed by is added to the SMBH feedback energy reservoir, E kin,res .Here ϵ f,kin = 0.008 is the kinetic AGN feedback efficiency.Once the SMBH feedback energy reservoir reaches a given threshold, i.e.
the energy in the reservoir is released to the surrounding N ngb gas particles according to their SPH kernel weights.Here, f thr = 20 is the user-specified threshold parameter controlling the feedback strength, Mgas is the total gas mass within the SMBH smoothing length, and σDM is the one-dimensional velocity dispersion estimated from 40 nearby dark matter particles.The feedback energy, E kin,i , is coupled to the i-th gas neighbour in a kinetic form, i.e. this gas particle receives a kick velocity with the magnitude of where mgas,i is the mass of the i-th gas neighbour.In the nonbinary phase, the direction of the kick velocity has an equal probability (50 per cent chance) of being either parallel or anti-parallel to the direction of the gas particle's angular momentum, ri × vi, where ri and vi are the position and velocity vectors of the gas particle with respect to the SMBH.In the binary phase, the kick direction has an equal 50 per cent probability of being parallel or anti-parallel to the direction of the binary orbital angular momentum.
Note that the thermal AGN feedback is implemented in a continuous approach, while the kinetic AGN feedback is in a pulsed manner.For more details of the galaxy formation subgrid model used in this study and the code implementations, we refer the interested reader to L23.

Initial conditions
In this study, we consider both disc-disc and elliptical-elliptical galaxy mergers.To generate the initial conditions, we first prepare the single progenitor galaxies.For the disc progenitor galaxy, we use the D1 galaxy from L23 but with an initial SMBH mass of MBH = 7.5×10 7 M ⊙ , which matches the observed SMBH massstellar velocity dispersion (MBH-σ⋆) relation from Kormendy & Ho (2013), instead of a seed SMBH mass (105 M ⊙ ) 5 .The SMBH spin is set to zero, and we do not model the evolution of the SMBH spin due to gas accretion in this study.The D1 disc progenitor is created following the method in Springel et al. (2005).It consists of a Hernquist (1990) dark matter halo (with a concentration of c200 = 9 and a spin parameter of λ = 0.033), an exponential disc (with a scale length r D1 ), and a central SMBH.The initial stellar ages and stellar/gaseous metallicities are set following Lahén et al. (2018).The density profiles of the different components of D1 are shown in the left panel of Fig. 1; the component masses and particle numbers are summarized in Table 1.We refer the reader to L23 for more details.
The elliptical progenitor galaxy, which is denoted as E1 here, is designed to be a 'fraternal twin' of D1, i.e. they have similar masses/sizes in different components but differ in galaxy morphology.Such controlled numerical setups enable us to study the impact of galaxy morphology on SMBH mergers.Following the setup of D1, we use the same virial velocity V200 = 200 km s −1 to generate E1, which leads to the same virial mass (M200 = 2.62×10 12 M ⊙ ) and virial radius6 (R200 = 282 kpc).The dark matter halo of E1 is similar to that of D1.For the stellar component, E1 only has a Hernquist bulge, and its total stellar mass is set to be the sum of the stellar disc and bulge masses of D1, and its projected effective radius is set to r E1 b,proj = r D1 d = 4.32 kpc, which corresponds to a 3D scale radius of r E1 b = 2.38 kpc.The ages and metal abundances of star particles are initialized in the same manner as for the D1 bulge component.For the gaseous component, E1 has a hot gas halo with the mass within R200 equal to the gas disc mass of D1, i.e.M gas 200 = 2.15 × 10 10 M ⊙ .The hot gas halo follows the β-model profile, which provides a good fit to the hot Initial elliptical galaxy E1 Figure 1.Spherically averaged density profiles of the initial disc (D1, left) and elliptical (E1, right) galaxies.The dark matter, gaseous, and stellar components are plotted with blue, orange, and green, respectively.For D1, the density profiles of the stellar disc and stellar bulge are shown with the dashed and dotted curves, respectively.The two initial galaxies, D1 and E1, are designed to be 'fraternal twins', i.e. they have identical virial masses, have similar masses in the dark matter, stellar, and gas components, but they differ in their morphology.
Table 1.Properties of single progenitor galaxies.Both D1 and E1 galaxies are generated with a virial velocity of V 200 = 200 km s −1 , which is equivalent to a virial mass of M 200 = 2.62 × 10 12 M ⊙ and a virial radius of R 200 = 282 kpc.Both galaxies contain an SMBH in the centre with a mass of M BH = 7.5 × 10 7 M ⊙ .From the left, the name of the progenitor galaxy, the dark matter mass M DM , the stellar disc mass M ⋆,disc , the stellar bulge mass M ⋆,bulge , the gas disc mass M gas,disc , the gas halo mass M gas,halo , the dark matter particle number N DM , the stellar disc particle number N ⋆,disc , the stellar bulge particle number N ⋆,bulge , the gas disc particle number N gas,disc , and the gas halo particle number N gas,halo .
gas around elliptical galaxies (Forman et al. 1985).We adopt an outer slope parameter β = 2/3 (Jones & Forman 1984) and a core radius of rc = 0.22rs = 6.9 kpc (Makino et al. 1998), where rs ≡ R200/c200 = 31.3kpc is the scale radius of the dark matter halo.The cutoff radius of the gas halo, beyond which there are no gas particles, is set to rcut = 50rc = 345 kpc.The total mass of the gas halo is M gas,halo = 2.64 × 10 10 M ⊙ .We assume no angular momentum transfer between the dark matter halo and the gas halo, and thus the gas halo has the same specific angular momentum as the dark matter halo.The radial metallicity profile of the gas halo is assumed to follow the stellar one.Following D1, E1 also contains a central SMBH with a mass of MBH = 7.5 × 10 7 M ⊙ and a spin of zero.Note that Eisenreich et al. (2017) adopted a similar approach to set up their initial elliptical galaxy.The density profiles of E1 are shown in the right panel of Fig. 1; the masses and particles numbers of different components are summarized in Table 1.
We have compared the properties of D1 and E1 to the observational data, and confirmed that both galaxies agree reasonably well with the observations; see Section 3 for details.The D1 and E1 progenitor galaxies are then used to create the initial conditions of equal-mass galaxy mergers.For both disc-disc and elliptical-elliptical galaxy mergers, two identical galaxies are set to approach each other on a parabolic orbit with an initial separation of dsep,ini = R200 and a pericentric distance of dperi,ini = 2r D1 d .We adopt the G5 retrograde galaxy orbit geometry from Naab & Burkert (2003) for the disc-disc galaxy merger, resulting in a relatively modest starburst during the merger and making the simulation computationally efficient.For the G5 orbit, the inclination of the disc with respect to the orbital plane, i, and the argument of pericentre, ω, of the two galaxies are i1 = −109 • , ω1 = −60 • , i2 = 180 • , ω2 = 0 • .Note that here the definitions of i and ω follow the convention in Toomre & Toomre (1972).Following the philosophy of fraternal twin simulations, for the elliptical-elliptical galaxy merger, as the rotation of the progenitor galaxy (e.g. the dark matter halo and the gas halo) introduces a special direction, we set the initial galaxy orbit also as G5 by replacing the disc orientation with the galaxy spin direction and the disc plane with the halo equatorial plane.
Following the naming convention of galaxy mergers in L23, i.e. 'progenitor galaxy types-galaxy mass ratio-orbit geometry', we denote our disc-disc and elliptical-elliptical galaxy mergers as DD-Table 3. Simulation sets including different physical processes.The symbols '✓' (Yes) and '✗' (No) are used to denote whether a physical process is included in a simulation.Here, the hydrodynamics is solved using the SPH method.The radiative cooling process takes metal-dependent cooling into account.Star formation criteria are based on the density and temperature thresholds.Stellar feedback includes the feedback from SNIa, SNII, and AGB stars.SMBH accretion is modelled using the BHL accretion and the preferential accretion from circumbinary discs.The thermal AGN feedback is performed by increasing the internal energy of the surrounding gas particles, while the kinetic AGN feedback is implemented by adding velocity kicks to the ambient gas.Further details regarding various physical processes are provided in Section 2.1.

Name
11-G5 and EE-11-G5, respectively.The particle masses and the total particle numbers of the simulation initial conditions are summarized in Table 2.

Simulations
To study how different physical processes (especially AGN feedback) affect the properties of merging galaxies and the orbital decay of SMBH binaries, for both DD-11-G5 and EE-11-G5 galaxy mergers, we perform five simulation sets including different physical processes and model implementations: (i) NoGas.This simulation set is gravity-only.In the initial condition, all gas particles are converted into stars by inheriting the positions, velocities, and masses.This set of simulations is an analogue to the collisionless gas-free elliptical galaxy mergers in previous studies.Note that here we do not refer to these simulations as 'collisionless' runs because they do resolve the SMBH-SMBH and the SMBH-star collisional interactions.
(ii) NoCool.Compared to the NoGas set, this set of simulations includes not only gravity but also gas hydrodynamics with an adiabatic equation of state.These runs do not incorporate the radiative cooling process, and as a result, the processes of star formation, as well as stellar and AGN feedback, are not taken into account.Note that similar runs have been referred to as 'adiabatic' or 'non-radiative' in the literature.
(iii) CoolStarNoAGN.Compared to the NoCool set, these simulations further include the subgrid models of gas cooling, star formation, and stellar feedback.However, SMBH accretion and AGN feedback processes are not considered in these simulations.
(iv) CoolStarThmAGN: Compared to the CoolStarNoAGN set, this set of simulations includes in addition the SMBH binary accretion model and the thermal AGN feedback mode.
(v) CoolStarKinAGN: These simulations are similar to the CoolStarThmAGN ones, but use the kinetic AGN feedback mode instead of the thermal mode.
A summary of the physical processes included in the different simulation sets is given in Table 3.All simulations in this work start from the cosmic time of t0 = 10.7 Gyr, with most evolving for 3 Gyr.However, the SMBH binaries in the DD-11-G5 NoGas and the EE-11-G5 CoolStarThmAGN runs take longer to merge, so these two simulations are evolved for 4 Gyr.For a straightforward comparison with the EE-11-G5 CoolStarThmAGN run, the EE-11-G5 CoolStarKinAGN run is also evolved for 4 Gyr, despite the fact that the SMBH merger takes place at ∼2.95 Gyr.
Following L23, the Plummer-equivalent gravitational soften-ing lengths for dark matter, gas, star, and SMBH particles are set to 100 pc, 20 pc, 5 pc, and 5 pc, respectively.Note that only the SMBH-star and SMBH-SMBH interactions are unsoftened in our simulations, the SMBH-gas and SMBH-dark matter gravitational interactions are still softened.The radius of the KETJU region around each SMBH is 15 pc, i.e. three times the stellar/SMBH softening length.Compared to the runs in L23 which start with a seed SMBH mass and the KETJU integration is only switched on when the SMBH-to-star particle mass ratio reaches the threshold MBH/m⋆ = 300, the simulations in this work have an initial mass ratio of MBH/m⋆ = 750, and the KETJU integration is switched on from the very beginning.In this paper, the simulation time, t, is defined as the time span after t0.The CoolStarThmAGN and CoolStarKinAGN simulations utilize the same numerical setup and parameters as in L23, thus we refer the interested reader to L23 for further simulation details.

Overview of galaxy evolution
We start by looking at the overall galaxy merger process.In prior to the first pericentre passage), while the extended dark matter haloes of the two galaxies are already interacting with each other, the central regions of the galaxies, including the central SMBHs, remain tightly bound and experience minimal tidal interactions.As a result, they continue to follow the input parabolic orbits and move ahead of their own dark matter haloes.When the galaxy centres get closer and go through the first pericentre passage, the intense tidal interactions between discs deform their appearances and produce elongated and narrow tails and bridges (Toomre & Toomre 1972).The gaseous tidal features tend to be narrower and better defined compared to the stellar ones, reflecting that radiative cooling in gas helps to damp the random motions.At the same time, the strong gravitational pulls from the lagging dark matter haloes slow down the galaxy centres.At t ∼ 1 Gyr, the galaxy centres reach their first orbital apocentre and fall back.During this process, the tidal forces convert the orbital energy and angular momentum into the internal degrees of freedom in each galaxy, and consequently the galaxy centres deviate from the initial parabolic orbits.The centres of galaxies undergo a few passages, and with each cycle, the associated tidal interactions increasingly stimulate the internal motions of the inner regions.As a result of the continuous loss of orbital energy and angular momentum, the two galaxies/SMBHs eventu-ally merge together and the galaxy merger remnant continues to undergo relaxation.
The overall orbital evolution of the EE-11-G5 case is qualitatively similar as the aforementioned disc system.At the early stage, the galaxy centres also follow the initial parabolic orbits before their first pericentre passage.However, unlike the cold disc system, the spherical elliptical galaxies do not exhibit narrow tidal tails or bridges, instead, we only observe broad features, i.e. some stars are excited to wider orbits by the tidal interactions and the overall stellar distribution is more extended.The physical reason behind these broad tidal features is that the random motions of individual stars in such dispersion-supported systems are comparable to the velocities imparted by the tidal forces.The initial spherical hot gas haloes also undergo deformation due to tidal torques and particularly the additional gas pressure, leading to the dispersal of their high-density cores.

Star formation and SMBH growth histories
In Fig. 3, we present the star formation and SMBH growth histories for the CoolStarNoAGN, CoolStarThmAGN, and CoolStarKi-nAGN runs, for both disc and elliptical galaxy mergers.Note that the NoGas and NoCool runs do not include the processes of star formation and SMBH mass growth, hence they are excluded from this figure.
In the top panels, the total star formation rates (SFRs) in the simulations as a function of time are plotted.The CoolStarNoAGN runs consistently exhibit the highest SFRs throughout the entire simulation for both the DD-11-G5 and EE-11-G5 mergers.However, once thermal AGN feedback is included, star formation is suppressed, particularly after t ∼ 1.5 Gyr, when the SMBHs exhibit high accretion rates and produce strong AGN feedback.Notably, the implementation of kinetic AGN feedback proves to be more effective in suppressing star formation compared to thermal AGN feedback.In the case of the DD-11-G5 merger, the Cool-StarKinAGN run exhibits a lower total SFR compared to the Cool-StarThmAGN run throughout the entire simulation period.In the last ∼0.5 Gyr, the SFR of the CoolStarKinAGN merger remnant becomes approximately two orders of magnitude lower than that of the CoolStarThmAGN remnant.In the EE-11-G5 merger, the CoolStarKinAGN run results in a complete quenching of the SFR in both merging galaxies after the first pericentre passage (t ∼ 0.7 Gyr).In contrast, the CoolStarThmAGN run maintains an SFR ≳ 0.1 M ⊙ yr −1 for most of the simulation duration.These results are in agreement with previous simulation studies which also found that kinetic AGN feedback is more efficient in removing gas from the galaxy centre and thus quenching galaxies (e.g.Choi et al. 2012Choi et al. , 2014;;Barai et al. 2014;Eisenreich et al. 2017; see also section 5.1 of Costa et al. 2020 for a detailed discussion).
The impact of AGN feedback on gas can be further illustrated in the gas phase diagram.In Fig. 4, we show the gas temperaturedensity phase diagrams at t = 2 Gyr, when the SMBHs either have recently merged or have formed a bound binary.The horizontal and vertical solid lines mark the star formation temperature and density thresholds, respectively.The gas particles in the bottom-right corner are star-forming.Compared to the no AGN case, the thermal AGN feedback heats up more gas particles in the dense region (i.e.ρgas ≳ 10 −24 g cm −3 ).In contrast, the kinetic AGN feedback kicks these gas particles with high velocity, causing them to leave the dense region and collide with gas particles at larger radii, thereby heating the outer diffuse gas.Consequently, the CoolStarK-inAGN run tends to have more gas particles in the hot and diffuse phase compared to the no AGN feedback and the thermal AGN feedback cases.In the EE-11-G5 remnant, the kinetic AGN feedback is so effective that no star-forming gas is present, leading to a fully quenched galaxy remnant.
The evolution of the total SMBH accretion rates and the total SMBH masses are displayed in the middle and bottom panels of Fig. 3, respectively.In the DD-11-G5 case, the thermal and kinetic AGN feedback runs show similar overall trends in SMBH accretion histories.However, the CoolStarKinAGN run, with its pulsed feedback, displays more bursty SMBH accretion (or equivalently higher AGN variability) compared to the continuous thermal feedback.Following the merger of the two SMBHs, we observe a slight decrease in the total SMBH mass, attributed to the energy loss through GW emission (Zlochower & Lousto 2015).In the EE-11-G5 case, after the first pericentre passage, the CoolStarKinAGN run exhibits SMBH accretion rates approximately two orders of magnitude lower than those of the CoolStarThmAGN run.Consequently, the total SMBH mass in the CoolStarKinAGN run shows minimal growth, reaching only half of that in the CoolStarThmAGN run by the end of the simulation (t = 4 Gyr).This low SMBH accretion rate in the kinetic AGN feedback run can be understood as follows: unlike the cold gas reservoir in disc galaxies, the gas from the extended hot gas halo in elliptical galaxies must first cool down and reach the galaxy centre to be accreted by the SMBH.Since kinetic AGN feedback effectively maintains the gas heated to a large distance, preventing the formation of cooling flows and keeping a low gas density in the galaxy centre, the central SMBH experiences a considerably lower accretion rate.
We have conducted a resolution study to address the reliability of our results by running simulations with both higher and lower mass resolutions.This study confirms that the star formation and SMBH growth histories presented above have reached a satisfactory resolution convergence.See Appendix A for details.
From the comparison among the CoolStarNoAGN, Cool-StarThmAGN, and CoolStarKinAGN runs, we can clearly see that the presence and specific numerical implementation of AGN feedback have significant impact on star formation and SMBH growth.The kinetic AGN feedback is necessary to efficiently quench galaxy merger remnants, which typically become elliptical galaxies.In Appendix B, we further test the AGN feedback implementations in maintaining a quiescent elliptical galaxy in isolation, and the results again emphasize the importance of kinetic AGN feedback in generating a red and quiescent early-type galaxy.

Galaxy scaling relations
In Fig. 5, we provide a comparison of the galaxy properties between the merger remnants from the CoolStarNoAGN, Cool-StarThmAGN, and CoolStarKinAGN simulations, as well as the properties of observed galaxies.The properties of the initial galaxy are also included for reference (shown in cyan).From Panels (a) to (f), we present the MBH-σ⋆ relation, the galaxy size-stellar mass relation, the hot gas X-ray luminosity, the stellar metallicity-stellar mass relation, the stellar mass-halo mass relation, and the dark matter mass fractions within the projected stellar effective radius.
By comparing the properties of our initial galaxies with the properties of the observed galaxies, we find that both our disc and elliptical progenitor galaxies match the observations well.It is worth noting that in Panel (c), which computes the soft X-ray luminosity (0.3-5 keV) emitted by the hot gas, the initial disc galaxy is not shown as it only has a cold gas disc, resulting in an X-ray luminosity value of zero.
We first study the merger remnants from the DD-11-G5 runs, which are represented by filled circles in all panels.In the absence of AGN feedback, the CoolStarNoAGN run exhibits an excessive level of star formation in the galaxy centre, which leads to an elevated σ⋆, the line-of-sight stellar velocity dispersion within the projected half stellar mass radius.Combined with the lack of SMBH accretion, this consequently leads to the deviation from the best- In Panel (c), the hot gas X-ray luminosity of the CoolStarKinAGN remnant is approximately one order of magnitude lower than the other two runs.This is attributed to the enhanced efficiency of kinetic AGN feedback in expelling the gas from the centre, resulting in reduced central gas density and a more extended hot gas halo (see Choi et al. 2014, for similar conclusions).From Panel (d) to (f), the three merger remnants exhibit more similar properties.This similarity arises because the stellar properties within the 30 kpc aperture are less affected by the specific distribution in the very central region.
The merger remnants from the EE-11-G5 runs are plotted with filled diamonds.Unlike the aforementioned DD-11-G5 runs, the stellar properties from all three EE-11-G5 runs are remarkably similar.This similarity comes from the lower overall star formation in elliptical mergers, resulting in negligible differences among the different runs.The most prominent difference among the three runs is found in the hot gas X-ray luminosity, as shown in Panel (c).Specifically, the X-ray luminosity of the CoolStarKinAGN remnant is nearly two orders of magnitude lower than that of the Cool-  The colour represents the number of gas particles in the different temperature-density bins, and the colour bar for each row is displayed on the right.The horizontal and the vertical solid lines mark the temperature and the density thresholds for star formation, respectively.The dashed line marks the Jeans mass resolution limit, M J,lim = N ngb mgas.For gas particles below this dashed line, the Jeans mass is not resolved and they can collapse under self-gravity due to numerical noise (see Eisenreich et al. 2017).Compared to the cases of no AGN feedback and thermal AGN feedback, the kinetic AGN feedback model results in distinct gas distributions in the phase diagram, i.e. less gas particles in the star-forming phase and more gas in the hot and diffuse phase.
StarThmAGN remnant, further indicating that kinetic AGN feedback is more effective in redistributing the gas from the galaxy centre and thus reducing X-ray luminosity (see also Choi et al. 2015).
It is evident from Fig. 5 that for both the DD-11-G5 and EE-11-G5 mergers, the runs with AGN feedback included (i.e. both CoolStarThmAGN and CoolStarKinAGN) generate merger remnants which closely resemble the observed galaxy properties.In Panel (b), some galaxy remnants (e.g. the DD-11-G5 ones) tend to lie below the best-fitting observed relations.This behaviour is anticipated in idealized merger simulations, as the galaxy remnants in our study have undergone a single major merger event.In contrast, real galaxies in the Universe often experience additional minor mergers, which can contribute to an increase in their sizes (Naab et al. 2009;Johansson et al. 2012).In Panel (d), the galaxy remnants tend to slightly exceed the observed median stellar metallicity-stellar mass relation.This can be attributed in part to the initial galaxy setup and is further explained by the absence of fresh gas inflow from the intergalactic medium in idealized merger simulations, which would otherwise dilute the metallicity of the interstellar medium and affect the stars formed.In Panel (e), our galaxy remnants tend to be positioned above the stellar mass-halo mass relations (e.g.Moster et al. 2013;Kravtsov et al. 2018), which again partially originates from the initial galaxy setup, but the remnants still fall within the range of two times the intrinsic scatter.
We can conclude that our simulations with AGN feedback reproduce reasonably well the observed galaxy properties.The agreement between the simulations and observations demonstrates that the SMBH subgrid model developed in L23 for disc galaxy mergers can be applied to elliptical galaxy mergers as well.

SMBH ORBITAL DECAY
In Section 3, we demonstrated that simulations incorporating different physical processes (e.g. with versus without SMBH accretion and AGN feedback) or employing different numerical implementations (e.g.thermal versus kinetic AGN feedback) yield different galaxy properties, including different central stellar properties, gas distributions, and SMBH masses.Such variations in galaxy properties are expected to result in different rates of orbital decay for the merging SMBHs, influenced by distinct dynamical friction and gravitational potential environments.In this section, we quantitatively compare the SMBH orbital decays across different simulation runs.
Following Merritt (2013), in the subsequent analyses, we define the SMBH influence radius as the radius where the circular velocity around an SMBH equals the one-dimensional (or line-ofsight) stellar velocity dispersion, When the separation of two SMBHs reaches roughly the influence radius of the more massive SMBH, ∆dBH ∼ R infl , the two SMBHs form a bound binary.Usually, this is regarded as the end of the dynamical friction-dominated phase, after which the interactions between star particles and the SMBH binary (i.e.three-body interactions) start to play a role in the SMBH orbital decay.The hard binary separation is defined as (Merritt 2013) where M red = MBH,1MBH,2/M bin is the reduced mass.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.When the semimajor axis of the SMBH binary reaches the hard binary separation, a ∼ R hard , the SMBH Figure 5.Comparison of galaxy properties in simulations and observations.In all panels, the initial isolated galaxies are plotted with cyan, while the merger remnants in the CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are shown with green, red, and purple, respectively.The DD-11-G5 and EE-11-G5 runs are plotted with filled circles and diamonds, respectively.(a) M BH -σ⋆ relation.The error bars in the line-of-sight stellar velocity dispersion show the standard deviations computed from 50 random projected orientations.The solid line plots the best-fitting M BH -σ⋆ relation from van den Bosch (2016) and the grey-shaded regions show the intrinsic one-sigma scatter.The dashed lines give the best-fitting relations from Kormendy &Ho (2013) andTremaine et al. (2002).(b) Galaxy size-stellar mass relation.Here the galaxy size is characterized using the projected half stellar mass radius Re.The solid lines and shaded regions show the observed relations for early-type galaxies (ETGs) and late-type galaxies (LTGs) from Shen et al. (2003) and Newman et al. (2012).(c) The luminosity of soft X-rays (0.3-5 keV) from hot and diffuse gas versus the galaxy stellar mass.The X-ray luminosity is computed following the method detailed in L23.Note that the initial disc galaxy is not shown in this panel because it only has a cold gas disc and its hot gas X-ray luminosity is thus zero.The observational data from the MASSIVE (Ma et al. 2014) and ATLAS3D (Cappellari et al. 2011) 2020).Overall, the galaxy merger remnants in the simulations with AGN feedback included agree well with the observations.orbital decay undergoes a complete transition from dynamical friction to being entirely dominated by three-body interactions.In our simulations, typically R infl are ∼10 pc while R hard are ∼1 pc.

Overall evolution of SMBHs and their embedded galaxy centres
In Fig. 6, we present the time-evolution of the SMBH separation in the first 2 Gyr for all simulation runs.In the left panels, we mark the SMBH influence radius, R infl , and the hard binary separation, R hard , from different runs.
We first focus on the DD-11-G5 mergers, which are summarized in the top panels.The evolution of the different runs during the first ∼1 Gyr are quite similar as the plotted lines are almost on top of each other.As discussed in Section 3.1, initially the two SMBHs and the bound galaxy centres that they are embedded in move together following the input parabolic orbits until they experience strong tidal interactions.During the first pericentre passage, the tidal effects are not identical for galaxies in different runs   as they have distinct galaxy compositions (i.e.NoGas versus other runs) and they have experienced distinct star formation and SMBH growth histories (i.e.runs except NoGas and NoCool; see Fig. 3) which lead to different central galaxy properties.As a result, the time-evolution of the SMBH separation in different runs start to deviate from each other and the differences become more noticeable after the first apocentre (e.g. at t ∼ 1 Gyr).After the first apocentre, the NoGas run exhibits the slowest orbital decay, followed by the NoCool run, the CoolStarKinAGN run, the CoolStarThmAGN run, and finally the CoolStarNoAGN run which shows the most rapid orbital decay.
The comparison between the NoGas and the NoCool runs showcases the impact of gas, i.e. the gas pressure assists in the conversion of orbital energy and angular momentum into internal motions, resulting in a more rapid SMBH orbital decay.Similar phenomenon was observed in Barnes & Hernquist (1996) for their collisionless and gaseous paired runs of a retrograde disc galaxy merger.
Compared to the NoCool run, the runs with gas cooling and star formation (i.e.CoolStarKinAGN, CoolStarThmAGN, and CoolStarNoAGN) display even more rapid orbital decay, which is an outcome of the more bound galaxy centres.The enclosed total mass profiles (including SMBH, star, gas, and dark matter) at t = 1 Gyr from different runs are shown in the top-left panel of Fig. 7.We can easily see that the central matter distributions tend to be more concentrated in those runs including gas cooling and star formation.The enclosed mass profiles for the new stars and gas are plotted in the bottom-left panel of Fig. 7. Here, the new stars are defined as the stars formed after the start of the simulation, i.e. stars with a formation time tSF ⩾ t0.The CoolStarNoAGN run has a significant amount of new stars formed in the very inner region (i.e.r ≲ 1 kpc), while the AGN feedback suppresses the star formation in the very inner region in the CoolStarThmAGN and the CoolStarKinAGN runs.Compared to the NoCool run, the gas distributions in the CoolStarThmAGN and the CoolStarKinAGN runs are less concentrated due to two factors, i.e. some inner gas particles have already been converted into stars and more importantly the AGN feedback has heated and pushed the gas to the outer region.
The more concentrated galaxy centres are more bound, and thus they can resist the tidal stripping effect more effectively and sustain higher masses for a longer time.As a result, they experience stronger dynamical friction and their SMBH orbits decay faster.This explains why the CoolStarNoAGN run exhibits the most rapid SMBH orbital decay, followed by the CoolStarThmAGN and the CoolStarKinAGN runs, and why these three runs show more rapid decays than the NoCool and the NoGas runs.This picture is in line with the work of Callegari et al. (2009), who found that in minor mergers of disc galaxies, incorporating radiative cooling and star formation creates more bound satellite galaxy centres and consequently enhances the SMBH pair formation, in contrast to dry merger simulations without gas.The recent work of Chen et al. (2023) has also suggested the important role of high central stellar density in assisting the orbital decay of SMBHs.
The EE-11-G5 mergers are summarized in the bottom panels of Fig. 6.Similar to the DD-11-G5 case, during the first ∼1 Gyr, the evolution of SMBH separations from different EE-11-G5 runs are almost identical, since the SMBHs and their embedded galaxy centres are moving together and the bound masses are very similar across all runs.In addition, we notice that the evolution during the first ∼1 Gyr are fairly similar between the DD-11-G5 and the EE-11-G5 runs, as a result of the designed fraternal twin initial conditions.
After the first apocentre, we start to observe some noticeable differences among the different runs.However, overall the differences among the different runs are smaller than those in the DD-11-G5 runs, because the star formation and SMBH growth are much weaker in the EE-11-G5 runs and the enclosed mass profiles are more similar for the different runs (see the right panels of Fig. 7).
In the EE-11-G5 runs, the largest difference in the SMBH separation evolution is between the NoGas run and the other runs that include gas.Interestingly, in contrast to the DD-11-G5 case, here the SMBHs in the NoGas run reaches their second pericentre passages earlier than the runs including gas.This is caused by the fact that we convert all gas particles in the hot gas halo into star particles in the initial condition for the NoGas run, and the evolution of this 'stellar halo' is distinct from that of the gas halo, creating different gravitational drags on the galaxy centres.We have performed a test run by simply removing the gas particles in the hot gas halo to create an initial condition, which is dubbed 'NoGas test' and is plotted with grey in the bottom panel of Fig. 6.We can see that without this 'stellar halo', the SMBHs reach their second pericentre passage later compared to the runs including gas, in line with the DD-11-G5 case.For the sake of completeness, we also plot the 'NoGas test' run for the DD-11-G5 merger in the top panel of Fig. 6, illustrating how the least concentrated galaxy centre results in the slowest orbital decay rate.
Among all EE-11-G5 runs, the CoolStarNoAGN run has the most star formation in the galaxy centre, creating the most bound central regions and consequently having its SMBH separation reaching R infl first (i.e.exhibiting the most rapid orbital decay).In Section 4.2, we will show that the gas dynamical friction also contributes to the SMBH orbital decay in the CoolStarNoAGN case.The CoolStarThmAGN run has the second most rapid orbital decay, followed by the NoCool and the CoolStarKinAGN runs.As the star formation is completely quenched after the first pericentre passage and kinetic AGN feedback is effective in maintaining a hot gas halo, the CoolStarKinAGN run displays relatively similar SMBH orbital decay as the NoCool run.
Overall, the orbital decay of SMBHs in the EE-11-G5 case is slower than that in the DD-11-G5 case.This is because compared to the DD-11-G5 runs, the EE-11-G5 runs exhibit less star formation and less SMBH growth, resulting in less bound galaxy centres, weaker dynamical friction, and slower SMBH orbital decays.

Dynamical friction on SMBHs: stars versus gas
As the separation between the SMBHs decreases, the matter initially bound to the galaxy centre continues to be tidally stripped away.During the last stage of the dynamical friction-dominated phase, most of the matter that was initially bound to the SMBHs has been stripped away, leaving the two SMBHs effectively as isolate massive point particles moving within the stellar and gaseous background.
To address the relative importance of stellar and gaseous dynamical friction, we first look at the stellar (solid lines) and gas (dashed lines) density profiles when the SMBH separations are approximately twice the SMBH influence radius (i.e.∆dBH ∼ 20 pc).These density profiles are presented in the top panels of Fig. 8.Note that here the profiles are centred on the CoM position of the two SMBHs.We can clearly see that in all simulations including gas, the inner mass distributions are dominated by stars.In both the DD-11-G5 and EE-11-G5 cases with AGN feedback, the SMBHs push gas away, leading to a significant reduction in the gas density within ∼100 pc.The kinetic AGN feedback in the EE-11-G5 case is particularly effective in achieving this.Without AGN feedback, one would expect gas to cool down and condense in the galaxy centre, resulting in a higher gas central density, as observed in the EE-11-G5 CoolStarNoAGN run.However, in DD-11-G5 CoolStarNoAGN, the unrealistically massive stellar centre produces strong supernova feedback which pushes the gas to larger radii, resulting in a gas density profile similar to the CoolStarTh-mAGN one.
The central density has an important impact on the dynamical friction.To further provide a quantitative comparison, we estimate the stellar and gaseous dynamical frictions from our simulations.The stellar dynamical friction can be modelled as (Chandrasekhar 1943;Binney & Tremaine 2008) where ρ⋆ is the central stellar density, v rel,⋆ is the magnitude of the relative velocity between the SMBH and the stellar background, and Here, ln Λ is the Coulomb logarithm, X = v rel,⋆ /( √ 2σ⋆) is the ratio between the SMBH-star relative velocity and the stellar velocity dispersion, and erf(X) is the error function.On the other hand, the gas dynamical friction can be estimated by (Ostriker 1999) where ρgas is the central gas density, v rel,gas is the magnitude of the relative velocity between the SMBH and the gas background, and the velocity-dependent factor, Igas, has different forms in the subsonic and the supersonic regimes, i.e.
Here, M = v rel,gas /cs is the Mach number.Therefore, the ratio between the stellar and the gaseous dynamical frictions is It was shown in Ostriker (1999) that, if ρgas = ρ⋆, v rel,gas = v rel,⋆ , and cs = σ⋆, then, in the supersonic (subsonic) regime, the gaseous dynamical friction is stronger (weaker) than the stellar one.
To estimate the force ratio, in our calculations, we use the star particles within the radius7 of 100 pc centring on the SMBH to compute the related stellar quantifies, i.e. ρ⋆ is estimated as the mean stellar density within 100 pc, v rel,⋆ is the relative velocity between the SMBH and the CoM velocity of the used star particles, and σ⋆ is their velocity dispersion.The gas related quantities, ρgas, v rel,gas , and cs, are estimated at the SMBH position using the SPH approach.The Coulomb logarithm is estimated as (Binney & Tremaine 2008) where V 2 200 = GM200/R200 is the squared virial velocity, and we have adopted M200 ≈ 3 × 10 12 M ⊙ and MBH ≈ 10 8 M ⊙ .The detailed results are summarized in Table 4.
In most runs, except for DD-11-G5 CoolStarKinAGN and EE-11-G5 CoolStarNoAGN, the gas density is significantly lower than the stellar density and the SMBH motion is subsonic (i.e.M ∼ 0.4), therefore the stellar dynamical friction dominates over the gaseous one.In DD-11-G5 CoolStarKinAGN, even though the SMBH moves in the supersonic regime (M ∼ 1.7), the stellar density is much higher than the gas density, leading to the dominance of stellar dynamical friction.Interestingly, in EE-11-G5 CoolStarNoAGN, due to weak stellar feedback and the absence of AGN feedback, the central gas density becomes comparable to the stellar density.Furthermore, the SMBH moves in the supersonic regime, resulting in the dominance of gas dynamical friction.Note that we have performed a similar analysis in later snapshots (with ∆dBH ∼ 5-10 pc), and the results are qualitatively similar.
From these comparisons, we can conclude that when AGN feedback is included, the gas density around the SMBHs is drastically reduced, resulting in a decrease in gas dynamical friction on the SMBHs.Consequently, stellar dynamical friction becomes the dominant force affecting the SMBH motion.Importantly, this conclusion is applicable to both the thermal and kinetic AGN feedback implementations.
Our conclusions echo the previous works of e.g.Mayer et al. (2007) and Souza Lima et al. (2017).Mayer et al. (2007) demonstrated that in a coplanar and prograde merger of two equal-mass disc galaxies, gas funnelled by tidal torques forms a massive nuclear disc in the remnant centre, leading to SMBH orbital decay driven primarily by gas dynamical friction instead of stellar dynamical friction.We attribute the main cause of this to the lack of AGN feedback in their simulations, resulting in excessive gas accumulation in the remnant centre, and the lack of feedback from supernova explosions which has further exacerbated this effect.In contrast, Souza Lima et al. (2017) showed that once AGN feedback is included, the spherically distributed feedback can generate a hot bubble around the SMBH and consequently shuts off gas dynamical friction, a phenomenon named 'wake evacuation' by the authors.Recently, Bollati et al. (2023b) also concluded from their simulations that the inclusion of radiative AGN feedback can lead to inefficient or even reversed gas dynamical friction, delaying the orbital decay of SMBH pairs compared to the no feedback case.Our results also agree with those of Pfister et al. (2017), who found that dynamical friction from stars predominantly drives the orbital decay of SMBHs in their galaxy merger simulations incorporating AGN feedback, whilst the contribution from gas is negligible.
It is worth emphasizing that our comparison of stellar and gaseous dynamical frictions above is based on the dynamical friction models from Binney & Tremaine (2008) and Ostriker (1999).In our simulations, the dynamical friction acting on SMBHs from stars is accurately accounted for, thanks to the precise computation of non-softened gravity between the SMBHs and star particles using a regularized approach (Karl et al. 2015;Rantala et al. 2017).In contrast, the gravity between the SMBHs and gas is softened in the KETJU code, which can potentially lead to the underestimation of gas dynamical friction since the gas particles with impact factors lower than 2.8 times the softening length are not properly modelled.However, this does not affect our conclusion that the stellar dynamical friction dominates over the gaseous one in our runs with AGN feedback, as AGN feedback plays a more significant role in affecting the gas density and thus the gas friction.In addition, we note that the accretion of gas mass and momentum into SMBHs in CoolStarThmAGN and CoolStarKinAGN runs also contribute to the SMBH orbital decay.
Approximately when the separation of the two SMBHs reaches the influence radius of the more massive SMBH, they form a bound binary.As can be seen from Fig. 6, in all runs, the SMBH separations drop very rapidly in the regime of R hard ≲ ∆dBH ≲ R infl , as a result of the combined effects from dynamical friction and the gravitational slingshot interactions between stars and the SMBH binary.When the separation (or the binary semimajor axis) reaches ∼R hard , a hard SMBH binary forms and its subsequent orbital decay is discussed in the following subsection.

SMBH binary hardening and GW emission phases
In Fig. 9, we plot the time-evolution of the PN-corrected orbital parameters (Memmesheimer et al. 2004;Mannerkoski et al. 2019), specifically the semimajor axis a and the eccentricity e, of the SMBH binaries from the different runs.The filled circles mark the time t = t hard when a binary becomes hard, i.e. when the semimajor axis of an SMBH binary reaches approximately the hard binary separation R hard ∼ 1 pc.
After t hard , during the gravitational slingshot interaction phase, the central stellar properties have an important impact on driving the binary hardening (e.g.Quinlan 1996;Quinlan & Hernquist 1997).The density profiles centred on the CoM position of the SMBH binary are plotted in the bottom panels of Fig. 8.For both disc and elliptical galaxy mergers, the central stellar densities in the CoolStarNoAGN runs are significantly higher than those in other runs, therefore, it is unsurprising to find that the SMBH binaries in both the DD-11-G5 and EE-11-G5 CoolStarNoAGN runs exhibit the shortest merger time-scales.Especially, for the DD-11-Table 4. Comparison between stellar and gaseous dynamical friction.As the results for the two SMBHs are very similar, here we only present the computed properties of one SMBH.The information is derived from the snapshots shown in the top panels of Fig. 8 (i.e. when ∆d BH ≈ 2R infl ).From the left, the galaxy merger name, the simulation setup, the gas density at the SMBH position, the SMBH-gas relative velocity, the sound speed at the SMBH position, the Mach number, the mean stellar density within 100 pc, the relative velocity between the SMBH and the CoM velocity of star particles within 100 pc, the stellar velocity dispersion computed from star particles within 100 pc, and the ratio between the estimated stellar and gaseous dynamical friction forces.

Galaxy
Simulation  Figure 10.Axis ratios of the galaxy merger remnants as a function of distance from the CoM of the SMBH binary at t = t hard .The upper and lower panels show the results from the DD-11-G5 and EE-11-G5 simulations, respectively.From left to right, the shape profiles of the stellar, dark matter, and gas components are plotted.For each component, the top and bottom panels show the intermediate-to-major axis ratios (b/a) and the minor-to-major axis ratios (c/a), respectively.To ensure a robust estimation of the shapes, we plot the profiles only for radii where the enclosed particle number exceeds 300 (Bett et al. 2007).The merger remnants in all runs exhibit triaxiality in the different components all the way to the centre.The inclusion of gas cooling, stellar physical processes, and AGN feedback has a significant impact on the triaxiality of the matter distribution, especially in disc galaxy mergers and for the gas component.
G5 CoolStarNoAGN case, the SMBH merger time-scale is as short as ∼10 Myr.
The departure from the spherical potential can also lead to rapid binary hardening, which has been regarded as a possible solution to the final parsec problem (Milosavljević & Merritt 2003), i.e. stellar orbits in a triaxial potential can refill the loss cone more efficiently compared to the case of a spherical potential, resulting in efficient hardening of SMBH binaries (Merritt & Poon 2004;Berczik et al. 2006;Khan et al. 2011;Preto et al. 2011).In Fig. 10, we show the axis ratio (including the intermediate-to-major axis ratio, b/a, and the minor-to-major axis ratio, c/a) profiles of different components for our galaxy remnants at t hard .The axis ratios are computed following the iterative method described in Katz (1991).For a given radius of r, we first use all particles within r to compute the moment of inertia tensor and obtain the axis ratios from the tensor's eigenvalues.Then, new axis ratios are computed using only the particles enclosed within the ellipsoidal volume with the previously determined axis ratios.This procedure is repeated for the ellipsoidal volume with the new axis ratios, and iterated until the axis ratios converge.
From Fig. 10, we can see that the merger remnants in all runs exhibit triaxiality all the way to the centre.The NoGas and No-Cool runs show quite similar axis ratio profiles in the stellar and dark matter components.Conversely, in the other runs, the galaxy formation processes (especially AGN feedback) show a significant impact on the triaxial matter distributions, especially in disc galaxy mergers and for the gas component.With weaker AGN feedback, the stellar and gas components tend to show more triaxial distributions at smaller radii, as a result of more efficient cooling and star formation.At larger radii (i.e.≳ 1 kpc), compared to the NoGas and NoCool cases, the energetic galaxy formation processes have resulted in more dynamically hot and thus less triaxial distributions in the stellar and dark matter components.Overall, the disc galaxy merger remnants tend to be more triaxial compared to the elliptical ones.
The eccentricity of an SMBH binary when it becomes hard, e hard , is another crucial factor affecting the merger time-scale.During the gravitational slingshot interaction phase, the eccentricity only shows marginal evolution, i.e. usually it tends to grow slightly according to the three-body scattering experiments (e.g.Mikkola & Valtonen 1992;Quinlan 1996;Sesana et al. 2006), which is also observed in our runs.According to Peters & Mathews (1963) and Peters (1964), the transition to the GW emission-driven phase and the efficiency of the GW emission depend sensitively on the binary eccentricity.Specifically, during the GW emission-driven phase, the decaying rates in the semimajor axis and the eccentricity can be approximately (at the 2.5PN level) described by (Peters 1964) A higher eccentricity value closer to 1 will significantly accelerate the GW emission phase due to the factor of (1−e 2 ) −7/2 .This is the dominant mechanism in driving the rapid SMBH merger in the DD-11-G5 NoCool run, i.e. the SMBH binary has an e hard ∼ 0.995 and the GW emission rapidly dominates the orbital decay over other processes, resulting in a very short merger time-scale after t hard (i.e.∼20 Myr).Similarly, the low e hard ∼ 0.5 is partly responsible for why the NoGas (CoolStarThmAGN) run of the DD-11-G5 (EE-11-G5) merger has a much longer merger times-scale compared to other runs with higher e hard .We note that in gas-free N -body simulations, e hard has been found to suffer from stochasticity 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).Specifically, pure N -body simulations of galaxy mergers starting from different realizations of initial conditions can lead to distinct e hard for their hard binaries.Previous works have suggested that the scatter of e hard due to this stochasticity can be reduced when the simulations are run with higher mass resolutions (or equivalently more particles; see e.g.Quinlan & Hernquist 1997;Berczik et al. 2005;Nasim et al. 2020).For the simulations including gas in this study, e hard can have higher degrees of stochasticity compared to the gas-free case, as the galaxy formation subgrid model introduces additional stochastic processes.Another uncertainty in the binary eccentricity comes from the torque interaction between the binary and its circumbinary disc.Recently, some small-scale hydrodynamical simulations (e.g.D 'Orazio & Duffell 2021;Zrake et al. 2021;Siwek et al. 2023) suggested that the binary-disc torque interaction might drive the eccentricity toward two 'attractor' solutions, i.e. an eccentric and equal-mass binary with e ≳ 0.1 (e ≲ 0.1) evolves to an equilibrium value of eeq ∼ 0.5 (eeq ∼ 0).
Although the detailed SMBH merger time-scales of the different simulations are affected by varying e hard , in Fig. 9, we can still see an overall trend that the SMBH binaries in disc galaxy mergers tend to merge more rapidly than those in elliptical mergers.Specifically, the SMBHs in the DD-11-G5 runs tend to cluster around a lower e hard ∼ 0.6 and merge within a few hundred Myr (i.e.≲ 500 Myr) after they become hard (except the NoGas run which has a low central stellar density and a low e hard ).In contrast, the SMBHs in the EE-11-G5 runs tend to cluster around a higher e hard ∼ 0.9 but merge at time-scales ≳ 1 Gyr after t hard (except the CoolStarNoAGN run which has the highest central stellar density among all elliptical merger runs).In Appendix A, we demonstrate that the contrast in SMBH merger time-scales between disc and elliptical galaxy mergers remains consistent across runs with different mass resolutions.Such a difference in SMBH merger time-scales aligns with the higher central stellar densities and more triaxial potentials observed in disc galaxy merger remnants compared to elliptical cases.
For a direct comparison of SMBH merger time-scales across runs including different physical processes, it is necessary to fix e hard to similar values to mitigate the significant impact from varying eccentricities.In Paper I, we rerun the simulations from t hard by resetting e hard to the same value, which allow us to perform a more direct and more quantitative comparison of the merger timescales.In Paper I, we demonstrate that the torque interaction between the binary and the circumbinary disc only plays a minimal role in affecting the shrinking of the orbit in the G5 galaxy merger.

CONCLUSIONS
In this work, we presented the RABBITS simulations and investigated how the presence and the specific implementation (i.e.thermal versus kinetic) of AGN feedback influence the properties of merging galaxies and the orbital decay of their SMBHs.We considered both disc-disc and elliptical-elliptical galaxy mergers, and performed systematically controlled simulation sets including different physical processes (i.e.NoGas, NoCool, CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN; see Table 3 for details).
Our key findings are summarized as follows.
AGN feedback plays a critical role in shaping the galaxy properties: (i) Simulations with cooling and stellar processes (i.e.star formation and stellar feedback), but without AGN feedback, exhibit excessive star formation in the galaxy centre, resulting in galaxy remnants which deviate from some observed galaxy scaling relations (e.g. the MBH-σ⋆ relation).
(ii) Compared to thermal AGN feedback, kinetic AGN feedback is more effective in removing gas from the galaxy centre, preventing cooling flows, quenching star formation, and maintaining a red and quiescent galaxy remnant.
Different central galaxy properties caused by different AGN models further affect the SMBH orbital decaying process: (iii) Galaxies with higher star formation and higher SMBH masses (attributed to increased accretion) tend to possess denser central regions which are more resistant to tidal stripping, experience larger dynamical friction, and consequently form bound SMBH binaries earlier.
(iv) In the dynamical friction phase, as AGN feedback effectively reduces the gas density around SMBHs, stellar dynamical friction dominates over the gaseous friction in shrinking the SMBH orbit.
(v) In the binary hardening phase, compared to ellipticalelliptical galaxy mergers, disc-disc mergers tend to have more rapid SMBH hardening and shorter SMBH merger time-scales (i.e.≲ 500 Myr for most disc mergers versus ≳ 1 Gyr for most elliptical mergers), as a result of the higher central stellar density (attributed to increased star formation) and the more triaxial potential.
(vi) The exact SMBH merger time-scale in a galaxy merger is sensitive to the eccentricity when the SMBHs become hard.Given the results above, which demonstrate the sensitivity of the SMBH coalescing process to the galaxy properties, we can conclude that to improve the modelling of SMBH merger time-scales in galaxy mergers, it is critical to thoroughly consider AGN feedback and its subgrid implementations.
In this study, we have explored the impact of pure thermal and pure kinetic AGN feedback.In future work, we plan to take this a step further by considering a two-mode or more sophisticated AGN feedback model (see discussions in e.g.Naab & Ostriker 2017).In addition, it will be beneficial to expand our exploration to a broader parameter space by considering galaxy mergers with other properties of the initial galaxies (e.g.virial mass, SMBH mass, gas fraction, etc.) and various orbital configurations.

APPENDIX A: RESOLUTION STUDY
To study the impact of the numerical resolution, we perform a resolution study with the kinetic AGN feedback model (i.e.CoolStarK-inAGN), which was recently introduced into the KETJU code by L23, for both the DD-11-G5 and the EE-11-G5 galaxy mergers.
The runs with the fiducial resolution, which are presented and analysed in the main text, are named 'Fiducial-Np' here.We further increase the resolution by a factor of 2 ('Hi-res-2Np') and decrease the resolution by a factor of 2 ('Low-res-Np/2') and 4 ('Low-res-Np/4'), to create simulations with different resolutions.The gravitational softening lengths are set according to ϵ = ϵ fid /f 1/3 , where ϵ fid is the softening length of the fiducial run and f ≡ Np/N p,fid is the particle number ratio between the resolution test run and the fiducial run.The detailed particle numbers, particle masses, and softening lengths of the different runs are summarized in Table A1.Only the softening lengths and accordingly the KETJU region radii are changed across the different resolution runs, all other simulation parameters remain identical to the fiducial ones.
The star formation and SMBH growth histories are plotted in Fig. A1.Overall, all the runs with different resolutions exhibit quite similar star formation and SMBH growth histories for both the DD-11-G5 and the EE-11-G5 mergers.The Low-res-Np/4 ones tend to have lower SFRs and higher SMBH accretion rates compared to the runs with higher resolutions.But the fiducial run shows results which are very similar to the Hi-res-2Np run, suggesting a relatively good resolution convergence.
The time-evolution of the orbital parameters are summarized in Fig. A2.As discussed in the main text, the exact SMBH merger time-scale of an individual simulation is strongly influenced by e hard .However, regardless of the numerical resolution, the disc galaxy mergers generally result in shorter SMBH merger timescales compared to elliptical galaxy mergers, which is attributed to the higher central stellar densities in the disc galaxy merger remnants.
To conclude, our runs with the kinetic AGN feedback exhibit a relatively good convergence across different resolution tests.

APPENDIX B: INFLUENCE OF AGN FEEDBACK ON ISOLATED ELLIPTICAL GALAXIES
It was showed in Eisenreich et al. (2017) that the details of the AGN feedback implementations have a significant influence on the evolution and properties of isolated elliptical galaxies.In this appendix, we perform similar simulations to study how our AGN feedback models affect the quiescence and the circumgalactic medium metal enrichment of isolated elliptical galaxies.
We adopt the E1 galaxy as the initial condition and evolve it for 3 Gyr using the CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN simulation setups.The simulation parameters used in these isolated galaxy runs are identical to those of the galaxy merger runs as detailed in Section 2.3.
The star formation and SMBH growth histories from different runs are plotted in Fig. B1.Without AGN feedback, the SFR settles to a value of ∼1 M ⊙ yr −1 , or equivalently the specific SFR (sSFR ≡ SFR/M⋆) is maintained at ∼10 −11 yr −1 , making the galaxy a green valley galaxy, which lies between the blue star-forming and red quiescent sequences (Salim 2014).Star formation occurs within the circumnuclear disc, which forms as gas cools and condenses from the rotating hot gas halo (i.e.cooling flows).In the CoolStarThmAGN run, the energy from the thermal AGN feedback heats the gas close to the SMBH, but the thermal energy quickly radiates away due to the efficient cooling in the circumnuclear disc.Consequently, thermal AGN feedback only affects the star formation in a small central region, and the total SFR as a function of time is quite close to that of the CoolStarNoAGN run, resulting in a green valley galaxy.However, in the CoolStarKinAGN run, we observe that for the majority of the simulation period, the total SFR remains below 0.1 M ⊙ yr −1 , which results in the galaxy being categorized as a red and quiescent galaxy.The strong pulsed kinetic AGN feedback effectively expels gas particles from the galaxy centre, causing collisions with gas particles at larger radii and thus heating the surrounding gas.This mechanism efficiently removes gas from the centre, preventing cooling flows and promoting a quiescent galaxy state.The SMBH accretion history demonstrates clear self-regulated behaviour, wherein an increased gas availability for SMBH accretion triggers stronger kinetic AGN feedback, Table A1.Details of CoolStarKinAGN simulations with different resolutions.From left to right, the simulation name, the star particle number, the gas particle number, the dark matter particle number, the star (or gas) particle mass, the dark matter particle mass, the Plummer equivalent star softening length, the gas softening length, and the dark matter softening length.From top to bottom rows, the total SFRs, the total SMBH accretion rates, and the total SMBH masses are plotted.The filled circles mark the SMBH merger events.Note that the SFRs and the SMBH accretion rates have been averaged over 10 Myr here.The simulations with different resolutions are distinguished using different colours.For both DD-11-G5 and EE-11-G5 mergers, the fiducial simulations (green), which are presented in the main text, show similar SFR and SMBH accretion histories as the high resolution runs (red), thus suggesting a relatively good resolution convergence.rapidly removing the excess gas and reducing the accretion rate.Consequently, the SMBH mass exhibits step-like growth patterns.

Simulations
The gas metallicity profiles of the evolved galaxies at t = 3 Gyr from the different runs are displayed in Fig. B2.Compared to the initial metallicity profile, the CoolStarNoAGN profile tends to be more concentrated as a result of gas cooling and the enrichment from the very central new stars.Once the thermal AGN feedback is included, the profile becomes slightly more extended than the CoolStarNoAGN one.On the other hand, the kinetic AGN feedback is much more effective in lowering the central gas metallicity and flattening the profile, which agrees better with the profiles of the observed isolated elliptical galaxies plotted in Fig. B2.
To summarize, compared to thermal AGN feedback, kinetic AGN feedback is more effective in driving gas out of the galaxy centre, quenching star formation, and flattening the gas metallicity profile.Overall, the results of our isolated elliptical galaxy tests are in line with the findings in Eisenreich et al. (2017).
This paper has been typeset from a T E X/L A T E X file prepared by the author.
Fig. 2, we use the CoolStarKinAGN runs of both the DD-11-G5 and EE-11-G5 mergers as examples to illustrate the SMBH trajectories and the galaxy interactions from our simulations.The left plot of Panel (a) and the plots in Panel (b) of Fig. 2 depict the evolution of the DD-11-G5 galaxy merger.Initially the two galaxies are set on parabolic orbits (i.e. the dashed lines in Panel (a)) with a separation of R200.Within the first ∼650 Myr (i.e.

Figure 2 .
Figure 2. SMBH trajectories and the evolution of the stellar and gaseous components in the CoolStarKinAGN simulations.(a) SMBH trajectories for both the DD-11-G5 (left) and EE-11-G5 (right) runs.The incoming parabolic orbits are plotted with dashed curves while the actual SMBH paths are shown with solid curves.The black circles mark the initial SMBH positions.The coloured circles mark the positions at six later snapshots, with smaller sizes representing later time.The six corresponding snapshots of the galaxy evolution are plotted in Panels b and c.(b) Evolution of the projected stellar (upper) and gaseous (lower) densities from the DD-11-G5 simulation.During the merger of disc galaxies, stellar/gaseous tails and bridges launch as a result of the tidal response.The colour bars for the projected densities are given at the top-right corner of the figure.(c) Similar to Panel b, but for the EE-11-G5 run.The same colour bars as in the DD-11-G5 case have been adopted here.Unlike the disc galaxy mergers, no narrow tidal tails or bridges are observed in the elliptical galaxy mergers.

Figure 3 .
Figure3.Star formation and SMBH accretion histories in the DD-11-G5 (left) and EE-11-G5 (right) mergers.From top to bottom, the time evolution of the total SFRs, the total SMBH accretion rates, and the total SMBH masses in the simulations are plotted.The CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are displayed with green, red, and purple lines, respectively.The filled circles mark the SMBH merger events in the different runs.Note that the SFRs and the SMBH accretion rates are averaged over 10 Myr here.The kinetic AGN feedback implementation is more effective in suppressing the star formation in galaxies and it results in more bursty SMBH accretion rates.

Figure 4 .
Figure 4. Gas phase diagrams at t = 2 Gyr.The DD-11-G5 and the EE-11-G5 runs are shown in the top and bottom panels, respectively.From left to right, the CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are plotted.The colour represents the number of gas particles in the different temperature-density bins, and the colour bar for each row is displayed on the right.The horizontal and the vertical solid lines mark the temperature and the density thresholds for star formation, respectively.The dashed line marks the Jeans mass resolution limit, M J,lim = N ngb mgas.For gas particles below this dashed line, the Jeans mass is not resolved and they can collapse under self-gravity due to numerical noise (seeEisenreich et al. 2017).Compared to the cases of no AGN feedback and thermal AGN feedback, the kinetic AGN feedback model results in distinct gas distributions in the phase diagram, i.e. less gas particles in the star-forming phase and more gas in the hot and diffuse phase.
Figure5.Comparison of galaxy properties in simulations and observations.In all panels, the initial isolated galaxies are plotted with cyan, while the merger remnants in the CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are shown with green, red, and purple, respectively.The DD-11-G5 and EE-11-G5 runs are plotted with filled circles and diamonds, respectively.(a) M BH -σ⋆ relation.The error bars in the line-of-sight stellar velocity dispersion show the standard deviations computed from 50 random projected orientations.The solid line plots the best-fitting M BH -σ⋆ relation from van den Bosch (2016) and the grey-shaded regions show the intrinsic one-sigma scatter.The dashed lines give the best-fitting relations fromKormendy &Ho (2013) andTremaine et al. (2002).(b) Galaxy size-stellar mass relation.Here the galaxy size is characterized using the projected half stellar mass radius Re.The solid lines and shaded regions show the observed relations for early-type galaxies (ETGs) and late-type galaxies (LTGs) fromShen et al. (2003) andNewman et al. (2012).(c) The luminosity of soft X-rays (0.3-5 keV) from hot and diffuse gas versus the galaxy stellar mass.The X-ray luminosity is computed following the method detailed in L23.Note that the initial disc galaxy is not shown in this panel because it only has a cold gas disc and its hot gas X-ray luminosity is thus zero.The observational data from the MASSIVE(Ma et al. 2014) and ATLAS3D(Cappellari et al. 2011) surveys are fromGoulding et al. (2016).(d) Stellar metallicity-stellar mass relation.The black line and the grey-shaded region show the median observed relation and the 16th/84th percentiles fromGallazzi et al. (2005).Here, we adopt Z ⊙ = 0.0127 for the solar metallicity.(e) Stellar mass-halo mass relation.The black curve shows the best-fitting relation fromKravtsov et al. (2018) and the grey-shaded regions plot the scatter in stellar mass (i.e.0.2 dex).The dashed line gives the best-fitting relation fromMoster et al. (2013).(f) Dark matter mass fraction within the projected effective radius as a function of galaxy stellar mass.The crosses plot the dark matter fractions of 149 MaNGA ETGs fromJin et al. (2020).Overall, the galaxy merger remnants in the simulations with AGN feedback included agree well with the observations.

Figure 6 .
Figure 6.Time evolution of the SMBH separation for both the disc (top) and the elliptical (bottom) galaxy mergers.The NoGas, NoCool, CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are plotted with blue, orange, green, red, and purple lines, respectively.We also plot the additional NoGas test runs which exclude the gas particles in the initial condition (grey) for comparison.The horizontal line segments in the left panels mark the influence radius R infl (solid) and the hard binary separation R hard (dashed).

Figure 7 .
Figure7.Enclosed mass profiles centring on one of the SMBHs in DD-11-G5 (left) and EE-11-G5 (right) runs at the first apocentre (t ≈ 1 Gyr).In this comparison, the SMBHs from different simulations are matched using their respective particle IDs.The upper panels show the enclosed total mass (including SMBH, stars, gas, and dark matter) profiles, whereas the lower panels plot the enclosed mass profiles for new stars (solid) and gas (dotted).Here the new stars are those formed after the start of the simulation.The colour convention used to denote the different simulations is given in the top-left panel.Note that the profiles centred on the other SMBH are very similar to the profiles plotted here.

Figure 8 .
Figure8.Spherically averaged density profiles centred on the CoM position of the two SMBHs for the galaxy merger remnants in the DD-11-G5 (left) and the EE-11-G5 (right) runs.The density profiles computed from all star particles, new star particles only, and all gas particles are plotted with solid, dotted, and dashed curves, respectively.Only bins with at least ten particles are shown.Different simulations are distinguished by colours according to the convention outlined in the top-left panel.The upper panels display the density profiles when the SMBH separation is roughly twice the SMBH influence radius (i.e. the SMBH orbital decay is dominated by the dynamical friction from surrounding stars and gas).The lower panels show the density profiles when the semimajor axis of the SMBH binary reaches approximately the hard separation (i.e. the onset of the SMBH binary hardening phase which is dominated by the three-body interactions between the SMBH binary and the individual stars).

Figure 9 .
Figure 9. Time-evolution of the PN-corrected orbital parameters, semimajor axis a and eccentricity e, in the DD-11-G5 (upper) and the EE-11-G5 (lower) runs.The filled circles mark the time when the SMBH binaries become hard.Different colours are adopted to distinguish different simulations as given in the upper semimajor axis panel.

Figure A1 .
Figure A1.Star formation and SMBH accretion histories from the CoolStarKinAGN runs with different resolutions.The DD-11-G5 and EE-11-G5 galaxy mergers are displayed in the left and right columns, respectively.From top to bottom rows, the total SFRs, the total SMBH accretion rates, and the total SMBH masses are plotted.The filled circles mark the SMBH merger events.Note that the SFRs and the SMBH accretion rates have been averaged over 10 Myr here.The simulations with different resolutions are distinguished using different colours.For both DD-11-G5 and EE-11-G5 mergers, the fiducial simulations (green), which are presented in the main text, show similar SFR and SMBH accretion histories as the high resolution runs (red), thus suggesting a relatively good resolution convergence.

Figure A2 .
Figure A2.Time-evolution of the SMBH binary orbital parameters from the CoolStarKinAGN runs with different resolutions.The DD-11-G5 and EE-11-G5 runs are shown in the upper and lower panels, respectively.From low to high resolutions, the runs are plotted with blue, orange, green, and red lines.The filled circles mark the time when the two SMBHs become hard in the different runs.Regardless of the numerical resolution, the disc galaxy mergers generally result in shorter SMBH merger time-scales compared to the elliptical galaxy mergers, with the exact time-scales being influenced by the eccentricities at t hard .

Figure B1 .Figure B2 .
Figure B1.Star formation and SMBH accretion histories in the isolated elliptical galaxy simulation.From top to bottom panels, the SFR, SMBH accretion rate, and SMBH mass as a function of simulation time are plotted.The CoolStarNoAGN, CoolStarThmAGN, and CoolStarKinAGN runs are shown with green, red, and purple, respectively.

Table 2 .
Initial conditions of galaxy mergers.
The fiducial simulations are the ones presented in the main text.