# Signatures of Massive Black Hole Merger Host Galaxies from Cosmological
Simulations II:
Unique Stellar Kinematics in Integral Field Unit Spectroscopy

Jaeden Bardati Department of Physics & Astronomy, Bishop’s University,
Sherbrooke, QC J1M 1Z7, Canada  The Division of Physics, Mathematics and
Astronomy, California Institute of Technology, Pasadena, CA 91125, USA  John
J. Ruan Department of Physics & Astronomy, Bishop’s University, Sherbrooke, QC
J1M 1Z7, Canada  Daryl Haggard McGill Space Institute and Department of
Physics, McGill University, 3600 rue University, Montreal, Quebec, H3A 2T8,
Canada  Michael Tremmel School of Physics, University College Cork, Cork,
Ireland  Patrick Horlaville Department of Physics & Astronomy, Bishop’s
University, Sherbrooke, QC J1M 1Z7, Canada

###### Abstract

Secure methods for identifying the host galaxies of individual massive black
hole (MBH) binaries and mergers detected by gravitational wave experiments
such as LISA and Pulsar Timing Arrays are currently lacking, but will be
critical to a variety of science goals. Recently in Bardati et al. (2024,
Paper I), we used the Romulus25 cosmological simulation to show that MBH
merger host galaxies have unique morphologies in imaging, due to their
stronger bulges. Here, we use the same sample of simulated MBH merger host
galaxies to investigate their stellar kinematics, as probed by optical
integral field unit (IFU) spectroscopy. We perform stellar population
synthesis and dust radiative transfer to generate synthetic 3D optical
spectral datacubes of each simulated galaxy, and produce mock stellar
kinematic maps. Based on a linear discriminant analysis of a combination of
kinematic parameters derived from these maps, we show that this approach can
identify MBH binary and merger host galaxies with accuracies that increase
with chirp mass and mass ratio. For mergers with high chirp masses
($\gtrsim$108.2 $M_{\odot}$) and high mass ratios ($\gtrsim$0.5), the
accuracies reach $\gtrsim$85%, and their host galaxies are uniquely
characterized by slower rotation and stronger stellar kinematic misalignments.
These kinematic properties are commonly associated with massive early-type
galaxies that have experienced major mergers, and naturally act as signposts
for MBH binaries and mergers with high chirp masses and mass ratios. These
results suggest that IFU spectroscopy should also play a role in telescope
follow-up of future MBH binaries and mergers detected in gravitational waves.

black holes – galaxies: structure – gravitational waves – N-body simulations –
radiative transfer

††software:  SKIRT: Camps & Baes (2015); Pynbody: Pontzen et al. (2013);
Tangos: Pontzen & Tremmel (2018); FSPS: Conroy & Gunn (2010); ppXF: Cappellari
(2023); StatMorph: Rodriguez-Gomez et al. (2019); sklearn: Pedregosa et al.
(2011); ChaNGa: Menon et al. (2015); Amiga Halo Finder: Knollmann & Knebe
(2009); astropy: Astropy Collaboration et al. (2018)

##  1 Introduction

Mergers of massive black holes (MBHs, with masses of $M_{\mathrm{BH}}\gtrsim
10^{5}M_{\odot}$) at the centers of galaxies are inevitable in our standard
picture of hierarchical galaxy formation (Begelman et al., 1980; Volonteri et
al., 2003; Volonteri & Natarajan, 2009). Following a merger of two galaxies,
their two central MBHs will sink to the gravitational potential minimum due to
dynamical friction, and form a gravitationally-bound MBH binary (Colpi et al.,
1999; Yu, 2002; Pfister et al., 2017). This binary will then harden in the
nuclear region of the galaxy, due to some combination of three-body
interactions with stars (e.g., Quinlan, 1996; Merritt & Poon, 2004; Berczik et
al., 2006; Vasiliev et al., 2015; Sesana & Khan, 2015) or torques from gas
(e.g., Armitage & Natarajan, 2002; Escala et al., 2004; Dotti et al., 2007;
Mayer et al., 2007; Cuadra et al., 2009; Chapon et al., 2013; Fontecilla et
al., 2019). Eventually, angular momentum loss becomes dominated by
gravitational wave emission, which drives the MBH binary to coalescence
(Wyithe & Loeb, 2003; Sesana et al., 2004).

MBH binaries and mergers are already beginning to be detected in low-frequency
gravitational waves. Current pulsar timing array (PTA; Jenet et al., 2004)
experiments such as the International Pulsar Timing Array (Verbiest et al.,
2016) have recently detected the nHz stochastic gravitational wave background
(Agazie et al., 2023a), produced by the ensemble population of MBH _binaries_
in nearby galaxies (Agazie et al., 2023b). Although these experiments have not
yet detected an individual MBH binary system (Agazie et al., 2023c, d),
forecasts predict that they will soon achieve sufficient sensitivity for a
first detection (Sesana et al., 2009; Ravi et al., 2015; Mingarelli et al.,
2017; Kelley et al., 2018). Future mHz gravitational wave experiments such as
the _Laser Interferometer Space Antenna_ (_LISA_ ; Amaro-Seoane et al., 2017)
will detect _mergers_ of MBH with mass $10^{4}M_{\odot}\lesssim
M_{\mathrm{BH}}\lesssim 10^{8}M_{\odot}$, to high redshifts of $\gtrsim$10
(Sesana et al., 2005; Banks et al., 2022). These discoveries will represent a
breakthrough in our ability to probe the formation and growth of MBHs over
cosmic time.

Once an individual MBH binary or merger is detected in gravitational waves,
identifying its host galaxy will be challenging. Identifying the exact host
galaxy of each MBH binary or merger will enable a variety of key science
goals, such as unveiling the structure of accretion flows around binary MBH
(e.g., Armitage & Natarajan, 2002; Cuadra et al., 2009; Noble et al., 2012;
Duffell et al., 2020), constraining cosmological parameters (Schutz, 1986;
Holz & Hughes, 2005), and probing the co-evolution of galaxies and their
central MBHs (Magorrian et al., 1998; Ferrarese & Merritt, 2000; Gebhardt et
al., 2000). However, gravitational wave detections will have relatively poor
sky localizations that contain an overwhelmingly large number of possible host
galaxies in follow-up telescope observations. Goldstein et al. (2019) and
Petrov et al. (2024) estimate that the error volumes of individual MBH
binaries detected by PTAs will contain of order $\sim$102 candidate host
galaxies, after both redshift and mass cuts (assuming empirical galaxy scaling
relations) based on the gravitational wave information. Similarly, Lops et al.
(2023) show that the error volumes of MBH mergers detected by _LISA_ will
contain of order $\sim$102-3 candidate host galaxies, for mergers of
relatively-massive MBHs ($\gtrsim$107 $M_{\odot}$) at relatively-low redshifts
($z\lesssim 2$), which are most amenable to telescope follow-up. Clearly,
additional approaches to selecting the best sources for follow-up will be
needed to identify the exact host.

Many strategies identifying the host galaxies of MBH binaries and mergers have
been suggested, although secure methods are still lacking. Previous studies
have often focused on transient or time-variable electromagnetic signatures
from gas surrounding MBH binaries (e.g., D’Orazio et al., 2015; Graham et al.,
2015; Kelley et al., 2019; Charisi et al., 2022; Milosavljević & Phinney,
2005; Megevand et al., 2009; O’Neill et al., 2009; Corrales et al., 2010;
Rossi et al., 2010) that can be observed in telescope follow-up (see recent
reviews by Bogdanović et al., 2022; D’Orazio & Charisi, 2023). However, these
time-domain signatures are highly uncertain, and may be applicable to only a
small fraction of MBH binaries and mergers if the majority of their
surrounding environments are too gas-poor to support accretion (Izquierdo-
Villalba et al., 2023b; Dong-Páez et al., 2023), or are heavily enshrouded in
dust that obscures the electromagnetic emission from the accretion flow (Koss
et al., 2018).

Recently in Bardati et al. (2024, hereafter Paper I), we developed an
alternative approach to identifying the host galaxies of MBH binaries and
mergers based on galaxy morphology in telescope imaging, thus demonstrating
that the host galaxy stellar emission contain unique signatures. Specifically,
we used simulated MBH merger host galaxies from a cosmological simulation, and
performed dust radiative transfer to produce synthetic broadband ultraviolet,
optical, and infrared images. Using a variety of morphological parameters
extracted from these synthetic images, we trained a Linear Discriminant
Analysis (LDA) predictor to distinguish between MBH merger host galaxies and a
control sample. We showed that the accuracy of this method increases with both
chirp mass and mass ratio of the MBH binary, reaching $\gtrsim$80% for mergers
with high chirp masses ($\gtrsim$108.2 $M_{\odot}$) and high mass ratios
($\gtrsim$0.5). We showed that this is because the primary signatures of these
host galaxies are prominent classical bulges, which are built from major
mergers of massive galaxies; these galaxy mergers naturally lead to the
formation of MBH binaries and mergers with high mass ratios and high chirp
masses. These results are broadly consistent with other studies that use
cosmological simulations to show that the host galaxies of MBH binaries
detectable by PTA experiments (which preferentially have high mass ratios and
high chirp masses) are massive early-type galaxies with low star formation
rates and high metallicities (e.g., Saeedzadeh et al., 2024; Cella et al.,
2024), while the host galaxies of MBH mergers detectable by LISA (which have
low chirp masses) are not morphologically distinct (Izquierdo-Villalba et al.,
2023a). Although our results from Paper I imply that there are indeed
observable signatures of MBH binaries and mergers with high mass ratios and
high chirp masses in their host galaxy starlight, it is unclear whether unique
host galaxy morphologies in imaging is actually the most robust observational
signature.

Here, we investigate the unique signatures of MBH binaries and mergers in
their host galaxy _stellar kinematics_ (in contrast to morphology in Paper I),
as probed by optical integral field unit (IFU) spectroscopy with medium
spectral resolution. Although IFU spectroscopy is significantly more resource-
intensive than imaging, a stellar kinematics approach may have unique
advantages over morphology. For example, if the primary distinguishing
characteristic of MBH merger host galaxies is indeed the presence of a
prominent bulge, then this feature will likely be more distinctive in stellar
kinematic maps than morphology in imaging. Furthermore, Nevin et al. (2021)
use simulations to show that stellar kinematic signatures of galaxy mergers
persist for several Gyrs, after morphological disturbances have disappeared
and when the MBH binary or merger would be detected in gravitational waves.
Thus, it is possible that the stellar kinematics of MBH binary and merger host
galaxies can contain unique signatures from _both_ permanent features such as
bulges, as well as kinematic disturbances from their preceding galaxy merger.

To assess whether stellar kinematics can identify MBH binary and merger host
galaxies, we use a sample of 201 simulated MBH merger host galaxies from a
cosmological simulation, as well as a mass- and redshift-matched control
sample for comparison. The samples we use here are the exact same samples we
previously used in Paper I. However, the post-processing and analysis of these
simulated galaxies here are distinct from Paper I, since we are producing and
analyzing synthetic 3D spectral datacubes instead of synthetic images. We
perform stellar population synthesis and dust radiative transfer to generate
synthetic spectral datacubes of each simulated galaxy, similar to observations
from optical integral field unit (IFU) spectrographs. We then perform full-
spectrum fitting of these spectral datacubes to produce 2D stellar kinematic
maps, and extract a variety of kinematic parameters for each simulated galaxy.
By training a LDA predictor on these kinematic parameters, we find a linear
combination of kinematic parameters that can optimally distinguish between the
MBH merger host galaxies and the control sample with high accuracies. Our
results thus motivate follow-up IFU observations (or analysis of archival IFU
data) of galaxies lying in the localization regions of MBH binaries and
mergers detected in gravitational waves.

The outline of this paper is as follows. In Section 2, we briefly describe the
cosmological simulation we use, and our selection of samples of simulated
galaxies. In Section 3, we describe the radiative transfer simulations we
perform to produce 3D spectral datacubes, as well our fitting of these
datacubes to produce 2D stellar kinematic maps. In Section 4, we describe the
various kinematic parameters we extract, and our linear discriminant analysis.
In Section 5, we present our main results, and use them to build a consistent
picture of the unique signatures of MBH binary and merger host galaxies. We
briefly summarize and conclude in Section 6.

##  2 Cosmological Simulation and Galaxy Samples Selection

In this section, we describe the cosmological simulation, and our selection of
samples of simulated galaxies. We note that the simulation and galaxy samples
we use are identical to those used in Bardati et al. (2024, Paper I). This
enables a direct comparison between previous results from the galaxy
morphology approach in Paper I to our stellar kinematics approach here. We
thus only briefly describe the most salient points below, and refer to Paper I
for more details.

###  2.1 Romulus25

For our analysis, we use the Romulus25 cosmological simulation of galaxy
formation (Tremmel et al., 2017). The simulation is run to $z=0$ using the
N-Body + Smooth Particle Hydrodynamics code ChaNGa (Menon et al., 2015), with
a 25 cMpc per side volume and a standard $\Lambda\mathrm{CDM}$ cosmology
consistent with Planck Collaboration (2016).

Romulus25 introduces new implementations of MBH formation, growth, and
dynamics. MBHs in Romulus25 are seeded from direct collapse with masses of
$M_{BH}=10^{6}M_{\odot}$. MBH formation is directly tied to the properties of
the surrounding gas particles, such as metallicity, gas density, and
temperature, which better reproduces MBH seeding at high redshift than methods
tied to halo mass. After formation, the MBHs grow both via accretion of
surrounding gas at a rate determined by a modified Bondi-Hoyle prescription,
and via numerical mergers with other MBHs when their mutual distance is less
than two softening lengths ($\lessapprox$700 pc). These MBH binaries are
tracked to sub-kpc scales using a sub-grid model of dynamical friction
(Tremmel et al., 2015), in contrast to the more common approach of
artificially pinning the MBHs to their host galaxy potential minimum at each
timestep. This allows MBHs to decouple from their host galaxy dynamics, thus
enabling MBH dynamics to be tracked down to $\lessapprox$700 pc before merging
numerically below the resolution limit of the simulation. This feature of
Romulus25 leads to shorter time-delays between the numerical merger and the
physical merger, when gravitational waves would be emitted. We refer to
Tremmel et al. (2017) for more details on Romulus25.

Jung et al. (2022) investigated the properties of Brightest Group Galaxies
(BGGs) in Romulus25, and found that their stellar kinematics are an excellent
match to observations, especially for BGG in lower-mass groups with halo
masses of $<$1013 $M_{\odot}$. These lower-mass BGGs have halo masses similar
to our subsample of MBH merger host galaxies (and the corresponding control
galaxy sample) with chirp masses of $\gtrsim$108.2 $M_{\odot}$ that we
highlight below in Sections 4.2 and 5. Thus, the results of Jung et al. (2022)
suggests that the stellar kinematics of simulated galaxies from Romulus25 are
realistic and robust.

###  2.2 MBH Merger Host Galaxies Sample

We select a sample of 201 MBH numerical mergers in Romulus25, consisting of
all MBH merger events that occur at a redshift of $z\leq 2$, with a mass ratio
of $\nicefrac{{M_{2}}}{{M_{1}}}>0.1$ (where $M_{1}>M_{2}$) and with MBH masses
of $M_{\text{BH}}\gtrsim 10^{7}M_{\odot}$. These cuts select relatively
massive and low-$z$ galaxies, suitable for follow-up telescope observations.

We use the Amiga Halo Finder (Knollmann & Knebe, 2009) to identify the host
galaxy of each MBH merger at the nearest available simulation snapshot in
time, subsequent to the numerical merger. Since the time-delay between
numerical and physical merger is poorly-understood, we also track these
galaxies for $\sim$1 Gyr after their MBH merger to investigate the time
evolution of our LDA results in Section 5.3 below. In Paper I, we showed that
this sample of MBH merger host galaxies follows the typical scaling relation
between central MBH mass ($M_{\text{BH}}$) and stellar mass ($M_{*}$),
consistent with other studies which investigate MBH mergers using cosmological
simulations (e.g., DeGraf et al., 2021).

###  2.3 Control Galaxy Sample

We select a mass- and redshift-matched control sample of galaxies in Romulus25
that do not host MBH mergers, for comparison to our MBH merger host sample. To
produce this control sample, we first bin the simulation snapshot timesteps
such that each time bin contains a minimum of 30 MBH mergers, sufficient to
produce a MBH mass distribution. For each timestep bin, we then select
galaxies for a control sample from the non-MBH-merger galaxy mass
distribution, weighted by the ratio of the number of MBH merger galaxies to
the number of non-MBH merger galaxies for each mass bin. The resulting mass-
and redshift-matched control sample avoids biases in our results stemming from
galaxy mass and redshift.

Figure 1: Stellar kinematic maps of a representative simulated galaxy from our
MBH merger hosts sample (top panels), as well as from our control sample
(bottom). The top left panel for each of the two galaxies displays the stellar
line of sight velocity ($v_{\mathrm{los}}$) and the top right panel shows the
stellar velocity dispersion ($\sigma_{\star}$). The bottom panels show an
example synthetic spectrum from the indicated spaxel in the datacube (black
line), as well as the fitted model from ppXF (red line). It is difficult to
visually classify these galaxies from these kinematic maps, and thus we
extract a variety of kinematic parameters in Section 4.1 to input into a
linear discriminant analysis predictor.

##  3 Generating 3D Spectral Datacubes and Stellar Kinematic Maps

In this section, we describe the dust radiative transfer simulation we perform
to generate synthetic 3D spectral datacubes of the simulated galaxies from
Sections 2.2 and 2.3 above, and our fitting of these datacubes to produce mock
stellar kinematic maps. We emphasize that in both these steps, we do not
attempt to carefully match the specifications of any particular telescope or
instrument. Instead, our work here is a first step to gauge the efficacy of
IFU spectroscopy for identifying MBH merger host galaxies, and we aim to make
our results broadly applicable to a wide range of instruments. Our inherent
assumptions in generating the spectral datacubes (e.g., medium-resolution
spectroscopy and nearly diffraction-limited seeing with $0\farcs 1$ pixel
scales) are already viable with current IFU instruments such as the Multi Unit
Spectroscopic Explorer (MUSE; Bacon et al., 2010) on the Very Large Telescope,
and achievable with relatively short exposure times using future instruments
such as High Angular Resolution Monolithic Optical and Near-infrared
Instrument (HARMONI; Thatte et al., 2020) on the Extremely Large Telescope. We
discuss some practical considerations for applying our results to new and
archival IFU observations in Section 5.5, and defer a more careful matching to
instrument specification (e.g., using the HSIM3 simulator for HARMONI
observations; Zieleniewski et al. 2015) for future work.

###  3.1 Dust Radiative Transfer with SKIRT

We use the SKIRT 3D Monte Carlo dust radiative transfer software (Baes et al.,
2011; Camps & Baes, 2015) to generate synthetic spectral datacubes for each
simulated galaxy. SKIRT performs stellar population synthesis on each star
particle of a simulated galaxy, and uses the resulting SEDs to propagate
photon packets through the galaxy’s interstellar medium to produce spatially-
resolved galaxy SEDs. Critically, the current version 9 of SKIRT includes
stellar and gas kinematics in the radiative transport by accounting for
Doppler shifts due to particle motion in the photon packets (Camps & Baes,
2020). This enables us to accurately generate line-of-sight velocity
dispersion maps in Section 3.2.

To first perform the stellar population synthesis, we use the Flexible Stellar
Population Synthesis (FSPS; Conroy & Gunn, 2010) package included in SKIRT,
which models each star particle as a single-age simple stellar population. In
FSPS, we assume a Kroupa (2001) initial mass function, and we use the MIST
evolutionary tracks (Choi et al., 2016). To produce a spectrum from the single
stellar population of each star particle, we use the MILES (Sánchez-Blázquez
et al., 2006) spectral libraries (Vazdekis et al., 2010). Since we are
interested primarily in the stellar kinematics, we neglect stellar and nebular
emission lines, and we do not include emission from an active galactic
nucleus.

For SKIRT to perform dust radiative transport, the properties of the dust must
be modeled. We thus assume that the spatial distribution of the dust follows
that of the gas particles in each galaxy. We use a THEMIS dust mix (Jones et
al., 2017) to represent our dust grain size and composition, and assume a
proportionality factor of 0.2 on the grain size distribution function. In our
tests, we find that our results are not strongly dependent on our choice of
these dust parameters.

We use SKIRT to produce output spectra over the rest-frame wavelength range of
3400 – 7000 Å, with a spectral resolution of $R=4000$. This wavelength range
includes prominent stellar absorption lines (such as the H Balmer series and
the Ca II H and K lines) that are commonly used to produce stellar light-of-
sight velocity and velocity dispersion maps in IFU observations. These spectra
also enable us to produce synthetic $u$-, $g$-, and $r$-band images, which are
used in our calculation of some of the kinematic parameters described in
Section 4.1. For our synthetic spectral datacubes, we assume a pixel scale of
$0\farcs 1$ pixel-1 to match the specifications used in Paper I for direct
comparison. At the median redshift ($z\sim 1$) of the simulated galaxies in
our sample, each pixel corresponds to a physical scale of approximately 0.8
kpc, larger than the gravitational softening scale of Romulus25.

We compute uncertainties on the flux density values in each synthetic spectrum
using reliability statistics output from SKIRT, following prescriptions from
Camps & Baes (2018) and Camps & Baes (2020). In addition to measuring the
accumulated weights $w_{i}$ of the $N$ individual photon packets launched by
the Monte Carlo radiative transfer simulation in each pixel
$\sum_{i}^{N}w_{i}$, we also track the higher-order moments
$W_{k}\equiv\sum_{i}^{N}w_{i}^{k}$ for $k\in\\{2,3,4\\}$. With these moments,
we compute the relative error

$R=\sqrt{\frac{W_{2}}{W_{1}^{2}}-\frac{1}{N}},$

as well as the variance of the variance

$VOV=\frac{W_{4}-4W_{1}W_{3}/N+8W_{2}W_{1}^{2}/N^{2}-4W_{1}^{4}/N^{3}-W_{2}^{2}/N}{(W_{2}-W_{1}^{2}/N)^{2}},$

which is a measure of the uncertainty in $R$. We use these quantities as a
measure of the uncertainty on the flux in each pixel, and to mask bad pixels
(see Section 3.2).

To minimize the effects of galaxy orientation on our results, we generate
spectral datacubes along four different isotropically-oriented viewing angles
for each galaxy. Specifically, for each simulated galaxy, we generate spectral
datacubes from the perspective of each of the four vertices of a regular
tetrahedron. Since each galaxy is randomly oriented, the tetrahedron is also
effectively randomly oriented, and thus the resultant datacubes sample a
variety of different viewing angles.

###  3.2 Producing Stellar Kinematic Maps with ppXF

To produce a stellar kinematics map of each galaxy, we fit its spectral
datacube using the penalized-pixel fitting software ppXF (Cappellari, 2023).
ppXF uses single stellar population templates to perform full-spectrum fitting
of galaxy spectra. The stellar templates are convolved with a line-of-sight
velocity distribution to fit kinematic properties from the observed absorption
lines. To include the effects of a Point Spread Function (PSF), we convolve
the monochromatic image at each wavelength slice of the spectral datacubes
with a $0\farcs 1$ full-width at half-maximum (FWHM) Gaussian. Since we are
not directly comparing our synthetic spectra datacubes to observed data, we do
not include the effects of a line-spread function. Before fitting, we mask out
all spectra in each datacube that have $>$10% of pixels with relative error
$R>0.1$ and variance of variance $VOV>0.1$; SKIRT spectra satisfying these
conditions are known to be of insufficient quality and unreliable for fitting
(Camps & Baes, 2015).

For each datacube, we fit the unmasked spectra in Voronoi bins to obtain a
line-of-sight velocity $v_{\mathrm{los}}$ and velocity dispersion
$\sigma_{\star}$ for the stellar component. Following observational
approaches, we first Voronoi bin each spectral datacube using VorBin
(Cappellari & Copin, 2003) to achieve a minimum signal-to-noise ratio of 50 in
each bin at a fiducial wavelength of 5000 Å. We then use ppxf to fit the
combined spectra in each Voronoi bin, while masking out the pixels in each
combined spectrum that have relative error $R>0.1$ and variance of variance
$VOV>0.1$. Since our SKIRT spectra do not include nebular emission, we do not
include any emission line components in our spectral fitting. We fit the
spectra over the 3450 – 6000 Å wavelength range, using the empirical MILES
(Sánchez-Blázquez et al., 2006) library of single stellar population template
spectra (Vazdekis et al., 2010), and include the two-parameter dust model from
Cappellari (2023). In our fits, we include the 1$\sigma$ uncertainties on the
flux densities in the spectra from photon noise in the SKIRT radiative
transfer. We correct the final $\sigma_{\star}$ values for differences in
spectral resolution between the templates and simulated spectra, following
Cappellari (2017).

Figure 1 shows the resulting $v_{\mathrm{los}}$ and $\sigma_{\star}$ stellar
kinematic maps of an example galaxy from our MBH merger sample, as well as a
similar galaxy from the control sample. Figure 1 also displays an example
spectrum from a central spaxel in both datacubes, along with the ppXF fit. In
all of the Voronoi bins, the full-spectrum fits are robust, with
$\chi^{2}/dof\approx 1$. We will use these these stellar kinematic maps to
extract a variety of kinematic parameters in Section 4, and use a combination
of these parameters to identify MBH merger host galaxies.

Figure 2: Our fitting of the line of sight velocity $v_{\mathrm{los}}$ map
(left) and stellar velocity dispersion $\sigma_{\star}$ map (right) to extract
stellar kinematic properties for the example MBH merger host galaxy from
Figure 1. The bottom left panel shows the $v_{\mathrm{los}}$ map produced by
ppXF, along with the kinematic position angle fitted using the
$v_{\mathrm{los}}$ map (blue dashed line), as well as the photometric position
angle fitted using the synthetic $r$-band image (red dashed line). The
difference between these two position angles is used to calculate the
$\Delta$PA parameter. The bottom right panel shows the $\sigma_{\star}$ map
produced by ppXF, along with the kinematic center fitted using the
$\sigma_{\star}$ map (blue cross) and morphological center fitted using the
synthetic $r$-band image (red cross). The top left and top right panels show
histograms of the distributions of $v_{\mathrm{los}}$ and $\sigma_{\star}$
values in the kinematic maps, respectively, which are used to compute the
statistical moments $\mu_{v}$ and $\mu_{\sigma}$ described in Section 4.1.

##  4 Analysis of the Stellar Kinematic Maps

###  4.1 Extracting Kinematic Parameters

From the kinematic maps of $v_{\mathrm{los}}$ and $\sigma_{\star}$ produced by
ppxf, we extract 12 stellar kinematic measurements, described in detail below.
To calculate these kinematic parameters, we follow the definitions from Nevin
et al. (2019), and use the Radon Transform from Stark et al. (2018) and the
kinemetry package (Krajnović et al., 2006). Some of these kinematic parameters
also use parameters from broadband images, such as the $r$-band effective
radius and ellipticity $\varepsilon$. We thus also use our spectral datacubes
to produce synthetic $r$-band images, and use the StatMorph package
(Rodriguez-Gomez et al., 2019) to compute these parameters.

We first calculate the statistical moments of the line-of-sight velocity
($\mu_{2,v}$, $\mu_{3,v}$, $\mu_{4,v}$) and velocity dispersion maps
($\mu_{1,\sigma}$, $\mu_{2,\sigma}$, $\mu_{3,\sigma}$, $\mu_{4,\sigma}$). The
subscripts 1 through 4 denote the mean, dispersion, skew, and kurtosis of the
1-dimensional distributions. We use the mean line-of-sight velocity
($\mu_{1,v}$) to normalize the line-of-sight velocity maps, and thus it is not
included as a parameter in the LDA. Since we found in Paper I that MBH merger
host galaxies have more prominent bulges in imaging, we expect their mean
velocity dispersion ($\mu_{1,\sigma}$) to be higher.

To compute the kinematic asymmetry parameters $A$ and $A_{2}$, as well as the
kinematic center (used in our other kinematic measurements), we follow the
technique described in Nevin et al. (2019). Specifically, we use the bounded
Absolute Radon Transform from Stark et al. (2018), defined as

$R_{\mathrm{AB}}(r,\theta)=\int_{0}^{R_{e}}|v(x,y)-\langle v(x,y)\rangle|dl,$

where for a given polar coordinate $(r,\theta)$, the line-of-sight velocity
map $v(x,y)$ is integrated over the line segment passing through $(r,\theta)$
and perpendicular to the vector from the origin to $(r,\theta)$, subtracted by
the mean line-of-sight velocity along that line segment $\langle
v(x,y)\rangle$. The line segment is bounded to a width that we set to the
$r$-band effective (i.e., half-light) radius $R_{e}$, as calculated by
StatMorph. Since the Radon Transform is sensitive to the location of the
origin used, we follow the method described by Stark et al. (2018) to find the
kinematic center. We use this kinematic center in all further kinematic
calculations that require a center point. We also use the Radon Transform
calculation to find the asymmetry parameters $A$ and $A_{2}$, as defined by
Nevin et al. (2019), but these asymmetry parameters are excluded from the LDA
because we find that they are not significantly different between the MBH
merger and control galaxy samples.

We use kinemetry (Krajnović et al., 2006) to calculate the difference between
photometric and kinematic position angles $\Delta$PA from the line-of-sight
velocity map, centered around the kinematic center calculated earlier. The
kinematic position angle $\mathrm{PA}_{\mathrm{kin}}$ is calculated by
kinemetry by modelling the line-of-sight velocity map $v(x,y)$ as a series of
elliptical rings, whose finite harmonic expansion is

$v(r,\theta)\approx
A_{0}(r)+\sum_{n=1}^{N}\left[A_{n}(r)\sin(n\theta)+B_{n}(r)\cos(n\theta)\right].$

The best fit model is then found by minimizing $\chi^{2}$ with respect to the
parameters $A_{n}$ and $B_{n}$. We then compare the resulting kinematic
position angle $\mathrm{PA}_{\mathrm{kin}}$ from the best-fit model to the
position angle $\mathrm{PA}_{\mathrm{img}}$ calculated from the $r$-band image
(using StatMorph) to obtain
$\Delta\mathrm{PA}\equiv|\mathrm{PA}_{\mathrm{kin}}-\mathrm{PA}_{\mathrm{img}}|$.
Galaxies with larger position angle differences $\Delta$PA display stronger
stellar kinematic misalignments. We also calculate the residual between the
model and actual galaxy velocity maps ($\mathrm{resid}$) as defined by Nevin
et al. (2019), but ultimately exclude it since it does not significantly
discriminate between the MBH merger and control galaxy samples.

We also calculate the spin parameter $\lambda_{R_{e}}$ (Emsellem et al.,
2007), which measures the specific angular momentum of the galaxy. We define
it as

$\lambda_{R_{e}}=\frac{\sum_{n=1}^{N}F_{n}r_{n}|v_{n}|}{\sum_{n=1}^{N}F_{n}r_{n}\sqrt{v_{n}^{2}+\sigma_{n}^{2}}},$

where $F_{n}$ is the $r$-band flux, $r_{n}$ is the radius to the kinematic
center, $v_{n}$ is the line-of-sight velocity map, $\sigma_{n}$ is the
velocity dispersion map, and $N$ is the number of pixels. We also compute the
morphological ellipticity $\varepsilon$ from the $r$-band image (using
StatMorph) to compute the slow rotator discriminant, defined by Cappellari
(2016) as

$\lambda_{R_{e}}<0.08+\varepsilon/4\;\;\;\mathrm{and}\;\;\;\varepsilon<0.4.$

Slow rotator galaxies (i.e., with small $\lambda_{R_{e}}$) are observed to be
dispersion-dominated, and are more likely to have disturbed kinematics
(Emsellem et al., 2011).

Finally, we calculate morphological properties ($Gini_{\sigma}$,
$M_{20,\sigma}$, $C_{\sigma}$, $A_{\sigma}$, $S_{\sigma}$, $A_{S,\sigma}$)
directly on the $\sigma_{\star}$ maps. We use the same calculations as in
Paper I, which are based on the definitions from Lotz et al. (2004) and Pawlik
et al. (2016).

Figure 2 displays a visualization that summarizes the calculation of a few
kinematic measurements for an example MBH merger host galaxy, including
histograms of $\sigma_{\star}$ and $v_{\mathrm{los}}$ used to compute the
statistical moments, the kinematic and morphological centers, and the position
angles used to calculate $\Delta$PA. Figure 3 displays normalized histograms
of the most significant kinematic parameters for a subsample of MBH merger
host galaxies with high ($\nicefrac{{M_{2}}}{{M_{1}}}\geqslant 0.5$) and high
chirp masses ($M_{\text{chirp}}\geqslant 10^{8.2}M_{\odot}$), in comparison to
the corresponding control sample of galaxies. This subsample of MBH merger
host galaxies with high chirp mass and mass ratio contains a total of 66
spectral datacubes, which includes datacubes generated along different lines
of sight to account for inclination effects. The histograms in Figure 3
suggest that the best individual kinematic predictors for the presence of a
MBH merger are $\lambda_{R_{e}}$, $\Delta$PA, and $\mu_{1,\sigma}$. In
particular, $\lambda_{R_{e}}$ alone can discriminate between MBH merger host
galaxies and control galaxies with $\sim$80% accuracy. In Section 4.2, we will
show that using a linear combination of parameters results in even higher
accuracies of $\gtrsim$85%.

Figure 3: Normalized histograms of the kinematic parameters of the subsample
of MBH mergers host galaxies (red) with high mass ratios
($\nicefrac{{M_{2}}}{{M_{1}}}\geqslant 0.5$) and high chirp masses
($M_{\text{chirp}}\geqslant 10^{8.2}M_{\odot}$), in comparison to a
corresponding mass- and redshift-matched control sample (blue). We only
include kinematic parameters that pass a two-sample Kolmogorov–Smirnov (K-S)
test with a $>$$3\sigma$ significance. Although these histograms show
significant differences between the MBH merger host galaxies and the control
sample, a higher accuracy of discrimination between the two classes can be
obtained by combining multiple kinematic parameters, thus motivating our use
of LDA in Section 4.2.

Figure 4:  Histogram of the LDA predictor (left) for MBH merger hosts (red)
and control galaxies (blue), for the same high mass and high mass ratio
subsample as shown in Figure 3. The LDA predictor finds a linear combination
of stellar kinematic parameters (Equation 1) that maximally discriminates
between these two classes. The decision boundary is indicated by the dashed
black line, with positive LDA predictor value indicating a MBH merger host,
and negative LDA predictor value indicating a control galaxy. The confusion
matrix (right) shows that for this subsample, the LDA predictor can classify
galaxies with an accuracy of 86.0% and a precision of 94.3%.

###  4.2 Linear Discriminant Analysis

We use a LDA to find a linear combination of kinematic parameters that
optimally discriminates between our MBH merger host galaxies and control
sample. Our LDA approach here is the same as that used in Paper I for
morphological parameters; we summarize the most salient points below, and
refer the reader to Paper I for more details. The LDA finds a hyperplane in
parameter space that optimally separates the two classes in the training data,
and the linear equation perpendicular to this hyperplane (the ‘LDA predictor’)
provides a signed distance between each datapoint to the hyperplane that can
be used for classification. The efficacy of this predictor is summarized in a
confusion matrix, in which the sum of the diagonal of the confusion matrix is
the accuracy (the ratio of training data predicted by the LDA correctly). We
use the sklearn Python package (Pedregosa et al., 2011) to perform the LDA.

In addition to the stellar kinematic parameters we discussed in Section 4.1,
we also include interaction terms in the LDA, which are multiplicative
combinations of morphological measurements (e.g.,
$\lambda_{R_{e}}\times\Delta\mathrm{PA}$). The inclusion of these interaction
terms account for cross-correlations between parameters, thus improving the
predictive ability of the LDA (James et al., 2013), and enabling direct
interpretation of the LDA coefficients in Section 5.1.

To prevent any particular parameter from unphysically dominating the LDA, we
first whiten (i.e., normalize) our stellar kinematic parameters. Specifically,
we subtract the mean of the sample from each data point, and then divide by
the standard deviation of the sample.

In training the LDA predictor, we remove stellar kinematic parameters that are
uninformative in discriminating between MBH merger host galaxies and the
control sample. Specifically, we only use parameters that pass a two-sample
Kolmogorov–Smirnov test with a $p$-value of $\leqslant$ 0.05 ($2\sigma$) to
train the LDA predictor, rejecting the hypothesis that the MBH merger host
galaxy sample is drawn from the same distribution as the control sample in the
histograms in Figure 3.

To estimate uncertainties on the classification accuracies from the LDA
predictor, we perform repeated stratified $k$-fold cross validation.
Specifically, we divide the sample of MBH merger host galaxies and control
sample into $k$ bins of equal size, where the first $k-1$ bins are used for
training the LDA, and the last bin is used as a test set to test the resultant
LDA. This is repeated a total of $n=10$ times, with a different set of $k$
bins. The LDA accuracy and uncertainty is thus the median and $1\sigma$ spread
across all runs. In our tests, we find that $k=5$ provides a good balance
between the LDA accuracy and uncertainties, and that the LDA is robust to the
assumed $k$ or $n$.

To reduce the dimensionality of the LDA predictor and improve its
interpretability, we reduce the number of parameters by training the LDA using
forward stepwise selection. We first begin with a LDA model that has only a
single parameter, and choose the parameter that maximizes the total LDA cross
validation accuracy (i.e., minimizes the number of misclassifications). We
then add additional parameters one at a time and retrain the LDA using
$k$-fold cross validation iteratively, until the accuracy no longer increases.
To confirm that our results are not affected by overfitting, we randomly
shuffle the morphological parameters between the MBH merger host galaxies and
and control samples, and verify that the LDA accuracies are consistent with
$\sim$50%.

Figure 4 shows the median-scoring trained LDA predictor for a sample of MBH
merger host galaxies with high chirp masses ($M_{\text{chirp}}\geqslant
10^{8.2}M_{\odot}$) and high mass ratios
($\nicefrac{{M_{2}}}{{M_{1}}}\geqslant 0.5$), as well as its associated
confusion matrix. As discussed below in Section 5.2, we specifically focus on
this sample because it is the host galaxies of these MBH mergers that have the
most unique stellar kinematic signatures, leading to high LDA accuracies. We
find that the LDA discriminates between this MBH merger host galaxy sample and
its corresponding control sample with a median accuracy of 85.7 $\pm$ 4.5% and
median precision of $93.1\pm 6.7$%. The linear equation of the LDA predictor
in Figure 4 is:

$\begin{split}\text{LDA
predictor}=0.51\log\Delta\mathrm{PA}-2.81\lambda_{R_{e}}+0.04\end{split}$ (1)

The cross validation also provides uncertainties on the coefficients in
Equation 1. For this subsample, the 1$\sigma$ uncertainties are $\pm 0.16$ and
$\pm 0.32$, for the $\log\Delta\mathrm{PA}$ and $\lambda_{R_{e}}$ parameters,
respectively. We discuss these results in Section 5.

Figure 5: The accuracy of the LDA predictor as a function of halo mass (upper
left), chirp mass (upper right), MBH mass ratio (lower left), and redshift
(lower right) of the MBH merger in the Romulus25 simulation. For visual
clarity, the halo mass, chirp mass, and redshift plot only contain MBH mergers
with high mass ratio ($\nicefrac{{M_{2}}}{{M_{1}}}\geqslant 0.5$), while the
mass ratio plot only contains MBH merger with high chirp masses of
$M_{\text{chirp}}\geqslant 10^{7.5}M_{\odot}$. The error bars are the
1$\sigma$ spread from cross validation. The accuracy of our stellar kinematics
based approach increases as a function of halo mass, chirp mass, and mass
ratio, and does not have a significant trend with redshift.

##  5 Results and Discussion

###  5.1 LDA Coefficient Analysis

The most distinct kinematic parameters identified by the LDA in Equation 1 are
the spin parameter $\lambda_{R_{e}}$, and difference between kinematic and
photometric position angles $\Delta\mathrm{PA}$. Specifically, the histograms
in Figure 3 show that MBH merger host galaxies have systematically lower
$\lambda_{R_{e}}$ (slower spin) and higher $\Delta\mathrm{PA}$ (stronger
kinematic misalignments), relative to the control sample. Formally, a two-
sample K-S test shows that the distributions of $\lambda_{\mathrm{R_{e}}}$ and
$\Delta\mathrm{PA}$ for MBH merger host galaxies is statistically distinct
from the control sample in Figure 3 at a $\gtrsim$$10\sigma$ and 6.1$\sigma$
level, respectively.

Figure 3 also shows that MBH merger host galaxies have systematically higher
mean velocity dispersion $\mu_{1,\sigma}$, at high significance of 6.4$\sigma$
in a two-sample K-S test. However, $\mu_{1,\sigma}$ was not chosen by the LDA
for inclusion in Equation 1 during training, because it is correlated with
$\lambda_{\mathrm{R_{e}}}$, and thus its inclusion in Equation 1 would not add
additional discriminating power. This latent correlation is observed in IFU
surveys of nearby galaxies, which show that the most massive galaxies have
lower $\lambda_{\mathrm{R_{e}}}$ (Emsellem et al., 2011; Falcón-Barroso et
al., 2019).

###  5.2 Trends in Discriminating Power

We investigate how the accuracy of the LDA predictor in Equation 1 depends on
properties such as the halo virial mass, chirp mass, mass ratio, and redshift.
To do this, we divide both our MBH merger host galaxies and corresponding
control samples into subsamples (e.g., in bins of chirp mass). For each
subsample, we retrain the LDA predictor to calculate its accuracy. The
resulting trends are shown in Figure 5.

Figure 5 shows that the LDA accuracy increases with halo mass, chirp mass, and
mass ratio, and is relatively insensitive to redshift. These trends from
stellar kinematics are similar to those for morphological signatures of these
same MBH merger host galaxies that we previously found Paper I. We discuss the
origin of these trends in Section 5.4 below.

###  5.3 Dependence of Results on Time-Delay

We test whether the accuracy of the LDA predictor changes as a function of
time-delay between the numerical merger and the physical merger of MBHs. For
each binary in our simulation, the time-delay between the numerical merger in
the simulation (at binary separations of $\sim$700 pc) and the physical merger
(when a gravitational wave chirp would be emitted) is highly uncertain, and
can range from 0.1 to 10 Gyrs (e.g., Volonteri et al., 2020; Li et al., 2022).
Instead of carefully estimating this time-delay for each merger in our
samples, we remain agnostic to the detailed physics below the resolution limit
of the simulation, and instead test whether the LDA accuracy changes up to 1
Gyr after the numerical merger. Figure 6 shows the LDA accuracy as a function
of time-delay, for the same high mass ($\gtrsim$108.2 $M_{\odot}$) and high
mass ratio ($\gtrsim$0.5) subsample from Figure 4. The accuracy of our trained
LDA predictor remains high for at least $\sim$1 Gyr after numerical merger
(see Section 5.4).

Figure 6: The accuracy of the LDA predictor as a function of time-delay since
numerical merger in the simulation, for a high mass ratio
($\nicefrac{{M_{2}}}{{M_{1}}}\geqslant 0.5$) and high chirp mass
($M_{chirp}>10^{8}M_{\odot}$) MBH merger subsample and corresponding control
sample. The LDA accuracy remains high for at least $\sim$1 Gyr, suggesting
that the unique stellar kinematic signatures of MBH merger host galaxies are
likely to be permanent. Thus, the accuracy of the stellar kinematic approach
to identifying MBH merger host galaxies will likely be high at the time of
gravitational wave detection.

Figure 7:  Right panel: A comparison of the spin parameter $\lambda_{R_{e}}$
and morphological ellipticity $\varepsilon$, for MBH merger host galaxies and
the control sample. MBH merger host galaxies have systematically slower
rotation (smaller $\lambda_{R_{e}}$) and are more round (smaller
$\varepsilon$). Left panel: A comparison of the difference between kinematic
and photometric position angles $\Delta\mathrm{PA}$ and morphological
ellipticity $\varepsilon$, for the same MBH merger host galaxies and control
sample. MBH merger host galaxies are systematically more kinematically
misaligned (larger $\Delta\mathrm{PA}$). Galaxies with all of these kinematic
properties are thought to have undergone major mergers, which causes kinematic
misalignments and slower rotation. For visual clarity, both panels show the
subsample of high mass ratio MBH merger host
galaxies($\nicefrac{{M_{2}}}{{M_{1}}}\geqslant 0.5$) with high chirp masses
($M_{\text{chirp}}\geqslant 10^{8.2}M_{\odot}$) (red contours), with a
corresponding mass- and redshift-matched control sample (blue contours). The
contours enclose 15.9% and 50% of each sample. The dashed black line is the
slow rotator discriminant from Cappellari (2016), in which galaxies with small
values of $\lambda_{R_{e}}$ and $\varepsilon$ below the line are classified as
slow rotator galaxies.

###  5.4 Interpretation and Comparison to Results from Morphology in Imaging

Our stellar kinematic results support the view that the most distinctive
characteristics of MBH binary and merger host galaxies are permanent features,
rather than transient disturbances stemming from the preceding galaxy mergers.
Previously in Paper I, we used the same MBH merger host galaxy and control
galaxy samples from Romulus25 that we use here to study their morphologies in
synthetic broadband images produced from dust radiative transfer. In that
study, we extracted morphological parameters from these synthetic images and
trained an LDA predictor, similar to our use of kinematic parameters to train
an LDA predictor in Section 4 above. Paper I showed that the accuracy of the
resultant morphology-based LDA predictor does not decrease as a function of
time-delay, which suggests that the unique morphological signatures of MBH
merger host galaxies are permanent, instead of transient features such as
short-lived morphological disturbances stemming from the preceding galaxy
merger. This is consistent with our results shown in Figure 6, which
demonstrates that the accuracy of our LDA predictor trained on the stellar
kinematic parameters of these galaxies also does not decrease over time for at
least 1 Gyr after the numerical merger. Thus, both morphology and stellar
kinematics indicate that MBH merger host galaxies have unique and long-lived
characteristics.

The stellar kinematics of our simulated MBH merger host galaxies also suggest
that their most distinctive characteristics are permanent kinematic properties
associated with major mergers of massive galaxies, consistent with our
morphological results from Paper I. The two most informative stellar kinematic
parameters in our LDA predictor for MBH merger host galaxies are a
systematically slower spin $\lambda_{R_{e}}$, and a larger difference between
the kinematic and photometric position angles $\Delta\mathrm{PA}$, as
discussed in Section 5.1 above and shown in Figure 3. Galaxy IFU surveys have
revealed galaxies with both slow spin (small $\lambda_{R_{e}}$; Emsellem et
al., 2011; Graham et al., 2018; Ene et al., 2020) and strong misalignment
between their kinematic and photometric position angles (large
$\Delta\mathrm{PA}$; Krajnović et al., 2011; Cortese et al., 2016; Ene et al.,
2018). These properties are more clearly shown for our simulated galaxies in
Figure 7, which compares our MBH merger host galaxies to the control sample in
the $\lambda_{R_{e}}$ – $\varepsilon$ and $\Delta\mathrm{PA}$ – $\varepsilon$
planes, as commonly plotted for IFU observations. Numerical simulations show
that these galaxy properties result from major galaxy mergers, which cause
irregular rotation, and also decrease the net rotation of galaxies (e.g.,
Jesseit et al., 2009; Bois et al., 2011; Khochfar et al., 2011; Naab et al.,
2014; Penoyre et al., 2017; Lagos et al., 2018). Furthermore, we also find in
Figure 3 that MBH merger hosts also have systematically higher mean stellar
velocity dispersion $\sigma_{*}$, although this parameter was not selected by
the LDA for inclusion in Equation 1 because it is correlated with
$\lambda_{\mathrm{R_{e}}}$, and thus its predictive information is redundant.
Higher $\sigma_{*}$ is indicative of more prominent bulges, which are also a
natural consequence of major galaxy mergers. Overall, our interpretation that
the unique stellar kinematic signatures of MBH merger host galaxies stem from
major mergers of massive galaxies is highly consistent with our previous
results in Paper I, which showed that the unique morphological signatures of
MBH merger host galaxies is the presence of a strong classical bulge, built
through major galaxy mergers.

Our interpretation also naturally explains the trends in the LDA accuracy that
we observe as a function of chirp mass and mass ratio. A key result from this
study is that the stellar kinematics of MBH merger host galaxies are only
distinct for MBH binaries and mergers with high chirp masses and high mass
ratios. These trends are clearly observed in Figure 5, and naturally arise
from a scenario in which the unique kinematic signatures owe to major mergers
of massive galaxies. Numerical simulations demonstrate that galaxies with slow
spin and large kinematic misalignments arise from major mergers of massive
gas-poor galaxies (e.g., Jesseit et al., 2009; Bois et al., 2011; Naab et al.,
2014). Since it is specifically these galaxy mergers that produce MBH binaries
and mergers with high chirp mass and mass ratio, this naturally causes the LDA
accuracy to increase with chirp mass and mass ratio, as observed in Figure 5.
Since chirp mass is correlated with halo mass, it is thus also expected that
the LDA accuracy is high for the most massive halos in Figure 5. Although the
LDA accuracy from our previous results from host galaxy morphology decreased
mildly with redshift, our LDA accuracy from stellar kinematics in Figure 5
does not appear to display strong trends with redshift. This may owe to the
fact that bulges are more difficult to resolve and detect in the morphology of
higher-redshift galaxies, while their stellar kinematic signatures in IFU
spectroscopy remain relatively distinct.

A consistent picture thus emerges from our study of the unique stellar
kinematics of MBH merger host galaxies, in broad agreement with our previous
results from morphology in Paper I, as well as other studies. As massive
galaxies undergo major mergers, they form prominent classical bulges that can
be observed in their morphologies in broadband imaging (e.g., through the
bulge discriminant statistic $F(Gini,M_{20})$ described in Paper I) and
directly through larger mean stellar velocity dispersions $\sigma_{*}$ in IFU
spectra. The history of major mergers in these galaxies also cause their
stellar kinematics to display misalignments (large $\Delta\mathrm{PA}$) and
slower rotation (lower $\lambda_{R_{e}}$). The MBHs in these massive merging
galaxies will then go on to form binaries with high chirp masses and mass
ratios, and eventually merge due to gravitational wave losses. It is thus
specifically in MBH binaries and mergers with high chirp masses and mass
ratios that both the morphological and stellar kinematic signatures of the
host galaxies are most unique. These signatures can thus be used to identify
the host galaxies of MBH binaries and mergers detected in gravitational waves
with high accuracies. Lastly, our results are also broadly consistent with
other studies that use cosmological simulations to investigate the host
galaxies of MBH binaries detectable by PTA experiments (which preferentially
have high mass ratios and high chirp masses). In particular, Saeedzadeh et al.
(2024) show that PTA-detectable MBH binaries in Romulus25 are massive early-
type galaxies with low star formation rates residing at the centers of galaxy
groups, while Cella et al. (2024) show that PTA-detectable MBH binaries in the
Illustris simulations have redder colors and higher metallicities. These
properties are all expected from major mergers of massive galaxies, similar to
the stellar kinematic and morphological signatures that we identify.

###  5.5 Practical Considerations for Observations

We caution that our trained LDA predictor may not be directly applicable to
observations from IFU instruments with very different specifications than
assumed here. In producing our synthetic spectral datacubes, we assumed
specifications such as medium spectral resolution of $R=4000$, pixel scales of
$0\farcs 1$/pixel, and Gaussian $0\farcs 1$ FWHM PSFs. The LDA predictor in
Equation 1 can thus be directly applied to observed spectral datacubes
obtained by similar instruments. However, our results may no longer hold true
when applied to spectral datacubes from telescopes and instruments with vastly
different specifications. We defer a more thorough test of these effects to
future work.

Despite our caution above, we suggest that in practice, our results may
actually be more broadly applicable even to more modest instruments and
observing conditions. The fact that both morphological and stellar kinematic
host galaxy signatures are only distinct for MBH binaries and mergers with
high chirp masses and mass ratios suggest that these approaches will be most
useful for gravitational wave detections of individual MBH binaries by PTA
experiments. These MBH binaries will be in relatively nearby galaxies, at
significantly lower redshifts than our sample of simulated host galaxies from
Romulus25. Thus, current IFU instruments with more modest capabilities may
already be able to achieve results comparable to our findings, by resolving
physical scales in each galaxy that are similar to our assumptions. For
example, at the median redshift ($z\sim 1$) of simulated galaxies in our
sample, the assumed $0\farcs 1$ pixels resolve physical scales of
approximately 0.8 kpc, while this same physical scale can be resolved with
$2\farcs$0 pixels for a nearby galaxy at 90 Mpc. This suggests that our
trained LDA predictor could be directly applied to IFU spectral datacubes of
nearby galaxies that have significantly worse angular resolutions, provided
that they probe similar physical scales in the host galaxy targets.
Furthermore, archival IFU spectra of a significant percent of these nearby
galaxies already exist from large IFU surveys such as SAURON (Bacon et al.,
2001), ATLAS3D (Cappellari et al., 2011), CALIFA (Sánchez et al., 2012),
MASSIVE (Ma et al., 2014), SAMI (Bryant et al., 2015), and MANGA (Bundy et
al., 2015). Thus, efforts to identify MBH merger host galaxies in
gravitational wave localization volumes using our stellar kinematics approach
can leverage existing archival data, minimizing the need to obtain expensive
new IFU spectroscopy.

The main disadvantage of using a stellar kinematics approach to identifying
MBH binary and merger host galaxies is that it will only be useful if the
number of candidate host galaxies lying in typical gravitational wave error
volumes is small. For a gravitational wave error volume, selection cuts based
on redshift and stellar mass (assuming empirical scaling relations) must be
used first to limit the number of candidate galaxies, before additional
approaches such as galaxy stellar kinematics or morphology can be applied.
Despite the high accuracies of $\gtrsim$85% we find using our stellar
kinematics approach, it will only be useful if there are a small number
($\lesssim$10) of candidate galaxies remaining in the gravitational wave error
volume. Such a scenario is possible for PTA detections of individual binaries
in optimistic scenarios (Petrov et al., 2024), but in practice a combination
of several approaches will be required to identify the host galaxy.

###  5.6 Caveats

Although we tested the accuracy of our LDA predictor in Figure 6 up to 1 Gyr
after each numerical merger, the time-delay to physical merger can be up to
several Gyrs. For each MBH merger, the time-delay between the numerical merger
and physical merger is highly uncertain, and is dependent on the poorly-
understood dynamical friction and stellar scattering below the resolution
limit of the simulation. More detailed calculations have suggested that this
time-delay can range between 0.1–10 Gyr (Volonteri et al., 2020; Li et al.,
2022). In Section 5.3, we remain agnostic to the exact time-delay of each MBH
merger, and instead show in Figure 6 that the LDA accuracy does not decrease
for at least $\sim$1 Gyr. This choice was primarily driven by the
computational costs of performing radiative transfer simulations. The constant
accuracy of our LDA predictor over $\sim$1 Gyr after numerical merger is
consistent with (and an expected consequence of) our result that the most
unique characteristics of MBH merger host galaxies are slower spin, stronger
kinematic misalignment, and larger mean velocity dispersions, since these are
relatively permanent features of galaxies. Nevertheless, future work should
explicitly test these results beyond 1 Gyr, to understand whether these unique
signatures eventually fade.

##  6 Conclusions

We demonstrated that the host galaxies of MBH binaries and mergers detected in
gravitational wave experiments have unique stellar kinematics in IFU spectra.
Using a sample of simulated galaxies hosting MBH mergers selected from a
cosmological simulation, we performed stellar population synthesis and dust
radiative transfer simulations to produce optical 3D spectral datacubes that
are similar to data from IFU spectrographs. We performed full-spectrum fitting
on these spectral datacubes to generate kinematic maps of $\sigma_{\star}$ and
$v_{\mathrm{los}}$ for each host galaxy. We then extracted a variety of
stellar kinematic parameters from these maps, and compared their properties to
a mass- and redshift-matched control sample of simulated galaxies that do not
host MBH mergers. Specifically, we trained a LDA predictor on the stellar
kinematic parameters, to produce a linear equation that optimally
distinguishes MBH merger host galaxies from the control sample. Our main
findings are:

  1. 1.

The stellar kinematic properties of MBH merger host galaxies are unique,
suggesting that IFU spectra can be used to identify the host galaxies of MBH
binaries and mergers detected in gravitational waves. The accuracy of this
approach increases with chirp mass and mass ratio; for mergers with
particularly high chirp masses ($\gtrsim$108.2 $M_{\odot}$) and high mass
ratios ($\gtrsim$0.5), the accuracies reach $\gtrsim$85%. These accuracies are
slightly higher than approaches using galaxy morphology in broadband imaging,
with similar trends. However, our results also suggest that in practice, both
stellar kinematic and morphology-based approaches to identifying host galaxy
counterparts will only be useful for the most massive MBH binaries and
mergers, such as those that will be detected in nearby galaxies by PTA
experiments.

  2. 2.

The most unique stellar kinematic properties of MBH merger host galaxies are
systematically lower specific angular momentum $\lambda_{R_{e}}$, larger
difference between kinematic and photometric position angles
$\Delta\mathrm{PA}$, and higher stellar velocity dispersions $\sigma_{*}$.
This implies that that MBH merger host galaxies have systematically slower
rotation, more irregular stellar kinematics, and stronger bulges.
Observationally, these properties are commonly associated with massive early-
type galaxies that have undergone major mergers, which would naturally host
MBH binaries and mergers. Our results are thus consistent with a scenario in
which MBH binaries and mergers with high masses and mass ratios form in major
mergers of massive galaxies, which have unique morphologies and stellar
kinematics as probed by broadband imaging and IFU spectroscopy.

The high accuracies of $\gtrsim$85% achieved by our stellar kinematics
approach suggests that IFU spectroscopy can play a role in telescope follow-up
of MBH binaries and mergers detected in gravitational waves. These
$\gtrsim$85% accuracies are slightly higher than the $\gtrsim$80% accuracies
we find in Paper I using host galaxy morphology in broadband imaging. However,
we emphasize that we find these high accuracies specifically for MBH binaries
and mergers with high chirp masses and mass ratios, and thus our stellar
kinematics approach may only be relevant for PTA detections of massive MBH
binaries. Nevertheless, these PTA detections are limited to nearby massive
galaxies, where IFU spectroscopy of galaxies in the localization region will
be most feasible, and a plethora of archival IFU spectra already exists.

Our results suggest that prospects are bright for multi-messenger science with
MBH binaries and mergers that will be detected in low-frequency gravitational
waves. When a MBH binary or merger with high chirp mass and mass ratio is
detected, initial redshift and stellar mass selection cuts will reduce the
candidate galaxies in the gravitational wave error volume by orders of
magnitude. If only a small number ($\lesssim$10) of candidates remain,
additional broadband imaging and IFU spectroscopy may uniquely identify the
host galaxy based on morphology and stellar kinematics, respectively.
Crucially, these approaches will be effective even if there are no other
electromagnetic signatures, such as if there is no active accretion onto the
MBH binary or merger, or if there is heavy dust obscuration. By combining
insights from a variety of host galaxy identification techniques, the promise
of the diverse array of science goals based on multi-messenger observations of
MBH binaries and mergers can be realized.

We thank Thomas R. Quinn, Jessie Runnoe, Stephen Taylor, Kelly Holley-
Bockelmann, and Maria Charisi for insightful discussions. This work made
extensive use of the Cedar cluster of Digital Research Alliance of Canada at
Simon Fraser University. Romulus25 is part of the Blue Waters sustained-
petascale computing project, which is supported by the National Science
Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois.
Blue Waters is a joint effort of the University of Illinois at Urbana-
Champaign and its National Center for Supercomputing Applications. J.B.
acknowledges support from the NSERC Summer Undergraduate Research Award
program. J.J.R. and D.H. acknowledge support from the Canada Research Chairs
program, the NSERC Discovery Grant program, and the FRQNT Nouveaux Chercheurs
Grant program. J.J.R. acknowledges funding from the Canada Foundation for
Innovation, and the Québec Ministère de l’Économie et de l’Innovation.

## References

  * Agazie et al. (2023a) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8, doi: 10.3847/2041-8213/acdac6
  * Agazie et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.16220, doi: 10.48550/arXiv.2306.16220
  * Agazie et al. (2023c) —. 2023c, ApJ, 951, L50, doi: 10.3847/2041-8213/ace18a
  * Agazie et al. (2023d) —. 2023d, ApJ, 956, L3, doi: 10.3847/2041-8213/acf4fd
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786, doi: 10.48550/arXiv.1702.00786
  * Armitage & Natarajan (2002) Armitage, P. J., & Natarajan, P. 2002, ApJ, 567, L9, doi: 10.1086/339770
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bacon et al. (2001) Bacon, R., Copin, Y., Monnet, G., et al. 2001, MNRAS, 326, 23, doi: 10.1046/j.1365-8711.2001.04612.x
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 773508, doi: 10.1117/12.856027
  * Baes et al. (2011) Baes, M., Verstappen, J., De Looze, I., et al. 2011, ApJS, 196, 22, doi: 10.1088/0067-0049/196/2/22
  * Banks et al. (2022) Banks, S., Lee, K., Azimi, N., et al. 2022, MNRAS, 512, 6007, doi: 10.1093/mnras/stac831
  * Bardati et al. (2024) Bardati, J., Ruan, J. J., Haggard, D., & Tremmel, M. 2024, The Astrophysical Journal, 961, 34, doi: 10.3847/1538-4357/ad055a
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307, doi: 10.1038/287307a0
  * Berczik et al. (2006) Berczik, P., Merritt, D., Spurzem, R., & Bischof, H.-P. 2006, ApJ, 642, L21, doi: 10.1086/504426
  * Bogdanović et al. (2022) Bogdanović, T., Miller, M. C., & Blecha, L. 2022, Living Reviews in Relativity, 25, 3, doi: 10.1007/s41114-022-00037-8
  * Bois et al. (2011) Bois, M., Emsellem, E., Bournaud, F., et al. 2011, MNRAS, 416, 1654, doi: 10.1111/j.1365-2966.2011.19113.x
  * Bryant et al. (2015) Bryant, J. J., Owers, M. S., Robotham, A. S. G., et al. 2015, MNRAS, 447, 2857, doi: 10.1093/mnras/stu2635
  * Bundy et al. (2015) Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7, doi: 10.1088/0004-637X/798/1/7
  * Camps & Baes (2015) Camps, P., & Baes, M. 2015, Astronomy and Computing, 9, 20, doi: 10.1016/j.ascom.2014.10.004
  * Camps & Baes (2018) —. 2018, ApJ, 861, 80, doi: 10.3847/1538-4357/aac824
  * Camps & Baes (2020) —. 2020, Astronomy and Computing, 31, 100381, doi: 10.1016/j.ascom.2020.100381
  * Cappellari (2016) Cappellari, M. 2016, ARA&A, 54, 597, doi: 10.1146/annurev-astro-082214-122432
  * Cappellari (2017) —. 2017, MNRAS, 466, 798, doi: 10.1093/mnras/stw3020
  * Cappellari (2023) —. 2023, MNRAS, 526, 3273, doi: 10.1093/mnras/stad2597
  * Cappellari & Copin (2003) Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345, doi: 10.1046/j.1365-8711.2003.06541.x
  * Cappellari et al. (2011) Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 413, 813, doi: 10.1111/j.1365-2966.2010.18174.x
  * Cella et al. (2024) Cella, K., Taylor, S. R., & Kelley, L. Z. 2024, arXiv e-prints, arXiv:2407.01659, doi: 10.48550/arXiv.2407.01659
  * Chapon et al. (2013) Chapon, D., Mayer, L., & Teyssier, R. 2013, MNRAS, 429, 3114, doi: 10.1093/mnras/sts568
  * Charisi et al. (2022) Charisi, M., Taylor, S. R., Runnoe, J., Bogdanovic, T., & Trump, J. R. 2022, MNRAS, 510, 5929, doi: 10.1093/mnras/stab3713
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Colpi et al. (1999) Colpi, M., Mayer, L., & Governato, F. 1999, ApJ, 525, 720, doi: 10.1086/307952
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833
  * Corrales et al. (2010) Corrales, L. R., Haiman, Z., & MacFadyen, A. 2010, MNRAS, 404, 947, doi: 10.1111/j.1365-2966.2010.16324.x
  * Cortese et al. (2016) Cortese, L., Fogarty, L. M. R., Bekki, K., et al. 2016, MNRAS, 463, 170, doi: 10.1093/mnras/stw1891
  * Cuadra et al. (2009) Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, M. C. 2009, MNRAS, 393, 1423, doi: 10.1111/j.1365-2966.2008.14147.x
  * DeGraf et al. (2021) DeGraf, C., Sijacki, D., Matteo, T. D., et al. 2021, Monthly Notices of the Royal Astronomical Society, 503, 3629, doi: 10.1093/mnras/stab721
  * Dong-Páez et al. (2023) Dong-Páez, C. A., Volonteri, M., Beckmann, R. S., et al. 2023, A&A, 676, A2, doi: 10.1051/0004-6361/202346435
  * D’Orazio & Charisi (2023) D’Orazio, D. J., & Charisi, M. 2023, arXiv e-prints, arXiv:2310.16896, doi: 10.48550/arXiv.2310.16896
  * D’Orazio et al. (2015) D’Orazio, D. J., Haiman, Z., & Schiminovich, D. 2015, Nature, 525, 351, doi: 10.1038/nature15262
  * Dotti et al. (2007) Dotti, M., Colpi, M., Haardt, F., & Mayer, L. 2007, MNRAS, 379, 956, doi: 10.1111/j.1365-2966.2007.12010.x
  * Duffell et al. (2020) Duffell, P. C., D’Orazio, D., Derdzinski, A., et al. 2020, ApJ, 901, 25, doi: 10.3847/1538-4357/abab95
  * Emsellem et al. (2007) Emsellem, E., Cappellari, M., Krajnović, D., et al. 2007, MNRAS, 379, 401, doi: 10.1111/j.1365-2966.2007.11752.x
  * Emsellem et al. (2011) —. 2011, MNRAS, 414, 888, doi: 10.1111/j.1365-2966.2011.18496.x
  * Ene et al. (2020) Ene, I., Ma, C.-P., Walsh, J. L., et al. 2020, ApJ, 891, 65, doi: 10.3847/1538-4357/ab7016
  * Ene et al. (2018) Ene, I., Ma, C.-P., Veale, M., et al. 2018, MNRAS, 479, 2810, doi: 10.1093/mnras/sty1649
  * Escala et al. (2004) Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2004, ApJ, 607, 765, doi: 10.1086/386278
  * Falcón-Barroso et al. (2019) Falcón-Barroso, J., van de Ven, G., Lyubenova, M., et al. 2019, A&A, 632, A59, doi: 10.1051/0004-6361/201936413
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9, doi: 10.1086/312838
  * Fontecilla et al. (2019) Fontecilla, C., Haiman, Z., & Cuadra, J. 2019, MNRAS, 482, 4383, doi: 10.1093/mnras/sty2972
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13, doi: 10.1086/312840
  * Goldstein et al. (2019) Goldstein, J. M., Sesana, A., Holgado, A. M., & Veitch, J. 2019, MNRAS, 485, 248, doi: 10.1093/mnras/stz420
  * Graham et al. (2015) Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74, doi: 10.1038/nature14143
  * Graham et al. (2018) Graham, M. T., Cappellari, M., Li, H., et al. 2018, MNRAS, 477, 4711, doi: 10.1093/mnras/sty504
  * Holz & Hughes (2005) Holz, D. E., & Hughes, S. A. 2005, ApJ, 629, 15, doi: 10.1086/431341
  * Izquierdo-Villalba et al. (2023a) Izquierdo-Villalba, D., Colpi, M., Volonteri, M., et al. 2023a, A&A, 677, A123, doi: 10.1051/0004-6361/202347008
  * Izquierdo-Villalba et al. (2023b) Izquierdo-Villalba, D., Sesana, A., & Colpi, M. 2023b, MNRAS, 519, 2083, doi: 10.1093/mnras/stac3677
  * James et al. (2013) James, G., Witten, D., Hastie, T., & Tibshirani, R. 2013, An Introduction to Statistical Learning (Springer New York), doi: 10.1007/978-1-4614-7138-7
  * Jenet et al. (2004) Jenet, F. A., Lommen, A., Larson, S. L., & Wen, L. 2004, ApJ, 606, 799, doi: 10.1086/383020
  * Jesseit et al. (2009) Jesseit, R., Cappellari, M., Naab, T., Emsellem, E., & Burkert, A. 2009, MNRAS, 397, 1202, doi: 10.1111/j.1365-2966.2009.14984.x
  * Jones et al. (2017) Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46, doi: 10.1051/0004-6361/201630225
  * Jung et al. (2022) Jung, S. L., Rennehan, D., Saeedzadeh, V., et al. 2022, MNRAS, 515, 22, doi: 10.1093/mnras/stac1622
  * Kelley et al. (2018) Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2018, MNRAS, 477, 964, doi: 10.1093/mnras/sty689
  * Kelley et al. (2019) Kelley, L. Z., Haiman, Z., Sesana, A., & Hernquist, L. 2019, MNRAS, 485, 1579, doi: 10.1093/mnras/stz150
  * Khochfar et al. (2011) Khochfar, S., Emsellem, E., Serra, P., et al. 2011, MNRAS, 417, 845, doi: 10.1111/j.1365-2966.2011.19486.x
  * Knollmann & Knebe (2009) Knollmann, S. R., & Knebe, A. 2009, The Astrophysical Journal Supplement Series, 182, 608, doi: 10.1088/0067-0049/182/2/608
  * Koss et al. (2018) Koss, M. J., Blecha, L., Bernhard, P., et al. 2018, Nature, 563, 214, doi: 10.1038/s41586-018-0652-7
  * Krajnović et al. (2006) Krajnović, D., Cappellari, M., de Zeeuw, P. T., & Copin, Y. 2006, MNRAS, 366, 787, doi: 10.1111/j.1365-2966.2005.09902.x
  * Krajnović et al. (2011) Krajnović, D., Emsellem, E., Cappellari, M., et al. 2011, MNRAS, 414, 2923, doi: 10.1111/j.1365-2966.2011.18560.x
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Lagos et al. (2018) Lagos, C. d. P., Stevens, A. R. H., Bower, R. G., et al. 2018, MNRAS, 473, 4956, doi: 10.1093/mnras/stx2667
  * Li et al. (2022) Li, K., Bogdanović, T., Ballantyne, D. R., & Bonetti, M. 2022, ApJ, 933, 104, doi: 10.3847/1538-4357/ac74b5
  * Lops et al. (2023) Lops, G., Izquierdo-Villalba, D., Colpi, M., et al. 2023, MNRAS, 519, 5962, doi: 10.1093/mnras/stad058
  * Lotz et al. (2004) Lotz, J. M., Primack, J., & Madau, P. 2004, AJ, 128, 163, doi: 10.1086/421849
  * Ma et al. (2014) Ma, C.-P., Greene, J. E., McConnell, N., et al. 2014, ApJ, 795, 158, doi: 10.1088/0004-637X/795/2/158
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285, doi: 10.1086/300353
  * Mayer et al. (2007) Mayer, L., Kazantzidis, S., Madau, P., et al. 2007, Science, 316, 1874, doi: 10.1126/science.1141858
  * Megevand et al. (2009) Megevand, M., Anderson, M., Frank, J., et al. 2009, Phys. Rev. D, 80, 024012, doi: 10.1103/PhysRevD.80.024012
  * Menon et al. (2015) Menon, H., Wesolowski, L., Zheng, G., et al. 2015, Computational Astrophysics and Cosmology, 2, doi: 10.1186/s40668-015-0007-9
  * Merritt & Poon (2004) Merritt, D., & Poon, M. Y. 2004, ApJ, 606, 788, doi: 10.1086/382497
  * Milosavljević & Phinney (2005) Milosavljević, M., & Phinney, E. S. 2005, ApJ, 622, L93, doi: 10.1086/429618
  * Mingarelli et al. (2017) Mingarelli, C. M. F., Lazio, T. J. W., Sesana, A., et al. 2017, Nature Astronomy, 1, 886, doi: 10.1038/s41550-017-0299-6
  * Naab et al. (2014) Naab, T., Oser, L., Emsellem, E., et al. 2014, MNRAS, 444, 3357, doi: 10.1093/mnras/stt1919
  * Nevin et al. (2019) Nevin, R., Blecha, L., Comerford, J., & Greene, J. 2019, The Astrophysical Journal, 872, 76, doi: 10.3847/1538-4357/aafd34
  * Nevin et al. (2021) Nevin, R., Blecha, L., Comerford, J., et al. 2021, ApJ, 912, 45, doi: 10.3847/1538-4357/abe2a9
  * Noble et al. (2012) Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51, doi: 10.1088/0004-637X/755/1/51
  * O’Neill et al. (2009) O’Neill, S. M., Miller, M. C., Bogdanović, T., Reynolds, C. S., & Schnittman, J. D. 2009, ApJ, 700, 859, doi: 10.1088/0004-637X/700/1/859
  * Pawlik et al. (2016) Pawlik, M. M., Wild, V., Walcher, C. J., et al. 2016, Monthly Notices of the Royal Astronomical Society, 456, 3032, doi: 10.1093/mnras/stv2878
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Penoyre et al. (2017) Penoyre, Z., Moster, B. P., Sijacki, D., & Genel, S. 2017, MNRAS, 468, 3883, doi: 10.1093/mnras/stx762
  * Petrov et al. (2024) Petrov, P., Taylor, S. R., Charisi, M., & Ma, C.-P. 2024, arXiv e-prints, arXiv:2406.04409, doi: 10.48550/arXiv.2406.04409
  * Pfister et al. (2017) Pfister, H., Lupi, A., Capelo, P. R., et al. 2017, MNRAS, 471, 3646, doi: 10.1093/mnras/stx1853
  * Planck Collaboration (2016) Planck Collaboration. 2016, Astronomy & Astrophysics, 594, A13, doi: 10.1051/0004-6361/201525830
  * Pontzen et al. (2013) Pontzen, A., Roškar, R., Stinson, G. S., et al. 2013, pynbody: Astrophysics Simulation Analysis for Python
  * Pontzen & Tremmel (2018) Pontzen, A., & Tremmel, M. 2018, ApJS, 237, 23, doi: 10.3847/1538-4365/aac832
  * Quinlan (1996) Quinlan, G. D. 1996, New A, 1, 35, doi: 10.1016/S1384-1076(96)00003-6
  * Ravi et al. (2015) Ravi, V., Wyithe, J. S. B., Shannon, R. M., & Hobbs, G. 2015, MNRAS, 447, 2772, doi: 10.1093/mnras/stu2659
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez, V., Snyder, G. F., Lotz, J. M., et al. 2019, MNRAS, 483, 4140, doi: 10.1093/mnras/sty3345
  * Rossi et al. (2010) Rossi, E. M., Lodato, G., Armitage, P. J., Pringle, J. E., & King, A. R. 2010, MNRAS, 401, 2021, doi: 10.1111/j.1365-2966.2009.15802.x
  * Saeedzadeh et al. (2024) Saeedzadeh, V., Mukherjee, S., Babul, A., Tremmel, M., & Quinn, T. R. 2024, MNRAS, 529, 4295, doi: 10.1093/mnras/stae513
  * Sánchez et al. (2012) Sánchez, S. F., Kennicutt, R. C., Gil de Paz, A., et al. 2012, A&A, 538, A8, doi: 10.1051/0004-6361/201117353
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703, doi: 10.1111/j.1365-2966.2006.10699.x
  * Schutz (1986) Schutz, B. F. 1986, Nature, 323, 310, doi: 10.1038/323310a0
  * Sesana et al. (2004) Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004, ApJ, 611, 623, doi: 10.1086/422185
  * Sesana et al. (2005) —. 2005, ApJ, 623, 23, doi: 10.1086/428492
  * Sesana & Khan (2015) Sesana, A., & Khan, F. M. 2015, MNRAS, 454, L66, doi: 10.1093/mnrasl/slv131
  * Sesana et al. (2009) Sesana, A., Vecchio, A., & Volonteri, M. 2009, MNRAS, 394, 2255, doi: 10.1111/j.1365-2966.2009.14499.x
  * Stark et al. (2018) Stark, D. V., Bundy, K. A., Westfall, K., et al. 2018, MNRAS, 480, 2217, doi: 10.1093/mnras/sty1991
  * Thatte et al. (2020) Thatte, N. A., Bryson, I., Clarke, F., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Ground-based and Airborne Instrumentation for Astronomy VIII, ed. C. J. Evans, J. J. Bryant, & K. Motohara, 114471W, doi: 10.1117/12.2562144
  * Tremmel et al. (2015) Tremmel, M., Governato, F., Volonteri, M., & Quinn, T. R. 2015, MNRAS, 451, 1868, doi: 10.1093/mnras/stv1060
  * Tremmel et al. (2017) Tremmel, M., Karcher, M., Governato, F., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 1121, doi: 10.1093/mnras/stx1160
  * Vasiliev et al. (2015) Vasiliev, E., Antonini, F., & Merritt, D. 2015, ApJ, 810, 49, doi: 10.1088/0004-637X/810/1/49
  * Vazdekis et al. (2010) Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2010, MNRAS, 404, 1639, doi: 10.1111/j.1365-2966.2010.16407.x
  * Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267, doi: 10.1093/mnras/stw347
  * Volonteri et al. (2003) Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559, doi: 10.1086/344675
  * Volonteri & Natarajan (2009) Volonteri, M., & Natarajan, P. 2009, MNRAS, 400, 1911, doi: 10.1111/j.1365-2966.2009.15577.x
  * Volonteri et al. (2020) Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219, doi: 10.1093/mnras/staa2384
  * Wyithe & Loeb (2003) Wyithe, J. S. B., & Loeb, A. 2003, ApJ, 590, 691, doi: 10.1086/375187
  * Yu (2002) Yu, Q. 2002, MNRAS, 331, 935, doi: 10.1046/j.1365-8711.2002.05242.x
  * Zieleniewski et al. (2015) Zieleniewski, S., Thatte, N., Kendrew, S., et al. 2015, MNRAS, 453, 3754, doi: 10.1093/mnras/stv1860