On the Formation of an Eccentric Nuclear Disk following the Gravitational Recoil Kick of a Supermassive Black Hole

The double nucleus of the Andromeda galaxy is well-explained by a massive disk of apse-aligned stellar orbits about a supermassive black hole (Tremaine 1995), a configuration we call an eccentric nuclear disk (END). Madigan et al. (2018) showed that ENDs are stabilized by inter-orbit, gravitational torques that limit differential apsidal precession and Gruzinov et al. (2020) demonstrated that they are a thermodynamic equilibrium configuration for rotating nuclear stellar clusters. These eccentric, lopsided stellar disks may turn out to be fairly common in galactic nuclei.

Tidal disruption events (TDEs) occur when stars are unbound by a black hole's tidal gravity (Hills 1975). TDEs are observed in post-starburst galaxies at higher rates than other galaxy types (Arcavi et al. 2014; French et al. 2016; Stone & van Velzen 2016) and these starbursts are likely linked to galaxy mergers (Barnes & Hernquist 1991; Cales & Brotherton 2015; Hammerstein et al. 2021). Stellar populations in ENDs can produce TDE rates 10³–10⁴ times that of a spherical, isotropic configuration (Madigan et al. 2018). We propose a formation mechanism for ENDs, which links galaxy mergers to enhanced TDE rates.

Binary supermassive black holes can form as a result of galaxy mergers (Begelman et al. 1980). Gravitational waves are emitted during their inspiral and eventual coalescence (Einstein 1916). In an unequal-mass binary system, the less massive black hole has a higher orbital velocity and is thus more effective at forward-beaming its gravitational wave radiation (Wiseman 1992). This results in anisotropic gravitational waves, which carry away linear momentum and the remnant receives a recoil kick (Bekenstein 1973). The same effect is produced for spinning black hole binaries (Herrmann et al. 2007). In rare cases, the kick velocity may exceed the escape velocity of the galaxy leading to black hole ejection (Gualandris & Merritt 2008), but even for smaller kick velocities, there are important dynamical consequences (Blecha & Loeb 2008). Here we show that ENDs can form as a result of a recoil kick, explaining the enhanced TDE rates we observe in post-merger galaxies.

In Section 2, we introduce a toy model to analytically explain apsidal alignment that may result from a black hole recoil kick. In Section 3, we confirm this formation mechanism using N-body simulations and discuss the orbital properties of the resultant ENDs. We summarize our findings and conclude in Section 4.

Here we set up a toy problem that allows us to predict the orbital structure of a disk of stars following a kick to the central black hole. We consider the simple case of a circular disk formed from a self-gravitating accretion disk (Paczynski 1978; Goodman 2003) and focus on an in-plane recoil kick motivated by Bogdanović et al. (2007; but see also Lodato & Gerosa 2013). To measure apsidal alignment we use the eccentricity vector,

$$\begin{eqnarray}&&\vec{e}=\displaystyle \frac{\vec{v}\times \vec{j}}{{{GM}}_{\mathrm{BH}}}-\displaystyle \frac{\vec{r}}{r},\end{eqnarray} \tag{ 1 }$$

where $\vec{r}$ is the radius vector, $\vec{v}$ is the velocity vector, $\vec{j}$ is the angular momentum vector, and M_BH is the mass of the central black hole. The eccentricity vector points from the apoapsis to the periapsis of an orbit and its magnitude is equal to the scalar eccentricity.

We consider a circular orbit in the x–y plane with positive angular momentum (∣j∣ = j_z ) centered on the black hole located at the origin. Pre-kick, the star's initial position and velocity are given by

$$\begin{eqnarray}&&\vec{r}=r\cos ({\rm{\Omega }}t)\ \hat{x}+r\sin ({\rm{\Omega }}t)\ \hat{y},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\vec{v}=-{v}_{\mathrm{circ}}\sin ({\rm{\Omega }}t)\ \hat{x}+{v}_{\mathrm{circ}}\cos ({\rm{\Omega }}t)\ \hat{y},\end{eqnarray} \tag{ 3 }$$

where ${\rm{\Omega }}\equiv {\left(\tfrac{{{GM}}_{\mathrm{BH}}}{{r}^{3}}\right)}^{1/2}$ is the angular velocity, and v_circ ≡ Ωr is the circular velocity at this radius. Without loss of generality, we give the black hole a kick in the positive x-direction with magnitude v_BH at t = t_kick. We define the initial anomaly of the star at this instant as θ ≡ Ωt_kick. Immediately following the kick, the position of the star remains unchanged but in the reference frame of the black hole, its velocity becomes

$$\begin{eqnarray}&&\vec{v}=(-{v}_{\mathrm{circ}}\sin (\theta )-{v}_{\mathrm{BH}})\ \hat{x}+{v}_{\mathrm{circ}}\cos (\theta )\ \hat{y}.\end{eqnarray} \tag{ 4 }$$

Setting the star's pre-kick eccentricity vector to zero, and using trigonometric identities, the post-kick eccentricity vector becomes

$$\begin{eqnarray}\begin{array}{rcl}\vec{e} & = & \displaystyle \frac{1}{{{GM}}_{\mathrm{BH}}}\left(\displaystyle \frac{1}{2}{{rv}}_{\mathrm{circ}}{v}_{\mathrm{BH}}\sin (2\theta )\ \hat{x}\ +({{rv}}_{\mathrm{BH}}^{2}\sin (\theta )\right.\\ & & \left.+{{rv}}_{\mathrm{circ}}{v}_{\mathrm{BH}}{\sin }^{2}(\theta )+{{rv}}_{\mathrm{circ}}{v}_{\mathrm{BH}})\ \hat{y}\right).\end{array}\end{eqnarray} \tag{ 5 }$$

To approximate a population of stars in an initially circular disk, we assume that the initial anomalies are evenly distributed in [0, 2π) radians. The average eccentricity vector then becomes

$$\begin{eqnarray}\begin{array}{rcl}\langle \vec{e}\rangle & = & \displaystyle \frac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }\vec{e}(\theta )\ d\theta \\ & = & \displaystyle \frac{3}{2}\displaystyle \frac{{v}_{\mathrm{BH}}}{{v}_{\mathrm{circ}}}\ \hat{y}.\end{array}\end{eqnarray} \tag{ 6 }$$

Figure 1(a) shows a schematic diagram representing two stars in this toy problem. We show the initial circular orbit in black with the kick direction indicated by the red arrow labeled v_BH. At the time of the kick, the top (purple) star's instantaneous circular velocity is pointed in the negative x-direction. In the reference frame of the black hole the star takes on a larger velocity as indicated by the top (purple) arrow, and moves onto a larger, eccentric orbit shown by the large (purple) ellipse. The bottom (green) star's instantaneous circular velocity is in the positive x-direction in the same direction as the kick. Thus, in the frame of the black hole the star slows down, takes on a smaller velocity as shown by the bottom (green) arrow, and moves onto a smaller, eccentric orbit indicated by the small (green) ellipse. Both ellipses have eccentricity vectors that point toward the positive y-direction, in agreement with Equation (6). This is analogous to the Hohmann transfer orbit, an orbital maneuver that is often used to transfer a spacecraft between two circular orbits of different radii (Palmore 1984). Black hole kicks induce apsidal alignment in nearby circular stellar orbits and, perhaps counterintuitively, the eccentric disk forms in such a way that the net eccentricity vector direction is orthogonal to the kick direction. We confirm our analytic approach with N-body simulations in the next section.

Zoom In Zoom Out Reset image size

3.1. Numerical Setup

We use the open-source, N-body simulation package REBOUND with the IAS15 integrator, a high accuracy nonsymplectic integrator with adaptive time-stepping (Rein & Liu 2012; Rein & Spiegel 2015). We use code units of G = 1, M_BH = 1, and the inner edge of the disk a_in = 1 such that the period of a circular orbit at the inner edge of the disk is P(a_in) = 2π.

We initialize stars in an axisymmetric, thin disk with N=400 stars in the range a = 1–2 with surface density profile Σ ∝ a⁻², eccentricities uniformly distributed in e = 0–0.1, inclination Rayleigh distributed with scale parameter σ = 3°, longitude of periapsis (ϖ) and mean anomaly uniformly distributed in [0–2π), and disk mass M_disk = 10⁻² M_BH. At t_kick = 10, we kick the black hole in the plane of the disk with varying magnitudes of v_BH = (0.1–1.0) v_circ(a_out) in increments of 0.1 v_circ(a_out), where v_circ(a_out) is the circular velocity of stars at the outer edge of the disk (a_out = 2). Scaling our code units loosely to the M31 system (Lauer et al. 1993; Bender et al. 2005) with M_BH = 10⁸ M_⊙ and a_in = 1 pc means t = 1 translates to∼1490 yr. The simulation in Figure 2 then runs until ∼75 Myr with a black hole kick of v_BH = 0.4 v_circ(a_out) ≈ 186 kms⁻¹.

Zoom In Zoom Out Reset image size

3.2. Formation of an END

In Figure 1(b), we show the stellar orbits projected into the x-y plane before and after the recoil kick for a simulation with v_BH = 0.3 v_circ(a_out). The black hole is indicated by the red cross and the kick direction with a blue arrow. Immediately following the kick, the orbits become more eccentric and the eccentricity vectors point in the positive y-direction as predicted.

We quantify the net apsidal alignment using the average unit eccentricity vector,

$$\begin{eqnarray}&&\langle \hat{e}\rangle =\displaystyle \frac{{\sum }_{i=1}^{{N}_{\mathrm{bound}}}{\hat{e}}_{i}}{{N}_{\mathrm{bound}}},\end{eqnarray} \tag{ 7 }$$

where N_bound is the total number of bound stellar particles, and ${\hat{e}}_{i}$ denotes the unit eccentricity vector of the ith particle. In Figure 2, we show the time evolution of the magnitude of the average unit eccentricity vector of the disk stars for a simulation with v_BH = 0.4 v_circ(a_out). The black hole receives a kick at t_kick = 10. The noise floor, calculated as $1/\sqrt{{N}_{\mathrm{bound}}}$, is shown in the (orange) shaded region. The average unit eccentricity vector lies below the noise floor before the kick. After the kick, it jumps to a magnitude of close to one indicating near perfect apsidal alignment. It then undergoes a coherent, large-amplitude oscillation. This behavior is well-explained by the alignment/misalignment of stars separated in semimajor axis space, as explained below.

In Figure 3, we show a polar plot of the eccentricity vector components (e_x , e_y ) after the kick, with the semimajor axis indicated with the color bar. This figure explains the origin of the oscillation in apsidal alignment seen in Figure 2. Immediately following the kick (a), stellar orbits become eccentric and aligned in the positive y-direction. Stars with small semimajor axes precess much faster than those with larger semimajor axes, and misalignment (b) and realignment (c) between these groups occur periodically.

Zoom In Zoom Out Reset image size

3.3. Eccentricity Gradient

Specifying v_BH, r, and its corresponding v_circ, the orbital eccentricity of a star post-kick depends only on the anomaly of the star at the time of the kick, θ, which we see from taking the magnitude of the eccentricity vector in Equation (5). The semimajor axis of the post-kick orbit can be determined from the vis-viva equation,

$$\begin{eqnarray}&&{v}^{2}={{GM}}_{\mathrm{BH}}\left(\displaystyle \frac{2}{r}-\displaystyle \frac{1}{a}\right).\end{eqnarray} \tag{ 8 }$$

Substituting Equation (4) we obtain,

$$\begin{eqnarray}&&\begin{array}{l}a(\theta )={{GM}}_{\mathrm{BH}}r\cdot \left[2{{GM}}_{\mathrm{BH}}\right.-\ \left({\left(-{v}_{\mathrm{circ}}\sin (\theta )-{v}_{\mathrm{BH}}\right)}^{2}\right)\\ \quad \qquad {\left.+\left({\left({v}_{\mathrm{circ}}\cos (\theta )\right)}^{2}\right)r\right]}^{-1}.\end{array}\end{eqnarray} \tag{ 9 }$$

In Figure 4(a), we show the post-kick eccentricities as a function of semimajor axes in an N-body simulation with v_BH = 0.4 v_circ(a_out). In solid (blue) and dashed–dotted (red) lines, we show the analytic expectation for circular orbits that begin at the inner edge (r = 1) and the outer edge (r = 2) of the disk, respectively. The eccentricity profile has a distinctive "tick-mark" shape. There is a minimum at the radius corresponding to the initial semimajor axis of the stellar orbits and eccentricity values increase both inward or outward, such that the eccentricity gradient is negative for stars with smaller post-kick semimajor axes and positive for stars with larger post-kick semimajor axes. While the orbits diffuse due to two-body relaxation, this pattern is maintained in our simulations for at least t = 50,000.

Zoom In Zoom Out Reset image size

3.4. Spiral in Mean Anomaly

We observe a clear pattern in the distribution of mean anomalies of the stars in the disk following the kick. We can see why a pattern might emerge from Figure 1(a). The top (purple) star, which takes on a large semimajor axis, is positioned such that it is located at the periapsis of its new orbit. On the other hand, the bottom (green) star, which takes on a small semimajor axis, is positioned at the apoapsis of its new orbit. This suggests that there is a pattern in anomalies following the kick that is correlated with the stars' semimajor axes. We calculate true anomaly as,

$$\begin{eqnarray}\nu =\left\{\begin{array}{cl}\arccos \displaystyle \frac{\vec{e}\cdot \vec{r}}{| \vec{e}| | \vec{r}| }, & \mathrm{if}\ \vec{r}\cdot \vec{v}\gt 0\\ 2\pi -\arccos \displaystyle \frac{\vec{e}\cdot \vec{r}}{| \vec{e}| | \vec{r}| }, & \mathrm{if}\ \vec{r}\cdot \vec{v}\lt 0.\end{array}\right.\end{eqnarray} \tag{ 10 }$$

We then convert to mean anomaly, M, by a standard series expansion (Smart 1953). Figure 4(b) shows results of an N-body simulation with v_BH = 0.4 v_circ(a_out) of mean anomalies as a function of semimajor axis immediately post-kick (orange circles) and shortly after at t = t_kick + 10 (purple crosses). Once again, the expected distribution for initially circular orbits at r = 1 and r = 2 are shown by solid (blue) and dashed–dotted (red) lines, respectively. A spiral structure emerges that continues to wind for several tens of orbital periods and dissipates in less than a precession period (≈M_BH/M_disk P).

In this Letter we propose a new formation mechanism for eccentric nuclear disks (ENDs) via the gravitational wave recoil kick of the central black hole. We show that in-plane kicks cause surrounding circular stellar orbits to take on eccentricities orthogonal to the direction of the recoil kick. Post-kick, a transient spiral structure emerges in mean anomaly as a function of semimajor axis. The orbits display a distinctive and long-lived "tick-mark" pattern in their eccentricity as a function of semimajor axis. On secular timescales we show that the average unit eccentricity vector undergoes a coherent oscillation due to the alignment/misalignment of stars separated in semimajor axis space (Figure 2).

In future work we will explore the stability of these aligned disks as a function of kick magnitude. Weak kicks will not result in significant apsidal alignment (Equation (6)), while extreme kicks will leave few stars bound to the black hole. A recoil kick with a significant in-plane component is required to generate an END from a circular disk of stars. In future work we will consider recoil kicks with out-of-plane components, as expected for spinning Kerr black holes (Baker et al. 2008). We will explore other astrophysically relevant stellar configurations including a scoured nuclear core (Begelman et al. 1980) and consider the effects of gas in the system.

Stone & Loeb (2011, 2012) show that recoil kicks temporarily enhance TDE rates up to 0.1 yr⁻¹ gal⁻¹ by instantaneously filling the central black hole's loss cone (see also Komossa & Merritt 2008). Our preliminary results suggest that if an END forms, elevated TDE rates are actually maintained for a long period of time (at least several precession periods). Such high rates could have a significant impact on the lifetime of ENDs formed via recoil kicks. However, much like bars in disk galaxies, which are also composed of aligned stellar orbits, ENDs will evolve dynamically over time. They may dynamically attract new stellar orbits and grow in mass as they age. To this end, we will study the long-term evolution of ENDs formed in this manner.

Finally, this work applies to many other astrophysical scenarios, for example, when the recoil kick is due to anisotropic mass loss from a star or its nonspherical collapse to a stellar mass black hole or neutron star (which leads us to contemplate an END of planets). In the neutron star case, typical kicks are on the order of ∼380 km s⁻¹ (Faucher-Giguère & Kaspi 2006), which means that only the most tightly bound planets will be retained by the neutron star after an in-plane kick. White dwarf kicks, on the other hand, are of much smaller velocity. Kicks of order ∼1 km s⁻¹ (Stone et al. 2015) may lead to END formation beyond the planetary regime at ∼1000 au. In our solar system, this relatively low mass region is populated by scattered disk objects on high eccentricity orbits. This low-kick-velocity, high-orbital-eccentricity configuration is one we will explore in future work.

We thank Angela Collier for the short title suggestion and the anonymous referee for their suggestions that greatly improved the quality of our paper. A.M. gratefully acknowledges support from the David and Lucile Packard Foundation. This work utilized resources from the University of Colorado Boulder Research Computing Group, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University.

Software: REBOUND (Rein & Liu 2012).