text:
# Collisional dynamics around binary black holes in galactic centers

## 1 Introduction

Currently the standard picture of galaxy formation involves the collapse of baryonic matter in hierarchically clustering dark matter halos and the subsequent building of big galaxies from small ones via merging processes e.g., (Peebles, 1993; Diaferio et al., 1999; Kauffmann et al., 1999a, b). While recent cosmological simulations can adequately reproduce many global properties of galaxies and their correlations, the details are still very much dependent on the gas physics and stellar feedback involved (see e.g., Navarro and Steinmetz (2000)). Additionally, most, if not all, galaxies harbor supermassive black holes in their center (Magorrian et al., 1998; Richstone et al., 1998; Kormendy and Richstone, 1995). Correlations have been recently detected between black hole masses, galaxy masses, and central velocity dispersions in galaxies (Ferrarese and Merritt, 2000; Gebhardt et al., 2000). These correlations are strong evidence that black holes in galactic nuclei are linked to the dynamical history of their host galaxies. Haehnelt and Kauffmann (2000) and Kauffmann and Haehnelt (2000) demonstrate how this is consistent with the framework of semi-analytic models that follow the formation and evolution of galaxies in a cold dark matter-dominated universe. They assume supermassive black holes are formed and fueled during major mergers, qualitatively explaining many aspects of the observed evolution of galaxies, including the observed relation between bulge luminosity, velocity dispersion, and central black hole mass. As already discussed by Begelman et al. (1980), such a scenario requires the formation of galactic nuclei containing at least two black holes, depending on the black hole merger rate relative to the galaxy merger rate. However, there is very little observational evidence for massive black hole binaries (Lehto and Valtonen, 1996; Halpern and Eracleous, 2000). This conflict between theory and observations has become known as the “sinking black hole problem”. As an alternative to minimally impacting stellar dynamical processes, Gould and Rix (2000) and Armitage and Natarajan (2002) have proposed mechanisms which lead to rapid decay of massive black hole orbits and subsequent black hole mergers in galactic centers. Also, Begelman et al. (1980) offered the solution that gas accretion could dominate the orbital decay in the intermediate phase of the sinking black hole problem when dynamical friction becomes inefficient. However, as we will discuss later, dynamical friction, as laid out by Chandrasekhar (1943), is not sufficiently effective by itself to lead to rapid coalescence of black hole binaries.

If there are no quick mergers, multiple black hole nuclei could lose black holes through slingshot ejections (Valtonen et al., 1994). Once a binary system becomes hard, the high orbital velocities of the black holes allow further hardening through close encounters and three-body interactions with stars. Such processes will evacuate field stars from the surroundings of the binary, therefore it can be argued that the stellar scatterings cannot produce rapid coalescence. The preceding argument assumes that the center of mass of the binary does not move with respect to the stellar system. However, we will show that even with a fairly symmetrical initial setup the binary gains some linear momentum. This introduces a wandering motion which exceeds the expectations from equipartition. The wandering of the binary guarantees an adequate supply of stars for binary hardening and rapid coalescence through purely stellar dynamical processes.

Our new computational method allows us to study in detail three-body interactions of a black hole binary with field stars. Although one may argue that the perturbing mass of the field stars is small compared to the black hole mass and should have negligible impact, there are many stars, and each encounter can lead to changes in binding energy and eccentricity of the black hole binary. In fact, our models show that the black hole binary keeps a rather high eccentricity due to the encounters. Thus high eccentricity will speed up gravitational radiation mergers very efficiently, and is, as noted by Gould and Rix (2000) and Armitage and Natarajan (2002), a way to expedite massive black hole mergers in a purely stellar dynamical way.

The correct theoretical prediction of the frequency of black hole mergers in galactic environments will be important in the search for gravitational waves. The merging of supermassive black holes of $`3\times 10^4`$ to $`3\times 10^7\mathrm{M}_{}`$ in the nuclei of merging galaxies and protogalaxies can be detected with high signal-to-noise at redshifts from $`0<z<100`$ (Phinney, 2000) by the Laser Interferometer Space Antenna (LISA) (Danzmann, 2000).

Previous attempts to quantify this prediction have been made by either solving the perturbed two and three-body problem in simplified models (Mikkola and Valtonen, 1992), direct $`N`$-body models (Makino et al., 1993; Makino, 1997), or a combination of the two (Merritt and Quinlan, 1998; Quinlan and Hernquist, 1997). Simulating binary black hole hardening is extremely challenging, algorithmically and computationally. Since the mass differences between the black holes and the stars is so large, high particle numbers are required in order to model the relaxation processes around the black holes accurately. The simulations have used softened particles on special purpose computers (Makino et al., 1993; Makino, 1997) or a hierarchical hybrid code in which all forces involving the black hole particles are Keplerian (Merritt and Quinlan, 1998; Quinlan and Hernquist, 1997). These schemes used particle numbers in the order of $`10^4`$.

In this paper, we describe a new hybrid field-particle code which treats all particles with orbits crossing the central regions of the system with a high precision direct $`N`$-body method appropriate for collisional stellar dynamics. All other particles are integrated using a field method. In order to adapt both parts of the hybrid code to each other, the field method (approximating the potential exerted by a set of particles by a series expansion, referred to here as “SCF”) had to be upgraded to a fourth order Hermite integrator. This integration also uses the time derivative of the potential, as in modern direct $`N`$-body codes.

In the following sections some details of the sinking black hole problem are introduced. Section 2 introduces the integration software used for the numerical experiments described in this paper. Section 3 is devoted to a comparison between the new collisional code with a well used workhorse simulator in this field called *NBODY6* (Aarseth, 1993, 1996, 1999), using its parallel implementation *NBODY6++* (Spurzem and Baumgardt, 2001; Spurzem, 1999). The application of the code to the sinking binary black hole problem is reported in section 4.

## 2 Collisional stellar dynamics with EuroStar

A numerical simulation of the hardening phase (until the massive black holes start to radiate gravitational waves) must be able to accurately follow three-body encounters. For this reason, the Keplerian potential should not be softened in the denser parts of the system. Computationally, the central part of the system would best be treated using a collisional integrator. The code must be able to integrate encounters leading to large angle deflections in an efficient way, while requiring neither too much computing time nor introducing energy errors. The overall $`N`$-body integration does not need to be symplectic, but should keep the energy error as low as possible. On the other hand, in a system showing a core halo structure, the bulk of the stars in the halo move under the influence of the mean field of the whole cluster. The halo part of the central galactic cluster can be integrated with a mean field method.

In the new method (which we will refer to as EuroStar), both the collisional code NBODY6++ (Aarseth, 1993, 1999; Spurzem and Baumgardt, 2001) and the SCF method (Hernquist and Ostriker, 1992; Hernquist et al., 1995; Zhao, 1996; Sigurdsson et al., 1997; Holley-Bockelmann et al., 2001) are merged to optimize large-$`N`$ collisional $`N`$-body simulations. The star cluster, which is assumed to be in equilibrium, is divided into two sections by applying a critical angular momentum criterion. As shown in Figure 1, this allows for distinction between particles orbiting solely in the halo and ones which have trajectories leading through the core of the system.

In a stationary gravitational point mass system, two-body relaxation leads to an exchange between halo particles and core particles in a divided cluster. In a system of more than $`10^4`$ particles, only a few particles cross the core halo border per dynamical time. This is very fortunate because it allows us to integrate the orbits of the halo particles with a collisionless method and the core particles with a collisional code. Exchanges of particles cause energy conservation problems since the contribution of a particle to the main potential would be changed from Keplerian to a sample point of a mean field in EuroStar. This means that switching of particles from the core to the halo part of the integrator and vice versa is not permitted.

NBODY6++ integrates trajectories of point masses in the core of the system. It is an Aarseth-type direct force integrator applying the Hermite integration scheme. NBODY6++ gains its efficiency by implementing an Ahmad-Cohen neighbor scheme and individual block time steps (Ahmad and Cohen, 1973; Aarseth, 1999). Close interactions between particles are treated by regularization of the equations of motion (Kustaanheimo and Stiefel, 1965). NBODY6++ scales well on parallel computer systems and can also be used with the GRAPE special purpose computer (Spurzem and Kugel, 1999; Sugimoto et al., 1995). The Hermite scheme requires one to compute $`𝐅_i`$ and $`\dot{𝐅}_i`$ at each time step where,

$$𝐅_i=\underset{ji}{}\frac{Gm_j𝐫_{ij}}{r_{ij}^3},$$
(1)
$$\dot{𝐅}_i=\underset{ji}{}Gm_j\left[\frac{𝐯_{ij}}{r_{ij}^3}+\frac{3(𝐯_{ij}𝐫_{ij})𝐫_{ij}}{r_{ij}^5}\right].$$
(2)
The relative distance between particles $`i`$ and $`j`$ is given by $`𝐫_{ij}=𝐫_i𝐫_j`$. Accordingly, the relative velocity is $`𝐯_{ij}=𝐯_i𝐯_j`$. The extra effort of computing two direct force quantities allows one to approximate the particle’s orbit to fourth order. By storing $`𝐅_i`$ and $`\dot{𝐅}_i`$ from the previous time step it is possible to interpolate the next two higher derivatives and to apply a predictor/corrector scheme (Aarseth, 1996).

The SCF method qualifies for the collisionless part of a spherical system since then the basis functions are given analytically. This allows one to implement the Hermite scheme for SCF, which makes SCF an ideal far force extension to NBODY6++. Its drawback, however, is that this method restricts the input systems to have an approximately spherical particle distribution around the coordinate center. The possible asphericity depends on the number of spherical harmonics used for the potential expansion. In order to have better convergence, one-parameter basis functions for $`\rho _{nlm}(𝐫)`$ and $`\mathrm{\Phi }_{nlm}(𝐫)`$ are used (Zhao, 1996):

$$\rho _{nlm}(𝐫)=\sqrt{4\pi }\frac{K_{nl}r^lC_n^{(\omega )}(\xi )}{r^{(2\frac{1}{\alpha })}(1+r^{\frac{1}{\alpha }})^{2+\alpha (2l+1)}}Y_{lm}(\vartheta ,\phi ),$$
(3)
$$\mathrm{\Phi }_{nlm}(𝐫)=\frac{\sqrt{4\pi }r^lC_n^{(\omega )}(\xi )}{(1+r^{\frac{1}{\alpha }})^{\alpha (2l+1)}}Y_{lm}(\vartheta ,\phi ).$$
(4)
The $`C_n^{(\omega )}`$ are called ultraspherical or Gegenbauer polynomials. The spherical harmonics are given by the $`Y_{lm}(\vartheta ,\phi )`$. Once the $`A_{nlm}`$ for a certain set of particles are known, an analytic expression of the potential and the density is found. Due to the truncation of the expansion they represent the mean field and mean density. This means that the force at each position and its derivative can be computed using the following expressions,

$$F(𝐫)=\underset{nlm}{}A_{nlm}\mathrm{\Phi }_{nlm}(𝐫),$$
(5)
$$\dot{F}(𝐫)=\frac{d}{dt}\left(\underset{nlm}{}A_{nlm}\mathrm{\Phi }_{nlm}(𝐫)\right).$$
(6)
Since Equations (5) and (6) lead to a significant modification of the SCF scheme, we provide detailed form of these expressions in the Appendix.

## 3 Testing the hybrid code

As a first test for the new method, we have followed the last stages of an ongoing merger between two galaxies, each containing a supermassive black hole. For the initial setup, it is assumed that the stellar systems have already arranged themselves into a spherical system. The two, formerly central, supermassive black holes are moving through the stellar component with a speed on the order of the relative velocity between the two initial galaxies.

In the present simulations, the stellar component is a realization of a Plummer model. The density and potential of the spherically symmetric Plummer model are given by (Plummer, 1911):

$$\rho (r)=\frac{3M}{4\pi }\frac{R^2}{(R^2+r^2)^{5/2}},$$
(7)
$$\mathrm{\Phi }(r)=GM\frac{1}{(R^2+r^2)^{1/2}}.$$
(8)
The quantity $`M`$ describes the total mass of the system and $`G`$ is the gravitational constant. With the Plummer radius chosen to be $`R=3\pi /16`$, the half mass radius of this system is at a radius of $`r_\mathrm{h}0.78`$ in the model units of the simulations. The total mass of the system $`M`$ is set to unity. The gravitational constant $`G`$ is set to unity as well, conforming to the model units described by Heggie & Mathieu (1986). The stellar system is centered around the origin.

The black hole particles contain 1% of the system’s total mass. Their initial positions are at $`x=\pm 0.5`$, and their initial velocities are 13.6% of the circular velocity at their initial radii.

The black hole orbits are analyzed during the simulation assuming the orbit can be approximated by the classical two-body problem. The binding energies and eccentricities of the black hole orbits are computed from their relative distances and velocities, assuming a Keplerian potential. Once the black holes become bound, their two-body attraction is the most important force. The eccentricity of the binary $`ϵ`$ and the binding energy $`h`$ are computed as follows, using the definition $`𝐫=𝐫_a𝐫_b`$. The vectors $`𝐫_a`$ and $`𝐫_b`$ denote the position vectors of the black holes $`a`$ and $`b`$. From this it follows that (Boccaletti and Pucacco, 1996),

$$h=\frac{1}{2}m_{\mathrm{red}}\dot{r}^2+\frac{|𝐥|^2}{2m_{\mathrm{red}}|𝐫|^2}\frac{m_am_b}{|𝐫|},$$
(9)
$$a=\frac{m_am_b}{2h},$$
(10)
$$ϵ=\sqrt{1+\frac{2h|𝐥|^2}{m_{\mathrm{red}}(m_am_b)^2}},$$
(11)
where $`𝐯=\dot{𝐫}`$, and
$$\dot{r}=\frac{1}{|𝐫|}𝐯𝐫,$$
(12)
$$m_{\mathrm{red}}=\frac{m_am_b}{m_a+m_b},$$
(13)
$$𝐥=m_{\mathrm{red}}𝐫\times 𝐯.$$
(14)
This method of analysis provides a sensitive measure for the moment when the black holes become bound to one other. Furthermore, this way of analyzing the data also offers a precise tool for following the hardening of the binary.

Since this sinking binary problem is the first application of EuroStar, its results are compared with those of the fully collisional code NBODY6++. For the comparison runs, 16384 particles were simulated. Figure 2 shows the results for the two comparative runs. The plot on the left hand side in Figure 2 shows the eccentricity of the binary as a function of the simulated time in model units. The plot on the right hand side shows the two-body binding energy as a function of time. The binary becomes bound after 10 time units, in both cases.

The fully collisional method and the hybrid code EuroStar show slightly different sinking rates for the binary at the beginning of the simulation. These differences result from the different density of collisional particles around the black holes in both codes. While the black holes in the fully collisional run suffer small angle encounters with every particle in the system, this is not possible in the hybrid code. Naturally, all particles treated by the mean field method can only interact with the system through changes in the mean potential. After the binary has become bound, the hardening process is driven by the stars which have a small enough impact parameter such that they have an encounter timescale smaller than the orbital timescale of the massive binary. Therefore, the hardening depends more on the two-body encounters with neighboring particles, which have large orbital velocities.

With increasing simulation time, the binary locks into an oscillating motion around the center of mass of the stellar component. This motion does not extend beyond of the dense galactic core. Effectively, the differences in the density of the collisional particles between the two methods vanish after the binary has become bound. Figure 2 reflects this by showing a parallel evolution of eccentricity and binding energy for times larger than 10 time units in the simulation.

## 4 Hardening of a massive binary

This new code is intended for raising the total particle number for collisional simulations of spherical $`N`$-body systems. Hence the evolution of one, two or several massive bodies in a dense stellar cluster appears to be an ideal problem for EuroStar. This is why we are addressing here the problem of a sinking massive black hole binary in galactic centers. Another useful potential application os the dynamics of globular clusters.

### 4.1 Initial conditions

The particles representing the stellar component are distributed according to Plummer’s model with $`R=3\pi /16`$. The total mass of the stars is fixed at $`0.98`$, while the black hole particles carry $`0.01`$ each, so $`M=1.0`$. This is a fairly high mass for the black holes compared to the total mass of the stellar system, since Ferrarese and Merritt Ferrarese and Merritt (2000) found the black hole mass in bulges to be smaller than that. However our simulations start at a situation resembling the final stage of a galactic merger, which means we are concentrating on the innermost part of the allready spherical system.

The black hole particles are initially placed symmetrically about the center of mass of the stellar component. Their initial radii are $`r0.64r_\mathrm{h}`$, their initial velocities are 13.6% of the circular velocity at this radius. In the given model units, this represents starting points for the black holes at $`x=\pm 0.5`$ and $`v_y=\pm 0.1`$. The center of mass of the stellar component is at the origin. The mass factor between a stellar particle and a black hole particle is: $`1338.5`$ for 131072 particles, $`669.7`$ for 65536 particles, and $`335.4`$ for 32768 particles.

In order to have a statistical basis for analysis, we compare the results from five runs with 32768 particles, two runs with 65536 particles, and three runs with 131072 particles. Not all runs reached the $`60`$ time unit mark due to time step scheduling problems caused by accuracy problems in very close encounters between stars and a black hole particle. The regularization methods implemented in EuroStar are identical to the ones suitable for open or globular cluster simulations. The extreme situation in the late stages of the sinking binary black hole problem may cause the chain algorithm to fail (Mikkola and Aarseth, 1993). This problem can be solved by applying different regularization methods. However, up to the point of failure, the simulations conserved the total energy with relative errors below $`10^4`$. In all runs, the binary becomes bound at approximately $`10`$ time units.

The parameters of the hybrid code have been adjusted in the following way: The SCF part uses the parameter $`\alpha =0.5`$ for the basis functions. With this choice, the basis functions represent a Plummer model to zeroth order, which is in accordance with the models used by Clutton-Brock (1973). Also, this choice ensures an optimal representation of the actual potential by the expansion method. To allow flexibility in the expansion, seven basis functions are used for the radial direction and five ($`l=[0\mathrm{}5]`$, $`m=[5\mathrm{}5]`$) for the angular expansions. The NBODY6++ part uses $`\eta _i=0.01`$ for the irregular time steps and $`\eta _r=0.02`$ for the regular time steps. Furthermore, the Ahmad-Cohen neighbor scheme (Ahmad and Cohen, 1973) has been modified in such a way that the search radius for neighbors is enhanced by a factor of $`7.7`$, $`14.4`$, and $`27.7`$ for interactions with the black holes in the runs with 32768, 65536, and 131072 particles respectively.

### 4.2 The motion of the massive bodies

In order to compare the runs, all data have been binned by the parameter $`t`$, which represents the integrated time of the system in $`N`$-body time units (Heggie and Mathieu, 1986). Table 1 gives the number of sample points for the orbital data of the black hole particles. For technical reasons, the runs with 131072 particles could not be continued to 60 time units. When binning the total particle number $`N_{\mathrm{tot}}`$, table 1 shows the number of samples in the row labeled “total”. In the following, we present the results for the motion of the black hole binary within the stellar system.

#### 4.2.1 Sinking rate of the binary

Figure 3 shows the evolution of the quantity $`1/a`$ as the average from the runs above. Equation (10) allows us to compute the semi-major axis $`a`$ from the orbital data of the runs. The data is binned for comparison according to the prescription above. From the data in each bin we evaluate the average $`1/a`$ and the standard deviation. In order to find the hardening rate, we fit a line to the averages, plotted as circles in Figure 3. The standard deviations of the data points, given in the plot as the error-bars, supply the weighting factors.

For the average over all runs plotted in Figure 3, the regression line has a slope of $`8.7\pm 0.4`$. The dependency on the particle number can be deduced from the data shown on the left side. The slopes are $`9.6\pm 0.5`$ for the 32768 particle simulation, $`8.6\pm 0.2`$ for the 65536 particle simulation, and $`6.8\pm 0.2`$ for the 131072 particle simulation. There is clearly a dependence of the results for the sinking rate on the particle number. Compared with other quantities we analyze in this study, the noise level in the data for $`1/a`$ is low. We observe strong interactions between the stellar and the black hole particles in runs with 32768 particles. For this reason, $`1/a`$ shows strong steplike changes in both directions. This is most likely due to the small particle number.

#### 4.2.2 Evolution of the eccentricity

We are only studying the evolution of the eccentricity after the binary became bound. Because after 20 time units the eccentricity evolves relatively smoothly for each run, we are concentrating our analysis on the time range between 20 and 60 time units.

Figure 4 shows the mean eccentricity binned in time slots with a width of five time units. The symbols represent the averages in these bins. 32768 data are plotted with open triangles, 65536 data with open squares, and 131072 data with circles. While the eccentricities settle at values between $`ϵ=0.6`$ and $`ϵ=0.9`$ for the runs with 32768 particles, the runs with higher particle numbers show a fairly parallel evolution. The averages of the 32768 particle runs agree very much with the averages from 131072 data. However, the averages for the 65536 data are clearly higher.

With our initial conditions, the binary evolves in a highly eccentric orbit, which is around $`ϵ=0.85`$. The system retains this high eccentricity until the end of our simulations.

#### 4.2.3 Evolution of the angular momentum

In order to study the evolution of the angular momentum of the bound binary we plot the angle $`\theta `$ between the $`z`$-coordinate axis and $`𝐥/l`$ versus time units in Figure 5. $`\theta `$ is zero initially. As with $`1/a`$, all data from the simulations are binned and averaged. The open circles in Figure 5 represent the averages, while the error bars are the standard deviations in the data.

Once the binary becomes bound, the $`\theta `$ changes only slightly in all simulations. Averaged over the time between the first bound orbit of the massive particles and the end of the simulations, the average value of $`\theta `$ becomes $`0.5\pm 0.3`$ for the 32768 runs, $`0.3\pm 0.1`$ for the 65536 runs, and $`0.5\pm 0.1`$ for the 131072 runs.

The results in Figure 5 can be fitted by a straight line. The slope of this line is $`0.003\pm 0.003`$. When we group the simulations according to particle number, the fitting lines have slopes of $`0.016\pm 0.004`$ for the runs with 32768 particles, $`0.0044\pm 0.0006`$ for the 65536 runs, and $`0.006\pm 0.002`$ for the 131072 runs. Though small, these slopes are all significantly nonzero and distinct from one other. Torques clearly act on the binary system throughout the simulations. The data for the runs with 32768 particles and with the small mass ratio between black holes and stellar particles is very noisy and shows steplike changes in $`\theta `$.

While $`\theta `$ evolves in an ordered way until the binary becomes bound, the angle $`\varphi `$ between the $`x`$-coordinate axis and the normalized angular momentum vector behaves more randomly. Until the massive particles become bound, $`\varphi `$ changes rapidly reaching all values between $`0`$ and $`2\pi `$. However, once the binary becomes bound, $`\varphi `$ settles to a single value for each run. The changes in $`\varphi `$ are subsequently of the same order of magnitude as for $`\theta `$.

#### 4.2.4 Wandering motion of the binary

Studying the wandering of the binary using the quantity $`r_{\mathrm{com}}^2`$, we can compare the observed motion to the expected Brownian motion in the system. $`r_{\mathrm{com}}`$ is the distance from the center of mass of the black hole binary to the center of mass of the stellar system. Figure 6 implies that the mean motion is not constant with time. However, since the slope of the fitting line is $`(1.0\pm 1.1)\times 10^5`$, the behavior is constant within $`1\sigma `$ uncertainty. For the individual particle number groups the situation is as follows: For 32768 particles we find a slope of $`(0.6\pm 1.5)\times 10^5`$, for 65536 particles a slope of $`(1.0\pm 0.7)\times 10^5`$, and for 131072 particles a slope of $`(0.9\pm 7.8)\times 10^6`$. Compared to its mean value over the whole simulation, the evolution of $`r_{\mathrm{com}}^2`$ with time introduces changes of not more than 10%. For this reason, we assume $`r_{\mathrm{com}}^2`$ to be constant for the analysis of the Brownian motion. Figure 7 shows the mean squared distance between the center of mass of the black hole system and the stellar system as a function of the total particle number of the simulations. The slope of the fitting line is $`(4.5\pm 5.6)\times 10^9`$. Given our small sample of runs we cannot determine a dependency of $`r_{\mathrm{com}}^2`$ on the particle number.

Figures 9 and 8 show the evolution of $`v_{\mathrm{com}}^2`$ as a function of time and total particle number $`N_{\mathrm{tot}}`$. The quantity $`v_{\mathrm{com}}`$ is the velocity of the center of mass of the black holes relative to the velocity of the center of mass of the stellar system. For the time dependence of $`v_{\mathrm{com}}^2`$, we find a slope of $`(7.6\pm 4.9)\times 10^5`$ for the fitting line in Figure 8. For differing total particle numbers this slope is $`(3.5\pm 7.0)\times 10^5`$ for 32768 particles, $`(2.6\pm 2.4)\times 10^5`$ for 65536 particles, and $`(4.2\pm 4.0)\times 10^5`$ for 131072 particles. The slope for the dependence of $`v_{\mathrm{com}}^2`$ on particle number in Figure 9 is $`(1.7\pm 3.9)\times 10^8`$.

#### 4.2.5 Connection between the wandering and the orbital decay

Figure 10 shows the ratio of the wandering and the semi-major axis of the binary orbit $`a^2`$ as a function of time. The evolution of this ratio has a strong dependence on the particle number, as the wandering is dependent on the simulation size. However, all simulations show the same trend in that wandering becomes more important with time for the binary. As the right plot in Figure 10 shows, a fitting line with a slope of $`1.7\pm 0.6`$ can fit the data. However, the data suggests a nonlinear behavior which should be roughly quadratic, since $`1/a`$ increases linearly and $`r_{\mathrm{com}}^2`$ is roughly constant.

#### 4.2.6 The effect of dynamical friction

To study the influence of dynamical friction on the decay of the binary orbit we analyze the behavior of its orbital angular momentum as a function of time. As Figure 11 shows, the decay shows a two mode evolution. Between 0 and 20 time units, linear regression for $`\mathrm{lg}(l)`$ gives a slope of $`(6.8\pm 0.9)\times 10^2`$. The line with the more shallow slope $`(1.3\pm 0.2)\times 10^2`$ represents the behavior between 20 and 60 time units.

### 4.3 Reaction of the stellar system

The stellar system reacts to the motion of the black hole in a generic fashion. We find that our statistical basis is too small for finding a clear dependency of the results on the total number of particles in the simulations. Hence we present only the averages from all of our runs. Figures 12 and 13 show the evolution of the density and the velocity dispersion respectively for particles within a radius of $`r_{\mathrm{csp}}=0.032`$ averaged over all runs as a function of time. While the black hole binary becomes bound at $`10`$ time units, the density has a maximum at $`18`$ time units, and the velocity dispersion is highest at $`23`$ time units. A linear fit ($`y=a+bx`$) has been applied to the evolution of $`\rho `$ and $`\sigma `$ between 20 and 60 time units as plotted in Figures 12 and 13. For $`\rho `$, we find $`a=2.6\pm 0.7`$ and $`b=0.028\pm 0.016`$, for $`\sigma `$ we find $`a=2.8\pm 0.6`$ and $`b=0.030\pm 0.013`$.

## 5 Discussion

### 5.1 Hardening rate

Following Hills (1992), and Quinlan (1996) the hardening rate $`H`$ of a massive binary floating in a sea of light stars is given by,

$$\frac{d}{dt}\frac{1}{a}=H\frac{G\rho }{\sigma }.$$
(15)
With $`G=1`$ and the assumption that the averages for $`\rho `$ and $`\sigma `$ evolve in the same way between 20 and 60 time units, which would render the ratio between $`\rho `$ and $`\sigma `$ constant, we find $`H=8.7\pm 0.4`$. This is significantly smaller than the values given by Hills (1992) ($`H=13.5`$) and Quinlan (1996) ($`H18`$).

Our smaller hardening rate compared to the results of Quinlan (1996) is caused by the lower central density and the core type radial density profile of our Plummer model. Quinlan (1996) uses Jaffe models for his simulations which allow rapid transfer of orbital energy into the dense cusp through tidal interactions. This is also represented in the destruction of the cusp Quinlan (1996) observes, while our simulations show a much weaker change for the central density.

Hills (1992) models the shrinking of the binary orbit through three body encounters. His greater value of $`H`$ is consistent with our simulations. We observe steplike changes of the binding energy at later times of the simulations, which is less pronounced with increasing particle numbers. Because the granularity of the potential is higher in low $`N_{\mathrm{tot}}`$ runs, three body interactions with the black hole binary become more likely. As shown in Figure 2 such three body encounters can enhance the decay of the orbit. Thus, our small $`H`$ indicates that in our simulations shrinking of the black hole orbits is mainly caused by dynamical friction and not so much by tidal destruction of cusps or three body encounters.

### 5.2 Brownian motion

If the black holes reach equipartition of kinetic energy with the stars, their expected mean square velocity follows from,

$$v_{\mathrm{equ}}^2=\frac{m_{}}{m_{\mathrm{com}}}v_{}^2.$$
(16)
$`v_{\mathrm{equ}}^2`$ is the mean square velocity we expect for a particle with mass $`m_{\mathrm{com}}`$, which is the combined mass of the two black holes. $`m_{}`$ is the mass of the stars and $`v_{}^2`$ their mean square velocity. Since our setup involves a Plummer model, we are expecting the binary to move in a harmonic potential in later stages of a simulation and for $`r_{\mathrm{com}}^2v_{\mathrm{com}}^2`$.

While the individual black holes do not reach equipartition, equation (16) can describe the Brownian motion of the system. The sum of the black hole masses is represented by $`m_{\mathrm{com}}`$ and the center of mass motion of the binary by $`v_{\mathrm{com}}`$. The mass ratios between the individual stars and the black hole binary are $`1.49\times 10^3`$, $`7.47\times 10^4`$, and $`3.74\times 10^4`$ in the runs with 32768, 65536, and 131072 particles, respectively. However, equation (16) does not describe the behavior of the center-of-mass motion correctly. While the velocity dispersion drops after 20 time units, the mean square velocity of the black holes $`v_{}^2`$ increases. For this reason, we compare the measured average $`v_{\mathrm{com}}^2`$ for the datasets with the mean expectation from the right side of equation (16). For $`v_{}^2`$, we take the average over all simulations, which is $`1.5`$, neglecting the variability over time. Using this, we can estimate $`v_{\mathrm{equ}}^2`$ and compare it with the measured $`v_{\mathrm{com}}^2`$,

| $`N_{\mathrm{tot}}`$ | $`v_{\mathrm{com}}^2`$ | $`v_{\mathrm{equ}}^2`$ | $`v_{\mathrm{com}}^2/v_{\mathrm{equ}}^2`$ |
| --- | --- | --- | --- |
| 32768 | 0.0072 | 0.0022 | 3.2 |
| 65536 | 0.0026 | 0.0011 | 2.3 |
| 131072 | 0.0022 | 0.0006 | 3.9 |

This means we find a center of mass motion for the binary which exceeds the expected value from Brownian motion. However, the motion is enhanced by larger factors than proposed by Merritt (2000b) or by Chatterjee et al. (2002). A more detailed discussion of this result remains for future work.

### 5.3 Dynamical friction

Following Begelman et al. (1980), dynamical friction becomes inefficient as the driving forces behind the binary black hole orbit decay after it becomes hard. In order to put constraints on this, we estimate the influence of dynamical friction on the decay of the binary in the simulation. From Binney and Tremaine (1987), we take the following expression, which is derived in Chandrasekhar (1943).

$$\frac{d𝐯_{}}{dt}=16\pi ^2G^2\mathrm{log}(\mathrm{\Lambda })M_{}(M_{}+M_{})\frac{\underset{0}{\overset{v_{}}{}}v_{}^2f(v_{},𝐫)𝑑v_{}}{v_{}^3}𝐯_{}.$$
(17)
Equation (17) describes the deceleration of a particle with mass $`M_{}`$ under the influence of weak encounters with surrounding particles having a uniform mass $`M_{}`$. In the special case of a Plummer model, the integral in equation (17) can be expressed in terms of the escape velocity $`v_{\mathrm{esc}}`$ (Aarseth et al., 1974),

$$\underset{0}{\overset{v_{}}{}}v_{}^2f(v_{},𝐫)𝑑v_{}=\frac{n(𝐫)}{C}\underset{0}{\overset{q_{}}{}}q^2(1q^2)^{\frac{7}{2}}𝑑q,$$
(18)
where,
$$C=\underset{0}{\overset{1}{}}q^2(1q^2)^{\frac{7}{2}}𝑑q.$$
(19)
The quantity $`n(𝐫)`$ defines the number density of the stars at the position $`𝐫`$ and $`q=v/v_{\mathrm{esc}}`$. Taking the limit of a continuous system with $`M_{}M_{}`$, the term $`n(𝐫)M_{}(M_{}+M_{})`$ becomes $`M_{}\rho (𝐫)`$. The integral over $`q`$ in equation (18) can then be solved in a closed form.

If the motion of the black holes is determined by their self interaction plus a frictional force term, this friction can be linked to the decay of the angular momentum $`l`$ as follows:

$$𝐚=𝐚_r+a_{\mathrm{df}}\frac{𝐯}{v},$$
(20)
$$\dot{l}=\frac{ma_{\mathrm{df}}}{v}(𝐫\times 𝐯)=\frac{la_{\mathrm{df}}}{v},$$
(21)
where $`𝐚_r`$ is the radial acceleration of the two body motion of the black holes, $`a_{\mathrm{df}}`$ is the dynamical friction acting on each black hole, and $`𝐯`$ is the two-body velocity of the black holes. With equations (17) and (18) we can evaluate the impact of dynamical friction on the orbital angular momentum according to

$$\frac{\dot{l}}{l}=\frac{16\pi ^2G^2}{Cv_{}^3}\mathrm{log}(\mathrm{\Lambda })M_{}\rho (𝐫)\underset{0}{\overset{q_{}}{}}q^2(1q^2)^{\frac{7}{2}}𝑑q.$$
(22)
The gravitational constant $`G`$ is unity in our model units, and the black holes have mass $`M_{}=0.01`$ each. We use the mean orbital velocities for $`v_{}`$. In order to evaluate $`\rho (r)`$, we use the mean separation of the black holes for $`r`$, which introduces only a very small error in a Plummer model. Assuming a linear behavior for $`\mathrm{lg}(l)=a+bt`$, we find $`\dot{l}/l=b\mathrm{ln}(10)`$. Using this to estimate the angular momentum from the slopes $`b`$ of the linear fits in Figure 11, we find $`\mathrm{log}(\mathrm{\Lambda })0.15`$. This result shows that the usual assumption of large $`\mathrm{\Lambda }`$ does not hold.

Both the possibility of a linear fit for the evolution of $`l`$ and the small hardening rate $`H`$ indicate that mainly dynamical friction causes the shrinking of black hole orbits in our simulations.

## 6 Physical units

As stated before, the collisional simulations which include black hole particles do not reach the observed mass contrast in galactic nuclei. In order to transform simulation units to physical units, a system size in parsec or a stellar mass in units of solar mass has to be chosen. Setting the gravitational constant $`G`$, all remaining units can be rescaled (Heggie and Mathieu, 1986).

In the following, a run with 65536 particles is scaled to a physical stellar system. Since this work focuses on the dynamics of galactic nuclei, the physical mass of the supermassive objects motivates the following choices,

$$M_{}1.00015\times 10^7\mathrm{M}_{},$$
(23)
$$M_{}1.49536\times 10^4\mathrm{M}_{}.$$
(24)
The masses are chosen so that the total mass of the system is $`M_{\mathrm{tot}}=10^9\mathrm{M}_{}`$ and the mean mass of a particle is $`\overline{M}=1.52588\times 10^4\mathrm{M}_{}`$. This choice means that every stellar particle with mass $`M_{}`$ represents a compact star cluster with the order of $`10^4`$ particles. The chosen mass for the black hole particle has approximately the same mass as the central black hole of M31 (Magorrian et al., 1998).

The conversion between physical units and $`N`$-body units follows $`x_{\mathrm{phys}}=X_{\mathrm{conv}}x_{\mathrm{sim}}`$ for simulated quantities. By choosing the central velocity dispersion to be $`110\mathrm{km}/\mathrm{s}`$, we find $`T_{\mathrm{conv}}`$ and $`R_{\mathrm{conv}}`$,

$$R_{\mathrm{conv}}=355.39\mathrm{pc},$$
(25)
$$M_{\mathrm{conv}}=10^9\mathrm{M}_{},$$
(26)
$$T_{\mathrm{conv}}=3.1590\times 10^6\mathrm{y},$$
(27)
$$V_{\mathrm{conv}}=110\mathrm{km}/\mathrm{s}.$$
(28)
In a Plummer model, the half mass radius $`r_h`$ is related to the Plummer radius $`R`$ by $`r_h=1.30R`$ (Spitzer, 1987). Since $`R=3\pi /16`$, the half mass radius in model units is $`r_h=0.766`$. Therefore, the initial model for the simulated decay of a black hole in the galactic center is a Plummer sphere with a half mass radius of 272.23 pc. The initial distance between the black holes is 355.39 pc. They become bound after approximately 40 million years. The total simulated time is approximately 190 million years. At the end of the simulation the black hole distances vary from 1 pc at apocenter to 0.2 pc at pericenter. The semi-major axis of the first bound orbit is 21 pc.

This scaling allows us to compare our results with Begelman et al. (1980). We find that our smallest average orbits at the end of the simulation are not yet small enough that gravitational radiation, according to their estimates, would dominate the evolution time scale.

However, at the end of our simulations evolution is still dominated by dynamical friction and not by long evolution time scales for hard binaries as proposed by Begelman et al. (1980). Their estimate for the gravitational radiation shrinking time scale assumes circular orbits for the binary. With eccentricities of roughly 0.85 for the black hole binaries in our runs, we expect gravitational radiation to be efficient and coalescence in roughly $`10^8`$ years after our simulations stopped.

## 7 Conclusions

We have created a new $`N`$-body hybrid code by merging a high accuracy direct Hermite integrator of the standard type (Aarseth, 1999; Spurzem and Kugel, 1999) with a collisionless $`N`$-body method which approximates the potential of a given particle distribution by a series expansion (Hernquist and Ostriker, 1992; Zhao, 1996). The SCF method has been completely rewritten to include a computation of the time derivative of the gravitational force and a fourth order Hermite integrator. We have used this code to model a galactic nucleus containing two massive black holes with up to 128k single particles. The evolution of the binary black hole is followed from an initial phase, to a phase driven by standard dynamical friction where the binary is bound, and then further hardened by three-body encounters with single stars. In that hardening phase, we take full advantage of the regularized three-body integration developed by Mikkola and Aarseth (1996) and Mikkola and Aarseth (1998). The method proves to work well, and reproduces standard expectations, such as the Chandrasekhar dynamical friction in the initial phase. In the final hardening phase due to three-body encounters, we find that the eccentricity of the black hole binary maintains a fairly large value (around 0.85). This is very interesting because it decreases the time scale for gravitational radiation merger of binary black holes dramatically, thus increasing our chances of detecting gravitational radiation from such events with LISA. Due to computational limitations, however, our particle numbers are still not large enough to fully describe the real physical situation. Any further scaling is problematic, and so further work with improved hardware and software must be done.

We study in detail the motion of a black hole binary in the center of a galaxy. We find that the wandering motion does not decay with increasing particle number as expected. The mechanism exciting these anomalous motions is unclear. If they exist in simulations with realistic particle numbers, they will solve the problem of feeding the black holes with fresh stellar dynamical material raised by Gould and Rix (2000).

## 8 Acknowledgements

The authors would like to thank S. Aarseth, D. Heggie, W. Sweatman, C. Theis, C. Boily, D. Merritt, M. Milosavljević, F. Cruz, H. Baumgardt, G. Hensler, L. Hernquist, H. S. Zhao, P. Ghavamian and E. Barnes for fruitful help and discussion. This project is funded by *Deutsche Forschungsgemeinschaft* (DFG) project Sp 345/9-1,2 and Sonderforschungsbereich (SFB) 439 funded at the University of Heidelberg, NSF grant AST 00-71099, NASA grants NAG 5-7019, NAG 5-6037, and NAG 5-9046. Technical help and computer resources are provided by *NIC* in Jülich, *HLRS* in Stuttgart, *TRACS* and *EPCC* in Edinburgh, *ZIB* in Berlin, *SSC* in Karlsruhe, University of Heidelberg, Rutgers University, and University of Kiel, and by the Pittsburgh Supercomputer Center and the San Diego Supercomputer Center. The authors thank the Aspen Center for Physics, the Institute of Astronomy and the Lorenz center at Leiden University for hospitality. The sources for EuroStar are available from the authors or via `http://www.physics.rutgers.edu/~marchems/`

## Appendix A Recurrence relations for Ultraspherical polynomials

Throughout the computation for the forces and force derivatives in our SCF-scheme several special functions have to be tabulated. Recurrence relations provide a very efficient means of calculating these functions. The following recursion relations have been applied to compute ultraspherical polynomials and their derivatives. As starting values for $`n0,1`$ the recurrence formulae for the Gegenbauer or ultraspherical polynomials obey the relation

$$C_n^{(\alpha )}(\xi )=\{\begin{array}{cc}1\hfill & \text{if}n=0,\hfill \\ 2\alpha \xi \hfill & \text{if}n=1.\hfill \end{array}$$
(A1)
The expressions for higher values are given by:

$$C_{n+1}^{(\alpha )}(\xi )=\frac{2(n+\alpha )\xi C_n^{(\alpha )}(\xi )(n+2\alpha 1)C_{n1}^{(\alpha )}(\xi )}{(n+1)}$$
(A2)
From that the first derivative can be computed as: (Abramowitz and Stegun (1972), equation (22.7.22) and table 22.7.)

$$C_{n1}^{(\alpha +1)}(\xi )=\frac{(n+2\alpha )\xi C_n^{(\alpha )}(\xi )(n+1)C_{n+1}^{(\alpha )}(\xi )}{2\alpha (1\xi ^2)}$$
(A3)
For practical reasons and higher accuracy the second derivative polynomial is computed using equation (A2):

$$C_{n+1}^{(\alpha +2)}(\xi )=\frac{2(n+\alpha +2)\xi C_n^{(\alpha +2)}(\xi )(n+2(\alpha +2)1)C_{n1}^{(\alpha +2)}(\xi )}{(n+1)}$$
(A4)
Because the particle track is approximated by using the Hermite scheme, one has to find forces and the first force derivative simultaneously. An approximation using two timesteps for the first force derivative introduces errors to the second and third derivative of the forces. All particles move within a time dependent potential; therefore, the first derivative has a term describing the change of the potential and a term describing the change of force depending on the particle’s orbit.

$$\frac{d}{dt}𝐚(t,𝐫)=\frac{𝐚(t,𝐫)}{t}+\frac{𝐫}{t}\frac{𝐚(t,𝐫)}{𝐫}.$$
(A5)
With the help of the orbit integration for the single particle in a given static potential case, equation (A5) evaluates to:

$$\begin{array}{cc}\hfill \frac{d}{dt}𝐚(t,𝐫)& =(\frac{a_r}{t}+\frac{a_r}{r}\dot{r}+\frac{a_r}{\vartheta }\dot{\vartheta }+\frac{a_r}{\phi }\dot{\phi }a_\vartheta \dot{\vartheta }a_\phi \dot{\phi }\mathrm{sin}\vartheta )𝐞_r\hfill \\ & +(\frac{a_\vartheta }{t}+\frac{a_\vartheta }{r}\dot{r}+\frac{a_\vartheta }{\vartheta }\dot{\vartheta }+\frac{a_\vartheta }{\phi }\dot{\phi }+a_r\dot{\vartheta }a_\phi \dot{\phi }\mathrm{cos}\vartheta )𝐞_\vartheta \hfill \\ & +(\frac{a_\phi }{t}+\frac{a_\phi }{r}\dot{r}+\frac{a_\phi }{\vartheta }\dot{\vartheta }+\frac{a_\phi }{\phi }\dot{\phi }+a_r\dot{\phi }\mathrm{sin}\vartheta +a_\vartheta \dot{\phi }\mathrm{cos}\vartheta )𝐞_\phi \hfill \end{array}$$
(A6)
The evalutation of the first term on the right hand side of equation (A5) is given in section B. The derivatives with respect to the spatial coordinates in equation (A6) can be found in section C.

## Appendix B Time-dependency of the potential

Because all positions and velocities of the dataset are time-dependent, the partial derivatives with respect to $`t`$ apply only to the coefficients $`A_{nlm}`$. These are implemented as the variables $`C_{lm}(r)`$, $`D_{lm}(r)`$, $`E_{lm}(r)`$, and $`F_{lm}(r)`$, from which the partial derivative can be formed:

$$\frac{C_{lm}(r)}{t}=N_{lm}\underset{n=0}{\overset{\mathrm{}}{}}\stackrel{~}{A}_{nl}\stackrel{~}{\mathrm{\Phi }}_{nl}(r)\underset{k}{}m_k\frac{}{t}\left(\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)P_{lm}(\mathrm{cos}(\vartheta _k))\mathrm{cos}(m\phi _k)\right),$$
(B1)
$$\frac{D_{lm}(r)}{t}=N_{lm}\underset{n=0}{\overset{\mathrm{}}{}}\stackrel{~}{A}_{nl}\stackrel{~}{\mathrm{\Phi }}_{nl}(r)\underset{k}{}m_k\frac{}{t}\left(\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)P_{lm}(\mathrm{cos}(\vartheta _k))\mathrm{sin}(m\phi _k)\right),$$
(B2)
$$\frac{E_{lm}(r)}{t}=N_{lm}\underset{n=0}{\overset{\mathrm{}}{}}\stackrel{~}{A}_{nl}\frac{d\stackrel{~}{\mathrm{\Phi }}_{nl}(r)}{dr}\underset{k}{}m_k\frac{}{t}\left(\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)P_{lm}(\mathrm{cos}(\vartheta _k))\mathrm{cos}(m\phi _k)\right),$$
(B3)
$$\frac{F_{lm}(r)}{t}=N_{lm}\underset{n=0}{\overset{\mathrm{}}{}}\stackrel{~}{A}_{nl}\frac{d\stackrel{~}{\mathrm{\Phi }}_{nl}(r)}{dr}\underset{k}{}m_k\frac{}{t}\left(\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)P_{lm}(\mathrm{cos}(\vartheta _k))\mathrm{sin}(m\phi _k)\right).$$
(B4)
With:
$$\begin{array}{cc}\hfill \frac{}{t}\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)& =\frac{r_k}{t}\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)\hfill \\ & \left[\frac{l}{r_k}\frac{r_k^{\frac{1}{\alpha }}}{r_k}\frac{2l+1}{1+r_k^{\frac{1}{\alpha }}}+\frac{4r_k^{\frac{1}{\alpha }}}{r_k}\frac{\alpha (2l+1)+\frac{1}{2}}{\alpha (1+r_k^{\frac{1}{\alpha }})^2}\frac{C_{n1}^{(\omega +1)}(\xi _k)}{C_n^{(\omega )}(\xi _k)}\right],\hfill \end{array}$$
(B5)
$$\frac{}{t}P_{lm}(\mathrm{cos}(\vartheta _k))=\frac{\vartheta _k}{t}\mathrm{sin}(\vartheta _k)\frac{P_{lm}(\mathrm{cos}(\vartheta _k))}{\mathrm{cos}(\vartheta _k)},$$
(B6)
$$\frac{}{t}\mathrm{cos}(m\phi _k)=m\frac{\phi _k}{t}\mathrm{sin}(m\phi _k),$$
(B7)
$$\frac{}{t}\mathrm{sin}(m\phi _k)=m\frac{\phi _k}{t}\mathrm{cos}(m\phi _k).$$
(B8)
in the coefficient computation section the standard leap frog integrator provided by Hernquist and Ostriker (1992) is extended by two additional variables, which are computed by using the recursion relations in section A.

## Appendix C Orbit dependency of the force derivative

In order to account for the change of force due to the particle orbit one has to calculate the nine partial derivatives in equation (A6). These nine derivatives will now be listed. In order to save some space the second derivative of $`\stackrel{~}{\mathrm{\Phi }}_{nl}(r)`$ is given first:

$$\begin{array}{cc}\hfill \frac{^2}{r^2}\stackrel{~}{\mathrm{\Phi }}_{nl}(r)& =\stackrel{~}{\mathrm{\Phi }}_{nl}(r)[(\frac{l}{r}\frac{r^{\frac{1}{\alpha }}}{r}\frac{(2l+1)}{(1+r^{\frac{1}{\alpha }})})^2\frac{l}{r^2}\hfill \\ & +\frac{(2l+1)}{(1+r^{\frac{1}{\alpha }})^2}\frac{r^{\frac{1}{\alpha }}}{\alpha r^2}(\alpha 1+\alpha r^{\frac{1}{\alpha }})\hfill \\ & +[\mathrm{\hspace{0.17em}8}\left(\frac{r^{\frac{1}{\alpha }}}{r}\frac{\omega }{\alpha (1+r^{\frac{1}{\alpha }})^2}\right)(\frac{l}{r}\frac{r^{\frac{1}{\alpha }}}{r}\frac{(2l+1)}{(1+r^{\frac{1}{\alpha }})})\hfill \\ & +\frac{4\omega r^{\frac{1}{\alpha }}}{r^2\alpha ^2(1+r^{\frac{1}{\alpha }})^3}(1\alpha (\alpha +1)r^{\frac{1}{\alpha }})]\frac{C_{n1}^{(\omega +1)}(\xi )}{C_n^{(\omega )}(\xi )}\hfill \\ & +16\left(\frac{r^{\frac{1}{\alpha }}}{\alpha r}\right)^2\frac{\omega (\omega +1)}{(1+r^{\frac{1}{\alpha }})^4}\frac{C_{n2}^{(\omega +2)}(\xi )}{C_n^{(\omega )}(\xi )}\hfill \end{array}$$
(C1)
The nine derivatives can be implemented as follows:

### C.1 Derivatives with respect to $`r`$

The radial derivative for the radial acceleration becomes:

$`{\displaystyle \frac{a_r}{r}}`$ $`={\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}P_{lm}(\mathrm{cos}(\vartheta ))\left[G_{lm}(r)\mathrm{cos}(m\phi )+H_{lm}(r)\mathrm{sin}(m\phi )\right],`$ (C2)
with:
$`G_{lm}`$ $`=N_{lm}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\stackrel{~}{A}_{nl}{\displaystyle \frac{^2}{r^2}}\stackrel{~}{\mathrm{\Phi }}_{nl}(r){\displaystyle \underset{k}{}}m_k\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)P_{lm}(\mathrm{cos}(\vartheta _k))\mathrm{cos}(m\phi _k),`$ (C3)
$`H_{lm}`$ $`=N_{lm}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\stackrel{~}{A}_{nl}{\displaystyle \frac{^2}{r^2}}\stackrel{~}{\mathrm{\Phi }}_{nl}(r){\displaystyle \underset{k}{}}m_k\stackrel{~}{\mathrm{\Phi }}_{nl}(r_k)P_{lm}(\mathrm{cos}(\vartheta _k))\mathrm{cos}(m\phi _k).`$ (C4)
The radial deriavtive for the acceleration in $`\vartheta `$ direction becomes:

$$\begin{array}{cc}\hfill \frac{a_\vartheta }{r}& =\mathrm{sin}(\vartheta )\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\frac{P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{cos}(\vartheta )}\hfill \\ & \times [(\frac{1}{r^2}C_{lm}(r)\frac{1}{r}E_{lm}(r))\mathrm{cos}(m\phi )\hfill \\ & +(\frac{1}{r^2}D_{lm}(r)\frac{1}{r}F_{lm}(r))\mathrm{sin}(m\phi )].\hfill \end{array}$$
(C5)
The radial deriavtive for the acceleration in $`\phi `$ direction becomes:

$$\begin{array}{cc}\hfill \frac{a_\phi }{r}& =\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\frac{mP_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{sin}(\vartheta )}\hfill \\ & \times [(\frac{1}{r^2}D_{lm}(r)\frac{1}{r}F_{lm}(r))\mathrm{cos}(m\phi )\hfill \\ & (\frac{1}{r^2}C_{lm}(r)\frac{1}{r}E_{lm}(r))\mathrm{sin}(m\phi )].\hfill \end{array}$$
(C6)
### C.2 Derivatives with respect to $`\vartheta `$

The derivative with respect to $`\vartheta `$ for the radial acceleration becomes:

$$\frac{a_r}{\vartheta }=\mathrm{sin}(\vartheta )\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\frac{P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{cos}(\theta )}\left[E_{lm}(r)\mathrm{cos}(m\phi )+F_{lm}(r)\mathrm{sin}(m\phi )\right].$$
(C7)
The derivative with respect to $`\vartheta `$ for the acceleration in $`\vartheta `$ direction becomes:

$$\begin{array}{cc}\hfill \frac{a_\vartheta }{\vartheta }& =\frac{1}{r}\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\left(\mathrm{cos}(\vartheta )\frac{P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{cos}(\vartheta )}\mathrm{sin}^2(\vartheta )\frac{^2P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{cos}(\vartheta )^2}\right)\hfill \\ & \times [C_{lm}(r)\mathrm{cos}(m\phi )+D_{lm}(r)\mathrm{sin}(m\phi )].\hfill \end{array}$$
(C8)
The derivative with respect to $`\vartheta `$ for the acceleration in $`\phi `$ direction becomes:

$$\begin{array}{cc}\hfill \frac{a_\phi }{\vartheta }& =\frac{1}{r}\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}m\left(\frac{P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{cos}(\vartheta )}+\mathrm{cos}(\vartheta )\frac{P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{sin}^2(\vartheta )}\right)\hfill \\ & \times [D_{lm}(r)\mathrm{cos}(m\phi )C_{lm}(r)\mathrm{sin}(m\phi )].\hfill \end{array}$$
(C9)
### C.3 Derivatives with respect to $`\phi `$

The derivative with respect to $`\phi `$ for the acceleration in radial direction becomes:

$$\frac{a_r}{\phi }=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}mP_{lm}(\mathrm{cos}(\vartheta ))[F_{lm}(r)\mathrm{cos}(m\phi )E_{lm}(r)\mathrm{sin}(m\phi )].$$
(C10)
The derivative with respect to $`\phi `$ for the acceleration in $`\vartheta `$ direction becomes:

$$\frac{a_\vartheta }{\phi }=\frac{\mathrm{sin}(\vartheta )}{r}\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}m\frac{P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{cos}(\vartheta )}[D_{lm}(r)\mathrm{cos}(m\phi )C_{lm}(r)\mathrm{sin}(m\phi )].$$
(C11)
The derivative with respect to $`\phi `$ for the acceleration in $`\phi `$ direction becomes:

$$\frac{a_\phi }{\phi }=\frac{1}{r}\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}m^2\frac{P_{lm}(\mathrm{cos}(\vartheta ))}{\mathrm{sin}(\vartheta )}[C_{lm}(r)\mathrm{cos}(m\phi )+D_{lm}(r)\mathrm{sin}(m\phi )].$$
(C12)