© The Authors 2024

1 Introduction

Supermassive black holes (SMBHs) reside at the centres of massive galaxies and are considered to affect their evolution profoundly. Numerous studies provided evidence of the relation between the mass of a SMBH and the properties of its host galaxy (e.g. Kormendy et al. 1993; Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Merritt & Ferrarese 2001; Ferrarese et al. 2001; Haering & Rix 2004; Gültekin et al. 2009; Graham & Driver 2007; McConnell & Ma 2013; Gaspari et al. 2019). The most widely accepted explanation is that during SMBH growth by gas accretion, a small fraction of the enormous amount of the released gravitational energy couples with the surrounding environment, regulating the galaxy star formation via various possible mechanisms (e.g. Silk & Rees 1998; Granato et al. 2004; Hopkins et al. 2005; Bower et al. 2006; Cattaneo et al. 2009; Gitti et al. 2012).

Given the influence the SMBHs play in shaping the environment where they reside, it is essential to trace their dynamics correctly. Massive BHs are thought to be formed at some early epochs through mechanisms ranging from direct collapse of primodial gas clouds or as the end stage of very massive Population III stars (e.g. Bromm & Loeb 2003; Volonteri & Bellovary 2012; Mayer & Bonoli 2018), and to subsequently grow in the dense cores of galaxies. Nonetheless, recent studies have consistently shown cases of SMBHs exhibiting substantial displacements from their host galaxies (e.g. Webb et al. 2012; Menezes et al. 2014; Combes et al. 2019; Reines et al. 2020). The dynamical behaviour of SMBHs is significantly affected by the dynamical friction (DF) force (Chandrasekhar 1943; Binney & Tremaine 2008) exerted by the matter distribution surrounding them. This drag force in general prevents an SMBH from escaping the centre of its host galaxy, is responsible for the migration of BH to the galaxy centre, and drives the initial stages of a merger event between two SMBHs, ultimately leading to the formation of a close pair (Begelman et al. 1980).

Cosmological simulations represent ideal tools for following the evolution of structures given the highly nonlinear astrophysical phenomena governing the interactions of baryonic matter – including gas and stars – and also interactions with dark matter (DM). N-body simulations describe the gravitational instability of a collisionless fluid, which is sampled by a discrete set of ‘macro particles’ (e.g. Borgani & Kravtsov 2011; Springel 2016, for reviews). For this reason, tracking the orbits of single collisionless particles has little physical relevance, as it is their ensemble properties that carry the most importance. A BH particle, on the other hand, is introduced in cosmological simulations as an individual collisionless particle, and has a specific physical meaning; its presence is capable of significantly impacting the global properties of the galaxy it belongs to. Unlike the surrounding macro-particles, BHs are not entities whose motion can be interpreted solely on a global scale. Instead, their motion mirrors the effective motion of an astrophysical object. The contrast between the nature of BH particles in these simulations and the interpretation of the surrounding ones presents a primary conceptual obstacle when tracking the trajectories of BH particles.

Specifically, scattering interactions between a BH and its surrounding particles can ‘heat’ the BH orbit. A spurious, numerical displacement of a BH from the centre of the host galaxy is a major consequence of this heating, and can eventually lead to the formation of unwanted ‘wandering’ BHs. Furthermore, this displacement also negatively impacts the capability of simulations to describe BH–BH merger events, ultimately leading to an incorrect description of AGN feedback, which in turn strongly affects the predictions of the simulations. For instance, Ragone-Figueroa et al. (2018) found that a better centring of the SMBH within the host galaxy in cluster simulations was key to predicting the masses of the brightest cluster galaxies (BCGs; confirmed by comparison with observations), which were overpredicted by a factor of a few in their previous work (Ragone-Figueroa et al. 2013).

The issues described above stem from numerical simulations failing to recover the DF force. As the magnitude of the drag due to DF depends on the interactions between BHs and the surrounding particles, any limitation in reconstructing the correct gravitational interactions at the N-body resolution level results in an inaccurate representation of this effect. In this context, the first question that we attempt to answer here pertains to whether it is feasible to introduce a correction to the gravitational acceleration provided by the N-body solver that accounts for the unresolved DF.

Instead of relying on some numerical ‘tricks’ to control the dynamics of the BHs, such as artificially repositioning the BHs at the position of a local minimum of the gravitational potential (Springel et al. 2005; Di Matteo et al. 2008; Sijacki et al. 2015; Davé et al. 2019; Ragone-Figueroa et al. 2018; Bassini et al. 2019; Bahé et al. 2022), or using a boosted dynamical mass to enhance the effect of the resolved DF (DeBuhr et al. 2011; Bassini et al. 2020), some authors have already suggested addressing this problem by introducing an explicit correction for the unresolved DF (Hirschmann et al. 2014; Tremmel et al. 2015; Bird et al. 2022; Chen et al. 2022; Ma et al. 2023). For instance, Hirschmann et al. (2014) proposed the application of a DF correction given by the Chandrasekhar DF formula. In their application, the maximum impact parameter entering in the DF correction is the half-mass radius of the substructure hosting the BH. In addition, Tremmel et al. (2015) stated that, under the hypotheses of a sufficiently shallow potential surrounding the BH, only unresolved interactions within the softening length require correction. Still, Chen et al. (2022) proposed that the application of the DF correction should be coupled with a boosted dynamical mass in order to account for those cases in which a BH has a rather small mass, close to its value at seeding, and is located in a poorly resolved halo. The absence of a consensus on the possibility, use, and actual computation of a DF correction to improve the description of the dynamics of BHs in cosmological simulations is what motivated this study.

In this paper, we propose a novel implementation of the DF correction, which we realise using the OpenGadget3 code for cosmological N-body and hydrodynamical simulations. In this implementation, a correction to DF force is computed by explicitly accounting for the contributions of numerical particles whose gravitational interactions with a BH particle – as provided by the N-body solution – are directly affected by force softening. As we extensively discuss in the following, the primary advantage of our approach is that it is less affected by the assumptions on which the derivation of the Chandrasekhar DF formula is based. In more detail, we aim to address the following questions: Does our approach provide an adequate description of the DF force acting on BHs? How does it compare against the other numerical ah hoc prescriptions (i.e. repositioning and dynamical mass) introduced to mimic the effect of DF on BH particles? To answer these questions, we simulate a group-sized halo and a cluster-sized halo, with initial conditions from the DIANOGA set (Bonafede et al. 2011; Bassini et al. 2019, 2020), along with a cosmological box with a comoving size of 16 Mpc per side (16 cMpc). The simulations were carried out three times, maintaining identical settings for all parameters, except for the sub-resolution prescription governing the BH dynamics: continuous repositioning on a local potential minimum, boosted dynamical mass, and our novel model to correct the unresolved DF. This approach enables us not only to focus on systems densely populated and rich in interactions, such as galaxy clusters and groups, but also to carry out a statistical analysis of the BH population in a cosmological volume.

The paper is organised as follows. In Sect. 2 we present the new DF model and compare it with both previous DF corrections and numerical prescriptions to constrain the BH dynamics. In Sect. 3 we describe the implementation of SMBH evolution and the ensuing AGN feedback in the OpenGadget3 code. After introducing the details of our test simulations in Sect. 4, we present the results of our analysis of the general properties of the SMBH population in Sect. 5. In Section 6 we zoom in, limiting our analysis to single BHs and merger episodes, and compare the small-scale effects of our sub-resolution model of DF.

2 Dynamics of BH particles

Due to the finite mass and force resolution of cosmological simulations, the effect of DF exerted on a BH particle by surrounding particles is always underestimated and, in general, subject to discreteness noise. Such limitations can lead to a grossly incorrect description of the orbits of BH particles, which leads in turn to an incorrect description of the ensuing AGN feedback and of the predictions on the SMBH population. We refer to Sect. 3 for details on the BH seeding, accretion and feedback mechanisms in cosmological hydrodynamical simulations. In this section we present different approaches, including our new one, which have been introduced to overcome this limitation. In Sect. 2.1, we present our new description of a correction to the gravitational acceleration to account for the unresolved DF. Previous works already proposed approaches to correct DF for finite-resolution effects, and some details of our implementation revisit their arguments. For this reason, we highlight in Sect. 2.2 the conceptual differences between our method and such previous approaches proposed in the literature. Finally, Sect. 2.2.1 and Sect. 2.2.2 describe other methods, based on ad hoc prescriptions, to correct BH dynamics for the unresolved DF.

2.1 A new prescription for dynamical friction

To clearly describe in detail the DF correction proposed in this work, we review the basic steps of the DF expression originally derived by Chandrasekhar (1943), which holds under the assumption of specific hypotheses. The starting point is a system of two particles moving on a Keplerian orbit around each other. The velocity variation ∆v_M of a particle of mass M caused by the interaction with another particle of mass m and velocity v can be expressed in terms of their impact parameter b and relative velocity, v₀ = v − v_M , when the two particles were initially at large distance. The components of ∆v_M parallel and perpendicular to the direction of u0 can be expressed as (Binney & Tremaine 2008): $${\left|\text{Δ}{v}_{\text{M}}\right|}_{\Vert }=\frac{2m\left|v_0\right|}{(M+m)}{\left[1+\frac{b^2{\left|v_0\right|}^4}{G^2{(\text{M}+m)}^2}\right]}^{-1};$$(1) $${\left|\text{Δ}{v}_{\text{M}}\right|}_{\perp }=\frac{2mb{\left|v_0\right|}^3}{G{(\text{M}+m)}^2}{\left[1+\frac{b^2{\left|v_0\right|}^4}{G{(\text{M}+m)}^2}\right]}^{-1}.$$(2)

If the particle M moves in a ‘sea’ of particles of mass m, then the DF force acting on the former arises as the sum of the contributions of the velocity variations given by eqs. (1) and (2), due to the interactions with all the surrounding particles. In the derivation by Binney & Tremaine (2008), the mass m of the ‘sea’ particles is assumed to be the same for all such particles. Under the assumption that the distribution of particles is uniform around the BH, the perpendicular contributions of Eq. (2) sum to zero. The rate of encounters having impact parameter in the range [b, b + db] is then given by: 2π b db |∆v_M |_‖f (v)d³v, where f(v) is the phase-space number density of stars. Integrating this rate over the impact parameter from 0 to b_max ¹, we have that the DF force acting on a BH of mass M from particles having mass m and velocities in the range (v, v + d³v) is: $${\frac{d{v}_{\text{M}}}{dt}|}_v=2\pi \ln \left[1+\text{Λ}{(m,v)}^2\right]G^{2}m(\text{M}+m)f(v)d^{3}v\frac{\left(v-v_{\text{M}}\right)}{{\left|v-v_{\text{M}}\right|}^3},$$(3)

where $$\text{Λ}(m,v)=\frac{b_{\max }{\left(v-v_{\text{M}}\right)}^2}{G(\text{M}+m)}.$$(4)

Here, b_max is the maximum impact parameter. In general, b_max is interpreted as the largest distance from the target BH particle, which contains all the particles contributing to the DF exerted on the BH itself. In the above expression, we have made explicit the dependence of Λ on mass m and velocity v.

The maximum allowed value of the impact parameter should in principle be set by the size of the system containing all the particles contributing to the DF exerted on a BH particle. In fact, different assumptions for the values of b_max have been introduced in the literature, all having a degree of arbitrariness. In the next section, we provide an extended comparison between different choices. As in Tremmel et al. (2015), we assume b_max to be given by the gravitational softening length of the BH particle, є_BH , meaning that above such length, the N-body solver is assumed to provide already a correct description of DF. In the context of our simulations, the BH is surrounded by particles which trace the underlying continuous density field. Each of these particles has its mass, m_i, and its velocity v_m,i . The phase-space number density of particles surrounding the BH can be then expressed as a sum of delta-functions, each with a mass m_i and a velocity v_m,i: $$f(v)=\frac{3}{4\pi {ϵ}_{\text{BH}}{}^3}{\displaystyle \sum _{i=1}^{N\left(<{ϵ}_{\text{BH}}\right)}\delta }\left(v-v_{m,i}\right),$$(5)

where the sum is over all the N(< є_BH) particles lying at a distance from the BH smaller than its gravitational softening scale.

Equation (5) relies on the hypothesis that the particles within the softening length provide an adequate sampling of the velocity field of the ‘sea’ surrounding the BHs. To validate this assumption, we performed several tests presented in Appendix A. For well-resolved sub-halos, where the number of particles within the softening length can be as large as ~ 10⁴, this assumption is valid.

In the simulations performed in this work, we only use DM and star particles to compute the correction to DF, while we defer to a forthcoming analysis the inclusion of gas particles (see for example Dubois et al. 2014). Using Eqs. (3) and (5), we can then derive the total DF force term, by integrating over the surrounding particles’ velocities: $$\frac{d{v}_{\text{M}}}{dt}=\frac{3G^2}{2{ϵ}_{\text{BH}}{}^3}{\displaystyle \sum _{i=1}^{N\left(<{ϵ}_{\text{BH}}\right)}\ln }\left[1+\text{Λ}{\left(m_i\right)}^2\right]m_i\left(\text{M}+m_i\right)\frac{\left(v_{m,i}-v_{\text{M}}\right)}{{\left|v_{m,i}-v_{\text{M}}\right|}^3}.$$(6)

The gravitational acceleration of the BH is then corrected by the DF acceleration a_df given by Eq. (6), so that the total acceleration acting on a single BH particle is: $$a_{\text{BH}}=a_{\text{g}}+a_{\text{df}},$$(7)

where a_g is the acceleration provided by the N-body solver. It is important to remark that the DF correction hence obtained is a contribution correcting for the softened interactions between the BH and the surrounding particles, but not for the absence of information on the sub-softening structure of phase-space.

In order to validate the performance of this model in a controlled numerical experiment, without the complexity of a full-physics simulation in a cosmological context, we carry out several tests for the sinking timescale of a BH initially placed on a circular orbit and infalling toward the centre of an isolated DM halo, and compare them to analytical predictions. The results of these tests are reported in Appendix B. They show in general that simulations which include our correction for unresolved DF are in very good agreement with analytical predictions, and the convergence with resolution to this analytical result is faster than when DF correction is neglected.

We also point out that our approach for the computation of the DF correction presents some advantages with respect to the models proposed in the literature. We discuss this in the next section, providing an overview of the implementations adopted nowadays and underlining their differences with our approach.

2.2 Previous approaches

Correcting the unresolved DF force acting on BH particles using a physically motivated approach instead of resorting to ad hoc prescriptions is clearly highly desirable. This is the reason why, the implementation of a correction of the DF force acting on BH particles, based on the original derivation by Chandrasekhar (1943), has been explored in various studies since the very first implementations (e.g. Hirschmann et al. 2014; Tremmel et al. 2015). A common aspect of all such approaches is that they start from Eq. (3), and implicitly assume that the BH is surrounded by a homogeneous and infinite distribution of particles all having the same mass. However, different approaches make different assumptions of the physical size of the region where to correct for the unresolved DF, i.e. on the actual value of b_max. For instance, in the work by Hirschmann et al. (2014), b_max is defined as the typical size of the system hosting the BH and, as such, it is set to the half-mass radius of the sub-halo². On the other hand, Tremmel et al. (2015) argued that DF is correctly computed by the N-body solver at scales larger than the gravitational softening length of the BH, є_BH. Accordingly, a correction term to the DF should be added only on scales smaller than є_BH , so that b_max = є_BH . In line with this approach, Pfister et al. (2019) accounted for the DF force from particles within a sphere centred on the BH position and with radius which is a multiple of the adaptive grid mesh size on which the gravitational force is computed (Teyssier 2002). Finally, the DF model recently presented by Chen et al. (2022) and Bird et al. (2022) assumes a constant value of b_max = 10h⁻¹ ckpc and 20 kpc, respectively. As for the velocity distribution function of the sea of particles around the BH, Hirschmann et al. (2014), Chen et al. (2022) and Bird et al. (2022) adopt the same standard hypothesis, originally formulated by Binney & Tremaine (2008), of a local Maxwellian velocity distribution. Under the further hypothesis that M_BH≫ m, where m is the mass of the particles around the BH, the DF force F_DF can be cast in the form $$F_{\text{DF}}=-4\pi \rho {\left(\frac{G{\text{M}}_{\text{BH}}}{v_{\text{BH}}}\right)}^2 \text{F}(x)\ln (\text{Λ}){\widehat{v}}_{\text{BH}}.$$(8)

Here, ρ is the smoothed density at the position of the BH, contributed by stellar and DM particles, using the BH smoothing length. Furthermore,${\widehat{v}}_{\text{BH}}$ is the versor of the BH velocity relative to the ‘sea’ of surrounding particles, while $$\text{F}(x)=\mathrm{erf}(x)-\frac{2x}{\sqrt{\pi }}{e}^{-x^2};x=\frac{v_{\text{BH}}}{{\sigma }_v},$$(9)

with σ_v the velocity dispersion of the surrounding particles. We assume that stars and DM particles exert a DF force on the BHs. Ostriker (1999) computed the contribution to DF from gas, lately included in simulations as an additional numerical corrective term by Chen et al. (2022). However, in their analysis, Chen et al. (2022) found that the DF correction is in fact dominated by the collisionless component.

Rather than assuming a specific expression for the velocity distribution function of the particles around a BH, Tremmel et al. (2015) proposed to incorporate the mass m of the surrounding particles within the integral over velocities, thus moving the uncertainty on the velocity distribution function to an uncertainty on the mass density of the surrounding particles. The approach by Tremmel et al. (2015) has been then applied also in Bellovary et al. (2018). They found that the effect of DF correction is efficient for BHs having a mass at least three times larger than that of the surrounding particles. This result justifies the choice made by Chen et al. (2022) to include both a DF correction and a boosted dynamical mass for BH particles (see Sect. 2.2.2 below).

Finally, a scheme that stands out from the others is the one proposed by Ma et al. (2023). In their approach, a discrete N-body correction, similar to the one proposed in this paper, is taken into account, but still acting only on scales above the gravitational force resolution of simulations. Differently from such an approach, the model that we described in Sect. 2.1 explicitly intends to correct for the interactions that take place below the BH softening scale which, by definition, are not correctly described by the N-body solver. Overall, the main differences between the DF correction proposed in this work and the previous ones are the absence of any a priori assumption on the velocity distribution of the surrounding particles and the relaxation of the hypothesis that the surrounding particles’ mass is negligible compared to the BH one. Our approach accounts for single scatters between the BH and each particle within the BH softening, each contributing with its own velocity and mass. In this way, we try to reduce the dynamical heating of BHs as a consequence of the noisy background potential, taking place whenever their mass is comparable to or even lower than the surrounding particles’ one. In summary, Table 1 lists the main differences between our novel DF correction introduced in Sect. 2.1 and those proposed by Hirschmann et al. (2014), Tremmel et al. (2015), and Chen et al. (2022).

2.2.1 Repositioning BH particles

Among the major issues encountered when introducing BH particles in N-body simulations is that they can escape from the centres of the host galaxies. One of the most widely used methods to avoid this consists of repositioning, at each time step, the BH particle at the position of the most bound particle among its neighbours. Different implementations of this method feature different choices for the search radius. Moreover, such alternative implementations often adopt additional constrains (e.g. on their relative velocity) for the selection of the neighbour particle on which to relocate the BH. In this way, the BH particle is generally forced to remain at the centre of its host sub-halo (see e.g. Di Matteo et al. 2008; Booth & Schaye 2009; Vogelsberger et al. 2013; Sijacki et al. 2015; Schaye et al. 2014; Pillepich et al. 2017; Ragone-Figueroa et al. 2018; Bahé et al. 2022).

As pointed out by Tremmel et al. (2015), this method may have major shortcomings during mergers or high-speed close encounters between galaxies. During a close encounter between two galaxies, one of the two BHs may select the most bound neighbour particle as a particle belonging to the other galaxy. In this case, at the next time step the BH is suddenly and unphysically relocated to the neighbouring galaxy, thus leaving its original host galaxy without a central BH. In addition, this BH, which has typically moved to an outer region of the galaxy, will be quickly repositioned closer to the centre of the new host galaxy, where another BH is located, within a few time steps. In this way, BH-BH mergers will become faster and more frequent.

To prevent the occurrence of these spurious behaviours, different definitions of the neighbours, over which to search for the most bound particle, have been employed. The original radius was set as the SPH smoothing radius of the BH particle, or some kernel radius associated with a different hydro solver (e.g. Di Matteo et al. 2008; Vogelsberger et al. 2013; Davé et al. 2019). Other authors preferred to search the most bound particle within the BH gravitational softening or a small multiple of it (e.g. Booth & Schaye 2009; Schaye et al. 2014; Ragone-Figueroa et al. 2018) because the smoothing length can be much larger, thus exacerbating the problem mentioned above.

A further condition usually introduced to search for the most bound particle is that its velocity relative to the BH has to be smaller than a threshold in the attempt to ensure that it belongs to the same galaxy. The most commonly used threshold is a fraction of the local sound speed, as originally introduced by Di Matteo et al. (2008). However, Ragone-Figueroa et al. (2018) found that this criterion is not effective, besides having a not-so-clear physical basis (see also Bahé et al. 2022). Their results improved by imposing a maximum velocity of the order of the typical motions of particles within galaxies, namely 100–200 km/s. These values and the smaller search radius limited the unphysical transfer of a BH to another galaxy in their zoom-in massive cluster simulations since the typical orbital velocities of galaxies in clusters are much larger. While including such additional criteria improves the performance of the method, in the simulation presented in this work we adopt a version of the repositioning closer to the original one, where the most bound neighbour particle is selected within the SPH smoothing length, provided that its relative velocity is smaller than 25 % of the local sound speed. The motivations behind this choice lies in the increased resolution of the simulations presented here, which provide more complex merger dynamics, and the purpose to highlight the limitations of this commonly used repositioning technique. However, we keep the BH velocity calculation unchanged, regardless of the specific sub-resolution approach used for BH dynamics.

2.2.2 Boosting the dynamical mass

An alternative approach to account for the limited mass resolution when the mass of the BH particles is close to its seeding mass is to increase the BH dynamical mass at seeding artificially. In this approach, once seeded, two different masses are assigned to the BH: the real mass, m_r, which grows continuously by the Eddington-limited Bondi-like prescription (see Sect. 3.2.2 below), and the dynamical mass, m_d, which enters in the computation of gravitational force. At seeding, the latter is set at a relatively large value, typically equal to the mass of DM particles, while the former is a few orders of magnitude smaller. As long as m_r < m_d, the value of m_d does not increase by gas accretion, which only affects the value of m_r. Once m_r ≥ m_d, the real mass increases in a continuous way, while the dynamical mass increases according to a stochastic prescription for the swallowing of neighbour gas particles (e.g. Springel et al. 2005). From then on, the two masses remain similar, differing only because of the stochastic swallowing of gas particles.

This artificial boost of the BH dynamical mass at seeding is intended to amplify the effect of the numerically resolved DF, eventually preventing the BH from escaping the host halo soon after it is seeded.

Clearly, this method is also prone to spurious effects. For instance, Tremmel et al. (2015) pointed out that initialising the BH mass with a value hundreds of times higher than its real value can affect the mass of the host galaxy, thus unavoidably impacting its subsequent evolution.

In the present section, we review the most commonly used techniques to deal with BH dynamics, mainly to provide a background to the simulations presented in this work. We refer the reader to Di Matteo et al. (2023) for a more comprehensive discussion.

3 Super-massive black holes in cosmological simulations

In this section, we discuss the description of SMBH evolution and of the ensuing AGN feedback in our simulations. In Sect. 3.1, we briefly give an overview of the OpenGadget3 code, within which we implemented our model to correct for unresolved DF. Given the finite force and mass resolution of cosmological simulations, sub-resolution models are needed to describe the processes of birth, accretion, and feedback of SMBHs. Furthermore, the merger events between BH pairs cannot be followed down to the final inspiraling of their orbits. Therefore, we need to include also some criteria to establish when to merge two BHs. In Sect. 3.2 we review the approach to treat BHs in cosmological simulations, by focusing on the seeding criterion (Sect. 3.2.1), on gas accretion (Sect. 3.2.2), and on the conditions allowing BH-BH mergers (Sect. 3.2.3).

3.1 The OpenGadget3 simulation code

Our simulations are based on the OpenGadget3 code (Dolag et al., in prep.; see also Groth et al. 2023), which represents an evolution of the GADGET-3 code (which, in turn, is an improvement of the previous GADGET-2 code by Springel 2005). OpenGadget3 solves gravity with the Tree-PM method (see also Ragagnin et al. 2016). In the simulations presented here, hydrodynamics is described by the SPH formulation presented by Beck et al. (2016), which overcomes several of the limitations of the original SPH formulation of GADGET-3.

OpenGadget3 is parallelised using a hybrid MPI/OpenMP/ OpenACC scheme (Ragagnin et al. 2020). Adopting a limited number of MPI tasks per node allows us to reduce the ‘communication surface’, while efficiently using OpenMP inside a single shared-memory node. Load-balancing is achieved using a domain decomposition based on a space-filling Peano-Hilbert curve, whose fragmentation into segments (each assigned to an MPI task) guarantees a very good computational balance, at the expense of some memory imbalance (Springel et al. 2005).

OpenGadget3 includes a sub-resolution description of a range of astrophysical processes relevant for the simulations presented here: metallicity-dependent radiative cooling (e.g. Wiersma et al. 2009), an effective model for star formation from a sub-resolution description of the multi-phase structure of the interstellar medium (Springel & Hernquist 2003), a model for stellar evolution and chemical enrichment from AGB stars and type-Ia and II supernovae (Tornatore et al. 2007; Bassini et al. 2020), and a model to follow the evolution of SMBHs and the ensuing AGN feedback (see below for details). As for the latter, we remind that the aim of this paper is to present an improved implementation of a sub-resolution description of the effect of DF on the dynamics of BH particles (see Sect. 2 for details).

3.2 Black holes in cosmological simulations

3.2.1 Black hole seeding

In our simulations, BHs are described by collisionless sink particles which are initially seeded within a halo hosting a ‘bona fide’ galaxy. The halo is identified through a Friend-of-Friend (FoF) algorithm (with linking length equal to 0.2 times the mean separation of DM particles). For a BH particle to be seeded, we require the host halo to fulfill few conditions, so as to guarantee that it is well resolved and that star formation already took place within it. Following Hirschmann et al. (2014), we added to the halo mass threshold criterion introduced by Springel et al. (2005) additional conditions for star and gas fraction of the halo.

In detail, in the simulations presented in this work, the following seeding conditions must all be met: (i) the DM mass of the halo exceeds the value of 6.94 × 10¹⁰ M_⊙; (ii) the stellar mass is at least 2 per cent of the total mass and 5 per cent the DM mass of the halo; (iii) the gas mass reaches a value of 10 per cent of the stellar mass; (iv) the halo does not contain any other BH particle. If a halo fulfils these conditions, the most bound star particle of the halo is converted into a BH. The mass of a BH particle at seeding, M_BH,seed, is not fixed, but scales with the amount of stars in the FoF group according to: $$M_{\text{BH},\text{seed}}=M_0\frac{M_{*, \text{h}}}{f_{*}{M}_{\text{DM},\text{seed}}},$$(10)

where M_*,h is the stellar mass assigned to a FoF halo, f_* and M_DM,seed are the input parameters for the fraction of stellar mass and the minimum DM mass for a FoF halo where to seed a BH, respectively. As for M₀ , it is the minimum seeding mass of a BH particle, when the seeding conditions are just met. In the simulations presented in this work, the mass of star particles (see Sect. 4 and Table 2) is larger than M₀ . Therefore, the assumption that the BH mass is larger than the one of the surrounding objects, on which the derivation of the Chandrasekhar DF formula is based, is not met at seeding due to limited mass resolution, and for the above condition to hold BHs need to grow by a large factor. Therefore, it is not surprising that, despite the initial BH position is at the location of the most bound particle, two-body scatterings with neighboring particles can easily cause the BH to be scattered outside its host galaxy.

3.2.2 Gas accretion onto BHs

During its evolution, the mass of a BH increases through two channels: accretion of surrounding gas and merging with other BHs. As for the former, we adopt the accretion model originally implemented by Springel et al. (2005). BH accretion rate is calculated according to Bondi-Hoyle formula (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952) as $${\dot{M}}_{\text{Bondi }}=\frac{4\pi G^2{M}_{\text{BH}}^2\alpha \rho }{{\left(c_s^2+v^2\right)}^{3/2}},$$(11)

where ρ and c_s are the density and the sound speed of the surrounding gas computed at the position of the BH particle, v is the relative velocity of the BH with respect to the surrounding gas particles, and α is a ‘boost’ factor introduced to account for the limited resolution with which gas density in the surroundings of the BHs is reconstructed. Following Steinborn et al. (2015), we distinguish between accretion from the hot (T > 10⁵ K) and the cold gas (T < 10⁵ K), using α = 10 and α = 100 respectively for the hot and the cold gas.

The accretion is always limited to the Eddington accretion rate, $${\dot{M}}_{\text{Edd}}=\frac{4\pi G{M}_{\text{BH}}{m}_p}{{\eta }_r{\sigma }_{T}c},$$(12)

where m_p is the proton mass, σ_T is the Thompson cross-section. The parameter η_r is the radiative efficiency and represents the fraction of the accreted rest-mass energy which is converted in radiation (Novikov & Thorne 1973; Noble et al. 2011). In our simulations, we use η_r = 0.1. In addition, we allowed the black hole to swallow gas particles via stochastic accretion, as originally proposed by Springel et al. (2005). Following Fabjan et al. (2010), gas particles are not swallowed entirely, but are sliced into three parts, so as to have a more continuous description of stochastic accretion. In this way, we can assign to each BH two masses that in general slightly differ: a ‘true’ mass, which grows in a continuous way according to the above Eddington-limited Bondi-like accretion model, and a dynamical mass, which is varied each time that a portion of a gas particle is stochastically selected for the swallowing.

3.2.3 Mergers

The finite force resolution set by the gravitational softening sets the scale below which gravitational interactions, including mergers between BHs, cannot be properly followed. The simplest criterion for defining when a BH-BH merging event occurs is to impose a limiting BH-BH distance, d_merg, below which the two BHs could immediately merge. In our simulations, we adopt a value of d_merg = 5 · ϵ_BH. However, during fly-by encounters between galaxies, it could happen that this criterion produces an unwanted behaviour, with the two BHs forced to merge, even if their relative velocities are large enough to make them gravitationally unbound.

To circumvent this potential problem, we include in our simulations two additional criteria that must be fulfilled for a merger to happen. Following similar arguments as in Di Matteo et al. (2008), we first require that $$v_{\text{rel}}\text{<}0.5c{}_s.$$(13)

Furthermore, we impose that the two BHs are gravitationally bound by requiring: $$\frac{\left|{\text{Φ}}_{{\text{BH}}_1}-{\text{Φ}}_{{\text{BH}}_2}\right|}{a}<0.5\cdot c_s^2-v_{\text{rel}}^2,$$(14)

where ${\text{Φ}}_{{\text{BH}}_1}$ and ${\text{Φ}}_{{\text{BH}}_2}$ are the values of the gravitational potential at the positions of the two merging black holes, c_s is the local sound speed, υ_rel their relative velocity and a is the scale factor, normalised to unity at redshift z = 0. During a merger event involving two BHs, the BH having a lower value of the potential swallows the other BH. Consequently, after the merger, the surviving BH retains its original position and velocity, while its mass increases by the mass of the swallowed BH.

3.2.4 Feedback energy

As a BH accretes gas, it injects energy in the surrounding region and a fraction of the energy radiated during gas accretion is thermally coupled in the form of feedback energy to the surrounding medium (Wurster & Thacker 2013, for an extensive comparison of different implementations of AGN feedback in simulations). Within the simulations presented in this work, we consider a purely thermal mechanism of feedback, whose energy rate is calculated according to $$\dot{E}={\eta }_r{\eta }_f{\dot{M}}_{\text{BH}}{c}^2.$$(15)

In the above expression, Ṁ_BH = min(Ṁ_Bondi, Ṁ_Edd) is the BH mass accretion rate, and η_f is the fraction of the radiated energy that is thermally coupled to the surrounding gas. Furthermore, we emulate a transition between the radio and quasar mode of BH feedback by varying the parameter η_f, as described by Sijacki & Springel (2006) and Fabjan et al. (2010). Whenever the accretion rate of the BH is one-hundredth of the Eddington limit, we increased the fraction of energy thermally coupled with the surrounding gas by a factor of four. We adopt a value of η_f = 0.05 during the quasar mode, increasing to η_f = 0.2 during the radio mode.

For the sake of clarity, we summarise in Table 2 the information on the force resolution and main characteristics of BH model implemented in the simulations performed for this work, which are described in the next section.

4 Simulations

To assess the performance of the novel DF correction and to compare it with other prescriptions (Sect. 2.2), we carried out a set of two zoom-in simulations (a group-sized and a cluster-sized halo), and a cosmological box (with side length of 16 comoving-Mpc), all at the same resolution. This allowed us to test the new implementation in different environments. All simulations are performed assuming a flat ΛCDM cosmology with Ω_m = 0.24, Ω_b = 0.0375 for total matter and baryon density parameters, h = 0.72 for the Hubble parameter, n_s = 0.96 for the primordial spectral index, σ₈ = 0.8 for the power spectrum normalisation. We report in Table 3 a description of the main characteristics of these simulations. The selected zoom-in regions belong to the Dianoga simulation set introduced by Bonafede et al. (2011) and Bassini et al. (2020). We refer to such papers for a detailed description of this set of zoom-in initial conditions for simulations of galaxy clusters, as well as for details on the model of star formation and chemical enrichment adopted. The D9 cluster has a mass M₂₀₀ ≃ 1.53 × 10¹⁴ M_⊙ at z = 0, while the other refers to a lower-mass group of M₂₀₀ ≃ 1.38 × 10¹³ M_⊙, namely Cl 13³. Both objects have been simulated at the same mass and spatial resolution (see Table 2).

As for the cosmological box, namely CosmoBox, it has been simulated by adopting exactly the same cosmological model and the same implementation of all physical processes as the zoom-in simulations. For each initial condition of D9, Cl13, and CosmoBox, we carried out three simulations using the same settings while changing only the sub-resolution technique to cope with the BH dynamics. We compared the repositioning scheme, the adoption of a boosted dynamical mass and the new implementation of DF introduced in Sect. 2. In the following sections, we refer to these schemes as, respectively, REPOS, DYNMASS and DYNFRIC. Comparing the results of these simulations allows us to assess the effect of our new implementation of DF with respect to some previously adopted ad hoc prescriptions to account for this effect.

In the following sections, we present our results using two different approaches. In Sect. 5, we present a statistical analysis to show how the method used to track BH orbits impacts the properties of their population. Section 6 is dedicated to the study of the evolutionary dynamic histories of individual BHs, comparing their path when they sink into the potential well of the host sub-halos and during mergers when governed by different sub-resolution prescriptions.

5 Properties of the BH population

The minimum aim of our DF model is to avoid spurious numerical effects that other methods may produce, as discussed in Sect. 2. Thus, in this section, we present an analysis of the overall properties of the population of SMBHs in our simulation, focusing on the ability of our DF model to: place and hold the BH at the centre of the host galaxies (Sect. 5.1), prevent BH particles from spuriously wandering outside the host galaxies (Sect. 5.2), reproduce the co-evolution of BHs and host galaxies (Sect. 5.4).

5.1 Centring the BHs within host galaxies

To assess the ability of each model to correctly locate the BHs at the centre of the host sub-halos, we select all the main halos and sub-halos having a BH within their DM half-mass radius (hereafter R_HMS). We associated each BH to the closest (sub-)halo, as this procedure allows us to assess how BHs are kept at the centre of their host or are off-centred. Sub-halos are identified by the SubFind algorithm (Springel et al. 2001; Dolag et al. 2009) and are requested to feature at least 20 particles. In the zoom-in simulations, we excluded from the analysis all the halos that contain within their R₂₀₀ at least one low-resolution DM particle spuriously scattered within the high-resolution region. We then calculated the distance between the (sub-)halo centres and the closest BH. The sub-halo centres are identified with the minimum of the gravitational potential occurring among the member particles of that sub-halo. Their values are retrieved from the SubFind output.

Figure 1 shows the histogram of the probability density function (PDF) of the number of (sub-)halos having a BH within their R_HMS, versus the distance ∆r from the associated BH. The sample of sub-halos in Fig. 1 is obtained by summing up all the (sub-)halos from D9, Cl13 and CosmoBox. However, we verified that the individual analysis of each of these simulations keep the overall results unchanged. In each column, we compare different prescriptions for BH dynamics with the DF model, shown in blue, and we report the results at z = 3, 1 and 0, from top to bottom. For reference, the dashed-dot vertical line marks the value of the Plummer-equivalent softening length of BH particles. Not surprisingly, the REPOS prescription predicts sub-halos that host generally well-centred BHs. This is the most consistent advantage of the repositioning technique, as well as the most predictable, since it explicitly forces the BHs to be relocated at each time-step at the centre of their host substructures. As for the DYNMASS, it shows the less pronounced tail for higher distances (∆r > 1 ckpc) at z = 3. In fact, at high redshift, the impact of a large dynamical mass is relatively more important given that a larger number of BHs is expected to have relatively smaller true mass and live in less resolved hosting halos. Interestingly, in this regime, our DF model performs quite well in keeping BHs at the centre of their host halos, without resorting to ad hoc prescription. The comparable performance of the two methods is actually a nontrivial and encouraging result for our DF implementation; in fact, accurate centring of the BHs in sub-halos is obtained as a result of an implementation of the DF correction instead of an artificial increase of BH masses, which have the side effect of altering the structure of the sub-halos gravitational potential, thereby impacting on the resulting galaxy formation process (e.g. Chen et al. 2022). Table 4 reports the percentage of sub-halos with which each simulation contributes to the differential distribution of Fig. 1. The D9 simulation provides the most numerous population of centred sub-halos, with at least 70% of the total distribution.

Fig. 1 Probability density distribution of the distances between sub-halos identified by SubFind and the closest BH particle within the RHMS of each sub-halo of Cl 13, D9 regions, and Cosmobox. The rows show the results obtained at different redshifts: z = 3 (up), z = 1 (central), and z = 0 (bottom). The column report the results comparing DYNFRIC and REPOS on the left, and DYNFRIC and DYNMASS on the right side. We include a dashed-dot line in each plot indicating the softening length of the BH as a reference for the spatial resolution of the simulation.

5.2 The population of wandering BHs

Another crucial aspect to consider when judging the reliability of a DF model is its capability of preventing BHs from escaping the host sub-halos because of two-body interactions with larger-mass particles or during merger events. We denote the population of these kicked-off BHs as wandering BHs. We classify a BH as wandering whenever its distance from the closest sub-halo is larger than twice the RHMs of the closest sub-halo. With this definition of WBHs we aim to focus on those BHs that have been expelled from their parent host galaxies as a consequence of spurious numerical effects. Thus, they should not be confused with the ‘classical’ population of off-centre WBHs as, for instance, investigated by Tremmel et al. (2018), Ricarte et al. (2021), which instead are expected to become wandering as a consequence of numerically resolved gravitational dynamics. It is worth noting that our definition can misidentify wandering BHs during close encounters between sub-halos, when a non-WBH beloging to a sub-halo can temporary become a wander of the closer sub-halo passing by. We verified that limiting the search only to those BHs with a distance exceeding twice the R_HMS of the closer among the two sub-halos, the results reported below do not vary considerably.

Throughout this section, we only present the results from the D9 simulation. Its larger population of WBHs enables a more accurate statistical analysis.

Figure 2 shows the cumulative fraction of wandering BHs to the total number of BHs as a function of BH mass at z = 0 in the left panel and as a function of redshift in the right panel, for all the considered sub-resolution models. In the left panel, for each BH mass on the x-axis we see the cumulative fraction of wandering BHs on the y-axis. We observe that, irrespective of the method employed to keep the BH centred, the ratio of the wandering BHs to the total number of BHs increases as we progressively consider larger BH masses at the numerator. Each curve extends up to the value of the mass of the most massive wandering BH found in the considered simulation. The REPOS simulation, in particular, exhibits the higher fraction of wandering BHs. The percentage of wandering BHs in the simulation using the repositioning scheme rises for massive BHs (M_BH > 10⁷M_⊙). Figure 3 shows the number of WBHs for each mass bin found in the D9 simulation at z = 0. For every simulations, the majority of wandering WBHs’ masses lies below 10⁸. The REPOS simulations shows a peak at M_BH ~ 10⁷M_⊙. Interestingly, this implementation shows a second lower peak for 10⁷ < M_BH < 10⁸. Again, this suggests that the repositioning technique adopted in this work can produce a significant displacement or BH ‘teleporting’ of even massive BHs. The peak of the distribution in the DYNMASS simulation is shifted toward larger masses compared to REPOS. Moreover, DYNMASS reproduces lower WBHs for M_BH < 10⁷ compared to DYNFRIC and more for M_BH > 10⁷. On one hand, BHs in the low-mass range are boosted in mass, preventing them from escaping the surrounding sub-halo. On the other, once their mass grows up to M_DM, no correction acts on the dynamics of the BHs. Using the DYNFRIC technique, the peak of the distribution flattens compared to the other simulations. The number of WBHs with M_BH < 10⁷ is slightly higher than for REPOS and DYNMASS, but it reduces in the high-mass end.

The increase of wandering BHs using the repositioning mechanism is surprising: these massive BH should reside in sub-halos having a deeper potential well. Therefore, they should be less prone to wander outside their host galaxy. The origin for the excess of wandering BHs in the REPOS scheme can be ascribed both to close encounters between sub-structures and, at decreasing redshift, to the presence of large-scale potential gradients. To illustrate how these mechanisms work, we show several wandering BHs in the two panels of Fig. 4.

Each panel of this figure shows, on the central plot, the projected stellar density around two wandering BHs, BH 2 in the top panel, and BH 4 in the bottom panel. We also show the position of the other BHs in the field, marking wandering BHs with cyan crosses and the BHs centred on their host as green crosses. Furthermore, the centres of neighboring sub-halos are denoted by shaded white dots, with circles corresponding to the sub-halo R_HMS enclosing them. The side plots display the gravitational potential Φ of the stars along the x-axis (on the top) and y-axis (on the right). In the top panel of Fig. 4, BH 2 wanders due to a close encounter between two sub-halos, one having a deeper potential well. BH 2 then migrates to the location of the most bound neighbouring particles, thus moving along the gradient of the potential of the central sub-halo, and leaving its initial host sub-halo. It is worth noticing that, as defined here, a BH lying within twice the R_HMS of a sub-halo can be a wander: our definition classifies as not wandering only those BHs located within twice the R_HMS of the closest sub-halo. The bottom panel of Fig. 4 shows the position at z = 0.06 of the supermassive wandering BH 4 with a mass of 10¹⁰ M_⊙. Being embedded in the large-scale potential gradient of the most massive halo, depicted in the side plots of the panel, BH 4 escaped from its original sub-halo. This demonstrates that for the repositioning scheme implemented in this work, even the most massive BHs can move away from their host galaxies, thereby completely changing their gas accretion and the ensuing release of feedback energy.

In the simulations using the DF correction and DYNMASS, wandering BHs instead originate from two-body or even three-body scattering. When the BH is not repositioned, the dynamics during these events can be quite complex; this aspect is analysed in more detail in Sect. 6.

Finally, Fig. 5 displays the position of all the BHs in a projected stellar density map within 1 Mpc from the centre of the most massive halo in the D9 simulation at redshift z = 0 for REPOS (top panel), DYNMASS (bottom left panel) and DYNFRIC (bottom right panel). We plot BHs with crosses of different colours depending on their distance ∆r from the closest sub-halo, comparing it to the R_HMS of this sub-halo. Dark-blue crosses indicate BHs with ∆r < R_HMS, green crosses are for R_HMS ≤ ∆r < 2 × R_HMS, and cyan crosses mark the wandering BHs, i.e. those with ∆r ≥ 2 × R_HMS. We also plot the circles for the ten most massive sub-halos in the region, each having radius that is suitably scaled according to the value of its R_HMS. Each circle is centred on the sub-halo position and has a radius proportional to its R_HMS. The radius of the yellow circle corresponds to the actual physical value of the R_HMS of the BCG. For the purpose of readability, we adopted a different scaling for marking the size of the other sub-halos, as indicated in the figure legend.

Figure 5 shows that the DYNFRIC simulation exhibits a less numerous population of wandering BHs compared to DYNMASS (8 WBHs against 13 within the main halo R_HMS), mainly located close to multiple BH systems. Most of the DYNMASS wandering BHs occupy the central region of the halo. However, crucial differences arise between REPOS and the outcomes of the other simulations: the overall population of BHs using the repositioning scheme, significantly decreases both in the core and in the outskirt of the halo. Most of the sub-halos, identified by stellar density peaks or with circles for the more massive ones, lacks a central BH.

Fig. 2 Cumulative number of wandering BHs to the total number of BHs found in the D9 region as a function of the BH mass at z = 0 (left panel), and as a function of redshift for z = 3, z = 1, and z = 0 (right panel). In the left panel, no mass threshold is adopted to select the BHs, and hence all the BHs with mass above the seeding mass are included (see Table 2). In both panels, the simulation using the DF model, marked with a blue line, is compared with REPOS in red and DYNMASS in green. The wandering BHs are defined as those having a distance larger than two times the RHMS from the closest sub-halo.

Fig. 3 Histogram of the number of WBHs as a function of their mass for the D9 simulation at z = 0. We overplot the distribution of DYNFRIC and REPOS on the left and DYNFRIC and DYNMASS on the right.

5.3 Merger events and multiple BH systems

The occurrence and timing of merger events are highly sensitive to the prescription to control the BH dynamics. We observe that the REPOS simulations facilitate mergers when two BHs approach a distance d_merg (see Table 2), making them merge on shorter timescales than when the other prescriptions are adopted (see Sect. 6).

Figure 6 displays the cumulative (upper part) and the differential (lower part) distributions of the number of merger events as a function of the redshift, for the three simulations. DYNFRIC and DYNMASS simulations predict similar results, with DYNFRIC showing slightly more mergers in denser regions (Cl 13 and D9) and less mergers in the CosmoBox simulation. By contrast, all the simulations based on REPOS feature a significantly higher number of mergers, more than twice the number of mergers of both DYNFRIC and DYNMASS in the zoom-in simulations and more than three times in the CosmoBox simulation. Interestingly, the density peak of mergers for the repositioning model occurs at z ≃ 1.5 in all the simulations.

Before the merging events occur, one would expect to find structures hosting systems of two or more BHs. However, the increase of mergers in the repositioning scheme is not counterbalanced by an increase of multiple systems, i.e. of structures containing two or more BHs.

Figure 7 shows the ratio between the number of sub-halos hosting more than one BH within the R_HMS to total number of sub-halos hosting at least one BH for z = 3, 1, 0 in the D9 simulations. While the percentage of sub-halos hosting multiple BHs reaches more than the 8% in DYNFRIC (blue) and DYNMASS (green) simulation, the REPOS shows a consistently lower percentage. The concurrence of a larger number of mergers and the nearly absence of multiple systems found for the REPOS simulations can be ascribed to extremely fast mergers, with multiple BHs coexisting within the same halo only for a rather short time.

We note that the number of merger events is influenced both by the adopted model to trace the BH dynamics and by the specific seeding prescription adopted. In principle, a higher number of merging events may be associated with a more frequent seeding. This situation may arise when, following a merging event between BHs but not between their corresponding halos, one of the two halos remains orphan of its BH, while still matching the seeding conditions (see Table 2). In this case, a new BH would be seeded at the centre of the orphan halo, thereby possibly contributing to increase the BH-BH merger rate. To quantify how this ‘seeding bias’ affects the merging predictions of from each model, we report in Table 5 the number of seeded BHs in each simulation.

We observe that the increase in the total number of seeded BHs in the REPOS simulations is only marginal. The absence of a proportional rise in the number of seeding black holes despite the increase in mergers in the REPOS simulations can be addressed on the concurrence of two particular seeding conditions.

On the one hand, the seeding occurs exclusively within the main halos identified by the FoF halo finder, and not within the sub-halos identified by SubFind. On the other hand, BHs are not seeded in FoF groups which already contain a BH particle. Instead, the more frequent scenario in the REPOS simulations is the merging of two BHs initially belonging to two sub-structures contained within the same FoF halo. In this case, no further seeding takes place within the sub-halo which eventually remain orphan of its BH.

Moreover, the simulations adopting the dynamical mass implementation predict less mergers in the D9 and Cl13 simulations (Fig. 6) and a slightly higher fraction of multiple BH systems at z = 0 (Fig. 7) compared to DYNFRIC, even if the latter exhibits a higher number of seeded BHs (Table 5). This circumstances already suggest that the dynamical mass model fails to reproduce some merger events, something that will be explored in detail in Sect. 6.

Fig. 4 Stellar density projection along the z-axis with depths of 100 kpc (top) and 300 kpc (bottom), centred on the position of a wandering BH at redshiſt z = 3 (tagged as BH 2; left panel), and of a wandering BH at z = 0.06 (tagged as BH 4; right panel) in the D9 simulation using the REPOS scheme. In both panels, the cyan crosses identify the wandering BHs, while the green crosses identify the BHs centred in their host. We show sub-halo centres as white shaded dots, with the white circles indicating the corresponding RHMS. On the top and on the right of each panel we show the gravitational potential Φ of all the star particles (light white dots) and of the BHs (crosses) in the field, projected along the two orthogonal directions.

Fig. 5 Stellar density maps along the z-axis with a depth of 1 Mpc centred on the most massive halo of D9 at z = 0 in the REPOS (top), DYNMASS (bottom left), and DYNFRIC (bottom right) simulations. The panels are all 1 Mpc on a side. In each panel, we plot with dark-blue crosses the BHs located within the RHMS of the associated sub-halo (the same criterion adopted in Sect. 5.1). BHs lying between RHMS and 2RHMS of the closest sub-halo are indicated with green crosses, while wandering BHs (defined as those located beyond 2RHMS of the closest sub-halo) are shown as light-blue crosses. The values of RHMS of the ten most massive structures correspond to the radii of the circles, each centre on the position of the corresponding sub-halo. The RHMS of the BCG, in yellow, corresponds to the physical size of the yellow circle. The legend in the upper left panel of the plot shows the scaling size of the other sub-structures in the region, marked in white.

Fig. 6 Cumulative and differential distributions of the number of merger events as a function of redshift. Cl13, D9, and CosmoBox are in the first, second, and third rows, respectively. The REPOS simulation results are marked with a dash-dotted red line (with a shaded area marking the differential distributions), the DYNFRIC simulations with a blue solid line and the DYNMASS simulations with a dashed green line. The curves stop at the redshift corresponding to the occurrence of last merger event.

Fig. 7 Ratio between the number of sub-halos hosting more than one BH and the total number of sub-halos hosting at least one BH in D9 simulations. We consider DYNMASS (dashed green line), DYNFRIC (solid blue line) and REPOS (dot-dashed red line) at z = 3, 1, 0.

5.4 The M_*−M_BH relation

A first diagnostic to investigate how the processes of star formation and galaxy evolution are intertwined with the evolution of the population of SMBHs is to look at the relationship between SMBH masses and stellar masses of the bulges of the host galaxies, the so-called Kormendy-Magorrian relation (Kormendy et al. 1993; Magorrian et al. 1998). Figure 8 shows this relation obtained by varying the implementation for the BH dynamics in all the simulations for the three initial conditions considered in our analysis. In detail, for each halo with a mass larger than a threshold value M_thr = 10¹⁰ M_⊙, we distinguish between the central galaxy, to be identified with BCG or Brightest Group Galaxy (BGG; marked with diamond symbols in Fig. 8), and the satellite galaxies, which are hosted within the substructures identified by SubFind (marked with circles). Since the SubFind algorithm does not split the diffuse stellar component of the Intra-Cluster light to the one bound to the BCG, at least for large-mass halos we adopted a fixed physical aperture within which measuring the BCG/BGG stellar masses (see for instance Marini et al. 2020). A fixed aperture is also commonly adopted for a fair comparison with observational measurements of BCG/BGG stellar masses (Bassini et al. 2019; Ragone-Figueroa et al. 2018). For the BCGs, we calculate the stellar mass of the stars belonging to the BCG/BGG within 70 kpc from the centre. As for the satellite galaxies, their stellar mass is computed by considering all the star particles that SubFind assigns to the substructure. For the central galaxies, we associate the most massive BH within 70 kpc from the centre of the structure. For the satellite galaxies, we link the closest BH within their R_HMS . Symbols are colour- coded according to the distance between the galaxy centres and the associated BHs. For reference, in each panel we also show the relation obtained by McConnell & Ma (2013) and the data for BCGs and BGGs as measured by Gaspari et al. (2019). We note in general that all the simulations reproduce rather well the observational relation, possibly with a slightly larger normalisation. Besides proving the good calibration of the AGN feedback parameters, this agreement also demonstrates that this scaling relation is relatively insensitive to the details of the model adopted to account for unresolved DF on BH dynamics. Table 6 presents data on stellar mass, associated BH mass, and the distance of the BH from the halo centre for the most massive halo within each simulation. The DYNFRIC simulations predict better centred BHs compared to the DYNMASS simulations, although all the distance values are below the resolution limit of the simulations.

As for Cl13, we note that the DYNFRIC, REPOS and DYNMASS simulations all produce very similar results. In all the three cases, the stellar mass of the BGG, and the mass of the hosted SMBHs are quite similar at z = 0.

In the D9 region, the larger statistics of sub-halos helps to better understand what happens in different scenarios. Still, the DYNFRIC and the DYNMASS simulations produce comparable results, with the DYNFRIC simulation even further reducing the offset of the most massive BHs from the centre of their host galaxies. The results for the D9 region using REPOS are in good agreement with the observational data for M_* < 10¹¹ M_⊙.. However, this simulation predicts BHs that are more massive compared to the other implementations, again due to excess of merging episodes, as already discussed. We note that the BH having the highest mass in the D9 region results from the merging between the BH 3 and the massive wandering BH 4 in Fig. 4 that reached the centre of the BCG from z = 0.06 to z = 0.

Finally, the CosmoBox results confirm the substantial agreement between our predictions and observations. The DF model further demonstrates its increased efficiency in centring the BHs compared to DYNMASS. Nonetheless, both simulations reveal three massive sub-halos with a stellar mass > 10^10.5 M_⊙ , whose BHs have a significant displacement. The visual representation of these three occurrences is depicted in Fig. 9, which shows the projection of the stellar density centred around these pathological sub-halos. In the first row, two analogous situations are presented: a close encounter between two substructures hosting BHs of different masses. The displacement obtained in Fig. 8 follows from the association of the most massive BH to the central halo, marked in Fig. 9 with a red circle. This feature, rather than being a real off-centring of the BH, is a consequence of the method that we adopted to associate a BH to a given structure.

The third sub-halo, shown in the bottom right panel of Fig. 9, hosts a truly off-centre BH. This is the outcome of a merger event between two substructures which took place at redshift z = 0.19 (see bottom left panel). Since then, the BH belonging to the merging substructure did not have the time to sink to the centre of the merged sub-halos, hence the displacement at z = 0 between the centre identified by SubFind and the BH.

Fig. 8 Relationship between BH mass and stellar mass of the host galaxies in the Cl13 (upper panels), D9 (central panels), and CBox simulations (lower panels). Each column contains the results obtained with a different BH dynamics prescription: DYNFRIC (left), DYNMASS (central), REPOS (right). The diamonds refer to the BHs associated with the main halos while circles correspond to sub-halos belonging to halos more massive than 1010 M⊙ Diamonds and circles are colour-coded according to the distance between the sub-halo centres and the associated BHs. For comparison, we also plot observational data from Gaspari et al. (2019) with black dots with errorbars, and the relation obtained by McConnell & Ma (2013) with a green and shaded area.

6 Analysis of individual events

Besides the analysis of the properties that emerge from the population of BHs in each simulation, we perform a detailed study of the effect of different sub-resolution prescriptions on the dynamics of single BHs.

To better understand the BH dynamics, due to the different prescriptions to follow it, we adopt the following approach. We freeze the configuration of BHs and of their host galaxies at two snapshots from the Cl 13-DYNMASS simulation at z = 3, and z = 1.26. Using these snapshots as initial conditions, we run simulations using: the DYNFRIC model introduced in Sect. 2.1, the DYNMASS method based on boosting dynamical mass at seeding (see Sect. 2.2.2), the REPOS method based on relocating the BH on the local potential minimum (see Sect. 2.2.1), and a fourth simulation not correcting to the BH position is introduced (NOCORR in the following). The simulations having as initial conditions the snapshot at z = 3 evolve until z = 2, while the ones starting at z = 1.3 reach z = 0.95. In this way, we trace the histories of BHs that, at the beginning of the simulations have the same position, mass and velocity, and are located in substructures with the same characteristics. Therefore, any difference between their subsequent orbits is purely driven by the different tracing methods adopted.

To ensure that the results were reproducible and marginally affected by the possible chaotic nature of a simulation, we carried out each of them twice, obtaining results which show marginal differences in the timings, while leaving the general results qualitatively unchanged. Clearly, restarting DYNFRIC and REPOS simulations from an output produced at a given intermediate redshift by the DYNMASS run drives to a structure evolution and dynamics that are different from those produced by using DYNFRIC and REPOS since the initial redshift, as described in the previous sections. In particular, our analysis focuses on the evolution of binary BH systems, typically characterised by a massive BH located at the centre of a substructure, indicated as BH_cen, and a second displaced BH, the ‘satellite’, labelled as BH_sat. We focus in the following on three events which represent three very different values of the central-to-satellite mass ratio between the BHs, f_m = m_BH,cen/m_BH,sat; they are labelled as Event 1, Event 2 and Event 3 in the following. Event 1 consists of a binary system of two BHs initially displaced by 4 kpc from each other. The initial BH mass ratio is f_m = 1.1, and both the BHs have a mass smaller than m_DM (see Table 7). For that reason, both masses are boosted in the DYNMASS simulation. Event 2 involves a binary system of two BHs, initially at a distance of 1.5 kpc, with a higher mass ratio (f_m = 50) compared to Event 1, but at the same redshift. Still, the two masses are both increased in the DYNMASS run. Event 3 consists of a binary system of two BHs initially separated by 10 kpc atz = 1.26. The BHs have a large mass ratio f_m = 313. This time, only BH_sat, whose mass is less than m_DM, undergoes the dynamical mass correction in the DYNMASS simulation.

The characteristics of these three events are summarised in Table 7. Figure 10 provides a graphical representation of the stellar density map at the redshift reported in the legend. The figure also includes the trajectory of the BH_sat, colour-coded according to redshift, within the substructure hosting that event. Each row displays a single event, and results from the DYNFRIC, DYNMASS, REPOS and NOCORR simulations are presented from left to right in each column.

Figure 11 displays on the left column, for the different events in each row, the distance between BH_cen and BH_sat, namely ∆r, during the event. The dashed green line, the solid blue line, the dash-dotted red line and the short-dashed orange line indicate the DYNMASS, DYNFRIC, REPOS and NOCORR runs, respectively. The horizontal grey solid line represents the distance threshold which is necessary for the merger event d_merg to happen (see Table 2). The right side of Fig. 11 focuses on the results of the DYNFRIC runs. In particular, we plot ∆r in the top panel, and the ratio between the DF force and the gravitational force (including the contribution of the DF correction), both for BH_cen (light-blue) and BH_sat (light green) in the bottom panel. In the following subsections, we study each event separately.

Fig. 9 Stellar density maps along the z-axis with depth of 350 kpc (upper row) and 260 kpc (bottom row) of the massive sub-halos having a large separation between sub-halo centre and the associated BH in the CosmoBox simulation, using the DYNFRIC model (blue diamonds in the corresponding plot in Fig. 8, for log10(M*/M⊙) > 10), at redshift z = 0. BHs are marked by green crosses. Upper panels: two close encounters between structures. The halo centre is identified as a red circle. The arrows link to the values of the BH masses (in units as specified in the label). In both the left and right panels, the most massive BH belongs to the off-centre substructure. Bottom panels: a merger event between two substructures at redshift z = 0.19 (left panel) and the resulting merged halo at z = 0 (right panel). The plot on the right side captures an off- centre BH while sinking toward the halo centre.

6.1 Event 1

Looking at the upper row of Fig. 11, the simulation adopting the DYNFRIC model exhibits oscillations with gradually decreasing amplitude and gently drives the two BHs toward the merger. In the DYNMASS simulation, instead, the distance between the BHs exhibits persistent oscillations, which do not decrease in amplitude. Rather than driving the BHs to form a close pair, the enhanced dynamical mass intensifies their mutual gravitational attraction, causing collisions resulting in sustained relative distance fluctuations. The simulation that does not employ any correction shows, after an initial gradual decrease of the distance, a second phase during which the two BHs keep oscillating with respect to each other, with a nearly constant amplitude. The merger between the two BHs, which occurs right after z ~ 2.5, is preceded by a sudden decrease in distance. The inset in the upper right side zooms on the results of the REPOS simulation. We observe a merger event occurring very rapidly, with the distance between the two BHs decreasing through discrete ‘jumps’ as large as almost 1 kpc.

The figure shows that the two BHs can be closer than d_merg , yet without merging. We remind that, this represents a necessary, but not sufficient condition for the merger to happen. Indeed, merging requires the BHs to fulfill all the conditions listed in Table 2. Whenever they are closer than d_merg , they merge if their relative velocity is less than half of their surrounding sound speed and if the gravitational binding criterion defined in eq. (14) is also satisfied. From the visual representation of Fig. 10 we can infer that DYNFRIC, DYNMASS and NOCORR show similar paths for BH_sat until it first crosses the denser region of the sub-halo. Then, DYNFRIC and NOCORR simulations bound BH_sat in the core of the host. Besides, DYNMASS reproduces constant oscillation at rather fixed apocentric distances. A detailed analysis of the reason why the merger does not occur in the DYNMASS simulation shows that the merger failure is due to the enhanced relative velocity and to the difference of the gravitational potential between BH_sat and BH_cen . While crossing the distance threshold imposed by the merging criteria at the pericentre, the higher relative velocity and the larger potential difference compared to DYNFRIC prevent both merging conditions from being satisfied. The sound speed surrounding the BHs remains comparable in the two simulations. These considerations are still valid to explain why the DYNMASS simulation fails to reproduce the merger in Event 2 and Event 3.

On the other hand, the path of BH_sat in the REPOS simulation is discontinuous and short compared to the other: through successive repositioning steps, the BH drops at the centre of the sub-halo. The right panel of the plot highlights the features arising from the application of the DF correction in the DYNFRIC model. In this case, the DF force responds to a local increase of stellar density. The ratio between the gravitational and the DF forces oscillates for both BHs with oscillation of comparable amplitude, indicating that they are orbiting around a central denser zone. Indeed, looking at the upper row of Fig. 10 we observe that the final merger will take place in the central, denser region on the substructure.

6.2 Event 2

The second rows of Fig. 10 and of Fig. 11 show the evolution of the distances during the Event 2. Again, throughout the event, REPOS and DYNFRIC lead to immediate and gradual BH coalescences of the BHs, respectively.

DYNMASS and NOCORR do not produce any merger, thus leaving BH_sat swing around the central one. Thus, neither an ad hoc increase of the BH mass nor the N-body gravity solver without any correction is sufficient to dump the oscillations of the smaller BH. Furthermore, the visual representation in Fig. 10 demonstrate that DYNFRIC, DYNMASS and NOCORR simulations all produce an initial tilt of the orbit of BH_sat in the direction of the elongated denser region of the sub-halo. Then, in the DYNFRIC simulation, the contribution of DF leads the BH toward the centre. The right panel in Fig. 11 illustrates the F_DF /F_grav trend, showing that for this Event 2 the peaks of the contribution of the DF on BH_sat usually take place in correspondence of the minimum distance between the two BHs: as BH_sat approaches BH_cen , it crosses the central region where the density of stellar particles is higher, thus leading to an increase of the DF acting on it. On the other hand, the DF acting on BH_cen is relatively stronger during all the duration of the event, consistent with the more central position that it occupies.

Fig. 10 Stellar density map along the z-axis with depth of 10 kpc in the top and central row and 70 kpc in the bottom row, centred on the central BH (BHcen) of Event 1 in the top row, Event 2 in the central row and Event 3 in the bottom row. The maps of stellar density refer to the redshift indicated on the left panel of each row. In each panel, the curve show the orbit of the satellite BH (BHsat) involved in the merger event, colour-coded according to the redshift. The crosses indicate the last position of the satellite BH before either the merging event or the end of the simulations (see Sect. 6). The columns display the BH evolution using different sub-resolution prescriptions, from the left: using DF (first column), using a boosted dynamical mass (second), adopting the repositioning scheme (third) and finally without any correction for BH dynamics (fourth).

6.3 Event 3

Finally, the bottom rows of Fig. 10 and of Fig. 11 show the evolution during the Event 3. In this case, at the lower redshift reached by restarting the simulations at z = 1.26 to z = 0.96, the mass of BH_cen is more than 300 times higher than that of BH_sat. The latter is initially revolving around the central BH with large oscillations of approximately 30 kpc in amplitude. Once again, REPOS allows the two BHs, initially displaced by 10 kpc, to merge nearly instantaneously, with a ‘long-range teleporting’ of BH_sat from the outskirt to the centre of the substructure, reported in the topright inset in the bottom left panel of Fig. 10. Numerical details of the repositioning scheme can affect the amplitude of the jump: its size could be smaller should different or additional constraints be taken into account (e.g. Ragone-Figueroa et al. 2018, as also discussed in Sect. 2.2.1).

Distance oscillation in DYNFRIC, DYNMASS, and NOCORR nearly coincide with each other. However, while DYNFRIC and NOCORR enable the merger to occur after a short time, this does not happen for DYNMASS. The reason for this lies in the different ways in which these simulations match the merging conditions: the distance threshold criterion is satisfied by all three simulations, but the other merging conditions remains unsatisfied for the DYNMASS simulation. BH_sat, having its dynamical mass boosted, keeps revolving around BH_cen performing wide oscillations of 30 kpc in amplitude. In any case, the similarity of the oscillations described by DYNFRIC and NOCORR suggests that the DF correction in the former should be subdominant in this case. The orbits of BHsat further demonstrate that DYNFRIC, DYNMASS and NOCORR produce almost identical results until z ≃ 1.1. In the DYNMASS simulation, the subsequent oscillations of the BH seem to be results of ‘kicks’ that BH_sat receives while crossing the denser core of the sub-halo. In the right panel, we note that F_df /F_grav for BH_sat oscillates in phase with the distance between the BHs. We verified that this feature is driven by wide oscillations of the gravitational force, which increases when approaching the galactic core where BH_cen resides. Furthermore, the value of F_df/F_grav for BH_cen is significantly higher than for BH_sat . This happens both because of the higher mass of the central BH and due to the strong DF correction due to the concentration of stellar particles in the central sub-halo region.

In summary, these three events highlight the importance of adding a correction to the gravitational force onto BHs contributed by the unresolved DF, as demonstrated in Sect. 6. Furthermore, the description of such examples of merging events also reveals the significant limitations in correctly describing the dynamics of a system of two BHs of both the REPOS (in the worst-case scenario considered here) and the DYNMASS schemes. Indeed, our simulations based on REPOS generate ‘jumps’ of the BHs, while DYNMASS, by artificially boosting the BH mass, produces two-body scattering effects, preventing the formation of a bound system. Lastly, we point out that using an additional constraint on the relative velocity between the BH and the selected neighbor particles in the REPOS model (see Sect. 2.2.1) does not sensibly change merging and sinking timescales when the two BHs are located within the same sub-halo, as in Event 1 and Event 2.

Fig. 11 Left column: evolution of the distance between the central and the satellite BH in Event 1 (top), Event 2 (central), Event 3 (bottom). The dashed green line displays the evolution using the dynamical mass scheme (DYNMASS), the blue solid line is for the DF scheme (DYNFRIC), the dash-dotted red line for the repositioning (REPOS), and the densely dashed orange line is for the case without any sub-resolution prescription (NOCORR). The horizontal grey solid line represents the distance threshold which is necessary for the merger event dmerg to happen. Right column: evolution of the distance between the BHs for the DYNFRIC case (top panel) and ratio between the DF and gravitational forces during each event.

7 Conclusions

We introduce a novel approach to partially correcting for the effect of unresolved dynamical friction (DF) on the orbits of black hole (BH) particles in cosmological hydrodynamical simulations. We implement this correction using the OpenGadget3 code. The main motivation of this study is the need for a more physically motivated method for preventing BH particles in simulations from spuriously leaving the host galaxies as a result of finite numerical resolution. The second aim of this study is to improve the description of BH-BH mergers. The DF correction implemented here has been extensively tested both in zoomedin and fully cosmological simulations. Specifically, we ran two zoomed-in simulations, one for a group-sized halo and another for a cluster-sized halo, and a cosmological volume with a box size of (16 comoving Mpc)³ .

We assessed the performance of the new model by carrying out each simulation with other prescriptions for the BH dynamics: the repositioning scheme and the adoption of a large dynamical mass to reinforce the unresolved DF contribution. To distinguish the effects of variations in the prescription of BHs dynamics, all simulations were run with identical resolution and sub-resolution physics. We summarise the main conclusions of our analysis as follows:

• Offset: The simulations using the continuous repositioning of BH particles on the local potential minimum (REPOS scheme) exhibit the smallest offset between the BH and the centre of the host sub-halo. Our model to correct unresolved DF (DYNFRIC scheme) provides an accurate centring, even outperforming, at z < 1, the scheme based on boosting the dynamical mass of BH particles at seeding (DYNMASS scheme). See Sect. 5.1 and Fig. 1.

• Wandering BHs: The DYNFRIC simulations produce the smallest population of wandering BHs, which are predominantly found in multiple BH systems. The REPOS prescription features the most numerous and massive set of wandering BHs, as close encounters between galaxies and large-scale potential environments favour the spurious repositioning of BHs outside their host sub-halos. See Sect. 5.2, Fig. 2, Fig. 4, and Fig. 5.

• Mergers: The overestimate of merger events in the REPOS simulations in denser environment leads to an excess of galaxies deprived of a central BH, an effect that is more pronounced in the denser environment of a massive galaxy cluster. See Sect. 3.2.3, Fig. 6, and Fig. 5.

• M_* − M_BH relation: The good agreement with the observational data, that is, the ability of all the simulations to reproduce the observed M_* − M_BH scaling relation, demonstrates a relative insensitivity of this diagnostic to the particular prescription adopted, while it is sensitive to the choice of the parameters regulating BH accretion and the ensuing AGN feedback efficiency. See Sect. 5.4 and Fig. 8.

To delve deeper into the details of how different prescriptions affect the orbits of BH particles, we focused on specific BH– BH interaction events. To disentangle the possible diversity in the substructure evolution between different simulations, we restarted the simulations from two snapshots of the group-sized halo at z = 3 and at z = 1.26 to explore how BHs respond to different methodologies governing their dynamics in the same environment. In this analysis, we also carried out simulations without any sub-resolution prescription to correct the orbits of BH particles. The results of this further analysis, introduced in Sect. 6, Fig. 10, and Fig. 11, can be summarised as follows:

• DF: Our novel model, adding the DF correction, predicts damped orbits of the satellite BHs, which gradually approach the centre of a host galaxy, and eventually form a close BH–BH pair.

• Dynamical mass: The large BH dynamical masses in the DYNMASS scheme can lead to spuriously strong interactions between BHs, which could delay or even prevent merger events. Moreover, such interactions eventually inhibit the BHs from satisfying the gravitational boundedness criteria for merging. The overall results, in terms of the quality of the description of orbital decay in BH pairs, can be worse than those obtained without any prescription to correct BH dynamics for unresolved DF.

• Repositioning: Our implementation of the REPOS method allows extremely rapid mergers preceded by large and sudden movements of the satellite BH, which promptly reaches the central region of the sub-halo with a few ‘jumps’; these jumps can eventually span several kiloparsecs.

In summary, an extensive analysis of the BH population arising in different simulations demonstrates that our novel implementation of the correction for unresolved DF force acting on BH particles introduced in this paper provides a robust and reliable description of the DF exerted on BHs by their surroundings. This model achieves at least the same performance as other ad hoc numerical prescriptions in terms of centring BHs in their host halos, while significantly reducing the population of wandering BHs, and overcoming the limitations of this prescription in terms of its ability to describe BH-BH merger events. Thanks to the recent work by Sala et al. (2024), the OpenGadget3 code is now equipped with a spin-evolution model that, coupled to the DF correction, is now supported by a more precise description of the BH dynamics. The compelling estimate of the sinking timescales of BHs onto host halos (see Appendix B) arising from the application of the DF correction proposed in this work merits further investigation (see Taffoni et al. 2003), which we will present in an upcoming work (Damiano et al. in prep.). Such an analysis is crucial for further studies of BH merger rates based on cosmological hydrodynamical simulations (e.g. DeGraf et al. 2024), and will serve as a powerful tool for fully exploiting the potential of gravitational wave astrophysics. Furthermore, examining the interaction between BH dynamics and galaxy evolution may offer fresh insight. As discussed in the present study, adding a DF correction alters both the demography of the BH population and their merger timescales, which have an impact on where and how they interact with the surrounding gas through accretion and heating.

Acknowledgements

We are grateful to the referees for their thoughtful feedback, which has significantly improved our work. We thank Gian Luigi Granato for a number of useful discussions and comments on this paper. Simulations have been carried out at the CINECA Supercomputing Center (Bologna, Italy), with computing time assigned through ISCRA-B calls and through CINECA- INAF and CINECA-UNITS agreements, and at the computing center of INAF – Astronomical Observatory of Trieste (Bertocco et al. 2020; Taffoni et al. 2020). We acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support. This paper is supported by: the Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2, within the activities of Spoke 3, Astrophysics and Cosmos Observations; by the PRIN 2022 PNRR project (202259YAF) “Space-based cosmology with Euclid: the role of High-Performance Computing”, and by a INAF Grant within the “Astrofisica Fondamentale” funding scheme. We acknowledge partial financial support from the INFN Indark Grant and by the INAF project “CONNECTIONS” (COllaboratioN oN codE development for future Cosmological simulaTIONS). Klaus Dolag acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme grant agreement ERC-2019-AdG 882679 as well as support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311 We acknowledge the European Union’s HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Sklodowska-Curie grant agreement number 101086388 – Project acronym: LACEGAL.

Appendix A Velocity distribution surrounding the BHs

One of the hypotheses on which our model to account for unresolved DF is based, is that the velocity distribution function in the vicinity of a BH particle can be expressed as in Eq. (5). $$\text{f}(\text{v})=\frac{3}{4\pi {ϵ}_{\text{BH}}^3}{\displaystyle \sum _i^{\text{N}\left(<{ϵ}_{\text{BH}}\right)}\delta }\left(\text{v}-{\text{v}}_{\text{BH}}\right)$$(A.1)

meaning that the local velocity distribution is proportional to the PDF within the softening length. The possibility to recover, through this formula, the velocity distribution surrounding the BH strongly depends on density and resolution. To address the question of how many particles are necessary to accurately sample the background velocity distribution, we investigated whether we are able to recover the velocity distribution of the particles surrounding the BHs by considering only the particles within the BH softening length. We report here few examples of typical frequent conditions that we encountered. From snapshots of the D9 and Cl13 simulations, by visual inspection we identified BHs whose host galaxies were not undergoing mergers or tidal disruptions. For those cases, we verified that a Maxwellian distribution of all the stars belonging to such galaxies properly describes the velocity distribution. Then we selected star particles within different distances from each BH, and compared the distribution of the relative velocities between star particles and BHs. We considered the following limiting distances: 30 kpc, 5 ϵ_BH, 2 ϵ_BH and ϵ_BH. Table A.1 reports the BH mass and the stellar mass enclosed within 30 kpc from the halo centre for the examples we describe below.

Fig. A.1 Probability distribution functions of the relative velocity between a massive BH (see A.1) in the Cl13 simulation using DF at redshift z = 2.5, and the surrounding stars within 30 kpc (enclosing all the stars in the sub-halo) (red line), 5 ϵBH (yellow line), 2 ϵBH (green line) and ϵBH (blue line). In the left panel, we compare the best fitting Maxwellian PDFs for these four cases. The right plots superimpose the PDF measured from the simulation with the Maxwellian fit for the 30 kpc case (top) and the ϵBH case (bottom).

Figure A.1 compares the velocity probability distribution function (PDF) for a selected BH hosted at the centre of a galaxy having more than 16000 star particles. The dynamical BH mass corresponds to 1.53 times the DM particle mass. While less than 10% of the star particles lies within the softening length (see legend), their velocity distribution can be quite accurately described by a Maxwellian distribution, as shown in the left plot. Interestingly, the distributions that best fit in the surroundings of the BH (2 ϵ_BH, ϵ_BH) shift their peaks down to lower relative velocities. Particles with lower relative velocity compared to the BH contribute the most to the computation of the DF correction, according to Eq. (6). The same trend can be observed when the host sub-halo is resolved with fewer star particles.

Fig. A.2 PDF fitting of the velocity distribution as Fig. A.1 but for a smaller BH at z = 2.5 in the Cl13 simulation.

Fig. A.3 PDF fitting of the velocity distribution as in Fig. A.1 for a BH in the D9 simulation using the DF correction at z = 2.5. The relative velocity PDF around the BH is locally different from the global distribution.

For the case presented in Fig.A.2, the particle number within the entire galaxy is 1265 and only 49 of them are found within ϵ_BH, thus contributing to our DF correction. Moreover, the BH mass reaches only 6% of the mass of the DM particles. Even if the system is less resolved than that described in Fig. A.1, the dominant contribution to the PDF of velocity in the vicinity of the BH comes from particles having a velocity relative to the BH smaller than 200 km/s, whose contribution is again dominant in the correction of the DF. Finally, when the BH local velocity distribution deviates from the host galaxy velocity distribution profile, this technique can trace the local velocity distribution. Although the DF force is a global effect, whose ‘long-range’ contribution is already captured by the N-body solver, the added correction should have a local nature.

For instance, the velocity distributions in Fig. A.3 describe a peculiar situation. As the BH is displaced from the sub-halo centre of 4 ckpc/h, the velocity distribution in the BH’s surroundings differs from that of the sub-halo. Thus, enclosing only particles in the local vicinity of the BH leads to velocity distributions different from the one of the entire system.

Figures A.2, A.1, and A.3 are only three examples that are representative of the vast majority of the cases we obtained when performing tests on the validity of Eq. A.1. Thus, we conclude from this analysis that for well-resolved systems the velocity PDF from particles within the softening length reproduces fairly well the PDF that characterises the BH underlying velocity field. Although for less-resolved systems the velocity dispersion is less constrained, the contribution from particles with low velocity relative to the BH is dominant. These particles are those that provide the dominant contribution in the DF correction of Eq. 6.

Appendix B Validating the DF correction in an isolated halo

To validate the DF model introduced in this work, we test its performance in a controlled experiment in which a BH particle is placed on an initially radial orbit within an isolated DM halo, which follows a Navarro, Frenk and White (NFW, Navarro et al. 1996) mass density profile. We then compute the sinking timescale of a BH falling inwards to the centre of the halo, when the halo is described at increasing resolution. In this idealised scenario, the mass density and velocity distribution of the halo are fully specified, thus allowing for a direct comparison with analytical results from the original DF formula from Chandrasekhar (1943). Furthermore, the limited computational cost of these simulations enables one to perform multiple simulations at different resolution. This allows us to validate the DF correction of Eq. (2.1), which depends in fact on the numerical resolution.

We generate the initial conditions of a NFW DM-only halo using the MAKEGALAXY code (Springel & White 1999) For the tests shown in this section, we adopted the same halo parameters of Genina et al. (2024). The halo has a virial radius of R_vir = 350.54 kpc, a virial mass of M_vir = 10¹³ M_⊙ and concentration c = 4.38. Having verified that the halo density stabilses into the expected NFW mass density profile and the velocity distribution defines a Maxwellian distribution with dispersion σ(r), we seed a BH at 20 kpc from the centre with a mass of 10⁹ M_⊙, setting it onto a circular orbit around the halo centre. We repeat this seeding process by progressively increasing mass and force resolution by sampling the halo with 10⁵,10⁶,10⁷,5 ⋅ 10⁷ particles within R_vir. We adopt Plummer-equivalent softenings derived from Power et al. (2003). Once seeded, the BH has the same softening length as the surrounding particles. We report in Table B.1 the number of particles, the mass of DM particles, and their softening for each resolution. Every simulation is carried out twice, with and without the DF correction introduced in Sect. 2.1.

Fig. B.1 Sinking timescales for a BH orbiting in a NFW DM halo. We plot the distance to the halo centre found through the shrinking sphere algorithm (Power et al. 2003) as a function of time. We plot the results from numerical simulation without any correction to the BH dynamics (top panel) and with the DF implementation introduced in Sect. 2.1 (bottom panel). The colourmap indicates the number of particles sampling the halo such that darker lines refer to higher numerical resolution. The dashed-dotted vertical lines are indicative sinking timescales and horizontal lines report the softening lengths, both are colour-coded with the simulation they refer to. The dotted lines are the analytical preditions for an infalling BH into a NFW halo, choosing the same BH and halo parameter of the numerical simulations, but varying the maximum impact paramter from bmax = 20 (red) and bmax = 80 (green).

Given this ideal setup, where the density profile is known and the velocity distribution is Maxwellian, we can compare the sinking timescale ts obtained from numerical simulation with the analytical expectation provided by Eq. (8). Following a similar procedure to that presented in Rodriguez et al. (2018), it is possible to compute analytically the path of a BH initially seeded on a circular orbit at ${\text{r}}_{{\text{BH}}_{\text{i}}}$ and with initial velocity ${\text{v}}_{{\text{BH}}_{\text{i}}}$ assuming that its orbit remain circular. A detailed description of these calculations as well as a more extended version of the tests shown here below will be presented in Damiano et al. (in prep).

Figure B.1 shows the distance of the BH from the centre of the halo, determined using the shrinking sphere algorithm (Power et al. 2003), across all simulations performed. The upper and lower panels display the results without and with DF correction, respectively. In both panels, higher resolution corresponds to darker lines, as indicated by the corresponding colour maps. The softening lengths, obtained by multiplying the corresponding Plummer-equivalent values by a factor of 2.8, are marked by horizontal lines, colour-coded according to the corresponding resolution. The vertical dashed-dotted lines indicate approximate sinking timescales which correspond to the time at which the orbits stop further shrinking. The dotted green and red curves correspond to the analytical predictions of the BH decaying orbital radius, which are obtained by setting ${\text{b}}_{\text{min}}\text{\hspace{0.17em}=\hspace{0.17em}}{\text{GM}}_{\text{BH}}/\left({\text{v}}_{\text{BH}}^2+2/3{\sigma }^2\right)$ (Just et al. 2011) for the minimum value of the impact parameters appearing in the Coulomb logarithm and choosing ${\text{b}}_{\max }={\text{r}}_{{\text{BH}}_{\text{i}}}=20$ kpc (red), corresponding to the initial BH orbit radius, and b_max = r_s = 80 kpc (green), corresponding to the scale radius of the NFW halo.

In the test at the lowest resolution, applying the DF correction enables the BH to spiral inward toward the halo centre. Without any correction, the BH remains at a roughly constant distance near the softening length, showing persistent oscillations. As the resolution increases, without the addition of any DF correction (NODF case), the sinking timescale t_s decreases from t_s = 8 Gyr when sampling the halo with 10⁶ particles, to t_s = 6.5 Gyr when using 5 ⋅ 10⁷ particles. On the other hand, the DF correction significantly shortens t_s compared to the NODF case, at fixed resolution. For example, when sampling the halo with 10⁶ particles, t_s ≃ 6.5 Gyr, which is only reached at the highest resolution in the NODF case. In addition, when using the DF correction, the higher-resolution simulations (using 10⁷ and 5 ⋅ 10⁷ particles) converge toward similar sinking timescales. Comparing these results with the analytical predictions, we observe that without the DF correction, the highest-resolution simulation approaches the prediction based on the smaller value of the impact parameter. Using the DF correction, the simulations with a number of particles from 10⁶ to 5 ⋅ 10⁷ show timescales that are fully consistent with the analytical ones. In conclusion, we validate our novel DF correction by setting idealised initial conditions for a halo with an NFW profile, adopting parameters similar to those used by Genina et al. (2024). These controlled experiments enable a reliable comparison with the analytical results derived from theoretical calculations. Our results show that applying the DF correction proposed in Sect. 2.1 drives to sinking timescales which are fully consistent with those predicted from the DF formula in Eq. (8).

References

All Tables

All Figures

Fig. 1 Probability density distribution of the distances between sub-halos identified by SubFind and the closest BH particle within the RHMS of each sub-halo of Cl 13, D9 regions, and Cosmobox. The rows show the results obtained at different redshifts: z = 3 (up), z = 1 (central), and z = 0 (bottom). The column report the results comparing DYNFRIC and REPOS on the left, and DYNFRIC and DYNMASS on the right side. We include a dashed-dot line in each plot indicating the softening length of the BH as a reference for the spatial resolution of the simulation.

In the text

Fig. 2 Cumulative number of wandering BHs to the total number of BHs found in the D9 region as a function of the BH mass at z = 0 (left panel), and as a function of redshift for z = 3, z = 1, and z = 0 (right panel). In the left panel, no mass threshold is adopted to select the BHs, and hence all the BHs with mass above the seeding mass are included (see Table 2). In both panels, the simulation using the DF model, marked with a blue line, is compared with REPOS in red and DYNMASS in green. The wandering BHs are defined as those having a distance larger than two times the RHMS from the closest sub-halo.

In the text

	Fig. 3 Histogram of the number of WBHs as a function of their mass for the D9 simulation at z = 0. We overplot the distribution of DYNFRIC and REPOS on the left and DYNFRIC and DYNMASS on the right.
In the text

Fig. 4 Stellar density projection along the z-axis with depths of 100 kpc (top) and 300 kpc (bottom), centred on the position of a wandering BH at redshiſt z = 3 (tagged as BH 2; left panel), and of a wandering BH at z = 0.06 (tagged as BH 4; right panel) in the D9 simulation using the REPOS scheme. In both panels, the cyan crosses identify the wandering BHs, while the green crosses identify the BHs centred in their host. We show sub-halo centres as white shaded dots, with the white circles indicating the corresponding RHMS. On the top and on the right of each panel we show the gravitational potential Φ of all the star particles (light white dots) and of the BHs (crosses) in the field, projected along the two orthogonal directions.

In the text

Fig. 5 Stellar density maps along the z-axis with a depth of 1 Mpc centred on the most massive halo of D9 at z = 0 in the REPOS (top), DYNMASS (bottom left), and DYNFRIC (bottom right) simulations. The panels are all 1 Mpc on a side. In each panel, we plot with dark-blue crosses the BHs located within the RHMS of the associated sub-halo (the same criterion adopted in Sect. 5.1). BHs lying between RHMS and 2RHMS of the closest sub-halo are indicated with green crosses, while wandering BHs (defined as those located beyond 2RHMS of the closest sub-halo) are shown as light-blue crosses. The values of RHMS of the ten most massive structures correspond to the radii of the circles, each centre on the position of the corresponding sub-halo. The RHMS of the BCG, in yellow, corresponds to the physical size of the yellow circle. The legend in the upper left panel of the plot shows the scaling size of the other sub-structures in the region, marked in white.

In the text

Fig. 6 Cumulative and differential distributions of the number of merger events as a function of redshift. Cl13, D9, and CosmoBox are in the first, second, and third rows, respectively. The REPOS simulation results are marked with a dash-dotted red line (with a shaded area marking the differential distributions), the DYNFRIC simulations with a blue solid line and the DYNMASS simulations with a dashed green line. The curves stop at the redshift corresponding to the occurrence of last merger event.

In the text

	Fig. 7 Ratio between the number of sub-halos hosting more than one BH and the total number of sub-halos hosting at least one BH in D9 simulations. We consider DYNMASS (dashed green line), DYNFRIC (solid blue line) and REPOS (dot-dashed red line) at z = 3, 1, 0.
In the text

Fig. 8 Relationship between BH mass and stellar mass of the host galaxies in the Cl13 (upper panels), D9 (central panels), and CBox simulations (lower panels). Each column contains the results obtained with a different BH dynamics prescription: DYNFRIC (left), DYNMASS (central), REPOS (right). The diamonds refer to the BHs associated with the main halos while circles correspond to sub-halos belonging to halos more massive than 1010 M⊙ Diamonds and circles are colour-coded according to the distance between the sub-halo centres and the associated BHs. For comparison, we also plot observational data from Gaspari et al. (2019) with black dots with errorbars, and the relation obtained by McConnell & Ma (2013) with a green and shaded area.

In the text

Fig. 9 Stellar density maps along the z-axis with depth of 350 kpc (upper row) and 260 kpc (bottom row) of the massive sub-halos having a large separation between sub-halo centre and the associated BH in the CosmoBox simulation, using the DYNFRIC model (blue diamonds in the corresponding plot in Fig. 8, for log10(M*/M⊙) > 10), at redshift z = 0. BHs are marked by green crosses. Upper panels: two close encounters between structures. The halo centre is identified as a red circle. The arrows link to the values of the BH masses (in units as specified in the label). In both the left and right panels, the most massive BH belongs to the off-centre substructure. Bottom panels: a merger event between two substructures at redshift z = 0.19 (left panel) and the resulting merged halo at z = 0 (right panel). The plot on the right side captures an off- centre BH while sinking toward the halo centre.

In the text

Fig. 10 Stellar density map along the z-axis with depth of 10 kpc in the top and central row and 70 kpc in the bottom row, centred on the central BH (BHcen) of Event 1 in the top row, Event 2 in the central row and Event 3 in the bottom row. The maps of stellar density refer to the redshift indicated on the left panel of each row. In each panel, the curve show the orbit of the satellite BH (BHsat) involved in the merger event, colour-coded according to the redshift. The crosses indicate the last position of the satellite BH before either the merging event or the end of the simulations (see Sect. 6). The columns display the BH evolution using different sub-resolution prescriptions, from the left: using DF (first column), using a boosted dynamical mass (second), adopting the repositioning scheme (third) and finally without any correction for BH dynamics (fourth).

In the text

Fig. 11 Left column: evolution of the distance between the central and the satellite BH in Event 1 (top), Event 2 (central), Event 3 (bottom). The dashed green line displays the evolution using the dynamical mass scheme (DYNMASS), the blue solid line is for the DF scheme (DYNFRIC), the dash-dotted red line for the repositioning (REPOS), and the densely dashed orange line is for the case without any sub-resolution prescription (NOCORR). The horizontal grey solid line represents the distance threshold which is necessary for the merger event dmerg to happen. Right column: evolution of the distance between the BHs for the DYNFRIC case (top panel) and ratio between the DF and gravitational forces during each event.

In the text

Fig. A.1 Probability distribution functions of the relative velocity between a massive BH (see A.1) in the Cl13 simulation using DF at redshift z = 2.5, and the surrounding stars within 30 kpc (enclosing all the stars in the sub-halo) (red line), 5 ϵBH (yellow line), 2 ϵBH (green line) and ϵBH (blue line). In the left panel, we compare the best fitting Maxwellian PDFs for these four cases. The right plots superimpose the PDF measured from the simulation with the Maxwellian fit for the 30 kpc case (top) and the ϵBH case (bottom).

In the text

	Fig. A.2 PDF fitting of the velocity distribution as Fig. A.1 but for a smaller BH at z = 2.5 in the Cl13 simulation.
In the text

	Fig. A.3 PDF fitting of the velocity distribution as in Fig. A.1 for a BH in the D9 simulation using the DF correction at z = 2.5. The relative velocity PDF around the BH is locally different from the global distribution.
In the text

Fig. B.1 Sinking timescales for a BH orbiting in a NFW DM halo. We plot the distance to the halo centre found through the shrinking sphere algorithm (Power et al. 2003) as a function of time. We plot the results from numerical simulation without any correction to the BH dynamics (top panel) and with the DF implementation introduced in Sect. 2.1 (bottom panel). The colourmap indicates the number of particles sampling the halo such that darker lines refer to higher numerical resolution. The dashed-dotted vertical lines are indicative sinking timescales and horizontal lines report the softening lengths, both are colour-coded with the simulation they refer to. The dotted lines are the analytical preditions for an infalling BH into a NFW halo, choosing the same BH and halo parameter of the numerical simulations, but varying the maximum impact paramter from bmax = 20 (red) and bmax = 80 (green).

In the text