11institutetext: Dipartimento di Fisica dell’Università di Trieste, Sez. di
Astronomia, via Tiepolo 11, I-34131 Trieste, Italy 22institutetext: INAF –
Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131, Trieste, Italy
33institutetext: IFPU, Institute for Fundamental Physics of the Universe, Via
Beirut 2, 34014 Trieste, Italy 44institutetext: INFN, Instituto Nazionale di
Fisica Nucleare, Via Valerio 2, I-34127, Trieste, Italy 55institutetext: ICSC
- Italian Research Center on High Performance Computing, Big Data and Quantum
Computing, via Magnanelli 2, 40033, Casalecchio di Reno, Italy
66institutetext: INAF-Osservatorio di Astrofisica e Scienza dello Spazio di
Bologna, Via Piero Gobetti 93/3, I-40129 Bologna, Italy 77institutetext:
Dipartimento di Fisica e Astronomia ”Augusto Righi”, Alma Mater Studiorum
Università di Bologna, via Gobetti 93/2, I-40129 Bologna, Italy
88institutetext: Instituto de Astronomía Teórica y Experimental (IATE),
Consejo Nacional de Investigaciones Científicas y Técnicas de la
República Argentina (CONICET), Universidad Nacional de Córdoba, Laprida 854,
X5000BGR, Córdoba, Argentina 99institutetext: Universitäts-Sternwarte München,
Scheinerstr. 1, D-81679, München, Germany 1010institutetext: Max-Plank-
Institut für Astrophysik, Karl-Schwarzschild Strasse 1, D-85740 Garching,
Germany

# Dynamical friction and evolution of black holes in cosmological simulations:
a new implementation in OpenGadget3

Alice Damiano 11224455 Milena Valentini 11224455 Stefano Borgani 1122334455
Luca Tornatore 2255 Giuseppe Murante 223355 Antonio Ragagnin 336677 Cinthia
Ragone-Figueroa 8822 Klaus Dolag 991010

###### Abstract

Aims. We introduce a novel sub-resolution prescription to correct for the
unresolved dynamical friction onto black holes (BHs) in cosmological
simulations, to describe BH dynamics accurately and to overcome spurious
motions induced by numerical effects.

Methods. We implement a sub-resolution prescription for the unresolved
dynamical friction onto BHs in the OpenGadget3 code. We carry out cosmological
simulations of a volume of $\mathrm{16\ cMpc^{3}}$ and zoom-ins of a galaxy
group and of a galaxy cluster. The advantages of our new technique are
assessed in comparison to commonly adopted methods to hamper spurious BH
displacements, i.e. repositioning onto a local minimum of the gravitational
potential and ad-hoc boosting of the BH particle dynamical mass. We inspect
variations in BH demography in terms of offset to the host sub-halo centers,
wandering population of BHs, BH-BH merger rates, and occupation fraction of
sub-halos. We also analyze in detail the impact of the different prescriptions
on individual events of BH interactions.

Results. The newly-introduced dynamical friction correction provides centering
of BHs on host halos which is at least comparable with the other techniques,
becoming gradually more effective as the redshift decreases. It predicts half
as many merger events with respect to the repositioning prescription, with the
advantage of being less prone to leave sub-structures without any central BH.
Simulations featuring our dynamical friction prescription produce a smaller
(by up to $\mathrm{\sim 50\%}$ with respect to repositioning) population of
wandering BHs and final BH masses in good agreement with observations. As for
individual BH-BH interactions, our dynamical friction model captures the
gradual inspiraling of orbits before the merger occurs. By contrast, the
repositioning scheme, in its most classical renditions considered here,
describes extremely fast mergers, while the dynamical mass misrepresents the
black holes’ dynamics, introducing numerical scattering between the orbiting
BHs.

Conclusions. Given its performances in describing the centering of BHs within
host galaxies and the orbiting of BH pair before their merging, our dynamical
friction correction opens interesting applications for an accurate description
of the evolution of BH demography within cosmological simulations of galaxy
formation at different cosmic epochs and within different environments.

##  1 Introduction

Supermassive black holes (SMBHs) reside at the center of massive galaxies and
are considered to affect their evolution profoundly. Numerous studies provided
evidence of relations between the SMBH mass and properties of the host
galaxies (e.g. Kormendy et al., 1993; Magorrian et al., 1998; Ferrarese &
Merritt, 2000; Gebhardt et al., 2000; Merritt & Ferrarese, 2000; Ferrarese et
al., 2001; Haering & Rix, 2004; Gultekin 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. In fact, 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 SMBHs from escaping the centers of the host galaxies, is responsible
to migrate BHs to the galaxy centers, and drives the initial stages of the
merger event between two SMBHs, ultimately leading to the formation of a close
pair (Begelman et al., 1980).

Cosmological simulations represent ideal tools to follow the evolution of
structures given the highly non-linear astrophysical phenomena governing the
interactions of baryonic matter, including gas and stars, in addition to 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 orbits of single collisionless particles has little physical
relevance, since what matters is their ensemble properties. A BH particle, on
the other hand, is introduced in cosmological simulations as an individual
collisionless particle which has a specific physical meaning and whose
presence is capable to significantly impact on 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 BHs in these simulations and the interpretation
given to surrounding particles presents a primary conceptual obstacle when
tracking the trajectories of BH particles. Specifically, scattering
interactions between a BH and surrounding particles can ”heat” the BH orbit. A
spurious, numerical displacement of a BH from the center of the host galaxy is
a major consequence of this ”heating”, which can eventually lead to the
formation of unwanted ”wandering” BHs. Furthermore, it also means a poor
predicting capability of describing 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 the key to predict masses of the Brightest Cluster Galaxies
(BCGs) in good agreement with observations, which instead were overpredicted
by a factor of a few in their previous work (Ragone-Figueroa et al., 2013).

These issues stem from numerical simulations failing to recover the DF force.

Since the magnitude of DF drag 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, our first query in
this paper is 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, e.g. by artificially repositioning the BHs at the position of a local
minimum of the gravitational potential (Springel et al., 2005b; 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, which is contributed by all
the particles within the half-mass radius of the substructure hosting the BH.
Besides, Tremmel et al. (2015) stated that, under the hypotheses of a
sufficient shallow potential surrounding the BH, only unresolved interactions
within the softening length need correction. Still, Chen et al. (2022)
proposed that the application the DF correction should be coupled with a
boosted dynamical mass, to account for those cases in which a BH has a rather
small mass, close to its value at seeding, and it 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 motivates this study.

In this paper, we propose a novel implementation of DF correction, that we
implement in 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 will extensively
discuss in the following, the primary benefit 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 a adequate description of the DF force
acting on BHs? How does it compare against the other numerical ah-hoc
prescriptions (i.e., repositioning, dynamical mass) introduced to mimic DF’s
effect on BH particles? To answer these questions, we will simulate a group-
sized and a cluster-sized halos, with initial conditions from the DIANOGA set
(Bonafede et al., 2011; Bassini et al., 2019, 2020), along with a cosmological
box having a comoving size of 16 cMpc a side. The simulations have been
carried out three times, maintaining identical settings for all parameters,
except for the sub-resolution prescription governing BH dynamics: continuous
repositioning on a local potential minimum, boosted dynamical mass, and our
novel model to correct unresolved dynamical friction. This approach will
enable us not only to focus on systems densely populated and rich in
interactions, such as galaxy clusters and groups but also to undertake a
statistical analysis of the BH population in a cosmological volume.

The paper is organized 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 present the implementation of this
model in the OpenGadget3  code of the description of the SMBH evolution and of
the ensuing AGN feedback. After introducing the details of our test
simulations in Sect. 4, we present in Sect. 5 the results of our analysis of
the general properties of the SMBH population. Sect. 6 will ”zoom” on single
BHs and merger episodes, by comparing in details the small-scale effect of our
sub-resolution model of DF.

##  2 Dynamics of black-hole 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.

In this section, we will present different approaches, including our new one,
that 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 $\Delta\@vec{v}_{\rm M}$ of a
particle of mass $\mathrm{M}$ caused by the interaction with another particle
of mass $m$ and velocity $\@vec{v}$ can be expressed in terms of their impact
parameter $b$ and relative velocity, $\@vec{v}_{0}=\@vec{v}-\@vec{v}_{\rm M}$,
when the two particles were initially at large distance. The components of
$\Delta\@vec{v}_{\rm M}$ parallel and perpendicular to the direction of
$\@vec{v}_{0}$ can be expressed as (Binney & Tremaine, 2008):

$|\Delta\@vec{v}_{\rm
M}|_{\parallel}=\frac{2m|\@vec{v}_{0}|}{({M}+m)}\left[1+\frac{b^{2}|\@vec{v}_{0}|^{4}}{G^{2}({\rm
M}+m)^{2}}\right]^{-1}\,;$ (1) $|\Delta\@vec{v}_{\rm
M}|_{\perp}=\frac{2mb|\@vec{v}_{0}|^{3}}{G({\rm
M}+m)^{2}}\left[1+\frac{b^{2}|\@vec{v}_{0}|^{4}}{G({\rm
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\pi\,b\,db\
{{|\Delta\@vec{v}_{M}|}_{\parallel}}f(\@vec{v})d^{3}\@vec{v}$, where
$f(\@vec{v})$ is the phase-space number density of stars. Integrating this
rate over the impact parameter from $0$ to $b_{\rm max}$ 111In principle this
integration should be performed from a minimum value of the impact parameter
equal to the Schwarzschild radius of the BH. It is unnecessary to consider
this contribution to derive an expression for our DF correction. Therefore, we
will assume the minimum value of $b$ to be zero., we have that the DF force
acting on a BH of mass M from particles having mass $m$ and velocities in the
range $(\@vec{v},\@vec{v}+d^{3}\@vec{v})$ is:

$\frac{d\@vec{v}_{\rm
M}}{dt}\bigg{|}_{\@vec{v}}=2\pi\ln\left[{1+{\Lambda(m,\@vec{v})}^{2}}\right]G^{2}m({\rm
M}+m)f(\@vec{v})d^{3}\@vec{v}\frac{(\@vec{v}-\@vec{v}_{\rm
M})}{|\@vec{v}-\@vec{v}_{\rm M}|^{3}}\,,$ (3)

where:

$\Lambda(m,\@vec{v})=\frac{b_{\rm max}(\@vec{v}-\@vec{v}_{\rm M})^{2}}{G({\rm
M}+m)}.$ (4)

Here, $b_{\rm max}$ is the _maximum impact parameter_. In general, $b_{\rm
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 $\Lambda$ on
mass $m$ and velocity $\@vec{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_{\rm max}$ have been introduced in the literature, all having a degree of
arbitrariness. In the next section, we will provide an extended comparison
between different choices. As in Tremmel et al. (2015), we will assume $b_{\rm
max}$ to be given by the gravitational softening length of the BH particle,
$\epsilon_{\rm 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 $\@vec{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 $\@vec{v}_{{m},i}$:

$f(\@vec{v})=\frac{3}{4\pi{\epsilon_{\rm BH}}^{3}}\sum_{i=1}^{N(<\epsilon_{\rm
BH})}\delta(\@vec{v}-\@vec{v}_{m,i})\,,$ (5)

where the sum is over all the $N(<\epsilon_{\rm BH})$ particles lying at a
distance from the BH smaller than its gravitational softening scale. 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\@vec{v}_{\rm M}}{dt}=\frac{3}{2{\epsilon_{\rm
BH}}^{3}}\sum_{i=1}^{N(\epsilon)}\ln\left[{1+{\Lambda}(m_{i})^{2}}\right]m_{i}({\rm
M}+m_{i})\frac{(\@vec{v}_{m,i}-\@vec{v}_{\rm
M})}{|\@vec{v}_{m,i}-\@vec{v}_{\rm M}|^{3}}\,.$ (6)

The gravitational acceleration of the BH is then corrected by the DF
acceleration $\@vec{a}_{\rm df}$ given by Eq. (6), so that the total
acceleration acting on a single BH particle is:

$\@vec{a}_{\rm BH}=\@vec{a}_{\rm g}+\@vec{a}_{\rm df}\,.$ (7)

where $\@vec{a}_{\rm 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. This approach for the computation of the DF
correction presents some advantages with respect to the models proposed in the
literature. We will 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_{\rm max}$.
For instance, in the work by Hirschmann et al. (2014), $b_{\rm 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-halo222This parametrization required an
on-the-fly execution of the SubFind algorithm (Springel et al., 2001a; Dolag
et al., 2009) to identify the substructures hosting the BHs and to compute
their half-mass radius.. 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, $\epsilon_{\rm BH}$. Accordingly, a
correction term to the DF should be added only on scales smaller than
$\epsilon_{\rm BH}$, so that $b_{\rm max}=\epsilon_{\rm BH}$. In line with
this approach, Pfister et al. (2019) accounted for the DF force from particles
within a sphere centered 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^{-1}$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_{\rm
BH}\gg m$, where $m$ is the mass of the particles around the BH, the dynamical
friction force $\@vec{F}_{\rm DF}$ can be cast in the form

$\@vec{F}_{\rm DF}=-4\pi\rho\left(\frac{G{\rm M_{\rm BH}}}{{v}_{\rm
BH}}\right)^{2}\mathrm{F}(x)\ln(\Lambda)\widehat{v}_{\rm BH}\,.$ (8)

Here $\rho$ is the smoothed density at the position of the BH, contributed by
stellar and DM particles, using the BH smoothing length. Furthermore,
$\widehat{v}_{\rm BH}$ is the versor of the BH velocity relative to the ”sea”
of surrounding particles, while

$\mathrm{F}(x)={\rm erf}(x)-\frac{2x}{\sqrt{\pi}}e^{-x^{2}}\ \ ;\ \
x=\frac{v_{\rm BH}}{\sigma_{v}},$ (9)

with $\sigma_{v}$ the velocity dispersion of the surrounding particles. We
will assume that stars and DM particles exert a dynamical friction 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.
In our approach, without any a-priori assumption on the velocity dispersion of
the surrounding particles, BH motion correction arises from summing over the
single contributions of scattering events taking place below the softening
length. Consequently, each scattering particle will contribute to the
correction with its specific velocity and mass. In this sense, our approach
relaxes the approximation of negligible mass of the particles belonging to the
”sea”, which is needed to derive Eq. (8).

####  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 centers 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 to 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. The choice of this limiting relative
velocity becomes a parameter that should be tuned to the typical merger events
and the resolution of the simulation.

In the simulations presented in this work, we adopted a version of the
repositioning close to the original implementation, 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.

We emphasize that the main purpose of our analysis is to highlight the
potential dangers in the use of an ad-hoc repositioning model.

####  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}\geq 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., 2005a). 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, also this method is prone to spurious effects. For instance, Tremmel
et al. (2015) pointed out that initializing the BH mass to a value hundreds of
times higher than its real value can affect the mass of the host galaxy, thus
unavoidably impacting on its subsequent evolution.

In this Section, we reviewed 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 (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 will 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
preparation; 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 parallelized 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., 2005b).

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. (2005b) 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\times 10^{10}$ M⊙; _(ii)_ the stellar mass is at least 2 per cent of the
total mass and 0.05 times the DM mass of the halo; _(iii)_ the gas mass
reaches a value of 10 per cent of the stellar mass; _(vi)_ 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_{\rm BH,seed}$, is not fixed, but scales with the
amount of stars in the FoF group according to:

$M_{\rm BH,seed}=M_{0}\,\frac{M_{\rm*,h}}{f_{*}\,M_{\rm DM,seed}}\,,$ (10)

where $M_{\rm*,h}$ is the stellar mass assigned to a FoF halo, $f_{*}$ and
$M_{{\rm 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_{0}$, 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 1 ) is larger than $M_{0}$. 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. (2005b).
BH accretion rate is calculated according to Bondi-Hoyle formula (Hoyle &
Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952) as

${\dot{M}_{\rm Bondi}}=\frac{4\pi^{2}M^{2}_{\rm
BH}\alpha\rho}{(c_{s}^{2}+v^{2})^{3/2}}\,,$ (11)

where $\rho$ 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
$\alpha$ 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 ($\mathrm{T>10^{5}}$ K) and the cold gas ($\mathrm{T<10^{5}}$ K), using
$\alpha=10$ and $\alpha=100$ respectively for the hot and the cold gas.

The accretion is always limited to the Eddington accretion rate,

$\dot{M}_{\rm Edd}=\frac{4\pi M_{\rm BH}m_{p}}{\eta_{r}\sigma_{T}c}\,,$ (12)

where $m_{p}$ is the proton mass, $\sigma_{T}$ is the Thompson cross-section.
The parameter $\eta_{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
$\eta_{r}=0.1$. In addition, we allowed the black hole to swallow gas
particles via stochastic accretion, as originally proposed by Springel et al.
(2005b). 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_{\rm merg}$, below
which the two BHs could immediately merge. In our simulations, we adopt a
value of $d_{\rm merg}=5\cdot\epsilon_{\rm 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, we follow Di Matteo et al. (2008) to include in our
simulations the following criterion, that should guarantee that the BHs are
actually gravitationally bound:

$\frac{|\Phi_{\rm BH_{1}}-\Phi_{\rm BH_{2}}|}{a}<0.5\cdot c_{s}^{2}-v_{\rm
rel}^{2}\,,$ (13)

where $\Phi_{\rm BH_{1}},\Phi_{\rm 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, $v_{\rm rel}$ their relative velocity and
$a$ is the scale factor, normalized 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, 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}_{\rm BH}c^{2}\,.$ (14)

In the above expression, $\dot{M}_{\rm BH}=\min(\dot{M}_{\rm
Bondi},\dot{M}_{\rm Edd})$ is the BH mass accretion rate, and $\eta_{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 $\eta_{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 $\eta_{f}=0.05$ during the quasar mode, increasing
to $\eta_{f}=0.2$ during the radio mode.

For the seek of clarity, we summarize in Table 1 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.

Table 1: Summary of the relevant equations and parameters that characterize the sub-resolution model of BH evolution and AGN feedback. See Sect.3 for additional details. Column 1: values of the masses of the different particle species in the high-resolution regions of Dianoga clusters and Cosmobox: $m_{\rm DM}$, $m_{\rm gas}$ and $m_{*,\rm}$ are the masses of DM, gas and star particles, respectively; $\epsilon_{\rm DM},\epsilon_{\rm gas},\epsilon_{*}$, and $\epsilon_{\rm BH}$ are the Plummer-equivalent softening lengths of the different particle species. For BH, gas and star particles, the softenings are fixed in comoving units. For DM particles, we fix the softening in comoving units until $z=2$, and then in physical units until $z=0$. Column 2: criteria for BH seeding: $M_{\rm BH,seed}$ is the seeding mass of a BH; $M_{\rm 0}$ is the minimum BH mass seed; $M_{\rm DM,h}$, $M_{\rm gas,h}$ and $M_{\rm*,h}$ are the halo masses in the DM, gaseous and stellar components; $M_{\rm h}$ is the total halo mass and $M_{\rm h,seed}$ is the minimum halo mass required to seed a BH. Column 3: BH mass accretion rate: $\dot{M}_{\rm Bondi}$ and $\dot{M}_{\rm Edd}$ are defined in Eqs.11, 12. Column 4: criteria for BH-BH merging: $\Delta r_{\rm BH}$ is the distance between the merging BHs, $d_{\rm merg}$ is the minimum threshold distance required for merging, $\Phi_{\rm BH_{1}}$ and $\Phi_{\rm BH_{2}}$ are the values of the gravitational potentials of the two BHs, $v_{\rm rel}$ the module of their relative velocity, $c_{\rm s}$ the local sound speed and $a$ is the scale factor. Column 5: expression of the rate of the feedback energy released by the BH and thermally coupled to the surrounding gas; $c$ is the light speed, $\eta_{\rm r}$ and $\eta_{\rm f}$ are the radiative and coupling efficiency, respectively. Resolution |  Seeding |  Accretion |  Merger |  Feedback
---|---|---|---|---
Mass $m_{\rm DM}=4.69\times 10^{7}$ M⊙ $m_{\rm gas}\ =8.67\times 10^{6}$ M⊙ $\ \ m_{{\rm*}}\ \ =8.7\times 10^{6}$ M⊙ Force softening $\epsilon_{\rm DM}=4.13$ ckpc / $1.38$ kpc $\epsilon_{\rm gas}=1.38$ ckpc $\epsilon_{\rm*}\ \ =0.35$ ckpc $\epsilon_{\rm BH}\ =0.35$ ckpc |  $M_{{\rm BH,seed}}=M_{\rm 0}\frac{M_{\rm*,h}}{f_{*}M_{\rm DM,seed}}$ $M_{\rm 0}=2.7\times 10^{5}$ M⊙ $M_{{\rm h,seed}}=6.94\times 10^{10}$ M⊙ $\frac{M_{\rm*,h}}{M_{\rm{h}}}=2\%$ $\frac{M_{\rm*,h}}{M_{\rm DM,h}}=5\%$ $\frac{M_{\rm gas,h}}{M_{{*,h}}}=10\%$ |  $\min(\dot{M}_{\rm Bondi},\dot{M}_{\rm Edd})$ |  Threshold distance: $\Delta r_{\rm BH}<d_{\rm merg}$ $d_{\rm merg}=5\cdot\epsilon_{\rm BH}$ Gravitationally bound: $\frac{|\Phi_{\rm BH_{1}}-\Phi_{\rm BH_{2}}|}{a}<0.5\cdot c_{\rm s}^{2}-v_{\rm rel}^{2}$ |  $\dot{E}_{\rm feed}=\eta_{\rm r}\eta_{\rm f}\dot{M}_{\rm BH}c^{2}$ $\eta_{\rm r}=0.1$ Quasar : $\eta_{\rm f}=0.05$ Radio: $\eta_{\rm f}=0.2$

##  4 Simulations

Table 2: Set of the simulations. Column 1: simulation name; Column 2: method to account for unresolved DF; Column 3 and 4: virial mass and virial radius, respectively, at $z=0$, of the most massive halo in the simulations. Simulation | Name |  $M_{\rm vir}$ [$10^{13}$ M⊙] |  $\rm R_{\rm vir}$ [kpc]
---|---|---|---
| REPOS | 1.27 | 610.89
Cl13 | DYNMASS | 1.26 | 609.33
| DYNFRIC | 1.26 | 609.75
| REPOS | 18.76 | 1497.17
D9 | DYNMASS | 19.02 | 1504.28
| DYNFRIC | 18.89 | 1499.24
| REPOS | 2.57 | 751.01
CosmoBox | DYNMASS | 2.58 | 751.55
| DYNFRIC | 2.65 | 760.15

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 cMpc), all at the same resolution. This allowed us to
test the new implementation in different environments. All simulations are
performed assuming a $\Lambda$CDM cosmology with $\Omega_{m}=0.24$,
$\Omega_{b}=0.0375$ for the total matter and baryon density parameters,
$h=0.72$ for the Hubble parameter, $n_{s}=0.96$ for the primordial spectral
index, $\sigma_{8}=0.8$ for the power spectrum normalization. We report in
Table 2 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_{\rm
200}\simeq 1.53\times 10^{14}$ M⊙ at $z=0$, while the other refers to a lower-
mass group of M${}_{\rm 200}\simeq 1.38\times 10^{13}$ M⊙, namely Cl13333We
indicate with $M_{\rm 200}$ the total mass contained within a sphere of radius
$\rm R_{\rm 200}$ which contains an overdensity equal to $200\rho_{c}(z)$,
where $\rho_{c}(z)=3H(z)^{2}/8\pi G$ is the critical cosmic density at
redshift $z$. In Table 2 we use the virial radius $\rm R_{\rm vir}$ defined as
the radius within which the mean halo density corresponds to the prediction of
spherical collapse (e.g. Eke et al., 1996). Accordingly, the virial mass $\rm
M_{\rm vir}$ is the mass of a sphere enclosed within a radius $\rm R_{\rm
vir}$.. Both objects have been simulated at the same mass and spatial
resolution (see Table 1).

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 will 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 will present our results using two different
approaches. In Sect. 5, we will conduct a statistical analysis to show how the
method to track BH orbits impacts on the properties of their population. Sect.
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
analyse 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
center 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 Centering 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 $\rm R_{\rm HMS}$). Sub-halos are
identified by the SubFind algorithm (Springel et al., 2001b; 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 $\rm R_{200}$ at least one low-resolution DM particle spuriously
scattered within the high-resolution region. We then calculated the distance
between the (sub-)halo centers and the closest BH. Figure 1 shows the
histogram of the probability density function (PDF) of the number of
(sub-)halos having a BH within their $\rm R_{\rm HMS}$, versus the distance
$\Delta 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. 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-centered 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 center of their host substructures. As for the DYNMASS,
it shows the less pronounced tail for higher distances ($\Delta 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 center of their host halos, without resorting to ad-hoc prescription.
The comparable performance of the two methods is actually a non-trivial and
encouraging result for our DF implementation; in fact, accurate centering 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 3 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 centered sub-halos, with at least
$70\%$ of the total distribution.

Figure 1: Probability density distribution of the distances between sub-halos identified by SubFind and the closest BH particle within the $\rm R_{\rm HMS}$ of each sub-halo of Cl13, 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. Table 3: Percentages of sub-halos from each simulation that contribute to the distributions shown in Fig.1. For each implementation, we indicate the percentage of sub-halos for each simulation: Cl13, D9 and Cosmobox, in the left, central and right columns, respectively. Percentages are listed for different redshifts (as in Fig.1): $z=3$, $z=1$, and $z=0$.  | REPOS | DYNMASS | DYNFRIC
---|---|---|---
| Cl13 | D9 | CosmoBox | Cl13 | D9 | CosmoBox | Cl13 | D9 | CosmoBox
$z=3$ | 11 | 82 | 7 | 12 | 83 | 5 | 12 | 83 | 5
$z=1$ | 10 | 78 | 12 | 9 | 80 | 11 | 10 | 78 | 12
$z=0$ | 13 | 70 | 17 | 10 | 75 | 15 | 10 | 75 | 15

###  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 $\rm R_{\rm HMS}$ of the closest sub-halo.

Figure 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$, in the left panel,
and as a function of redshift for $z=3$, $z=1$, $z=0$ in the right panel. In
the left panel, no mass threshold is adopted to select the BHs, hence all the
BHs with mass above the seeding mass are included (see Table 1). 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 $\rm R_{\rm HMS}$ from
the closest sub-halo.

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 centered, 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 ($\rm M_{\rm BH}>10^{7}\rm M_{\rm\odot}$). This is somewhat
surprising since massive BHs are those for which DF should be better
numerically resolved and, therefore, less prone to wander outside their host
galaxy. Moreover, the repositioning mechanism predicts the highest percentage
of wandering BHs, substantially increasing from $z=1$ to $z=0$. 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. 3.

Figure 3: Stellar density projection along the $z$-axis with depth of $100$
kpc (top) and $300$ kpc (bottom), centered on the position of a wandering BH
at redshift $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 centered in their host. We show sub-halo
centers as white shaded dots, with the white circles indicating the
corresponding $\rm R_{\rm HMS}$. On the top and on the right of each panel we
show the gravitational potential $\Phi$ of all the star particles (light white
dots) and of the BHs (crosses) in the field, projected along the two
orthogonal directions.

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 centered on their host as green
crosses. Furthermore, the centers of neighboring sub-halos are denoted by
shaded white dots, with circles corresponding to the sub-halo $\rm R_{\rm
HMS}$ enclosing them. The side plots display the gravitational potential
$\Phi$ of the stars along the $x$-axis (on the top) and $y$-axis (on the
right). In the top panel of Fig.3, 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 neighbour 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 $\rm R_{\rm HMS}$ of a sub-halo can be a wander: our definition
classifies as not wandering only those BHs located within twice the $\rm
R_{\rm HMS}$ of the closest sub-halo. The bottom panel of Fig. 3 shows the
position at $z=0.06$ of the supermassive wandering BH 4 with a mass of $\rm
10^{10}\ M_{\odot}$. 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, and will
be analyzed in more detail in Sect. 6.

Figure 4: Stellar density maps along the $z$-axis with depth of $1$ Mpc
centered 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 $\rm R_{\rm HMS}$ of the associated sub-halo (the same criterion
adopted in Sect. 5.1). BHs lying between $\rm R_{\rm HMS}$ and $2\ rmR_{\rm
HMS}$ of the closest sub-halo are indicated with green crosses, while
wandering BHs (defined as those located beyond $2\rm R_{\rm HMS}$ 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
centered 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.

Finally, Fig. 4 displays the position of all the BHs in a projected stellar
density map within 1 Mpc from the center 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
colors depending on their distance $\Delta r$ from the closest sub-halo,
comparing it to the $\rm R_{\rm HMS}$ of this sub-halo. Dark-blue crosses
indicate BHs with $\Delta r<\rm R_{\rm HMS}$, green crosses are for $\rm
R_{\rm HMS}\leq\Delta r<2\times\rm R_{\rm HMS}$, and cyan crosses mark the
wandering BHs, i.e. those with $\Delta r\geq 2\times\rm R_{\rm HMS}$. We also
plot the circles the ten most massive sub-halos in the region, each having
radius that is suitably scaled according to the value of its $\rm R_{\rm
HMS}$. Each circle is centered on the sub-halo position and has a radius
proportional to its $\rm R_{\rm HMS}$. The radius of the yellow circle
corresponds to the actual physical value of the RHMS 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 4 shows that the DYNFRIC simulation exhibits the less numerous
population of wandering BHs, 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.

###  5.3 Merger events and multiple BH systems

Figure 5: Cumulative and differential distributions of the number of merger
events as a function of redshift. Cl13, D9 and CosmoBox are on the first,
second and third row, 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.

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_{\rm merg}$ (see Table
1), making them merge on shorter timescales than when the other prescriptions
are adopted (see Sect. 6).

Fig. 5 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 (Cl13
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\simeq 1.5$ in all the simulations.

Figure 6: The ratio between the number of sub-halos hosting more than one BH
to 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$.

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. Fig. 6 shows the ratio
between the number of sub-halos hosting more than one BH within the $\rm
R_{\rm 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 1). In this case, a new BH would be seeded at the center
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 4 the number of seeded BHs in each
simulation.

Table 4: Number of seeded BHs in each simulation for every sub-resolution prescription adopted. Simulation | REPOS | DYNFRIC | DYNMASS
---|---|---|---
Cl13 | 206 | 191 | 184
D9 | 2239 | 2206 | 2169
CosmoBox | 328 | 314 | 323

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.

###  5.4 The $M_{\rm*}$–$M_{\rm 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 7 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_{\rm
thr}=10^{10}M_{\odot}$, we distinguish between the central galaxy, to be
identified with BCG or Brightest Group Galaxy (BGG; marked with diamond
symbols in Fig. 7), and the satellite galaxies, which are hosted within the
substructures identified by SubFind (marked with circles). For the BCGs, we
calculate the stellar mass of the stars belonging to the BCG/BGG. within 70
kpc from the center, to avoid contributions from the intra-cluster light. 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 center of
the structure. For the satellite galaxies, we link the closest BH within their
$\rm R_{\rm HMS}$. Symbols are color-coded according to the distance between
the galaxy centers 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 normalization. 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.

Figure 7: The relationship between BH mass and stellar mass of the host
galaxies for 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 to the main halos while circles
correspond to sub-halos belonging to halos more massive than $10^{10}$ M⊙.
Diamonds and circles are color-coded according to the distance between the
sub-halo centers 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.  Table 5: Stellar mass (left column), associated BH mass (central
column) and its distance from the halo centers (right column) for the most
massive halo identified within each simulation.

$\rm log_{\rm 10}{(M_{\rm*}/M_{\odot})}$ $\rm log_{\rm 10}{(M_{\rm
BH}/M_{\odot})}$ Distance (kpc) Cl13 REPOS 11.28 9.13 0.03 DYNMASS 11.27 9.01
0.10 DYNFRIC 11.25 9.07 0.07 D9 REPOS 11.75 10.0 0.10 DYNMASS 11.89 9.52 0.12
DYNFRIC 11.83 9.47 0.07 CosmoBox REPOS 11.36 9.35 $1.25\times 10^{-3}$ DYNMASS
11.39 9.16 0.16 DYNFRIC 11.47 9.4 0.08

Table 5 presents data on stellar mass, associated BH mass, and the distance of
the BH from the halo center for the most massive halo within each simulation.
The DYNFRIC simulations predict better centered BHs compared to the DYNMASS
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 center of
their host galaxies. The results for the D9 region using REPOS are in good
agreement with the observational data for $M_{\rm*}<10^{11}$ 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. 3 that reached the center 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 centering 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. 8, which shows the projection of
the stellar density centered 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. 7 follows from the association of the most massive BH to the central
halo, marked in Fig.8 with a red circle. This feature, rather than being a
real off-centering 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.8, hosts a truly
off-center 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
center of the merged sub-halos, hence the displacement at $z=0$ between the
center identified by SubFind and the BH.

Figure 8:  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 center and the associated BH in the CosmoBox
simulation, using the DYNFRIC model (blue diamonds in the corresponding plot
in Fig.7, for $\text{log}_{10}(M_{\rm*}/M_{\odot})>10$), at redshift $z=0$.
BHs are marked by green crosses. Upper panels: two close encounters between
structures. The halo center 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-center
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-center BH while sinking
toward the halo center.

##  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 Cl13-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
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. In
particular, our analysis focuses on the evolution of binary BH systems,
typically characterized by a massive BH located at the centre of a
substructure, indicated as $\mathrm{BH_{\rm cen}}$, and a second displaced BH,
the ”satellite”, labelled as $\mathrm{BH_{\rm 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_{\rm m}=\mathrm{m_{\rm
BH,{cen}}}/\mathrm{m_{\rm 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 $\mathrm{4}$ kpc from each other. The initial BH mass
ratio is $f_{\rm m}=1.1$, and both the BHs have a mass smaller than $m_{\rm
DM}$ (see Table 6). 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_{\rm 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 at $\mathrm{z=1.26}$. The BHs have a large mass ratio
$f_{\rm m}=313$. This time, only $\rm BH_{\rm sat}$, whose mass is less than
$m_{\rm DM}$, undergoes the dynamical mass correction in the DYNMASS
simulation.

Table 6: Column 1 and 2: initial true masses of the central ($m_{BH,{\rm cen}}$) and satellite ($m_{BH,{\rm sat}}$) BHs in the events described in Sect.6, at the initial redshifts $z=3$ (for Events 1 and 2) and $z=1.26$ (for Event 3). In the Events 1 and 2, the BH masses of both $\rm BH_{\rm cen}$ and $\rm BH_{\rm sat}$ are boosted in the DYNMASS runs. On the other hand, we only boost the satellite BH mass in the Event 3. Column 3: mass of the sub-halo hosting the events. Column 4: initial distance between the two BHs. |  $m_{\rm BH,{cen}}$ (M⊙) |  $m_{\rm BH,{sat}}$ (M⊙) |  $M_{sub}$ (M⊙) |  $\Delta r$ (kpc)
---|---|---|---|---
Event 1 |  $3.23\times 10^{5}$ |  $3.63\times 10^{5}$ |  $6.96\times 10^{10}$ |  4.08
Event 2 |  $3.60\times 10^{6}$ |  $7.20\times 10^{5}$ |  $8.11\times 10^{10}$ |  1.48
Event 3 |  $6.39\times 10^{8}$ |  $2.04\times 10^{6}$ |  $4.50\times 10^{11}$ |  10
Figure 9: 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, centered on the central BH
($\rm BH_{\rm cen}$) of Event 1 in the top row, Event 2 in the central row and
Event 3 in the bottom row. The maps refer to the final redshift at which each
event is followed. In each panel, the line is the path of the satellite BH
($\rm BH_{\rm sat}$) involved in the merger event, color-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).

The characteristics of these three events are summarized in Table 6. Figure 4
provides a graphical representation of the stellar density map at the redshift
marking the conclusion of the simulation or the BH merger, whichever occurs
first. The figure also includes the trajectory of the $\rm BH_{\rm sat}$,
color-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 10 displays on the left column, for the different events in each row,
the distance between $\mathrm{BH_{\rm cen}}$ and $\mathrm{BH_{\rm sat}}$,
namely $\Delta 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_{\rm merg}$ to happen (see Table 1). The right side of Fig.10 focuses
on the results of the DYNFRIC runs. In particular, we plot $\Delta 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 $\mathrm{BH_{\rm
cen}}$ (light-blue) and $\mathrm{BH_{\rm sat}}$ (light green) in the bottom
panel. In the following subsections, we will study each event separately.

###  6.1 Event 1

Figure 10: 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 $d_{\rm merg}$ 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.

Looking at the upper row of Fig. 10, 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\sim 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_{\rm 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 1. Whenever they are closer
than $d_{\rm merg}$, they merge if the gravitational binding criterion is also
satisfied. From the visual representation of Fig.9 we can infer that DYNFRIC,
DYNMASS and NOCORR show similar paths for $\rm BH_{\rm sat}$ until it first
crosses the denser region of the sub-halo. Then, DYNFRIC and NOCORR
simulations bound $\rm BH_{\rm sat}$ in the core of the host. Besides, DYNMASS
reproduces constant oscillation at rather fixed apocentric distances. On the
other hand, the path of $\rm BH_{\rm sat}$ in the REPOS simulation is
discontinuous and short compared to the other: through successive
repositioning steps, the BH drops at the center 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. 9 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. 9 and of Fig.10 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 $\rm BH_{\rm 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.9 demonstrate that DYNFRIC, DYNMASS and NOCORR simulations all produce an
initial tilt of the orbit of $\rm BH_{\rm 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 center. The right panel in Fig.10
illustrates the $F_{\rm DF}/F_{\rm grav}$ trend, showing that for this Event 2
the peaks of the contribution of the DF on $\rm BH_{\rm sat}$ usually take
place in correspondence of the minimum distance between the two BHs: as $\rm
BH_{\rm sat}$ approaches $\rm BH_{\rm 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 $\rm BH_{\rm cen}$
is relatively stronger during all the duration of the event, consistent with
the the more central position that it occupies.

###  6.3 Event 3

Finally, the bottom rows of Fig. 9 and of Fig. 10 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 $\mathrm{BH_{\rm cen}}$ is
more than 300 times higher than that of $\mathrm{BH_{\rm 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 $\mathrm{10}$ kpc, to merge nearly instantaneously,
with a ”long-range teleporting” of $\rm BH_{\rm sat}$ from the outskirt to the
center of the substructure, reported in the top-right inset in the bottom left
panel of Fig. 9. 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
gravitational binding criterion remains unsatisfied for the DYNMASS
simulation. $BH_{\rm sat}$, having its dynamical mass boosted, keeps revolving
around $\mathrm{BH_{\rm 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 $\rm BH_{\rm sat}$ further demonstrate
that DYNFRIC, DYNMASS and NOCORR produce almost identical results until
$z\simeq 1.1$. In the DYNMASS simulation, the subsequent oscillations of the
BH seem to be results of ”kicks” that $\rm BH_{\rm sat}$ receives while
crossing the denser core of the sub-halo. In the right panel, we note that
$F_{\rm df}/F_{\rm grav}$ for $\rm BH_{\rm sat}$ oscillates in phase with the
distance between the BHs. We verified that this feature is driven by wide
oscillations of the gravitational force, that increases when approaching the
galactic core where $\mathrm{BH_{\rm cen}}$ resides. Furthermore, the value of
$F_{\rm df}/F_{\rm grav}$ for $\rm BH_{\rm cen}$ is significantly higher than
for $\rm BH_{\rm 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 by Sect. 6.2. 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.

##  7 Conclusions

We introduced a novel approach to partially correct for the effect of
unresolved dynamical friction (DF) on the orbits of black hole (BH) particles
in cosmological hydrodynamical simulations. We implemented it in the
OpenGadget3 code. The main motivation for this study is to introduce a more
physically motivated method to prevent BH particles in simulations to
spuriously leave the host galaxies, as a result of finite numerical
resolution, and to improve the description of BH-BH mergers. This DF
correction has been extensively tested both in zoom-in and fully-cosmological
simulations. Specifically, we ran two zoom-in simulations, one for a group-
sized and another for a cluster-sized halo, and a cosmological volume having a
box size of 16 $\mathrm{cMpc}^{3}$.

We assessed the new model’s performance for unresolved DF 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 sort out the effect of variations in the prescription for
BHs dynamics, all simulations were run with identical resolution and sub-
resolution physics.

We summarize here below the main conclusions of our analysis.

  * •

Offset: The simulations using the continuous repositioning of BH particles on
the local potential minimum (REPOS scheme) exhibit the smallest offset between
BHs and host sub-halo centers. Our model to correct unresolved DF (DYNFRIC
scheme) provides an accurate centering, even outperforming at $\mathrm{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.3 and Fig.4.

  * •

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

  * •

$\mathbf{M_{\rm\textbf{*}}-M_{\rm\textbf{BH}}}$ relation: The good agreement
with the observational data in reproducing the observed $M_{\rm*}-M_{\rm BH}$
scaling relation in all the simulations demonstrates a relative insensitivity
of this diagnostic to the particular prescription adopted, while it is
sensitive to the choice of the parameters regulating BH accretions and the
ensuing AGN feedback efficiency. See Sect. 5.4 and Fig. 7.

To delve deeper in the details of how different prescriptions affect the
orbits of BH particles, we focused on specific events of BH-BH interactions.
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. 9 and Fig. 10, are summarized
here below.

  * •

DF: Our novel model to correct DF predicts dumped orbits of the satellite BHs,
which gradually approach the center of a host galaxy, thus gently driving to
the formation of 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 merging events. In fact, such interactions eventually inhibit the BHs
from satisfying the gravitational boundedness criteria for merging. The
overall results in terms of correct description of orbit 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 for extremely
rapid mergers, preceded by large and sudden movements of the satellite BH,
that promptly reaches the central region of the sub-halo with few ”jumps”,
that can eventually span several kpcs.

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 reaches at least the same performances of other ad-
hoc numerical prescriptions in terms of centering BH in their host halos,
while sensibly reducing the population of wandering BHs, and overcoming the
limitations of such prescription in the description of BH-BH merger events.
Thanks to the recent work by Sala et al. (2023), the OpenGadget3 code is now
equipped with a spin-evolution model that coupled to the DF correction can
rely on a more precise description of BH dynamics. Also, the introduction of a
DF model is crucial for further investigations on the BH merger rates from
cosmological hydrodynamical simulations (e.g. DeGraf et al., 2024), serving as
a powerful tool to fully exploit the new window offered by gravitational wave
astrophysics.

###### Acknowledgements.

We thank Gian Luigi Granato for a number of useful discussions and comments on
this paper, that contributed to improve the presentation of the results.
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 program 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.

## References

  * Bahé et al. (2022) Bahé , Y. M., Schaye, J., Schaller, M., et al. 2022, Monthly Notices of the Royal Astronomical Society, 516, 167
  * Bassini et al. (2020) Bassini, L., Rasia, E., Borgani, S., et al. 2020, A&A, 642, A37
  * Bassini et al. (2019) Bassini, L., Rasia, E., Borgani, S., et al. 2019, A&A, 630, A144
  * Beck et al. (2016) Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
  * Bellovary et al. (2018) Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2018, Monthly Notices of the Royal Astronomical Society
  * Bertocco et al. (2020) Bertocco, S., Goz, D., Tornatore, L., et al. 2020, in Astronomical Society of the Pacific Conference Series, Vol. 527, Astronomical Data Analysis Software and Systems XXIX, ed. R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter, 303
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bird et al. (2022) Bird, S., Ni, Y., Matteo, T. D., et al. 2022, Monthly Notices of the Royal Astronomical Society, 512, 3703
  * Bonafede et al. (2011) Bonafede, A., Dolag, K., Stasyszyn, F., Murante, G., & Borgani, S. 2011, Monthly Notices of the Royal Astronomical Society, 418, 2234
  * Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195
  * Bondi & Hoyle (1944) Bondi, H. & Hoyle, F. 1944, MNRAS, 104, 273
  * Booth & Schaye (2009) Booth, C. & Schaye, J. 2009, Monthly Notices of the Royal Astronomical Society, 398, 53
  * Borgani & Kravtsov (2011) Borgani, S. & Kravtsov, A. 2011, Advanced Science Letters, 4, 204
  * Bower et al. (2006) Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645
  * Bromm & Loeb (2003) Bromm, V. & Loeb, A. 2003, ApJ, 596, 34
  * Cattaneo et al. (2009) Cattaneo, A., Faber, S. M., Binney, J., et al. 2009, Nature, 460, 213, arXiv:0907.1608 [astro-ph]
  * Chandrasekhar (1943) Chandrasekhar, S. 1943, ApJ, 97, 255
  * Chen et al. (2022) Chen, N., Ni, Y., Tremmel, M., et al. 2022, MNRAS, 510, 531
  * Combes et al. (2019) Combes, F., García-Burillo, S., Audibert, A., et al. 2019, A&A, 623, A79
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, Monthly Notices of the Royal Astronomical Society, 486, 2827
  * DeBuhr et al. (2011) DeBuhr, J., Quataert, E., & Ma, C.-P. 2011, Monthly Notices of the Royal Astronomical Society, 420, 2221
  * DeGraf et al. (2024) DeGraf, C., Chen, N., Ni, Y., et al. 2024, MNRAS, 527, 11766
  * Di Matteo (2023) Di Matteo, Daniel Angles-Alcazar, F. S. 2023 [arXiv:2304.11541]
  * Di Matteo et al. (2008) Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33
  * Dolag et al. (2009) Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497
  * Dubois et al. (2014) Dubois, Y., Volonteri, M., & Silk, J. 2014, Monthly Notices of the Royal Astronomical Society, 440, 1590–1606
  * Eke et al. (1996) Eke, V. R., Cole, S., Frenk, C. S., & Navarro, J. F. 1996, MNRAS, 281, 703
  * Fabjan et al. (2010) Fabjan, D., Borgani, S., Tornatore, L., et al. 2010, MNRAS, 401, 1670
  * Fabjan et al. (2010) Fabjan, D., Borgani, S., Tornatore, L., et al. 2010, Monthly Notices of the Royal Astronomical Society, 401, 1670
  * Ferrarese & Merritt (2000) Ferrarese, L. & Merritt, D. 2000, The Astrophysical Journal, 539, L9, arXiv:astro-ph/0006053
  * Ferrarese et al. (2001) Ferrarese, L., Pogge, R. W., Peterson, B. M., et al. 2001, Supermassive Black Holes in Active Galactic Nuclei. I. The Consistency of Black Hole Masses in Quiescent and Active Galaxies, arXiv:astro-ph/0104380
  * Gaspari et al. (2019) Gaspari, M., Eckert, D., Ettori, S., et al. 2019, The Astrophysical Journal, 884, 169, arXiv:1904.10972 [astro-ph, physics:physics]
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, The Astrophysical Journal, 539, L13, arXiv:astro-ph/0006289
  * Gitti et al. (2012) Gitti, M., Brighenti, F., & McNamara, B. R. 2012, Advances in Astronomy, 2012, 1
  * Graham & Driver (2007) Graham, A. W. & Driver, S. P. 2007, The Astrophysical Journal, 655, 77, arXiv:astro-ph/0607378
  * Granato et al. (2004) Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580
  * Groth et al. (2023) Groth, F., Steinwandel, U. P., Valentini, M., & Dolag, K. 2023, arXiv e-prints, arXiv:2301.03612
  * Gultekin et al. (2009) Gultekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, The M-sigma and M-L Relations in Galactic Bulges and Determinations of their Intrinsic Scatter, arXiv:0903.4897 [astro-ph]
  * Haering & Rix (2004) Haering, N. & Rix, H.-W. 2004, The Astrophysical Journal, 604, L89, arXiv:astro-ph/0402376
  * Hirschmann et al. (2014) Hirschmann, M., Dolag, K., Saro, A., et al. 2014, Monthly Notices of the Royal Astronomical Society, 442, 2304
  * Hopkins et al. (2005) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2005, A Unified, Merger-Driven Model for the Origin of Starbursts, Quasars, the Cosmic X-Ray Background, Supermassive Black Holes and Galaxy Spheroids, arXiv:astro-ph/0506398
  * Hoyle & Lyttleton (1939) Hoyle, F. & Lyttleton, R. A. 1939, Proceedings of the Cambridge Philosophical Society, 35, 405
  * Kormendy et al. (1993) Kormendy, J., Beckman, J., Colina, L., & Netzer, H. 1993, Consejo Superior de Investigaciones Cientificas (Madrid, 1993)
  * Ma et al. (2023) Ma, L., Hopkins, P. F., Kelley, L. Z., & Faucher-Giguère, C.-A. 2023, Monthly Notices of the Royal Astronomical Society, 519, 5543
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, The Astronomical Journal, 115, 2285, arXiv:astro-ph/9708072
  * Mayer & Bonoli (2018) Mayer, L. & Bonoli, S. 2018, Reports on Progress in Physics, 82, 016901
  * McConnell & Ma (2013) McConnell, N. J. & Ma, C.-P. 2013, The Astrophysical Journal, 764, 184, arXiv:1211.2816 [astro-ph]
  * Menezes et al. (2014) Menezes, R. B., Steiner, J. E., & Ricci, T. V. 2014, The Astrophysical Journal, 796, L13
  * Merritt & Ferrarese (2000) Merritt, D. & Ferrarese, L. 2000, The M(BH)-Sigma Relation for Supermassive Black Holes, arXiv:astro-ph/0008310
  * Noble et al. (2011) Noble, S. C., Krolik, J. H., Schnittman, J. D., & Hawley, J. F. 2011, The Astrophysical Journal, 743, 115
  * Novikov & Thorne (1973) Novikov, I. D. & Thorne, K. S. 1973, Black holes (Les astres occlus), 1, 343
  * Ostriker (1999) Ostriker, E. C. 1999, The Astrophysical Journal, 513, 252
  * Pfister et al. (2019) Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, Monthly Notices of the Royal Astronomical Society, 486, 101
  * Pillepich et al. (2017) Pillepich, A., Springel, V., Nelson, D., et al. 2017, Monthly Notices of the Royal Astronomical Society, 473, 4077
  * Ragagnin et al. (2020) Ragagnin, A., Dolag, K., Wagner, M., et al. 2020, arXiv e-prints, arXiv:2003.10850
  * Ragagnin et al. (2016) Ragagnin, A., Tchipev, N., Bader, M., Dolag, K., & Hammer, N. J. 2016, in Advances in Parallel Computing, 411–420
  * Ragone-Figueroa et al. (2018) Ragone-Figueroa, C., Granato, G. L., Ferraro, M. E., et al. 2018, Monthly Notices of the Royal Astronomical Society
  * Ragone-Figueroa et al. (2013) Ragone-Figueroa, C., Granato, G. L., Murante, G., Borgani, S., & Cui, W. 2013, MNRAS, 436, 1750
  * Reines et al. (2020) Reines, A. E., Condon, J. J., Darling, J., & Greene, J. E. 2020, The Astrophysical Journal, 888, 36
  * Sala et al. (2023) Sala, L., Valentini, M., Biffi, V., & Dolag, K. 2023, Supermassive black hole spin evolution in cosmological simulations with OpenGadget3
  * Schaye et al. (2014) Schaye, J., Crain, R. A., Bower, R. G., et al. 2014, Monthly Notices of the Royal Astronomical Society, 446, 521
  * Sijacki & Springel (2006) Sijacki, D. & Springel, V. 2006, Monthly Notices of the Royal Astronomical Society, 366, 397
  * Sijacki et al. (2015) Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, Monthly Notices of the Royal Astronomical Society, 452, 575
  * Silk & Rees (1998) Silk, J. & Rees, M. J. 1998, Quasars and Galaxy Formation, arXiv:astro-ph/9801013
  * Springel (2005) Springel, V. 2005, Monthly Notices of the Royal Astronomical Society, 364, 1105
  * Springel (2016) Springel, V. 2016, in Saas-Fee Advanced Course, Vol. 43, Saas-Fee Advanced Course, ed. Y. Revaz, P. Jablonka, R. Teyssier, & L. Mayer, 251
  * Springel et al. (2005a) Springel, V., Di Matteo, T., & Hernquist, L. 2005a, Monthly Notices of the Royal Astronomical Society, 361, 776
  * Springel & Hernquist (2003) Springel, V. & Hernquist, L. 2003, MNRAS, 339, 289
  * Springel et al. (2005b) Springel, V., Matteo, T. D., & Hernquist, L. 2005b, Monthly Notices of the Royal Astronomical Society, 361, 776
  * Springel et al. (2001a) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001a, MNRAS, 328, 726
  * Springel et al. (2001b) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001b, MNRAS, 328, 726
  * Steinborn et al. (2015) Steinborn, L. K., Dolag, K., Hirschmann, M., Prieto, M. A., & Remus, R.-S. 2015, MNRAS, 448, 1504
  * Taffoni et al. (2020) Taffoni, G., Becciani, U., Garilli, B., et al. 2020, in Astronomical Society of the Pacific Conference Series, Vol. 527, Astronomical Data Analysis Software and Systems XXIX, ed. R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter, 307
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Tornatore et al. (2007) Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050
  * Tremmel et al. (2015) Tremmel, M., Governato, F., Volonteri, M., & Quinn, T. R. 2015, Monthly Notices of the Royal Astronomical Society, 451, 1868
  * Vogelsberger et al. (2013) Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, Monthly Notices of the Royal Astronomical Society, 436, 3031
  * Volonteri & Bellovary (2012) Volonteri, M. & Bellovary, J. 2012, Reports on Progress in Physics, 75, 124901
  * Webb et al. (2012) Webb, N., Cseh, D., Lenc, E., et al. 2012, Science, 337, 554
  * Wiersma et al. (2009) Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99
  * Wurster & Thacker (2013) Wurster, J. & Thacker, R. J. 2013, Monthly Notices of the Royal Astronomical Society, 431, 2513