# Forming off-center massive black hole binaries in dwarf galaxies through Jacobi capture

Thibaut L. François<fnref n="1" />, Christian M. Boily<fnref n="1" />, Jonathan Freundlich<fnref n="1" />, Simon Rozier<fnref n="2" />, and Karina Voggel<fnref n="1" />

Received 13 November 2023 / Accepted 26 April 2024

## ABSTRACT

It is well established that black holes reside in the central regions of virtually all types of known galaxies. Recent observational and numerical studies however challenge this picture, suggesting that intermediate-mass black holes in dwarf galaxies may be found on orbits far from the center. In particular, constant-density cores minimize orbital energy losses due to dynamical friction, and allow black holes to settle on stable off-center orbits. Using controlled simulations, we study the dynamics of off-center black holes in dwarf galaxies with such cores. We propose a new scenario to describe off-center mergers of massive black holes, starting with a Jacobi capture. We focus on initially circular and co-planar black hole orbits and explore a large parameter space of black hole masses and orbital parameters. We find that Jacobi captures are a complex and chaotic phenomenon that occurs in about 13% of cases in this simplified setup, and we quantify how the likelihood of capture depends on the simulation parameters. We note that this percentage is likely an upper limit of the general case. Nevertheless, we show that Jacobi captures in cored dwarf galaxies can facilitate the formation of off-center black hole binaries, and that this process is sufficiently common to have a substantial effect on black hole growth. While our setup only allows for temporary captures, we expect dissipative forces from baryons and post-Newtonian corrections to maintain the captures over time and to lead to the formation of stable binary systems. This motivates future studies of the effectiveness of such dissipative forces, within stripped nuclei or globular clusters, in forming stable bound systems.

**Key words.** chaos – gravitation – methods: numerical – galaxies: dwarf – galaxies: kinematics and dynamics

## 1. Introduction

Over the last decades, observations have uncovered massive black holes (BHs) at the center of massive galaxies (see Kormendy & Richstone 1995). Their mass correlates with the properties of the host galaxy (Ferrarese & Merritt 2000; Gebhardt et al. 2000; van den Bosch 2016), which suggests a common past history. Recent studies with different probes (e.g., optical, X-ray, Radio or infrared) have brought to light the presence of such BHs within dwarf galaxies (see Reines et al. 2013; Pardo et al. 2016; Mezcua et al. 2018, 2019; Zaw et al. 2020; Cann et al. 2020). And although a clear consensus has yet to emerge due to limitations in observed kinematics and dynamical modeling, studies reveal indications of intermediate-mass black holes (IMBHs) at the centers of some globular clusters (see Noyola et al. 2008; van der Marel & Anderson 2010; McNamara et al. 2012; Lützgendorf et al. 2013a,b; Lanzoni et al. 2013). Major galaxy mergers, accretion events and in-spiraling of globular clusters can lead to multiple of these BHs within the same galaxy (see Schweizer et al. 2018; Nguyen et al. 2019; Kollatschny et al. 2020; Shen et al. 2021; Stemo et al. 2021; Voggel et al. 2022; Pechetti et al. 2022).

These BHs can then merge through a three-stage process (Begelman et al. 1980): first, massive BHs are driven to the center of the galaxy, due to the loss of angular momentum from dynamical friction (Chandrasekhar 1943; Binney & Tremaine 2008, hereafter BT08, Sect. 8.1, Eqs. (8.1)–(8.4)). When the separation between two BHs is comparable to the gravitational influence radius of the larger one, they form a bound pair and

a second phase takes place. The binary gives away some of its angular momentum via three-body interactions to stars with lower angular momentum (see Heggie 1975). In scenarios where a sufficient number of stars are available to facilitate angular momentum exchange, the binary system enter the relativistic regime, emitting gravitational waves before merging. These three distinct phases constitute the conventional framework for the merger of massive BHs through the assembly of galaxies (depicted in the top panel of Fig. 1).

The merging scenario relies mainly on dynamical friction to bring BHs together at the center. However, the time required to bring IMBH ($10^3-10^6 M_\odot$) to the center of a typical dwarf galaxy can be much longer than a Hubble time in the case of cored dark matter density profiles and/or low stellar densities<fnref n="1" /> (see Bellovary et al. 2019, 2021; Pfister et al. 2019). As a consequence, BHs are left on off-center orbits, with their guiding radius reducing only very slightly over 14 Gyr. According to simulations and observations, it is estimated that half of the IMBHs are off-center in dwarf galaxies (see Bellovary et al. 2019, 2021; Reines et al. 2020). Furthermore, their electromagnetic detection is highly challenging due to their low accretion luminosities, and their dynamical detection is equally problematic as they have a very small impact on the dynamics of the stars within their sphere of influence.

<footnote n="1">¹ Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France
e-mail: thibaut.francois@proton.me</footnote>

<footnote n="2">² Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="1">¹ With Eq. (8.1) of BT08, we compute a time ~22 Gyr for a $10^6 M_\odot$ BH to fall from 5 kpc to the center of a cored dwarf spheroidal galaxy with a virial mass of $10^{11} M_\odot$. For a $10^3 M_\odot$ BH falling from 10 kpc, it takes ~85 000 Gyr.</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Central versus off-center BH merger scenarios. Top row: classical massive BH merger scenario in the galactic center. Phase 1: the massive BHs lose angular momentum through dynamical friction with the galaxy and halo, which causes them to sink toward the center. Phase 2: once both BHs have reached the center, they may form a bound pair and continue to lose energy via dissipative forces. Phase 3: the binary eventually enters the relativistic regime and loses energy through gravitational waves emission, allowing the two BHs to merge. Bottom row: proposed scenario for an off-center merger. Phase 1: the BHs are on off-center orbits e.g., as a result of long dynamical time for their radius to decrease through dynamical friction or core stalling acting against dynamical friction. If both BHs are sufficiently close energetically, they may form an off-center bound pair through Jacobi capture. This is the scenario explored in the present paper.</figcaption></figure>

In certain instances, however, the long infall time argument alone does not provide a satisfactory explanation for the existence of off-center objects, such as in the case of the Fornax dwarf galaxy. This galaxy harbors five globular clusters orbiting at approximately 1 kpc from its center (Mackey & Gilmore 2003), and based on conventional assumptions, these clusters should have migrated from their current locations within a few billion years (see Goerdt et al. 2006, though this picture has been challenged, given the modeling freedom of the initial GC formation radius, and uncertainties arising from projection effects, see e.g., Borukhovetskaya et al. 2022). This is, in part, due to dynamical processes which cannot be modeled using Chandrasekhar dynamical friction alone: numerical simulations have revealed instances of the perturber's fall coming to a halt, known as "core stalling" (see Read et al. 2006; Inoue 2011; Petts et al. 2015), as well as the occurrence of an outward push, or "dynamical buoyancy" (see Cole et al. 2012). These phenomena are observed in constant-density cores (as those observed at the center of dark matter halos in low-mass galaxies, see Binney & Evans 2001; Borriello & Salucci 2001; de Blok et al. 2001), and are not accounted for in Chandrasekhar's description. Recent analytical models provide an explanation in terms of resonances (see Boily et al. 2008; Kaur & Sridhar 2018; Banik & van den Bosch 2021): when a perturber enters the core region, significant resonances responsible for dynamical friction are suppressed, leading to a cessation in the fall. Within a core, these resonances can even give energy to the perturber, causing it to migrate outward.

In light of these cases where IMBHs are found off-center, either because their infall time greatly exceeds the dynamical time of the system or because it is halted due to the presence

of a constant density core, the question arises whether a IMBH merger is possible outside the galactic center. Indeed, the aforementioned three-phase scenario must be revised. Firstly, there is no guarantee that the BHs may get sufficiently close in phase space to form a binary, as opposed to the classical case where the galactic center acts as an attractor, allowing an agglomeration at a single point. In this context, specific phase-space conditions must be met for them to form bound configurations away from the center.

Our proposed scenario of off-center Jacobi capture may be modeled, as a first approximation, as a circular-restricted three-body problem. This is a variation of the three-body problem in which two very light objects orbit around a much more massive body. This problem, known as Hill's problem, has been extensively studied to describe the movements of satellites around a planet (Petit & Hénon 1986, 1990; Henon & Petit 1986) as well as gravitational wave capture of stellar-mass BHs around a supermassive BH (Boekholt et al. 2023; Li et al. 2022). The principal distinctions in our case lies in the fact that a galaxy is an extended massive body rather than point-like, so the mass of the central body depends on the orbital radius of the BHs and the shape of the potential is different, featuring a density core at the center. For a merger to occur, it is first necessary for the two light bodies to meet a suitable configuration where their gravitational attraction becomes comparable to the force generated by the massive body.

In this case, they bond for a certain period of time, known as Jacobi capture, and produce several close encounters. This is the first assumed step in an off-center merger scenario. Then, it has been observed that dissipative forces play a crucial role in the hardening of the binary (e.g., Fabian et al. 1975;
Goldreich et al. 2002). In Li et al. (2022), Boekholt et al. (2023) this energy dissipation is facilitated by the emission of gravitational waves, while in Tagawa et al. (2018, 2020), Li et al. (2022), Rowan et al. (2023), DeLaurentiis et al. (2023), it is dynamical friction with the surrounding gas that tightens the binary. We anticipate similar behaviors for the case of two IMBHs: the bound stellar populations that typically surround them can generate dynamical friction on the other BH, eventually leading to the last step of the classical scenario (see bottom panel on Fig. 1).

In this paper, we focus on the first step of this process: the Jacobi capture. We aim to analyze which initial configurations, in terms of BH masses and kinematics, can allow Jacobi capture to occur within a Hubble time. In order to isolate the capture step, we perform idealized simulations of two point-like BHs in a smooth background gravitational potential, varying the initial BH masses and kinematics. This work allows to pinpoint regions in phase space where Jacobi captures are most likely to occur, paving the way for a comprehensive investigation of the subsequent merger stages through N-body simulations.

This paper is organized as follows. In Sect. 2, we outline our methodology, including the gravitational potential, parameter selection, core stalling radius and influence radius computation, capture criterion and an examination of the chaotic behavior of Jacobi captures. In Sect. 3, we present our findings. Initially, we examine the general case involving five parameters, wherein no assumptions are made regarding the cause of the offset of BHs. Subsequently, we explore the scenario where BHs are positioned at their core stalling radius. Lastly, we determine the probability of capture by employing a Monte Carlo generation of merger trees to select BH masses. Our outcomes are summarized in Sect. 4 and we discuss our results.

## 2. Method

### 2.1. Numerical setup

In this study, we aim to explore the influence of different initial properties of BHs on the occurrence of a Jacobi capture. We also aim to provide physical insights on the capture process itself, occurring off-center in an extended galactic potential. As a consequence, we are reduced to a number of simplifying assumptions. Because we want to sample a large parameter space of initial BH masses and kinematics, we cannot afford to simulate the host galaxy as a live N-body component. Indeed, the problem we are considering requires high computational accuracy along two lines. First, we need to resolve the dynamics of the BH binary, implying a high time resolution of the integration. Second, an extreme number of galactic particles is required to resolve processes such as dynamical friction<fnref n="2" />. As a workaround, we consider the galaxy as a background static potential, in a simplified setup that enables us to execute thousands of simulations within a reasonable timeframe.

We integrate the orbits of two BHs in an external gravitational potential representing the dark matter halo of a dwarf galaxy. Simulations in this paper made use of the REBOUND N-body code (Rein & Liu 2012). The simulations were integrated using IAS15, a 15th order Gauss-Radau integrator

<figure><img src="image_2.png" /><figcaption>Fig. 2. Density profile and integrated mass profile of the dark matter halo (background potential). The two horizontal lines represent the limits of the mass range we consider for the BHs. The shaded area denotes the scope of this study, and the x-axis at the top of the figure is expressed in units of the core radius of the potential $r_c$, defined by its inflection point.</figcaption></figure>

(Rein & Spiegel 2015). This integrator incorporates a step-size control mechanism to automatically select an optimal timestep and is adept at handling close encounters and orbits with high eccentricity. With systematic errors maintained well below machine precision, IAS15 demonstrates adherence to Brouwer’s law, where energy errors exhibit random walk behavior.

We choose a dwarf galaxy halo as various studies indicate that this type of galaxy hosts a substantial population of off-center BHs (Bellovary et al. 2019, 2021). Furthermore, in the standard $\Lambda$CDM paradigm dwarf galaxies are dominated by dark matter, hence we can initially approximate them by considering only their halo. We consider a galaxy with a $10^{11} M_\odot$ cored dark matter halo with a virial radius $r_{200}$ = 100 kpc (this corresponds to a stellar mass of $\sim 2 \times 10^8 M_\odot$). This halo mass range is expected to develop cored profiles more efficiently (see Tollet et al. 2016; Freundlich et al. 2020). The dark matter component of our dwarf galaxy model is described by a cored $(\alpha, \beta, \gamma)$ density profile (Hernquist 1990; Zhao 1996)

$$ \rho = \rho_0 \left(\frac{r}{a}\right)^{-\gamma} \left[1 + \left(\frac{r}{a}\right)^{\alpha}\right]^{\frac{\gamma-\beta}{\alpha}}, \quad (1) $$

with $\alpha = 1$, $\beta = 3$ and $\gamma = 0$ respectively the outer, transitional and inner slope of the density profile. The scale radius is $a = 6.7$ kpc. The density profile as well as the integrated mass profile are depicted in Fig. 2.

### 2.2. Toy model

We adopted a simplified model where the two BHs are on circular orbits in the same plane. This configuration closely resembles Hill’s problem (Petit & Hénon 1986, 1990; Henon & Petit 1986; Boekholt et al. 2023). However, our massive body is extended rather than point-like, so its mass depends on the orbital radius of the BHs, and the potential shape differs, featuring a density core at the center. If our BHs are close to the center (which is the case for core stalling), the mass of the central body is comparable to the mass of the BHs. When the proximity of the two

<footnote n="2">² We simulate IMBHs of $10^3-10^6 M_\odot$ in a DM halo of $10^{11} M_\odot$. In order to resolve dynamical friction, a mass gap of $\sim 1/1000$ is required between the galactic particles and the BHs. Consequently, the mass of the galactic particles must not exceed $1 M_\odot$. We therefore need to simulate $10^{11}$ particles which is well outside the technical capabilities at our disposal in terms of memory and computing time.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Parameters of an initial configuration, $r_i$ the initial radius of the inner BH, $\Delta r_i$ the radial separation between BHs, $\Delta\varphi_i$ the phase difference and $M_1, M_2$ the two masses.</figcaption></figure>

<table><caption>**Table 1.** Parameter ranges for the five-parameter case.</caption><thead><tr><th></th><th>Minimum value</th><th>Maximum value</th><th>Number of values</th><th>Spacing</th></tr></thead><tbody><tr><td><i>M</i><sub>1</sub>/<i>M</i><sub>☉</sub></td><td>10<sup>3</sup></td><td>10<sup>6</sup></td><td>8</td><td>Log</td></tr><tr><td><i>M</i><sub>2</sub>/<i>M</i><sub>☉</sub></td><td>10<sup>3</sup></td><td>10<sup>6</sup></td><td>8</td><td>Log</td></tr><tr><td><i>r</i><sub>i</sub>/kpc</td><td>0.1</td><td>6</td><td>14</td><td>Log</td></tr><tr><td>Δ<i>r</i><sub>i</sub>/kpc</td><td>0.08</td><td>1.2</td><td>30</td><td>Log</td></tr><tr><td>Δ<i>φ</i><sub>i</sub></td><td>0</td><td>2π</td><td>9</td><td>Lin</td></tr></tbody></table>

**Notes.** For a definition of the parameters in use, see Fig. 3.

BHs leads to gravitational forces comparable to those exerted by the galaxy, temporary captures occur. These are known as Jacobi captures and are marked by multiple close encounters between BHs.

Our aim is to probe the phase space of the initial conditions to determine the percentage of cases in which the BHs produce a Jacobi capture. A priori, we have 14 parameters to describe an initial configuration: 6 phase space coordinates per BH ($r, v_r, \theta, v_\theta, \varphi, v_\varphi$) and the mass of each object ($M_1$ and $M_2$). We decide to simplify the problem by considering the 2D case of two BHs orbiting in the same plane, which reduces our number of parameters to 10. We place the BHs on circular orbits in the $\theta = \pi/2$ plane which fixes $v_r$ and $v_\varphi$ once the position is known. Finally, we do not need each BH's phase ($\varphi_1$ and $\varphi_2$) but only the phase difference $\Delta\varphi$ because of cylindrical symmetry. We thus end up with five initial parameters: $r_1, r_2, \Delta\varphi_i, M_1$ and $M_2$. However, it is more convenient to use $r_i, \Delta r_i, \Delta\varphi_i, M_1$ and $M_2$ with $r_i$ the initial radius of the inner BH and $\Delta r_i$ the initial radial separation between the two BHs (these parameters are shown in Fig. 3). Without loss of generality we impose that $M_1$ is always the inner BH (this means that $\Delta r_i$ is always positive). We quantify the effect of the initial parameters on Jacobi capture and the percentage of cases in which such an event occurs. Parameter ranges are summarized in Table 1.

## 2.3. Core stalling

The general five-parameter case helps us emphasize the influence of masses and kinematics on Jacobi captures. Yet, to derive a probability of capture, we must delve deeper into a physically

<table><caption>**Table 2.** Parameter ranges for the three-parameter case.</caption><thead><tr><th></th><th>Minimum value</th><th>Maximum value</th><th>Number of values</th><th>Spacing</th></tr></thead><tbody><tr><td><i>M</i><sub>1</sub>/<i>M</i><sub>☉</sub></td><td>10<sup>3</sup></td><td>10<sup>6</sup></td><td>20</td><td>Log</td></tr><tr><td><i>M</i><sub>2</sub>/<i>M</i><sub>☉</sub></td><td>10<sup>3</sup></td><td>10<sup>6</sup></td><td>20</td><td>Log</td></tr><tr><td>Δ<i>φ</i><sub>i</sub></td><td>2π/1000</td><td>2π</td><td>400</td><td>Lin</td></tr></tbody></table>

motivated parameter space with a reduced number of degrees of freedoms. To achieve this, we conduct a second simulation, assuming that BHs are off-center at their core stalling radius. This fixes the radii once the masses are defined, effectively reducing the number of free parameters to three. The stalling radius (noted $r_{cs}$) is directly related to the mass of the BH through (see Kaur & Sridhar 2018)

$$ \Omega_p(r_{cs}) = \Omega_b \quad \text{with} \quad \Omega_p(r) = \sqrt{G \frac{M_b(<r) + M_p}{r^3}}, \qquad (2) $$

where $\Omega_p$ is the BH's orbital frequency, $\Omega_b$ the orbital frequency of stars at the very center of the cored potential and $M_b(<r)$ the enclosed mass of our background density. Numerical integration of Eq. (2) leads to a near-power solution which we find is accurately recovered from

$$ \frac{r_{cs}}{\text{kpc}} \approx A \left( \frac{M}{M_{\odot}} \right)^{1/4} \qquad (3) $$

with $A \approx 1/80$ a dimension-less constant. We note that $r_{cs}$ increases with the mass of the perturber, which is in agreement with Read et al. (2006). We end up with only three parameters: $M_1, M_2$ and $\Delta\varphi_i$ (parameter ranges used are summarized in Table 2).

We can impose $M_2 \ge M_1$ without loss of generality, ensuring that $M_1$ consistently represents the inner BH, as in the previous section. Not imposing this condition would mean simulating all possible cases twice, merely swapping the BH labels. This approach would yield mass-symmetric outcomes, offering limited insights into the impact of BH positioning on captures. Thus, we consider a total of 84 000 initial configurations.

## 2.4. Black hole radius of influence

The volume of influence, also known as the Hill volume, defines a spatial domain surrounding the BH where its gravitational influence dominates. Traditionally, this volume is approximated as a sphere, characterized by the Hill radius (or radius of influence). In the classical restricted three-body problem, this radius is derived by determining the distance between the body of interest and the Lagrange point L1, situated along the line connecting the two massive bodies, between them. Under the distant-tide approximation, assuming a Hill radius much smaller than the distance from the body to the center and a significant mass hierarchy, an analytical expression for the Hill radius emerges: $R_H = (m/3M)^{1/3}R$, where $m$ represents the mass of the body of interest, $R$ its distance from the center, and $M$ the mass of the central body (BT08, Sect. 8.3.1, Eq. (8.91)).

To address the extended galactic case, using the previous analytical expression for the Hill radius, which involves replacing $M$ by $M_b(<R)$ the enclosed mass of our background density at the radius of the BH, is not desirable. Complications arise due to several reasons: (i) our investigation deviates from previous assumptions, particularly when we examine the dynamics
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Effective potential of the galaxy + BH system plotted in the reference frame rotating at the angular frequency of the BH orbit in the (x, y) plane. The coordinates (x, y) = (0, 0) denote the center of mass, calculated for the sum of the BH mass and the mass of the host galaxy integrated up to the BH orbit. Hence the galactic center is shifted significantly from (0, 0) in some cases. The Lagrange points, along with the positions of the black hole (denoted BH) and galactic center (denoted GC), are highlighted in red. R = f × r<sub>c</sub> is the orbital radius of the BH, where f is a dimensionless factor and r<sub>c</sub> represents the core radius of the potential, defined by its inflection point. We observe distinct differences in the topology of the effective potential, as well as in the positioning and number of Lagrange points, between the cored galactic potential (on the left) and the Keplerian case (on the right) where we approximate the galaxy as a point mass equal to the enclosed galactic mass up to the BH orbit.</figcaption></figure>

of BHs in close proximity to the center ($R_H \sim R$), where the galaxy's integrated mass is comparable in order to that of the BH ($M_b(<R) \sim M_{BH}$); (ii) in density profiles featuring a central core, distinct Lagrange point structures emerge, differing from those observed in the Keplerian case (Banik & van den Bosch 2022); (iii) the density profile employed lacks an analytical for-

mula for the potential. For these reasons, we have conducted numerical computations to determine the Lagrange points and the radius of influence.

Figure 4 shows the effective potential of the galaxy + BH system, plotted in the reference frame rotating with the BH. Each row in Fig. 4 corresponds to a distinct position of the BH
within the core. The left-hand column is computed for the cored galactic potential, while the right-hand column corresponds to the Keplerian case, where we approximate the galaxy as a point mass equal to the enclosed galactic mass up to the BH orbit. Near the galactic center, the topology of the effective potential in the $(\alpha, \beta, \gamma)$ case differs from the Keplerian case. As the BH approaches the center and reaches a critical radius $R_{crit}$, a bifurcation occurs and the Lagrange points L1, L3 and L0 merge into a single point. In the last two rows of Fig. 4, we observe the scenario where the BH's radius exceeds this critical value, resulting in the same number of Lagrange points as in the Keplerian case, albeit at different positions. Conversely, the first row illustrates the case where the BH's radius is smaller than the critical value, leading to the merging of Lagrange points L1, L3, and L0³. As the BH's orbital radius increases, the galactic potential gradually approaches the Keplerian model.

Figure 5 shows the Hill radius as a function of the BH orbital radius. It is computed in two ways, by taking the distance between the BH and the L1 Lagrange point and the distance between the BH and L2. Four background potentials were considered: a Keplerian potential and three $(\alpha, \beta, \gamma)$ profiles with different inner density slopes.

First, we observe that the distance between the BH and L1 matches that between the BH and L2 around and beyond the core radius. However, this symmetry does not persist within the core. Here, the distant-tide approximation, which enforces precise L1/L2 symmetry by eliminating asymmetry terms in the calculations, breaks down. The exception is the NFW profile ($\gamma = 1$), where symmetry holds across all radii.

In our study profile with a core ($\gamma = 0$), we discern a discontinuity in the BH-L1 distance at a specific radius, indicative of the earlier mentioned bifurcation. Beyond this critical radius, the L1 Lagrange point lies between the BH and the galactic center. However, as the BH's orbital radius diminishes, the L1 point converges toward the galactic center. Below the critical radius, the distance between the BH and the galactic center is measured, accounting for the linearity of the curve.

The Keplerian profile closely aligns with the Hill radius predicted by the distant-tide approximation formula, maintaining this consistency even within the core despite L1/L2 asymmetry. It is only in very close proximity to the center that significant divergence occurs. In contrast, the extended profiles $(\alpha, \beta, \gamma)$ exhibit clear discrepancies with the distant-tide approximation, particularly within the core, notably for inner slopes $\gamma = 0$ and $\gamma = 1$. Within our study profile ($\gamma = 0$), we see a difference of a factor of 2 between the distant-tide approximation's analytical formula and the precise computation. This discrepancy, coupled with the distinct shape of the Hill radius curve compared to the Keplerian case, suggests that the galaxy cannot be adequately treated as a radius-dependent point mass. Moreover, the disparities among the three extended density profiles highlight the significant dependence of the Hill radius in our study area on the inner density slope. Consequently, we anticipate an equally pronounced sensitivity of the capture process to this parameter.

Figures 4 and 5 show that the dynamics and likelihood of Jacobi captures are significantly influenced by the galactic context. Firstly, the zone of influence of BHs is different from the Keplerian case, and secondly, the altered topology of the effective potential gives rise to novel types of BH orbits, including

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

Fig. 5. Hill radius ($R_H$) as a function of the black hole (BH) orbital radius $R$ in units of the core radius for a mass $M_{BH} = 10^6 M_\odot$. The Hill radius is computed as the BH-L1 distance (dashed lines) and the BH-L2 length (solid lines). Four background potentials were considered: a Keplerian potential (shown as light-blue lines); and three $\alpha, \beta, \gamma$ profiles with $\alpha = 1, \beta = 3$ and inner slopes $\gamma = 0$ (red), $1/2$ (green) and $1$ (yellow) respectively. The thick dark blue solid line is for the distant-tide approximation, which goes to a constant value as $R \to 0$. On the figure, the shaded area denotes the scope of this study. The core stalling radius and galactic core radius are indicated by vertical lines. The separation of the dashed and solid lines indicates an asymmetry between L1 and L2 on either side of the black hole when entering the core region. Note that for the case of an NFW profile ($\gamma = 1$, yellow) the L1 and L2 solutions very nearly overlap inside the core. The discontinuity of the dashed red line shows a bifurcation of the Lagrange points in the core profile. It is striking that the Keplerian potential remains close to the distant-tide approximation well inside the core of the galaxy, before diverging significantly at $\log R/r_c \le -1.3$. By comparison, extended profiles may show differences with this approximation exceeding a factor 2 anywhere inside the core for $\gamma = 0$ (red).

Pacman orbits (see Fig. C.1), which are pivotal for core stalling (Banik & van den Bosch 2022).

## 2.5. Capture criterion

We aim to contextualize Jacobi captures within an observable framework to evaluate the likelihood of forming and detecting off-center BH binaries. To achieve this, we conduct simulations spanning a timeframe of 14 Gyr. Consequently, some captures may be truncated by halting the simulation, while others that are theoretically feasible may not have sufficient time to occur. For instance, when the orbital radii of BHs are extremely close, the time required to catch-up their phase difference may exceed 14 Gyr. Likewise, in order to deduce a probability from initial conditions, we incorporate conditions that lead to capture at $t = 0$ Gyr. These cases are indeed conditions allowing BH binding off-center, and due to the ergodic nature of this study, any stretch of 14 Gyr should be equivalent to another. Put differently, a time shift in the simulation range would interchange initial and dynamical captures (for instance, initiating the simulation later would render some dynamical captures as initial and vice versa). Lastly, although we simulate 14 Gyr of evolution during which an initial condition can result in several captures, we solely focus on the first, as higher-order effects can significantly influence the sequence of events in reality (such as dynamical friction and post-Newtonian corrections).

<footnote n="3">³ A detailed investigation is warranted to confirm whether this convergence indeed represents a single point or falls below our resolution limit. This aspect will be addressed in future research endeavors.</footnote>
To determine a capture, two conditions must be met. First, the binding energy between the BHs must be negative. This is given by

$$E_b = \frac{1}{2}\mu v_{\text{rel}}^2 - \frac{GM_1M_2}{d_{\text{rel}}}, \quad (4)$$

with $\mu = M_1M_2/(M_1+M_2)$ the reduced mass, G the gravitational constant, $d_{\text{rel}}$ the relative distance between the two BHs and $v_{\text{rel}}$ their relative velocity. Secondly, during the period of negative binding energy, the distance between the BHs must fall below the binary Hill radius. This radius is defined as the average distance between BH-L1 (or BH-L0 if L1 no longer exists) and BH-L2 for a total mass $M_{\text{tot}} = M_1+M_2$ located at $\tilde{r} = (M_1\cdot r_1+M_2\cdot r_2)/(M_1+M_2)$. This additional condition prevents mislabeling due to coincidental alignment of velocity vectors, which could lead to a negative energy reading. It is important to bear in mind that classical two-body dynamics do not allow for gravitational capture, as BHs get accelerated upon approach due to energy conservation. However, in the three-body system, the presence of a third body enables the binding energy to dip below zero and later return to positive values, allowing for unbinding (Petit & Hénon 1986; Renaud et al. 2011; Peñarrubia 2023). The time spent in this sub-zero state, during which the BHs can revolve around each other a few times, is challenging to predict due to the chaotic nature of the three-body problem, occurring over variable timescales. However, dissipative forces like dynamical friction, three-body scattering or post-Newtonian corrections can enhance temporary captures when considered (Samsing et al. 2018; Li et al. 2022; Boekholt et al. 2023; Rowan et al. 2023; DeLaurentiis et al. 2023). Consequently, since our focus in this study is on Jacobi captures, our aim is to identify sub-zero crossings independently of their duration.

Although flybys with a single close encounter are traditionally distinguished from captures with more than one close encounter, we opt to label both scenarios as captures in this article: (i) unlike other distance-based criteria (Boekholt et al. 2023), Eq. (4) is more stringent as it incorporates a velocity condition for BHs. This refinement ensures that the surviving flybys consist of BHs close in energy and with mean eccentricities below 1 (see Fig. A.1). Consequently, we anticipate that these flybys may be more easily stabilized in the future (for a comparison with another, less restrictive criterion yielding numerous hyperbolic flybys, refer to Sect. B); (ii) it is conceivable that even a single close encounter could lead to eventual stabilization, especially when dissipative forces, such as those arising from the post-Newtonian corrections, play a significant role at the pericenter of eccentric binaries (cf. Fig. A.1). The question of the conditions necessary for stabilization remains a subject for future investigation.

Figure 6 illustrates six instances of Jacobi captures characterized by varying numbers of close encounters. The top row displays the relative coordinates between the BHs during these captures, with the reference frame centered on the most massive BH and rotating with it. The binary Hill radius is depicted by the dotted circle. The bottom row shows the binding energy and the relative distance between the BHs.

## 2.6. Chaos

To underscore the chaotic behavior of Jacobi captures in this context, we choose a specific capture as a point of reference and examine how small changes in initial conditions affect the number of close encounters between the BHs. In the left-hand panel of Fig. 7, this reference capture is shown, characterized

by 5 close encounters ($N_c$) between the BHs. We set the initial conditions ($M_1, M_2, r_i, \Delta r_i, \Delta \varphi_i$) for this capture and conducted two refined simulations: one varying ($M_1, M_2$) more finely and another varying the radial parameters ($r_i, \Delta r_i$). For both sets of parameters, we sampled 400 linearly distributed values, varying them by plus or minus 1% relative to the original initial conditions. The central and right-hand panels of Fig. 7 present the number of close encounters between the BHs for each case. In these regions of phase space, the number of close encounters fluctuates between 5 to 33 for the central panel and 2 to 34 for the right-hand panel. The upper limits are due to the fixed integration time, truncating the captures beyond 14 Gyr. The significant irregularities in the count of close encounters when both pairs of parameters are varied illustrate the high sensitivity to initial conditions. On the central panel, the patterns of irregularity align along lines where $M_1 + M_2 \sim C$, with C being a constant. This occurs because at a larger radius (here, the internal BH is 6 kpc from the galactic center), dynamics is governed by the binary Hill radius, which depends on the total mass.

# 3. Results

## 3.1. Bound black holes in the general case

For the general case where we made no assumption regarding the cause of the offset from the center, we integrated $8 \times 8 \times 14 \times 30 \times 9 = 241920$ initial configurations (cf. Table 1) over 14 Gyr. We find a capture in 7.1% of cases, of which 26.2% are flybys. These captures can be further categorized into two distinct types: those formed at the beginning of the simulations ($t = 0$ Gyr, hereafter referred to as initial captures), comprising 66.3% of the total captures, and those emerging during the integration at later times ($t > 0$ Gyr, hereafter referred to as dynamical captures), making up the remaining 33.7% of the captures. It is essential to bear in mind that in this section our focus is on the initial conditions that initiate the capture process i.e., the transition from two separate two-body problems to a three-body problem. Consequently, the results presented here are not subject to chaos. We examine the initiation of the process rather than the subsequent three-body interaction (such as capture duration or number of close encounters) which typically exhibit chaotic behavior (cf. Sect. 2.6).

Figure 8 shows the percentage of initial captures in the space of different parameter pairs, marginalized over the others. The right-hand panels show this percentage in ($M_1, \Delta r_i$) space (top) and in ($M_2, \Delta r_i$) space (bottom). In both panels, the proportion of captures increases as the radial separation between the BHs decreases and their masses increase. These conditions elevate the potential term of the binding energy (see Eq. (4)), consequently reducing the binding energy of the binary. Beyond a certain threshold of radial separation, captures cease to occur. This limit does not precisely match the maximum binary Hill radius (for the parameter range explored in this study, $\log R_H^{\text{max}} = -0.35$), as it is possible that initially, the BHs are slightly beyond the sphere of influence radius but their relative velocity may be low enough to initiate the capture process. Nevertheless, the maximum radial separation allowing initial capture is rather close to $R_H^{\text{max}}$.

The top left panel of Fig. 8 shows the percentage of initial captures as a function of ($M_1, r_i$), while the lower left panel shows it as a function of ($M_2, r_i$). Regarding the $r_i$ dependence, we observe a sharp increase in the percentage of capture at low radii ($\log r_i/kpc < -0.4$) for higher masses, which then stabilizes for higher radii. The constant part at higher radii is attributed to the nearly constant Hill radius within this range
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Examples of Jacobi captures with different numbers of close encounters between BHs. Top row: relative coordinates between BHs during Jacobi captures. The reference frame is centered on the most massive BH and rotates with it. The binary Hill radius is marked by the dotted circle. Bottom row: the black line represents the binding energy ($E_b$) per unit mass of the binary ($M_{\text{tot}} = M_1 + M_2$), with negative values shown in the colored area. The black marker indicates the point when the energy becomes negative. The blue curve represents the relative distance ($d_{\text{rel}}$) between the BHs. The binary Hill radius is indicated by a dashed blue line. The number of close encounters, $N_e$, is given in the box. To contextualize captures within observational framework, we extend our simulations over a period of 14 Gyr. Consequently, we observe instances where captures are halted prematurely, as illustrated by cases with $N_e = 20$ and $N_e = 77$. Additionally, we retain flybys due to our energy-based criterion, which is more stringent compared to distance-based criteria. This ensures that surviving flybys comprise black holes with similar energy levels, potentially facilitating stabilization through dissipative forces.</figcaption></figure>

(refer to Fig. 5 for the case of a $10^6 M_{\odot}$ BH), the potential term is therefore dominant in this zone (for further details, see Fig. B.1 and Appendix B). At low radii, the probability increase cannot be explained by Hill's radius, which decreases in this zone. Instead, the kinetic component of the binding energy plays a significant role in these captures. In this low-radius zone, the gradient of circular velocities decreases due to the presence of the density core ($\gamma = 0$), leading to a reduction in the relative velocity between BHs and consequently lowering the binding energy. Another contributing factor is the marginalization over other parameters. At low radii, there are more values of $\Delta\varphi_i$ that enable the BHs to remain within their sphere of influence, despite its decrease.

Figure 9 shows the percentage of dynamical captures (i.e., happening at $t > 0$ Gyr) as a function of ($M_1, M_2$). The percentage of captures increases as the mass of the two BHs increases,

and we notice that the percentage is higher when the mass of the outer BH is greater than the mass of the inner one ($M_2 > M_1$). This is of particular interest within the context of core stalling, which predicts a larger stalling radius for more massive BHs (cf. Eq. (3)). This outcome can be attributed to the interaction between BHs, which imparts eccentricity to the initially circular orbits over time. This eccentricity gain depends on the BH mass: when one BH is much more massive, it remains nearly circular while the lighter one becomes eccentric. Consequently, if the more massive BH is positioned closer to the center, the lighter one will approach it more centrally, resulting in higher kinetic energy and a more challenging capture. On the contrary, when the more massive BH occupies the outer position, the lighter BH approaches it at a larger radius, resulting in reduced kinetic energy and thus facilitating the capture process. In Fig. 9, we also overplotted contour lines representing the average time it
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Influence of the initial conditions on the number of close encounters ($N_e$). Left panel: reference Jacobi capture ($M_1^0 = 3.7276 \times 10^5 M_\odot$, $M_2^0 = 3.7276 \times 10^5 M_\odot$, $r_1^0 = 6.000$ kpc, $\Delta r_1^0 = 0.6242$ kpc, $\Delta\varphi_1^0 = 3.4907$, with significant figures up to the resolution of the central and right panels). The top row shows the relative coordinates between the BHs during the capture. The reference frame is centered on the most massive black hole and rotates with it. The binary Hill radius is indicated by the dotted circle. In the bottom row, the black line represents the binding energy ($E_b$) per unit mass of the binary ($M_{\text{tot}} = M_1 + M_2$), with negative values depicted in the colored area. The black marker indicates the point at which the energy becomes negative. The blue curve represents the relative distance ($d_{\text{rel}}$) between the BHs. The binary Hill radius is shown by a dashed blue line. The reference capture involves 5 close encounters between the black holes. Central panel: $N_e$ variations as a function of initial conditions ($M_1, M_2$). The ($M_1, M_2$) values are distributed plus or minus 1% around the initial condition of the reference capture, denoted as ($M_1^0, M_2^0$). Right panel: $N_e$ variations as a function of initial conditions ($r_1, \Delta r_1$). The ($r_1, \Delta r_1$) values are distributed plus or minus 1% around the initial condition of the reference capture, denoted as ($r_1^0, \Delta r_1^0$). The reference capture is located at (0, 0) in each panel. The strong variations of $N_e$ with initial conditions is suggestive of chaotic motion.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Percentage of initial captures as a function of different parameter pairs, marginalising each time over the three others. Top-left panel: as a function of the mass of the inner BH ($M_1$) and its initial radius ($r_1$); bottom-left panel: as a function of the mass of the outer BH ($M_2$) and $r_1$; top-right panel: as a function of $M_1$ and the initial radial separation between the BHs ($\Delta r_1$); bottom-right panel: as a function of $M_2$ and $\Delta r_1$.</figcaption></figure>

took for the capture to occur since the start of the simulation. The average capture time shows an increase as the total BH mass diminishes.

Figure 10 represents the percentage of captures as a function of ($r_i, \Delta r_i$), when we fix the masses $M_1, M_2$ and marginalize over $\Delta\varphi_i$ only. We show the results for two different pairs of masses,
<figure><img src="image_10.png" /><figcaption>Fig. 9. Percentage of dynamical captures as a function of inner ($M_1$) and outer ($M_2$) BH mass. The contour lines represent the average time it took for the capture to occur since the start of the simulation, measured in Gyr. The percentage of captures increases as the mass of the two BHs increases, and the mass of the outer BH is greater than the mass of the inner one ($M_2 > M_1$).</figcaption></figure>

on the left-hand panel $log(M_{tot}/M_{\odot}) = 4.8$ and on the right panel $log(M_{tot}/M_{\odot}) = 6.1$. Red areas represent initial captures and blue areas dynamical captures. Note the presence of blue bands that shift to the right ($\Delta r_i$ larger) as the total mass increases. To the left of these bands ($\Delta r_i$ smaller), we notice the presence of initial captures on a shallow red background. The plot can be divided into three regions: to the left, within, and to the right of the bands.

To the right of the band ($\Delta r_i$ larger), there is no capture: the BHs are too distant initially to reach a configuration where their gravitational attraction is comparable to the force generated by the external potential. To the left of the band ($\Delta r_i$ smaller), captures occur at $t = 0$ (in red). Indeed, two possible scenarios arise: (i) captures form at $t = 0$ if BHs have a sufficiently low initial phase shift $\Delta\phi_i$, (ii) the initial phase shift is too large for an initial capture, and $\Delta r_i$ is too small to allow phase catch-up within the simulation timeframe<fnref n="4" />. Also note the darker red area in the lower left of the right panel of Fig. 10. In this region the BHs are positioned very near the galactic center ($r_i$ small) and are in close radial proximity to each other. Under these conditions, the shorter distances for the same angular separation allow for more flexibility in terms of the initial phase shifts that result in Jacobi capture. In other words, when $r_i$ decreases larger initial phase shifts can lead to a capture, which explains the increased capture percentage.

## 3.2. Probability of capture

The previous section provides valuable insights into the impact of various parameters on Jacobi captures. In this section, our objective is to reduce the number of free parameters by making two assumptions concerning the origin and decentering of BHs. This allows us to calculate a probability of capture.

### 3.2.1. Bound black holes from the core stalling radii

In order to reduce the number of free parameters, we first make the assumption that the BHs are off-center due to core stalling, which sets $r_i$ and $\Delta r_i$ (cf. Sect. 2.3), leaving only

three free parameters (parameter ranges remaining unchanged). By employing this approach, we analyzed 84 000 simulations, 11 712 of which show a Jacobi capture, corresponding to 13.9% of cases. There is an increase in the percentage of captures, meaning that core stalling allows BHs to be placed on radii that favor capture.

The left panel of Fig. 11 shows the percentage of captures as a function of the mass of the inner ($M_1$) and outer ($M_2$) BH. The upper left region of the figure is left empty because of the condition $M_2 \ge M_1$, guaranteeing that BH2 consistently holds the outer position (cf. Sect. 2.3). The right-hand panel depicts the percentage of captures as a function of the initial phase difference $\Delta\phi_i$. The left panel reveals a higher percentage of captures when both $M_1$ and $M_2$ exhibit higher values, and when the mass gap between BHs is modest. The right-hand panel illustrates that captures predominantly occur with small initial phase difference ($\Delta\phi_i$ close to 0).

Core stalling positions the BHs on very similar radii, as can be seen in Fig. 10, where the black lines indicate the boundaries within which core stalling operates for the studied mass range of $10^3-10^6 M_{\odot}$. This region is in the regime of small radial separation discussed in Sect. 3.1. In that region, the vast majority of captures (95.4%) occur initially, when the BHs are sufficiently close in both initial radius and initial phase shift. Furthermore, to understand why there is a higher percentage of capture when masses are high, we can once again refer to Fig. 10. In each panel, the black cross indicates the positioning of BHs by core stalling. Notably, in the right-hand panel with higher masses, this cross falls within a darker region. This zone shifts from left to right as the mass of the BHs increases. In the right panel of Fig. 11, captures are observed at low angular separations, as we simulate 14 Gyr of evolution to evaluate captures within a Hubble time. When the angular separation is too large, BHs are unable to catch up with their phase delay within the given timeframe.

### 3.2.2. Black holes from a cosmological sequence of mergers

It is possible to reduce further the parameter set by employing a statistical sampling approach for BH masses, based on a cosmological model. This leaves a single parameter $\Delta\phi_i$ assumed to be uniform, allowing the derivation of a probability of capture. The selection of BH masses hinges on the assumption that they originate from the nuclei of two galaxies that underwent major mergers, ultimately forming our target galaxy for study. We also assume that within these nuclei, the BHs were accreting matter from the galaxy, which makes it possible to correlate their mass with that of their host galaxy.

First, we used dark matter merger trees obtained within the Semi-Analytic Satellite Generator (SatGen) introduced by Jiang et al. (2021). These merger trees are generated using an algorithm (Parkinson et al. 2008) based on the extended Press-Schechter formalism (Lacey & Cole 1993). We generated 4000 such merger trees for the target halo of mass $M = 10^{11} M_{\odot}$ at $z = 0$, and assess the associated stellar masses of the two progenitors according to the stellar-to-halo mass relation of Behroozi et al. (2013). Given the 0.16 dex scatter of the relation, we drew 10 realizations of the stellar mass per halo mass, following a log-normal distribution. In the context of accreting BH, mass accumulation is tied to the stellar mass. We assumed that the BH masses follow the $M_{BH} - M_{\star}$ relation of Greene et al. (2020)<fnref n="5" />:

$$ \log(M_{BH}/M_{\odot}) = \alpha + \beta \log(M_{\star}/M_{0}) \pm \epsilon \quad (5) $$

<footnote n="4">⁴ When $\Delta r_i$ is small, the BHs share similar orbital frequencies, so their phase difference evolves slowly.</footnote>

<footnote n="5">⁵ Low-mass galaxies pose challenges as they lack $\sigma_{\star}$ measurements and have diminished bulge fractions (MacArthur et al. 2003). Conse-</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Percentage of captures as a function of initial inner BH radius ($r_i$) and initial radial separation ($\Delta r_i$) for two mass pairs. On the left panel, $\log(M_1/M_\odot) = 3.9$ and $\log(M_2/M_\odot) = 4.7$, on the right panel, $\log(M_1/M_\odot) = 5.6$ and $\log(M_2/M_\odot) = 6$. Red areas represent initial captures while blue areas dynamical captures. The black lines indicate the limits within which core stalling can position our BHs, given the mass range studied. The black crosses indicate where core stalling would position the BHs in each case.</figcaption></figure>

with $M_0 = 3 \times 10^{10} M_\odot$, $\alpha = 7.43$, $\beta = 1.61$ and $\epsilon = 0.81$ the intrinsic scatter. We drew a series of 10 BH masses per stellar mass assuming again a lognormal distribution given the scatter of the relation. Following a selection process that involves filtering major mergers, requiring progenitors' halo mass ratio to be above 1/3, we end up with 4367 pairs of BH masses ($M_1, M_2$). A visualization of this distribution is plotted in Fig. 12 and is also reported as contour lines in Fig. 11. We then convolved the merger probability (Fig. 11) with the mass distribution derived from the semi-analytic model (Fig. 12). For each pair of masses ($M_1, M_2$) generated by the monte carlo algorithm, we calculated the capture proportion using our results from Sect. 3.2.1 as

$$p_i(M_1, M_2) = \frac{n(\text{bound} | M_1, M_2)}{n(\text{tot} | M_1, M_2)} \quad (6)$$

with $p_i(M_1, M_2)$ the capture proportion for a given pair of masses, $n(\text{bound}|M_1, M_2)$ the number of captures for the pair of masses ($M_1, M_2$) and $n(\text{tot}|M_1, M_2)$ the total number of simulations conducted for this pair of masses. We finally estimated the total probability:

$$p_{\text{tot}} = \frac{1}{N} \sum_{i=1}^{N} p_i \quad (7)$$

where the sum relates to the cases generated by SatGen, $N$ being the number of BH pairs generated. We obtain: $p_{\text{tot}} = 13.2\%$. This result, which is not negligible, carries significant implications for the assembly of BH masses, especially if these captures can be sustained over time through dissipative forces. Further investigations into the effectiveness of these forces, particularly within stripped nuclei or globular clusters, are warranted.

# 4. Conclusions

## 4.1. Summary

As a result of prolonged dynamical friction timescales and core stalling in cored dark matter haloes, a significant population of off-center intermediate-mass BHs is expected in dwarf galaxies. In this article, we investigate the possibility for such BHs to merge outside the galactic center. Dwarf galaxies being the most common type of galaxy in the Universe, off-center mergers could have a major impact on the assembly of BH mass and gravitational wave emission. We focus on the first step of an

off-center merger of IMBHs: Jacobi capture. In order to isolate the capture step, we perform idealized simulations of two BHs in coplanar circular orbits in a smooth background gravitational potential. Given the cored background potential and the fact that the BHs may not be far from the center, the Lagrange points and the zone of influence differ significantly from the Keplerian case (cf. Sect. 2.4). We find that the captures are highly sensitive to initial conditions, suggesting chaotic dynamics (cf. Sect. 2.6 and Fig. 7).

First in Sect. 3.1, we aim to analyze which configurations, in terms of BH masses and kinematics, favor Jacobi captures. This work allows to pinpoint regions in phase space where Jacobi captures are most likely to occur, paving the way for a comprehensive investigation of the subsequent merger stages through N-body simulations. We find that:

- In 7.1% of our cases of study (see Table 1), Jacobi captures occur. These captures are categorized into two distinct types: those formed at the beginning of the simulations, referred to as initial captures, and those emerging during the integration at later times, which we call dynamical captures. However, due to the ergodic nature of this study, any stretch of 14 Gyr should be equivalent to another and a time shift in the simulation range would interchange initial and dynamical captures.

- Initial captures represent 66.3%. They are favored at high radii due to the predominance of the potential component in our capture criterion over energy (Eq. (4)). Conversely, at low radii, the decrease in the circular velocity shear, attributed to the core within the density profile, contributes to their occurrence (see Fig. 8).

- Dynamical captures make up the remaining 33.7% of the captures. They are favored when the mass of the outer BH exceeds that of the inner one (see Fig. 9), which is of particular interest within the context of core stalling. We find capture bands in $\Delta r_i$, i.e., a range of radial separation values leading to capture. The capture band moves to higher $\Delta r_i$ as the binary mass increases. To the right of the band ($\Delta r_i$ larger), the BHs are too far away from each other to produce a capture. To the left of the band ($\Delta r_i$ smaller), dynamical capture is infeasible due to the impossibility of phase catch-up within a Hubble time.

Then, in Sect. 3.2 we reduce the number of free parameters in order to estimate a probability of capture. To reach this goal, we made assumptions about the origin of the BHs and the nature of their shift from the center.

- Our first assumption is that BH would be placed on their orbit by core stalling, which fixes their initial radius. We find that in 13.9% of our cases of study (see Table 2), Jacobi captures

<footnote>quently, Greene et al. (2020) revisited the $M_{BH} - M_\star$ correlation, focusing on dynamical studies involving BHs with $M_{BH} < 10^6 M_\odot$.</footnote>
<figure><img src="image_12.png" /></figure>

**Fig. 11.** Percentage of captures in the three-parameter case. Left panel: as a function of inner ($M_1$) and outer ($M_2$) BH mass when positioned on their core stalling radius. The contour lines depict the probability (in percentage) of obtaining a mass pair (derived from a Monte Carlo sampling using SatGen, see Sect. 3.2.2). Right panel: as a function of the initial phase difference.

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Monte Carlo sampling of BH mass using SatGen.</figcaption></figure>

occur. The increase in the capture percentage compared to the previous scenario suggests that core stalling positions BHs at radii that favor capture.

– Given the mass range under study, core stalling positions the BHs within a regime characterized by close radial separations. In this context, the vast majority of captures (95.4%) occur initially, if the BHs are sufficiently close. This condition necessitates $M_1$ and $M_2$ to be close, corresponding to a small radial separation, and a low initial phase shift $\Delta\phi_i$ (see Fig. 11).

– Our second assumption is that these BHs originate from the central regions of two galaxies that underwent a major merger, ultimately forming our target galaxy for study. This enables to recover a statistical distribution of BH masses, employing a Monte Carlo sampling algorithm to construct merger trees and using scaling relations between stellar and BH masses (see Fig. 12). By doing this, we obtain a probability of capture of 13.2%.

## 4.2. Caveats

The primary goal of our work is to gain a better qualitative understanding of Jacobi capture on a case-by-case basis. On the way, we had to make two particular assumptions that will need to be alleviated in the future, as they may significantly impact our picture for BH growth: we restricted the problem to specific initial orbits, and we assumed a static galactic potential.

The first restrictive assumption on BH orbits concerns the geometry of the problem. Indeed, we considered only coplanar BHs, while in all generality their orbits could be randomly oriented. It may seem that this assumption artificially boosts our probability of Jacobi capture, because coplanar orbits imply that the BHs are, on average, closer both spatially, and in velocity. However, in the scenario where off-center BHs originate from accretion events, it is expected that such occurrences demonstrate spatial correlation. Consequently, the BH orbital planes have correlated orientations, in a way that favors coplanar orbits. The second assumption on the BH orbits is that they are initially circular. But this is actually quite realistic, because dynamical friction tends to circularize a wandering BH's orbit, all the more that it has been orbiting for a long time. Finally, we assumed all orbits to revolve in the same sense, which effectively boosts the probability of capture by a factor 2. Similarly, this factor could be mitigated by the correlated origin of accreted BHs.

The assumption of a static background galaxy may also be a concern. Indeed, in some of our cases of study, we expect the reflex response from the galaxy to be significant. This is particularly the case for a stalling BH, where the very definition of stalling implies that the galaxy is responding. This motion of the background may induce intricate modifications of the dynamics of individual setups, even at the level of a single BH orbiting a galaxy. The consequences could be even stronger in the case of two BHs. Then, dynamical friction and buoyancy exerted on each individual BH could couple, in a way that would even modify the stalling radii of each BH, or suppress the stalling process altogether. A more detailed treatment of such cases with live N-body runs is required to assess the influence of the galactic response on Jacobi capture.

In Sect. 3.2.2 we have selected the BH masses from a semi-analytical model, it is important to note that it is contingent upon the assumption that they are remnants of central BHs from past major mergers. The application of a correlation between stellar mass and BH mass hinges on the assumption of BH accretion of
galactic matter. For scenarios involving wandering BHs situated outside the galactic center, refraining from accretion, reference to Fig. 11 is required on a case-by-case basis for each BH mass pair.

### 4.3. Implications for black holes

In this study we investigate if and under which circumstances two BHs can become bound while outside the center of a galaxy. We do find that such captures do happen at substantial fractions in our simplified simulations. In recent years, IMBHs were successfully detected outside the centers of galaxies (Reines et al. 2020), as well as single relic nuclear clusters at kpc from the center of the galaxy have been confirmed (Seth et al. 2014; Ahn et al. 2018). Simulations and model predictions show that these relics of the accretion event should be extremely common (Tremmel et al. 2018; Voggel et al. 2019) and mergers have been shown to be the reason for the M-σ relation (Volonteri & Natarajan 2009).

If the hidden off-center population of BHs merge repeatedly before reaching the center they would increase their mass with smaller merging events and be more massive once they reach the center of the main galaxy. The existence of off-center BH mergers also implies that three massive BHs simultaneously descending toward the galaxy's center becomes more rare. This would cause an unstable three-body system ensues, which can lead to the ejection of one BH. However, if two of the three BHs merge during their inspiral path, only two BHs reach the center, enabling a more stable situation and thus a more likely subsequent merger that yields a more massive BH than in the previous scenario. The history of BH mergers being altered by the possibility of off-center mergers, this should change the expected gravitational wave signals.

As this type of mergers can have major implications on the mass growth of BHs but has not been explored much in the literature it is crucial to improve upon our idealized simulations. The next step in the investigation of the off-center merger phenomena could be to determine whether the inclusion of stellar populations associated with BHs, such as stripped galactic nuclei or stars within globular clusters changes the Jacobi captures. Such an inclusion is necessary to be more realistic, as typically BHs are surrounded by nuclear star clusters (Seth et al. 2008). In particular it will be important to see whether dissipative forces such as dynamical friction between the stellar populations can help sustain these captures.

**Acknowledgements.** We thank Raphaël Errani, Florent Renaud, Paolo Bianchini, Katarina Kraljic, Jorge Peñarrubia, Anna Lisa Varri, Marta Volonteri, Benoit Famaey, Gary Mamon, Karamveer Kaur and Julien Montillaud for the helpful discussion and the anonymous referee for the constructive suggestions. The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsac calcul/Big Data. This research made use of the following software: Agama (Vasiliev 2018), Jupyter (Granger & Perez 2021), Matplotlib (Hunter 2007), Numpy (Harris et al. 2020), Python (Van Rossum & Drake 1995), SatGen Jiang et al. (2021) and Scipy (Virtanen et al. 2020). Simulations in this paper made use of the REBOUND N-body code (Rein & Liu 2012). The simulations were integrated using IAS15, a 15th order Gauss-Radau integrator (Rein & Spiegel 2015).

## References

* Ahn, C. P., Seth, A. C., Cappellari, M., et al. 2018, ApJ, 858, 102

* Banik, U., & van den Bosch, F. C. 2021, ApJ, 912, 43

* Banik, U., & van den Bosch, F. C. 2022, ApJ, 926, 215

* Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307

* Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

* Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2019, MNRAS, 482, 2913

* Bellovary, J. M., Hayoune, S., Chafla, K., et al. 2021, MNRAS, 505, 5129

* Binney, J. J., & Evans, N. W. 2001, MNRAS, 327, L27

* Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press)

* Boekholt, T. C. N., Rowan, C., & Kocsis, B. 2023, MNRAS, 518, 5653

* Boily, C. M., Padmanabhan, T., & Paiement, A. 2008, MNRAS, 383, 1619

* Borriello, A., & Salucci, P. 2001, MNRAS, 323, 285

* Borukhovetskaya, A., Errani, R., Navarro, J. F., Fattahi, A., & Santos-Santos, I. 2022, MNRAS, 509, 5330

* Cann, J. M., Satyapal, S., Bohn, T., et al. 2020, ApJ, 895, 147

* Chandrasekhar, S. 1943, ApJ, 97, 255

* Cole, D. R., Dehnen, W., Read, J. I., & Wilkinson, M. I. 2012, MNRAS, 426, 601

* de Blok, W. J. G., McGaugh, S. S., Bosma, A., & Rubin, V. C. 2001, ApJ, 552, L23

* DeLaurentiis, S., Epstein-Martin, M., & Haiman, Z. 2023, MNRAS, 523, 1126

* Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS, 172, 15

* Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9

* Freundlich, J., Jiang, F., Dekel, A., et al. 2020, MNRAS, 499, 2912

* Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13

* Goerdt, T., Moore, B., Read, J. I., Stadel, J., & Zemp, M. 2006, MNRAS, 368, 1073

* Goldreich, P., Lithwick, Y., & Sari, R. 2002, Nature, 420, 643

* Granger, B. E., & Perez, F. 2021, Comput. Sci. Eng., 23, 7

* Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

* Heggie, D. C. 1975, MNRAS, 173, 729

* Henon, M., & Petit, J. M. 1986, Cel Mech., 38, 67

* Hernquist, L. 1990, ApJ, 356, 359

* Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

* Inoue, S. 2011, MNRAS, 416, 1181

* Jiang, F., Dekel, A., Freundlich, J., et al. 2021, MNRAS, 502, 621

* Kaur, K., & Sridhar, S. 2018, ApJ, 868, 134

* Kollatschny, W., Weilbacher, P. M., Ochmann, M. W., et al. 2020, A&A, 633, A79

* Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581

* Lacey, C., & Cole, S. 1993, MNRAS, 262, 627

* Lanzoni, B., Mucciarelli, A., Origlia, L., et al. 2013, ApJ, 769, 107

* Li, J., Lai, D., & Rodet, L. 2022, ApJ, 934, 154

* Lützgendorf, N., Kissler-Patig, M., Gebhardt, K., et al. 2013a, Mem. Soc. Astron. It., 84, 129

* Lützgendorf, N., Kissler-Patig, M., Neumayer, N., et al. 2013b, A&A, 555, A26

* MacArthur, L. A., Courteau, S., & Holtzman, J. A. 2003, ApJ, 582, 689

* Mackey, A. D., & Gilmore, G. F. 2003, MNRAS, 345, 747

* McNamara, B. J., Harrison, T. E., Baumgardt, H., & Khalaj, P. 2012, ApJ, 745, 175

* Mezcue, M., Civano, F., Marchesi, S., et al. 2018, MNRAS, 478, 2576

* Mezcue, M., Suh, H., & Civano, F. 2019, MNRAS, 488, 685

* Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2019, ApJ, 872, 104

* Noyola, E., Gebhardt, K., & Bergmann, M. 2008, in Dynamical Evolution of Dense Stellar Systems, eds. E. Vesperini, M. Giersz, & A. Sills, 246, 341

* Pardo, K., Goulding, A. D., Greene, J. E., et al. 2016, ApJ, 831, 203

* Parkinson, H., Cole, S., & Helly, J. 2008, MNRAS, 383, 557

* Peñarrubia, J. 2023, MNRAS, 519, 1955

* Pechetti, R., Seth, A., Kamann, S., et al. 2022, ApJ, 924, 48

* Petit, J. M., & Hénon, M. 1986, Icarus, 66, 536

* Petit, J. M., & Hénon, M. 1990, in Dynamics and Stochastic Processes Theory and Applications, eds. R. Lima, L. Streit, & R. Vilela Mendes, 355, 225

* Petts, J. A., Gualandris, A., & Read, J. I. 2015, MNRAS, 454, 3778

* Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101

* Read, J. I., Goerdt, T., Moore, B., et al. 2006, MNRAS, 373, 1451

* Rein, H., & Liu, S. F. 2012, A&A, 537, A128

* Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424

* Reines, A. E., Greene, J. E., & Geha, M. 2013, ApJ, 775, 116

* Reines, A. E., Condon, J. J., Darling, J., & Greene, J. E. 2020, ApJ, 888, 36

* Renaud, F., Gieles, M., & Boily, C. M. 2011, MNRAS, 418, 759

* Rowan, C., Boekholt, T., Kocsis, B., & Haiman, Z. 2023, MNRAS, 524, 2770

* Samsing, J., MacLeod, M., & Ramirez-Ruiz, E. 2018, ApJ, 853, 140

* Schweizer, F., Seitzer, P., Whitmore, B. C., Kelson, D. D., & Villanueva, E. V. 2018, ApJ, 853, 54

* Seth, A., Agüeros, M., Lee, D., & Basu-Zych, A. 2008, ApJ, 678, 116

* Seth, A. C., van den Bosch, R., Mieske, S., et al. 2014, Nature, 513, 398
Shen, Y., Chen, Y.-C., Hwang, H.-C., et al. 2021, Nat. Astron., 5, 569

Stemo, A., Comerford, J. M., Barrows, R. S., et al. 2021, ApJ, 923, 36

Tagawa, H., Saitoh, T. R., & Kocsis, B. 2018, Phys. Rev. Lett., 120, 261101

Tagawa, H., Haiman, Z., & Kocsis, B. 2020, ApJ, 898, 25

Tollet, E., Macciò, A. V., Dutton, A. A., et al. 2016, MNRAS, 456, 3542

Tremmel, M., Governato, F., Volonteri, M., Pontzen, A., & Quinn, T. R. 2018, ApJ, 857, L22

van den Bosch, R. C. E. 2016, ApJ, 831, 134

van der Marel, R. P., & Anderson, J. 2010, ApJ, 710, 1063

Van Rossum, G., & Drake, Jr., F. L. 1995, Python Reference Manual (Centrum voor Wiskunde en Informatica Amsterdam)

Vasiliev, E. 2018, ArXiv e-prints [arXiv:1802.08255]

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Voggel, K. T., Seth, A. C., Baumgardt, H., et al. 2019, ApJ, 871, 159

Voggel, K. T., Seth, A. C., Baumgardt, H., et al. 2022, A&A, 658, A152

Volonteri, M., & Natarajan, P. 2009, MNRAS, 400, 1911

Zaw, I., Rosenthal, M. J., Katkov, I. Y., et al. 2020, ApJ, 897, 111

Zhao, H. 1996, MNRAS, 278, 488
# Appendix A: Binary eccentricity during Jacobi captures

The histogram in the left panel of Fig. A.1 displays the mean eccentricity of the binary during captures. Jacobi captures tend to result in rather eccentric binaries, a factor that can aid in stabilization through energy dissipation caused by dissipative forces. Post-Newtonian corrections can particularly influence highly eccentric orbits, especially at the pericenter (Samsing et al. 2018). Notably, our criterion ensures that the mean eccentricity during captures remains below 1. On the right panel of Fig. A.1, the standard deviation of eccentricity during captures is depicted. These variations may permit a very high eccentricity during the capture, despite the lower mean eccentricity.

<figure><img src="image_14.png" /></figure>

**Fig. A.1.** Eccentricity during Jacobi captures. Left panel: histogram in percent of the mean eccentricity of the binary during captures (initial captures in red and dynamical captures in blue). The dark blue line represents only dynamical captures with a single close encounter (i.e., flybys). Jacobi captures produce rather eccentric binaries. A significant proportion are in the very high eccentricity range ($e \ge 0.9$), for which Post-Newtonian corrections can be significant at the pericenter. Right panel: percentage histogram of standard deviation of eccentricity during captures. These variations may permit a very high eccentricity during the capture, despite the lower mean eccentricity.

# Appendix B: Alternative capture criterion

In this section, we revisit some of the findings presented in the article using an alternative capture criterion based solely on the radius of influence (i.e., only on gravitational potential). Here, we define BHs as captured when their separation is less than the Hill radius of the binary. This approach allows us to emphasize the role of the kinetic component of the criterion utilized throughout the paper (Eq. 4).

When the separation between the BHs is smaller than the binary Hill radius, the potential generated by the BHs competes with that of the galaxy, resulting in a significant perturbation of the BH trajectories. Therefore, it is natural to consider that this criterion alone may suffice to define a capture. We examine here the differences that this approach entails.

Applying this new criterion to the general five-parameter case, we find captures in 16.4% of the simulated instances. This is more than double the 7.1% obtained with the energy criterion. There are 5.4% initial captures (compared to 4.7% with the energy criterion) and 11.0% dynamical captures (compared

<figure><img src="image_15.png" /></figure>

**Fig. B.1.** Percentage of initial captures as a function of different parameter pairs, marginalising each time over the three others. Top-left panel: as a function of the mass of the inner BH ($M_1$) and its initial radius ($r_1$); bottom-left panel: as a function of the mass of the outer BH ($M_2$) and $r_1$; top-right panel: as a function of $M_1$ and the initial radial separation between the BHs ($\Delta r_1$); bottom-right panel: as a function of $M_2$ and $\Delta r_1$. This figure is analogous to Fig. 8 with a different capture criterion. BHs are considered bound if their relative distance falls below the binary's Hill radius.

to 2.4%). We observe a significant increase in the number of dynamical captures.

Figure B.1 serves as the counterpart to Fig. 8, depicting the percentage of initial captures across various parameter pairs, marginalized over the remaining three. Unlike the latter, this figure does not consider the relative velocity between BHs. In the top-left panel, the percentage of captures with respect to $r_1$ closely mirrors the Hill radius curve depicted in Figure 5, albeit with a leftward shift as the BH mass decreases.

The white area in the top-left corner of Fig. B.1 appears dark on the figure 8 due to cases where the BHs initially maintain a relative distance slightly larger than the Hill radius, yet exhibit very low relative velocities. They are energetically close and plunge into the Hill sphere after a few time steps. Under this criterion, they are classified as dynamical captures, while their low relative velocity categorizes them as initial captures according to the energy criterion.

Figure B.2 serves as the counterpart to Fig. 10, depicting the percentage of dynamical captures relative to the initial radius of the inner BH ($r_1$) and the initial radial separation between the BHs ($\Delta r_1$). Comparing it with Figure 10, we observe that the dynamical binding bands extend further down in $r_1$. To understand the nature of these dynamical captures introduced by the alternative criterion, we can turn to Fig. B.3, mirroring Figure A.1. The histogram of mean eccentricities during captures on the left-hand panel extends beyond 1. Consequently, these newly added dynamical captures are mostly hyperbolic flybys close to the galactic center and for which BHs have high relative velocities.
<figure><img src="image_16.png" /></figure>

**Fig. B.2.** Percentage of dynamical captures as a function of initial inner BH radius ($r_1$) and initial radial separation ($\Delta r_1$) for two mass pairs, On the left panel, $\log(M_1/M_\odot) = 3.9$ and $\log(M_2/M_\odot) = 4.7$, on the right panel, $\log(M_1/M_\odot) = 5.6$ and $\log(M_2/M_\odot) = 6$. This figure is analogous to Fig. 10 with a different capture criterion. BHs are considered bound if their relative distance falls below the binary's Hill radius. We observe that the binding bands of dynamical captures extend closer to the galactic center. The majority of these captures involve hyperbolic flybys, indicating that the BHs involved have high relative velocities.

<figure><img src="image_17.png" /></figure>

**Fig. B.3.** Eccentricity during Jacobi captures. Left panel: histogram in percent of the mean eccentricity of the binary during captures (initial captures in red and dynamical captures in blue). The dark blue line represents only dynamical captures with a single close encounter (i.e., flybys). Right panel: percentage histogram of standard deviation of eccentricity during captures. This figure is analogous to Fig. A.1 with a different capture criterion. BHs are considered bound if their relative distance falls below the binary's Hill radius. Note that the alternative criterion adds captures for which the mean eccentricity is greater than 1.

## Appendix C: Pacman orbit and number of close encounters

The existence of a central core within the density profile facilitates the emergence of a novel orbital type known as Pacman orbits (Banik & van den Bosch 2022). These orbits play a key role in the phenomenon of dynamical buoyancy within cores, which drives perturbers outward to an equilibrium radius with the dynamical friction, ultimately resulting in core stalling. During capture events, BHs experience these unique dynamics. Figure C.1 illustrates an example of a Pacman capture.

<figure><img src="image_18.png" /></figure>

**Fig. C.1.** Example of Pacman orbit. Top row: relative coordinates between BHs during the Jacobi capture. The reference frame is centered on the most massive BH and rotates with it. The binary Hill radius is marked by the dotted circle. Bottom row: the black line represents the binding energy ($E_b$) per unit mass of the binary ($M_{\text{tot}} = M_1 + M_2$), with negative values shown in the colored area. The black marker indicates the point when the energy becomes negative. The blue curve represents the relative distance ($d_{\text{rel}}$) between the BHs. The binary Hill radius is indicated by a dashed blue line.

This example serves as a perfect illustration of how the labeling of the number of close encounters between BHs depends, in some borderline cases like this one, on the criterion used to label a capture. In the bottom panel of Fig. C.1, the binding energy and relative distance between the BHs during the capture are depicted. With the energy-based capture criterion (Eq. 4), we would label two captures with a single close encounter, as the energy returns to positive values between the encounters. However, if we utilize a criterion solely based on distance, as outlined in Section B (where BHs are labeled as captured when the relative distance is less than the binary Hill radius), we would label only a single capture with two close encounters. This case underscores the complexity of selecting an appropriate capture criterion, as finding a perfect one remains a challenging task.