# Statistical analysis of kicked black holes from TNG300 simulation

M. Smole  Astronomical Observatory, Volgina 7, 11060 Belgrade, Serbia  [ M.
Micic  Astronomical Observatory, Volgina 7, 11060 Belgrade, Serbia

###### Abstract

Asymmetric emission of gravitational waves during mergers of black holes (BHs)
produces a recoil kick, which can set a newly formed BH on a bound orbit
around the center of its host galaxy, or even completely eject it. To study
this population of recoiling BHs we extract properties of galaxies with
merging BHs from Illustris TNG300 simulation and then employ both analytical
and numerical techniques to model unresolved process of BH recoil.

This comparative analysis between analytical and numerical models shows that,
on cosmological scales, numerically modeled recoiling BHs have a higher escape
probability and predict a greater number of offset active galactic nuclei
(AGN). BH escaped probability $>40~{}\%$ is expected in 25 $\%$ of merger
remnants in numerical models, compared to 8$\%$ in analytical models. At the
same time, the predicted number of offset AGN at separations $>5$ kpc changes
from 58 $\%$ for numerical models to 3 $\%$ for analytical models. Since BH
ejections in major merger remnants occur in non-virialized systems, static
analytical models cannot provide an accurate description. Thus we argue that
numerical models should be used to estimate the expected number density of
escaped BHs and offset AGN.

###### keywords:

black hole physics – gravitational waves – galaxies: interactions – methods:
numerical

M<EMAIL_ADDRESS>
##  1 Introduction

The concept of hierarchical growth of structures suggests that galaxy mergers
play an important role in their evolution. Alongside building their masses,
galaxy mergers can reshape their structure, enhance star formation and
potentially trigger active galactic nuclei (AGN). Both theory and observations
reveal the potential of major galaxy mergers to fuel enough gas to the central
galaxy regions and increase the central black hole (BH) accretion rate (e.g.
cox (13); lambas (32); patton (44); Weston (63); Ellison (15); montero (50);
Byrne (6)).

When two galaxies merge, the central black holes (BHs) gradually lose energy
and move toward the center of the merged galaxy due to dynamical friction.
This results in the formation of a binary BH, which can continue to harden
through interactions with stars and gas. The formation of a close BH binary
system is significantly affected by the morphology and mass ratio of
progenitor galaxies. Major mergers of gas-rich galaxies are expected to result
in effective hard BH binary formation (Escala1 (17, 18); Mayer (39); roskar
(51); delValle (14); Goicovic (20); Khan_2016 (28)), while some minor mergers
do not necessarily lead to central BHs merger (Callegari_2009 (8, 7); Tremmel
(60)).

In the final stages of BH merger, the emission of gravitational waves will
effectively remove angular momentum and energy from the binary system, leading
to a rapid collision Begelman (2). Those gravitational waves could be directly
detected with the upcoming space-based interferometer LISA (Laser
Interferometer Space Antenna, Amaro (1)), whose focus will be on the mergers
of BHs in the mass range $10^{4}-10^{7}\>{\rm M_{\odot}}$. Any asymmetry in
the binary system, such as the presence of BHs with unequal masses or spins,
can cause the asymmetric emission of gravitational radiation and BH recoil.
The emission of gravitational waves has a preferred direction due to a non-
zero net linear momentum, and this causes the center of mass of the binary to
move in the opposite direction Redmount (48). During this process, the newly
formed BH receives a recoil whose amplitude depends on the mass ratio of the
merging BH, the magnitude and orientation of their spins with respect to the
binary’s orbital plane, and the eccentricity of the orbit. Mergers of fast
rotating BHs with spin vectors in the orbital plane and orientated in opposite
directions result with super-kicks as high as $~{}4000\>{\rm~{}km}\,{\rm
s}^{-1}$, while efficient spin alignment during BH merger reduces kick
velocities (gonz (21); camp (9); lousto (34)).

Depending on the amplitude of the kick velocity, a recoiling BH could be
either placed on a bound orbit around the galaxy center or be completely
ejected from its host. Bound BHs can spend up to several Gyr oscillating with
low amplitude around the galaxy center (e.g., Gualandris (22); Komossa (29)),
while this time can be significantly reduced in gas-rich systems (e.g.,
Blecha08 (4); Guedes (23); Sijacki (52)). In the case that merged BH receives
a recoil velocity higher than the escape velocity from the host galaxy, it
will be completely ejected. Those wandering BHs are less likely to merge with
other BHs, thus gravitational wave recoil could have a negative influence on
merger-driven BH growth (e.g., Haiman (24); Merritt (40); Volonteri (61)). If
active at the time of the merger, the accretion disk will be ejected along
with a recoiling BH, during which time it can be observed as spatially or
kinematically offset AGN (e.g. Madau (36); Loeb (33); Blecha08 (4); Blecha11
(3)). There is a large number of observed offset AGN candidates in the
literature (e.g. Komossa1 (30); Robinson (49); Jonker (25); Civano10 (12, 11);
Eracleous (16); Koss (31); Markakis (38); Chiaberge (10) ;Kalfountzou (26)).
More recently, Stemo21 (59) composed a catalog of 204 offset and dual AGN,
using the Advanced Camera for Surveys Active Galactic Nuclei Catalog Stemo20
(58) of AGNs observed by Hubble Space Telescope. However, none of these cases
have been confirmed since alternative explanations, such as a binary BH system
with only one active BH, cannot be ruled out.

Even though cosmological simulations provide a powerful tool to investigate
the population of merging galaxies, the actual process of a BH merger cannot
be simulated due to the resolution limits. Instead, the results of
cosmological simulations are often combined with analytical models for
unresolved processes. Blecha (5) investigated trajectories of recoiling BHs
using Illustris cosmological simulations (Vogelsberger et al. 2014a,b ; illu3
(19); illu4 (53)). By extracting properties of merging BHs and progenitor
galaxies directly from the simulation, authors constructed analytical models
of merger remnant galaxies to estimate escape velocities, BH kick velocities
and offset AGN distribution.

On the other hand, higher mass resolution achieved by simulations of isolated
galaxies can provide more detailed insights into processes within the merger
remnant galaxy. smole19 (54) showed that numerical models provide a more
realistic representation of merging galaxies compared to analytical models.
Mass redistribution during interaction can lead to a 25 $\%$ lower escape
velocities in numerical models of major merger remnant galaxies (for more
details we refer to Methods).

In this work, we combine the results of a cosmological simulation on one side
and isolated galaxy simulations on the other side. We apply the similar
analysis described by Blecha (5) to state-of-the-art The Next Generation
Illustris Simulations111https://www.tng-project.org/ (Illustris TNG) project
(Springel18 (57); Naiman (42); Pillepich (46); Nelson (43); Marinacci (37)).
We employ the simulation with the largest volume, TNG300, to extract
properties of galaxies with merging BHs. In addition to analytical models, we
used the results of numerical simulation of isolated galaxy mergers in order
to estimate escape velocities smole19 (54). We investigate the difference
between the predicted number of recoiling BHs in analytical and numerical
models on cosmological scales and estimate the probability for complete BH
ejections during major galaxy mergers.

In Methods we introduce the employed techniques, together with data sets used
in this work, and we describe their implementation in our methodology. Our
findings are presented and discussed in Results. In Conclusions we summarise
the main findings of this work and draw conclusions.

##  2 Methods

###  2.1 Illustris TNG300

In order to calculate the statistics of recoiling BHs in numerical and
analytical models, we use publicly available data from Illustris TNG project
(Springel18 (57); Naiman (42); Pillepich (46); Nelson (43); Marinacci (37)).
Illustris TNG cosmological hydrodynamical simulations of galaxy formation were
performed using the arepo code (Springel10 (56)). The simulations were
initialized at redshift $z=127$ using planck (47) cosmological parameters:
matter density $\Omega_{\rm{m}}=0.3089$, baryon density
$\Omega_{\rm{b}}=0.0486$, dark energy density $\Omega_{\Lambda}=0.6911$,
Hubble constant $H_{\rm{0}}=0.6774$ $\rm{km}/\rm{s}/\rm{Mpc}$, power spectrum
normalization $\sigma_{\rm{8}}=0.8159$ and a primordial spectral index
$n_{\rm{s}}=0.9667$. For the analysis employed in this work we choose the
simulation with the largest cosmological box $(\sim 300~{}\rm{Mpc})^{3}$,
TNG300. Mass resolution for TNG300 simulation is $5.9\times 10^{7}\>{\rm
M_{\odot}}$ for DM particles and $1.1\times 10^{7}\>{\rm M_{\odot}}$ for
baryonic component.

Results of Illustris TNG300 simulation are stored in 100 snapshots, from
redshift $z=20$ to $z=0$. Snapshot files contain information about each
particle in the simulation volume. In TNG300 simulation a BH particle with a
seed mass of $1.18\times 10^{6}\>{\rm M_{\odot}}$ is placed in each galaxy
more massive than $7.8\times 10^{10}\>{\rm M_{\odot}}$, which do not already
have a BH, and then allowed to grow via gas accretion and by mergers. BH
accretion is modeled according to Eddington-limited Bondi accretion rate:

$\dot{M}_{\textrm{Bondi}}=\frac{4\pi G^{2}M_{\text{BH}}^{2}\rho}{c_{s}^{3}},$
(1a)

$\dot{M}_{\textrm{Edd}}=\frac{4\pi
GM_{\textrm{BH}}m_{p}}{\epsilon_{r}\sigma_{\rm T}}c,$ (1b)

$\dot{M}=\min\left(\dot{M}_{\textrm{Bondi}},\dot{M}_{\text{Edd}}\right),$ (1c)

where $M_{\text{BH}}$ is BH mass, $G$ is the gravitational constant, $\rho$
the kernel-weighted ambient density around the BH, $c_{s}$ the kernel-weighted
ambient sound speed including the magnetic signal propagation speed, $m_{p}$
the proton mass, $c$ the speed of light, $\epsilon_{r}=0.2$ the BH radiative
efficiency and $\sigma_{\rm T}$ the Thompson cross-section. We refer to
Weinberger (62) for more details about BH growth and the feedback processes
included.

An additional supplementary data catalog is available providing details of
each BH-BH merger in the simulation, otherwise not available from the
snapshots alone. This file contains information about both merging BH ID
numbers (unique throughout the simulation), their masses immediately preceding
the merger and snapshot during which, or immediately following, this merger
event occurred.

The total number of BH merger events in Illustris TNG300 simulation is 590328.
First, we exclude all multiple mergers, that is all BHs that participate in
more than one merger event per snapshot. In a simplified case during triple
mergers, the lowest mass BH will be ejected while the remaining BHs will form
a binary system. However, the dynamic of multiple BH systems requires complex
treatment which is beyond the scope of this paper. Even though this excludes
144048 merger pairs ($\sim~{}24\%$), our main goal is to obtain the
comparative statistics between analytical and numerical predictions. Thus, we
limit our analysis to a smaller sample of single BH mergers for which we have
well-studied analytical and numerical models.

The next step was to link all of the merging BHs from BH mergers supplementary
catalog with the associated subhaloes, i.e. the host galaxies. First, we
identify the host galaxies of each BH participating in the merger event at the
snapshot preceding the merger, while BHs still occupy separate systems. Next,
we follow backward the evolution of each host galaxy, from the snapshot
preceding the BH merger to the snapshot when the galaxies first appeared in
the simulation. This step is necessary in order to accurately determine the
mass ratio of the merging galaxies, since in the snapshot preceding the BH
merger host galaxies have already exchanged mass as a consequence of galaxy
interaction. The mass ratio of the merging galaxies is calculated at the
snapshot where the secondary, i.e. less massive galaxy, has the maximum mass,
thus prior mass loss due to galaxy merger. That snapshot is taken as the last
snapshot before the galaxy interaction started. However, since a galaxy merger
is a prolonged process, it is possible to imagine a scenario where one of the
progenitor galaxies undergoes a new merger between the last snapshot before
galaxy interaction and the snapshot of the BH merger. That merger would change
the total mass and other properties of the progenitor galaxy. We limit our
analysis to single mergers, thus we eliminate from our sample not only
galaxies that undergo multiple mergers per snapshot but also galaxies that
experience multiple mergers during galaxy interaction.

The additional selection criterion is imposed to ensure that each host galaxy
is well resolved. Following Blecha (5) and Kelley (27) we require that each
progenitor galaxy, at the last snapshot before the galaxy interaction, has
$M_{\textrm{DM}}>10^{10}\>{\rm M_{\odot}}$ and $M_{\textrm{star}}>10^{8}\>{\rm
M_{\odot}}$. In addition, we also limit the upper mass of the host galaxy to
$M_{\textrm{DM}}<8\times 10^{13}\>{\rm M_{\odot}}$ to exclude galaxy clusters
(Paul (45)).

As we will discuss in the following subsection, smole19 (54) showed that
escape velocities from minor merger remnants in analytical and numerical
models do not differ significantly. In addition, since BH mass is expected to
scale with the host galaxy mass, minor mergers would lead to the formation of
low mass ratio BH binaries, and the expected kick velocity in such systems is
not enough to eject the central BH (see the following subsections). Thus we
limit our analysis to major merger remnants. Our final sample contains 46031
major merger remnants or $\sim 8~{}\%$ of the total number of all BH mergers
in TNG300 simulation.

###  2.2 BH escape velocities in analytical and numerical galaxy models

smole19 (54) investigated differences between recoiling BHs in analytical and
numerical potential using isolated N-body simulations performed with GADGET-2
code gadget (55). The authors constructed numerical models for a set of
progenitor galaxies with different total masses and mass distributions within
galaxy components and simulated galaxy mergers with various galaxy mass
ratios. In addition, the authors constructed analytical models of merger
remnant galaxies with the same characteristics as their numerical
counterparts. The main goal was to compare escape velocities from numerical
and analytical merger remnants. In numerical models, BH was represented as one
massive particle placed at the merger remnant center and its trajectory was
followed directly from the simulation. In analytical models BH trajectory was
numerically integrated, using leapfrog integration. Both analytical and
numerical techniques result in the same value of escape velocity for isolated
and virialized galaxies. However, galaxy mergers cause mass redistribution
within merger remnants which reduces escape velocities in numerical models.
During galaxy mergers energy of individual particles is not conserved. In the
process called violent relaxation Lynden-Bell (35) weakly bound particles can
escape the host galaxy potential, which results in mass loss during mergers.
This process is not depicted in static analytical potential, making the
evolving numerical model a more realistic description of dynamical processes
in galaxies with merging BHs. As a consequence, smole19 (54) showed that
analytical models overestimate BH escape velocities. BH escape velocities in
numerical major merger remnant galaxies can be up to 25 $\%$ lower compared to
those in analytical models. However, this effect is predominantly limited to
major mergers. During minor mergers, the secondary galaxy will experience
significant mass redistribution, but its contribution to the total potential
is not sufficient to influence escape velocities, leading to the similar
escape velocities in numerical and analytical minor merger remnants.

Here, we use the results of the above work to make comparative statistics
between recoiling BHs in TNG300 simulation, assuming analytical versus
numerical models. We extract properties of merger remnant galaxies directly
from the simulation. For the given merger remnant mass, its central BH mass,
and the mass ratio of the merging galaxies, analytical and numerical escape
velocities calculated by smole19 (54) are extrapolated and assigned to each
galaxy from our TNG300 sample.

###  2.3 BH kick velocities

After assigning escape velocities to each major merger remnant galaxy, the
next step is to calculate kick velocities for each newly formed BH, adopting
the method described by micic11 (41). Kick velocity depends on the merging BH
mass ratio, their spin amplitude, and the alignment to the orbital angular
momentum:

$V_{\mathrm{k}}=[(V_{m}+V_{\perp}\rm{cos}{\xi})^{2}+(V_{\perp}\sin{\xi})^{2}+(V_{\parallel})^{2}]^{1/2},$
(2a)

where

$V_{m}=A\frac{q^{2}\left(1-q\right)}{\left(1+q\right)^{5}}\left[1+B\frac{q}{\left(1+q\right)^{2}}\right],$
(2b)
$V_{\perp}=H\frac{q^{2}}{\left(1+q\right)^{5}}\left(\alpha_{2}^{\parallel}-q\alpha_{1}^{\parallel}\right),$
(2c)

and

$V_{\parallel}=K\cos\left(\Theta-\Theta_{0}\right)\frac{q^{2}}{\left(1+q\right)^{5}}\left(\alpha_{2}^{\perp}-q\alpha_{1}^{\perp}\right).$
(2d)

The fitting constants are A$=1.2\times 10^{4}\>{\rm~{}km}\,{\rm s}^{-1}$,
B$=-0.93$, H$=(7.3\pm 0.3)\times 10^{3}\>{\rm~{}km}\,{\rm s}^{-1}$ and
K$\cos\left(\Theta-\Theta_{0}\right)=(6.0\pm 0.1)\times 10^{4}\;$; $q$ is the
mass ratio of the merging BHs, $\alpha_{i=1,2}=\frac{S_{i}}{M_{i}}$ is the
reduced spin parameter, $S_{i}$ is the spin angular momentum of BH, and the
orientation of the merger is determined with angles $\Theta$ i $\xi$.
Following micic11 (41) we calculate kick velocities for two distinct spin
distributions. In both models, BH spin amplitudes are taken from a uniform
distribution, while BH spin orientation depends on a model choice. In the
first model, BH spin orientations are also taken from a uniform distribution
(’random’ model), while in the second model, BH spins are always aligned with
the orbital angular moment of the binary (’aligned’ model), which results in
lower kick velocities. The mass ratio of the merging BHs is taken directly
from the BH merger supplementary catalog. BH kick velocities are sampled from
one of those distributions using 10000 realizations for each merging BH pair
from our sample.

###  2.4 Escaped BHs versus Offset AGN statistics

Next, we compare the sampled kick velocities to the escape velocity from the
merger remnant. Throughout this work we will refer to realizations with
$v_{\textrm{k, i}}>v_{\textrm{esc}}$ ($i\in[1,...,10000]$) as escaped BHs. The
escape probability of a host galaxy is defined as the fraction of escaped BHs
out of 10000 realizations. Recoiling BHs with $v_{\textrm{k,
i}}<v_{\textrm{esc}}$ will not be ejected from the galaxy, but they can spend
an extended period of time oscillating around the galaxy center. We refer to
those as BHs on bound orbits. We use the term offset AGN to refer to both
escaped BHs and BHs on bound orbits, thus all BHs with kick velocities large
enough to displace BH from the center of the galaxy.

In this work, we separately investigate the statistics of escaped BHs and
offset AGN. Our main goal is to calculate the comparative statistics between
analytical and numerical predictions on cosmological scales, using a sample of
well-resolved major merger remnants from TNG300 simulation, that reside
outside of galaxy clusters.

##  3 Results

###  3.1 Major merger remnants in TNG300

Figure 1: Total mass of progenitor galaxy as a function of the central BH
mass.

Fig. 1 shows the total mass of each progenitor galaxy as a function of the
central BH mass. The color bar indicates the estimated 2D density of data
points. Our sample is dominated by lower mass galaxies,
$M_{\textrm{gal}}<5\times 10^{11}\>{\rm M_{\odot}}$, that host central BHs
with masses $M_{\textrm{BH}}<10^{7}\>{\rm M_{\odot}}$.

Fig. 2 shows the distribution of merging BH mass ratio, divided by the total
BH mass, for all mergers (solid lines) and for mergers occurring at $z<1$
(dashed lines). A lower fraction of high mass ratios at $z<1$ indicates that
mergers of equal mass BHs were more common at high redshifts before multiple
mergers and accretion led to BH growth. The distribution of all BHs (blue
solid line), shows an increase toward equal mass mergers, while massive BHs
show a flatter distribution. Dotted lines show merging BH mass ratios in
Illustris simulation, calculated by Blecha (5). This distribution is flatter
and with a higher fraction of low mass ratios. Assuming that central BH mass
scales with the host galaxy mass, choosing only major galaxy mergers, as we
did in our sample, would lead to higher BH mass ratios. A lower fraction of
the most massive BHs reflects the imposed higher limit for total galaxy mass
of $M_{\textrm{DM}}<8\times 10^{13}\>{\rm M_{\odot}}$, so we exclude galaxy
clusters.

Figure 2: The distribution of merging BH mass ratio, for different total BH
mass bins. Dashed lines represent the same distribution calculated for BH
mergers at $z<1$. Dotted lines show the distribution in Illustris simulation,
calculated by Blecha et al. 2016.  Figure 3: Numerically calculated escape
velocity as a function of redshift.

BH mergers at high redshifts are expected to produce more escaped BHs due to
generally lower masses of merger remnant galaxies, and thus lower escape
velocities. Fig. 3 shows escape velocities as a function of redshift,
calculated using numerical models. As predicted by the hierarchical growth of
structures, the number of high-mass galaxies with high escape velocities
decreases with redshift as galaxies build their masses and thus potential
wells deep enough to retain a recoiling BH.

For the extracted galaxy mass, its central BH mass, and the mass ratio of the
merging galaxies we estimate the escape velocity for each galaxy in our sample
using the results obtained by smole19 (54). Fig. 4 shows escape velocities
from galaxies in analytical models as a function of their escape velocities in
numerical models. The color bar indicates the estimated 2D density of data
points on a logarithmic scale. The black dashed line represents a linear fit
to our data, while the solid red line denotes $v_{\textrm{esc
(an.)}}=v_{\textrm{esc (num.)}}$. The majority of galaxies from our sample
occupy the region above the red line, thus on a statistical level numerical
models of major merger remnants predict lower escape velocities compared to
analytical models.

Figure 4: Escape velocities from galaxies in analytical models as a function
of their escape velocities in numerical models. The black dashed line
represents a linear fit to our data and the solid red line denotes
$v_{\textrm{esc (an.)}}=v_{\textrm{esc (num.)}}$.

In this work, we adopted the method described by micic11 (41) to calculate
kick velocities of merged BHs. Fig. 5 shows the distribution of kick
velocities from our data set, as a function of the merging BH mass ratio.
Values of kick velocities are the result of 10000 realizations per merger
event, in which BH spins are sampled from random or from aligned distribution.
If BH spins prior BH merger are drawn from a uniform distribution (random
model, black vertical lines) the resulting kick velocities can be
$>2000~{}\textrm{km/s}$, which is enough to eject central BH even from the
most massive galaxies in our sample. On the other hand, if BH spins are always
aligned with the orbital angular moment of the binary (aligned model, red
horizontal lines), the highest possible kick velocities are $\sim
250~{}\textrm{km/s}$ and recoiling BH are expected to be rare. Both
distributions yield the highest kick velocities for mergers of BHs with
comparable masses.

Figure 5: Distribution of kick velocities as a function of the mass ratio of
merging BHs. Black vertical lines represent random spin model, while red
horizontal lines denote to aligned spin model.

###  3.2 Statistics of escaped BHs

Figure 6: Probability function for a BH to escape for random spin model of BH
kick velocity distribution. Different colors indicate numerical (solid line,
pink) or analytical (dashed line, cyan) values for escape velocity
calculation.

Next, we calculate the probability that each BH merger pair from our sample
produces a kick large enough for complete BH removal. For each merger event,
BH kick velocities are sampled from either random or aligned spin distribution
using 10000 realizations. Escape probability, or the probability for each
merger event to completely eject BH, is calculated by comparing each of the
sampled kick velocities $v_{\textrm{k, i}}$ ($i\in[1,...,10000]$) to the
escape velocity from the merger remnant. Realizations with $v_{\textrm{k,
i}}>v_{\textrm{esc}}$ are referred to as escaped BHs.

Fig. 6 shows BH escape probability adopting random spin model for BH kick
velocity distribution. Different colors correspond to different models used to
estimate the escape velocity from the merger remnant galaxy, numerical (solid
line, pink) and analytical (dashed, cyan). If the BH spin parameters are taken
from a uniform distribution, the maximal probability for escaped BH is $\sim
67~{}\%$ and $\sim 55~{}\%$ for numerical and analytical models, respectively.
Both models predict that the majority of merger events from our sample would
not lead to BH removal, however, both distributions show a local maximum value
around $\sim~{}40\%$ (numerical) and $\sim 35~{}\%$ (analytical). Thus, the
plot shows that the model with numerically estimated escape velocities results
in a higher probability for escaped BHs compared to the analytical model.

Under the assumption that BH spins are always aligned with the orbital angular
moment of the binary (aligned spin model), resulting kick velocities are lower
(Fig 5). In fact, our sample does not contain any escaped BHs for the
analytical model for escape velocity, while numerical model predicts the
highest probability of $\sim 2~{}\%$ for a merger event to remove BH. Thus, if
BH spin alignment during BH mergers is efficient, ejection of a newly formed
BH from its host potential well is not excepted.

Further, we explore how the probability that a BH merger will produce an
escaped BH evolves with time. Fig. 7 shows BH probability to escape for
different redshift bins. At redshifts $z>5$ escaped BH probability function
peaks close to 40 $\%$, followed by a bimodal distribution at redshifts
$1<z<4$, with the second peak at zero probability. At low redshifts the
function shows one peak and our sample is dominated by BH mergers that will
not lead to BH removal. Again, for each redshift bin escaped BH probability
function is shifted toward higher values when numerical models are used.

Figure 7: Same as Fig. 6, for different redshift bins.

Table 1 shows the total number of BH mergers that have escape probability
$>20~{}\%$ and $>40~{}\%$, assuming analytical versus numerical potential.
Random spin model predicts that over $50\%$ of BH mergers from our sample,
produce a BH with escape probability $>20~{}\%$ for both numerical and
analytical models. However, at higher escape probabilities the difference
between our models becomes significant. Numerical models predict that $25\%$
of BH mergers leave a remnant BH with escape probability $>40~{}\%$, compared
to $8~{}\%$ in analytical models. Under the assumption that BH spins are
aligned, escape probabilities $>2~{}\%$ are not expected.

Similarly, Table 2 shows the total number of BH mergers with escape
probability $>20~{}\%$ and $>40~{}\%$, per redshift bin. The results presented
here refer to random model of spin distribution. At each redshift bin,
numerical models produce a higher fraction of BH mergers with escape
probability exceeding $40~{}\%$. Moreover, the gap between the fraction of
escaped BHs predicted analytically versus numerically widens as redshift
increases. This is a consequence of generally more turbulent galaxy evolution
at higher redshifts, with numerous mergers occurring before galaxies acquire
sufficient mass to retain the majority of recoiled BHs. Numerical models
better capture these processes in merging galaxies and thus provide a more
accurate estimate of the number of escaped BHs. If an offset AGN is defined as
a kicked BH with escape probability $>40~{}\%$, numerical models predict a
substantially more significant number of BH ejected from the host galaxy
centre, compared to analytical models. However, the gap between the two models
becomes less pronounced when considering an escape probability $>20~{}\%$.

Table 1: Total number of escaped BHs spin | escape | total number of escaped BHs
---|---|---
model | probability | analytical | numerical
| >20% | 24456 (53%) | 26510 (57%)
random | >40% | 3741 (8%) | 11580 (25%)
| >20% | 0 | 0
aligned | >40% | 0 | 0
Table 2: Total number of escaped BHs per redshift bin redshift |  escape | total number of escaped BHs
---|---|---
| probability | analytical | numerical
$0<z\leq 1$ | >20% | 5064 (34%) | 6106 (40%)
| >40% | 257 (2%) | 1594 (11%)
$1<z\leq 2$ | >20% | 11410 (58%) | 12126 (62%)
| >40% | 1907 (10%) | 5680 (29%)
$2<z\leq 3$ | >20% | 6170 (69%) | 6404 (72%)
| >40% | 1324 (15%) | 3427 (39%)
$3<z\leq 4$ | >20% | 1442 (70%) | 1491 (72%)
| >40% | 214 (10%) | 715 (35%)
$4<z\leq 5$ | >20% | 308 (74%) | 321 (77%)
| >40% | 32 (7%) | 136 (33%)
$z>5$ | >20% | 39 (85%) | 39 (85%)
| >40% | 4 (9%) | 19 (41%)
Figure 8: Median escape probability as a function of redshift. Kick velocities
are calculated assuming random spin distribution.

Fig. 8 shows the median escape probability as a function of redshift, for
numerical (red) and analytical (blue) models. The median probability that a BH
merger pair would produce an escaped BH increases with redshift. At high
redshifts a greater percentage of BH mergers are expected to result in kick
velocities high enough to eject central BH.

###  3.3 Statistics of offset AGN

Figure 9: The distribution of $\log(v_{\textrm{k,median}}/v_{\textrm{esc}})$
ratio for random (left panel) and aligned (right panel) spin distributions.
Solid and dashed lines denote numerically and analytically calculated escape
velocities, respectively. Different colors indicate the total merging BH mass
bins. The dotted line represents the distribution obtained by Blecha et al.
(2016). For kick velocities $v_{\textrm{k,median}}<v_{\textrm{esc}}$ our
models predict a greater number of offset AGN.

Further, we explore the distribution of offset AGN, thus we also take into
account recoiling BHs with kick velocities lower than the escape velocity, but
potentially still high enough to produce a spatial offset. Since Blecha (5)
provided a similar distribution for the Illustris sample, we reproduce their
figure 4 and compare our results. Fig. 9 shows the distribution of median kick
velocity for a BH merger pair, scaled to the host escape velocity in numerical
(solid line) and analytical (dashed line) models. The distribution is
separated by total BH mass, indicated by different colors. Blue line
($\log(M_{\textrm{BH,tot}}/\>{\rm M_{\odot}})>6$) shows all BH mergers in our
sample. A greater fraction of escaped BHs is predicted for low-mass BHs, while
complete BH ejections of most massive BHs are not expected. However, the plot
shows the median kick velocity obtained after 10000 realizations. Assuming
random spin distribution even the most massive BHs occasionally get kick
velocities large enough for complete ejection from their hosts. If BH spins
are aligned, kick velocities
$\log(v_{\textrm{k,median}}/v_{\textrm{esc}})>-0.5$ are possible only for the
lowest mass BHs. Regardless of the total BH mass and the chosen model for BH
spin distribution, numerical models always predict a greater number of offset
AGN.

Dotted lines on Fig. 9 represent the distribution obtained by Blecha (5),
assuming random and aligned spin distributions. Our models predict a greater
fraction of BHs with kick velocities $v_{\textrm{k,median}}<v_{\textrm{esc}}$.
Recoiled BHs in range $\log v_{\textrm{k,
median}}/v_{\textrm{esc}}\in(-0.5,1)$ are expected to produce offset AGN that
can spend a significant time on bound orbits, thus our models predict a
greater number of offset AGN compared to the findings of Blecha (5). However,
we note that our results cannot be compared directly, since our sample
contains only major galaxy mergers from TNG300 simulation, while Blecha (5)
shows the distribution of galaxy mergers with all mass ratios from lower
volume cube of Illustris simulation. In fact, distributions calculated in this
work predict a higher fraction of offset AGN, using both analytical and
numerical estimates for the escape velocity. This indicates that the values of
kick velocities in our sample are higher. Alongside BH spin parameters, the
mass ratio of merging BHs plays a key role in the resulting kick velocity
amplitude. Fig. 5 shows that the highest kick velocities can be achieved in
mergers of BHs with comparable masses. In contrast, mergers of low mass ratio
binaries result in low kick velocities and escaped BHs in such systems are
expected to be rare. As shown in Fig. 2 the distribution of merging BH mass
ratios from our sample has a lower fraction of low mass ratios compared to
Blecha (5) sample. By including only major mergers in our sample, we have
introduced a preference for higher kick velocities. Again, a lower fraction of
most massive BHs in our sample reflects the imposed upper limit for host
galaxy mass of $M_{\textrm{DM}}<8\times 10^{13}\>{\rm M_{\odot}}$.

Next, we explore the maximal distance from a galaxy center reached by recoiled
BH, assuming the median value of kick velocity. Fig. 10 shows the maximal
separation from a galaxy center as a function of $v_{\textrm{k,
median}}/v_{\textrm{esc}}$ ratio. Maximal distances from the host center that
BHs could reach in our analytical models (left panel) are in agreement with
Blecha (5) predictions (red lines). Recoiled BHs from both TNG300 (this work)
and Illustris (Blecha (5)) samples will reach comparable separations if their
trajectories are calculated analytically. This agreement is anticipated as we
employ similar analytical models. In numerical models (right panel), due to
the violent relaxation process, redistribution of mass within the merger
remnant galaxies will allow recoiled BHs to reach greater separations. BHs on
larger galactocentric distances also spend more time on bound orbits outside
of the galactic center, thus for the given $v_{\textrm{k,
median}}/v_{\textrm{esc}}$ ratio numerical models predict more spatially
offset AGNs.

Table 3 shows the total number of recoiled BHs at separations $>5$ kpc and
$>20$ kpc per redshift bin, assuming numerical versus analytical models. For a
median value of kick velocity majority of recoiled BHs will reach separations
$>5$ kpc in numerical models, in contrast to analytical models where $\geq
90~{}\%$ BHs will stay bound to central galaxy regions. The total number of
offset AGN, as well as the difference between analytical and numerical models,
increases with redshift, which reflects the hierarchical growth of structures
and a greater number of low-mass galaxies in the early Universe.

Figure 10: Maximal separation from a galaxy center reached by a recoiled BH as a function of $v_{\textrm{k, median}}/v_{\textrm{esc}}$ ratio for analytical (left) and numerical (right) models. Red lines represent the upper and lower limits of the distribution obtained by Blecha et al. (2016).  Table 3: Total number of offset AGN redshift |  $R_{\textrm{max}}$ (kpc) | analytical | numerical
---|---|---|---
all $z$ | >5 | 1550 (3%) | 26722 (58%)
| >20 | 208 (0.5%) | 16408 (36%)
$0<z{\scriptstyle{\leq}}1$ | >5 | 78 (0.5%) | 6094 (41%)
| >20 | 9 (0.06%) | 2695 (18%)
$1<z{\scriptstyle{\leq}}2$ | >5 | 777 (4%) | 12246 (68%)
| >20 | 103 (0.5%) | 7899 (40%)
$2<z{\scriptstyle{\leq}}3$ | >5 | 595 (7%) | 6465 (73%)
| >20 | 77 (0.9%) | 4568 (51%)
$3<z{\scriptstyle{\leq}}4$ | >5 | 86 (4%) | 1524 (74%)
| >20 | 18 (0.9%) | 1000 (49%)
$4<z{\scriptstyle{\leq}}5$ | >5 | 10 (3%) | 329 (79%)
| >20 | 1 (0.2%) | 202 (49%)
$z>5$ | >5 | 2 (4%) | 41 (90%)
| >20 | 0 | 31 (67%)

##  4 Conclusions

Static analytical models of merger remnant galaxies have been widely used to
study recoiling BHs (e.g. Blecha (5)). However, smole19 (54) suggested that
numerical models represent more realistic approach since galaxy mergers can
lead to a decrease in galaxy mass via violent relaxation process that is not
depicted in static analytical models. The authors showed that escape
velocities in numerical models can be up to 25 $\%$ lower compared to those in
analytical models.

Here, we extended the above work and investigated the comparative statistics
between recoiling BHs in analytical and numerical models of galaxies extracted
from TNG300 simulation. Our sample consisted of 46031 well-resolved major
merger remnant galaxies outside galaxy clusters. For the given merger remnant
galaxy mass, central BH mass, and mass ratio of progenitor galaxies we
estimated analytical and numerical escape velocities, extrapolating the escape
velocities calculated by smole19 (54). Kick velocities are calculated using
the method described by micic11 (41). For each merging BH pair from our
sample, BH spins are sampled from random or aligned spin distributions, using
10000 realizations. Comparing kick velocities to the escape velocity from the
host galaxy, we are able to distinguish between recoiled BHs on bound orbits
and escaped BHs.

If BH spins are taken from random spin distribution, numerical models predict
that 25 $\%$ of merging BHs from our sample will have escape probability
$>40~{}\%$, compared to $8~{}\%$ for analytical models. On the other hand,
aligned spin distribution does not yield kick velocities large enough to
produce escaped BHs. As redshift increases, the disparity between analytical
and numerical models becomes more prominent, highlighting the turbulent
evolution of galaxies at high redshifts. Processes in merging galaxies cannot
be described using static analytical models alone, thus numerical models
provide a better estimate of the number of escaped BHs.

High escape probabilities predicted by numerical models can have a negative
influence on BH growth through mergers. Major merges of gas reach galaxies can
fuel gas into central galaxy regions, trigger episodes of gas accretion and
exponential BH growth. At redshifts $z>5$, $\sim 40~{}\%$ of major mergers in
numerical models will produce a recoiling BH with escape probability
$>40~{}\%$. This makes merger-driven BH growth at high redshifts challenging.
However, if a kicked BH remains bound, future mergers are less likely to eject
it since escape probability decreases at lower redshifts, as galaxies build
their mass. For aligned BH spin model, kick amplitudes are not high enough to
halt BH growth through mergers, since neither numerical nor analytical models
produce recoiling BHs with escape probability $>2~{}\%$.

Numerical models also predict a greater number of offset AGN. For any given
fraction of the escape velocity, recoiled BHs in numerical models are able to
reach greater distances from the host center, increasing the probability of
their detection. Assuming median values of kick velocity numerical models
predict that 58 $\%$ of BH mergers in our sample would produce kicks large
enough to displace BHs at separations $>5$ kpc in numerical models, while
analytical models predict only 3 $\%$ of BHs at those separations. Recoiled
BHs with larger spatial offsets also have longer return times, which prolongs
the period during which they could possibly be observed.

Thus, on cosmological scales, numerically modeled recoiling BHs have a higher
escape probability and predict a greater number of offset AGN. We conclude
that numerical models should be favored over analytical ones since BH
ejections take place in non-virialized merger remnants.

#### Funding Statement

This research was supported by the Ministry of Science, Technological
Development and Innovation of the Republic of Serbia (MSTDIRS) through
contract no. 451-03-47/2023-01/200002 made with Astronomical Observatory of
Belgrade. This research was partially supported by the computing
infrastructure acquired in the frame of the project BOWIE (PROMIS 6060916,
Science Fund of the Republic of Serbia).

#### Competing Interests

None

#### Data Availability Statement

The data generated in this study are available from the corresponding author,
M.S., upon reasonable request.

## References

  * (1) Pau Amaro-Seoane et al.  “Astrophysics with the Laser Interferometer Space Antenna”  In _arXiv e-prints_ , 2022, pp. arXiv:2203.06016  DOI: 10.48550/arXiv.2203.06016
  * (2) M.. Begelman, R.. Blandford and M.. Rees  “Massive black hole binaries in active galactic nuclei”  In _Nature_ 287.5780, 1980, pp. 307–309  DOI: 10.1038/287307a0
  * (3) Laura Blecha, Thomas J. Cox, Abraham Loeb and Lars Hernquist  “Recoiling black holes in merging galaxies: relationship to active galactic nucleus lifetimes, starbursts and the $M_{\textrm{BH}}$-$\sigma$ relation”  In _Monthly Notices of the RAS_ 412.4, 2011, pp. 2154–2182  DOI: 10.1111/j.1365-2966.2010.18042.x
  * (4) Laura Blecha and Abraham Loeb  “Effects of gravitational-wave recoil on the dynamics and growth of supermassive black holes”  In _Monthly Notices of the RAS_ 390.4, 2008, pp. 1311–1325  DOI: 10.1111/j.1365-2966.2008.13790.x
  * (5) Laura Blecha et al.  “Recoiling black holes: prospects for detection and implications of spin alignment”  In _Monthly Notices of the RAS_ 456.1, 2016, pp. 961–989  DOI: 10.1093/mnras/stv2646
  * (6) Shoshannah Byrne-Mamahit et al.  “Interacting galaxies in the IllustrisTNG simulations - IV: enhanced supermassive black hole accretion rates in post-merger galaxies”  In _Monthly Notices of the RAS_ 519.4, 2022, pp. 4966–4981  DOI: 10.1093/mnras/stac3674
  * (7) Simone Callegari et al.  “GROWING MASSIVE BLACK HOLE PAIRS IN MINOR MERGERS OF DISK GALAXIES”  In _The Astrophysical Journal_ 729.2 The American Astronomical Society, 2011, pp. 85  DOI: 10.1088/0004-637X/729/2/85
  * (8) Simone Callegari et al.  “PAIRING OF SUPERMASSIVE BLACK HOLES IN UNEQUAL-MASS GALAXY MERGERS”  In _The Astrophysical Journal_ 696.1 The American Astronomical Society, 2009, pp. L89  DOI: 10.1088/0004-637X/696/1/L89
  * (9) Manuela Campanelli, Carlos O. Lousto, Yosef Zlochower and David Merritt  “Maximum Gravitational Recoil”  In _Physical Review Letters_ 98 American Physical Society, 2007, pp. 231102  DOI: 10.1103/PhysRevLett.98.231102
  * (10) M. Chiaberge et al.  “The puzzling case of the radio-loud QSO 3C 186: a gravitational wave recoiling black hole in a young radio source?”  In _Astronomy and Astrophysics_ 600, 2017, pp. A57  DOI: 10.1051/0004-6361/201629522
  * (11) F. Civano et al.  “Chandra High-resolution observations of CID-42, a Candidate Recoiling Supermassive Black Hole”  In _Astrophysical Journal_ 752.1, 2012, pp. 49  DOI: 10.1088/0004-637X/752/1/49
  * (12) F. Civano et al.  “A Runaway Black Hole in COSMOS: Gravitational Wave or Slingshot Recoil?”  In _Astrophysical Journal_ 717.1, 2010, pp. 209–222  DOI: 10.1088/0004-637X/717/1/209
  * (13) T.. Cox et al.  “The effect of galaxy mass ratio on merger-driven starbursts”  In _Monthly Notices of the RAS_ 384.1, 2008, pp. 386–409  DOI: 10.1111/j.1365-2966.2007.12730.x
  * (14) L. del Valle et al.  “Supermassive Black Holes in a Star-forming Gaseous Circumnuclear Disk”  In _Astrophysical Journal_ 811.1, 2015, pp. 59  DOI: 10.1088/0004-637X/811/1/59
  * (15) Sara L. Ellison et al.  “A definitive merger-AGN connection at z $\sim$ 0 with CFIS: mergers have an excess of AGN and AGN hosts are more frequently disturbed”  In _Monthly Notices of the RAS_ 487.2, 2019, pp. 2491–2504  DOI: 10.1093/mnras/stz1431
  * (16) Michael Eracleous, Todd A. Boroson, Jules P. Halpern and Jia Liu  “A Large Systematic Search for Close Supermassive Binary and Rapidly Recoiling Black Holes”  In _Astrophysical Journals_ 201.2, 2012, pp. 23  DOI: 10.1088/0067-0049/201/2/23
  * (17) Andrés Escala, Richard B. Larson, Paolo S. Coppi and Diego Mardones  “The Role of Gas in the Merging of Massive Black Holes in Galactic Nuclei. I. Black Hole Merging in a Spherical Gas Cloud”  In _Astrophysical Journal_ 607.2, 2004, pp. 765–777  DOI: 10.1086/386278
  * (18) Andrés Escala, Richard B. Larson, Paolo S. Coppi and Diego Mardones  “The Role of Gas in the Merging of Massive Black Holes in Galactic Nuclei. II. Black Hole Merging in a Nuclear Gas Disk”  In _Astrophysical Journal_ 630.1, 2005, pp. 152–166  DOI: 10.1086/431747
  * (19) Shy Genel et al.  “Introducing the Illustris project: the evolution of galaxy populations across cosmic time”  In _Monthly Notices of the RAS_ 445.1, 2014, pp. 175–200  DOI: 10.1093/mnras/stu1654
  * (20) Felipe G. Goicovic, Alberto Sesana, Jorge Cuadra and Federico Stasyszyn  “Infalling clouds on to supermassive black hole binaries - II. Binary evolution and the final parsec problem”  In _Monthly Notices of the RAS_ 472.1, 2017, pp. 514–531  DOI: 10.1093/mnras/stx1996
  * (21) José A. González et al.  “Supermassive Recoil Velocities for Binary Black-Hole Mergers with Antialigned Spins”  In _Physical Review Letters_ 98.23, 2007, pp. 231101  DOI: 10.1103/PhysRevLett.98.231101
  * (22) Alessia Gualandris and David Merritt  “Ejection of Supermassive Black Holes from Galaxy Cores”  In _Astrophysical Journal_ 678.2, 2008, pp. 780–797  DOI: 10.1086/586877
  * (23) Javiera Guedes, Piero Madau, Lucio Mayer and Simone Callegari  “RECOILING MASSIVE BLACK HOLES IN GAS-RICH GALAXY MERGERS”  In _The Astrophysical Journal_ 729.2 The American Astronomical Society, 2011, pp. 125  DOI: 10.1088/0004-637X/729/2/125
  * (24) Zoltán Haiman  “Constraints from Gravitational Recoil on the Growth of Supermassive Black Holes at High Redshift”  In _The Astrophysical Journal_ 613.1, 2004, pp. 36  DOI: 10.1086/422910
  * (25) P.. Jonker et al.  “A bright off-nuclear X-ray source: a type IIn supernova, a bright ULX or a recoiling supermassive black hole in CXO J122518.6+144545”  In _Monthly Notices of the RAS_ 407.1, 2010, pp. 645–650  DOI: 10.1111/j.1365-2966.2010.16943.x
  * (26) E. Kalfountzou, M. Santos Lleo and M. Trichas  “SDSS J1056+5516: A Triple AGN or an SMBH Recoil Candidate?”  In _Astrophysical Journall_ 851.1, 2017, pp. L15  DOI: 10.3847/2041-8213/aa9b2d
  * (27) Luke Zoltan Kelley, Laura Blecha and Lars Hernquist  “Massive black hole binary mergers in dynamical galactic environments”  In _Monthly Notices of the RAS_ 464.3, 2017, pp. 3131–3157  DOI: 10.1093/mnras/stw2452
  * (28) Fazeel Mahmood Khan et al.  “SWIFT COALESCENCE OF SUPERMASSIVE BLACK HOLES IN COSMOLOGICAL MERGERS OF MASSIVE GALAXIES”  In _The Astrophysical Journal_ 828.2 The American Astronomical Society, 2016, pp. 73  DOI: 10.3847/0004-637X/828/2/73
  * (29) S. Komossa and David Merritt  “Gravitational Wave Recoil Oscillations of Black Holes: Implications for Unified Models of Active Galactic Nuclei”  In _The Astrophysical Journal_ 689.2, 2008, pp. L89  DOI: 10.1086/595883
  * (30) S. Komossa, H. Zhou and H. Lu  “A Recoiling Supermassive Black Hole in the Quasar SDSS J092712.65+294344.0?”  In _Astrophysical Journall_ 678.2, 2008, pp. L81  DOI: 10.1086/588656
  * (31) Michael Koss et al.  “SDSS1133: an unusually persistent transient in a nearby dwarf galaxy”  In _Monthly Notices of the RAS_ 445.1, 2014, pp. 515–527  DOI: 10.1093/mnras/stu1673
  * (32) D.. Lambas, S. Alonso, V. Mesa and A.. O’Mill  “Galaxy interactions. I. Major and minor mergers”  In _Astronomy and Astrophysics_ 539, 2012, pp. A45  DOI: 10.1051/0004-6361/201117900
  * (33) Abraham Loeb  “Observable Signatures of a Black Hole Ejected by Gravitational-Radiation Recoil in a Galaxy Merger”  In _Physical Review Letters_ 99 American Physical Society, 2007, pp. 041103  DOI: 10.1103/PhysRevLett.99.041103
  * (34) Carlos O. Lousto and Yosef Zlochower  “Hangup Kicks: Still Larger Recoils by Partial Spin-Orbit Alignment of Black-Hole Binaries”  In _Physical Review Letters_ 107 American Physical Society, 2011, pp. 231102  DOI: 10.1103/PhysRevLett.107.231102
  * (35) D. Lynden-Bell  “Statistical Mechanics of Violent Relaxation in Stellar Systems”  In _Monthly Notices of the RAS_ 136.1, 1967, pp. 101–121  DOI: 10.1093/mnras/136.1.101
  * (36) Piero Madau and Eliot Quataert  “The Effect of Gravitational-Wave Recoil on the Demography of Massive Black Holes”  In _Astrophysical Journall_ 606.1, 2004, pp. L17–L20  DOI: 10.1086/421017
  * (37) Federico Marinacci et al.  “First results from the IllustrisTNG simulations: radio haloes and magnetic fields”  In _Monthly Notices of the RAS_ 480.4, 2018, pp. 5113–5139  DOI: 10.1093/mnras/sty2206
  * (38) K. Markakis et al.  “Subaru and e-Merlin observations of NGC 3718. Diaries of a supermassive black hole recoil?”  In _Astronomy and Astrophysics_ 580, 2015, pp. A11  DOI: 10.1051/0004-6361/201425077
  * (39) L. Mayer et al.  “Rapid Formation of Supermassive Black Hole Binaries in Galaxy Mergers with Gas”  In _Science_ 316.5833, 2007, pp. 1874  DOI: 10.1126/science.1141858
  * (40) David Merritt et al.  “Consequences of Gravitational Radiation Recoil”  In _The Astrophysical Journal_ 607.1, 2004, pp. L9  DOI: 10.1086/421551
  * (41) Miroslav Micic, Kelly Holley-Bockelmann and Steinn Sigurdsson  “Growing massive black holes in a Local Group environment: the central supermassive, slowly sinking and ejected populations”  In _Monthly Notices of the RAS_ 414.2, 2011, pp. 1127–1144  DOI: 10.1111/j.1365-2966.2011.18444.x
  * (42) Jill P. Naiman et al.  “First results from the IllustrisTNG simulations: a tale of two elements - chemical evolution of magnesium and europium”  In _Monthly Notices of the RAS_ 477.1, 2018, pp. 1206–1224  DOI: 10.1093/mnras/sty618
  * (43) Dylan Nelson et al.  “First results from the IllustrisTNG simulations: the galaxy colour bimodality”  In _Monthly Notices of the RAS_ 475.1, 2018, pp. 624–647  DOI: 10.1093/mnras/stx3040
  * (44) David R. Patton et al.  “Galaxy pairs in the Sloan Digital Sky Survey - VI. The orbital extent of enhanced star formation in interacting galaxies”  In _Monthly Notices of the RAS_ 433.1, 2013, pp. L59–L63  DOI: 10.1093/mnrasl/slt058
  * (45) S. Paul, R.. John, P. Gupta and H. Kumar  “Understanding ‘galaxy groups’ as a unique structure in the universe”  In _Monthly Notices of the RAS_ 471.1, 2017, pp. 2–11  DOI: 10.1093/mnras/stx1488
  * (46) Annalisa Pillepich et al.  “First results from the IllustrisTNG simulations: the stellar mass content of groups and clusters of galaxies”  In _Monthly Notices of the RAS_ 475.1, 2018, pp. 648–675  DOI: 10.1093/mnras/stx3112
  * (47) Planck Collaboration et al.  “Planck 2015 results. XIII. Cosmological parameters”  In _Astronomy and Astrophysics_ 594, 2016, pp. A13  DOI: 10.1051/0004-6361/201525830
  * (48) I.. Redmount and M.. Rees  “Gravitational-radiation rocket effects and galactic structure.”  In _Comments on Astrophysics_ 14, 1989, pp. 165
  * (49) Andrew Robinson et al.  “SPECTROPOLARIMETRIC EVIDENCE FOR A KICKED SUPERMASSIVE BLACK HOLE IN THE QUASAR E1821+643”  In _The Astrophysical Journal Letters_ 717.2 The American Astronomical Society, 2010, pp. L122  DOI: 10.1088/2041-8205/717/2/L122
  * (50) Francisco Rodríguez Montero et al.  “Mergers, starbursts, and quenching in the SIMBA simulation”  In _Monthly Notices of the RAS_ 490.2, 2019, pp. 2139–2154  DOI: 10.1093/mnras/stz2580
  * (51) Rok Roškar et al.  “Orbital decay of supermassive black hole binaries in clumpy multiphase merger remnants”  In _Monthly Notices of the RAS_ 449.1, 2015, pp. 494–505  DOI: 10.1093/mnras/stv312
  * (52) Debora Sijacki, Volker Springel and Martin G. Haehnelt  “Gravitational recoils of supermassive black holes in hydrodynamical simulations of gas-rich galaxies”  In _Monthly Notices of the RAS_ 414.4, 2011, pp. 3656–3670  DOI: 10.1111/j.1365-2966.2011.18666.x
  * (53) Debora Sijacki et al.  “The Illustris simulation: the evolving population of black holes across cosmic time”  In _Monthly Notices of the RAS_ 452.1, 2015, pp. 575–596  DOI: 10.1093/mnras/stv1340
  * (54) Majda Smole, Miroslav Micic, Mitrašinović and Ana  “Recoiling supermassive black holes in analytical and numerical galaxy potential”  In _Monthly Notices of the RAS_ 488.4, 2019, pp. 5566–5579  DOI: 10.1093/mnras/stz2107
  * (55) Volker Springel  “The cosmological simulation code GADGET-2”  In _Monthly Notices of the RAS_ 364.4, 2005, pp. 1105–1134  DOI: 10.1111/j.1365-2966.2005.09655.x
  * (56) Volker Springel  “E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh”  In _Monthly Notices of the RAS_ 401.2, 2010, pp. 791–851  DOI: 10.1111/j.1365-2966.2009.15715.x
  * (57) Volker Springel et al.  “First results from the IllustrisTNG simulations: matter and galaxy clustering”  In _Monthly Notices of the RAS_ 475.1, 2018, pp. 676–698  DOI: 10.1093/mnras/stx3304
  * (58) Aaron Stemo et al.  “A Catalog of AGN Host Galaxies Observed with HST/ACS: Correlations between Star Formation and AGN Activity”  In _Astrophysical Journal_ 888.2, 2020, pp. 78  DOI: 10.3847/1538-4357/ab5f66
  * (59) Aaron Stemo et al.  “A Catalog of 204 Offset and Dual Active Galactic Nuclei (AGNs): Increased AGN Activation in Major Mergers and Separations under 4 kpc”  In _Astrophysical Journal_ 923.1, 2021, pp. 36  DOI: 10.3847/1538-4357/ac0bbf
  * (60) M Tremmel et al.  “Dancing to changa: a self-consistent prediction for close SMBH pair formation time-scales following galaxy mergers”  In _Monthly Notices of the RAS_ 475.4, 2018, pp. 4967–4977  DOI: 10.1093/mnras/sty139
  * (61) Marta Volonteri  “Gravitational Recoil: Signatures on the Massive Black Hole Population”  In _Astrophysical Journall_ 663.1, 2007, pp. L5–L8  DOI: 10.1086/519525
  * (62) Rainer Weinberger et al.  “Supermassive black holes and their feedback effects in the IllustrisTNG simulation”  In _Monthly Notices of the RAS_ 479.3, 2018, pp. 4056–4072  DOI: 10.1093/mnras/sty1733
  * (63) Madalyn E. Weston et al.  “Incidence of WISE -selected obscured AGNs in major mergers and interactions from the SDSS”  In _Monthly Notices of the RAS_ 464.4, 2017, pp. 3882–3906  DOI: 10.1093/mnras/stw2620


# Pairwise Adjusted Mutual Information

Denys Lazarenko
Technische Universität München
Germany
<EMAIL_ADDRESS>&Thomas Bonald
Institut Polytechnique de Paris
France
<EMAIL_ADDRESS>
###### Abstract

A well-known metric for quantifying the similarity between two clusterings is
the adjusted mutual information. Compared to mutual information, a corrective
term based on random permutations of the labels is introduced, preventing two
clusterings being similar by chance. Unfortunately, this adjustment makes the
metric computationally expensive. In this paper, we propose a novel adjustment
based on pairwise label permutations instead of full label permutations.
Specifically, we consider permutations where only two samples, selected
uniformly at random, exchange their labels. We show that the corresponding
adjusted metric, which can be expressed explicitly, behaves similarly to the
standard adjusted mutual information for assessing the quality of a
clustering, while having a much lower time complexity. Both metrics are
compared in terms of quality and performance on experiments based on synthetic
and real data.

##  1 Introduction

A well-known metric for quantifying the similarity between two clusterings of
the same data is the adjusted mutual information [5, 10]. Compared to mutual
information, this metric is adjusted against chance, meaning that the
similarity cannot be due to randomness but only to the structure of the
dataset, appearing in both clusterings. This is the reason why this metric is
widely used in unsupervised learning, see [13, 9, 8, 12, 11] for various
applications.

The standard way of adjusting mutual information against chance is through
random label permutations of one of the clusterings [10]. Unfortunately, this
adjustment makes the metric computationally expensive. Specifically, the time
complexity of the metric is in $O(\max(k,l)n)$, where $k,l$ are the numbers of
clusters in each clustering and $n$ is the number of samples [6]. As a
comparison, the time complexity of mutual information is equal to $O(kl)$
given the contingency matrix of the clusterings, i.e., the matrix counting the
number of samples in each cluster pair. The additional computational effort
required by adjustment is significant as the number of samples $n$ is
typically much larger than the numbers of clusters $k,l$.

In this paper, we propose a novel adjustment based on pairwise permutations.
That is, we consider permutations where only two samples, selected uniformly
at random, exchange their labels. We show that the corresponding adjusted
metric, we refer to as pairwise adjusted mutual information, is as efficient
as adjusted mutual information for assessing the quality of a clustering, with
a much lower time complexity. In particular, the time complexity is the same
as that of mutual information. The gain in complexity is significant, as the
computation time is now independent of the number of samples $n$, given the
contingency matrix.

The rest of the paper is organized as follows. We first provide the definition
and key properties of adjusted mutual information in the general setting of
information theory. We then introduce mutual information with pairwise
adjustement and explain why the exact same properties are satisfied by this
new notion of adjusted mutual information. The application of both notions of
adjustment to clustering, including the explicit expressions of the
corresponding metrics, is presented in section 4. Experiments on both
synthetic and real data are presented in section 5. Section 6 concludes the
paper.

##  2 Adjusted mutual information

Let $P$ be the uniform probability measure on $\Omega=\\{1,\ldots,n\\}$, for
some positive integer $n$. Let $X,Y$ be random variables on the probability
space $(\Omega,P)$. Without any loss of generality, we assume that $X$ and $Y$
are mapping from $\Omega$ to sets consisting of consecutive integers, starting
from 1. Denoting by $H$ the entropy, the mutual information between $X$ and
$Y$ is defined by [1]:

$I(X,Y)=H(X)+H(Y)-H(X,Y).$ (1)

This is the information shared by $X$ and $Y$, which is equal to 0 if $X$ and
$Y$ are independent. A distance between $X$ and $Y$ can then be defined by:

$d(X,Y)=H(X,Y)-I(X,Y)=H(X|Y)+H(Y|X).$

This distance, known as the variation of information, is a metric in the
quotient space of random variables under the equivalence relation $X\sim Y$ if
and only if there is some bijection $\varphi$ such that $X=\varphi(Y)$ [4].

### Adjusted mutual information.

The adjusted mutual information between $X$ and $Y$, corresponding to the
mutual information between $X$ and $Y$ adjusted against chance, is defined by:

$\Delta I(X,Y)=I(X,Y)-{\rm E}(I(X,Y_{\sigma})),$ (2)

where $Y_{\sigma}$ is the random variable $Y\circ\sigma$, for any permutation
$\sigma$ of $\\{1,\ldots,n\\}$, and the expectation is taken over all
permutations $\sigma$, chosen uniformly at random.

######  Remark 1 (Normalization).

It is frequent to also normalize adjusted mutual information, so as to get a
score between 0 and 1 [10, 6]. In this paper, we only focus on the adjustment
step. Note that normalization can be equally applied to both considered
notions of adjustment and thus be studied separately.

We have the equivalent definition:

$\displaystyle\Delta I(X,Y)$ $\displaystyle={\rm E}(H(X,Y_{\sigma}))-H(X,Y),$
$\displaystyle=\frac{1}{2}({\rm E}(d(X,Y_{\sigma}))-d(X,Y)).$ (3)

This equivalence follows from Proposition 1 and the fact that the definition
is symmetric in $X$ and $Y$. All proofs are deferred to the appendix.

######  Proposition 1.

We have for any random variables $X$ and $Y$:

$\displaystyle H(X)$ $\displaystyle={\rm E}(H(X_{\sigma})),$
$\displaystyle{\rm E}(H(X,Y_{\sigma}))$ $\displaystyle={\rm
E}(H(X_{\sigma},Y)),$ $\displaystyle{\rm E}(I(X,Y_{\sigma}))$
$\displaystyle={\rm E}(I(X_{\sigma},Y)).$

In view of (3), we expect $\Delta I(X,Y)$ to be positive if $X$ and $Y$ share
information, as $X$ is expected to be closer to $Y$ (for the distance $d$)
than to $Y_{\sigma}$, a randomized version of $Y$. There are specific cases
where $\Delta I(X,Y)=0$, as stated in Proposition 2; these cases will be
interpreted in terms of clustering in section 4.

######  Proposition 2.

We have $\Delta I(X,Y)=0$ whenever $Y$ (or $X$, by symmetry) is constant or
equal to some permutation of $\\{1,\ldots,n\\}$.

### Adjusted entropy.

Observing that $H(X)=I(X,X)$, we define similarly the adjusted entropy of $X$
by:

$\Delta H(X)=\Delta I(X,X)=H(X)-{\rm E}(I(X,X_{\sigma})).$

By (1), we get:

$\Delta H(X)={\rm E}(H(X,X_{\sigma}))-H(X)=\frac{1}{2}{\rm
E}(d(X,X_{\sigma})).$ (4)

Since $d$ is a metric, this shows that the adjusted entropy of $X$ is non-
negative.

######  Proposition 3.

We have $\Delta H(X)=0$ if and only if $X$ is constant or equal to some
permutation of $\\{1,\ldots,n\\}$.

Proposition 3 characterizes random variables with zero adjusted entropy.
Again, this result will be interpreted in terms of clustering in section 4.

##  3 Pairwise adjustment

In this section, we introduce pairwise adjusted mutual information. The
definition is the same as adjusted mutual information, except that the
permutation $\sigma$ is now restricted to the set of pairwise permutations.
Specifically, we consider permutations $\sigma$ for which there exists
$i,j\in\\{1,\ldots,n\\}$ such that $\sigma(i)=j$ and $\sigma(j)=i$, whereas
$\sigma(t)=t$ for all $t\neq i,j$. We consider the set of such permutations
$\sigma$ where the samples $i,j$ are drawn uniformly at random in the set
$\\{1,\ldots,n\\}$. We denote by ${\sigma_{\rm p}}$ such a random permutation.
Observe that ${\sigma_{\rm p}}$ is the identity with probability $1/n$ (the
probability that $i=j$).

### Pairwise adjusted mutual information.

We define the pairwise adjusted mutual information as:

$\Delta_{\rm p}I(X,Y)=I(X,Y)-{\rm E}(I(X,Y_{\sigma_{\rm p}})).$

This is exactly the same definition as the adjusted mutual information, except
for the considered permutations ${\sigma_{\rm p}}$. It can be readily verified
that the same properties apply, with the exact same proofs, a key property
being that the random permutations ${\sigma_{\rm p}}$ and ${\sigma_{\rm
p}}^{-1}$ have the same distributions. In particular, we have the analogue of
(3):

$\displaystyle\Delta_{\rm p}I(X,Y)$ $\displaystyle={\rm E}(H(X,Y_{\sigma_{\rm
p}}))-H(X,Y),$ $\displaystyle=\frac{1}{2}({\rm E}(d(X,Y_{\sigma_{\rm
p}}))-d(X,Y)).$ (5)

Moreover, $\Delta_{\rm p}I(X,Y)=0$ whenever $X$ or $Y$ is constant or equal to
some permutation of $\\{1,\ldots,n\\}$.

### Pairwise adjusted entropy.

We also define the pairwise adjusted entropy as:

$\Delta_{\rm p}H(X)=\Delta_{\rm p}I(X,X)=H(X)-{\rm E}(I(X,X_{\sigma_{\rm
p}})).$

We have $\Delta_{\rm p}H(X)\geq 0$, with equality if and only if $X$ is
constant or equal to some permutation of $\\{1,\ldots,n\\}$.

##  4 Application to clustering

Let $A=\\{A_{1},\ldots,A_{k}\\}$ and $B=\\{B_{1},\ldots,B_{l}\\}$ be two
partitions of some finite set $\\{1,\ldots,n\\}$ into $k$ and $l$ clusters,
respectively. Let $\Omega=\\{1,\ldots,n\\}$ and ${\rm P}$ be the uniform
probability measure over $\Omega$. Consider the random variables $X$ and $Y$
defined on $(\Omega,{\rm P})$ by $X^{-1}(i)=A_{i}$ for all $i=1,\ldots,k$ and
$Y^{-1}(j)=B_{j}$ for all $j=1,\ldots,l$. Note that $X(\omega)$ and
$Y(\omega)$ can be interpreted as the labels $i$ and $j$ of sample $\omega$ in
clusterings $A$ and $B$, for each $\omega\in\\{1,\ldots,n\\}$.

We denote by $a_{i}=|A_{i}|$ the size of cluster $A_{i}$, by $b_{j}=|B_{j}|$
the size of cluster $B_{j}$, and by $n_{ij}=|A_{i}\cap B_{j}|$ the number of
samples both in cluster $A_{i}$ and cluster $B_{j}$, for all $i=1,\ldots,k$
and $j=1,\ldots,l$. The matrix $(n_{ij})_{1\leq i\leq k,1\leq j\leq l}$ is
known as the contingency matrix. Note that $a_{i}$ and $b_{j}$ are the
respective sums of row $i$ and column $j$ of the contingency matrix.

### Adjusted mutual information.

A well-known metric for assessing the similarity $s(A,B)$ between clusterings
$A$ and $B$ is the adjusted mutual information111Recall that we don’t
normalize the metric, see Remark 1. $\Delta I(X,Y)$ between the corresponding
random variables $X$ and $Y$. In words, this is the common information shared
by clusterings $A$ and $B$ not due to randomness.

By Proposition 2, we have $s(A,B)=0$ whenever clustering $A$ (or $B$, by
symmetry) is trivial, that is, it consists of a single cluster or of $n$
clusters (one per sample). This is a key property, showing the interest of the
adjustment.

It is known that [10]:

$\displaystyle\begin{split}&s(A,B)=-\sum_{i=1}^{k}\sum_{j=1}^{l}\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}+\sum_{i=1}^{k}\sum_{j=1}^{l}\sum_{c=(a_{i}+b_{j}-n)^{+}}^{\min(a_{i},b_{j})}\\\
&\frac{a_{i}!b_{j}!(n-a_{i})!(n-b_{j})!}{n!c!(a_{i}-c)!(b_{j}-c)!(n-a_{i}-b_{j}+c)!}\frac{c}{n}\log\frac{c}{n},\end{split}$
(6)

with the notation $(\cdot)^{+}=\max(\cdot,0)$. The time complexity of this
formula, which is dominated by the second term, is in $O(\max(k,l)n)$ [6]. In
particular, it is linear in the number of samples $n$.

Interestingly, we can similarly assess the quantity of information $q(A)$
contained in clustering $A$ through the adjusted entropy $\Delta H(X)$ of the
corresponding random variable $X$. This is the information contained in $A$
not due to randomness. We have $q(A)\geq 0$ and, by Proposition 3, $q(A)=0$ if
and only if clustering $A$ is trivial, that is, it consists of a single
cluster or of $n$ clusters (one per sample).

Since $q(A)=s(A,A)$, it follows from (6) that:

$\displaystyle\begin{split}q(A)&=-\sum_{i=1}^{k}\frac{a_{i}}{n}\log\frac{a_{i}}{n}+\sum_{i,j=1}^{K}\sum_{c=(a_{i}+a_{j}-n)^{+}}^{\min(a_{i},a_{j})}\\\
&\frac{a_{i}!a_{j}!(n-a_{i})!(n-a_{j})!}{n!c!(a_{i}-c)!(a_{j}-c)!(n-a_{i}-a_{j}+k)!}\frac{c}{n}\log\frac{c}{n}.\end{split}$

The time complexity of this formula, also dominated by the second term, is in
$O(kn)$. Again, this complexity is linear in the number of samples $n$.

### Pairwise adjusted mutual information.

The main contribution of the paper is the following new measure of similarity
$s_{\rm p}(A,B)$ between clusterings $A$ and $B$, based on the pairwise
adjusted mutual information $\Delta_{\rm p}I(X,Y)$ between the corresponding
random variables $X$ and $Y$. We have an explicit expression for this
similarity:

######  Theorem 1.

We have:

$\displaystyle\begin{split}s_{\rm
p}(A,B)&=2\sum_{i=1}^{k}\sum_{j=1}^{l}\frac{n_{ij}(n-a_{i}-b_{j}+n_{ij})}{n^{2}}\\\
&\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}-1}{n}\log\frac{n_{ij}-1}{n}\right)\\\
&+2\sum_{i=1}^{k}\sum_{j=1}^{l}\frac{(a_{i}-n_{ij})(b_{j}-n_{ij})}{n^{2}}\\\
&\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}+1}{n}\log\frac{n_{ij}+1}{n}\right).\end{split}$
(7)

The time complexity of this formula is in $O(kl)$, like mutual information. It
is independent of the number of samples $n$, given the contingency matrix.
Corollary 1 shows that the time complexity reduces to $O(m)$ the number of
non-zero entries of the contingency matrix, provided the latter is stored in
sparse format.

######  Corollary 1.

We have:

$\displaystyle\begin{split}s_{\rm
p}(A,B)&=2\sum_{i,j:n_{ij}>0}\frac{n_{ij}(n-a_{i}-b_{j}+n_{ij})}{n^{2}}\\\
&\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}-1}{n}\log\frac{n_{ij}-1}{n}\right)\\\
&+2\sum_{i,j:n_{ij}>0}\frac{(a_{i}-n_{ij})(b_{j}-n_{ij})}{n^{2}}\\\
&\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}+1}{n}\log\frac{n_{ij}+1}{n}+\frac{1}{n}\log\frac{1}{n}\right)\\\
&-2\left(n^{2}-\sum_{i=1}^{k}a_{i}^{2}-\sum_{j=1}^{l}b_{i}^{2}+\sum_{i,j:n_{ij}>0}n_{ij}^{2}\right)\\\
&\times\frac{1}{n}\log\frac{1}{n}.\end{split}$

Similarly, we can define the quantity of information $q_{\rm p}(A)$ in
clustering $A$ through the pairwise adjusted entropy $\Delta_{\rm p}H(X)$ of
the corresponding random variable $X$. Again, $q_{\rm p}(A)\geq 0$, with
$q_{\rm p}(A)=0$ if and only if clustering $A$ is trivial.

######  Corollary 2.

We have:

$\displaystyle\begin{split}q_{\rm
p}(A)&=2\sum_{i=1}^{k}\frac{a_{i}(n-a_{i})}{n^{2}}\\\
&\times\left(\frac{a_{i}}{n}\log\frac{a_{i}}{n}-\frac{a_{i}-1}{n}\log\frac{a_{i}-1}{n}-\frac{1}{n}\log\frac{1}{n}\right).\end{split}$

Note that the time complexity of this formula in $O(k)$. It only depends on
the number of clusters $k$, and not on the number of samples $n$.

##  5 Experiments

In this section, we compare both notions of adjusted mutual information
through experiments involving synthetic and real data. The experiments are run
on a computer equipped with an AMD Ryzen Threadripper 1950X 16-Core Processor
and 32 GB of RAM, with a a Debian 10 OS. All codes and datasets used in the
experiments are available online222See
https://github.com/denyslazarenko/Pairwise-Adjusted-Mutual-Information.

### Synthetic data.

We start with the simple case of $n=100$ samples with clusters of even sizes,
consisting of consecutive samples. Specifically, we consider the set of
clusterings $A^{(s)}$, consisting of clusters of size $s$ (except possibly the
last one), for $s=1,2,\ldots,100$. In particular, both $A^{(1)}$ and
$A^{(100)}$ are trivial clusterings while $A^{(5)}$ consists of 20 clusters of
size 5.

Figure 1 gives the similarity between clusterings $A^{(10)}$ and $A^{(s)}$
with respect to $s$ in terms of adjusted mutual information, for both notions
of adjustment. We observe very close behaviors, suggesting that both notions
of adjustment tend to capture the same patterns in the clusterings. Note that
the maximum similarity is attained for $s=10$ in both cases, as expected. The
similarity is equal to 0 for $s\in\\{1,100\\}$ for both cases, in agreement
with Proposition 2. We also observe local peaks at $s=20,30,\ldots,90$, which
can be interpreted by the fact that clustering $A^{(10)}$ is a refinement of
clustering $A^{(s)}$ for these values of $s$; similarly, the local peak at
$s=5$ may be interpreted by the fact that clustering $A^{(5)}$ is a refinement
of clustering $A^{(10)}$.

(a) Full adjustment

(b) Pairwise adjustment

Figure 1: Comparison of metrics on synthetic data ($n=100$).

We now consider random clusterings. Specifically, we assign $n$ samples to $k$
clusters independently at random, according to some probability distribution
$p=(p_{1},\ldots,p_{k})$, which is itself drawn at random333Namely, $p\propto
U$ where $U=(U_{1},\ldots,U_{k})$ is a vector of $k$ i.i.d. random variables
uniformly distributed over $[0,1]$.. Consider three such random clusterings
$A$, $B$, $C$ (with the same parameters $n$ and $k$, but different probability
distributions $p$). We would like to know whether $A$ is “closer” to $B$ or to
$C$. In particular, we are interested in testing whether both notions of
adjusted mutual information give the same ordering in the sense that:

$(s(A,B)-s(A,C))(s_{\rm p}(A,B)-s_{\rm p}(A,C))\geq 0.$ (8)

We compute the average precision score (fraction of triplets $A,B,C$ for which
(8) is true) over 1 000 independent samples of $A,B,C$, for different values
of $n$ and $k$. We repeat the experiment 100 times to get the mean and
standard deviation. The results are given in Table 1. We observe a very high
precision score, always higher than $93\%$, showing that both notions of
adjusted mutual information tend to give the same ordering of these random
clusterings.

$n$ | $k$ | Precision score
---|---|---
100 | 2 | $0.972\pm 0.004$
100 | 5 | $0.952\pm 0.007$
100 | 10 | $0.943\pm 0.006$
100 | 20 | $0.955\pm 0.008$
500 | 20 | $0.936\pm 0.007$
1000 | 20 | $0.933\pm 0.006$
1000 | 50 | $0.949\pm 0.008$
Table 1: Precision score (mean $\pm$ standard deviation)

For the performance gain, we compare the computation times of both versions of
adjusted mutual information for the similarity between clusterings $A$ and
$B$, where $A$ consists of $k=10$ clusters of same size and $B$ is a random
clustering, drawn as in the previous experiment. Both versions of adjusted
mutual information are coded in Python, with the standard version imported
from scikit-learn. Figure 2 shows the computation time when the number of
samples $n$ grows from $10^{2}$ to $10^{7}$. The performance gain brought by
pairwise adjustement is significant. In particular, the computation time
becomes independent of the number of samples.

Figure 2: Computation time with respect to $n$ (mean $\pm$ standard
deviation).

### Real data.

For real data, we consider the 79 datasets of the benchmark suite [2]444See
https://github.com/gagolews/clustering_benchmarks_v1.. We apply to each
dataset each of the following clustering algorithms:

  1. 1.

$k$-means

  2. 2.

Affinity propagation

  3. 3.

Mean shift

  4. 4.

Spectral clustering

  5. 5.

Ward

  6. 6.

Agglomerative clustering

  7. 7.

DBSCAN

  8. 8.

OPTICS

  9. 9.

Birch

  10. 10.

Gaussian Mixture

We use the scikit-learn implementation of these algorithms, with the
corresponding default parameters555See https://scikit-
learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html.666Dimension
reduction is applied to the MNIST datasets, consisting of 70 000 images of
size $28\times 28$ each, see the supplementary material for details.. We get
10 clusterings per dataset. The quality of each clustering is assessed through
the similarity with the available ground-truth labels, using adjusted mutual
information with either full adjustment or pairwise adjustment. We then
compute the Spearman correlation of the corresponding similarities, a value of
1 meaning the exact same ordering of the 10 clusterings with full adjustment
and pairwise adjustment. The results are shown in Figure 3, together with the
speed-up in computation time due to pairwise adjustment. In both cases, the 79
datasets are ordered by the number of samples, ranging from 105 to 105 600
[2].

(a) Spearman correlation.

(b) Speed-up of pairwise adjustment.

Figure 3: Comparison of metrics on real data.

We first observe that the correlation is very high, suggesting again that both
notions of adjusted mutual information tend to provide the same results. For
65 datasets among 79, the Spearman correlation is higher than 95%. As for the
computation time, we observe a significant performance gain, by one order of
magnitude for the largest datasets.

##  6 Conclusion

We have proposed another way of adjusting mutual information against chance,
through pairwise label permutations. The novel metric, whose explicit
expression is given in Theorem 1, has a much lower complexity than the usual
adjusted mutual information. Interestingly, both metrics can also be used to
assess the quantity of information contained in a clustering, which the common
property of being equal to 0 if and only if the clustering is trivial, as
stated in Proposition 3; again, the pairwise adjusted entropy, given in
Corollary 2, has a much lower complexity. Experiments on synthetic and real
data show that pairwise adjusted mutual information tends to provide the same
results as the usual adjusted mutual information for comparing clusterings,
while involving much less computations.

For future work, we plan to extend this idea to other similarity metrics.
While the practical interest is less obvious for the Adjusted Rand Index [3],
due to the fact that the time complexity of this metric is already independent
of the number of samples, it would be worth considering other versions of
information theoretic measures, as those studied in [7].

## Appendix

## Proof of Proposition 1

The first equality follows from the fact that $X$ and $X_{\sigma}$ have the
same distribution. Specifically, we have for any positive integer $k$,

$\displaystyle{\rm P}(X_{\sigma}=k)$ $\displaystyle={\rm
P}((X\circ\sigma)^{-1}(k)),$ $\displaystyle={\rm P}(\sigma^{-1}(X^{-1}(k))),$
$\displaystyle={\rm P}(X^{-1}(k)),$ $\displaystyle={\rm P}(X=k).$

For the second, we observe that if $\sigma$ is a random permutation of
$\\{1,\ldots,n\\}$, chosen uniformly at random, so is $\sigma^{-1}$ which
implies:

${\rm E}(H(X,Y_{\sigma}))={\rm E}(H(X,Y_{\sigma^{-1}}))={\rm
E}(H(X_{\sigma},Y)),$

where we have used the first equality and the fact that:

$(X,Y_{\sigma^{-1}})\circ\sigma=(X_{\sigma},Y).$

The third equality is a direct consequence of the two first.

## Proof of Proposition 2

If $Y$ is constant, then $Y=Y_{\sigma}$ for all permutations $\sigma$ and the
result follows from (3). Now assume that $Y$ is a permutation of
$\\{1,\ldots,n\\}$. Then $H(X,Y)=H(Y)=\log(n)$ and $I(X,Y)=H(X)$, for any
random variable $X$. It then follows from (2) (and the symmetry in $X$ and
$Y$) that

$\Delta I(X,Y)=I(X,Y)-{\rm E}(I(X_{\sigma},Y))=0.$

## Proof of Proposition 3

If $\Delta H(X)=0$, then $d((X,X_{\sigma}))=0$ for all permutation $\sigma$.
In particular, there exists some bijection $f$ such that $X_{\sigma}=f(X)$.
Now assume that for some integer $i$, the event:

$A=\\{\omega:X(\omega)=i\\}$

is such that $1<|A| <n$. Then there exists some $j\neq i$ such that the event
$B=\\{\omega:X(\omega)=j\\}$ is not empty. Choose $a\in A$ and $b\in B$ and
define $\sigma$ as the permutation of $a$ and $b$. Then $X_{\sigma}(a)=X(b)=j$
while $X_{\sigma}(a^{\prime})=X(a^{\prime})=i$ for all $a^{\prime}\in
A\setminus\\{a\\}$. So $X_{\sigma}(a)\neq X_{\sigma}(a^{\prime})$ while
$X(a)=X(a^{\prime})$ for all $a^{\prime}\in A\setminus\\{a\\}$, which
contradicts the existence of some mapping $f$ that $X_{\sigma}=f(X)$. Thus for
each integer $i$, the cardinal of the event $A=\\{\omega:X(\omega)=i\\}$ is
$0,1$ or $n$. This implies that $X$ is constant or equal to some permutation
of $\\{1,\ldots,n\\}$.

## Proof of Theorem 1

Consider two items selected uniformly at random in $\\{1,\ldots,n\\}$. Let
$A_{i_{1}},B_{j_{1}}$ be the clusters of the first item, $A_{i_{2}},B_{j_{2}}$
be the clusters of the second item. In particular, these items belong
respectively to the sets $A_{i_{1}}\cap B_{j_{1}}$ and $A_{i_{2}}\cap
B_{j_{2}}$. The probability of this event is:

$\frac{n_{i_{1}j_{1}}n_{i_{2}j_{2}}}{n^{2}}.$

Now assume that these items exchange their labels for the first clustering, so
that the first item move to set $A_{i_{2}}$ while the second item move to the
set $A_{i_{1}}$. If $i_{1}=i_{2}$ or $j_{1}=j_{2}$, the new contingency matrix
remains unchanged; now if $i_{1}\neq i_{2}$ and $j_{1}\neq j_{2}$, the new
contingency matrix $n^{\prime}_{ij}$ remains unchanged except for the
following entries:

$n^{\prime}_{ij}=\left\\{\begin{array}[]{ll}n_{ij}-1&\text{for
}i,j=i_{1},j_{1}\text{ and }i_{2},j_{2},\\\ n_{ij}+1&\text{for
}i,j=i_{1},j_{2}\text{ and }i_{2},j_{1}.\end{array}\right.$

Using (5), we obtain the similarity between clusterings $A$ and $B$:

$\displaystyle s_{\rm p}(A,B)$ $\displaystyle=\sum_{i_{1}\neq i_{2},j_{1}\neq
j_{2}}\frac{n_{i_{1}j_{1}}n_{i_{2}j_{2}}}{n^{2}}$
$\displaystyle\times\left(\frac{n_{i_{1}j_{1}}}{n}\log\frac{n_{i_{1}j_{1}}}{n}-\frac{n_{i_{1}j_{1}}-1}{n}\log\frac{n_{i_{1}j_{1}}-1}{n}\right.$
$\displaystyle+\frac{n_{i_{2}j_{2}}}{n}\log\frac{n_{i_{2}j_{2}}}{n}-\frac{n_{i_{2}j_{2}}-1}{n}\log\frac{n_{i_{2}j_{2}}-1}{n}$
$\displaystyle+\frac{n_{i_{1}j_{2}}}{n}\log\frac{n_{i_{1}j_{2}}}{n}-\frac{n_{i_{1}j_{2}}+1}{n}\log\frac{n_{i_{1}j_{2}}+1}{n}$
$\displaystyle+\left.\frac{n_{i_{2}j_{1}}}{n}\log\frac{n_{i_{2}j_{1}}}{n}-\frac{n_{i_{2}j_{1}}+1}{n}\log\frac{n_{i_{2}j_{1}}+1}{n}\right),$

where by convention, $x\log x=0$ for any $x\leq 0$. Observing that for any
given $i_{1},j_{1}$,

$\sum_{i_{2}\neq i_{1},j_{2}\neq
j_{1}}n_{i_{1}j_{1}}n_{i_{2}j_{2}}=n_{i_{1}j_{1}}(n-a_{i_{1}}-b_{j_{1}}+n_{i_{1}j_{1}}),$

while for any given $i_{1},j_{2}$,

$\sum_{i_{2}\neq i_{1},j_{1}\neq
j_{2}}n_{i_{1}j_{1}}n_{i_{2}j_{2}}=(a_{i_{1}}-n_{i_{1}j_{2}})(b_{j_{2}}-n_{i_{1}j_{2}}),$

we get by symmetry:

$\displaystyle s_{\rm p}(A,B)$
$\displaystyle=2\sum_{i,j}\frac{n_{ij}(n-a_{i}-b_{j}+n_{ij})}{n^{2}}$
$\displaystyle\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}-1}{n}\log\frac{n_{ij}-1}{n}\right)$
$\displaystyle+2\sum_{i,j}\frac{(a_{i}-n_{ij})(b_{j}-n_{ij})}{n^{2}}$
$\displaystyle\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}+1}{n}\log\frac{n_{ij}+1}{n}\right).$

## Proof of Corollary 1

The proof follows on observing that the second sum in (7) can be written:

$\displaystyle S$ $\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\sum_{i,j}\frac{(a_{i}-n_{ij})(b_{j}-n_{ij})}{n^{2}}$
$\displaystyle\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}+1}{n}\log\frac{n_{ij}+1}{n}\right)$
$\displaystyle=\sum_{i,j:n_{ij}>0}\frac{(a_{i}-n_{ij})(b_{j}-n_{ij})}{n^{2}}$
$\displaystyle\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}+1}{n}\log\frac{n_{ij}+1}{n}\right)$
$\displaystyle-\sum_{i,j:n_{ij}=0}\frac{(a_{i}-n_{ij})(b_{j}-n_{ij})}{n^{2}}\frac{1}{n}\log\frac{1}{n}.$

Since

$\displaystyle\sum_{i,j:n_{ij}=0}$
$\displaystyle{(a_{i}-n_{ij})(b_{j}-n_{ij})}=$
$\displaystyle\sum_{i,j}(a_{i}-n_{ij})(b_{j}-n_{ij})$
$\displaystyle-\sum_{i,j:n_{ij}>0}{(a_{i}-n_{ij})(b_{j}-n_{ij})}$
$\displaystyle=n^{2}-\sum_{i}a_{i}^{2}-\sum_{j}b_{j}^{2}$
$\displaystyle+\sum_{i,j}n_{ij}^{2}-\sum_{i,j:n_{ij}>0}{(a_{i}-n_{ij})(b_{j}-n_{ij})},$

we get:

$\displaystyle S$
$\displaystyle=\sum_{i,j:n_{ij}>0}\frac{(a_{i}-n_{ij})(b_{j}-n_{ij})}{n^{2}}$
$\displaystyle\times\left(\frac{n_{ij}}{n}\log\frac{n_{ij}}{n}-\frac{n_{ij}+1}{n}\log\frac{n_{ij}+1}{n}+\frac{1}{n}\log\frac{1}{n}\right)$
$\displaystyle-\left(n^{2}-\sum_{i}a_{i}^{2}-\sum_{j}b_{j}^{2}+\sum_{i,j}n_{ij}^{2}\right)\frac{1}{n}\log\frac{1}{n}.$

## Proof of Corollary 2

The proof follows from (7) applied to the diagonal contingency matrix
$n_{ij}=a_{i}\delta_{ij}=b_{j}\delta_{ij},$ where $\delta_{ij}$ denotes the
Kronecker symbol.

## References

  * [1] T. M. Cover and J. A. Thomas.  Elements of Information Theory.  Wiley, 1991.
  * [2] M. Gagolewski.  Benchmark suite for clustering algorithms – version 1, 2020.
  * [3] L. Hubert and P. Arabie.  Comparing partitions.  Journal of classification, 2(1):193–218, 1985.
  * [4] M. Meilă.  Comparing clusterings by the variation of information.  In Learning theory and kernel machines, pages 173–187. Springer, 2003.
  * [5] X. V. Nguyen, J. Epps, and J. Bailey.  Information theoretic measures for clusterings comparison: is a correction for chance necessary?  In ICML, 2009.
  * [6] S. Romano, J. Bailey, V. Nguyen, and K. Verspoor.  Standardized mutual information for clustering comparisons: one step further in adjustment for chance.  In International Conference on Machine Learning, pages 1143–1151, 2014.
  * [7] S. Romano, N. X. Vinh, J. Bailey, and K. Verspoor.  Adjusting for chance clustering comparison measures.  The Journal of Machine Learning Research, 17(1):4635–4666, 2016\.
  * [8] A. A. Taha and A. Hanbury.  Metrics for evaluating 3d medical image segmentation: analysis, selection, and tool.  BMC medical imaging, 15(1):29, 2015.
  * [9] B. Thirion, G. Varoquaux, E. Dohmatob, and J.-B. Poline.  Which fmri clustering gives good brain parcellations?  Frontiers in neuroscience, 8:167, 2014.
  * [10] N. X. Vinh, J. Epps, and J. Bailey.  Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance.  The Journal of Machine Learning Research, 11:2837–2854, 2010.
  * [11] B. Wang, J. Zhu, E. Pierson, D. Ramazzotti, and S. Batzoglou.  Visualization and analysis of single-cell rna-seq data by kernel-based similarity learning.  Nature methods, 14(4):414–416, 2017.
  * [12] Z. Yang, R. Algesheimer, and C. J. Tessone.  A comparative analysis of community detection algorithms on artificial networks.  Scientific reports, 6:30750, 2016.
  * [13] J. Zhang, P. Kapli, P. Pavlidis, and A. Stamatakis.  A general species delimitation method with applications to phylogenetic placements.  Bioinformatics, 29(22):2869–2876, 2013.