# Anisotropic Star Clusters around Recoiling Supermassive Black Holes

Tatsuya Akiba JILA and Department of Astrophysical and Planetary Sciences, CU
Boulder, Boulder, CO 80309, USA<EMAIL_ADDRESS>Ann-Marie Madigan
JILA and Department of Astrophysical and Planetary Sciences, CU Boulder,
Boulder, CO 80309, USA

###### Abstract

Gravitational wave recoil kicks from merging supermassive black hole binaries
can have a profound effect on the surrounding stellar population. In this
work, we study the dynamic and kinematic properties of nuclear star clusters
following a recoil kick. We show that these post-kick structures present
unique signatures that can provide key insight to observational searches for
recoiling supermassive black holes. In Akiba & Madigan (2021), we showed that
an in-plane recoil kick turns a circular disk into a lopsided, eccentric disk
such as the one we observe in the Andromeda nucleus. Building on this work,
here we explore many recoil kick angles as well as initial stellar
configurations. For a circular disk of stars, an in-plane kick causes strong
apsidal alignment with a significant fraction of the disk becoming retrograde
at large radii. If initial orbits are highly eccentric, an in-plane kick forms
a bar-like structure made up of two anti-aligned lopsided disks. An out-of-
plane kick causes clustering in the argument of periapsis, $\omega$,
regardless of the initial eccentricity distribution. Initially isotropic
configurations form anisotropies in the form of a torus of eccentric orbits
oriented perpendicular to the recoil kick. Post-kick surface density and
velocity maps are presented in each case to highlight the distinct, observable
structures of these systems.

Galaxy nuclei - Stellar dynamics - Gravitational waves

††journal: ApJ††software: REBOUND (Rein & Liu, 2012)

##  1 Introduction

Binary supermassive black holes are a natural consequence of galaxy mergers
(Begelman et al., 1980), and gravitational waves are emitted during their
inspiral and eventual coalescence. In an asymmetric black hole binary (i.e.,
different masses and/or spins), anisotropic gravitational waves carry away
linear momentum in addition to angular momentum and energy. Most radiation
occurs near time of merger (González et al., 2007; Lousto & Zlochower, 2008;
Blanchet et al., 2005) and the remnant black hole receives a recoil kick
(Peres, 1962; Bekenstein, 1973; Wiseman, 1992; Campanelli et al., 2007;
Herrmann et al., 2007). Typical kick velocities of $\mathcal{O}(100)~{}\rm{km\
s^{-1}}$ result in damped oscillations of the merger remnant within the
galactic potential with its eventual return to the nucleus on
$\mathcal{O}$($10^{6}$–$10^{9}$ yr) timescales (Merritt et al., 2004; Madau &
Quataert, 2004; Komossa & Merritt, 2008a). In rare cases, the kick velocity
may exceed the escape velocity of the galaxy leading to black hole ejection
(Gualandris & Merritt, 2008). Detections of recoiling supermassive black holes
would provide key constraints on black hole binary mass and spin evolution, on
gravitational wave event rates, and the merger evolution of galaxies.

Observational searches for recoiling black holes have thus far focused on
active galaxies. If a recoiling black hole is accreting at the time of kick,
it carries with it its accretion disk and broad line region. The recoiling
remnant may then be observed as an offset active galactic nucleus (AGN) (Madau
& Quataert, 2004; Blecha & Loeb, 2008). Numerous such candidates have been
discovered either spectroscopically via Doppler-shifted optical broad lines
and/or via kpc-scale spatial offsets (e.g., Komossa et al., 2008; Chiaberge et
al., 2018; Kim et al., 2018; Hogg et al., 2021). However, most black holes are
not actively accreting as AGN (Kormendy & Ho, 2013) which greatly limits
candidate numbers. It is, furthermore, difficult to rule out alternate
possibilities such as gaseous outflows, AGN jet activity, and dual AGN (e.g.,
Comerford et al., 2017).

In this paper, we explore the effect that the gravitational wave recoil kick
has on the surrounding stellar population, as the basis for using the distinct
structures to detect recoiling systems. Komossa & Merritt (2008b) and Merritt
et al. (2009) first explored this topic finding that a recoiling black hole
carries with it a “hypercompact stellar system” (HCSS). This system may reveal
itself via tidal disruption flares or accretion of gas from stellar winds onto
the black hole. Furthermore, the properties of HCSSs contain information about
the merger histories of galaxies, binary black hole evolution, and the
distribution of gravitational-wave kicks.

In Akiba & Madigan (2021) we explored the post-kick orbits of the stars within
a HCSS. To make the problem analytically tractable, we focused on the simple
case of the stars initially being distributed in a circular disk, inspired by
the theory of circumbinary star formation in a self-gravitating accretion disk
(Paczynski, 1978; Goodman, 2003). We furthermore considered the effect of an
in-plane recoil kick motivated by Bogdanović et al. (2007). We showed that the
kick experienced by the merger remnant can directly result in the formation of
a lopsided stellar disk of eccentric, aligned orbits, a configuration we call
an eccentric nuclear disk (END). Eccentric disks are fairly common in our
local universe (Lauer et al., 2005). The double nucleus of M31 (Lauer et al.,
1993), for instance, is well-explained by an eccentric disk (Tremaine, 1995).
If the eccentric disk was formed via gravitational wave recoil, its dynamic
and kinematic properties can illuminate M31’s merger history.

To measure alignment of stellar orbits we use the eccentricity vector,

$\vec{e}=\frac{\vec{v}\times\vec{j}}{GM_{\bullet}}-\frac{\vec{r}}{r}\ ,$ (1)

where $\vec{r}$ is the radius vector, $\vec{v}$ is the velocity vector,
$\vec{j}$ is the angular momentum vector, and $M_{\bullet}$ 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.
In Akiba & Madigan (2021), we showed that a circular orbit in the $x$-$y$
plane about a central object that is given an in-plane recoil kick in the
$x$-direction will have a preferential direction of eccentricity vectors,

$\langle\vec{e}\rangle=\frac{3}{2}\frac{v_{\rm kick}}{v_{*}}\ \hat{y}\ ,$ (2)

where $v_{\rm kick}$ is the recoil kick speed and $v_{*}$ is the initial speed
of the star. That is, a lopsided disk of stars will form around the recoiling
black hole perpendicular to its direction of motion. We showed using $N$-body
simulations that the eccentric nuclear disk maintains its apsidal-alignment
over tens of thousands of orbits (the limit of our simulations) and that this
leads to a greatly enhanced rate of tidal disruption events.

Here we extend our original study in Akiba & Madigan (2021) to more general
initial stellar distributions and a range of kick directions. Our goal is to
systematically explore these parameter spaces and present a catalog of
asymmetric stellar distributions following a recoil kick imparted on the
central object. We note that while we focus our attention on stars surrounding
a recoiling supermassive black hole, the dynamics explored here are directly
applicable to other astrophysical contexts such as the orbits of planetary
material following a white dwarf natal kick or a supernova kick.

We begin in Section 2 with stars once again in a circular disk experiencing an
in-plane kick, and sequentially add complexity. In Section 3 we add an out-of-
plane component to the kick, and in Section 4 we give non-zero values of
eccentricity to the stars. In Section 5, we relax the constraint of the stars
being in a disk and expand our exploration to spherical distributions. In each
section we confirm analytic predictions with numerical simulations. We
summarize our findings and conclude in Section 6.

##  2 Circular Disk with an In-plane Kick Revisited

###  2.1 Analytics: circular disk, in-plane kick

In Akiba & Madigan (2021), we examined the kick magnitude required to induce
apsidal alignment in a disk with a very narrow radial range. In reality, we
are more interested in the inverse problem. The recoil kick is determined by
relativistic effects, which in theory can be calculated using the black hole
orbits, masses and spins. The surrounding stellar disk, which spans a wide
range of semi-major axes, reacts to this kick with certain radial ranges
responding more strongly than others. In Akiba & Madigan (2021), we
hypothesized that there are three distinct radial regimes of interest based on
comparing the circular speed of stars, $v_{*}$, and the recoil kick speed,
$v_{\rm kick}$. When $v_{*}\gg v_{\rm kick}$, the stellar orbits are
negligibly affected and the $\sim$circular disk is retained. When $v_{*}\sim
v_{\rm kick}$, the apsidal alignment is strong and an eccentric nuclear disk
forms. Finally, when $v_{*}\ll v_{\rm kick}$, stars become unbound from the
black hole. However, due to the narrow radial range investigated, we
overlooked a key feature — the presence of a predominantly retrograde disk at
large radii.

Figure 1: Circular disk with an in-plane recoil kick. A schematic diagram
illustrating four pre- and post-kick orbits. (Left:) Initial circular orbits
A, B, and C with markers indicating the stars’ positions at the time of kick,
and the arrow labeled $v_{\rm kick}$ is the black hole kick direction. The
stars initially move on counter-clockwise (prograde) orbits. (Right:) The
post-kick eccentric orbits where dashed lines indicate where the stars move on
clockwise (retrograde) orbits. Thus, we expect a retrograde population to
emerge above a certain radius in the post-kick stellar distribution.

A schematic diagram illustrating the emergence of a retrograde disk is shown
in Figure 1. In the left panel, three different initial circular orbits A, B,
and C are shown in solid lines with several possible star locations indicated
by markers. The recoil kick is shown by the arrow labeled $v_{\rm kick}$.
Orbit A is the inner region investigated in Akiba & Madigan (2021): the top
(pink triangle) star speeds up in the reference frame of the black hole moving
onto a larger semi-major axis orbit ($\rm{A}^{\prime}_{2}$) while the bottom
(pink diamond) star slows down and goes into a smaller semi-major axis orbit
($\rm{A}^{\prime}_{1}$), as shown in the right panel, with both orbits
remaining prograde. The triangle star will remain bound for kick velocities
$v_{\rm kick}<(\sqrt{2}-1)v_{*}$, which correspond to radii

$r<\frac{(\sqrt{2}-1)^{2}GM_{\bullet}}{v_{\rm kick}^{2}}\ .$

Apsidal alignment will be strongest where $\langle\vec{e}\rangle\sim 1$. Using
Equation 2, we define this point to be the characteristic radius of the
eccentric nuclear disk,

$r_{c}\equiv\frac{4GM_{\bullet}}{9v_{\rm kick}^{2}}\ .$ (3)

The diamond star will have a prograde post-kick orbit if $v_{\rm kick}<v_{*}$,
at radii

$r<\frac{GM_{\bullet}}{v_{\rm kick}^{2}}\ .$

Moving to larger radii, the stellar velocity can reverse and become retrograde
in the reference frame of the black hole, illustrated by orbits B and C in the
schematic. For orbit B, the (orange five-pointed) star moves onto a retrograde
orbit with smaller semi-major axis ($\rm{B}^{\prime}$). This happens when
$v_{*}<v_{\rm kick}<2v_{*}$, or at radii

$\frac{GM_{\bullet}}{v_{\rm kick}^{2}}<r<\frac{4GM_{\bullet}}{v_{\rm
kick}^{2}}\ .$

Even further out at orbit C, the (green six-pointed) star goes into a larger
semi-major axis orbit that is retrograde ($\rm{C}^{\prime}$). This occurs when
$2v_{*}<v_{\rm kick}<(\sqrt{2}+1)v_{*}$, corresponding to radii

$\frac{4GM_{\bullet}}{v_{\rm
kick}^{2}}<r<\frac{(\sqrt{2}+1)^{2}GM_{\bullet}}{v_{\rm kick}^{2}}\ .$

Beyond $r_{\rm{out}}\approx 13\ r_{c}$, no star remains bound to the recoiling
black hole in this set-up.

###  2.2 Numerical Set-Up

To explore the post-kick radial structures further, we use the open-source,
$N$-body simulation package REBOUND (Rein & Liu, 2012). We use code units of
$G=1$, $M_{\bullet}=1$, and characteristic radius $r_{c}=1$ such that the
period of a circular orbit at this radius is $P(r_{c})=2\pi$. For
astrophysical context, we will translate these units to $M_{\bullet}=4\times
10^{6}\ M_{\odot}$ and $r_{c}=0.765$ pc (given by Equation 3) where we choose
$v_{\rm kick}=100$ km/s.

To study the instantaneous post-kick stellar distribution, we initialize
$N=5\times 10^{4}$ stars in an axi-symmetric, thin disk spanning four orders
of magnitude in semi-major axis space, $a=0.005$–$50$, with surface density
profile $\Sigma\propto a^{-1}$. Inclination is Rayleigh distributed with scale
parameter $\sigma_{i}=3^{\circ}$. Longitude of periapsis, $\varpi$, and mean
anomaly, $\mathcal{M}$, are uniformly distributed in $[0,2\pi)$. Since the
characteristic radius is $r_{c}=0.765$ pc, the initial semi-major axis range
translates to $a=0.00382$–$38.2$ pc. The outer edge is beyond the radius of
influence for a black hole of this mass, but we note the sensitive dependence
of the characteristic radius on the kick velocity, $r_{c}\propto v_{\rm
kick}^{-2}$. If the kick velocity doubles, $v_{\rm kick}=200$ km/s, the
characteristic radius drops by a factor of four, $r_{c}=0.19$ pc, and so on.
As we are investigating the distribution of stars instantaneously following
the kick, the stars are massless in our simulations. Surface density and
velocity maps are generated by evolving stars on their post-kick Keplerian
orbits.

###  2.3 Numerical Results: circular disk, in-plane kick

####  2.3.1 Radial Regions and Apsidal Alignment

Figure 2: Circular disk with an in-plane recoil kick. (Top:) The post-kick
distribution of the $z$-component of angular momentum (normalized by the
circular angular momentum at a given radius) as a function of post-kick semi-
major axis scaled by the characteristic radius, with the color bar showing the
$y$-component of the eccentricity vector. We show radii corresponding to the
circular orbits illustrated in the schematic Figure 1, A, B, and C (black
dashed lines), and the corresponding post-kick distribution emanating from
those radii (pink, orange, and green solid lines). (Second from top:) Two
measures of post-kick apsidal alignment: the magnitude of the mean
eccentricity vector (blue solid line), $|\langle\vec{e}\rangle|$, and the mean
unit eccentricity vector (orange dashed line), $|\langle\hat{e}\rangle|$,
binned by semi-major axis. The shaded noise floor shows the maximum expected
deviation from 0 (at a $3\sigma$ level) for an isotropic distribution. (Second
from bottom:) The eccentricity profile as a function of semi-major axis color-
coded by orbital inclination. (Bottom:) The retrograde fraction of the disk as
a function of semi-major axis. The post-kick stellar orbits exhibit a rich
structure of apsidal alignment, eccentricity, and retrograde fraction each
with a sensitive dependence on semi-major axis.

In Figure 2 in the top panel, we show the post-kick distribution of the
$z$-component of angular momentum (normalized by the circular angular momentum
at a given radius), $j_{z}/j_{\rm{circ}}$, as a function of semi-major axis,
with the color bar showing the $y$-component of the eccentricity vector.
Prograde orbits are above ($j_{z}>0$) and retrograde orbits are below
($j_{z}<0$) the dashed black line. Circular orbits are farther from the dashed
line with $|j_{z}/j_{\rm{circ}}|\sim 1$ and orbits become highly eccentric as
$|j_{z}/j_{\rm{circ}}|\rightarrow 0$. Overlaid on top are the three radii
schematically illustrated in Figure 1. All stars that start at an inner orbit
A remain prograde, but some stars become eccentric with a smaller semi-major
axis ($\rm{A}_{1}^{\prime}$) whereas others take on a larger semi-major axis
($\rm{A}_{2}^{\prime}$) depending upon their initial position with respect to
the kick vector. Stars that start at B can become retrograde, once again
depending on their initial position, but stars that flip orientation move to
smaller semi-major axis orbits ($\rm{B}^{\prime}$). The stars that remain on
prograde orbits move to high eccentricity orbits. Stars beginning at C that
remain bound all flip to retrograde orientation and take on a larger semi-
major axis orbit ($\rm{C}^{\prime}$). Based on these results, we identify five
kinematically distinct post-kick radial regions instantaneously following the
kick:

  1. 1.

At $a<0.1$, $v_{*}\gg v_{\rm kick}$ and the circular orbits are mostly
unaffected.

  2. 2.

At $0.1<a<1$, we have a prograde eccentric disk in the $+y$-direction with
eccentricities generally increasing with increasing semi-major axis.

  3. 3.

At $1<a<10$, we have an eccentric disk in the same direction ($+y$) and an
emerging retrograde population. The prograde stars become more eccentric
whereas retrograde stars become more circular with increasing semi-major axis.

  4. 4.

At $10<a\lesssim 3\times 10^{3}$, we get a disk where all of the orbits are
eccentric and the prograde orbits are aligned in the same direction as before
($+y$), but the retrograde orbits are aligned in the opposite direction
($-y$). Both prograde and retrograde orbits become more eccentric with
increasing semi-major axis.

  5. 5.

Finally, at $a\gtrsim 3\times 10^{3}$, we do not observe many stars bound to
the black hole. For most reasonable choices of astrophysical scale, this
region would be outside of the black hole’s radius of influence.

In the second from the top panel of Figure 2, we show two measures of apsidal
alignment binned by semi-major axis. The mean eccentricity vector (solid blue
line) is given by

$|\langle\vec{e}\rangle|=\left|\frac{\sum_{i}^{N}\vec{e}_{i}}{N}\right|\ ,$

where $N$ is the number of bound particles and $\vec{e}_{i}$ is the
eccentricity vector of the $i$-th particle. The mean unit eccentricity vector
(dashed orange line) is similarly defined by replacing $\vec{e}$ with
$\hat{e}$, the unit eccentricity vector. A perfectly aligned eccentric nuclear
disk would have $|\langle\hat{e}\rangle|=1$ and
$|\langle\vec{e}\rangle|=\langle e\rangle$, the mean scalar eccentricity,
whereas a random distribution of eccentricity vectors gives both
$|\langle\hat{e}\rangle|$ and $|\langle\vec{e}\rangle|\sim 0$. The mean unit
eccentricity vector was used in Akiba & Madigan (2021) and measures the
alignment of eccentricity vectors independent of how eccentric the orbits are.
For an eccentric disk to be stable, however, both the magnitude and alignment
of eccentricity vectors are important for orbits to effectively
gravitationally torque each other (Madigan et al., 2018). The difference
between $|\langle\vec{e}\rangle|$ and $|\langle\hat{e}\rangle|$ is clearest at
small semi-major axes. While the eccentricity vectors are extremely aligned
with $|\langle\hat{e}\rangle|\sim 1$, the orbits remain $\sim$circular at
small semi-major axis, so apsidal alignment is unlikely to be maintained long-
term. $|\langle\vec{e}\rangle|$ captures this effect and gives us the expected
peak in apsidal alignment near $r_{c}$. The noise floor shows the maximum
expected $|\langle\vec{e}\rangle|$ or $|\langle\hat{e}\rangle|$ deviations
(within 3$\sigma$) if the stellar distribution were isotropic.
$\sigma=1/\sqrt{N}$ for $|\langle\hat{e}\rangle|$ and $\sigma=\sqrt{\langle
e^{2}\rangle/N}$ for $|\langle\vec{e}\rangle|$. Since
$1/\sqrt{N}>\sqrt{\langle e^{2}\rangle/N}$ for bound stars, we plot
$1/\sqrt{N}$ as the “noise floor”. The eccentricity vectors become most
aligned when $0.1<a<10$ as expected.

In the second from the bottom panel of Figure 2, we show the eccentricity
profile as a function of semi-major axis with the inclination distribution
shown with the color bar. The prograde population has a predominantly positive
eccentricity gradient, $de/da>0$, whereas the retrograde population hits a
minimum eccentricity at $a\sim 10$ with eccentricities increasing both below
and above. Finally, in the bottom panel of Figure 2, we show the retrograde
fraction as a function of semi-major axis. We see that the retrograde
population emerges for semi-major axis larger than $\sim r_{c}$, and the
fraction is $\gtrsim 40\%$ for $a>10$. This significant retrograde fraction
produces distinct density and velocity maps offering a unique opportunity to
identify recoiling black holes using their surrounding star clusters.

####  2.3.2 Post-Kick Density and Velocity Profile

Figure 3: Circular disk with an in-plane recoil kick. Pre- and post-kick
number density as a function of apoapsis distance, $r_{+}$. The initial pre-
kick number density is flat, but the post-kick number density extends to
higher apoapsis distances with a turn around from a flat region to a power-law
region around $r_{+}\approx 4$. At $r_{+}\gtrsim 4$, the number density is
roughly distributed as $N\propto r^{-0.93}$ according to our linear fit.

The stellar density profile changes as a result of the recoil kick. Since we
start with a circular disk of stars, the initial radius of a stellar orbit is
well-defined. Post-kick, however, many of the stellar orbits obtain
significant eccentricities. In this case, stars change their distance from the
central object dramatically over the course of their orbit. Since stars on
eccentric orbits spend most of their time near apoapsis, we assume that
characterizing the number density as a function of the apoapsis distance,
$r_{+}$, will be the most accurate in an orbit-averaged sense. In Figure 3, we
show how the number density as a function of apoapsis distance changes
following the kick. Before the kick, the number density is flat with apoapsis
distance, consistent with $\Sigma\propto a^{-1}$. Post-kick, the stellar
distribution extends to larger apoapsis distances and there appears a turnover
from a flat region to a power-law sloped region beyond an apoapsis distance
$r_{+}\approx 4$. Beyond this radius, the approximate number density profile
is $N\propto r^{-0.93}$, which corresponds to a cuspy surface density profile
$\Sigma\propto r^{-1.93}$. In other words, the number density remains
relatively unchanged at small radii and becomes significantly steeper at
larger radii.

Figure 4: Circular disk with an in-plane recoil kick. The face-on, post-kick
distribution of stars at (top left:) $t=0$, (top right:) $t=0.15$ Myr, and
(bottom:) $t=0.75$ Myr. Bound stars are separated into prograde (blue triangle
markers) and retrograde (orange “+” markers). The unbound stars are marked
with green dots. The “x” marker shows the center of the galaxy and the “o”
marker shows the recoiling black hole. The arrow labeled $v_{\rm kick}$
indicates the recoil kick direction. All positions are shown in the galactic
center frame. Length scale and time units assume a black hole mass of
$M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s.
The prograde, bound orbits form an eccentric nuclear disk in the
$+y$-direction. The unbound stars form an arc of over-density behind the
recoiling black hole. Figure 5: Circular disk with an in-plane recoil kick.
The post-kick surface density distribution of bound stars in the (top:) face-
on and (bottom:) edge-on orientations at (left:) $t=0$ and (right:) $t=0.1$
Myr, zoomed-in to sub-parsec scales. On this scale, all stars are on prograde
orbits. The recoil kick direction is marked in each orientation and labeled as
$v_{\rm kick}$. The “o” marker shows the recoiling black hole, and density
contours are computed in the reference frame of the recoiling black hole.
Length scale and time units assume a black hole mass of $M_{\bullet}=4\times
10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s. The post-kick
density contours reveal an eccentric nuclear disk. Figure 6:  Circular disk
with an in-plane recoil kick. The line-of-sight velocity maps in the frame co-
moving with the recoiling black hole (“o” marker) in the (top:) $x$-$y$ plane,
(center:) $x$-$z$ plane, and (bottom:) $y$-$z$ plane at (left:) $t=0$,
(center:) $t=0.1$ Myr, and (right:) $t=0.35$ Myr. In each orientation, the
color bar is adjusted such that approaching velocities are blue and receding
velocities are red. The recoil kick direction is marked and labeled as $v_{\rm
kick}$ for each row. Length scale and time units assume a black hole mass of
$M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s.
The edge-on orientations clearly show the presence of a retrograde disk at
large radii: there exists a mix of blue-shifted and red-shifted stars on
either side of the black hole.

While this radial number density profile is informative, the post-kick stellar
distribution is non-axisymmetric and so we need to explore the full two-
dimensional profile. We plot the post-kick distribution of stars in the face-
on orientation in Figure 4 at three different times following the kick. We
split the population into bound and unbound stars (green dots), and the bound
stars are further separated into ones on prograde (blue triangles) and
retrograde (orange “+” markers) orbits. The kick is in the $+x$-direction as
indicated by the arrow, the galactic center is marked by the “x” symbol and
the recoiling black hole is shown with the “o” symbol. Positions are computed
in the galactic center frame. To obtain physical units, we assume a black hole
mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ recoiling at $v_{\rm
kick}=100$ km/s. We see the asymmetry immediately emerge at $t=0$: stars are
more likely to remain bound to the black hole if initially located toward the
$-y$-direction, and one can clearly see the predominantly retrograde
population at large radii. As we move forward in time, the prograde population
reveals its apsidal alignment in the $+y$-direction, consistent with Figure 2.
The unbound stars also exhibit structure: at $t=0.15$ Myr, they form a large
spiral which leaves a dense arc of unbound stars behind the recoiling black
hole. We note that the stellar orbits do not interact with each other and
therefore do not precess. We also do not take into account an external
potential for the galaxy.

We investigate the density structure at smaller radii using a kernel density
estimation. The resulting surface density contours are shown in Figure 5 at
two different times following the kick for both the face-on (top) and edge-on
(bottom) orientations. The “o” marker shows the recoiling black hole, and all
contours are computed with respect to the black hole frame of reference.
Notice that we have zoomed-in to sub-parsec scales, so the structure we see
here concerns the bound, prograde population of stars. At $t=0$, the structure
is symmetric since the stars have not yet evolved on their new orbits. At
$t=0.1$ Myr, we can clearly see that the stars have arranged themselves into a
lopsided disk.

In Figure 6, we plot the line-of-sight velocity of the stars in the face-on
($x$-$y$) orientation as well as both edge-on ($x$-$z$ and $y$-$z$)
orientations at three different times following the kick. In each orientation,
the color bar is adjusted such that approaching velocities are blue and
receding velocities are red. All velocities are calculated in the black hole’s
reference frame. In the face-on orientation, no clear structure can be seen.
In the $x$-$z$ plane, the initial symmetry breaks down over time with
different radial ranges exhibiting different velocity signatures. While the
disk is initially entirely prograde with the left side blue-shifted and the
right side red-shifted, the retrograde population emerges over time, and at
$t=0.35$ Myr, the right side is blue-shifted and the left side is red-shifted
at large radii. Closer in, the prograde population dominates the velocity
signal: the left side is blue-shifted and the right side is red-shifted. The
$y$-$z$ plane line-of-sight velocities clearly show the emergence of the
retrograde orbits from the beginning: the stars that remain bound at large
radii are primarily located toward the $-y$-direction and are all red-shifted
since they are on retrograde orbits.

##  3 Introducing Out-of-plane Kicks

###  3.1 Analytics: circular disk, arbitrary kick direction

Now we expand our toy model in Akiba & Madigan (2021) by considering a recoil
kick with an out-of-plane component. For simplicity, we once again assume that
the circular orbit is initially in the $x$-$y$ plane, and without loss of
generality, we assume that the kick has some positive $x$-component as well as
some $z$-component. Say the magnitude of the kick is $v_{\rm kick}$ and it
makes an angle $\alpha$ with the $z$-axis, where $0\leq\alpha\leq\pi$. The
initial circular orbit has position and velocity given by

$\vec{r}=r\cos(\theta)\ \hat{x}+r\sin(\theta)\ \hat{y}\equiv x\ \hat{x}+y\
\hat{y}\ ,$ $\vec{v}=-v_{*}\sin(\theta)\ \hat{x}+v_{*}\cos(\theta)\
\hat{y}\equiv v_{x}\ \hat{x}+v_{y}\ \hat{y}\ .$

Post-kick, the velocity vector in the reference frame of the black hole
becomes

$\vec{v}^{\prime}=(v_{x}-v_{\rm kick}\sin(\alpha))\ \hat{x}+v_{y}\
\hat{y}-v_{\rm kick}\cos(\alpha)\ \hat{z}\ .$ (4)

Using vector identities and subtracting the initial eccentricity vector yields
the post-kick eccentricity vector components

$\displaystyle GM_{\bullet}e_{x}$
$\displaystyle=\frac{1}{2}r\sin(2\theta)v_{*}v_{\rm kick}\sin(\alpha)$
$\displaystyle+r\cos(\theta)v_{\rm kick}^{2}\cos^{2}(\alpha)\ ,$

$\displaystyle GM_{\bullet}e_{y}$ $\displaystyle=r\sin(\theta)v_{\rm
kick}^{2}$ $\displaystyle+r(1+\sin^{2}(\theta))v_{*}v_{\rm kick}\sin(\alpha)\
,$

$\displaystyle GM_{\bullet}e_{z}$ $\displaystyle=-r\cos(\theta)v_{\rm
kick}^{2}\sin(\alpha)\cos(\alpha)\ .$

Additionally, the post-kick angular momentum vector components are

$\displaystyle j_{x}=-r\sin(\theta)v_{\rm kick}\cos(\alpha)\ ,$ $\displaystyle
j_{y}=r\cos(\theta)v_{\rm kick}\cos(\alpha)\ ,$ $\displaystyle
j_{z}=rv_{*}+r\sin(\theta)v_{\rm kick}\sin(\alpha)\ .$

####  3.1.1 Limit: $\alpha=\pi/2$, an in-plane kick

When the kick is completely in the plane of the disk, $\alpha=\pi/2$, the
eccentricity and angular momentum vector components reduce to

$\displaystyle GM_{\bullet}e_{x}=\frac{1}{2}r\sin(2\theta)v_{*}v_{\rm
kick}\sin(\alpha)\ ,$ $\displaystyle GM_{\bullet}e_{y}=r\sin(\theta)v_{\rm
kick}^{2}+r(1+\sin^{2}(\theta))v_{*}v_{\rm kick}\ ,$ $\displaystyle
GM_{\bullet}e_{z}=0\ ,$ $\displaystyle j_{x}=0\ ,$ $\displaystyle j_{y}=0\ ,$
$\displaystyle j_{z}=rv_{*}+r\sin(\theta)v_{\rm kick}\ ,$

so the orbit remains in the $x$-$y$ plane. If we average over
$\theta:0\rightarrow 2\pi$ for the eccentricity vector components, we recover
Equation 2, the analytic result from Akiba & Madigan (2021).

####  3.1.2 Limit: $\alpha=0$ (or $\pi$), an out-of-plane kick

The other limiting case is if the kick is completely out-of-plane, $\alpha=0$
(or $\pi$). In this limit,

$\displaystyle GM_{\bullet}e_{x}=r\cos(\theta)v_{\rm kick}^{2}\ ,$
$\displaystyle GM_{\bullet}e_{y}=r\sin(\theta)v_{\rm kick}^{2}\ ,$
$\displaystyle GM_{\bullet}e_{z}=0\ $ $\displaystyle\implies$ $\displaystyle
GM_{\bullet}\vec{e}=v_{\rm kick}^{2}\ \vec{r}\ \rm{,\ and}$

$\displaystyle j_{x}=\mp r\sin(\theta)v_{\rm kick}\ ,$ $\displaystyle
j_{y}=\pm r\cos(\theta)v_{\rm kick}\ ,$ $\displaystyle j_{z}=rv_{*}$
$\displaystyle\implies$
$\displaystyle\vec{j}=\vec{j}_{\rm{init}}\pm\frac{rv_{\rm kick}}{v_{*}}\
\vec{v}\ ,$

so the circular orbit develops an eccentricity in the direction of the star’s
initial position, and the orbital angular momentum inclines in the direction
(opposite that) of the star’s initial velocity vector. The result is that all
of the orbits collectively roll over their major axes and the stars are at
periapsis of these new elliptical orbits. In addition, we see that the post-
kick eccentricity is simply given by $e=(v_{\rm kick}/v_{*})^{2}$. This means
that stars which initially orbit at radii smaller than

$r_{\rm{out}}=\frac{GM_{\bullet}}{v_{\rm kick}^{2}}$

will remain bound where as stars beyond this radius are all ejected. The
inclination of the orbit monotonically increases with radius as the ratio
$v_{\rm kick}/v_{*}$ increases, and at $r_{\rm{out}}$ the angular momentum
vector is simply given by

$\vec{j}=rv_{*}(\hat{z}+\hat{v})\ ,$

which means that the maximum possible inclination is 45∘.

###  3.2 Numerical Results: circular disk, out-of-plane kick

####  3.2.1 Distribution of Orbital Elements and Tilt Vector Alignment

Figure 7: Circular disk with an out-of-plane recoil kick. (Top:) Distribution
of the post-kick eccentricity as a function of semi-major axis with the
inclination on the color bar. (Bottom:) Two measures of $\omega$-clustering:
the magnitude of the mean tilt vector (blue solid line),
$|\langle\vec{t}\rangle|$, and the mean unit tilt vector (orange dashed line),
$|\langle\hat{t}\rangle|$, where the tilt vector is defined in Equation 5,
binned by semi-major axis. The shaded noise floor region shows the maximum
expected deviation from 0 (within $3\sigma$) if the distribution were
isotropic. Both eccentricity and inclination monotonically increase with semi-
major axis, and tilt vector alignment is statistically significant at most
semi-major axes. Figure 8: Circular disk with an out-of-plane recoil kick. A
histogram of the post-kick mean anomaly, $\mathcal{M}$, and argument of
periapsis, $\omega$. Mean anomaly is strongly clustered around $0$ since every
star begins its new orbit at periapsis, and $\omega$ is strongly clustered
around $\pi$ since the new periapsis corresponds to the orbit’s descending
node for all orbits. Figure 9: Circular disk with an out-of-plane recoil kick.
The edge-on, post-kick distribution of stars at (top left:) $t=0$, (bottom
left:) $t=0.15$ Myr, and (right:) $t=0.75$ Myr. Bound and unbound stars are
shown with blue triangle markers and green dots, respectively. The “x” marker
shows the center of the galaxy and the “o” marker shows the recoiling black
hole. The recoil is applied in the direction indicated by the arrow labeled
$v_{\rm kick}$. All positions are computed in the galactic center frame of
reference. Length scale and time units assume a black hole mass of
$M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s.
The post-kick distribution is a large ball of trailing stars with the bound
stars forming a distinct dome-shape behind the recoiling black hole. Figure
10: Circular disk with an out-of-plane recoil kick. The face-on, post-kick
distribution of stars at (left:) $t=0$, (center:) $t=0.15$ Myr, and (right:)
$t=0.75$ Myr. Bound and unbound stars are shown with blue triangle markers and
green dots, respectively. The “x” marker shows the center of the galaxy and
the “o” marker shows the recoiling black hole. The recoil kick direction is
out of the plane, marked and labeled as $v_{\rm kick}$. Length scale and time
units assume a black hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and
a kick of $v_{\rm kick}=100$ km/s. The distribution is always axi-symmetric,
but an under-density seems to emerge between the central region and the
surrounding ring which both radially expand over time. Figure 11: Circular
disk with an out-of-plane recoil kick. The post-kick surface density
distribution of bound stars in the (top:) face-on and (bottom:) edge-on
orientations at (left:) $t=0$, (center:) $t=0.01$ Myr, and (right:) $t=0.04$
Myr. The “o” marker shows the recoiling black hole, and these density contours
are calculated in the frame co-moving with the black hole. In each
orientation, the recoil kick direction is marked and labeled as $v_{\rm
kick}$. Length scale and time units assume a black hole mass of
$M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s.
The face-on orientation reveals the evolution of the central region and the
surrounding ring over-density. The edge-on orientation shows the development
of the characteristic dome-shape. Figure 12: Circular disk with an out-of-
plane recoil kick. The line-of-sight velocity maps in the frame co-moving with
the recoiling black hole (“o” marker) in the (top:) face-on and (bottom:)
edge-on orientations at (left:) $t=0$, (center:) $t=0.1$ Myr, and (right:)
$t=0.35$ Myr. In each orientation, the color bar is adjusted such that
approaching velocities are blue and receding velocities are red. The recoil
kick direction is marked and labeled as $v_{\rm kick}$ for each row. Length
scale and time units assume a black hole mass of $M_{\bullet}=4\times 10^{6}\
M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s. The strong correlation in
anomalies and arguments of periapsis causes clear, circular patterns in line-
of-sight velocity for the face-on orientation. The edge-on orientation shows
that the dome-shaped structure maintains coherent rotation.

We numerically study the case of a circular disk with the central black hole
receiving a purely out-of-plane kick. The set-up is identical to what is
described in Section 2.2 with the only difference being that the kick is
perpendicular to the disk plane. In the top panel of Figure 7, we show the
post-kick eccentricity distribution as a function of semi-major axis with the
inclination color-coded. Both eccentricity and inclination monotonically
increase with semi-major axis with the inclination plateauing at $\sim
45^{\circ}$ and eccentricity $e\sim 1$, as analytically predicted. In the
bottom panel of Figure 7, we introduce a new parameter — a vector that can be
used to measure clustering in the argument of periapsis, $\omega$, similar to
how the eccentricity vector, $\vec{e}$, is a vector quantity that lets us
measure clustering in the longitude of periapsis, $\varpi$. Note that
$\varpi=\Omega+\omega$ where $\Omega$ is the longitude of ascending node.
$\Omega$ is the angle between the reference direction ($+x$) and the ascending
node, the point at which an orbit crosses the reference ($x$-$y$) plane and
the star is moving upward (from $-z$ to $+z$). $\omega$, the argument of
periapsis, is the angle between the ascending node and the orbit’s periapsis.
Thus, in a disk, $\varpi$ measures the angle between the reference direction
and the periapsis. The new quantity which we call the tilt vector, $\vec{t}$,
has Cartesian components

$\displaystyle t_{1}=\cos(\Omega)e_{x}+\sin(\Omega)e_{y}\ ,$ (5)
$\displaystyle t_{2}=-\sin(\Omega)e_{x}+\cos(\Omega)e_{y}\ ,$ $\displaystyle
t_{3}=e_{z}\ ,$

where $e_{x}$, $e_{y}$, and $e_{z}$ are the $x$, $y$, and $z$-components of
the eccentricity vector, respectively. This is a rotation of the eccentricity
vector by an angle $-\Omega$, which allows us to calculate the position of the
eccentricity vector from the ascending node of the orbit rather than from the
reference direction. The ratio between $t_{1}$ and $t_{3}$ indicates the
degree to which the orbit is tilting over their major axis versus minor axis.
$t_{1}$ dominates when an orbit rolls over its major axis whereas $t_{3}$
dominates when an orbit pitches over its minor axis. It is important to note
that the tilt vector is sensitive to our choice of reference plane and
direction. Equation 5, for instance, gives the tilt vector components when we
choose the $x$-$y$ plane to be our reference plane with $+x$ being our
reference direction.

The tilt vector is closely related to the inclination angles $i_{a}$ and
$i_{b}$ as defined in Madigan & McCourt (2016) which measure the inclination
components that roll and pitch, respectively (see Appendix A). Using the tilt
vector, we can use the same measure of vector alignment we use for
eccentricity vectors. In the bottom panel of Figure 7, we show the mean tilt
vector magnitude (blue solid line),

$|\langle\vec{t}\rangle|\equiv\left|\frac{\sum_{i}^{N}\vec{t}_{i}}{N}\right|\
,$

as well as the similarly-defined mean unit tilt vector magnitude (orange
dashed line) binned by semi-major axis. We notice that there is significant
$\vec{t}$ and $\hat{t}$-alignment at most semi-major axis ranges, since all of
the orbits tilt coherently over their major axes. Every star begins its new
orbit at the descending node, which also corresponds to its new periapsis.
This creates a highly non-uniform post-kick distribution of the mean anomaly,
$\mathcal{M}$, in addition to the argument of periapsis, $\omega$, as shown in
Figure 8. Mean anomaly is clustered near 0 since stars are initially located
at their periapsis, and $\omega$ is clustered near $\pi$ because the ascending
node of the new orbit corresponds to the apoapsis.

####  3.2.2 Post-Kick Density and Velocity Profile

The post-kick stellar distribution exhibits unique density and velocity
profiles for the purely out-of-plane kick case as well. We show the edge-on
distribution of the stars at three different times in Figure 9. The bound
stars are denoted with blue triangles and the unbound stars are marked by
green dots, with the black hole recoiling in the $+z$-direction as indicated
by the arrow. The galactic center is shown with the “x” marker and the
recoiling black hole is shown with the “o” marker. All positions are
calculated in the galactic center frame. Length and time units assume a black
hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of $v_{\rm
kick}=100$ km/s. At $t=0$, the distribution is an axi-symmetric disk where
only stars within $r_{\rm{out}}\sim 2$ pc remain bound to the recoiling black
hole. At subsequent times, we see a large ball of trailing stars: the bound
stars form a “dome” shape behind the black hole and the unbound stars make up
the rest of the spherical distribution following the black hole. For the same
three times following the kick, we show the face-on distribution of stars in
Figure 10. Since the kick is purely out-of-plane, the distribution remains
axi-symmetric about the $z$-axis throughout, but the bound population expands
significantly in radius, and at $t=0.75$ Myr, there appears an interesting
density structure where there is a central over-density as well as a dense
ring of stars surrounding it, with an under-density in between.

We further explore the density structure closer-in, at $\sim$parsec scales,
using a kernel density estimation. We present both face-on (top) and edge-on
orientations (bottom) of these zoomed-in surface density maps in Figure 11 at
three different times following the kick. The “o” marker shows the recoiling
black hole, and the density contours are constructed in the frame co-moving
with the black hole. In the face-on orientation, we see that the central over-
density and the surrounding ring forms immediately following the kick with
both components expanding in radius over time. In the edge-on orientation, the
appearance is a dome-like shape: the stars that are most tightly bound to the
black hole forms a puffy disk at small radii, and the stars at slightly larger
radii form an arc that trails the recoiling black hole. The puffy disk and arc
in the edge-on orientation correspond to the central over-density and the
surrounding ring in the face-on orientation, respectively.

The alignment of the tilt vector (bottom panel of Figure 7) accompanied by the
strong mean anomaly correlation (Figure 8) produces coherent structures in the
line-of-sight velocity of stars. We show velocity maps in the reference frame
co-moving with the black hole in both the face-on and edge-on orientations
over time in Figure 12. In each orientation, the color bar is adjusted such
that approaching velocities are blue and receding velocities are red. In the
face-on orientation, the post-kick stellar distribution exhibits distinct
circular patterns in line-of-sight velocity due to the fact that the stars all
gain eccentricities and tilt over their major axes. With the recoil kick
pointed toward us, the outer edge of the ring is preferentially red-shifted
and the inner regions are preferentially blue-shifted. In the edge-on
orientation, we see that the dome-shaped, puffy disk still maintains the same
angular momentum axis since stars all remain prograde post-kick: the right
side is preferentially red-shifted and the left side is preferentially blue-
shifted. We note that this is very similar to Figure 2 in Lippai et al. (2008)
who show an aerial view of the density profile of a tightly-bound gaseous disk
a week after responding to an out-of-plane recoil kick.

##  4 Adding Initial Eccentricities

###  4.1 Analytics: disk with initial eccentricities, arbitrary kick direction

We now expand our analyses to incorporate initial eccentricities. The initial
position vector of a body on an eccentric orbit is given by

$\vec{r}=\frac{a(1-e^{2})}{1+e\cos(\theta-\varpi)}\ (\cos\theta\
\hat{x}+\sin\theta\ \hat{y})\ ,$

where $a$ is the semi-major axis, $e$ is the eccentricity, $\theta$ is the
true anomaly, and $\varpi$ is the longitude of periapsis. The initial velocity
vector is given by

$\displaystyle\vec{v}=\frac{1}{\sqrt{a(1-e^{2})}}$
$\displaystyle[-(\sin\theta+e\sin\varpi)\ \hat{x}$
$\displaystyle+(\cos\theta+e\cos\varpi)\ \hat{y}]\ .$

We move into the recoiling black hole’s reference frame by shifting the
velocity vector as we have done in Equation 4 before, and calculate the
eccentricity vector. We average over mean anomaly, $\mathcal{M}$, and
$\varpi$, so we obtain the analytic expectation starting with an axi-symmetric
disk of stellar orbits with a certain given initial eccentricity.

####  4.1.1 Limit: $e=0$, a circular orbit

When we set $e=0$, the initial position and velocity vectors above reduce to

$\vec{r}=a\cos\theta\ \hat{x}+a\sin\theta\ \hat{y}\ ,$

$\vec{v}=-\frac{1}{\sqrt{a}}\sin\theta\ \hat{x}+\frac{1}{\sqrt{a}}\cos\theta\
\hat{y}\ ,$

and since $a=r$ and $v_{*}=\sqrt{\frac{GM_{\bullet}}{r}}=\frac{1}{\sqrt{r}}$
for a circular orbit, we recover the analytics from Section 3.1.

####  4.1.2 Limit: $e\rightarrow 1$, a high eccentricity orbit

Stars on eccentric orbits spend a majority of their time near apoapsis. In the
limit of $e\rightarrow 1$, we can capture the disk’s behavior by replacing
$v_{*}$ in Equation 2 with $v_{\rm{apo}}$, the star’s speed at apoapsis. This
is given by $v_{\rm{apo}}=\sqrt{\frac{(1-e)GM_{\bullet}}{(1+e)a}}$.
Substituting this and setting $|\langle\vec{e}\rangle|\sim 1$ yields a
characteristic semi-major axis,

$a_{c}=\frac{1-e}{1+e}\ r_{c}\ .$ (6)

This effectively reduces the semi-major axis at which the strongest apsidal
alignment occurs. For an initial eccentricity $e=0.99$, for instance,
$v_{\rm{kick}}=100$ km/s gives $a_{c}=4\times 10^{-3}$ pc as opposed to
$r_{c}=0.8$ pc for a circular disk, a factor of $\sim$200 difference.

##### In-plane kick

Re-investigating kinematically distinct radial regions as in Section 2.3.1 but
for an initially high eccentricity orbit, we see a significant change due to
the fact that the initial apoapsis speed,
$v_{\rm{apo}}=\sqrt{\frac{(1-e)GM_{\bullet}}{(1+e)a}}\rightarrow 0$ is much
slower than the circular velocity at that radius,
$v_{\rm{circ}}=\sqrt{\frac{GM_{\bullet}}{(1+e)a}}\rightarrow\sqrt{\frac{GM_{\bullet}}{2a}}$.
We find that:

  * •

A star with its apoapsis located along the $y$-axis will remain bound if
$v_{\rm kick}<\sqrt{2}v_{\rm{circ}}\pm
v_{\rm{apo}}\rightarrow\sqrt{2}v_{\rm{circ}}$ which translates to
$a<\frac{GM_{\bullet}}{v_{\rm kick}^{2}}$.

  * •

A star with its apoapsis initially located toward the $-y$-direction will only
remain prograde for kicks $v_{\rm kick}<v_{\rm{apo}}\rightarrow 0$, so roughly
half of the stars will flip to retrograde at most radii.

  * •

A star along the $y$-axis will have the same apoapsis pre- and post-kick if
$v_{\rm kick}<v_{\rm{circ}}\pm v_{\rm{apo}}\rightarrow v_{\rm{circ}}$ or
$a<\frac{GM_{\bullet}}{2v_{\rm kick}^{2}}$.

  * •

The same star will move onto a larger semi-major axis orbit where the pre-kick
apoapsis becomes the post-kick periapsis if $v_{\rm{circ}}<v_{\rm
kick}<\sqrt{2}v_{\rm{circ}}$ or at $\frac{GM_{\bullet}}{2v_{\rm
kick}^{2}}<a<\frac{GM_{\bullet}}{v_{\rm kick}^{2}}$.

In other words, for initially high eccentricity orbits, the initial speed of
the star is negligible, and their post-kick evolution is dictated by the
recoil kick at most radii. Stars initially located toward the $-y$-direction
will flip to retrograde orientation with eccentricity vectors pointing in the
$+y$-direction at small semi-major axes and $-y$-direction at larger semi-
major axes. Stars initially in the $+y$-direction will have eccentricity
vectors pointing in the $-y$-direction at small semi-major axes and
$+y$-direction at larger semi-major axes. Hence, for a given radial range, we
expect a prograde eccentric nuclear disk and an anti-aligned retrograde
eccentric disk creating a linear, bar-like structure.

##### Out-of-plane kick

It is fairly straightforward to apply the same analytic limit for purely out-
of-plane kicks. As with the circular disk case, the recoil kick will tilt all
of the orbits coherently over their major axes, but since $v_{\rm{circ}}\gg
v_{\rm{apo}}$ for high eccentricity orbits, most stars will begin their new
orbit at the apoapsis which coincides with the descending node. Thus, while we
get $\omega$-clustering near $\pi$ for an initially circular disk, we expect
clustering near 0 when starting with a high eccentricity disk. In addition, at
larger semi-major axes, $v_{\rm kick}\gg v_{\rm{apo}}$, so the recoil kick
should cause significant inclinations of $\sim 90^{\circ}$.

###  4.2 Numerical Results: parameter study of initial eccentricities

Figure 13: Parameter study of an initially thin disk with varying
eccentricities and kick directions. (Top:) The mean $y$-component of the
eccentricity vector, $\langle e_{y}\rangle$, and (bottom:) the mean first
component of the tilt vector, $\langle t_{1}\rangle$, which measures roll, as
a function of the initial eccentricity with different lines showing different
out-of-plane kick components. We show $\langle e_{y}\rangle$ and $\langle
t_{1}\rangle$ because they dominate the mean eccentricity vector and the mean
tilt vector, respectively. The other components, $\langle e_{x}\rangle$,
$\langle e_{z}\rangle$, $\langle t_{2}\rangle$, and $\langle t_{3}\rangle$,
are all $\sim$0\. $\alpha=0^{\circ}$ is a completely out-of-plane recoil kick
whereas $\alpha=90^{\circ}$ is a completely in-plane kick. The solid lines are
analytic expectations and the “x” markers show the numerical results. We omit
the analytic expectation for the purely in-plane tilt vector, since the
longitude of ascending node is not well-defined for a two-dimensional orbit in
the $x$-$y$ plane. Eccentricity vector alignment becomes weaker with
increasing initial eccentricities and larger out-of-plane kick components.
Tilt vector alignment happens in all cases except when the kick is purely in-
plane.

We study the instantaneous alignments of the eccentricity vector and the tilt
vector as we change initial eccentricities and the out-of-plane kick component
with a suite of numerical simulations. Once again, the set-up is the same as
what is described in Section 2.2 except we incrementally change the initial
eccentricity of orbits and the kick direction. For simulations starting with
high eccentricities, we adjust the inner edge of the disk to be $a_{c}/2$
where $a_{c}$ is given by Equation 6 to capture the behavior near the
characteristic semi-major axis. For instance, for $e=0.99$, $a_{c}\approx
1/200\ r_{c}$, so the inner edge of the disk is set to $a_{c}/2\approx
0.0025$. In the top panel of Figure 13, we plot the numerical results and
corresponding analytic calculations of the mean $y$-component of the
eccentricity vector, $\langle e_{y}\rangle$, as a function of the initial
eccentricity of orbits for various kick angles. We plot $\langle e_{y}\rangle$
because it is the component that dominates the mean eccentricity vector:
$\langle e_{x}\rangle$ and $\langle e_{z}\rangle$ are both $\sim$0 for all
simulations. These results are averaged over the semi-major axis range,
$a=0.005$–$50$, in our simulations with the given surface density profile.
Numerical results are shown with the “x” markers whereas the analytic results
are shown in solid lines. $\alpha=90^{\circ}$ is an in-plane kick and
$\alpha=0^{\circ}$ is an out-of-plane kick. The circular, in-plane
($\alpha=90^{\circ}$) case induces the most apsidal alignment, and the
alignment becomes weaker as the out-of-plane kick component and initial
eccentricities increase. Nonetheless, the preferential direction for apsidal
alignment remains the same: if the in-plane recoil kick component is in the
$+x$-direction, the eccentricity vectors will preferentially become aligned in
the $+y$-direction, orthogonal to the in-plane kick. A statistically
significant ($>3\sigma$) alignment is induced in most cases except when the
kick is purely out-of-plane ($\alpha=0^{\circ}$).

In the bottom panel of Figure 13, we show the mean first component of the tilt
vector, $\langle t_{1}\rangle$, as a function of the initial eccentricity for
various kick directions. $\langle t_{1}\rangle$ is the component that
dominates the mean tilt vector which means that all of the stellar orbits are
coherently tilting over their major axes. $\langle t_{2}\rangle$ and $\langle
t_{3}\rangle$ are $\sim$0 for all simulation runs. The solid lines show our
analytic expectations and the “x” markers show our numerical results. We omit
the analytic estimate for the purely in-plane kick case, since the longitude
of ascending node is not well-defined for an orbit in the reference plane. As
we have shown before, the tilt vectors are aligned for a circular, out-of-
plane case, but $\langle t_{1}\rangle<0$ due to the fact that stars begin at
periapsis. Pointing in the direction of the eccentricity vector, each orbit
rolls to the left, so $t_{1}$ is preferentially negative. On the other hand,
there is also a very strong $\vec{t}$-alignment induced in the limit of
initially high eccentricity orbits provided that there is at least some out-
of-plane kick component. This is due to the $\omega$-clustering at 0 predicted
in Section 4.1.2, where orbits roll over their major axes but the stars begin
their new orbits preferentially at apoapsis. Pointing in the direction of the
eccentricity vector, each orbit rolls to the right, so $t_{1}$ is
predominantly positive. It is important to note that while we only capture the
average tilt alignment here, nearly every star begins at periapsis or apoapsis
and rolls over its major axis following the kick. In general, we see two
clusters at $\omega=0$ and $\pi$ in competition, so $\langle t_{x}\rangle\sim
0$ simply means that these two clusters are similar in size. It does not mean
that there is no tilt alignment. $\omega$ tends to cluster at $\pi$ for
circular orbits and at 0 for eccentric orbits, and the turning point when
$\langle t_{x}\rangle$ crosses 0 depends largely on $\alpha$. Statistically
significant ($>3\sigma$) alignment emerges in most cases except when the kick
is completely in-plane ($\alpha=90^{\circ}$).

###  4.3 Numerical Results: disk with initial $e=0.99$, in-plane kick

####  4.3.1 Radial Regions and Apsidal Alignment

Figure 14: Disk with high initial eccentricities ($e=0.99$) with an in-plane
recoil kick. (Top:) The post-kick distribution of the $z$-component of angular
momentum (normalized by the circular angular momentum at a given radius) as a
function of semi-major axis, with the color bar showing the $y$-component of
the eccentricity vector. (Second from top:) The $y$-component of the mean
eccentricity vector, $\langle e_{y}\rangle$, binned by semi-major axis and
separated into prograde ($j_{z}>0$) and retrograde ($j_{z}<0$) stars. The mean
for prograde orbits is shown in the blue solid line whereas the mean for
retrograde orbits is shown in the orange dashed line. The shaded noise floor
shows the maximum expected deviation from 0 (at a $3\sigma$ level) if the
distribution were isotropic. (Second from bottom:) The eccentricity profile as
a function of semi-major axis with the inclination color-coded. (Bottom:) The
retrograde fraction of the disk at a given semi-major axis range. Two anti-
aligned eccentric nuclear disks create a linear, bar-like structure with a
$\gtrsim 40\%$ retrograde fraction at most semi-major axis ranges. The
eccentricity gradient, $de/da$, is negative inward of $a\sim 3$ and positive
outward of $a\sim 5$.

We numerically study the case of a thin disk of high eccentricity orbits
receiving an in-plane kick. This set-up is not motivated by an astrophysically
realistic scenario. Rather we are building orbital complexity in a sequential
manner. In the top panel of Figure 14, we show the post-kick distribution of
the $z$-component of angular momentum (normalized by the circular angular
momentum at a given radius), $j_{z}/j_{\rm{circ}}$, as a function of the post-
kick semi-major axis for a simulation starting with $e=0.99$ orbits with the
black hole recoiling in-plane. The color bar shows the $y$-component of the
eccentricity vector. As predicted analytically, the prograde and retrograde
sub-populations are nearly mirror images of each other except at very small
semi-major axes ($a\lesssim 10^{-2}$) when $v_{\rm kick}\lesssim
v_{\rm{apo}}$. At these very small semi-major axes, there exists two anti-
aligned lopsided modes but both are prograde. When $10^{-2}\lesssim a\lesssim
4$, there is a prograde eccentric nuclear disk in the $-y$-direction and a
retrograde eccentric disk in the $+y$-direction. At $a\gtrsim 4$, there is a
prograde eccentric disk in the $+y$-direction and a retrograde eccentric disk
in the $-y$-direction.

Figure 15: Disk with high initial eccentricities ($e=0.99$) with an in-plane
recoil kick. The face-on, post-kick distribution of stars at (top left:)
$t=0$, (top right:) $t=0.15$ Myr, and (bottom:) $t=0.75$ Myr. Bound stars are
separated into prograde (blue triangle markers) and retrograde (orange “+”
markers). The unbound stars are marked with green dots. The “x” marker shows
the center of the galaxy and the “o” marker shows the recoiling black hole.
The arrow labeled $v_{\rm kick}$ indicates the recoil direction. Positions are
computed in the galactic center frame of reference. Length scale and time
units assume a black hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and
a kick of $v_{\rm kick}=100$ km/s. The prograde and retrograde populations are
nearly mirror images of each other about the $x$-axis. Prograde stars are
aligned in one direction and retrograde stars are anti-aligned with the
prograde orbits. Figure 16: Disk with high initial eccentricities ($e=0.99$)
with an in-plane recoil kick. The post-kick surface density distribution of
bound stars in the (top:) face-on and (bottom:) edge-on orientation at (left:)
$t=0$ and (right:) $t=0.5$ Myr, where the stars have been separated into
prograde (blue contours) and retrograde (orange contours). The “o” marker
shows the recoiling black hole, and density contours are calculated in the
frame co-moving with the black hole. In each orientation, the recoil kick
direction is marked and labeled as $v_{\rm kick}$. Length scale and time units
assume a black hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick
of $v_{\rm kick}=100$ km/s. The density contours reveal a prograde eccentric
disk and an anti-aligned retrograde disk forming a bar-like structure. Figure
17: Disk with high initial eccentricities ($e=0.99$) with an in-plane recoil
kick. The line-of-sight velocity maps in the frame co-moving with the
recoiling black hole (“o” marker) in the (top:) $x$-$y$, (center:) $x$-$z$,
and (bottom:) $y$-$z$ plane at (left:) $t=0$, (center:) $t=0.1$ Myr, and
(right:) $t=0.35$ Myr. In each orientation, the color bar is adjusted such
that approaching velocities are blue and receding velocities are red. The
recoil kick direction is marked and labeled as $v_{\rm kick}$ for each row.
Length scale and time units assume a black hole mass of $M_{\bullet}=4\times
10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s. In the $x$-$z$
plane, there appears a cluster of stars with high line-of-sight velocity
immediately behind the recoiling black hole. In the $y$-$z$ plane, different
radial regions become coherently blue-shifted over time as both prograde and
retrograde populations evolve along their orbits.

This can also be seen by plotting the $y$-component of the mean eccentricity
vector, $\langle e_{y}\rangle$, separated into prograde and retrograde orbits
as we have done in the second from the top panel of Figure 14. The mean
$e_{y}$ among prograde stars is shown in the blue solid line whereas the mean
for retrograde stars is shown in the orange dashed line. Statistically
significant apsidal alignment is seen in both anti-aligned eccentric disks
outward of $a\sim 10^{-2}$, and the flip in alignment occurs nearly
simultaneously: around $a\approx 3$ for the prograde population and $a\approx
5$ for the retrograde population (with $r_{c}=1$). We also plot the
eccentricity profile as a function of semi-major axis with the inclination
distribution color-coded in the second from the bottom panel of Figure 14. The
prograde and retrograde populations have very similar eccentricity
distributions, hitting a minimum at around $a\approx 3$ and $a\approx 5$,
respectively, and eccentricities increase both inward and outward. In the
bottom panel of Figure 14, we plot the retrograde fraction in the disk as a
function of semi-major axis. Compared to the retrograde profile we saw in
Figure 2, here the retrograde population emerges at much smaller semi-major
axes ($a\sim 10^{-2}$) and is $\gtrsim 40\%$ for $a>r_{c}$. This makes sense
since $a_{c}\sim 10^{-2}\ r_{c}$ for $e=0.99$, so orbits are affected much
closer in.

####  4.3.2 Post-Kick Density and Velocity Profile

For a disk with initially $e=0.99$ orbits with an in-plane kick to the black
hole, we plot the distribution of stars in the face-on orientation at three
different times post-kick in Figure 15. Length scale and time units assume a
black hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of
$v_{\rm kick}=100$ km/s. Once again the unbound (green dots) and bound
populations are separated and the bound population is further split into
prograde (blue triangles) and retrograde (orange “+” markers) orbits, with the
galactic center marked with the “x” symbol and the recoiling black hole shown
by the “o” symbol. The positions are in the galactic center frame of
reference. At this zoomed-out scale of $\sim 10$ pc, we mainly see stars on
large semi-major axis orbits. At $t=0$, the prograde orbits are located in the
$+y$-direction and the retrograde orbits in the $-y$-direction. The bound
population is split $\sim$equally about the $x$-axis in this radial range.
Over time, the prograde and retrograde populations switch places: at $t=0.15$
Myr and $t=0.75$ Myr, the prograde population is now located toward the
$-y$-direction and the retrograde population toward the $+y$-direction. This
is because at large semi-major axis, the stars begin their new orbits at
periapsis, which is located in the $+y$-direction for prograde orbits and
$-y$-direction for retrograde orbits. As the stars evolve along their
Keplerian orbits, they begin to show clustering near their collective
apoapses, which is in the opposite direction ($-y$-direction for prograde and
$+y$ for retrograde orbits).

We investigate the density structure with a kernel density estimation in
Figure 16, where we plot surface density contours in both the face-on and
edge-on orientations at two different times post-kick, zoomed-in to sub-parsec
scales. We separate the population into prograde (blue) and retrograde
(orange) orbits with the recoiling black hole shown with the “o” marker.
Density contours here are in the black hole frame. We see that in this region
of small semi-major axes, the prograde and retrograde populations are
initially located in the $+y$ and $-y$-directions, respectively, and after
$0.5$ Myr, the clustering is unchanged. This is because at small semi-major
axes, stars begin their new orbit at apoapsis.

Figure 18: Disk with high initial eccentricities ($e=0.99$) with an out-of-
plane recoil kick. (Top:) The post-kick eccentricity distribution as a
function of semi-major axis with the inclination on the color bar. (Bottom:)
The post-kick mean first component of the tilt vector, $\langle t_{1}\rangle$,
which measures roll, as a function of semi-major axis. The shaded noise floor
shows the maximum expected deviation from 0 (at a $3\sigma$ level) if the
distribution were isotropic. The eccentricity gradient, $de/da$, is negative
inward of $a\sim 2$ and positive outward, with most stars taking on a high
inclination $\sim 90^{\circ}$. The tilt alignment is statistically significant
at most semi-major axis ranges. Figure 19: Disk with high initial
eccentricities ($e=0.99$) with an out-of-plane recoil kick. The edge-on, post-
kick distribution of stars at (top left:) $t=0$, (bottom left:) $t=0.15$ Myr,
and (right:) $t=0.75$ Myr. Bound and unbound stars are shown with blue
triangle markers and green dots, respectively. The “x” marker shows the center
of the galaxy and the “o” marker shows the recoiling black hole. The arrow
labeled $v_{\rm kick}$ indicates the recoil kick direction. Positions are
computed in the galactic center frame of reference. Length scale and time
units assume a black hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and
a kick of $v_{\rm kick}=100$ km/s. Both bound and unbound stars form an
elongated ball around the recoiling black hole. Figure 20: Disk with high
initial eccentricities ($e=0.99$) with an out-of-plane recoil kick. The face-
on, post-kick distribution of stars at (left:) $t=0$, (center:) $t=0.15$ Myr,
and (right:) $t=0.75$ Myr. Bound and unbound stars are shown with blue
triangle markers and green dots, respectively. The “x” marker shows the center
of the galaxy and the “o” marker shows the recoiling black hole. The recoil
kick direction is out of the plane, marked and labeled as $v_{\rm kick}$.
Length scale and time units assume a black hole mass of $M_{\bullet}=4\times
10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s. Both bound and
unbound populations seem to radially expand over time. Figure 21: Disk with
high initial eccentricities ($e=0.99$) with an out-of-plane recoil kick. The
post-kick surface density distribution of bound stars in the (top:) face-on
and (bottom:) edge-on orientations at (left:) $t=0$, (center:) $t=0.01$ Myr,
and (right:) $t=0.04$ Myr. The “o” marker shows the recoiling black hole, and
the contours are computed in the frame co-moving with the black hole. In each
orientation, the recoil kick direction is marked and labeled as $v_{\rm
kick}$. Length scale and time units assume a black hole mass of
$M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s.
The edge-on orientation reveals an over-density behind the recoiling black
hole. Figure 22: Disk with high initial eccentricities ($e=0.99$) with an out-
of-plane recoil kick. The line-of-sight velocity maps in the frame co-moving
with the recoiling black hole in the (top:) face-on and (bottom:) edge-on
orientations at (left:) $t=0$, (center:) $t=0.1$ Myr, and (right:) $t=0.35$
Myr. In each orientation, the color bar is adjusted such that approaching
velocities are blue and receding velocities are red. The recoil kick direction
is marked and labeled as $v_{\rm kick}$ for each row. Length scale and time
units assume a black hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and
a kick of $v_{\rm kick}=100$ km/s. The face-on orientation shows circular
patterns where the edges of the distribution becomes more blue-shifted over
time. The edge-on orientation demonstrates that the elongated ball-structure
maintains coherent rotation as it evolves.

Finally, in Figure 17, we plot the line-of-sight velocity maps in the face-on
($x$-$y$) and both edge-on ($x$-$z$ and $y$-$z$) orientations at three
different times post-kick. Velocities are in the frame co-moving with the
recoiling black hole, and the color bar in each orientation is adjusted so
that approaching velocities are blue and receding velocities are red. We do
not see much velocity structure in the face-on orientation. In the $x$-$z$
plane, a cluster of stars with high line-of-sight velocity seems to form right
behind the recoiling black hole after some time ($t=0.35$ Myr). This is due to
high eccentricity, large semi-major axis stars that initially begin at
periapsis. Since their eccentricity vectors are aligned along the $y$-axis and
they have high eccentricity orbits, shortly after the kick they have
velocities almost purely in the line-of-sight with a magnitude close to its
high periapsis speed. This creates the collection of high line-of-sight
velocity stars behind the recoiling black hole. In the $y$-$z$ plane, we
clearly see the symmetry about the $x$-axis caused by the $\sim$equal split
between prograde and retrograde stars. At $t=0$, all stars are preferentially
red-shifted since stars toward the $+y$-direction are on prograde orbits and
stars toward the $-y$-direction are on retrograde orbits. At $t=0.1$ Myr, the
two populations begin to switch places and we see several blue-shifted stars.
Finally, at $t=0.35$ Myr, the edges are now preferentially blue-shifted since
the retrograde population has migrated toward the $+y$-direction and the
prograde population has moved toward the $-y$-direction.

###  4.4 Numerical Results: disk with initial $e=0.99$, out-of-plane kick

####  4.4.1 Distribution of Orbital Elements and Tilt Vector Alignment

In the top panel of Figure 18, we plot the post-kick eccentricity distribution
as a function of the semi-major axis with the inclination color-coded.
Circular orbits are found in a narrow region around $a\approx 2$ (in units
where $r_{c}=1$), and eccentricities grow both inward and outward. Since
$v_{\rm{circ}}\gg v_{\rm{apo}}$, the post-kick inclinations of these orbits
are generally very high, $i\sim 90^{\circ}$: one only finds low inclination
orbits at smaller semi-major axes. In the bottom panel of Figure 18, we plot
the first component of the mean tilt vector as a function of semi-major axis.
The high, positive $\langle t_{1}\rangle$ at small semi-major axes is due to
$\omega$-clustering at 0 since the stars begin at their apoapsis which also
corresponds to their descending node (so periapsis coincides with the
ascending node). On the other hand, the negative $\langle t_{1}\rangle$ at
large semi-major axes is due to $\omega$-clustering at $\pi$, since the
descending node now corresponds to the periapsis.

####  4.4.2 Post-Kick Density and Velocity Profile

We show the edge-on distribution of stars in Figure 19 at three different
times post-kick. Stars are split into bound (blue triangles) and unbound
(green dots) orbits, the galactic center is indicated by a “x” marker, and the
black hole with a “o” marker. The positions are computed in the galactic
center frame of reference. At $t=0$, the bound stars are confined to a compact
disk surrounding the recoiling black hole. As the system evolves, both bound
and unbound stars seem to form an elongated ball surrounding the black hole
that expands in size over time. This elongated appearance is due to the high
eccentricity, high inclination orbits with large semi-major axes we see in
Figure 18. In Figure 20, we show the same distribution of stars in the face-on
orientation. The structure is axi-symmetric at all times, as expected, and we
see the ball of bound and unbound stars expand in radius over time. In Figure
21, we show surface density contours in the face-on and edge-on orientations
zoomed-in to $\sim$parsec-scales at three different times following the kick,
using a kernel density estimation. The “o” marker shows the recoiling black
hole, and these density contours are computed in the frame co-moving with the
black hole. In the face-on orientation, the structure always remains axi-
symmetric in appearance. In the edge-on orientation, a puffy disk of tightly
bound stars surround the recoiling black hole as it lifts above the disk
plane, and the more loosely bound stars seem to cluster behind the black hole.

Finally, we plot the line-of-sight velocity profile in both face-on and edge-
on orientations in Figure 22. Velocities are computed in the frame co-moving
with the black hole, and in each orientation, the color bar is adjusted such
that approaching velocities are blue and receding velocities are red. In the
face-on orientation, the stars are initially all red-shifted since they all
start moving away from us on their new orbits; every star begins their new
orbit at its descending node. At $t=0.1$ Myr, there appears to be a circular
pattern in the line-of-sight velocity once again caused by the coherent roll
of the orbits over their major axes, but this pattern is less pronounced than
the circular case (Figure 12) and seems to disperse at earlier times ($t=0.35$
Myr). In the edge-on orientation, one can see that as the black hole begins to
recoil and the disk puffs up into an elongated ball, the rotation profile of
the initial disk is maintained. While the inclination distribution of the
cluster significantly changes, the right side is always red-shifted and the
left side blue-shifted.

##  5 Adding Initial Inclinations

###  5.1 Numerical Results: parameter study of initial inclinations

While we do not perform further calculations here, we can make predictions
based on our analyses thus far. As we increase the inclinations within the
disk, stellar orbital angular momentum vectors will have a wider range of
angles with respect to the recoil kick. This will reduce the overall alignment
of the eccentricity vector. However, we expect to see non-zero tilt vector
alignment even in a puffy disk. As shown in Figure 13, if there is some out-
of-plane component to the kick, the post-kick orbits have significant
$\omega$-clustering near 0 or $\pi$.

We numerically study the instantaneous alignment of the eccentricity vector
and the tilt vector as we change the initial inclination distribution and the
initial eccentricity distribution. Once again, we only slightly alter the set-
up described in Section 2.2. We vary the scale parameter, $\sigma_{i}$, that
dictates the inclination range and gradually puff-up the disk while
simultaneously changing initial eccentricities incrementally. We additionally
simulate an isotropic distribution with a $N(e)=2e\ de$ eccentricity profile
relevant for Keplerian orbits with isotropic velocities (Binney & Tremaine,
1987) as the most general case. In addition, we vary the initial density
profile. We first use the same number density as before which translates to a
steep volume density profile of $\rho\sim a^{-2}$, but we also run simulations
with the Bahcall-Wolf cusp density profile ($\rho\sim a^{-7/4}$) (Bahcall &
Wolf, 1976) following the simulations presented in Merritt et al. (2009).

Figure 23: Parameter study with disks of varying initial inclination
distributions. (Top:) The average $y$-component of the eccentricity vector,
$\langle e_{y}\rangle$, for an in-plane kick and (bottom:) the average first
component of the tilt vector, $\langle t_{1}\rangle$, for an out-of-plane
kick, as a function of the initial eccentricity with different lines showing
different inclination distributions taken from a suite array of numerical
simulations. Apsidal alignment becomes weaker with increased disk thickness,
but tilt alignment is expected regardless of the inclination distribution.

In the top panel of Figure 23, we show $\langle{e_{y}}\rangle$ as a function
of eccentricity for a variety of inclination distributions for an in-plane
kick. As expected, the apsidal alignment becomes weaker as we increase
$\sigma_{i}$, the scale parameter of the Rayleigh distribution. $\langle
e_{y}\rangle$ seems to show relatively weak apsidal alignment in puffy disks
(e.g. $\sigma_{i}=60^{\circ}$) and no alignment in the isotropic case, but as
we discuss in the next section, the mean eccentricity vector fails to capture
symmetric anisotropies, so the isotropic case warrants further investigation.
In the bottom panel of Figure 23, we show $\langle t_{1}\rangle$ as a function
of eccentricity for a variety of inclination distributions for an out-of-plane
kick. As expected, $\omega$-clustering is seen at various eccentricities
either at 0 or $\pi$ with the turning point having a dependence on the
inclination distribution. This makes sense since changing $\sigma_{i}$ changes
the distribution of kick angles with respect to the orbital angular momenta of
individual stars. Notably, $\langle t_{1}\rangle$ predicts tilt alignment even
in the isotropic case, but we emphasize that this measure is dependent on our
choice of reference plane: specifically, tilt alignment is expected in the
plane perpendicular to the recoil kick.

Figure 24: Parameter study with stellar distributions of varying initial
inclinations. (Top two rows:) The post-kick distribution of the $z$-component
of angular momentum (normalized by the circular angular momentum at a given
radius) as a function of semi-major axis, with the color bar showing the
$y$-component of the eccentricity vector for six different simulation runs.
Rows correspond to (top:) inclination Rayleigh-distributed with
$\sigma_{i}=30^{\circ}$ and (second from top:) isotropically-distributed
orbits. Columns correspond to (left:) circular, (center:) $e=0.99$, and
(right:) eccentricity profile $N(e)=2e\ de$ with a Bahcall-Wolf cusp density
profile ($\rho\sim a^{-1.75}$) (Bahcall & Wolf, 1976). The first two columns
use the same number density as before which translates to a volume density
$\rho\sim a^{-2}$. (Second from bottom:) The $y$-component of the mean
eccentricity vector separated into prograde and retrograde orbits (about the
$z$-axis) as a function of semi-major axis corresponding to the isotropic
simulations in the panels above. (Bottom:) The mean first component of the
tilt vector, $\langle t_{1}\rangle$, binned by semi-major axis where the plane
of reference is the plane perpendicular to the kick ($y$-$z$) and the
reference direction is $+y$. In this rotated coordinated system, $t_{1}$ still
measures roll over the major axis but now with respect to the new plane of
reference ($y$-$z$). For the bottom two rows, the shaded noise floor shows the
maximum expected deviation from 0 (at a $3\sigma$ level) for an isotropic
distribution. All of these simulations develop post-kick anisotropies. The
isotropic simulations form a torus-like structure of eccentric orbits that
preferentially align in the plane perpendicular to the recoil kick.

We compare the post-kick distribution of the $z$-component of angular momentum
(normalized by the circular angular momentum at a given radius) as a function
of semi-major axis for six different simulations in Figure 24. The top row
shows simulations where the inclination is initially Rayleigh-distributed with
$\sigma_{i}=30^{\circ}$, and the second from the top row shows simulations
with initially isotropically-distributed stellar orbits. Columns correspond to
different eccentricity distributions: the left column is initially circular
orbits, the center column is initially $e=0.99$, and the right column is an
eccentricity profile $N(e)=2e\ de$. The first two use the same number density
profile as before which is a volume density of $\rho\sim a^{-2}$. For the
right-most column, we change the density profile to that of a Bahcall-Wolf
cusp ($\rho\sim a^{-7/4}$) (Bahcall & Wolf, 1976). There is no net angular
momentum axis for an isotropic configuration, but we pick the $z$-axis around
which we determine whether the orbit is prograde or retrograde, and the kick
in this case is orthogonal to that axis ($+x$). This is a convenient way to
visualize how anisotropies manifest in different ways as we alter the initial
eccentricity and inclination distributions. Focusing on the top row, which is
a puffy disk with $\sigma_{i}=30^{\circ}$, the initially circular case is
fairly similar to the top panel of Figure 2. Inward of $a\sim 10$, the
majority of stars have a significant, positive $e_{y}$ component so we have a
single eccentric nuclear disk. At $a\gtrsim 10$, there is a retrograde
population with apsidal alignment in the opposite direction. The $e=0.99$ case
is similar to the top panel of Figure 14. Inward of $a\sim 4$, we have a
prograde eccentric nuclear disk in the $-y$-direction and a retrograde
eccentric disk in the $+y$-direction forming a bar-like structure. At
$a\gtrsim 4$, they flip so we have a prograde eccentric disk in the
$+y$-direction and a retrograde eccentric disk in the $-y$-direction. The
$N(e)=2e\ de$ simulation (right panel) contains eccentricities from the most
circular to the most eccentric, so it is effectively an intermediate between
the left two panels. Closer in ($a\lesssim 5$), the retrograde population
shows alignment in the $+y$-direction since both circular and eccentric cases
show coherent alignment in this quadrant. The prograde orbits, however, have
seemingly randomized since initially circular orbits align in the opposite
direction ($+y$) as initially eccentric orbits ($-y$). At $a\gtrsim 5$, both
circular and eccentric cases exhibit similar structure so we obtain the bar-
like structure again.

For the isotropic cases (second from top row), the prograde and retrograde
sub-populations must be mirror images of each other due to the axi-symmetry
about the kick vector. For a given eccentricity distribution, one can see how
the isotropic cases are obtained by combining the disk (top row) cases with
their mirror images. For instance for the circular cases (right column), the
“V” shape of the retrograde population around $a\sim 10$ in the disk case is
still very clearly seen in the isotropic case, with the prograde orbits now
showing the same mirror-image structure: an inverse “V” shape around $a\sim
10$. In a similar manner, the $e=0.99$ isotropic case is a mirror-image
version of the disk case. Finally, the $N(e)=2e\ de$ case can be understood as
an intermediate between the circular and $e=0.99$ isotropic cases. Inward of
$a\sim 10$, the alignment is opposite for the initially circular and eccentric
simulations, so the post-kick eccentricity vectors are seemingly random. At
$a\gtrsim 10$, however, the alignment in the circular and eccentric cases are
the same, so we get $+y$-alignment for prograde orbits and $-y$-alignment for
retrograde orbits. It is clear from these plots that the recoil kick induces
distinct orbital structures even when starting from an isotropic
configuration. One point of clarification is that the initially isotropic
simulations exhibit these anisotropies in an axi-symmetric fashion, so
statistically-significant apsidal alignment presents itself as a torus-like
structure of eccentric orbits rather than a bar-like structure.

####  5.1.1 Limitation of the Mean Eccentricity Vector

While the mean eccentricity vector is a convenient measure of lopsidedness
(i.e., if the entire disk is eccentric and apse-aligned), it fails to capture
symmetric anisotropies. As a quick proof of concept let’s consider the bar-
like structure of two anti-aligned eccentric nuclear disks. The mean
eccentricity vector of two anti-aligned eccentric disks of approximately the
same strength in apsidal alignment is $\sim$0 even though the stellar
distribution is clearly anisotropic. This poses a problem when starting with
an isotropic distribution of orbits. By design, the post-kick configuration
will be axi-symmetric about the kick vector, so the mean eccentricity vector
will always average out to $\sim$0.

To capture symmetric anisotropies, we use two measures. The first measure
separates orbits into prograde and retrograde populations about an arbitrary
axis perpendicular to the kick ($+z$). Then, the $y$-component of the mean
eccentricity vector binned by semi-major axis can detect axi-symmetric
anisotropies at different radii. The second measure looks at the tilt vector
alignment with respect to the plane perpendicular to the recoil kick
($y$-$z$). We choose $+y$ to be our reference direction, so one can obtain the
tilt vector components from Equation 5 in this new plane of reference by
simply rotating the axes once over: replacing $e_{x}$ with $e_{y}$, $e_{y}$
with $e_{z}$, and $e_{z}$ with $e_{x}$.

###  5.2 Spherically Isotropic Distribution

We show our two measures of anisotropy in the bottom two rows of Figure 24
with these panels corresponding to the isotropic simulations in the second
from top row. The circular case is in the left panel, the $e=0.99$ case is in
the center panel, and the $N(e)=2e\ de$ case is in the right panel. For both
measures, the shaded noise floor shows the maximum expected deviation from 0
(at a $3\sigma$ level) for an isotropic distribution. The second from the
bottom row shows the mean $y$-component of the eccentricity vector separated
into prograde and retrograde orbits about the $z$-axis and binned by semi-
major axis. When the eccentricity vector is separated in this way, we notice
that statistically-significant apsidal alignment is induced along the
$y$-axis. These panels are also consistent with the distribution of angular
momenta in the second from top row. For the circular case, the prograde orbits
are preferentially aligned in the $+y$-direction for most stars and the
retrograde orbits are preferentially aligned in the $-y$-direction.
Correspondingly, $\langle e_{y}\rangle$ is positive for prograde orbits and
negative for retrograde orbits. For the $e=0.99$ case, the flip in alignment
between prograde and retrograde stars at $a\sim 5$ is captured by $\langle
e_{y}\rangle$ as well. Finally, for the $N(e)=2e\ de$ case, we see alignment
emerging at large semi-major axes, once again consistent with the color-coded
angular momenta in the panel above. Keeping in mind that there is always axi-
symmetry along the kick vector ($+x$), the same eccentricity pattern exists
along the $z$-axis (or any axis perpendicular to the kick vector), so we can
understand these anisotropy patterns as a torus of eccentric orbits that is
perpendicular to the kick vector. We note, however, that this is not a typical
torus in which individual orbits have clustered angular momentum vectors. The
eccentricity vectors lie in the plane perpendicular to the kick, but the
angular momentum vectors are not clustered.

In the bottom row of Figure 24, we show $\langle t_{1}\rangle$, the tilt
vector component that measures roll, with respect to the plane perpendicular
to the recoil ($y$-$z$) binned by semi-major axis. We immediately see that
this measure detects the same anisotropic structure as above. This makes sense
based on how the tilt vector is computed. If we, for example, take the
circular, isotropic case (the bottom three panels of the left-most column),
the prograde orbits (with respect to the $z$-axis) have eccentricity vectors
coherently aligned in $+y$-direction. This means that the ascending node with
respect to the $y$-$z$ plane is preferentially in the $-y$-direction, and
since the eccentricity vector ($+y$) and the ascending node ($-y$) are
misaligned, $\langle t_{1}\rangle<0$ in this case. Another way to understand
this is that $\Omega$ (with respect to the $y$-$z$ plane) is distributed axi-
symmetrically following the kick, but since the tilt vector rotates every
eccentricity vector by $-\Omega$, the axi-symmetry is erased and the torus-
like anisotropies are uncovered.

####  5.2.1 Numerical Results: spherically isotropic, circular

Figure 25: Initially circular orbits in a spherically isotropic distribution;
recoil kick along the $x$-axis. A two-dimensional projection of the post-kick
eccentricity vectors in spherical coordinates where $\theta$ is the polar
angle and $\phi$ is the azimuthal angle. With respect to the $z$-axis,
prograde stars are marked with blue triangles and retrograde stars are marked
with orange “+” markers. Eccentricity vectors preferentially lie in the plane
perpendicular to the recoil kick ($y$-$z$), and this torus-like structure has
a clear division when separated into prograde and retrograde stars about a
certain axis within that plane ($z$-axis). Prograde stars have preferentially
positive $e_{y}$ and retrograde stars have preferentially negative $e_{y}$.
Figure 26: Spherically isotropic distribution with a $N(e)=2e\ de$
eccentricity profile and a Bahcall-Wolf cusp density profile. A histogram of
the post-kick (left:) mean anomaly, $\mathcal{M}$, and (right:) argument of
periapsis, $\omega$. Mean anomaly is clustered at 0 and $\omega$ is clustered
at $\pi$, even for an initially isotropic distribution.

To illustrate the torus-like anisotropy, we show a projection of the post-kick
eccentricity vectors in spherical coordinates for an initially isotropic
distribution of circular orbits in Figure 25, where we have separated prograde
(blue triangle) and retrograde (orange “+” marker) orbits with respect to the
$z$-axis. Just as we saw in Figure 24 in the bottom three panels of the left-
most column, the eccentricity vectors show coherent alignment preferentially
perpendicular to the kick, forming a torus of eccentric orbits. The torus has
a clear division between prograde and retrograde stars: prograde stars have
eccentricity vectors that preferentially point in the $+y$-direction and
retrograde stars have eccentricity vectors preferentially pointing in the
$-y$-direction, consistent with Figure 24. This intuitively makes sense
according to our in-plane and out-of-plane kick simulation results for
initially circular orbits. We have seen that a circular orbit in the $x$-$y$
plane preferentially develops eccentricities along the $y$-direction from
Section 2. Similarly a circular orbit in the $x$-$z$ plane should develop
eccentricities along the $z$-direction. Finally, for an out-of-plane kick with
an orbit initially in the $y$-$z$ plane, eccentricities develop along their
initial position vector which is also contained within the $y$-$z$ plane
(Section 3). It makes sense, then, that an isotropic distribution of circular
orbits will predominantly develop eccentricities in the plane perpendicular to
the kick vector.

####  5.2.2 Numerical Results: spherically isotropic, eccentric

As the most general case, we employ a $N(e)=2e\ de$ eccentricity distribution
and run a recoil kick simulation starting with an isotropic star cluster. As
shown in the bottom three panels of the right-most column in Figure 24, the
post-kick distribution seems to have no discernible structure before
$a=1$–$10$, and beyond $a\gtrsim 10$, a clear eccentric torus forms. As shown
in the bottom row of this figure, tilt vector alignment (with respect to the
plane perpendicular to the kick) is still strong in the isotropic cases. This
is because tilt vector alignment occurs even for a relatively small out-of-
plane kick component (see Figure 23). In Figure 26, we show the post-kick
histogram of mean anomalies, $\mathcal{M}$, in the left panel and arguments of
periapsis, $\omega$, on the right panel. While the distribution is much more
spread-out compared to what we saw for the circular, out-of-plane case in
Figure 8, both distributions are clearly anisotropic: mean anomalies skew
toward 0, which means stars are preferentially starting their new orbits at
periapsis, and the distribution of $\omega$ is skewed toward $\pi$, which
means that the stars periapses preferentially correspond to their descending
node (i.e., moving toward the $-x$-direction at periapsis).

##### Post-Kick Density and Velocity Profile

Figure 27: Spherically isotropic distribution with a $N(e)=2e\ de$
eccentricity profile and a Bahcall-Wolf cusp density profile. The post-kick
distribution of stars in the $x$-$y$ plane at (top left:) $t=0$, (top right:)
$t=0.15$ Myr, and (bottom:) $t=0.75$ Myr. Bound stars are separated into
prograde (blue triangles) and retrograde (orange “+” markers) about the
$z$-axis. The unbound stars are marked with green dots. The “x” marker shows
the center of the galaxy and the “o” marker shows the recoiling black hole.
The recoil kick direction is marked by the arrow labeled $v_{\rm kick}$. All
positions are computed in the galactic center frame of reference. Length scale
and time units assume a black hole mass of $M_{\bullet}=4\times 10^{6}\
M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s. Even for an isotropic
distribution, the prograde and retrograde stars seem to have a clear division
at large radii.

For the most general case of an initially isotropic star cluster with a
$N(e)=2e\ de$ eccentricity profile, we take a look at the density and velocity
profiles to explore any potentially observable structures. In Figure 27, we
show the distribution of stars separated into bound and unbound (green dots)
stars, with the bound stars further separated into prograde (blue triangles)
and retrograde (orange “+” markers) about the $z$-axis, at three different
times post-kick. The “x” marker shows the galactic center and the “o” marker
shows the recoiling black hole. All positions are computed in the galactic
center frame of reference. The evolution is quite similar to what was shown in
Figure 15. At $t=0$ the stars on prograde orbits start toward the
$+y$-direction and the stars on retrograde orbits start in the $-y$-direction.
As the stars move along their Keplerian orbits, the two populations switch
places such that at $t=0.15$ Myr and $t=0.75$ Myr, the prograde population is
now in the $-y$-direction and the retrograde orbits are in the $+y$-direction.
Of course, one key difference between the two simulations is that this
evolution only occurs in the $x$-$y$ plane for the disk case whereas the
initially isotropic case exhibits this structure axi-symmetrically about the
$x$-axis, so the same structure should be observed in the $x$-$z$ plane, for
instance.

Figure 28: Spherically isotropic distribution with a $N(e)=2e\ de$
eccentricity profile and a Bahcall-Wolf cusp density profile. The post-kick
surface density distribution of bound stars in the (top:) $x$-$y$ plane and
(bottom:) $y$-$z$ plane at (left:) $t=0$, (center:) $t=0.1$ Myr, and (right:)
$t=0.2$ Myr. The “o” marker shows the recoiling black hole, and the density
contours are computed in the frame co-moving with the black hole. In each
orientation, the recoil kick direction is marked and labeled as $v_{\rm
kick}$. Length scale and time units assume a black hole mass of
$M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of $v_{\rm kick}=100$ km/s.
The $x$-$y$ plane density contours reveal the “fan”-shape behind the recoiling
black hole found in Merritt et al. (2009). Figure 29: Spherically isotropic
distribution with a $N(e)=2e\ de$ eccentricity profile and a Bahcall-Wolf cusp
density profile. The line-of-sight velocities in the frame co-moving with the
recoiling black hole in the (top:) $x$-$y$ plane and (bottom:) $y$-$z$ plane
at (left:) $t=0$, (center:) $t=0.25$ Myr, and (right:) $t=0.5$ Myr. In each
orientation, the color bar is adjusted such that approaching velocities are
blue and receding velocities are red. The recoil kick direction is marked and
labeled as $v_{\rm kick}$ in each row. Length scale and time units assume a
black hole mass of $M_{\bullet}=4\times 10^{6}\ M_{\odot}$ and a kick of
$v_{\rm kick}=100$ km/s. In the $x$-$y$ plane, there appears a cluster of
stars with high line-of-sight velocities right behind the recoiling black
hole. In the $y$-$z$ plane, the edges become more blue-shifted over time.

In Figure 28, we show the density profiles in two orientations ($x$-$y$ and
$y$-$z$) using a kernel density estimation at three different times following
the recoil. The “o” marker shows the recoiling black hole, and these contours
are in the frame of reference co-moving with the black hole. In the top row,
we see that an initially isotropic distribution forms a “fan” shape behind the
black hole as seen at $t=0.1$ Myr and $t=0.2$ Myr, which is consistent with
density plots presented in Merritt et al. (2009). The density profile
projected onto the $y$-$z$ plane is always axi-symmetric, as expected, and
expands in size over time. Finally, we show the line-of-sight velocity
distribution in the $x$-$y$ and $y$-$z$ planes at three different times post-
kick in Figure 29. Velocities are calculated in the frame co-moving with the
recoiling black hole, and the color bar is adjusted such that approaching
velocities are in blue and receding velocities are in red. In the $x$-$y$
plane, the fan-shape behind the black hole seems to exhibit higher line-of-
sight velocities than the rest of the stellar distribution. This is caused by
extremely eccentric and large semi-major axis orbits that begin their new
orbits at their periapsis, similar to what we found in the center row of
Figure 17. They develop eccentricities perpendicular to the kick vector and
tilt over their major axes such that the periapsis is their descending node
($-x$-direction). Then, shortly after the kick, these stars form a cluster
right behind the black hole in this fan-shape where they exhibit line-of-sight
velocities close to their high periapsis speed. The $y$-$z$ plane shows unique
line-of-sight velocity structures as well. Initially, stars start at their
descending node, so every star is red-shifted. Over time, stars will switch to
being blue-shifted at particular radii depending on their post-kick orbital
period, so that the edges of the cluster are preferentially blue-shifted
($t=0.25$ Myr). This pattern gets weaker over time as the cluster expands in
size and the anomalies are randomized ($t=0.5$ Myr).

##  6 Discussion

In Akiba & Madigan (2021), we explored the post-kick orbits of stars that
remain bound to a recoiling black hole. For a narrow range of initial
conditions, we showed that a lopsided stellar disk emerges post-kick. In this
paper, we expand to more general initial stellar distributions and kicks. Our
findings can be summarized as follows:

  1. 1.

Circular disks which experience an in-plane kick of the black hole map to a
lopsided, eccentric disk with a rich, complex orbital structure. Maximum
orbital alignment occurs at a characteristic post-kick radius

$r_{c}\equiv\frac{4GM_{\bullet}}{9v_{\rm kick}^{2}}\ ,$ (7)

and significant sections of the post-kick disk are counter-rotating with
respect to the disk’s initial angular momentum. The fraction of retrograde
orbits within the bound population tends toward 50% for large post-kick semi-
major axes. The azimuthally-averaged post-kick density profile remains similar
to its pre-kick profile out to $\sim 4\ r_{c}$, then drops sharply.

  2. 2.

Circular disks which experience an out-of-plane kick map to a distinctive
dome-shape. Post-kick, all stars start at the periapsis of their new
elliptical orbits which results in a clustering of mean anomalies at 0. The
angular momentum vectors of the orbits collectively roll over their major axes
in the direction of their initial velocity vectors, which results in
clustering of arguments of periapsis at $\pi$. Eccentricities and inclinations
increase with semi-major axis from the black hole reaching a maximum of 1 and
$45^{\circ}$ for bound orbits, respectively. Stars within a narrow radial
range will have identical values for 5 out of their 6 Kepler elements
($a,e,i,\omega,\mathcal{M}$). Only longitude of ascending node, $\Omega$, will
differ.

  3. 3.

The critical post-kick semi-major axis in which apsidal alignment is strongest
is much smaller for high eccentricity initial conditions,

$a_{c}=\frac{1-e}{1+e}\ r_{c}\ .$ (8)

In disks of high eccentricity orbits which experience an in-plane kick,
roughly half of the stars will map to retrograde orientations post-kick,
except at the smallest semi-major axes. At these very small distances from the
black hole, there exist two anti-aligned lopsided prograde modes. At larger
semi-major axes, there is a prograde eccentric nuclear disk apsidally-aligned
in one direction and a retrograde eccentric disk aligned in the other. This
forms a bar-like structure. The orientations reverse around $a\sim 4$ (Figure
14).

  4. 4.

In disks of high eccentricity orbits which experience an out-of-plane kick,
all orbits tilt significantly. Stellar orbits at low semi-major axes cluster
with $\omega=0$ while those at large semi-major axes cluster with $\omega=\pi$
and gain inclinations of up to 90∘.

  5. 5.

Even for the most general case of stellar orbits being distributed
isotropically about the black hole before recoil, the post-kick distribution
has complex kinematic structure. In particular, irrespective of the initial
stellar orbital eccentricities, the recoil kick produces a torus of eccentric
orbits perpendicular to the kick vector at large semi-major axes (bottom three
rows of Figure 24 and Figure 26). Furthermore, mean anomalies and arguments of
periapsis are non-uniformly distributed.

Mastrobuono-Battisti et al. (2023) recently showed that mergers of
supermassive black holes embedded in nuclear star clusters can leave an
imprint on the shape, density profile, rotation, and velocity structure of the
post-merger nuclear star cluster. While the density profile is sensitive to
the initial merger conditions, the post-merger nuclear star cluster is always
found to be flattened and rotating (expected due to the conserved orbital
angular momentum of its progenitors). The post-merger cluster is also
tangentially anisotropic at its center. These are the initial conditions —
flattened, rotating, with central tangential anisotropy — that we might expect
for a cluster that then responds to a recoil kick when the supermassive black
holes merge at the center. The simulation in this paper that best matches this
initial condition is a low-eccentricity disk with $\sigma_{i}=30^{\circ}$
which corresponds to an oblate spheroid with an axial ratio between
$0.6$–$0.7$. As can be seen in Figure 23 and the top panel of Figure 24, these
initial conditions result in statistically significant apsidal and/or tilt
alignment of post-kick orbits. In particular, at large post-kick semi-major
axes we would expect anti-aligned prograde and retrograde disks of high
eccentricity orbits.

Disks of stars responding to in-plane components of recoil kicks end up with a
significant fraction of their outer orbits with retrograde orientation (see
bottom panels of Figures 2 and 14). This type of counter-rotation in an edge-
on stellar ring surrounding a central black hole has already been observed in
an early-type galaxy NGC 3706 (Gültekin et al., 2011).

In this paper, we have presented the instantaneous orbital distributions of
stars following a recoil kick of the central black hole. We leave as future
work the effects of mass loss due to gravitational wave emission (e.g. Penoyre
& Haiman, 2018) and the long-term evolution of these systems. We expect many
configurations to evolve dramatically with time once self-gravity between the
stars is included. For example, the anti-aligned eccentric nuclear disks we
find in Section 4.3 should gravitationally torque each other to extreme
eccentricities. Furthermore, counter-rotating disks are unstable to buckling
out of the plane (Toomre, 1964; Merritt & Sellwood, 1994). For the case in
which stars are initially on circular orbits, those that map to retrograde
orientations originate further from the black hole. If the stellar population
is initially mass segregated with heavier stars and compact objects located
closer to the black hole, we would expect the post-kick prograde population to
be composed of heavier bodies and the retrograde population to be lighter.

##  7 Acknowledgements

AM 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.

## References

  * Akiba & Madigan (2021) Akiba, T., & Madigan, A.-M. 2021, ApJ, 921, L12
  * Bahcall & Wolf (1976) Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
  * Bekenstein (1973) Bekenstein, J. D. 1973, ApJ, 183, 657
  * Binney & Tremaine (1987) Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton University Press)
  * Blanchet et al. (2005) Blanchet, L., Qusailah, M. S. S., & Will, C. M. 2005, ApJ, 635, 508
  * Blecha & Loeb (2008) Blecha, L., & Loeb, A. 2008, MNRAS, 390, 1311
  * Bogdanović et al. (2007) Bogdanović, T., Reynolds, C. S., & Miller, M. C. 2007, ApJ, 661, L147
  * Campanelli et al. (2007) Campanelli, M., Lousto, C. O., Zlochower, Y., & Merritt, D. 2007, Phys. Rev. Lett., 98, 231102
  * Chiaberge et al. (2018) Chiaberge, M., Tremblay, G. R., Capetti, A., & Norman, C. 2018, ApJ, 861, 56
  * Comerford et al. (2017) Comerford, J. M., Barrows, R. S., Greene, J. E., & Pooley, D. 2017, ApJ, 847, 41
  * González et al. (2007) González, J. A., Sperhake, U., Brügmann, B., Hannam, M., & Husa, S. 2007, Phys. Rev. Lett., 98, 091101
  * Goodman (2003) Goodman, J. 2003, MNRAS, 339, 937
  * Gualandris & Merritt (2008) Gualandris, A., & Merritt, D. 2008, ApJ, 678, 780
  * Gültekin et al. (2011) Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2011, ApJ, 741, 38
  * Herrmann et al. (2007) Herrmann, F., Hinder, I., Shoemaker, D., Laguna, P., & Matzner, R. A. 2007, ApJ, 661, 430
  * Hogg et al. (2021) Hogg, J. D., Blecha, L., Reynolds, C. S., Smith, K. L., & Winter, L. M. 2021, MNRAS, 503, 1688
  * Kim et al. (2018) Kim, D. C., Yoon, I., & Evans, A. S. 2018, ApJ, 861, 51
  * Komossa & Merritt (2008a) Komossa, S., & Merritt, D. 2008a, ApJ, 689, L89
  * Komossa & Merritt (2008b) —. 2008b, ApJ, 683, L21
  * Komossa et al. (2008) Komossa, S., Zhou, H., & Lu, H. 2008, ApJ, 678, L81
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
  * Lauer et al. (1993) Lauer, T. R., Faber, S. M., Groth, E. J., et al. 1993, AJ, 106, 1436
  * Lauer et al. (2005) Lauer, T. R., Faber, S. M., Gebhardt, K., et al. 2005, AJ, 129, 2138
  * Lippai et al. (2008) Lippai, Z., Frei, Z., & Haiman, Z. 2008, ApJ, 676, L5
  * Lousto & Zlochower (2008) Lousto, C. O., & Zlochower, Y. 2008, Phys. Rev. D, 77, 044028
  * Madau & Quataert (2004) Madau, P., & Quataert, E. 2004, ApJ, 606, L17
  * Madigan et al. (2018) Madigan, A.-M., Halle, A., Moody, M., et al. 2018, ApJ, 853, 141
  * Madigan & McCourt (2016) Madigan, A.-M., & McCourt, M. 2016, MNRAS, 457, L89
  * Mastrobuono-Battisti et al. (2023) Mastrobuono-Battisti, A., Ogiya, G., Hahn, O., & Schultheis, M. 2023, MNRAS, doi:10.1093/mnras/stad898
  * Merritt et al. (2004) Merritt, D., Milosavljević, M., Favata, M., Hughes, S. A., & Holz, D. E. 2004, ApJ, 607, L9
  * Merritt et al. (2009) Merritt, D., Schnittman, J. D., & Komossa, S. 2009, ApJ, 699, 1690
  * Merritt & Sellwood (1994) Merritt, D., & Sellwood, J. A. 1994, ApJ, 425, 551
  * Paczynski (1978) Paczynski, B. 1978, Acta Astronomica, 28, 91
  * Penoyre & Haiman (2018) Penoyre, Z., & Haiman, Z. 2018, MNRAS, 473, 498
  * Peres (1962) Peres, A. 1962, Physical Review, 128, 2471
  * Rein & Liu (2012) Rein, H., & Liu, S.-F. 2012, A&A, 537, A128
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217
  * Tremaine (1995) Tremaine, S. 1995, AJ, 110, 628
  * Wiseman (1992) Wiseman, A. G. 1992, Phys. Rev. D, 46, 1517

##  Appendix A Tilt Vector and $i_{a}$, $i_{b}$, and $i_{e}$ Relations

Two angles are required to orient the plane of an orbit in space. In Madigan &
McCourt (2016) we introduced two inclination angles, $i_{a}$ and $i_{b}$.
These angles represent rotations of an orbit with angular momentum vector
($\hat{j}$) about its semi-major axis ($\hat{a}$) and its semi-minor axis (or
semi-latus rectum) ($\hat{b}\equiv\hat{j}\times\hat{a}$), respectively. An
additional angle is needed to specify the direction of the orbit within the
orbital plane. We use $i_{e}$ which is the rotation of an orbit about the
angular momentum vector. The angles $i_{a}$, $i_{b}$, and $i_{e}$ are
equivalent to that of an aircraft’s roll, pitch and yaw.

Here we relate these angles to the tilt vector as introduced in Section 3.2.
The tilt vector results from the rotation of the eccentricity vector by an
angle $-\Omega$, so taking the $x$-$y$ plane to be our reference plane and
$+x$ to be our reference direction,

$\displaystyle t_{1}=\cos(\Omega)e_{x}+\sin(\Omega)e_{y}\ ,$ $\displaystyle
t_{2}=-\sin(\Omega)e_{x}+\cos(\Omega)e_{y}\ ,$ $\displaystyle t_{3}=e_{z}\ .$

This allows us to view the eccentricity vector from the ascending node of the
orbit rather than from a common reference direction. The ratio between $t_{1}$
and $t_{3}$ indicates the degree to which the orbit is tilting over their
major axis versus minor axis.

###  A.1 $i_{a}$ is related to $t_{1}$ and they measure roll

$i_{a}$ is defined as

$i_{a}=\arctan\left(\frac{\hat{b}_{z}}{\sqrt{1-\hat{b}_{z}^{2}}}\right)\ .$

We note that the longitude of ascending node, $\Omega$, is always $90^{\circ}$
counter-clockwise in the $x$-$y$ plane from $\vec{j}_{\perp}$, the projection
of $\vec{j}$ onto the $x$-$y$ plane. Since the tilt vector rotates the orbit
such that the ascending node coincides with the positive $x$-axis, this means
that the angular momentum vector becomes rotated such that it resides in the
$y$-$z$ plane. The angular momentum vector post-rotation is

$\vec{j}^{\prime}=\left\langle
0,-\sqrt{j_{x}^{2}+j_{y}^{2}},j_{z}\right\rangle\ .$

The definition of the tilt vector is the same rotation applied to the
eccentricity vector, so $\vec{e}^{\prime}=\vec{t}$. Then, we can calculate
$\hat{b}_{z}$ as

$\displaystyle\hat{b}_{z}$
$\displaystyle=\hat{j}^{\prime}_{x}\hat{t}_{2}-\hat{j}^{\prime}_{y}\hat{t}_{1}$
$\displaystyle=-\frac{j_{\perp}}{j}\frac{t_{1}}{e}\ ,$

so $i_{a}$ can be written in terms the first component of the tilt vector,
$t_{1}$.

###  A.2 $i_{b}$ is related to $t_{3}$ and they measure pitch

$i_{b}$ is given by

$i_{b}=\arctan\left(-\frac{\hat{a}_{z}}{\sqrt{1-\hat{a}_{z}^{2}}}\right)\ ,$

and since $\hat{a}=\hat{e}$, the unit eccentricity vector, $i_{b}$ can be
equivalently written as

$i_{b}=\arctan\left(-\frac{\hat{e}_{z}}{\sqrt{1-\hat{e}_{z}^{2}}}\right)\ .$

Furthermore, $\hat{e}_{z}=\hat{t}_{3}$, so

$i_{b}=\arctan\left(-\frac{\hat{t}_{3}}{\sqrt{1-\hat{t}_{3}^{2}}}\right)\ .$

###  A.3 $i_{e}$ is related to the ratio $t_{2}/t_{1}$ and they measure yaw

$i_{e}$ is defined as

$i_{e}=\arctan\left(\frac{\hat{e}_{y}}{\hat{e}_{x}}\right)\ ,$

and this measures yaw with respect to the $x$-axis. A similar definition of
$i_{e}$ using the tilt vector would be

$i_{e}=\arctan\left(\frac{\hat{t}_{2}}{\hat{t}_{1}}\right)\ ,$

which measures yaw with respect to the longitude of ascending node. This is
effectively a measure of yaw that is independent of our choice of reference
direction.