EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A231 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A231
Number of page(s) 15
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202453070
Published online 24 March 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Massive black holes (MBHs) are found at the centre of most massive galaxies (Kormendy & Richstone 1995; Magorrian et al. 1998). MBHs co-evolve over cosmic history with their host galaxies as inferred observationally from the mass of the MBH (M•), found to correlate with properties of its host galaxy, such as the mass and stellar velocity dispersion of the central bulge (Magorrian et al. 1998; Tremaine et al. 2002; Marconi & Hunt 2003; Häring & Rix 2004) and, with larger scatter, with the stellar mass of the galaxy (M*) (Reines & Volonteri 2015; Greene et al. 2020). Both the MBH and the stellar mass can grow efficiently from the dense galactic gas (e.g. Di Matteo et al. 2005; Springel et al. 2005; Dubois et al. 2014c), which is replenished by cosmic inflows and cooling (e.g. White & Rees 1978; Birnboim & Dekel 2003; Kereš et al. 2005) and disrupted by the injection of energy from stars and MBH accretion themselves (e.g. Dubois et al. 2014b, 2015; Habouzit et al. 2017; Anglés-Alcázar et al. 2017).

Mergers are another key ingredient in the co-evolution of MBH and their galactic hosts. The tidal forces from galaxy mergers can trigger a phase of rapid star formation and MBH accretion (Hernquist 1989; Barnes & Hernquist 1991; McAlpine et al. 2018; Lapiner et al. 2021). More importantly, a merger of two galaxies can lead to the merger of the two central MBHs, thus aggregating the masses of the two parent galaxies and MBHs in a correlated manner (Peng 2007; Kormendy & Ho 2013). The journey from a galaxy merger to an MBH merger is a complicated one. After two galaxies merge, the two central MBHs will dissipate energy to dynamical friction and sink to the centre of the new potential. If they form an MBH binary, the binary can continue to shrink due either to stellar hardening (star sling-shots extracting the orbital energy of the binary) or viscous gas drag from the circum-binary disc until, if the evolution is fast enough, dissipation via gravitational waves (GWs) becomes efficient and the MBHs coalesce (Begelman et al. 1980). The future Laser Interferometer Space Antenna (LISA) mission aims to measure the GWs in the final stages of MBH mergers (Amaro-Seoane et al. 2023).

In the last stage of coalescence, any asymmetries in the spins or masses of the binary lead to the final burst of GWs being emitted anisotropically and carrying away linear momentum. As a response, the remnant receives a velocity kick in the opposite direction (Bekenstein 1973). The magnitude of this kick can reach values of up to 5000 km s−1 for MBHs with high spins partly aligned in the direction normal to the orbital plane but pointing in opposite directions in the orbital plane (Lousto & Zlochower 2011). Even with other configurations, the velocity can reach several hundreds of km s−1. This means that the recoil velocity can exceed the escape velocity of haloes and displace MBH merger remnants from the nucleus or eject them from their hosts.

GW recoil can have a significant impact on the cosmic evolution of MBHs. The kicked MBH can see its growth severely reduced if it leaves the nuclear gas reservoir and has a high velocity relative to the interstellar medium (ISM) until it dissipates its orbital energy and returns to the nucleus (Blecha & Loeb 2008; Sijacki et al. 2011; Blecha et al. 2011, 2016). If accretion is reduced or the MBH is ejected from the nucleus, energy injection into the nucleus is also reduced, which can assist nuclear star formation and deplete gas available for future MBH growth (Blecha et al. 2011). For larger kicks that exceed the escape velocity of the host, the MBH can be ejected and the galaxy will be devoid of a central MBH until another MBH from another seeding event or galaxy merger is brought to the nucleus (Schnittman 2007; Volonteri et al. 2010; Dunn et al. 2020). The MBHs that escape the centre can become ‘wandering’ MBHs in the outskirts of the galaxy or in the galactic halo (Volonteri & Perna 2005; Micic et al. 2006, 2011; Untzaga et al. 2024)

The reduction of MBH growth and the ejection of MBHs due to GW recoil can leave observational imprints. Ejected MBHs can reduce the MBH occupation fraction of galaxies, the MBH merger rate, and decrease the normalisation and increase the scatter of the MBH–host mass correlations (Schnittman 2007; Volonteri 2007; Volonteri et al. 2010, 2011; Blecha et al. 2011; Dunn et al. 2020; Mannerkoski et al. 2022). GW recoil also decreases the MBH merger rate (e.g. Micic et al. 2006; Dunn et al. 2020). GW recoil can also leave a trace in the galactic morphology, acting in combination with MBH binary scouring to produce large stellar cores (Gualandris & Merritt 2008; Nasim et al. 2021; Rawlings et al. 2025). In some cases, the recoiling MBH could leave a trail of star-forming material (Ogiya & Nagai 2024) or it could even be directly detected (Komossa 2012).

Studying the cosmic effect of GW recoil is a difficult task. First, the magnitude of the kick, which regulates its effect on the future evolution of the MBH, is dependent on the mass ratio, spins, and orbital angular momentum of the MBH pair. Extreme mass ratios or equal masses and aligned configurations of the angular momenta yield small kicks, while some degree of misalignment leads to larger kicks (e.g. Campanelli et al. 2007; Lousto et al. 2012). These quantities depend complexly on the previous accretion, mergers and dynamics of the MBH pair (Dong-Páez et al. 2023b). Previous studies have found that high-density gas with coherent angular momentum aligns the orbital angular momentum and the spins leading to small recoil kicks, while in mergers with little gas content the directions are close to random and the kicks can have a larger effect (Bogdanović et al. 2007; Dotti et al. 2010; Volonteri et al. 2010; Blecha et al. 2016; Dunn et al. 2020). Moreover, the fate of a recoiling MBH depends on its dynamics in an often asymmetric and time-varying host potential in an expanding cosmology (Sijacki et al. 2009, 2011; Blecha et al. 2011; Choksi et al. 2017).

It follows that, to capture the cosmic effect of GW recoil and the complex co-dependence between GW recoil and MBH evolution and dynamics, one needs to run cosmological simulations that can capture the complexity of MBH environments and include realistic sub-grid models for accretion, spin evolution, and dynamics. Previous studies have used semi-analytic models relying on simplified prescriptions and often with unresolved dynamics, idealised simulations, or cosmological simulations that do not include all of these requirements. Often recoil kicks are estimated by assuming some distribution of spin magnitudes and directions at merger instead of calculating them from MBH spins computed self-consistently across the MBH accretion history. So far the only cosmological hydrodynamical simulations modelling MBH evolution with GW recoil kicks on-the-fly are those of Sijacki et al. (2009) and Mannerkoski et al. (2022). Sijacki et al. (2009) used a treatment of MBH spin evolution accounting for spin changes by mergers but neglected the role of MBH accretion for spin evolution, which is the dominant source of spin evolution for most grown MBHs (Dubois et al. 2014b; Bustamante & Springel 2019; Peirani et al. 2024; Sala et al. 2024).

In this paper, we present the first hydrodynamical cosmological simulation including on-the-fly GW recoil kicks calculated from self-consistently evolved MBH masses, spins (evolved from both accretion and mergers), and dynamics. The masses and spins are evolved due to gas accretion and MBH mergers. The simulations also include sub-grid models for the unresolved dynamical friction onto the MBHs. In this paper, we present the first hydrodynamical cosmological simulation including on-the-fly GW recoil kicks calculated from MBH masses, spins, and dynamics all directly evolved in the simulation, increasing the self-consistency with respect to previous investigations. The masses and spins are evolved according to gas accretion and MBH mergers. The simulations also include sub-grid models for the unresolved dynamical friction onto the MBHs. We use this simulation to study the interplay between GW recoil and MBH evolution and the dynamics of recoiling MBHs in a cosmological environment. In Sect. 2 we describe our simulations and our analysis method. In Sect. 3, we study the evolution of MBH in our simulations. In Sect. 4, we study the dynamics of recoiling MBHs and compare recoiling MBHs to the global population of wandering MBHs. We conclude in Sect. 5.

2. Method

2.1. The OBELISK-RECOIL and OBELISK-NORECOIL simulations

We have run two cosmological zoom-in hydrodynamical simulations, OBELISK-RECOIL and OBELISK-NORECOIL, using the adaptive mesh refinement code RAMSES (Teyssier 2002) down to z ∼ 4.4. Our simulations are similar to the OBELISK simulation (Trebitsch et al. 2021), the main difference being that OBELISK-RECOIL includes a sub-grid model for GW recoil kicks and both simulations add adjustments to other sub-grid models in OBELISK. We used these two simulations and OBELISK to cover a larger parameter space. Below we present a summary of the OBELISK-RECOIL and OBELISK-NORECOIL simulations, highlighting especially the differences with OBELISK. For a more detailed description of the simulation model, we refer the reader to Trebitsch et al. (2021).

2.1.1. Initial conditions

The initial conditions are based on the HORIZON-AGN simulation (Dubois et al. 2014a). We selected particles enclosed within a sphere centred on the most massive halo in HORIZON-AGN at z = 1.97 with radius 2.51 h−1 cMpc. We found the convex hull that encloses these particles in the initial conditions. This region was simulated with a succession of initial nested coarse grids down to maximum initial level of refinement ℓ = 12 corresponding to a high-resolution dark matter (DM) particle mass of 1.2 × 106 M⊙. The rest of the 100 h−1 cMpc box was simulated at lower resolution with a minimum level of refinement ℓ = 8.

We assumed a standard ΛCDM cosmology. The cosmological parameters were taken to be the best-fit parameters from the WMAP-7 analysis (Komatsu et al. 2011) – Hubble constant H0 = 70.4  km s−1 Mpc−1, dark energy density parameter ΩΛ = 0.728, total matter density parameter Ωm = 0.272, baryon density parameter Ωb = 0.0455, amplitude of the power spectrum σ8 = 0.81, and spectral index ns = 0.967.

2.1.2. Gravity and hydrodynamics

The gas was evolved using an unsplit second-order MUSCL–Hancock scheme (van Leer 1979). The inter-cell conservative variables were reconstructed using a total variation diminishing scheme with minmod slope limiter for the two states of the Riemann problem solved with the approximate Harten-Lax-van Leer-Contact Riemann solver. We assumed the gas is ideal and monoatomic with an adiabatic index of 5/3. The time step follows the courant condition with a two-fold adaptive time stepping between levels larger than ℓ = 13 and single time stepping between coarser levels.

Any cell that exceeds a mass of 8 times the initial mass resolution is refined down to a minimum cell size of Δx ≃ 35 pc (proper, hence, corresponding to level ℓ = 19 and 20 when passing redshifts 9 and 4 respectively) within the zoomed-in region. DM, stars, and MBHs were modelled as collisionless particles, using a cloud-in-cell interpolation which size corresponded to that of the underlying cell (down to 35 proper pc), except for DM particles that have a minimum cloud size of 540 comoving pc (i.e. level ℓ = 18 or 100 proper pc at z = 4.41) in order to reduce the shot noise effect of the mass resolution imbalance between star and DM particles. Gravitational acceleration is obtained solving the Poisson equation with a multigrid solver (Guillet & Teyssier 2011) on levels coarser than ℓ < 13 and a conjugate gradient otherwise. We further enforced the mesh to be refined to the minimum cell size of Δx ≃ 35 pc around MBHs to better capture their dynamics (Lupi et al. 2015).

2.1.3. Sub-grid modelling for gas and stars

Unlike OBELISK, our simulations do not include radiation. Instead of using a non-equilibrium thermo-chemistry model for hydrogen and helium, we assume collisional-ionisation equilibrium in a homogeneous UV background (Haardt & Madau 1996) below a redshift of reionisation of z = 8.5, similarly to the NEWHORIZON simulation (Dubois et al. 2021). Metal-enriched gas can cool further following the cooling functions from Dalgarno & McCray (1972) below a temperature of T < 104 K, and Sutherland & Dopita (1993) for T ≥ 104 K. We also account for the self-shielding against UV radiation of high-density gas (Rosdahl & Blaizot 2012).

A cell was considered to be star-forming if the density of gas exceeds 5 H cm−3 and the turbulent Mach number ℳ ≥ 2. Star formation follows a Schmidt (1959) law inversely proportional to the gas free-fall time with an efficiency parameter that depends on ℳ and αvir = 2Ek/Eg (where Ek and Eg are respectively the turbulent kinetic energy and the gravitational energy of the gas) following the multi-free fall model (Hennebelle & Chabrier 2011; Federrath & Klessen 2012) of Padoan & Nordlund (2011). Stars were modelled with 104 M⊙ particles that represent stellar populations with a Kroupa (2001) initial mass function. 5 Myr after birth, 20% of the stellar mass leads to supernova (SN) explosions, releasing metals with yield y = 0.075 and an energy of 1051 erg per 10 M⊙ individual SN into the gas, following the numerical scheme from Kimm & Cen (2014).

2.1.4. Sub-grid modelling for MBHs

We seeded MBHs with a mass of Mseed = 105 M⊙ in cells where both the stellar and gas densities exceed a threshold of ngas = 200 H cm−3. We note that our seeding model has a more stringent criterion and produces more massive seeds than in OBELISK, where Mseed = 3 × 104 M⊙ and ngas = 100 H cm−3. To avoid the spurious formation of multiple seeds in one galaxy, we only allowed a cell to form a MBH if there are no other MBHs within a radius of 50 ckpc.

The unresolved dynamical friction from collisionless particles and gas was modelled using the implementation by Dubois et al. (2013) and Pfister et al. (2019). The frictional force from gas follows the analytical expression from Ostriker (1999) but boosted by a factor (ρ/ρDF, th) if ρ > ρDF, th = 10 H cm−3. We note that the normalisation of the dynamical friction force from the gas is over a factor of 106 smaller than in OBELISK, arguably leading to more realistic MBH dynamics.

The MBHs in the simulation grow by gas accretion or MBH mergers. Gas accretion was modelled as spherical Bondi–Hoyle–Lyttleton accretion (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944),

M ˙ acc = 4 π G 2 M • 2 B ¯ . $$ \begin{aligned} \dot{M}_\mathrm{acc} =4\pi G^2M_\bullet ^2 \bar{B}. \end{aligned} $$(1)

The variable B is defined as B = ρ/(cs2 + vrel2)3/2, where ρ is the gas density, cs is the sound speed, and vrel is the gas velocity relative to the MBH. The bar denotes the average in a Gaussian kernel around the MBH with a characteristic scale comparable to the BHL radius RBHL = GM•/(cs2 + vrel2). We note that, compared to OBELISK, the average is computed on the composite quantity B instead of the individual quantities ρ, cs, and vrel.

The accretion rate was limited by the Eddington rate,

M ˙ Edd = 4 π G M • m p σ T ε r c ∼ 2 M ⊙ yr − 1 ( ε r 0.1 ) − 1 ( M • 10 8 M ⊙ ) , $$ \begin{aligned} \dot{M}_\mathrm{Edd} =\frac{4\pi GM_\bullet m_\mathrm{p} }{\sigma _\mathrm{T} \varepsilon _\mathrm{r} c}\sim 2\,M_\odot \,{yr}^{-1}\left(\frac{\varepsilon _\mathrm{r} }{0.1}\right)^{-1}\left(\frac{M_\bullet }{10^8\mathrm M_\odot }\right), \end{aligned} $$(2)

where mp is the proton mass, and σT is the Thompson cross-section, and εr is the radiated efficiency. The accretion flows was assumed to be radiatively efficient if fEdd = Ṁacc/ṀEdd ≥ fEdd,crit = 0.01, otherwise they are assumed to be fed by radiatively inefficient accretion flows (RIAFs). In the radiatively efficient regime, a fraction εr of the accreted mass is radiated and the remaining fraction is added to the MBH mass. The radiative efficiency εr depends on the MBH spin, ranging from 0.038 for counter-rotating material to 0.32 for co-rotating material around maximally spinning MBHs (a• = ±0.998). In the radiatively inefficient regime, the radiative efficiency was further decreased by a factor of fEdd/fEdd, crit (Benson & Babul 2009).

Every coarse time step of the simulation, two MBHs were allowed to merge if their separation became shorter than 4Δx. Any losses to GWs were neglected, and the mass of the remnant was taken to be the sum of the two parent MBHs.

The spins of MBHs were evolved with gas accretion and MBH mergers following Dubois et al. (2021) (similar to Dubois et al. 2014b; c) but with a modified model at fEdd < fEdd, crit). In the radiatively efficient regime the magnitude of the spin following an accretion episode was calculated as follows (Bardeen 1970),

a• = 1 3 r ∼ ( 4 − 3 r ∼ − 2 ) , $$ \begin{aligned} a_\bullet =\frac{1}{3}\sqrt{\tilde{r}}\left(4-\sqrt{3\tilde{r}-2}\right), \end{aligned} $$(3)

where r ∼ = r ISCO [ M• / ( M• + Δ M• ) ] 2 $ \tilde r=r_{\mathrm{ISCO}}[M_\bullet/(M_\bullet+\Delta M_\bullet)]^2 $. The radius of the innermost stable circular orbit in units of the gravitational radius, rISCO, is a function of the spin only. The spin magnitude evolved by accretion was limited to −0.998 < a• < 0.998 (Thorne 1974). The angular momentum of the gas in the outer accretion disc was assumed to align with the gas on the resolved scales. The inner disc aligned or anti-aligned with the MBH spin through the Bardeen & Petterson (1975) effect depending on the criterion of King et al. (2005) for misalignement. The direction of the MBH spin following an accretion episode was calculated by summing the angular momenta of the accreted gas and the MBH. Assuming the MBH accretes from a reservoir of gas with coherent angular momentum, the MBH spin aligns with the gas in a short timescale of approximately (Perego et al. 2009)

t al ∼ 0.1 a ∙ 5 / 7 ( M• 10 6 M ⊙ ) − 2 / 35 f Edd − 32 / 35 Myr . $$ \begin{aligned} t_\mathrm{al} \sim 0.1 a_\bullet ^{5/7} \left(\frac{M_\bullet }{10^6\mathrm M_\odot }\right)^{-2/35} f_\mathrm{Edd} ^{-32/35}\,\mathrm{Myr} . \end{aligned} $$(4)

In the RIAFs regime, the jets were assumed to be powered from the rotational energy of the MBH, and the spin decreased following the simulations of McKinney et al. (2012) with the fitting function provided in Dubois et al. (2021).

Following a MBH merger the spin was updated as (Rezzolla et al. 2008)

a f = 1 ( 1 + q ) 2 ( a 1 + a 2 q 2 + ℓ q ) , $$ \begin{aligned} \boldsymbol{a}_f = \frac{1}{(1+q)^2}\left(\boldsymbol{a}_1 +\boldsymbol{a}_2 q^2 +\boldsymbol{\ell } q \right), \end{aligned} $$(5)

where a1 and a2 are the spins of the primary (the most massive) and secondary (least massive) MBHs. The mass ratio is the ratio of the secondary to the primary mass, q = M2/M1, and ℓ is related to the orbital angular momentum, and calculated as

ℓ = s 4 ( 1 + q 2 ) 2 ( a 1 2 + a 2 2 q 4 + 2 a 1 . a 2 q 2 ) + s 5 η + t 0 + 2 1 + q 2 ( a 1 cos ϕ 1 + a 2 q 2 cos ϕ 2 ) + 2 3 + t 2 η + t 3 η 2 , $$ \begin{aligned} \ell =&\frac{s_4}{\left(1+q^2\right)^2}\left(a_1^2 + a_2^2 q^4 + 2 \boldsymbol{a}_1 . \boldsymbol{a}_2 q^2\right)\nonumber \\&+ \frac{s_5\eta + t_0 + 2}{1+q^2}\left(a_1 \cos \phi _1 + a_2 q^2 \cos \phi _2\right) \\ &+ 2\sqrt{3} + t_2\eta + t_3\eta ^2, \end{aligned} $$(6)

where ϕ1 and ϕ2 are the angle of a1 and a2 compared to ℓ, and the symmetric mass ratio is η = q/(1 + q)2. The constants are set to s4 = −0.129, s5 = −0.384, t0 = −2.686, t2 = −3.454, and t3 = 2.353. ℓ is taken to be aligned to the orbital angular momentum of the binary in the time step before the merger.

Active galactic nucleus (AGN) feedback was modelled using a dual-mode model. In the RIAF regime, the MBH injected kinetic energy in two cylindrical regions in the direction of its spin. This represents two relativistic jets of a radio-mode AGN. The fraction of energy into the jet is a function of the MBH spin and was calculated from a fit to the simulations of magnetically chocked accretion flows from McKinney et al. (2012). In the radiatively efficient regime, the MBH injected as thermal energy 5% of the accreted energy isotropically, similar to a radiative-mode AGN. We note that the feedback efficiency has been decreased compared to OBELISK, which assumes a fraction of 15%. This choice is supported by results from a similar high-resolution simulation, ‘NewCluster’, a galaxy cluster simulation with 70 pc resolution, that better matches the MBH-to-stellar mass relation at low redshift (Han et al., in prep.).

To assess our uncertainty in the accretion of MBHs, we ran two additional simulations where the accretion rate was enhanced by setting vrel = 0 in Eq. (1), one simulation with recoil and one without recoil. These simulations were run down to redshift 6.

2.1.5. Sub-grid GW recoil kicks following MBH mergers

As described above, MBHs were merged if their separation decreases below 4Δx. The remnant MBH was kicked in the coarse time step immediately after an MBH merger, which is an important simplification of the true MBH dynamics (e.g. Milosavljević & Merritt 2001; Khan et al. 2011; Rantala et al. 2017). Indeed, the final coalescence can take several hundred and up to several billions of years after they numerically connect at a scale of a few 100 pc (Katz et al. 2020; Volonteri et al. 2020; Chen et al. 2022; Li et al. 2022, 2024; Dong-Páez et al. 2023b). However, we defer to future work the self-consistent treatment of MBH ‘sub-grid’ dynamics (using e.g. the RAMCOAL model of Li et al. 2024) that leads to the final MBH binary coalescence and its GW recoil kick in a more consistent way. We modelled the magnitude and direction of the recoil velocity using a fitting formula to the results of numerical relativity simulations, as quoted by Lousto et al. (2012):

v recoil ( q , a 1 , a 2 , L ̂ orb ) = v m e 1 + v ⊥ ( cos ξ e 1 + sin ξ e 2 ) + v ‖ L ̂ orb , $$ \begin{aligned} \boldsymbol{v}_\mathrm{recoil} (q,\boldsymbol{a}_1, \boldsymbol{a}_2, \hat{\boldsymbol{L}}_\mathrm{orb} ) = v_m\boldsymbol{e}_1+v_\bot (\cos \xi \boldsymbol{e}_1 +\sin \xi \boldsymbol{e}_2) + v_\parallel \hat{\boldsymbol{L}}_\mathrm{orb} \, , \end{aligned} $$(7)

with

v m = A m η 2 ( 1 − q ) 1 + q [ 1 + B m η ] , v ⊥ = H η 2 1 + q ( a 1 , ‖ − q a 2 , ‖ ) , v ‖ = 16 η 2 1 + q [ V 1 , 1 + V A S ∼ ‖ + V B S ∼ ‖ 2 + V C S ∼ ‖ 3 ] | a 1 , ⊥ − q a 2 , ⊥ | × × cos ϕ , $$ \begin{aligned} v_m&=A_m\frac{\eta ^2(1-q)}{1+q}\left[1+B_m\eta \right],\nonumber \\ v_\perp&=H\frac{\eta ^2}{1+q}\left(a_{1,\parallel }-qa_{2,\parallel }\right) ,\\ v_\parallel&=16\frac{\eta ^2}{1+q}\left[V_{1,1}+V_A\tilde{S}_\parallel +V_B\tilde{S}_\parallel ^2+V_C\tilde{S}_\parallel ^3\right]\left|\boldsymbol{a}_{1,\perp }-q\boldsymbol{a}_{2,\perp }\right|\times \nonumber \\&\times \cos \phi ,\end{aligned} $$(8)

where e1 is the unit vector from the secondary to the primary, L ̂ orb $ \hat{\boldsymbol{L}}_{\mathrm{orb}} $ is the unit orbital angular momentum vector, e 2 = L ̂ orb × e 1 $ \boldsymbol{e}_2=\hat{\boldsymbol{L}}_{\mathrm{orb}}\times\boldsymbol{e}_1 $, η = q/(1 + q)2 is the symmetric mass ratio, and S ∼ $ \tilde{\boldsymbol{S}} $ is defined as S ∼ = 2 ( a 1 + q 2 a 2 ) / ( 1 + q ) 2 $ \tilde{\boldsymbol{S}}=2(\boldsymbol{a}_1+q^2\boldsymbol{a}_2)/(1+q)^2 $. The subscript ∥ and ⊥ refer to the parallel and perpendicular components relative to the orbital angular momentum. The constants were set to ξ = 145°, H = 6.9 × 103 km s−1, Am = 1.2 × 104 km s−1, Bm = −0.93, V1,1 = 3677.76 km s−1, VA = 2481.21 km s−1, VB = 1792.45 km s−1, VC = 1506.52 km s−1 (Lousto & Zlochower 2008, 2009; Zlochower et al. 2011; Lousto et al. 2012). ϕ is a phase we assumed to be random following Lousto & Zlochower (2011). Recoil velocity is zero for equal-mass and equal-spin binaries, up to ∼200 km s−1 for non-spinning black holes, and up to ∼5000 km s−1 for spinning black holes in particular configurations.

The spins, masses, and orbital angular momentum of the primary and the secondary were measured at the time-step before the numerical merger and the resulting recoil velocity was injected into the merger remnant at the time-step after the MBH merger.

2.2. MBH-galaxy assignment

Since MBHs were formed based on local – cell level – properties and they were free to move in the simulation, they have to be assigned to host galaxies. This assignment was performed on the catalogue of galaxies resulting from the simulation. Galaxies and DM haloes were identified using the same method presented in Trebitsch et al. (2021). We used a version of ADAPTAHOP (Aubert et al. 2004; Tweed et al. 2009) that runs on the total mass distribution of stars and DM. We only considered uncontaminated (with no low-resolution particles) galaxies and haloes with more than 100 star and DM particles. The centre of the galaxy was defined to be the peak of the stellar density distribution using a shrinking sphere approach (Dubois et al. 2021).

To further analyse the catalogue of simulated MBH and their host galaxies, we assigned MBHs from the simulation to galaxies identified using the structure-finding algorithm in post-processing using a similar approach to Dong-Páez et al. (2023a,b). We first assigned ‘main’ MBHs, defined as the most massive MBH inside the half-mass radius R50, to galaxies in descending order of stellar mass M*. Unassigned MBHs are then assigned as ‘satellite’ MBHs to the most massive galaxies that enclose them within R90. R90 is the radius encompassing 90% of the stellar mass. Finally, any unassigned MBHs was assigned as ‘halo’ MBHs to the most massive galaxy that encloses them within the virial radius of their DM halo (Rvir).

We also assigned MBH mergers to galaxies. Since recoil kicks can eject MBHs on timescales shorter than the typical time between two snapshots (15 Myr), we assigned MBH mergers using the information from the snapshot before the merger. The host was taken to be the most massive galaxy hosting either of the MBHs within r50. Mergers involving two MBHs located outside of R50 were considered to be spurious and discarded from the sample. The meaning of “spurious” here is that we considered that MBHs that merged numerically far from the galaxy centre would not actually merge in reality, when accounting for further dynamical evolution: on the one hand if they were to form a binary off-centre, stellar hardening would not be effective, as it can shrink binaries only at high stellar densities, typical of galaxy nuclei. On the other hand, if they were to first sink to the centre and then form a binary, their dynamical friction timescale would be long. We consider MBHs within R50 as “central MBHs”, as in Dong-Páez et al. (2023b). Nearly 30% of mergers are discarded.

2.3. Analysis of recoiling and wandering MBH orbits

The orbit of a recoiling MBH was followed after the kick. The position of the MBHs was recorded at every coarse time-step in the simulation (approximately every ∼0.1 Myr), while the position of each galaxy was obtained from the snapshots, recorded every ∼15 Myr. We obtained the MBH position relative to its host by interpolating the galaxy position between snapshots.

An MBH was considered to have returned to the centre of the galaxy after a recoil kick when there was an apocentre of its orbit that lied inside R50. An MBH was considered to have escaped the halo if it escaped Rvir. An MBH was considered to be escaping the halo if it was outside R90 and its velocity was larger than four times the circular velocity at Rvir. The factor of four was chosen by inspection of trajectories. An MBH was considered to be wandering if it remained inside the halo but had not returned to the centre and its velocity did not exceed four times the circular velocity. See Table 1 for a summary of this classification.

Table 1.

Dynamical classification for recoiling MBH that we used in this study.

The circularity of an orbit is defined as

λ = L• L circ ( E ∙ ) . $$ \begin{aligned} \lambda =\frac{L_\bullet }{L_\mathrm{circ} (E_\bullet )}. \end{aligned} $$(9)

We have calculated the escape velocity of an MBH from its host halo as v esc = − 2 [ Φ ( x • ) − Φ ( R vir , rec ) ] $ v_{\mathrm{esc}}=\sqrt{-2[\Phi(\boldsymbol{x}_\bullet)-\Phi(\boldsymbol{R}_{\mathrm{vir,rec}})]} $. Rvir, rec is the point in the surface of the virial radius in the direction of vrecoil. This is to make sure we account for any inhomogeneities on the Rvir surface. Φ(x) is obtained directly from the gravity solver of the simulation. For an MBH merger, vesc was calculated from the snapshot before the merger. Blecha et al. (2011) show that the escape velocity can vary rapidly with time after a galaxy merger. We tested for any time variation in vesc by computing as an upper limit the value of vesc from the centre of the galaxy after the merger. We find that the variation between our estimate and this upper limit is small – most differences are below 10% and the median difference is of 5%. This is probably because galaxy mergers in the simulations are more frequently minor mergers than major mergers, and therefore the potential is not dramatically deepened by the collision.

3. GW recoil and the cosmic evolution of MBH

Throughout this study, we used the population of mergers from the OBELISK-RECOIL and OBELISK-NORECOIL. In addition to this, we also considered the OBELISK simulation. The OBELISK simulation included a boosted sub-grid dynamical friction, that couples MBHs dynamically to the dense gas and increases the efficiency of accretion. We used OBELISK as a limiting case for efficient dynamical friction and accretion.

Figure 1 shows the MBH mass distribution in OBELISK-RECOIL and OBELISK-NORECOIL. The two simulations have similar mass distributions that only differ in the high-mass tail, where OBELISK-RECOIL has a dearth of MBHs. In the following section, we study the interplay between GW recoils sourced by MBH mergers and MBH evolution.

thumbnail Fig. 1.

Mass of the main MBHs in the OBELISK-RECOIL (red), OBELISK-NORECOIL (blue) and OBELISK simulations at z = 4.41.

3.1. Recoil and MBH growth

GW recoil can decrease the growth of MBH by ejecting and removing MBH mass from galaxies and by hindering gas accretion. In Fig. 2, we show the fraction of the mass gained through mergers fmergers as a fraction of the total mass gain, M• − Mseed for the three simulations analysed here, OBELISK (in green), OBELISK-RECOIL (red), and OBELISK-NORECOIL (blue), where Mseed is 3 × 104 M⊙ for the former and 105 M⊙ for the latter two. In OBELISK-RECOIL and OBELISK-NORECOIL, MBHs grow mostly through MBH-MBH mergers, while accretion is subdominant. In OBELISK, a significant fraction of the MBHs that have gained more than 105 M⊙ settle in a phase of sustained efficient accretion. OBELISK featured an ad hoc boost in the gas dynamical friction, resulting in MBHs being much more dynamically coupled to the dense and cold ISM. As a result, MBHs are better centred in OBELISK-RECOIL and its twin than in OBELISK, but they have a larger velocity relative to the gas, which suppresses the accretion rate.

thumbnail Fig. 2.

Fraction of mass gained through mergers against the total mass gained above the seed mass for OBELISK-RECOIL (red), OBELISK-NORECOIL (blue), and OBELISK (green).

Recoil reduces the growth of MBHs both from accretion and mergers. The kick increases the relative velocity with the ISM, reducing the cross-section of the MBH for accretion. For stronger kicks, the MBH escapes the nucleus of the galaxy to a lower-density region, with less gas available for accretion. This effect is shown in Fig. 3. The average accretion rate of remnant MBHs over the 5 Myr following the mergers tends to experience a dramatic decrease after a recoil kick. The accretion rates 5 Myr following the merger in OBELISK-RECOIL are significantly smaller than in OBELISK-NORECOIL. The MBH growth is suppressed until the MBH dissipates its kinetic energy and settles in the centre – we discuss the timescales for an MBH to return to the centre of their host in Sect. 4.1.

thumbnail Fig. 3.

Distribution of the accretion rate normalised by the Eddington limit for OBELISK-RECOIL and OBELISK-NORECOIL, averaged over the 5 Myr after the merger for the remnants. The first bin shows the number of MBHs below the lower limit of the plot.

The integrated impact of the reduction in gas accretion over the cosmic evolution of MBHs is also seen in Fig. 2. In OBELISK-NORECOIL there are a few large MBHs (M∙ − Mseed > 105 M⊙) that have gained most of their mass through gas accretion, fmergers < 0.5. In OBELISK-RECOIL there are no MBHs in this regime. The growth from mergers is also reduced. There are fewer MBH mergers in OBELISK-RECOIL (92 mergers) compared to OBELISK-NORECOIL (104 mergers): this is because recoil can remove MBHs from the centre of galaxies, decreasing the MBH occupation probability of galaxies. In total, there are 13 fewer mergers, a difference of ∼10%. We note that the number of mergers for the simulation with recoil is possibly overestimated since in our simulations new MBHs can be seeded if the original MBH has been ejected from a galaxy and there are no other satellite MBHs.

The total effect of recoil on the growth of the global MBH population is explored in Fig. 1. The difference in MBH mass between the two simulations is noticeable. The OBELISK-RECOIL lacks the high mass tail present in the OBELISK-NORECOIL. The distributions in the host galaxy mass and fEdd (not shown) are similar for the two simulations. Recoil does not have a significant effect on the growth of the global galaxy population. This is partly because the MBHs do not grow efficiently in our simulation and they do not release significant energy into the ISM.

Recoils also decrease the typical total mass of MBH mergers. The total mass of MBH mergers tends to be lower for OBELISK-RECOIL. This is because more massive MBHs tend to experience more mergers and recoil kicks. This means that their growth via accretion is suppressed, resulting in smaller masses, but also that they are preferentially ejected from galaxies and therefore missing from subsequent mergers.

Recoils produce also a noticeable decrease in the normalisation of the M• − M* relation, shown in Fig. 4. The relative difference between the two simulations remains comparable to the scatter of the relation (and so the error in the mean is smaller than the difference in the mean). We find MBHs close to the seed mass all across the galaxy mass spectrum in both simulations, but differences appear more clearly at M* > 1010 M⊙, where in OBELISK-NORECOIL more galaxies host MBHs with mass at least twice the seed mass, although the number of objects at the massive end is low because of the small underlying sample of massive galaxies.

thumbnail Fig. 4.

Correlation between the main MBH mass and the host galaxy stellar mass for OBELISK-RECOIL (red) and OBELISK-NORECOIL (blue). The mean in several galaxy mass bins with the standard deviation is shown in the lines and shaded regions.

Some of our results are qualitatively confirmed by the simulations in which accretion is boosted by setting artificially vrel = 0 in Eq. (1). In this case, recoil also reduces the mass of MBHs and the normalisation of the M• − M* relation, arguing in favour of the robustness of this result. For more details, we refer the reader to Appendix A.

In short, we find that recoils can hinder significantly MBH growth, both by reducing the number of mergers and the accretion rate. Our results agree qualitatively with previous studies using different approaches (Micic et al. 2006; Sijacki et al. 2009, 2011; Volonteri et al. 2011; Blecha et al. 2011; Dunn et al. 2020). This makes it even harder to explain the efficient growth of the MBH population inferred from the large population of very massive and active MBHs detected by JWST (Harikane et al. 2023; Maiolino et al. 2024; Matthee et al. 2024; Greene et al. 2024), unless MBHs do not experience mergers in low-mass galaxies.

3.2. Recoil velocities as a function of MBH evolution

Physically, the recoil velocity is a function of the spin magnitudes and directions, the direction of the orbital angular momentum and the mass ratio. Typically, higher recoil velocities are sourced by mergers with higher spin magnitudes, and higher degrees of misalignment between the two spin directions and between the spin directions and that of the orbital angular momentum. For non-zero spins, the recoil velocity tends to be maximal for mass ratios close to unity. Therefore, the magnitude of the recoil velocity is linked to the accretion and merger history of the MBHs, over which the spin and masses are built, and to the dynamics of the MBHs, which set the direction of the orbital angular momentum and determine the population of MBH pairs that undergo all dynamical stages from the galaxy merger to coalescence.

Figure 5 shows the magnitude of the recoil velocity, vrecoil and the ratio between the recoil velocity and the escape velocity vrecoil/vesc as a function of the host galaxy mass. This is shown for all the mergers in OBELISK-RECOIL. To enlarge the parameter space, we also plot the mergers in the OBELISK simulation, which reaches lower redshift and larger MBH and galaxy masses. For OBELISK we calculate recoil kick magnitudes in post-processing. For OBELISK, we select major mergers with q > 0.3 since below this mass ratio the recoil kicks get arbitrarily small and also mergers that fulfill M2 > 2Mseed and so the MBH has erased the model-dependent initial conditions. In this way, we aim to encompass two opposite physical situations – the MBHs in OBELISK-RECOIL have more realistic dynamics and less growth, while the MBHs in OBELISK are dynamically coupled to the gas and experience more efficient growth (and therefore spin coupling to the angular momentum of the gas). That is, OBELISK-RECOIL corresponds to more random alignment and OBELISK to more efficient alignment through accretion.

thumbnail Fig. 5.

Recoil velocity (top panel) and ratio of the recoil velocity to the escape velocity of the host halo (bottom panel) as a function of the galaxy stellar mass, for OBELISK-RECOIL (red) and OBELISK (green). The recoil kick magnitudes of OBELISK are calculated in post-processing. The mean in several galaxy mass bins with the standard deviation is shown in the lines and shaded regions.

We find that MBHs can receive large recoil kicks after a major merger across a wide range of galaxy masses, MBH masses and accretion rates (MBH masses and accretion rates not shown in the figure). The correlation of vrecoil with these parameters linked to the MBH evolution is generally weak. We find this to be the case in all our simulations – OBELISK-RECOIL, the simulation with the OBELISK-RECOIL model but more efficient accretion (see Appendix A), and OBELISK, which suggests that our result is robust for against several accretion models.

The amount of gas during a galaxy merger could have an impact on recoil kicks – typically, ‘wet’ gas-rich mergers are associated with high spins and aligned configurations as the MBHs couple to the dynamics and angular momentum of the surrounding gas, while ‘dry’ gas-poor mergers are associated with random configurations (Bogdanović et al. 2007; Volonteri 2007; Volonteri et al. 2010; Blecha et al. 2016; Sayeb et al. 2021). We find that in our simulation the distinction between high and low accretion rate mergers (relative to Eddington) is not sufficient to drive a bimodality in the magnitude of the kicks. We find only a very weak correlation between the accretion rate before the merger and the alignment of the orbital and spin angular momenta with the surrounding gas.

Firstly, the dynamical timescale of the MBH merger might not be long enough to align the angular momenta (e.g. Gerosa et al. 2015). However, even if the timescale is long enough, the angular momentum of the nuclear region need not be coherent over time. For low-mass galaxies (M* ≲ 109 M⊙), strong SN feedback prevents coherent rotationally supported structures from forming and the MBHs to accrete efficiently (Dubois et al. 2014b, 2015; Habouzit et al. 2017; Bower et al. 2017; Anglés-Alcázar et al. 2017; Peirani et al. 2024; Beckmann et al. 2025), favouring random spins. In higher-mass galaxies (M* ≳ 109 M⊙), which are the hosts of most of the mergers in our samples, coherent rotationally supported gas can produce efficient mass and spin growth. However, the gas angular momentum is rarely completely coherent, there are generally small fluctuations in the direction over short timescales at high accretion rates (e.g. see figures 1 and 2 in Dong-Páez et al. 2023b). These fluctuations can prevent the high pre-merger alignment of angular momenta needed for low recoil kicks. The bottom-centre panel in Fig. 6 shows the distribution of the alignment between the spins, showing that even if there is a preference for alignment in the simulations, most mergers do not have perfectly aligned spins. The timescale for spin alignment (Eq. (4)) is typically shorter than that for gas dynamical friction, which means that on top of the misalignment between the spins, the spins will also tend to be misaligned with the orbital angular momentum (centre-left panel, Fig. 6). Small misalignments and large spins can still yield large recoil kicks. For example, Lousto et al. (2012) show that rapidly spinning MBHs misaligned by just 15° with the orbital angular momentum and opposite in the azimuthal plane can produce recoil kicks over 1000 km s−1. In Sect. 4 we discuss an example of a merger embedded in a rotationally supported disc that nonetheless receives a large recoil kick.

thumbnail Fig. 6.

Distribution of the cosine of the angle in the time step before the merger between: the gas angular momenta around each MBH (top left), the orbital angular momentum and the gas angular momenta (top middle), the MBH spins and the gas angular momenta (top right), the MBH spins and the orbital angular momentum (bottom left) and the two MBH spins (bottom middle). The MBH spin magnitudes are shown in the bottom right panel. The distributions are shown for OBELISK-RECOIL (red), OBELISK (green), and the ‘hot’ and ‘cold’ disc models in the circumbinary disc simulations of Dotti et al. (2010) (orange and yellow). If both the quantities of the primary and the secondary are shown separately, the primary is shown in solid lines and the secondary in dotted lines.

We note that we are considering numerical mergers, that is, we are ignoring the sub-grid dynamical delay from tens of pc scales to coalescence and the additional alignment that could take place during this evolution. However, if the angular momentum of the circumnuclear discs around bound binaries is preserved from larger scales, the fluctuations in the angular momentum at such scales should also affect the smaller unresolved scales. Furthermore, the timescale for alignment can be longer than the binary evolution timescales (Lodato & Gerosa 2013) and idealised high-resolution simulations of MBHs in a circumnuclear disc by Dotti et al. (2010) show that, even if the angular momentum of the circumnuclear disc is coherent, the MBH spins need not fully align with the orbital angular momentum.

Even if large kicks are produced across all galaxy masses in our sample, the escape velocity increases with the galaxy mass. The ratio vrecoil/vesc of major mergers (Fig. 5) decreases slightly with increasing galaxy mass, meaning that in more massive galaxies recoil kicks for major mergers can have a somewhat more moderate but still significant effect.

Figure 6 investigates further the typical distribution of the spin magnitudes and the relative orientations between the orbital, MBH spin and the angular momentum of the gas around each MBH. The orbital angular momentum and the individual MBH spins tend to be aligned with the ambient gas (top panels). The alignment of each spin with respect to the orbital angular momentum, and the relative alignment between the two MBH spins, which are the terms that enter the calculation of vrecoil, tend to have a somewhat poorer alignment (lower left and middle panels).

Even in the optimistic model of efficient dynamical (and hence angular momentum) coupling to the gas in OBELISK the alignment is high but not perfect. The alignment is lower than in the idealised high-resolution simulations of lower scales (circumbinary discs) by Dotti et al. (2010). We calculated the distributions by drawing random spins and spin vectors from the analytical fits from (Lousto et al. 2012) to the results of (Dotti et al. 2010) in the ‘hot’ and ‘cold’. This is labelled as ‘L12’ in the figure. Finally, in the bottom-right panels we show that in OBELISK, where MBHs grow efficiently, both primaries and secondaries tend to have spin magnitudes close to maximal. As discussed above, this combination of large spin magnitudes and small misalignment can lead to large recoil velocities with vrecoil ≳ 103 km s−1. The spin magnitudes in OBELISK-RECOIL are lower, but alignment is also poorer, as evident comparing the green and red histograms in the lower panels of Fig. 6, resulting in similar kick velocities in the two simulations as reported in Fig. 5.

Overall, we find that high recoil velocities can occur in galaxies and MBHs all across the spectrum in masses and accretion rates found in the simulations. We explored two simulations with very different efficiencies in the coupling of MBH angular momenta with the ambient gas. Even if the coupling is efficient and the nuclear gas is rotationally supported there can be small degrees of misalignment that are enough to still produce a high fraction of high recoil velocities. This result is different from simplified models that assume highly aligned configurations (Bogdanović et al. 2007; Volonteri 2007; Volonteri et al. 2010; Blecha et al. 2016; Dunn et al. 2020; Sayeb et al. 2021). Blecha et al. (2016) and find that large kicks imply optimistic prospects for the detection of recoiling MBHs as offset AGN.

4. Dynamics of recoiling MBHs

Now we turn to study the orbits of recoiling MBHs. Since the MBHs in OBELISK-RECOIL experience limited growth, the sample is composed mostly of major mergers. This makes it suitable for investigating the effect of large recoil kicks with a good statistical sample.

4.1. Do recoiling MBHs escape their host galaxies?

Recoil velocities can be quite large, often exceeding the escape velocity of most MBH merger hosts. Early studies of the cosmological impact of GW recoil using semi-analytic models operated under the assumption that MBHs escape their host galaxies if vrecoil > vesc and remain in the galaxy otherwise (e.g. Haiman 2004; Volonteri & Perna 2005; Schnittman 2007). However, the fate of the MBH does not depend only on the depth of the gravitational well. Dynamical friction, especially from the high stellar densities in galactic cores (e.g. Madau & Quataert 2004), can remove significant energy from the orbits in the early stages of the recoiling trajectory, and at the subsequent pericentres if the orbit remains sufficiently radial. Moreover, realistic galaxies can be irregular and inhomogeneous – MBHs might scatter or interact with clumps, satellite galaxies, and other features (e.g. Fiacconi et al. 2013; Pfister et al. 2019). The halo also grows and increases its escape velocity with time (Choksi et al. 2017). This problem can only be tackled in its full complexity by using cosmological simulations.

We find that, in agreement with simplified prescriptions, MBHs tend to escape the halo rapidly if vrecoil > vesc and stay or return rapidly to the galactic centre otherwise. In Fig. 7, we show the recoil velocity as a function of the escape velocity of the host galaxy. The markers and colours indicate whether the MBH returns to R50 before the end of the simulation, stays within the halo as a wandering MBH or escapes the halo. As physically expected, MBHs with vrecoil < vesc do not escape the galaxy. Above vesc, recoiling MBHs escape their host halo, since dynamical friction is less efficient for fast MBHs that escape the dense regions of the galaxy (Choksi et al. 2017). In total, 37 out of the 92 mergers escape their host haloes.

thumbnail Fig. 7.

Recoil velocity as a function of the escape velocity of the host for the mergers in OBELISK-RECOIL. MBHs that return to R50 are shown in yellow points, wandering MBHs are shown in orange triangles, MBHs that are on their way to escaping are shown in red squares, and MBHs that escape their haloes are shown in purple crosses (see Table 1 for a summary of the definitions). The dotted line is the locus where vrecoil = vesc. Two outliers caused by spurious measurements have been removed.

The time taken for the MBH to return to the centre or to escape behaves asymptotically around vrecoil/vesc ∼ 1. Figure 8 shows the time taken to return to R50 (a lower limit in the case of wandering MBHs that have not returned by the end of the simulation) or the time to escape the halo (a lower limit in the case of clearly unbound MBHs that have not yet escaped). We observe that for vrecoil/vesc < 0.5, treturn is very small or zero – the MBH never escapes R50 or returns in only a few orbits. In the opposite region of the parameter space, for vrecoil/vesc > 2, the MBHs escape quickly, typically in less than 50 Myr.

thumbnail Fig. 8.

Top panel: time to return to R50, treturn, (for returning, in yellow points, and wandering MBHs, in orange triangles) and or time to escape the halo, tesc, (for escaping, in red squares, and escaped MBHs, in purple crosses) as a function of the ratio vrecoil/vesc. For wandering and escaping MBHs the time shown is just the time from the recoil kick to the end of the simulation, and should be taken as a lower limit to treturn or tesc. Bottom panel: radius of the first apocentre of the orbit, in units of R50 as a function of vrecoil/vesc.

Longer return or escape times can only be obtained in a ‘resonant’ regime spanning 0.5 ≲ vrecoil/vesc ≲ 2. Here, the time needed to return to the galaxy or to escape the halo can be quite long, of the order of several hundreds of Myr. A significant fraction of MBHs in this regime do not return to within R50 before the end of the simulation and remain as wandering MBHs. There is still a population of MBHs in this regime that return to R50 or escape in a short time.

This trend in treturn is a reflection of another asymptote in the first apocentric radius of the orbit (Rapo, 0) as a function of vrecoil/vesc. Again, only MBHs with vrecoil/vesc > 0.5 manage to escape the inner regions of the galaxy. In this regime, some of the MBHs manage to dissipate their energy shortly after the kick and stay close to R50, while other MBHs escape R50 and have difficulties sinking back inside. This is because the dynamical friction time for low-energy orbits is short since the MBH tends to stay in the higher-density regions inside the galaxy. For high-energy orbits, the timescale is long, especially if the MBH spends most of its time in the outskirts or outside of the galaxy in low-density regions. This resembles the analytical results by Madau & Quataert (2004) assuming a cored single isothermal sphere and the idealised simulations of galaxy mergers by Blecha et al. (2011) – the radius of the first apocentre diverges as the kick approaches the escape velocity, and so does the dynamical friction time for the MBH to return to the centre.

In Figs. 9, 10 and 11, we show examples of trajectories of MBHs with vrecoil/vesc > 0.5 that return to the centre, wander in the outskirts and escape the halo, respectively. The projections are integrated in the direction of the angular momentum of the galaxy. The galaxies are typically quite compact and dense in the centre, as is common in high-redshift galaxies. The central stellar and gas density tend to be above 102 and 1 M⊙ pc−3 respectively. The orbit is circularised with time. The MBH in Fig. 9 is ejected outside of the nucleus with vrecoil/vesc = 0.55 but remains inside of the galaxy and decays on a short timescale. As shown in the bottom panel it passes through the central high-density regions, dissipating energy progressively. The MBH in Fig. 10 is kicked with vrecoil/vesc = 0.77 to a higher-energy orbit and decays on a longer timescale, remaining as a wandering MBH by the end of the simulation. The orbit is radialised with time. Despite the high central densities of both stars and gas (see bottom panels), the MBH is unable to dissipate energy efficiently. The MBH in Fig. 11 experiences a larger kick with vrecoil/vesc = 1.07 and escapes its host halo. The densities of gas and stars near the MBH decay rapidly as the MBH exits the galaxy.

thumbnail Fig. 9.

Trajectory of a returning MBH, plotted against the projected stellar and gas density (top and middle panels). The colour lightens with time. The bottom panels show the distance to the centre and the density of stars, gas and DM along the trajectory. Light grey, grey, black dotted horizontal lines denote R50, R90, Rvir.

thumbnail Fig. 10.

Similar to Fig. 9 but for a recoiling wandering MBH.

thumbnail Fig. 11.

Similar to Fig. 9 but for an MBH as it escapes the host galaxy.

Dynamical friction is not the only relevant additional effect. The trajectories are often complex. The MBH might interact with clumps and galaxies on its way out – some recoiling MBHs show signs of gravitational interactions with other structures. An example of this is shown in Fig. 12. The MBH is kicked into a low-energy orbit but scatters later off a sub-halo and gains energy. In other cases, the MBH can even be ejected by such interactions. Bound orbits can also radialise or circularise with time, as discussed above.

thumbnail Fig. 12.

Trajectory of a wandering MBHs that interacts with the halo substructure, plotted against the projected stellar density.

In conclusion, in a realistic environment, kicked MBHs tend in most cases to either escape quickly or stay in the centre of the galaxy. There is also a population of MBHs whose fate is difficult to determine without considering the full evolution of the orbit in their complex environment. This is especially the case in the dense and irregular environment characteristic of high-redshift galaxies and haloes.

We highlight that many of our results are also applicable to other possible ejection mechanisms of central MBHs. For example, ejections arising from multiple MBH interactions are likely to be astrophysically relevant (e.g. Volonteri et al. 2003; Hoffman & Loeb 2007; Bonetti et al. 2019; Partmann et al. 2024; Sayeb et al. 2024). Finally, we note that our MBHs release little feedback into their environment, which could reduce the effect of dynamical friction or even reverse it (Sijacki et al. 2011; Park & Bogdanović 2017).

4.2. Recoiling MBHs and other wandering MBHs

We have shown that a significant fraction of recoiling MBHs remain as wandering MBHs. The population of MBHs wandering the galaxy or the halo is not limited to recoiling MBHs. Some other wandering MBHs can originate from former central MBHs that have been displaced (Bellovary et al. 2021), or that have merged into another galaxy and have been stripped of their host galaxies, but dynamical friction has not been efficient enough to sink them to the centre of the potential. Here we investigate whether recoiling MBHs can constitute a significant fraction of the wandering MBH population and whether there are any tell-tale kinematic features distinguishing recoiling MBHs.

In our simulations, recoiling MBHs are a significant fraction of the population of wandering MBHs. Figure 13, shows the distance of MBHs to the centre of the galaxy as a function of the galaxy mass, for the range in galaxy mass in which mergers occur in our simulation. We observe wandering MBHs at a wide range of distances. In general, we see that recoiling MBHs in haloes tend to pile up near Rvir. At lower distances where the dynamical friction times are lower, they represent a smaller fraction compared to wandering MBHs from other origins. In comparison with Untzaga et al. (2024), we find recoiling MBHs located at smaller distances. We do not consider MBHs outside the virial radius which can be defined in different ways in different models.

thumbnail Fig. 13.

Distance from the galactic centre as a function of the host stellar mass at the last output (z = 4.41) for MBHs that have previously experienced a recoil kick, which includes wandering (orange triangles), escaping (red squares), and escaped (purple crosses) MBHs. We also show the distribution of MBHs that have not experienced a merger that are main (dark blue empty circles), satellite (blue empty circles), or halo (light blue empty circles) MBHs. See the text for details of the sample selection. The distance corresponding to 4Δx is denoted with a dashed line.

As discussed in Fig. 8, recoiling MBHs are kicked to orbits with a wide range of apocentric distances. Those with small apocentres have a short lifetime as wandering MBHs since they can decay fast through dynamical friction, especially if they remain in radial orbits inside of the galaxy. Those with large apocentres have long dynamical friction times and it is likely to find them as wandering MBHs. Consequently, we find that most recoiling MBHs are either main MBHs that manage to sink back to the centre or MBHs ejected in the resonant regime identified in Fig. 8 in the outskirts of the halo. There are not many recoiling MBHs at intermediated distances, as satellite MBHs inside galaxies. There is also a population of unbound recoiling MBHs on their way to escape the halo. MBHs on high-energy orbits are likely to be found at larger distances, as they spend more time at larger distances where their velocities are smaller. There is also a population (purple crosses) that have exited the halo but then re-enter the same halo or a different one. OBELISK-RECOIL models an overdense region, making MBHs more likely to enter another halo upon leaving their original halo.

Since MBH mergers and recoil kicks occur in the centre of galaxies, one might expect that the orbits of recoiling MBHs are radial. This could be a clear kinematic property to distinguish recoiling MBHs from other kinds of wandering MBHs. However, in our simulations, we find that the picture is more complex. Figure 14 shows the circularity, calculated from the mass distribution in the simulation assuming a spherically symmetric potential (see Eq. (9)), as a function of the MBH–galaxy distance. The circularity varies from 0 for purely radial orbits to 1 for purely circular orbits. We only show MBHs that are bound in the sphericalised potential. We exclude recoiling MBHs that have escaped a halo.

thumbnail Fig. 14.

Orbit circularity as a function of the distance from the galactic centre for MBHs associated with galaxies with stellar mass larger than 109 M⊙. The MBH decomposition is similar to that of Fig. 13. MBHs that are unbound in a simplified sphericalised potential are not shown. We also exclude recoiling MBHs that have escaped a halo. Only MBHs at distances larger than 2Δx are shown.

Many wandering recoiling MBHs retain radial orbits. However, not all of them do. Triaxial or irregular potentials can circularise the orbit (Madau & Quataert 2004; Guedes et al. 2009). Less bound orbits can also be more susceptible to interacting with satellites and being deflected, especially in the violent environments at high redshifts. Even for the recoiling MBHs that preserve radial orbits, this is not a distinctive feature compared to other wandering MBHs. The orbits of the population of other wandering MBHs can also be radial at large distances since satellite galaxies and other objects tend to enter the haloes on radial orbits, as shown also by Fastidio et al. (2024).

Figure 15 shows the accretion rate as a function of the distance to the centre for the sample of satellite and halo MBHs. The accretion rates drop fast outside of the inner 100 pc with increasing distance. We find again that recoiling MBHs do not show any particular signature in the accretion rate that can distinguish them from other wandering MBHs.

thumbnail Fig. 15.

Eddington ratio as a function of the distance from the galactic centre for MBHs that have previously experienced a recoil kick and associated with in galaxies with stellar mass larger than 109 M⊙. The MBH decomposition is similar to that of Fig. 13. We exclude all main MBHs and recoiling MBHs that return to the centre of their galaxy.

In summary, wandering recoiling MBHs are more likely to be found at large distances, in the outskirts of the halo, compared to other wandering MBHs. Some recoiling MBHs preserve radial orbits, but there are other wandering MBHs with similarly radial orbits at the same distances, making it difficult to discriminate recoiling MBHs just based on dynamics and accretion rate.

5. Conclusions

The effect of GW recoil depends on the cosmic evolution and dynamics of MBHs, and vice versa. The dynamics of an MBH after a recoil kick are complex since MBHs inhabit asymmetric and time-varying potentials. Cosmological simulations are required to study this complex problem. We perform a high-resolution simulation (OBELISK-RECOIL) of an overdense region down to redshift z ∼ 4.4. The simulation includes a model for GW recoil kicks following MBH mergers coupled to detailed sub-grid models for MBH growth, spin evolution and feedback, dynamical friction, as well as other sub-grid models that are relevant for galaxy formation. We also run a twin simulation (Obelisk-noRecoil) without recoil for comparison. We summarise below our main findings.

  • GW recoil decreases the growth of MBH through accretion (Figs. 2 and 3) and ejects MBHs. This leads to a decrease in the distribution of MBH mass (Fig. 1), and a decrease in the merger rate. GW recoil also can also decrease the normalisation of the M• − M* by a factor comparable to the scatter in the relation (Fig. 4).

  • Nearly-equal mass MBH mergers can produce large kicks (vrecoil > 1000 km s−1) which have a large impact on dwarf galaxies, but also on galaxies with stellar masses of 1011 M⊙ or even higher (Fig. 5). Even when MBHs are fed and dynamically coupled to rotationally supported gas, in our simulations we find that the gas is never sufficiently coherent to produce high alignment and low recoil kicks (Fig. 6).

  • After a GW kick, recoiling MBH trajectories behaves asymptotically around vrecoil/vesc = 1. MBHs with low vrecoil/vesc never escape the centre of the galaxy or return quickly driven by the strong gravity and dynamical friction. MBHs with large vrecoil/vesc escape quickly. Only MBH is the ‘resonant’ regime 0.5 ≲ vrecoil/vesc ≲ 2 can take longer times (∼100 Myr) to return or escape (Figs. 7 and 8).

  • Some recoiling MBHs remain in the galaxy or the halo as wandering MBHs by the end of the simulation. Recoiling MBHs constitute a significant fraction of the wandering MBH population. Wandering recoiling MBHs tend to be located at large distances outside of the galaxy, where dynamical friction times are longer. Recoiling MBH orbits are complex – they can change eccentricity with time and get scattered by the substructure in the halo into higher energy orbits or escaping trajectories (Figs. 9, 10, 11, and 12). Some recoiling MBHs maintain their original radial orbits, but this is not enough to distinguish them dynamically from the global population of wandering MBHs, which can also have radial orbits (Fig. 14).

We conclude that the GW recoil can have a significant effect on the cosmic evolution of MBHs for a wide range of astrophysical environments. It is important to include its effect in cosmological simulations that aim at modeling the growth of MBHs. For example, if one uses the local M• − M* relation to calibrate an AGN feedback model (as is often done in the field, e.g. Dubois et al. 2014a) in a simulation that does not account for the suppression of growth from GW recoil, the efficiency of AGN feedback might be overestimated. In the future, we plan to explore the parameter space of MBH evolution models (such as MBH accretion) to understand the effect of MBH evolution on GW recoil, couple GW kicks with the recently developed RAMCOAL model to track subgrid MBH binary dynamics (Li et al. 2024), and include ejections from triple interactions.

Acknowledgments

We thank Francisco Rodríguez Montero, Luc Blanchet, Corentin Cadiou, Yosef Zlochower, Carlos O. Lousto, Eugene Vasiliev, Kunyang Li, Alberto Mangiagli, Nimatou Diallo, and David Trestini for useful discussions. This work was made possible by funding from the French National Research Agency (grant ANR-21-CE31-0026, project MBH_waves). This work was supported by the CNES for the space mission LISA. This work has made use of the Infinity Cluster hosted by Institut d’Astrophysique de Paris; we thank Stéphane Rouberol for running smoothly this cluster for us.

References

  1. Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2023, Liv. Rev. Rel., 26, 2 [NASA ADS] [CrossRef] [Google Scholar]
  2. Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017, MNRAS, 472, L109 [Google Scholar]
  3. Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376 [Google Scholar]
  4. Bardeen, J. M. 1970, Nature, 226, 64 [Google Scholar]
  5. Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65 [Google Scholar]
  6. Barnes, J. E., & Hernquist, L. E. 1991, ApJ, 370, L65 [Google Scholar]
  7. Beckmann, R. S., Dubois, Y., Volonteri, M., et al. 2025, MNRAS, 536, 1838 [Google Scholar]
  8. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
  9. Bekenstein, J. D. 1973, ApJ, 183, 657 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bellovary, J. M., Hayoune, S., Chafla, K., et al. 2021, MNRAS, 505, 5129 [NASA ADS] [CrossRef] [Google Scholar]
  11. Benson, A. J., & Babul, A. 2009, MNRAS, 397, 1302 [NASA ADS] [CrossRef] [Google Scholar]
  12. Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349 [Google Scholar]
  13. Blecha, L., & Loeb, A. 2008, MNRAS, 390, 1311 [NASA ADS] [Google Scholar]
  14. Blecha, L., Cox, T. J., Loeb, A., & Hernquist, L. 2011, MNRAS, 412, 2154 [NASA ADS] [CrossRef] [Google Scholar]
  15. Blecha, L., Sijacki, D., Kelley, L. Z., et al. 2016, MNRAS, 456, 961 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bogdanović, T., Reynolds, C. S., & Miller, M. C. 2007, ApJ, 661, L147 [CrossRef] [Google Scholar]
  17. Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [Google Scholar]
  18. Bonetti, M., Sesana, A., Haardt, F., Barausse, E., & Colpi, M. 2019, MNRAS, 486, 4044 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bower, R. G., Schaye, J., Frenk, C. S., et al. 2017, MNRAS, 465, 32 [NASA ADS] [CrossRef] [Google Scholar]
  20. Bustamante, S., & Springel, V. 2019, MNRAS, 490, 4133 [CrossRef] [Google Scholar]
  21. Campanelli, M., Lousto, C., Zlochower, Y., & Merritt, D. 2007, ApJ, 659, L5 [NASA ADS] [CrossRef] [Google Scholar]
  22. Chen, N., Ni, Y., Holgado, A. M., et al. 2022, MNRAS, 514, 2220 [NASA ADS] [CrossRef] [Google Scholar]
  23. Choksi, N., Behroozi, P., Volonteri, M., et al. 2017, MNRAS, 472, 1526 [Google Scholar]
  24. Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375 [Google Scholar]
  25. Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604 [NASA ADS] [CrossRef] [Google Scholar]
  26. Dong-Páez, C. A., Volonteri, M., Beckmann, R. S., et al. 2023a, A&A, 676, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Dong-Páez, C. A., Volonteri, M., Beckmann, R. S., et al. 2023b, A&A, 673, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Dotti, M., Volonteri, M., Perego, A., et al. 2010, MNRAS, 402, 682 [NASA ADS] [CrossRef] [Google Scholar]
  29. Dubois, Y., Pichon, C., Devriendt, J., et al. 2013, MNRAS, 428, 2885 [NASA ADS] [CrossRef] [Google Scholar]
  30. Dubois, Y., Pichon, C., Welker, C., et al. 2014a, MNRAS, 444, 1453 [Google Scholar]
  31. Dubois, Y., Volonteri, M., & Silk, J. 2014b, MNRAS, 440, 1590 [Google Scholar]
  32. Dubois, Y., Volonteri, M., Silk, J., Devriendt, J., & Slyz, A. 2014c, MNRAS, 440, 2333 [Google Scholar]
  33. Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502 [Google Scholar]
  34. Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Dunn, G., Holley-Bockelmann, K., & Bellovary, J. 2020, ApJ, 896, 72 [NASA ADS] [CrossRef] [Google Scholar]
  36. Fastidio, F., Gualandris, A., Sesana, A., Bortolas, E., & Dehnen, W. 2024, MNRAS, 532, 295 [NASA ADS] [CrossRef] [Google Scholar]
  37. Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [Google Scholar]
  38. Fiacconi, D., Mayer, L., Roškar, R., & Colpi, M. 2013, ApJ, 777, L14 [NASA ADS] [CrossRef] [Google Scholar]
  39. Gerosa, D., Veronesi, B., Lodato, G., & Rosotti, G. 2015, MNRAS, 451, 3941 [Google Scholar]
  40. Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
  41. Greene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39 [CrossRef] [Google Scholar]
  42. Gualandris, A., & Merritt, D. 2008, ApJ, 678, 780 [NASA ADS] [CrossRef] [Google Scholar]
  43. Guedes, J., Madau, P., Kuhlen, M., Diemand, J., & Zemp, M. 2009, ApJ, 702, 890 [NASA ADS] [CrossRef] [Google Scholar]
  44. Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756 [Google Scholar]
  45. Haardt, F., & Madau, P. 1996, ApJ, 461, 20 [Google Scholar]
  46. Habouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935 [NASA ADS] [CrossRef] [Google Scholar]
  47. Haiman, Z. 2004, ApJ, 613, 36 [Google Scholar]
  48. Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 39 [NASA ADS] [CrossRef] [Google Scholar]
  49. Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89 [Google Scholar]
  50. Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29 [Google Scholar]
  51. Hernquist, L. 1989, Nature, 340, 687 [Google Scholar]
  52. Hoffman, L., & Loeb, A. 2007, MNRAS, 377, 957 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hoyle, F., & Lyttleton, R. A. 1939, Proc. Camb. Philos. Soc., 35, 405 [CrossRef] [Google Scholar]
  54. Katz, M. L., Kelley, L. Z., Dosopoulou, F., et al. 2020, MNRAS, 491, 2301 [NASA ADS] [Google Scholar]
  55. Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2 [Google Scholar]
  56. Khan, F. M., Just, A., & Merritt, D. 2011, ApJ, 732, 89 [Google Scholar]
  57. Kimm, T., & Cen, R. 2014, ApJ, 788, 121 [NASA ADS] [CrossRef] [Google Scholar]
  58. King, A. R., Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2005, MNRAS, 363, 49 [NASA ADS] [CrossRef] [Google Scholar]
  59. Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [Google Scholar]
  60. Komossa, S. 2012, Adv. Astron., 2012, 364973 [NASA ADS] [CrossRef] [Google Scholar]
  61. Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [Google Scholar]
  62. Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581 [Google Scholar]
  63. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  64. Lapiner, S., Dekel, A., & Dubois, Y. 2021, MNRAS, 505, 172 [NASA ADS] [CrossRef] [Google Scholar]
  65. Li, K., Bogdanović, T., Ballantyne, D. R., & Bonetti, M. 2022, ApJ, 933, 104 [NASA ADS] [CrossRef] [Google Scholar]
  66. Li, K., Volonteri, M., Dubois, Y., Beckmann, R., & Trebitsch, M. 2024, A&A, submitted [arXiv:2410.07856] [Google Scholar]
  67. Lodato, G., & Gerosa, D. 2013, MNRAS, 429, L30 [NASA ADS] [CrossRef] [Google Scholar]
  68. Lousto, C. O., & Zlochower, Y. 2008, Phys. Rev. D, 77, 044028 [Google Scholar]
  69. Lousto, C. O., & Zlochower, Y. 2009, Phys. Rev. D, 79, 064018 [Google Scholar]
  70. Lousto, C. O., & Zlochower, Y. 2011, Phys. Rev. Lett., 107, 231102 [NASA ADS] [CrossRef] [Google Scholar]
  71. Lousto, C. O., Zlochower, Y., Dotti, M., & Volonteri, M. 2012, Phys. Rev. D, 85, 084015 [NASA ADS] [CrossRef] [Google Scholar]
  72. Lupi, A., Haardt, F., & Dotti, M. 2015, MNRAS, 446, 1765 [Google Scholar]
  73. Madau, P., & Quataert, E. 2004, ApJ, 606, L17 [CrossRef] [Google Scholar]
  74. Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285 [Google Scholar]
  75. Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024, A&A, 691, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Mannerkoski, M., Johansson, P. H., Rantala, A., et al. 2022, ApJ, 929, 167 [NASA ADS] [CrossRef] [Google Scholar]
  77. Marconi, A., & Hunt, L. K. 2003, ApJ, 589, L21 [Google Scholar]
  78. Matthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129 [NASA ADS] [CrossRef] [Google Scholar]
  79. McAlpine, S., Bower, R. G., Rosario, D. J., et al. 2018, MNRAS, 481, 3118 [NASA ADS] [CrossRef] [Google Scholar]
  80. McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]
  81. Micic, M., Abel, T., & Sigurdsson, S. 2006, MNRAS, 372, 1540 [Google Scholar]
  82. Micic, M., Holley-Bockelmann, K., & Sigurdsson, S. 2011, MNRAS, 414, 1127 [Google Scholar]
  83. Milosavljević, M., & Merritt, D. 2001, ApJ, 563, 34 [Google Scholar]
  84. Nasim, I. T., Gualandris, A., Read, J. I., et al. 2021, MNRAS, 502, 4794 [NASA ADS] [CrossRef] [Google Scholar]
  85. Ogiya, G., & Nagai, D. 2024, MNRAS, 527, 5503 [Google Scholar]
  86. Ostriker, E. C. 1999, ApJ, 513, 252 [Google Scholar]
  87. Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40 [Google Scholar]
  88. Park, K., & Bogdanović, T. 2017, ApJ, 838, 103 [NASA ADS] [CrossRef] [Google Scholar]
  89. Partmann, C., Naab, T., Rantala, A., et al. 2024, MNRAS, 532, 4681 [NASA ADS] [CrossRef] [Google Scholar]
  90. Peirani, S., Suto, Y., Beckmann, R. S., et al. 2024, A&A, 686, A233 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Peng, C. Y. 2007, ApJ, 671, 1098 [NASA ADS] [CrossRef] [Google Scholar]
  92. Perego, A., Dotti, M., Colpi, M., & Volonteri, M. 2009, MNRAS, 399, 2249 [NASA ADS] [CrossRef] [Google Scholar]
  93. Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101 [Google Scholar]
  94. Rantala, A., Pihajoki, P., Johansson, P. H., et al. 2017, ApJ, 840, 53 [Google Scholar]
  95. Rawlings, A., Keitaanranta, A., Mattero, M., et al. 2025, MNRAS, 537, 3421 [Google Scholar]
  96. Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82 [NASA ADS] [CrossRef] [Google Scholar]
  97. Rezzolla, L., Barausse, E., Dorband, E. N., et al. 2008, Phys. Rev. D, 78, 044002 [NASA ADS] [CrossRef] [Google Scholar]
  98. Rosdahl, J., & Blaizot, J. 2012, MNRAS, 423, 344 [Google Scholar]
  99. Sala, L., Valentini, M., Biffi, V., & Dolag, K. 2024, A&A, 685, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Sayeb, M., Blecha, L., Kelley, L. Z., et al. 2021, MNRAS, 501, 2531 [CrossRef] [Google Scholar]
  101. Sayeb, M., Blecha, L., & Kelley, L. Z. 2024, MNRAS, 527, 7424 [Google Scholar]
  102. Schmidt, M. 1959, ApJ, 129, 243 [NASA ADS] [CrossRef] [Google Scholar]
  103. Schnittman, J. D. 2007, ApJ, 667, L133 [NASA ADS] [CrossRef] [Google Scholar]
  104. Sijacki, D., Springel, V., & Haehnelt, M. G. 2009, MNRAS, 400, 100 [NASA ADS] [CrossRef] [Google Scholar]
  105. Sijacki, D., Springel, V., & Haehnelt, M. G. 2011, MNRAS, 414, 3656 [Google Scholar]
  106. Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776 [Google Scholar]
  107. Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [Google Scholar]
  108. Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
  109. Thorne, K. S. 1974, ApJ, 191, 507 [Google Scholar]
  110. Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2021, A&A, 653, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, ApJ, 574, 740 [NASA ADS] [CrossRef] [Google Scholar]
  112. Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Untzaga, J., Bonoli, S., Izquierdo-Villalba, D., Mezcua, M., & Spinoso, D. 2024, MNRAS, 535, 3293 [Google Scholar]
  114. van Leer, B. 1979, J. Comput. Phys., 32, 101 [Google Scholar]
  115. Volonteri, M. 2007, ApJ, 663, L5 [NASA ADS] [CrossRef] [Google Scholar]
  116. Volonteri, M., & Perna, R. 2005, MNRAS, 358, 913 [CrossRef] [Google Scholar]
  117. Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559 [Google Scholar]
  118. Volonteri, M., Gültekin, K., & Dotti, M. 2010, MNRAS, 404, 2143 [NASA ADS] [Google Scholar]
  119. Volonteri, M., Natarajan, P., & Gültekin, K. 2011, ApJ, 737, 50 [NASA ADS] [CrossRef] [Google Scholar]
  120. Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219 [Google Scholar]
  121. White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [Google Scholar]
  122. Zlochower, Y., Campanelli, M., & Lousto, C. O. 2011, CQG, 28, 114015 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Simulations with vrel = 0

We ran twin simulations of OBELISK-RECOIL and OBELISK-NORECOIL varying the accretion model, namely setting vrel = 0 in equation 1.

Figure A.1 shows the M• − M* relation for the vrel = 0 and the fiducial sets of simulations. In both cases, we find that the number of the most massive MBHs is reduced in the recoil compared to the no recoil simulations, in both the vrel = 0 and fiducial cases. Similarly, with vrel = 0, the normalisation of the M• − M* relation is reduced in the recoil simulations compared to the no recoil simulations.

Figure A.2 shows the recoil velocity as a function of the host galaxy stellar mass. This confirms the trend found for OBELISK-RECOIL and OBELISK – large recoil kicks can occur in a variety of environments and galaxy masses.

Overall, we find qualitatively similar results when the model is varied to make accretion more efficient. This suggests that the results presented in the main body of the paper are robust.

thumbnail Fig. A.1.

Similar to Fig. 4. The correlation between the main MBH mass and the host galaxy stellar mass for the vrel = 0 versions of OBELISK-RECOIL (red) and OBELISK-NORECOIL (blue) (top panel) and the fiducial versions (bottom panel) at z = 6.01.

thumbnail Fig. A.2.

Similar to Fig. 5. Recoil velocity as a function of the galaxy stellar mass, for the major mergers in the vrel = 0 version of OBELISK-RECOIL.

All Tables

Table 1.

Dynamical classification for recoiling MBH that we used in this study.

In the text

All Figures

thumbnail Fig. 1.

Mass of the main MBHs in the OBELISK-RECOIL (red), OBELISK-NORECOIL (blue) and OBELISK simulations at z = 4.41.
In the text
thumbnail Fig. 2.

Fraction of mass gained through mergers against the total mass gained above the seed mass for OBELISK-RECOIL (red), OBELISK-NORECOIL (blue), and OBELISK (green).
In the text
thumbnail Fig. 3.

Distribution of the accretion rate normalised by the Eddington limit for OBELISK-RECOIL and OBELISK-NORECOIL, averaged over the 5 Myr after the merger for the remnants. The first bin shows the number of MBHs below the lower limit of the plot.
In the text
thumbnail Fig. 4.

Correlation between the main MBH mass and the host galaxy stellar mass for OBELISK-RECOIL (red) and OBELISK-NORECOIL (blue). The mean in several galaxy mass bins with the standard deviation is shown in the lines and shaded regions.
In the text
thumbnail Fig. 5.

Recoil velocity (top panel) and ratio of the recoil velocity to the escape velocity of the host halo (bottom panel) as a function of the galaxy stellar mass, for OBELISK-RECOIL (red) and OBELISK (green). The recoil kick magnitudes of OBELISK are calculated in post-processing. The mean in several galaxy mass bins with the standard deviation is shown in the lines and shaded regions.
In the text
thumbnail Fig. 6.

Distribution of the cosine of the angle in the time step before the merger between: the gas angular momenta around each MBH (top left), the orbital angular momentum and the gas angular momenta (top middle), the MBH spins and the gas angular momenta (top right), the MBH spins and the orbital angular momentum (bottom left) and the two MBH spins (bottom middle). The MBH spin magnitudes are shown in the bottom right panel. The distributions are shown for OBELISK-RECOIL (red), OBELISK (green), and the ‘hot’ and ‘cold’ disc models in the circumbinary disc simulations of Dotti et al. (2010) (orange and yellow). If both the quantities of the primary and the secondary are shown separately, the primary is shown in solid lines and the secondary in dotted lines.
In the text
thumbnail Fig. 7.

Recoil velocity as a function of the escape velocity of the host for the mergers in OBELISK-RECOIL. MBHs that return to R50 are shown in yellow points, wandering MBHs are shown in orange triangles, MBHs that are on their way to escaping are shown in red squares, and MBHs that escape their haloes are shown in purple crosses (see Table 1 for a summary of the definitions). The dotted line is the locus where vrecoil = vesc. Two outliers caused by spurious measurements have been removed.
In the text
thumbnail Fig. 8.

Top panel: time to return to R50, treturn, (for returning, in yellow points, and wandering MBHs, in orange triangles) and or time to escape the halo, tesc, (for escaping, in red squares, and escaped MBHs, in purple crosses) as a function of the ratio vrecoil/vesc. For wandering and escaping MBHs the time shown is just the time from the recoil kick to the end of the simulation, and should be taken as a lower limit to treturn or tesc. Bottom panel: radius of the first apocentre of the orbit, in units of R50 as a function of vrecoil/vesc.
In the text
thumbnail Fig. 9.

Trajectory of a returning MBH, plotted against the projected stellar and gas density (top and middle panels). The colour lightens with time. The bottom panels show the distance to the centre and the density of stars, gas and DM along the trajectory. Light grey, grey, black dotted horizontal lines denote R50, R90, Rvir.
In the text
thumbnail Fig. 10.

Similar to Fig. 9 but for a recoiling wandering MBH.
In the text
thumbnail Fig. 11.

Similar to Fig. 9 but for an MBH as it escapes the host galaxy.
In the text
thumbnail Fig. 12.

Trajectory of a wandering MBHs that interacts with the halo substructure, plotted against the projected stellar density.
In the text
thumbnail Fig. 13.

Distance from the galactic centre as a function of the host stellar mass at the last output (z = 4.41) for MBHs that have previously experienced a recoil kick, which includes wandering (orange triangles), escaping (red squares), and escaped (purple crosses) MBHs. We also show the distribution of MBHs that have not experienced a merger that are main (dark blue empty circles), satellite (blue empty circles), or halo (light blue empty circles) MBHs. See the text for details of the sample selection. The distance corresponding to 4Δx is denoted with a dashed line.
In the text
thumbnail Fig. 14.

Orbit circularity as a function of the distance from the galactic centre for MBHs associated with galaxies with stellar mass larger than 109 M⊙. The MBH decomposition is similar to that of Fig. 13. MBHs that are unbound in a simplified sphericalised potential are not shown. We also exclude recoiling MBHs that have escaped a halo. Only MBHs at distances larger than 2Δx are shown.
In the text
thumbnail Fig. 15.

Eddington ratio as a function of the distance from the galactic centre for MBHs that have previously experienced a recoil kick and associated with in galaxies with stellar mass larger than 109 M⊙. The MBH decomposition is similar to that of Fig. 13. We exclude all main MBHs and recoiling MBHs that return to the centre of their galaxy.
In the text
thumbnail Fig. A.1.

Similar to Fig. 4. The correlation between the main MBH mass and the host galaxy stellar mass for the vrel = 0 versions of OBELISK-RECOIL (red) and OBELISK-NORECOIL (blue) (top panel) and the fiducial versions (bottom panel) at z = 6.01.
In the text
thumbnail Fig. A.2.

Similar to Fig. 5. Recoil velocity as a function of the galaxy stellar mass, for the major mergers in the vrel = 0 version of OBELISK-RECOIL.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.