# The supermassive black hole merger driven evolution of high-redshift red
nuggets into present-day cored early-type galaxies

Antti Rantala1, Alexander Rawlings2, Thorsten Naab1, Jens Thomas3, Peter H.
Johansson2
1Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748,
Garching, Germany
2Department of Physics, University of Helsinki, P.O. Box 64 (Gustaf
Hällströmin katu 2), FI-00014, University of Helsinki, Finland
3 Max-Planck-Institut für Extraterrestriche Physik, Giessenbach-Str. 1,
D-85748, Garching, Germany E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Very compact ($R_{\mathrm{e}}\lesssim 1$ kpc) massive quiescent galaxies (red
nuggets) are more abundant in the high-redshift Universe ($z\sim 2$–$3$) than
today. Their size evolution can be explained by collisionless dynamical
processes in galaxy mergers which, however, fail to reproduce the diffuse low-
density central cores in the local massive early-type galaxies (ETGs). We use
sequences of major and minor merger N-body simulations starting with compact
spherical and disk-like progenitor models to investigate the impact of
supermassive black holes (SMBHs) on the evolution of the galaxies. With the
KETJU code we accurately follow the collisional interaction of the SMBHs with
the nearby stellar population and the collisionless evolution of the galaxies
and their dark matter halos. We show that only models including SMBHs can
simultaneously explain the formation of low-density cores up to sizes of
$R_{\mathrm{b}}\sim 1.3$ kpc with mass deficits in the observed range and the
rapid half-mass size evolution. In addition, the orbital structure in the core
region (tangentially biased orbits) is consistent with observation-based
results for local cored ETGs. The displacement of stars by the SMBHs boost the
half-mass size evolution by up to a factor of two and even fast rotating
progenitors (compact quiescent disks) lose their rotational support after
$6$–$8$ mergers. We conclude that the presence of SMBHs is required for merger
driven evolution models of high redshift red nuggets into local ETGs.

###### keywords:

gravitation – methods: numerical – galaxies: elliptical and lenticular, cD –
galaxies: evolution – galaxies: kinematics and dynamics – quasars:
supermassive black holes

††pubyear: 2024††pagerange: The supermassive black hole merger driven
evolution of high-redshift red nuggets into present-day cored early-type
galaxies–C

##  1 Introduction

Cosmological simulations supported by observations have revealed that the
formation of massive early-type galaxy formation can be divided into two
distinct phases (e.g Naab et al. 2009; Feldmann et al. 2010; Oser et al. 2010,
2012; Johansson et al. 2012; Hilz et al. 2012; Hilz et al. 2013; Wellons et
al. 2015; Rodriguez-Gomez et al. 2016; Qu et al. 2017; Naab & Ostriker 2017;
Ene et al. 2020; Cannarozzo et al. 2023). At early times ($z\gtrsim 1.5$)
galaxies grow rapidly by in-situ star formation fueled by cold gas flows and
accretion of multiple star-bursting progenitors. The latter extended phase,
characterized by mass and especially size growth of the galaxies is primarily
driven by the accretion of stars originally formed in smaller satellite
galaxies. In addition, simple arguments based on the virial theorem can be
used to show that multiple minor mergers will produce a more extended galaxy
compared to accreting the same stellar mass in few major mergers (Naab et al.,
2009; Bezanson et al., 2009).

The effective radii $R_{\mathrm{e}}$ of massive early-type galaxies (ETGs)
have been observed to grow by a factor of $\sim 4$–$7$ since $z\sim 2$ (e.g.
Daddi et al. 2005; Trujillo et al. 2007; van Dokkum et al. 2008; Buitrago et
al. 2008; Franx et al. 2008; Bezanson et al. 2009; Cenarro & Trujillo 2009;
van Dokkum et al. 2009; Carrasco et al. 2010; van Dokkum et al. 2010; Cassata
et al. 2011; van der Wel et al. 2011; Szomoru et al. 2011, 2012; van der Wel
et al. 2014; Belli et al. 2014) driven mainly by accretion of gas-poor
satellites (van Dokkum et al., 2010, 2015), as expected in the two-phased
formation scenario. The $z\sim 2$ compact, quiescent early-type galaxies,
sometimes termed ‘red nuggets’, have typically large stellar masses of the
order of $M_{\star}\sim{10^{11}}\>\mathrm{M_{\odot}}$ while having small
effective radii, $R_{\mathrm{e}}\lesssim 1\>\mathrm{kpc}$.

The exact surface brightness profile shapes and rotation properties of red
nuggets remain debated. While the median Sérsic indices $n$ of $z=2$ quiescent
galaxies are comparable to the local massive ETGs with
$n_{\mathrm{median}}=3.7$ (Szomoru et al., 2012), there is also evidence for
surface brightness profiles that are closer to exponential. Selecting galaxies
at a constant number density of $n=2\times 10^{-4}\>\mathrm{Mpc}^{-3}$, van
Dokkum et al. (2010) found that massive galaxies become more exponential at
high redshifts, e.g. their mean Sérsic index decreases from $n\sim 6$ at $z=0$
down to $n\sim 2$ at $z=2$. Individual massive, quiescent, rotationally
supported disk galaxies above $z=2$ have also been discovered, as reported by
Toft et al. (2017) who found a disky (axis ratio
$0.59\begin{subarray}{c}+0.03\\\ -0.09\end{subarray}$) rapidly rotating system
at $z=2.1478$ with
$\log_{\mathrm{10}}(M_{\star}/\mathrm{M}_{\odot})=11.15\begin{subarray}{c}+0.23\\\
-0.23\end{subarray}$ and $R_{\mathrm{e}}=1.73\begin{subarray}{c}+0.34\\\
-0.27\end{subarray}\ \rm kpc$. The Sérsic index of the galaxy
$n=1.01\begin{subarray}{c}+0.12\\\ -0.06\end{subarray}$, is consistent with an
exponential system. The shape of the system is congruous with observations
that show that compact massive galaxies tend to have highly elliptic shapes
both at high redshift (van der Wel et al., 2011; Buitrago et al., 2013) and at
$z\sim 0.15$ (Trujillo et al., 2012).

Due to the stochastic nature of mergers in the hierarchical picture of
cosmological structure formation, it is possible for a subset of $z\sim 2$
massive compact galaxies to only experience a few mergers during their later
evolution, and thus remain relatively compact until the present day (Quilis &
Trujillo, 2013). Such massive, compact galaxies with uniformly old stellar
populations are termed local relic galaxies. Local relic candidates have been
identified with varying degrees of relicness quantified by their star
formation and accretion histories (e.g. Trujillo et al. 2009; Taylor et al.
2010; Ferré-Mateu et al. 2012; Trujillo et al. 2012; Poggianti et al. 2013;
van den Bosch et al. 2015; Saulder et al. 2015; Tortora et al. 2016). Local
($z\lesssim 0.1$) relic galaxies are relatively rare, with Quilis & Trujillo
(2013) estimating that fewer than $0.1\%$ of massive galaxies have accreted
less than $10\%$ of their stars after their formation, corresponding to a
number density of $n\sim 10^{-6}\>\mathrm{Mpc}^{-3}$ (Trujillo et al., 2014).

Despite their low expected numbers, local relics have been confirmed, with NGC
1277 (Trujillo et al., 2014), Mrk 1216 and PGC 032873 (Ferré-Mateu et al.,
2017) being the prime examples. The three local relic galaxies have stellar
masses ranging between ${1.2\times 10^{11}}\>\mathrm{M_{\odot}}\leq
M_{\star}\leq{2.9\times 10^{11}}\>\mathrm{M_{\odot}}$ and effective radii in
the range of $1.2\ \mathrm{kpc}\leq R_{\mathrm{e}}\leq 2.3\ \mathrm{kpc}$. The
stellar masses and effective radii with their uncertainties are displayed in
Table 1.

local relic | $M_{\star}$ | $R_{\mathrm{e}}$ | reference
---|---|---|---
galaxy | [${10^{11}}\>\mathrm{M_{\odot}}$] | [kpc] |
NGC 1277 | $1.2\pm 0.4$ | $1.2\pm 0.1$ | Trujillo et al. (2014)
Mrk 1216 | $2.0\pm 0.8$ | $2.3\pm 0.1$ | Ferré-Mateu et al. (2017)
PGC 032873 | $2.3\pm 0.9$ | $1.8\pm 0.2$ | Ferré-Mateu et al. (2017)
Table 1: The stellar masses $M_{\star}$ and effective radii $R_{\mathrm{e}}$
of three local relic galaxies NGC 1277, Mrk 1216 and PGC 032873. The galaxies
are described in more detail in the text.

Resembling the $z\sim 2$ red nuggets, NGC 1277, Mrk 1216, and PGC 032873 show
elongated, disky shapes, strong rotation and highly peaked velocity dispersion
profiles (Trujillo et al., 2014; Ferré-Mateu et al., 2017). A straightforward
comparison to the Sérsic profiles of red nuggets is challenging as multiple
Sérsic components are required for good fits of the surface brightness
profiles for both NGC 1277 and Mrk 1216 (Yıldırım et al., 2015). In addition
to local relics, various survey programs have also revealed the existence of
several tens of relic candidates at intermediate redshifts ($0.1\lesssim
z\lesssim 1.0$, e.g. Valentinuzzi et al., 2010; Damjanov et al., 2014;
Spiniello et al., 2021, 2024).

The most massive present-day ($z\sim 0$) early-type galaxies (ETGs) differ
considerably from $z\sim 2$ compact red quiescent galaxies and the local
relics in their structural and kinematic properties. The brightest
($M_{\mathrm{B}}\lesssim-20.5$) and most massive
($M_{\star}\gtrsim{10^{11}}\>\mathrm{M_{\odot}}$) early-type galaxies observed
in the local Universe are characterized by their extended sizes, shallow
central surface brightness profiles and slow rotation (e.g. Kormendy & Bender
1996; Faber et al. 1997; Quenneville et al. 2024). Massive galaxies also
ubiquitously host supermassive black holes in their nuclei (e.g. Kormendy &
Richstone 1995; Ferrarese & Ford 2005; Kormendy & Ho 2013). Many host galaxy
properties have been found to tightly correlate with the mass of their SMBHs,
including the bulge mass ($M_{\bullet}$-$M_{\star}$ relation; Häring & Rix
2004) and the velocity dispersion ($M_{\bullet}$-$\sigma$ relation; Ferrarese
& Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002) suggesting a co-
evolution of massive galaxies and their SMBHs.

Binaries of SMBHs formed in the aftermath of galaxy mergers have been
suggested to explain the observed faint, diffuse central cores of massive
early-type galaxies, and there are indications of core-SMBH relations linking
the SMBH mass to the amount of missing starlight in the core region (mass
deficit; Graham 2004; de Ruiter et al. 2005; Lauer et al. 2007; Kormendy &
Bender 2009) as well as the spatial extent of the cores (Rusli et al., 2013;
Thomas et al., 2016). The characteristic kinematic fingerprint of the SMBH
binary core scouring, the tangentially biased velocity anisotropy ($\beta<0$)
in the core region, has also been observed in massive, cored early-type
galaxies (Thomas et al., 2014).

Theoretical arguments and later detailed numerical simulations have
established a picture of the formation, evolution, and mergers of SMBH
binaries in the aftermath of mergers of massive gas-poor galaxies (e.g.
Begelman et al. 1980; Ebisuzaki et al. 1991; Quinlan 1996; Quinlan & Hernquist
1997; Milosavljević & Merritt 2001, 2003; Merritt & Milosavljević 2005;
Merritt 2006, 2013; Nasim et al. 2021). First, dynamical friction
(Chandrasekhar, 1943) from stars, gas and dark matter causes the two SMBHs to
sink to the center of the merger remnant, and the SMBHs form a gravitationally
bound binary. The binary subsequently interacts with the stars in the galactic
nucleus, losing energy to the surrounding stars in complex three-body
interactions while shrinking towards smaller separations. Only stars with a
sufficiently small amount of angular momentum (the so-called loss-cone
population) can interact with the SMBH binary, and thus the binary
preferentially ejects stars from radial orbits, leading to a slow build-up of
stars on circular orbits in the core region (e.g. Thomas et al. 2014). The
rapidly lowered central stellar density and the more slowly developed
tangentially biased velocity dispersion of the core can be viewed as the two
phases of the core formation, transforming the orbit structure of the stars in
the nucleus (Frigo et al., 2021)

The so-called final-parsec problem, i.e. the stalling of the SMBH binary at
parsec-scale separations due to insufficient loss-cone refilling, seems to be
avoided in realistic simulation models of axisymmetric, triaxial or rotating
nuclei (e.g. Berczik et al. 2006; Preto et al. 2011; Khan et al. 2011; Khan et
al. 2013; Holley-Bockelmann & Khan 2015; Khan et al. 2016; Gualandris et al.
2017; Rantala et al. 2017; Mirza et al. 2017), although questions related to
the potentially insufficient mass resolution of the current simulation models
still remain (Vasiliev et al., 2014). Finally, at centiparsec separations,
relativistic radiation-reaction effects become important (Poisson & Will,
2014), and the binary radiates its remaining energy and angular momentum
(Peters & Mathews, 1963; Peters, 1964) until merging into a single SMBH in a
burst of gravitational radiation.

Despite the general success of the theoretical and numerical models explaining
the evolution of compact quiescent galaxies at $z\sim 2$ into the present-day
population of massive ETGs, open questions still remain (see Naab & Ostriker
2017 for a review). Large cosmological box runs and zoom-in simulations that
probe the size growth and rotation of ETGs from high redshift to the present-
day typically ignore the small-scale SMBH binary dynamics (and thus core
formation) due to the instant merger procedure at large separations of the
order of $\sim 0.1$ kpc (see e.g. the discussion of the merger procedure in
Johansson et al. 2009).

However, recently Mannerkoski et al. (2021, 2022) performed the first
cosmological zoom-in simulations that include accurate, non-softened dynamics
with post-Newtonian equations of motion for SMBH interactions using the
GADGET-3-based regularized tree code KETJU (Rantala et al., 2017). While the
code can in principle follow the SMBHs and SMBH binary dynamics including core
scouring in cosmologically evolving early-type galaxies, the accurate SMBH
dynamics requires that the mass ratio between the SMBHs and stellar population
particles is large enough ($M_{\bullet}/m_{\star}\gtrsim 100$). This
corresponds to accurate small-scale SMBH dynamics at redshifts $z\lesssim 1$
in these particular simulations. However, improved versions of the code and
upcoming high performance computing hardware capabilities will ameliorate the
problem in the future allowing for accurate SMBH binary dynamics for a wider
range in redshift (Keitaanranta et al., in prep.).

The main question to be answered in this study is whether the $z=2$ compact
early-type galaxies can be the progenitors of present-day massive, cored and
slowly rotating early-type galaxies. We especially focus on the question
whether the sizes, cores and rotation properties can be simultaneously
explained by a single physical process, namely major and minor galaxy mergers
including SMBHs. As simulations in a full cosmological setting with accurate
SMBH dynamics below $z\lesssim 2$ is still beyond the current computational
capabilities, we choose here instead to run samples of controlled, idealized
gas-free galaxy mergers.

The article is structured as follows. After the introduction we review the
GADGET-4-based simulation code KETJU and the initial conditions for our
numerical simulations in Section 2. We present the results on galaxy size
growth and SMBH binary core scouring in Section 3 and Section 4. Next, in
Section 5 we describe how core scouring can affect the global properties of
ETGs, especially for large mass deficits. Finally, we discuss the shape,
velocity anisotropy and kinematic structure of our simulated galaxies in
Section 6 and conclude in Section 7.

##  2 Methods and numerical simulations

###  2.1 Simulation code

####  2.1.1 The KETJU code: summary

We use the recent publicly
available111https://www.mv.helsinki.fi/home/phjohans/ketju implementation
(Mannerkoski et al., 2023) of the KETJU simulation code (Rantala et al.,
2017). The novel code version is itself based on the tree / fast multipole
method (FMM) gravity solver of GADGET-4 (Springel et al., 2021). The main
functionality of KETJU is to add an algorithmically regularized region around
each SMBH in the simulation in which the gravitational dynamics can be
accurately calculated without the need to employ gravitational softening. The
size of the regularized regions typically ranges from a few parsecs up to tens
of parsecs chosen based on the stellar density around each SMBH so that the
maximum number of stellar particles within any regularized region does not
exceed $N_{\star}\sim 5000$–$10000$. The spatial size of the regularized
region places an upper limit for the softening length
$\epsilon=\epsilon_{\star}=\epsilon_{\bullet}$ of SMBHs and stellar particles
in GADGET-4. This is because KETJU requires
$r_{\mathrm{region}}>2.8\times\epsilon$ in order to have all SMBH-SMBH and
star-SMBH interactions non-softened. In this work only SMBH and star particles
participate in the regularization integration while all interactions including
dark matter (DM) particles are always softened. Star-star interactions within
the regularized regions are softened as well.

####  2.1.2 MSTAR: the KETJU integrator

The algorithmically regularized (e.g. Mikkola & Tanikawa 1999; Preto &
Tremaine 1999) integrator of KETJU is based on the MSTAR integrator of Rantala
et al. (2020). MSTAR utilizes three main ingredients to achieve accurate and
efficient orbital integration. First, the equations of motion of the particles
in the regularized regions are time-transformed. Together with the use of the
standard leapfrog integration, the time transformation avoids the issues of
the diverging Newtonian potential at small separations (Mikkola & Tanikawa,
1999; Mikkola & Aarseth, 2002; Mikkola & Merritt, 2006, 2008). Next, MSTAR
uses a specialized minimum spanning tree (MST) coordinate system to reduce the
numerical floating-point round-off error (Rantala et al., 2020). Finally, in
order to reach the user-desired integration accuracy, MSTAR employs the Gragg-
Bulirsch-Stoer (GBS) extrapolation method (Gragg, 1965; Bulirsch & Stoer,
1966). The computational cost of the GBS method is mitigated by efficient two-
fold MPI parallelization of the code, based on the fact that the different
sub-step divisions (e.g. Deuflhard 1983) required by the method are
independent of each other (Rantala et al., 2020).

Compared to the version of Rantala et al. (2020), the MSTAR version (the KETJU
integrator) included in the public KETJU code (Mannerkoski et al., 2023)
includes a dynamical order control for the GBS extrapolation procedure (Hairer
et al., 1993) and an optional star-star softening. In addition, the code
includes post-Newtonian (PN) equations of motion up to order PN3.5 with
optional spin terms for the SMBH-SMBH interactions (Thorne & Hartle, 1985;
Blanchet, 2014), and finally a prescription for the relativistic gravitational
wave (GW) recoil kicks at SMBH-SMBH mergers (Zlochower & Lousto, 2015). For
additional PN and GW implementation details, see Mannerkoski et al. (2023).

####  2.1.3 GADGET-4: galactic-scale dynamics

Outside the regularized regions, the gravitational dynamics in KETJU is
treated as in the standard GADGET-4 code. Compared to the previous KETJU
implementations (Rantala et al., 2017) based on GADGET-3 (e.g. Springel 2005),
the GADGET-4 (Springel et al., 2021) based KETJU (Mannerkoski et al., 2023)
has several advantages. First, GADGET-4 has an option to use the so-called
hierarchical integration scheme based on hierarchical Hamiltonian splitting
instead of the common block time-step scheme. The hierarchical scheme is more
efficient than the block step scheme for systems with extended time-step
hierarchies as the rapidly evolving parts of the simulation decouple from the
more slowly evolving regions. Importantly, hierarchical integration allows for
manifest momentum conservation if coupled with a pair-wise force calculation
method such as direct summation or the FMM, which is not possible using the
traditional block time-step scheme.

GADGET-4 includes an accurate fast multipole method (FMM) gravity solver
(Greengard & Rokhlin 1987; Dehnen 2000, 2002, 2014; Zhu 2020) in addition to
the traditional one-sided tree code (Appel, 1985; Barnes & Hut, 1986; McMillan
& Aarseth, 1993; Springel et al., 2001; Springel, 2005). The FMM force
calculations between particles (and tree nodes) are always pair-wise and the
corresponding particle kicks are always synchronous, together with the
hierarchical integration making the force calculation and time integration
algorithm manifestly momentum conserving. This could in principle be very
important for SMBH binary hardening phases in KETJU, as spurious loss-cone
filling originating from one-sided tree force errors might potentially affect
the SMBH binary hardening rates, and thus their merger time-scales (Gualandris
et al., 2017). However, as shown by Mannerkoski et al. (2023) in their Fig. 8
the loss cone filling rates are nearly identical when using the one-sided tree
compared to the FMM solver.

####  2.1.4 Simulation code parameters

We use the following options, user-given accuracy and error tolerance
parameters for the GADGET-4-KETJU code. For GADGET-4, we use hierarchical
integration and the FMM gravity solver with second-order multipoles. The
GADGET-4 force and integration error tolerances are set to $\alpha=0.005$ and
$\eta=0.002$. For the MSTAR parameters, the Gragg-Bulirsch-Stoer per-step
relative error tolerance parameter is set to $\epsilon_{\mathrm{GBS}}=10^{-7}$
and the output time relative tolerance parameter is
$\epsilon_{\mathrm{t}}=10^{-3}$. Post-Newtonian terms are included in the
SMBH-SMBH interactions up to PN3.5.

The gravitational interactions of the dark matter particles are softened using
a Plummer-equivalent softening length of $\epsilon_{\mathrm{dm}}=65$ pc. The
stellar component has a softening length of $\epsilon_{\star}=6.5$ pc. The
size of the regularized regions around the SMBHs is set to
$r_{\mathrm{region}}=20$ pc, fulfilling the necessary criterion of
$r_{\mathrm{region}}>2.8\times\epsilon_{\star}$ (Rantala et al., 2017),
resulting in typically a few hundred stellar particles in the regions. Star-
star softening is used in the regularized regions (Mannerkoski et al., 2021;
Mannerkoski et al., 2023). SMBH-SMBH and SMBH-star interactions are always
unsoftened. In this study the interactions between the SMBHs and the dark
matter particles are always softened, and dark matter particles are not
included in the regularized integration. We note that a scenario in which
SMBHs and dark matter can strongly scatter each other has been recently
explored by Partmann et al. (2023) using GADGET-3-KETJU.

###  2.2 Initial conditions for galaxy models

####  2.2.1 Isolated galaxy models: spherical and disky stellar components

We generate isolated multi-component early-type galaxy initial conditions
(ICs) later to be merged with each other in our merger simulation setups. We
use a modified version of the standard galaxy IC generation code
MakeDiskGalaxy (Springel et al., 2005), which constructs equilibrium galaxy
models consisting of a dark matter halo, a stellar component (an axisymmetric
thick disk or a spherical bulge) and a central SMBH. We construct three types
of isolated models: two host galaxy progenitor models (spherical, disky) and
lower-mass spherical satellite galaxy models.

The host galaxy stellar mass is set to
$M_{\mathrm{\star,host}}={10^{11}}\>\mathrm{M_{\odot}}$ while the satellites
have $M_{\mathrm{\star,sat}}={2\times 10^{10}}\>\mathrm{M_{\odot}}$, resulting
in a satellite-host mass ratio of
$M_{\mathrm{\star,sat}}/M_{\mathrm{\star,host}}=1/5$. The host galaxies have
initially a three-dimensional stellar half-mass radius of $r_{\mathrm{h}}=1\
\mathrm{kpc}$. We construct satellites with five different half-mass radii
ranging from $r_{\mathrm{h}}=0.97$ kpc to $r_{\mathrm{h}}=3.36$ kpc, following
the observed mass-size relation for early-type galaxies (van der Wel et al.,
2014) between redshifts $z=2$ and $z=0$. In minor mergers beyond generation
five, we use the same fifth ($z=0$) satellite model. The basic physical
properties of the isolated IC galaxy models are listed in Table 2. The stellar
components of the galaxy models are realized using $N_{\mathrm{\star}}=10^{5}$
per host galaxy and $N_{\mathrm{\star}}=2\times 10^{4}$ for the satellites.
The necessary condition (Rantala et al., 2017) for well resolved stellar
component for the SMBH-star interactions in KETJU
($M_{\bullet}/m_{\star}\gtrsim 100$ with $m_{\star}=M_{\star}/N_{\star}$) is
always fulfilled as $\min{(M_{\bullet})}/m_{\star}=400$ in the simulations of
this study.

The host galaxy stellar components are either spherical (S) or disky222We note
that ‘diskiness’ in this study refers to the thick, axisymmetric shape of the
stellar component and not to the traditional disky versus boxy isophotal shape
classification of ETGs. (D). The satellites galaxy models are always
spherical. The spherical stellar components follow the Hernquist profile
(Hernquist, 1990) with a central power-law density slope $(\rho(r)\propto
r^{-\gamma})$ of $\gamma=1$. The Hernquist stellar models are isotropic and
non-rotating. The disky models are rotationally supported and have exponential
surface density profiles. The axisymmetric disky models are thick in the
vertical direction with the minor-to-major axis ratio being $c/a=0.5$ within
$r_{\mathrm{h}}$.

####  2.2.2 Dark matter halos

The dark matter halos of the galaxies follow the Hernquist profile, just as
the spherical stellar components. The host halo mass of
$M_{\mathrm{dm}}={10^{12}}\>\mathrm{M_{\odot}}$ is chosen following the
simulation setup of Hilz et al. (2012); Hilz et al. (2013). In general, it is
thought that dark matter may play a minor role in the central parts of high-
redshift massive galaxies due to their early dissipative evolution (e.g. Toft
et al. 2012). In addition NGC 1277, a local relic galaxy is known to be dark
matter deficient (Comerón et al., 2023). Motivated by such arguments, the
Hernquist scale radii of the host halos are chosen so that the dark matter
fraction $f_{\mathrm{dm}}=M_{\mathrm{dm}}/(M_{\mathrm{dm}}+M_{\star})$ is low,
only $f_{\mathrm{dm}}=3\%$ within $1\times r_{\mathrm{h}}$ and
$f_{\mathrm{dm}}=21\%$ within $5\times r_{\mathrm{h}}$. The satellite galaxy
models always have $f_{\mathrm{dm}}=20\%$ within the half-mass radius, more
typical of present-day massive early-type galaxies (Cappellari et al., 2013;
Courteau & Dutton, 2015). The host galaxy halos are realized with
$N_{\mathrm{dm}}=10^{6}$ particles while the satellite halos consist of
$N_{\mathrm{dm}}=2\times 10^{5}$ particles.

####  2.2.3 SMBH properties

We construct two versions of each galaxy model, with and without a central
SMBH. The SMBHs are placed at rest at the center of the isolated galaxy model.
The host galaxy SMBH mass is set to $M_{\mathrm{\bullet,host}}={2\times
10^{9}}\>\mathrm{M_{\odot}}$ following the SMBH mass estimates of Graham et
al. (2016) for the local relic galaxy NGC 1277. The satellite SMBH mass is set
to $M_{\mathrm{\bullet,sat}}={4\times 10^{8}}\>\mathrm{M_{\odot}}$
corresponding to the satellite-host black hole mass ratio of $1/5$.

model | shape |  $\rho(r)$, $\Sigma(r)$ | rotation | $M_{\mathrm{dm}}$ | $M_{\mathrm{\star}}$ | $M_{\mathrm{\bullet}}$ | $r_{\mathrm{h,\star}}$ | reference $z$
---|---|---|---|---|---|---|---|---
|  |  |  | [$10^{11}\mathrm{M_{\odot}}$] | [$10^{11}\mathrm{M_{\odot}}$] | [$10^{9}\mathrm{M_{\odot}}$] | [kpc] |
host-S | spherical | Hernquist | $\times$ | $10.0$ | $1.0$ | 2.0 | $1.00$ | $2.00$
host-D | axisymmetric | exponential | $\checkmark$ | $10.0$ | $1.0$ | 2.0 | $1.00$ | $2.00$
sat-1 | spherical | Hernquist | $\times$ | $2.0$ | $0.2$ | 0.4 | $0.97$ | $2.00$
sat-2 | spherical | Hernquist | $\times$ | $2.0$ | $0.2$ | 0.4 | $1.25$ | $1.56$
sat-3 | spherical | Hernquist | $\times$ | $2.0$ | $0.2$ | 0.4 | $1.68$ | $1.13$
sat-4 | spherical | Hernquist | $\times$ | $2.0$ | $0.2$ | 0.4 | $2.48$ | $0.69$
sat-5 | spherical | Hernquist | $\times$ | $2.0$ | $0.2$ | 0.4 | $3.36$ | $0.25$
Table 2: The initial galaxy models of this study. The redshift $z$ is for the
reference of the mass-size relation as discussed in the main text. The disky
shape refers to an axisymmetric profile with axis ratios $b/a=1.0$ and
$c/a=0.5$ with an exponential stellar surface density profile. Two versions of
each profile are generated, one with a central SMBH and another without.

###  2.3 Simulation sequences

####  2.3.1 Merger orbits

We adopt the basic galaxy merger orbit configuration from the previous KETJU
ETG merger studies (Rantala et al., 2017, 2018, 2019). The initial separation
of the nuclei of the galaxies is set to $r_{\mathrm{sep}}=30\ \mathrm{kpc}$.
The two merging galaxies are placed on a nearly parabolic orbit (Khochfar &
Burkert, 2006; Naab et al., 2006) with the initial relative velocity scaled by
a factor of $f_{\mathrm{vel}}=M(<r_{\mathrm{sep}})/M_{\mathrm{tot}}$, ensuring
a relatively rapid galaxy merger and limiting the computational expense of the
simulations. For all the minor merger simulations in this study the initial
pericenter distance $r_{\mathrm{p}}$ is fixed to $r_{\mathrm{p}}=2$ kpc.

For the major merger orbits of the spherical models we use two different
initial pericenter distances ($r_{\mathrm{p,lo}}=1$ kpc and
$r_{\mathrm{p,hi}}=9$ kpc) resulting in a slowly rotating and a fast-rotating
major merger remnant in the case of the spherical progenitors. Somewhat
counter-intuitively for the spherical non-rotating progenitors, the major
merger orbit with $r_{\mathrm{p,lo}}=1$ kpc results in the fast-rotating
merger remnant and $r_{\mathrm{p,hi}}=9$ kpc in a slowly rotating remnant,
even though the former orbit has more angular momentum than the latter. The
reason for this lies in how the stellar orbital angular momentum is re-
distributed during the merger and how much of it is still remaining when the
stellar bulges merge. The process is explained in more detail in Appendix A.
For the major merger simulations between disky models we use
$r_{\mathrm{p,mi}}=3$ kpc.

####  2.3.2 Merger sequences: overview

Figure 1: The four basic types of merger sequences performed for this study.
The sat-5 satellite model is used in minor merger generations beyond
generation five. Top left: the D-MA sequence. Two disky ETGs are first merged,
after which minor mergers are repeated until the remnant reaches the $z=0$
mass-size relation for ETGs. Top right: the D-MI sequence. A disky progenitor
experiences a sequence of minor mergers without the initial major merger.
Bottom panels: the merger sequences S-MA and S-MI. The two merger sequences
correspond to the sequences on the top row, but now the progenitor galaxies
are spherical and non-rotating. This study includes in total ten merger
sequences. Each of the sequences D-MA, D-MI, S-MA and S-MI is run twice, with
and without SMBHs. In addition, the setup S-MA is repeated with an alternative
merger orbit, both with and without SMBHs, resulting in a fast-rotating
progenitor after the major merger and before the minor merger sequence. These
two sequences are labeled S-MA-SROT and S-MA-FROT depending on the rotation
after the major merger. Note that the sizes of the galaxies are just for
illustrative purposes and do not reflect the relative masses or spatial
extents of the stellar components.

We construct in total $10$ different merger sequences from the host galaxies
and satellites described in the previous section and in Table 2. We mainly
focus on varying the following simulation parameters. First, we study merger
sequences with different progenitor galaxy shapes and rotation properties,
i.e. whether the host galaxy is spherical (S) or disky (D). Next, we explore
the effect of SMBHs by running each merger sequence with (label BH) and
without central SMBHs (no label). Finally, we investigate how including a
single major merger (MA) or alternatively excluding it (MI) in the merger
sequence affects the evolution of the merger remnant. If a major merger is
included in the merger sequence, the sequence begins with it, followed by the
multiple minor mergers.

These three options (disky/spherical host, major merger or not, SMBHs or not)
yield $8$ different basic merger sequences, outlined in Fig. 1. The last $2$
merger sequences explore different galaxy major merger orbit configurations
leading to very different merger remnant rotation properties (slow or fast
rotation) from the same progenitors.

The satellite galaxies are merged into the host progenitors from isotropic
random directions but in such a manner that any Nth minor merger occurs from
the same direction in all the merger sequences. The minor mergers are
continued until the merger remnant galaxies reach the $z=0$ mass-size
($M_{\star}$-$R_{\mathrm{e}}$) relation for early-type galaxies. The required
size growth to reach the observed relation is substantial as the effective
radius of the host galaxy progenitors is only $R_{\mathrm{e}}\sim 0.75\
\mathrm{kpc}$, and a size growth by a factor of $\sim 5$–$10$ is needed to
reach the relation (Shen et al., 2003). We run $0$ or $1$ major mergers and
$7$–$8$ minor mergers for each ten merger sequences. In merger generation
$6$–$8$ we use the satellite model sat-5.

####  2.3.3 Merger sequences with disky progenitors

sequence | host | major merger | SMBH
---|---|---|---
D-MA-BH | disky | $\checkmark$ | $\checkmark$
D-MA | disky | $\checkmark$ | $\times$
D-MI-BH | disky | $\times$ | $\checkmark$
D-MI | disky | $\times$ | $\times$
Table 3: The four merger sequences from disky progenitors (D-MA-BH, D-MA,
D-MI-BH, D-MI) with ($\checkmark$) and without ($\times$) a major merger and
with and without central SMBHs.

The four merger sequences with a disky progenitor galaxy (D-MA-BH, D-MA, D-MI-
BH, D-MI) are listed in Table 3. The sequences D-MA-BH and D-MA begin with a
major galaxy merger while the sequences D-MI-BH and D-MI do not include a
major merger. Central SMBHs are included in sequences D-MA-BH and D-MI-BH. In
the disky major merger simulations one galaxy lies initially in the orbital
plane while the other is tilted by $\theta=45^{\circ}$.

####  2.3.4 Merger sequences with spherical progenitors

sequence | host | major | rotation after | SMBH
---|---|---|---|---
|  | merger | major merger |
S-MA-SROT-BH | spherical | $\checkmark$ | slow | $\checkmark$
S-MA-SROT | spherical | $\checkmark$ | slow | $\times$
S-MA-FROT-BH | spherical | $\checkmark$ | fast | $\checkmark$
S-MA-FROT | spherical | $\checkmark$ | fast | $\times$
S-MI-BH | spherical | $\times$ |  | $\checkmark$
S-MI | spherical | $\times$ |  | $\times$
Table 4: The six merger sequences beginning from spherical progenitors.

The mergers beginning from spherical progenitors are listed in Table 4. There
are $6$ merger sequences in total. In addition to the sequences with and
without a major merger and SMBHs we perform the major merger simulations with
two different merger orbit configurations by varying the initial pericenter
separation $r_{\mathrm{p}}$. The two major merger orbit configurations result
in a slow rotating merger remnant and a fast rotating merger remnant. The
minor merger sequences from these two progenitors are labeled as S-MA-SROT-BH,
S-MA-SROT (slow rotation, with and without SMBHs) and S-MA-FROT-BH, S-MA-FROT
(fast rotation, with and without SMBHs). The two sequences without a major
merger are S-MI-BH and S-MI depending if SMBHs is present in the minor mergers
or not.

####  2.3.5 SMBH binary eccentricities

Figure 2: The eccentricities of the SMBH binaries in the simulation sequences
at a binary semi-major separation of $a=2$ pc. At this separation the
gravitational wave emission has not yet circularized the binary orbit. Most of
the binaries are very eccentric with $e\gtrsim 0.9$. A number of less
eccentric binaries occur in the samples S-MA-FROT-BH and D-MI-BH.

The SMBH binary eccentricities in the minor merger simulations at binary semi-
major axis of $a=2$ pc are shown in Fig. 2. The particular semi-major axis was
chosen so that even the most eccentric SMBH binaries in our simulations have
not yet entered the dynamical regime dominated by gravitational wave emission.

Most of the SMBH binaries in the minor merger simulations are highly eccentric
with $e>0.9$, especially in the runs S-MA-SROT-BH, S-MI-BH and D-MA-BH. The
simulation sequence S-MA-FROT-BH includes four moderately eccentric binaries
with $0.5<e<0.8$ while the sequence D-MI-BH involves two binaries with low and
moderate eccentricities with $e=0.08$ and $e=0.67$. The in general high binary
eccentricities likely somewhat depend on the selected galaxy merger orbits,
although there is recent evidence that stochastic processes may play a
significant role in determining the initial SMBH binary eccentricities
(Rawlings et al., 2023).

##  3 The size growth of massive early-type galaxies via major and minor
mergers

###  3.1 Reaching the $z=0$ mass-size relation

Figure 3: The evolution of the mass-size relation in the simulation sequences including SMBH: S-MA-SROT-BH, S-MA-FROT-BH and S-MI-BH (left panel) and D-MA-BH with D-MI-BH in the right panel. Both panels also show observational comparison data from local relics NGC 1277, PGC 032873 and Mrk 1216 (Trujillo et al., 2014; Ferré-Mateu et al., 2017), intermediate-redshift ($0.1<z<0.5$) relic candidates (Spiniello et al., 2021, 2024) and high redshift ($1.5<z<2.5$) quiescent (sSFR<$0.3/t_{\mathrm{Hubble}}$) galaxies (Szomoru et al., 2013). The local SDSS $z=0$ mass-size relation for ETGs (Shen et al., 2003) is displayed for comparison purposes. In general, $5$–$8$ minor mergers are required to reach the $z=0$ relation starting from a compact progenitor galaxy. Major mergers result in size growth $R_{\mathrm{e}}\propto M_{\star}^{\alpha_{\bullet}}$ with $\alpha_{\bullet}\sim 1$, almost parallel to the relation, so more minor mergers are required to reach the local mass-size relation in the merger samples beginning with a major merger. For minor mergers with SMBHs the size growth is steeper with $\alpha_{\bullet}>2.3$. This is further discussed in the text. Figure 4: The evolution of the mass-size relation as in Fig. 3 but for simulation samples without SMBHs: S-MA-SROT, S-MA-FROT and S-MI (left) and D-MA with D-MI (right). The size growth in the merger sequences is less prominent than in the corresponding simulations with SMBHs. sequence | major+SMBH | major | minor+SMBH | minor
---|---|---|---|---
S-MA-SROT | 1.02 | 0.92 | 2.31 | 1.80
S-MA-FROT | 1.08 | 0.92 | 3.10 | 2.10
S-MI | - | - | 2.51 | 2.33
D-MA | 0.54 | 0.48 | 2.81 | 1.79
D-MI | - | - | 2.46 | 2.10
Table 5: The power-law exponents $\alpha$ of the mass-size relation
$R_{\mathrm{e}}\propto M_{\star}^{\alpha}$ for the major and minor merger
sequences.

First we focus on the evolution of the mass-size
($M_{\star}$-$R_{\mathrm{e}}$) relation of the initially compact early-type
galaxies in the merger sequences. The evolution of the mass-size relation in
the consecutive merger simulations is presented in Fig. 3 (including SMBHs)
and in Fig. 4 (without SMBHs). Both our spherical and disky progenitors at
$M_{\star}={10^{11}}\>\mathrm{M_{\odot}}$ have initial effective radii of
$R_{\mathrm{e}}=0.75$ kpc. Compared to the local SDSS ETG mass-size relation
(Shen et al., 2003) at the same mass, the progenitors are by a factor of $\sim
5.6$ smaller in spatial size.

In our simulations the initially compact galaxy models can reach the $z=0$
mass-size relation after $6$–$8$ generations of minor mergers. As already
expected from simple virial arguments (Naab et al., 2009; Bezanson et al.,
2009), minor mergers are more efficient than major mergers in increasing the
galaxy sizes. A major merger moves a galaxy in the mass-size plane in an
almost parallel direction to the $z=0$ relation, thus requiring more minor
mergers in the sequences beginning with a major merger to reach the relation.

The power-law size growth of galaxies in merger sequences is typically
characterized as $R_{\mathrm{e}}\propto M_{\star}^{\alpha}$ with $\alpha>2$
required to explain the observed size growth of ETGs since $z=2$ (van Dokkum
et al., 2010). The values of the galaxy size growth power-law index $\alpha$
from our mergers sequences are collected in Table 5. In our simulations,
equal-mass major mergers of spherical progenitors show an almost linear trend
in size growth with $0.92\lesssim\alpha\lesssim 1.08$, consistent with the
results in the literature (Boylan-Kolchin et al., 2005; Ciotti et al., 2007;
Bezanson et al., 2009; Nipoti et al., 2009; Hilz et al., 2012; Hilz et al.,
2013). Mergers including SMBHs result in somewhat larger values for $\alpha$
with $\alpha=0.92$ (without SMBHs) and $1.02\lesssim\alpha_{\bullet}\lesssim
1.08$ (with SMBHs). For the disky progenitors the size growth in a major
merger is more limited with $0.48\lesssim\alpha\lesssim 0.54$.

Studying our minor merger simulations, the power-law exponent $\alpha$ of the
galaxy mass-size relation is close to $\alpha\sim 1.8$ in merger setups S-MA-
SROT and D-MA, and $\alpha\gtrsim 2.1$ in all other minor merger setups,
consistent with the observed size growth of ETGs since $z=2$ (van Dokkum et
al., 2010). The size growth with $\alpha=2.33$ in the spherical minor merger
sample without SMBHs S-MI is very close to the results of Hilz et al. (2013)
who used a similar setup to ours, although with less compact progenitor
galaxies.

Including SMBHs in the merger sequences increases $\alpha$ in all minor merger
simulations. The effect is especially prominent in the samples S-MA-SROT
($\alpha_{\bullet}-\alpha=0.51$), S-MA-FROT ($\alpha_{\bullet}-\alpha=1.0$)
and D-MA ($\alpha_{\bullet}-\alpha=1.02$), each of which begins with a major
galaxy merger. For the simulations without a major merger
$\alpha_{\bullet}-\alpha=0.18$ for the merger sequence S-MI and
$\alpha_{\bullet}-\alpha=0.36$ for the simulations D-MI. We explore this
phenomenon in Section 5.1, where we argue that the reason for the additional
size growth is simply the displacement of stellar material due to core
scouring by high-mass SMBHs.

###  3.2 Comparison to observed compact ETGs

Next, we compare our simulations results to the observed masses and sizes of
three local relic galaxies (NGC 1277, PGC 032873 and Mrk 1216; Trujillo et al.
2014; Ferré-Mateu et al. 2017), relic candidates at $0.1<z<0.5$ (Spiniello et
al., 2021, 2024) and massive quiescent333We define a galaxy as quiescent if it
has a specific star formation rate sSFR$<0.3\>t_{\mathrm{Hubble}}$. galaxies
at $1.5<z<2.5$ (Szomoru et al., 2013). The observed masses and radii of the
galaxies are displayed alongside our simulation results in the mass-size plane
in Fig. 3 and Fig. 4.

Our initial isolated galaxy models have masses
($M_{\star}={10^{11}}\>\mathrm{M_{\odot}}$) comparable both to the high-
redshift quiescent galaxies and the intermediate-redshift relic candidates
while having somewhat smaller effective radii than most of the observed
compact galaxies. In the simulation samples without major mergers, the
location of the simulated galaxies in the mass-size plane remains consistent
with the observed high and intermediate redshift galaxies during the first
$N_{\mathrm{merg}}=2$-$3$ minor mergers. After $N_{\mathrm{merg}}\geq 3$ the
galaxies experiencing minor mergers have larger effective radii
$R_{\mathrm{e}}$ than the observed galaxies in the samples. Three of our
models with two minor mergers (S-MI-BH, D-MI-BH and S-MI) are consistent with
the stellar mass and effective radius of the local relic NGC 1277. Compared to
the two other local relics, PGC 032873 and Mrk 1216, our minor merger models
without a major merger have in general larger sizes for their stellar mass.

The merger samples beginning with a major merger already have larger stellar
masses than any of the galaxies in the high-redshift quiescent galaxy sample
after the major merger. Nevertheless, they have effective radii of $1.05\
\mathrm{kpc}\lesssim R_{\mathrm{e}}\lesssim 1.54\ \mathrm{kpc}$ comparable to
the most massive galaxies in the high-redshift quiescent sample, $1.2\
\mathrm{kpc}\lesssim R_{\mathrm{e}}\lesssim 1.6\ \mathrm{kpc}$. Comparing to
the intermediate redshift relic candidates, the simulated galaxies in the
samples including a major merger have two observational counterparts with
masses close to $\log_{\mathrm{10}}(M_{\star}/\mathrm{M}_{\odot})=11.4$. The
simulated galaxies have effective radii comparable to these galaxies $(1.55\
\mathrm{kpc}\lesssim R_{\mathrm{e}}\lesssim 1.71\ \mathrm{kpc})$ after a major
merger and $N_{\mathrm{merg}}=1$-$3$ minor mergers, depending on the
simulation sample. The local relic PGC 032873 has also a similar stellar mass
$\log_{\mathrm{10}}(M_{\star}/\mathrm{M}_{\odot})=11.36$ and effective radius
($R_{\mathrm{e}}=1.8$ kpc) as our simulated galaxies after a major merger and
$N_{\mathrm{merg}}=2$ minor mergers. The third local relic, Mrk 1216 lies
between the major and minor merger simulation samples in the mass-size plane.
These results highlight the possibility that local relic galaxies may have
experienced a few mergers, but not the full sequence to evolve into the
largest ETGs in the local Universe.

##  4 Evolution of galactic nuclei in major and minor mergers

###  4.1 Surface density profiles

Figure 5: The azimuthally averaged radial stellar surface density profiles
$\Sigma_{\star}(R)$ of the merger remnants in the simulation samples including
SMBHs (left) column and without SMBHs (right column). Each profile is averaged
over $N_{\mathrm{los}}=50$ random viewing orientations. The number of mergers
increases from $N_{\mathrm{merg}}=1$ to $N_{\mathrm{merg}}=7$–$9$ as indicated
in the legend. For comparison, we show the surface density profiles of
selected local core ETGs (Rusli et al., 2013), typical Virgo ETGs (Kormendy &
Bender, 2009) and $1.5<z<2.5$ quiescent galaxies (Szomoru et al., 2013). The
central parts of the of the high-redshift quiescent galaxies are not resolved
at radii below $r\lesssim 0.7$ kpc (Szomoru et al., 2012), indicated by a
triangle symbol at the top of each panel. Only models with SMBHs (excluding
S-MI-BH which is further discussed in the text) resemble the surface density
profiles of local, massive, cored ETGs after the merger sequences while
without SMBHs the final surface density profiles are mildly cuspy at the
centers with values for central $\Sigma_{\star}(R)$, both from spherical and
disky initial conditions.

We show the stellar surface density profiles of our simulated galaxy sample
after each merger in Fig. 5. Each surface density profile is averaged over
$N_{\mathrm{los}}=50$ random line-of-sight viewing directions. The profiles
are azimuthally (circularly) averaged and range from $R=0.01$ kpc in the inner
parts to $R=30$ kpc in the outer parts. The surface density in the outer parts
($R\gtrsim 1$ kpc) increases with each merger generation in the all simulation
samples.

For comparison to our simulation results, we show the surface density profiles
of high redshift ($1.5<z<2.5$) quiescent galaxies (Szomoru et al., 2013),
local ETGs in the Virgo cluster with masses comparable or higher than for
observed $z=2$ compact ETGs (Kormendy & Bender, 2009; Szomoru et al., 2012)
and a number of local massive cored ETGs (Rusli et al., 2013). The central
parts of the high-redshift quiescent galaxies are not resolved at radii below
$r\lesssim 0.7\ \rm kpc$ (Szomoru et al., 2012), and their surface density
profiles are consistent both with our initial compact spherical setups and the
disky setups after a single major or minor merger.

####  4.1.1 Spherical models including a major merger

The spherical initial galaxy model host-S has a central (at $R=0.01$ kpc)
stellar surface density of $\Sigma_{\star}\sim{1.3\times
10^{5}}\>\mathrm{M_{\odot}\>pc^{-2}}$. In the spherical setups beginning with
a major merger (S-MA-SROT-BH and S-MA-FROT-BH) the central surface density
monotonically decreases with each merger generation as the SMBH binaries scour
increasingly extended, flat low-density cores into the galaxies.

SMBH binaries in the S-MA-FROT-BH sample have on average lower eccentricities
when they form compared to their counterparts with slowly rotating host
galaxies S-MA-SROT-BH as shown in Fig. 2. Thus, they enter the regime
dominated by GW losses at smaller spatial separations, having lost comparably
more orbital energy in the three-body slingshot hardening phase, resulting in
larger and less dense cores. After one major merger and eight minor mergers,
the central surface density of the S-MA-FROT-BH setup has decreased by two
orders of magnitude.

On the other hand, in the spherical sequences with a major merger but without
SMBHs, the central surface density gradually increases with each merger
generation, as opposed to their counterparts which include SMBHs. This is due
to violent relaxation (e.g. Lynden-Bell 1967; Hilz et al. 2012; Hilz et al.
2013). The setups S-MA-SROT and S-MA-FROT have very similar mildly cuspy
central structure after the minor mergers with $\Sigma_{\star}\sim{5\times
10^{5}}\>\mathrm{M_{\odot}\>pc^{-2}}$ at $R=0.01$ kpc.

The final merger remnants of the simulation sequences S-MA-SROT-BH and S-MA-
FROT-BH are in very good agreement with the observed local massive ETGs (e.g.
Rusli et al. 2013) with an extended, flat low-density core. Without SMBHs the
merger remnants have too high central surface densities compared to the
observed core ETGs and instead rather resemble cuspy power-law galaxies
without cores (Ferrarese et al., 1994).

####  4.1.2 Spherical models with only minor mergers

As shown in Fig. 2, in the spherical sequence S-MI-BH including only minor
mergers the formed SMBH binaries typically always have high eccentricities
($e\gtrsim 0.9$) and thus merge rapidly driven by gravitational wave emission.
Interestingly, the central surface density profile remains very close to the
initial mildly cuspy profiles despite the SMBH mergers as only relatively
little energy is transferred from the binary to the stellar component. We
performed additional SMBH merger simulations in which the SMBH binaries
initially have circular orbits in the setup. In these runs, presented in
Appendix B.1, a low-density core was scoured.

Despite the relatively unaffected central surface density profile
$\Sigma_{\star}(R)$ the comparison of the simulation sample S-MI-BH to its
counterpart S-MI without SMBHs clearly demonstrates that SMBH core scouring
has an effect on the profile. Without SMBHs, the central surface density
increases by each merger, reaching $\Sigma_{\star}\sim{1\times
10^{6}}\>\mathrm{M_{\odot}\>pc^{-2}}$ at $R=0.01$ kpc after the sequence of
seven minor mergers. Thus, in the sample S-MI-BH the SMBHs do not scour a
flat, extended low-density core, but rather prevent the central surface
density from increasing by a factor of $\sim 7$ during the minor merger
sequence. As we define the mass deficit as the missing central mass in the
runs with SMBHs compared to their counterparts without SMBHs, it is possible
to have $M_{\mathrm{def}}>0$ even if the central $\Sigma_{\star}(R)$ would
remain unaffected.

Compared to local ETGs, the models in the sample S-MI-BH lack a flat, low-
density core, even though their central surface densities are lower than in
the corresponding sample without SMBHs. Simulated galaxy merger remnants in
both sequences S-MA-BH and S-MI resemble mildly cuspy ETGs with their
Hernquist-like centers.

####  4.1.3 Disky models

The disky progenitor model host-D has a central ($R=0.01$ pc) surface density
of $\Sigma_{\star}\sim{1.9\times 10^{4}}\>\mathrm{M_{\odot}\>pc^{-2}}$, a
factor of $\sim 7.5$ lower compared to the spherical progenitor host-S. In the
sequences D-MA and D-MI without SMBHs, the central surface density rapidly
increases via minor mergers as the in-falling spherical satellites with a
Hernquist density profile fill the initially close to flat exponential center
of the disky model. After the merger sequences, the central surface densities
have reached values of $\Sigma_{\star}\sim{4.4\times
10^{5}}\>\mathrm{M_{\odot}\>pc^{-2}}$ (D-MA) and $\Sigma_{\star}\sim{5.5\times
10^{5}}\>\mathrm{M_{\odot}\>pc^{-2}}$ (D-MI).

In the disky models including SMBHs, the central surface density initially
increases in the first few mergers, both in the sequences with (D-MA-BH) and
without (D-MI-BH) a major merger. As with the model S-MI, SMBH core scouring
still acts in the simulations as the central surface densities without SMBHs
are systematically higher. As the later satellites are less bound (larger
$r_{\mathrm{h}}$) than the early ones as described in Table 2, core scouring
can lead to increasingly lowered central density in later mergers. The surface
density of the disky progenitor is again reached after $N_{\mathrm{merg}}\sim
4$-$5$ minor mergers, and the final central surface densities after the entire
merger sample are $\Sigma_{\star}\sim{9.6\times
10^{3}}\>\mathrm{M_{\odot}\>pc^{-2}}$ and $\Sigma_{\star}\sim{4.4\times
10^{3}}\>\mathrm{M_{\odot}\>pc^{-2}}$ for the samples D-MA-BH and D-MI-BH,
respectively. The final central $\Sigma_{\star}$ is very similar in the merger
samples D-MI-BH and S-MA-FROT-BH despite their different (spherical or disky)
progenitor properties.

The initial outer surface density profile of the exponential disky progenitor
is somewhat steep compared to the Sérsic outer parts of observed ETGs. A
single minor or major merger is enough to transform the outer simulated
profile into a Sérsic-like mildly curving shape. In the inner parts, only
models with SMBHs resemble the flat, low-density central parts of the local,
massive cored ETGs. Simulated merger remnants without SMBHs on the other hand
resemble cuspy power-law galaxies due to violent relaxation, even though the
slope of the central stellar cusp is not very steep.

###  4.2 Mass deficits with and without SMBHs

Figure 6: The mass deficits $M_{\mathrm{def}}/M_{\bullet}$ relative to the SMBH masses $M_{\bullet}$ after $N_{\mathrm{merg}}$ merger generations. sequence | $M_{\bullet}$ | $M_{\mathrm{def}}$ | $M_{\mathrm{def}}/M_{\star}$ | $\kappa$
---|---|---|---|---
| [${10^{9}}\>\mathrm{M_{\odot}}$] | [${10^{10}}\>\mathrm{M_{\odot}}$] |  |
Merritt (2006) |  |  |  | 0.5
S-MA-SROT | $6.8$ | $5.0$ | $0.16$ | 1.17
S-MA-FROT | $6.8$ | $7.1$ | $0.22$ | 1.99
S-MI | $4.8$ | $2.1$ | $0.09$ | 0.75
D-MA | $6.8$ | $7.0$ | $0.22$ | 1.74
D-MI | $4.8$ | $4.2$ | $0.17$ | 1.27
Table 6: The mass deficits $M_{\mathrm{def}}$ for the five merger sequences
compared to the final SMBH mass $M_{\bullet}$ and relative to the final
stellar mass $M_{\star}$ after the minor merger generations. The constant
$\kappa$ for the mass deficits relative to $M_{\bullet}$ is defined as
$M_{\mathrm{def}}=\kappa N_{\mathrm{merg}}M_{\bullet}$ from Fig. 6. We
consistently find higher values for $\kappa$ up to $\kappa\sim 2.0$ compared
to the value $\kappa=0.5$ of Merritt (2006).

For the purposes of our analysis, we use a practical definition to determine
the mass deficits $M_{\mathrm{def}}$ of our simulated merger remnants
including SMBHs. First, we calculate the influence radius $r_{\mathrm{infl}}$
of the SMBH using the cumulative stellar mass profile $M_{\star}(r)$ as
$M_{\star}(r<r_{\mathrm{infl}})=M_{\bullet}$. Next, we define the mass deficit
$M_{\mathrm{def}}$ as the difference of the stellar mass within $r=3\times
r_{\mathrm{infl}}$ from the center of the galaxy in the simulation with and
without SMBHs. Thus, our definition of the mass deficit $M_{\mathrm{def}}$ is

$M_{\mathrm{def}}=M_{\mathrm{\star}}^{\mathrm{noBH}}(r<3\times
r_{\mathrm{infl}})-M_{\mathrm{\star}}^{\mathrm{BH}}(r<3\times
r_{\mathrm{infl}})$ (1)

in which $M_{\mathrm{\star}}^{\mathrm{noBH}}(r)$ is the cumulative stellar
mass profile of a simulated galaxy without SMBHs and
$M_{\mathrm{\star}}^{\mathrm{BH}}(r)$ is the stellar mass profile of a
corresponding simulation including SMBHs.

Work based on numerical simulations of sinking binary SMBHs in stellar bulges
has predicted that mass deficit of a simulated core galaxy should correlate
with the final SMBH mass $M_{\bullet}$ and the total number of scouring events
the galaxy has experienced (Merritt, 2006). The mass deficit
$M_{\mathrm{def}}$ can be thus expressed as

$M_{\mathrm{def}}=\kappa N_{\mathrm{merg}}M_{\bullet}.$ (2)

Merritt (2006) predicted that this $\kappa$ to be a constant ($\kappa=0.5$)
from idealized simulations including a single, spherical stellar bulge with
two SMBHs with the mass deficit defined using the difference between the
initial and final density profiles.

We calculate the mass deficits from the simulated merger remnants of our five
simulation sample pairs S-MA-SROT, S-MA-FROT, S-MI, D-MA and D-MI. The results
are displayed in Fig. 6 with the slopes $\kappa$ of the individual merger
sequence pairs are listed in Table 6. The largest mass deficits (compared to
$M_{\bullet}$) occur in samples S-MA-FROT, D-MA and D-MI with
$M_{\mathrm{def}}/M_{\mathrm{\bullet}}\gtrsim 10$ after $N_{\mathrm{merg}}=7$,
as is already evident from the surface density profiles in Fig. 5. The large
mass deficits in the initially disky setups D-MA and D-MI originate from the
fact that the comparison simulation models without SMBHs acquire high central
stellar densities as the cuspy satellites fill the almost flat exponential
centers of the models in the mergers. The spherical samples with very
eccentric binaries (S-MA-SROT and S-MI) have the smallest mass deficits as
more SMBH binary orbital energy ends up being emitted in gravitational wave
radiation rather than transforming the central stellar structure of the
galaxy. For comparison, the less eccentric SMBH binaries in the spherical
sample S-MA-FROT acquire the largest mass deficits of this study as more SMBH
binary orbital energy couples to the galaxy rather than escaping via GWs. The
mass deficits relative to the stellar masses can also reach relatively high
values with high-mass SMBHs and a large number of mergers. In our simulations
we find $0.09\leq M_{\mathrm{def}}/M_{\star}\leq 0.22$. Thus, up to
approximately $\sim 10$–$20\%$ of the stellar material of the galaxy can be
displaced by SMBH core scouring in repeated galaxy mergers.

We find consistently larger values for $\kappa$ than Merritt’s value of
$\kappa=0.5$, ranging from $\kappa\sim 0.75$ (S-MI) to $\kappa\sim 2.0$ (S-MA-
FROT). We have additionally tested a simplified setup similar to Merritt
(2006) and find $\kappa\sim 0.5$ for initially circular binaries in the
context of the simple setup. Compared to a multi-component galaxy merger with
SMBHs, the setup of Merritt (2006) consists of two SMBHs on a circular orbit
at the center of a spherical stellar bulge. Our results suggest that in more
realistic merger sequence simulations mass deficits can be up to a factor of
$\sim 4$ larger compared to the idealized core scouring experiments found in
the literature. Similar large mass deficits have also been recently reported
by Partmann et al. (2023).

###  4.3 Core-Sérsic profiles

We fit the core-Sérsic profile (Graham et al., 2003; Trujillo et al., 2004) to
the surface density profiles of our simulated galaxies in Fig. 5 to quantify
the evolution of their central parts. Below the break radius or core size
$R_{\mathrm{b}}$, the core-Sérsic profile replaces the inner parts of the
common Sérsic profile (Sérsic, 1963; Caon et al., 1993) with a typically
almost flat central power-law component with a slope of $\gamma$. Above the
break radius, the Sérsic profile remains unchanged. In the terms of surface
density $\Sigma(R)$ the core-Sérsic profile can be written as

$\Sigma(R)=\Sigma^{\prime}\left[1+\left(\frac{R_{\mathrm{b}}}{R}\right)^{\mathrm{\alpha}}\right]^{\mathrm{\gamma/\alpha}}e^{-\mathrm{b}\left[\left(R+R_{\mathrm{b}}^{\mathrm{\alpha}}\right)/R_{\mathrm{e}}^{\mathrm{\alpha}}\right]^{\mathrm{1/(\alpha
n)}}}$ (3)

in which $n$ is the Sérsic index and $\Sigma^{\prime}$ is the overall
normalization related to the surface density at the break radius. The
parameter $\alpha$ determines the sharpness of the transition between the core
and the Sérsic parts of the profile. The constant $b$ is a somewhat complex
function of the core-Sérsic parameters defined in such a way that
$R_{\mathrm{e}}$ will enclose half of the total mass of the profile (Trujillo
et al., 2004).

As in Rantala et al. (2018), we use the cmpfit444cmpfit:
https://pages.physics.wisc.edu/$\sim$craigm/idl/cmpfit.html package to perform
the core-Sérsic fits. The radial fitting range for the core-Sérsic profiles is
$0.03$ kpc $<R<30.0$ kpc.

###  4.4 Observed vs simulated 2D mass deficits

Figure 7: The 2D mass deficits $M_{\mathrm{def,2D}}$ of the simulated merger
remnants as a function of the SMBH mass $M_{\bullet}$ in the merger sequences
compared to the mass deficits observed in a sample of local core ETGs (Rusli
et al., 2013). The evolution of the 2D mass deficits is steeper than the
linear relation $M_{\mathrm{def,2D}}\propto M_{\bullet}$.

Observationally, core light deficits are determined by measuring the
difference between the actually observed light profile $I(R)$ (or a fit to it)
and an inward extrapolation of the Sérsic profile $I_{\mathrm{Sersic}}$ that
describes the main body of the galaxy (Kormendy et al., 2009; Rusli et al.,
2013). These light deficits can be transformed into mass deficits
$M_{\mathrm{def}}$ with an appropriate stellar mass-to-light ratio
$\Upsilon=M/L$:

$M_{\mathrm{def,obs}}=2\pi\Upsilon\int_{\mathrm{0}}^{\mathrm{\infty}}dRR\left[I_{\mathrm{Sersic}}(R)-I(R)\right].$
(4)

This Sérsic-based approach would be desirable for determining the mass
deficits of simulated core galaxy models for a consistent comparison with
actual observed local ETG mass deficits. However, with increasing number of
merger generations the surface density profiles $\Sigma(R)$ in the outer parts
approach a power-law, making simple Sérsic fits with reasonable Sérsic indices
$n<10$ difficult. Hence, a consistent Sérsic-based mass deficit determination
from the simulated galaxies is not straightforward.

In general, steep power-law outer parts result in unreasonably large (and even
diverging) mass deficits when extrapolating the $n>10$ core-Sérsic fits to the
inner parts of the models. In principle, the issue might be avoided by
considering non-Hernquist initial models for the outer parts, exploring
different merger orbit geometries, or by performing careful multi-component
fitting (e.g. section 4 of Rusli et al. 2013). We note that also the observed
local relic galaxies require multiple Sérsic components for good surface
brightness profile fits (Yıldırım et al., 2015). However, here we instead
adopt a simple approach for determining the mass deficits from the surface
density profiles following section 4.2. Namely, we define the two-dimensional
mass deficit as $M_{\mathrm{def,2D}}$

$M_{\mathrm{def,2D}}=2\pi\int_{\mathrm{0}}^{\mathrm{3R_{\mathrm{b}}}}dRR\left[\Sigma_{\mathrm{\bullet}}(R)-\Sigma(R)\right]$
(5)

in which $\Sigma_{\mathrm{\bullet}}(R)$ is the scoured model and $\Sigma(R)$ a
corresponding simulation profile without SMBHs. We note that this mass deficit
determination is not fully consistent with observational mass deficit
determination procedures as observers have to rely on a single, scoured core-
Sérsic profile for measuring $M_{\mathrm{def,obs}}$. Nevertheless, as shown in
Fig. 7, the 2D mass deficits $M_{\mathrm{def,2D}}$ are in the range of
observed mass deficits in local ETGs (Rusli et al., 2013).

The mass deficits $M_{\mathrm{def,2D}}$ as function of SMBH mass $M_{\bullet}$
in different simulation samples behave overall in very similar manner and
scale steeper than $M_{\mathrm{def,2D}}\propto M_{\bullet}$. After the merger
sequences the 2D mass deficits are in the range $2.2\lesssim
M_{\mathrm{def,2D}}/M_{\bullet}\lesssim 3.1$ in the samples S-MA-SROT, S-MA-
FROT, D-MA and D-MI, consistent with the mass deficits of the local core ETG
population (Hopkins & Hernquist, 2010).

Here the prediction of Merritt (2006) based on sinking SMBH binaries in
isolated stellar bulges would yield $N_{\mathrm{merg}}\sim 4$–$6$, a somewhat
lower number than the true value. The merger sample S-MI that is not scouring
a low-density core has $M_{\mathrm{def,2D}}/M_{\bullet}\sim 1.3$. The good
agreement with the simulated 2D and observed mass deficits highlights that
observed mass deficits reliably describe the amount of stellar mass displaced
from the galactic nuclei due to SMBH binary core scouring even though the true
non-scoured profile is observationally unavailable. On the other hand, the
good agreement also strengthens the argument for the SMBH binary core scouring
origins for the diffuse low-density centers of massive ETGs.

###  4.5 Core scaling relations

Figure 8: The break radii $R_{\mathrm{b}}$ of the core-Sérsic fits of the
surface density profiles of the simulated merger remnants. The left panel
shows the setups with a spherical progenitors (S-MA-SROT, S-MA-FROT, S-MI)
while the right panel displays the models with disky progenitors (D-MA, D-MI).
The sequences including SMBHs are marked using solid lines while the runs
without SMBHs are shown as dotted lines. The shaded region in the left panel
indicates the central profiles with Hernquist-like shallow power-law cusps
which do not have a distinct flat core despite formally having non-zero break
radii $R_{\mathrm{b}}$ in the core-Sérsic fits. Only models with SMBHs can
grow central flat cores up to sizes of $\sim 1$ kpc. The sequence S-MI does
not produce a flat core with its very eccentric SMBH binaries as discussed in
the main text. Without SMBHs, the almost flat exponential central regions of
the initially disky models are filled by stellar material from the in-falling
satellites, and their somewhat core-like central initial appearance is lost.
Figure 9: The core sizes $R_{\mathrm{b}}$ as a function of the SMBH mass
$M_{\bullet}$ of the simulated galaxy merger samples compared to the local
observed relation of Thomas et al. (2016). Selected local cored ETGs (Rusli et
al., 2013; Thomas et al., 2016) are shown as open circles. Models with SMBHs
(again excluding the model S-MI-BH) can grow large cores up to
$R_{\mathrm{b}}\sim 1$ kpc.

####  4.5.1 Evolution of the core sizes $R_{\mathrm{b}}$

We show the evolution of the core sizes $R_{\mathrm{b}}$ in the galaxy models
in our major and minor merger sequences in Fig. 8 as a function of the stellar
mass after each merger generation. We first note that the projected Hernquist
profile has a shallow power-law center ($\Sigma(R)\propto R^{-\gamma}$ with
$\gamma<1$) such that the core-Sérsic fit obtains a core size
$R_{\mathrm{b,H}}\sim 0.1$ kpc for our isolated initial setups. Hence, we will
not consider $R_{\mathrm{b}}<R_{\mathrm{b,H}}=0.1$ kpc as actual cores in this
study.

For the spherical models without SMBHs (S-MA-SROT, S-MA-FROT, S-MI) the break
radius $R_{\mathrm{b}}$ never exceeds $R_{\mathrm{b,H}}=0.1$ kpc in the merger
sequences. In addition, their $R_{\mathrm{b}}$ always monotonically decrease
starting with the third minor merger generation. Thus, without SMBHs, the
initially spherical galaxies will not grow cores in the merger sequences and
retain their mildly cuspy central structure.

The disky models D-MA and D-MI originate from the exponential progenitor setup
host-D with an exponential surface density profile which by definition does
not have a core in the core-Sérsic sense as there is no break in the profile.
After a single major merger without SMBHs, the setup D-MA has a small core
with $R_{\mathrm{b}}\sim 0.12$ kpc. Similarly, in the setup D-MI without SMBHs
the simulated galaxy has a break radius of $R_{\mathrm{b}}\sim 0.22$ kpc after
the first minor merger. However, in the absence of SMBH core scouring in the
subsequent mergers, the small cores are filled and $R_{\mathrm{b}}$ steadily
decreases towards the later merger generations. After the entire merger
sequence the two models have very small core sizes of $R_{\mathrm{b}}\sim
0.059$ kpc (D-MA) and $R_{\mathrm{b}}\sim 0.035$ kpc (D-MI), respectively.
Even though the disky progenitor galaxies have almost flat core-like
exponential centers, the core-like central feature is erased by violent
relaxation in repeated galaxy mergers if no SMBHs are present.

With SMBHs, the model S-MI-BH does not form a core (even though
$M_{\mathrm{def}}>0$) as discussed in the previous sections. The break radius
of the model is always close to $R_{\mathrm{b}}\sim R_{\mathrm{b,H}}=0.1$ kpc.
This is not surprising as the surface density profiles of the model after each
merger generation remain practically unchanged, as seen already in Fig. 5.

With SMBHs, the spherical sequences S-MA-SROT-BH and S-MA-FROT-BH including a
major merger grow large, diffuse, flat central cores reaching up to kpc in
size. After seven minor mergers, the galaxy models have core sizes of
$R_{\mathrm{b}}\sim 0.58$ kpc (S-MA-SROT-BH) and $R_{\mathrm{b}}\sim 1.34$ kpc
(S-MA-FROT-BH), respectively, consistent with massive ETGs with large cores in
the local Universe (e.g. Lauer et al. 2007; Thomas et al. 2016). Reaching the
sizes of the largest observed galaxy cores ($R_{\mathrm{b}}\sim$ a few kpc)
potentially requires a subsequent merger of two already cored galaxies
(Mehrgan et al., 2019; Dullo, 2019).

With SMBHs, the disky setups D-MA-BH and D-MI-BH initially have
$R_{\mathrm{b}}\sim R_{\mathrm{b,H}}$ after their first mergers. In the
subsequent mergers including SMBHs, the break radius increases with each
merger generation, reaching $R_{\mathrm{b}}\sim 0.73$ kpc (D-MA-BH) and
$R_{\mathrm{b}}\sim 0.62$ kpc (D-MI-BH) after seven minor mergers. The core
sizes are very similar compared to especially the spherical model S-MA-SROT-
BH.

####  4.5.2 The $M_{\bullet}$–$R_{\mathrm{b}}$ relation

We present the evolution of the simulated galaxies in the
$M_{\bullet}$–$R_{\mathrm{b}}$ plane in Fig. 9. For comparison to the
simulation data, we display the core sizes and SMBH masses of a number of
observed local ETGs (Rusli et al., 2013; Thomas et al., 2016). In addition, we
show the observed relation between the core sizes $R_{\mathrm{b}}$ and the
SMBH masses $M_{\bullet}$ which can be expressed as

$\log_{\mathrm{10}}\left(\frac{M_{\bullet}}{M_{\odot}}\right)=a+b\times\log_{\mathrm{10}}\left(\frac{R_{\mathrm{b}}}{\text{kpc}}\right).$
(6)

For the coefficients $a$ and $b$ of the relation Thomas et al. (2016) provides
values $a=10.27\pm 0.51$ and $b=1.17\pm 0.14$.

The observed relation between the core size and the SMBH mass is not
especially tight, and all models in which SMBHs scour cores (S-MA-SROT-BH,
S-MA-FROT-BH, D-MA-BH, D-MI-BH) are consistent with Thomas et al. (2016). For
the models S-MA-BH the small break radii $R_{\mathrm{b}}\sim R_{\mathrm{b,H}}$
are consistent with the mildly cuspy Hernquist centers and are not considered
as actual cores. All merger sequences forming cores have their break radii
$R_{\mathrm{b}}$ for a given SMBH mass $M_{\bullet}$ initially lying below the
observed relation. Importantly, in the simulated merger sequences with cores
the core size $R_{\mathrm{b}}$ grows steeper as a function of $M_{\bullet}$
than the observed relation. Thus, each sequence with a core will eventually
cross the relation, and finally be far above it after a large number of minor
mergers. In our simulations, the galaxy models reach the observed
$M_{\bullet}$–$R_{\mathrm{b}}$ relation after $3$–$5$ minor mergers. Only the
model S-MA-FROT-BH grows a core large compared to the observed relation, this
occurs after $6$–$7$ minor mergers.

We note that the evolution of the simulated galaxies in the
$M_{\bullet}$–$R_{\mathrm{b}}$ plane depends on the slope of the initial
density profile. In general, the progenitors of local massive ETGs are thought
to have density profiles steeper than the Hernquist profile up to $\gamma\sim
2.0$. The Hernquist model ($\gamma=1.0$) adopted for this study primarily for
computational reasons results in more rapid increase of $R_{\mathrm{b}}$ as a
function of $M_{\bullet}$ compared to initially cuspier profiles such as the
$\gamma=3/2$ model as shown by Rantala et al. (2018) in their Fig. 11. The
cuspy $\gamma=3/2$ models move closely along the observed local
$M_{\bullet}$–$R_{\mathrm{b}}$ relation after each merger, indicating no
redshift evolution of the core-SMBH relations for such models.

Our results indicate that it is possible for a $z=2$ compact nugget-like ETGs
to reach the $z=0$ mass-size relation and simultaneously form a core
consistent with the observed $M_{\bullet}$–$R_{\mathrm{b}}$ relation after
$6$–$7$ minor mergers. The required number of mergers is consistent with
empirical galaxy formation models which predict that the minor merger rate
hardly exceeds $\sim 1$ Gyr-1 at $z\lesssim 2$ (O’Leary et al., 2021). The
presence of SMBHs is crucial for the scenario as without SMBHs the galaxy
merger remnants do not form cores, even when starting from core-like
exponential disky setups.

##  5 Core scouring and the global structure of ETGs

###  5.1 Large mass deficits can increase the effective radii of galaxies

Figure 10: The relative increase of the effective radii
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}>1$ of the galaxies in the merger
sequences as a function of increasing stellar mass (merger generation)
comparing mergers with and without SMBHs. The additional size growth due to
SMBHs is especially prominent in merger samples S-MA-FROT and D-MA, reaching
values up to $R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.76$ after a single
major merger and seven minor mergers.

In Section 3.1 we pointed out that merger sequences involving SMBHs grow their
effective radius faster in the minor merger sequences than their counterparts
without SMBHs, i.e. $R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}>1$. In other words,
for the size growth relation $R_{\mathrm{e}}\propto M_{\star}^{\alpha}$, the
power-law index $\alpha$ is larger in simulations including SMBHs, thus
$\alpha_{\bullet}-\alpha>0$.

In Fig. 10 we show the increase of the galaxy size
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}$ in the simulations with SMBHs relative
to the simulations without SMBHs as a function of the increasing galaxy
stellar mass in the merger simulations. In the three major merger simulations
of this study the extra size growth due to SMBHs is mild with
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\lesssim 1.09$. However, in the minor
merger sequences the additional size growth can be substantial, especially in
the setups S-MA-FROT and D-MA, with the ratio of the effective radii reaching
values of up to $R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.76$. The effect
is moderate in the sequences S-MA-SROT
($R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.29$) and D-MI
($R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.34$). The sample S-MI which does
not form a core with SMBHs even though having $M_{\mathrm{def}}>0$ has the
smallest relative effective radius size increase of
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.12$ after seven minor mergers.

Figure 11: The correlation of the mass deficit $M_{\mathrm{def}}/M_{\bullet}$
and the additional galaxy size growth as measured by the growth index
difference $\alpha_{\bullet}-\alpha$ (left panel) and the ratio of the
effective radii $R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}$ (right panel). The
dashed lines in each panel show a simple linear fit between the quantities
with the fit coefficients are listed in the main text. The displayed results
are after seven generations of minor mergers.

Our working hypothesis is that the stellar material displaced by the SMBH core
scouring is responsible for the observed size growth effect. We perform a
simple test in Appendix B.2 in which we show that the increase of galaxy size
indeed stops if the SMBH binary is merged in the simulation by hand. If the
core scouring is not responsible for the observed galaxy size growth in the
merger sequences there should be no correlation between the mass deficits
$M_{\mathrm{def}}$ and the size growth diagnostics
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}$ and $\alpha_{\bullet}-\alpha$.
However, we find such correlations between $M_{\mathrm{def}}$,
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}$ and $\alpha_{\bullet}-\alpha$.

We present the size growth - mass deficit correlation results in Fig. 11. Both
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}$ and $\alpha_{\bullet}-\alpha$ scale
linearly with the final mass deficits of the galaxy models from Table 6. The
best-fit coefficients for the linear relations with (95% confidence intervals)
are

$\alpha_{\bullet}-\alpha=a+b\times\frac{M_{\mathrm{def}}}{M_{\bullet}}$ (7)

with $a=0.71\pm 0.15$ and $b=0.53\pm 0.10$ for the size growth index
difference and

$\frac{R_{\mathrm{e,\bullet}}}{R_{\mathrm{e}}}=c+d\times\frac{M_{\mathrm{def}}}{M_{\bullet}}$
(8)

with $c=-0.41\pm 0.21$ and $d=0.74\pm 0.14$ for the effective radius ratio. We
proceed to explain this linear relationship between the increase of
$R_{\mathrm{e}}$ and the mass deficit $M_{\mathrm{def}}$ in the next section.

###  5.2 A simple analytic model for size evolution due to SMBH scouring

In the previous section we established the linear dependence of the additional
effective radius growth $R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}$ in simulations
including SMBHs on the mass deficits $M_{\mathrm{def}}$. We present a simple
scouring model to qualitatively explain why the mass deficit can increase the
effective radius of a galaxy and why the dependence is linear. The simple
scouring model can be summarized with the following points.

  * •

A galaxy with a total mass $M$, a cumulative mass profile $M(r)$ and an
initial effective radius $R_{\mathrm{e}}$ corresponding to a half-mass radius
$r_{\mathrm{h}}$ is scoured by one or more SMBH binaries with the final mass
deficit of $M_{\mathrm{def}}$. Stellar mass is removed from the nucleus of the
galaxy in such a manner that the central density profile $\rho(r)$ is flat
within the core radius.

  * •

Case one: scoured stars remain bound. The material scoured from the core is
deposited in the outer parts ($r>r_{\mathrm{h}}$) of the galaxy. In this case,
a mass of $M_{\mathrm{def}}$ is removed from inside the original
$r_{\mathrm{h}}$ and added outside it. The new half-mass radius
$r_{\mathrm{h,new}}$ is determined by the criterion
$M(r_{\mathrm{h,new}}<r)=M/2+M_{\mathrm{def}}$.

  * •

Case two: scoured stars get unbound. The scoured material becomes completely
unbound of the galaxy, and the mass of the galaxy is decreased by
$M_{\mathrm{def}}$. The material ejected from within $r_{\mathrm{h}}$ is
permanently removed, and the new half-mass radius is calculated as
$M(r_{\mathrm{h,new}}<r)=M/2+M_{\mathrm{def}}/2$.

The two cases can be combined into a single one by introducing a constant
$\eta$ which has a value of $\eta=1$ if the scoured material is retained in
the outer parts of the galaxy (case one) or $\eta=1/2$ if the material becomes
unbound (case two). Denoting the increase of the half-mass radius by $\Delta
r_{\mathrm{h}}=r_{\mathrm{h,new}}-r_{\mathrm{h}}$, the relative change of the
half-mass radius of the galaxy due to core scouring can be written as

$\frac{\Delta R_{\mathrm{e}}}{R_{\mathrm{e}}}=\frac{\Delta
r_{\mathrm{h}}}{r_{\mathrm{h}}}=\frac{\eta
M_{\mathrm{def}}}{r_{\mathrm{h}}}\left[\left.\frac{\mathrm{d}M(r)}{\mathrm{d}r}\right|_{\mathrm{r=r_{\mathrm{h}}}}\right]^{-1}.$
(9)

Most importantly, the relative change depends linearly on the mass deficit
$M_{\mathrm{def}}$, and on the derivative of the cumulative mass profile
$M(r)$ at the original half-mass radius $r_{\mathrm{h}}$. This can be
intuitively understood as follows. If the density profile $\rho(r)$ declines
sharply near $r_{\mathrm{h}}$, the mass profile $M(r)$ at $r_{\mathrm{h}}$ is
also very flat. Thus, the half-mass radius must be considerably increased to
enclose an additional mass of $\eta M_{\mathrm{def}}$. The opposite also
holds: if the density profile is not very steep at $r_{\mathrm{h}}$, the mass
profile steeply increases there, and even a small increase in the half-mass
radius is enough to enclose the extra mass $\eta M_{\mathrm{def}}$.

Figure 12: A simple scouring model to illustrate the increase of the half-mass
radius of a galaxy for large mass deficits. We assume
$M_{\mathrm{def}}=2N_{\mathrm{merg}}M_{\bullet}$ based on the maximum mass
deficits in our simulations listed in Table 6. The half-mass radius of the
initial model is $r_{\mathrm{h}}=1.00$ kpc. We assume that the scoured
material remains bound and is redistributed in the outer parts of the galaxy,
so $\eta=1$ for Eq. (9). After $N_{\mathrm{merg}}=1$, $N_{\mathrm{merg}}=3$
and $N_{\mathrm{merg}}=7$ scouring events, the half-mass radius of the galaxy
model has reached $r_{\mathrm{h}}=1.14$ kpc, $r_{\mathrm{h}}=1.41$ kpc and
$r_{\mathrm{h}}=2.00$ kpc, respectively.

For mass profiles that are shallow near $r_{\mathrm{h}}$, even small mass
deficits can result in considerably large relative changes in the half-mass
radius due to the derivative term in Eq. (9). We demonstrate this in Fig. 12
using the simple scouring model. We assume a galaxy model closely resembling
our isolated progenitors: $r_{\mathrm{h}}=1$ kpc and
$M_{\star}={10^{11}}\>\mathrm{M_{\odot}}$. We assume that the SMBH of the
galaxy with a mass of $M_{\bullet}={10^{9}}\>\mathrm{M_{\odot}}$ was assembled
in $N_{\mathrm{merg}}=1$–$7$ earlier SMBH mergers, and that the inspirals have
removed a mass of $M_{\mathrm{def}}$ from the center of the galaxy. We
calculate the mass deficit using $M_{\mathrm{def}}=\kappa
N_{\mathrm{merg}}M_{\bullet}$ with $\kappa=2.0$ motivated by the mass deficits
in our simulation sequences listed in Table 6. For the fate of the scoured
stars we assume that they remain bound in the outer parts of the galaxy model,
i.e. $\eta=1$ in Eq. (9).

For the Hernquist profile Eq. (9) can be easily evaluated the result being

$\frac{\Delta
R_{\mathrm{e}}}{R_{\mathrm{e}}}=\left(2+\sqrt{2}\right)\frac{\eta
M_{\mathrm{def}}}{M_{\star}}\sim 3.41\times\frac{\eta
M_{\mathrm{def}}}{M_{\star}}.$ (10)

Inserting the maximum mass deficit $M_{\mathrm{def}}/M_{\star}=0.22$ from our
simulations (D-MA-BH) in Section 4.2, we have $\frac{\Delta
R_{\mathrm{e}}}{R_{\mathrm{e}}}\sim 0.75$ for $\eta=1.0$. This corresponds to
a value of $R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.75$ which is in
remarkably good agreement with the actual simulation result
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.76$. For the merger sequences
S-MA-SROT-BH, S-MA-FROT-BH and S-MI the simple scouring model predicts
$R_{\mathrm{e,\bullet}}/R_{\mathrm{e}}\sim 1.27$–$1.55$ (simulation result
$1.29$, S-MA-SROT-BH), $1.37$–$1.75$ ($1.67$, S-MA-FROT-BH), $1.29$–$1.58$
(1.33, D-MI-BH). The agreement between the simple model and the simulation
results is in general very good. The model somewhat overpredicts
($1.15$–$1.31$) the size increase for the sequence S-MI-BH ($1.12$) which did
not scour a flat, low-density core.

We emphasize that the spatial galaxy size increase due to mass displacement by
core scouring is in general only significant for SMBHs above the
$M_{\bullet}$–$M_{\star}$ relation, but can reach surprisingly large values
for over-massive SMBHs. For SMBHs near the relation ($M_{\bullet}\sim
1.4\times 10^{-3}M_{\star}$; Häring & Rix 2004) the relative size increase is

$\frac{R_{\mathrm{e,\bullet}}}{R_{\mathrm{e}}}\sim
1.0+0.01N_{\mathrm{merg}}\text{\hskip 85.35826pt(typical ETGs)}$ (11)

assuming the Hernquist profile, $\kappa=2$ for the mass deficit and $\eta=1$
for the scoured material. Thus, for typical ETGs the effect is at most at a
few percent level. However, for SMBHs above the $M_{\bullet}$–$M_{\star}$
relation the relative increase may reach larger values due to the linear
dependence of the size increase on the mass deficit $M_{\mathrm{def}}\propto
M_{\bullet}$. For SMBHs which lie a factor of $10$ above the
$M_{\bullet}$–$M_{\star}$ relation we have

$\frac{R_{\mathrm{e,\bullet}}}{R_{\mathrm{e}}}\sim
1.0+0.1N_{\mathrm{merg}}\text{\hskip 8.5359pt(SMBH by a factor of 10 over-
massive)}$ (12)

which can easily lead to galaxy size increases of tens of percents or more. An
over-massive SMBH by a factor of $20$ compared to the
$M_{\bullet}$–$M_{\star}$ relation only requires $N_{\mathrm{merg}}\sim 5$–$6$
past mergers for a size increase by a factor of $2$.

It is expected that the progenitors of present-day massive ETGs had central
density profiles steeper $\gamma$ than the Hernquist model. We generalize Eq.
(10) and show the results in Appendix C. The galaxy size increase due to core
scouring is enhanced in more cuspy models. However, the relative size increase
depends on the exact value of $\gamma$ only mildly in the range
$0\leq\gamma\lesssim 2.5$. For example, the Jaffe (1983) model ($\gamma=2$)
grows in size only $\sim 17\%$ more than the Hernquist model ($\gamma=1.0$).
Thus, the mass deficit $M_{\mathrm{def}}$ itself rather than the exact
progenitor profile is the main driver of the relative galaxy size increase due
to SMBH core scouring.

##  6 Shape and rotation properties of ETGs

###  6.1 Galaxy shapes

Figure 13: The spherical, axisymmetric or triaxial shapes of the initial
galaxy models (host-D, host-S) and merger remnants.

We characterize the three-dimensional shapes of the stellar components of our
galaxy models which include SMBHs after each merger generation. The axis
ratios $b/a$ and $c/a$ (with $a\geq b\geq c$) within $r_{\mathrm{h}}$ are
determined by calculating the shape (or inertia) tensor of the stellar
particle distribution and its eigenvalues (e.g. Jesseit et al. 2005; Binney &
Tremaine 2008) following the recipe of Zemp et al. (2007). The triaxiality
parameter $T$ can in turn be calculated from the axis ratios $b/a$ and $c/a$
as

$T=\frac{1-\left(b/a\right)^{2}}{1-\left(c/a\right)^{2}}.$ (13)

We show the shapes of our galaxy models after each merger generation in the
$b/a$–$c/a$ plane in Fig. 13. By construction, the spherical galaxy progenitor
host-S has $b/a=c/a=1$ and for the disky progenitor $b/a=1$ and $c/a=0.5$.

The model S-MA-SROT-BH has a shape of $b/a=0.81$, $c/a=0.75$ ($T=0.77$) after
its major merger. This merger remnant is the most triaxial galaxy of the
study. In the minor mergers, both $b/a$ and $c/a$ steadily increase for the
model up to $b/a=0.98$, $c/a=0.91$ ($T=0.20$). On the contrary, the setup
S-MA-FROT-BH is axisymmetric ($b/a=0.99$, $c/a=0.79$) rather than triaxial
after its major merger. In the subsequent minor mergers, the model becomes
somewhat more triaxial as $b/a$ decreases to a value of $b/a=0.93$ while $c/a$
changes even less with a final value of $c/a=0.83$ ($T=0.46$). The merger
remnants in the merger sequence S-MI-BH with spherical progenitors and without
a major merger are the most spherical galaxies of the simulations of this
study with $c/a>0.92$ for each remnant.

The models originating from disky progenitors (D-MA-BH, D-MI-BH) are all
axisymmetric as $b/a>0.95$ for all the merger remnants. The model D-MA-BH has
$c/a=0.75$ after the major merger and $c/a=0.89$ at the end of the minor
merger sequence. The setup D-MI-BH remains the most disk-like among the merger
remnants of the study with its axisymmetric shape with $c/a=0.74$ after the
minor mergers.

###  6.2 SMBH binary core scouring as a global process: the orbit structure of
scoured ETGs

Figure 14: The ratio of box orbits to tube orbits
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ as a function of radius, scaled to the
half-mass radius of each galaxy. The top row shows
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ for disk progenitors (left) and spherical
progenitors (right), both with a BH, and the bottom row shows
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ for disk progenitors (left) and spherical
progenitors (right), both without a BH. In all panels, the dashed line
represents $N_{\mathrm{box}}/N_{\mathrm{tube}}$ following the first merger of
the sequence, and the solid line $N_{\mathrm{box}}/N_{\mathrm{tube}}$
following the final merger of the sequence.

SMBH binaries preferentially eject stars on close-to-radial orbits with little
angular momentum. Stars with low angular momentum orbits do not spend all
their time within the central region, but instead spend substantial periods of
time further out in the galaxy. Thus, SMBH core scouring removing these orbits
can be in a sense viewed as a process globally changing the orbit structure of
the galaxy, not just altering the core within $R_{\mathrm{b}}$ (see Frigo et
al. 2021). In the previous section we demonstrated that core scouring by
overmassive SMBHs can alter the global properties of the host galaxies, such
as their effective radii $R_{\mathrm{e}}$. Next, we study this process in more
detail using orbit analysis techniques.

We follow the methodology presented in Frigo et al. (2021) to construct an
analytical potential in which we can integrate the stellar particles to
classify the particles into one of six orbital families. To determine which
family a particular stellar orbit belongs to, the frequency of the particle’s
motion is determined along each principal axis of the galaxy’s inertia tensor
through a Fourier transform. Briefly, these orbit families are:

  1. 1.

$\pi$-box orbits, which have no net angular momentum and no resonant motion
along the three principal axes,

  2. 2.

boxlet orbits, which also have no net angular momentum, but do show resonant
motion along the three principal axes,

  3. 3.

$x$-tube orbits, which rotate about the major axis of the galaxy,

  4. 4.

$z$-tube orbits, which rotate about the minor axis of the galaxy,

  5. 5.

rosette orbits, which are typically associated with a spherically-symmetric
potential and conserve each component of the angular momentum vector, and

  6. 6.

irregular orbits, where the particle is not bound by any integrals of motion.

Having classified the stellar particles into the above orbital families, we
group the families into two broad classes: orbits without angular momentum
($\pi$-box and boxlet orbits) which we term ‘boxes’, and orbits with some
degree of angular momentum (rosette, $x$-, and $z$-tube orbits) which we term
‘tubes’.

Particles are then binned radially into eight logarithmically-spaced bins
spanning $0.5r_{\mathrm{h}}$ to $12r_{\mathrm{h}}$, and we determine the ratio
of box orbits to tube orbits $N_{\mathrm{box}}/N_{\mathrm{tube}}$ for each
radial bin. This analysis is done for the first snapshot after the first
merger (this may be either a major or minor merger depending on the merger
sequence) and the first snapshot after the final merger (which by construction
is always a minor merger). We show the box-tube ratio in Fig. 14, with dashed
lines indicating the ratio following the first merger, and solid lines
indicating the ratio following the final merger. A general trend we observe is
that the simulated merger remnants tend to be significantly dominated by tube
orbits at small radii compared to box orbits following the final merger. We
discuss each simulated merger remnant in detail below.

Looking first at disky mergers with a BH, we see that for $r\leq
r_{\mathrm{h}}$, little change occurs for the ratio of box to tube orbits.
This is in agreement with the galaxy shapes in Fig. 13, where for both the
D-MA-BH and D-MI-BH sequences the $b/a$ fraction is consistent after the first
merger and after the final merger. For radii $r>r_{\mathrm{h}}$, both D-MI-BH
and D-MA-BH increase the fraction of box orbits during their merger sequence.
The increase in $N_{\mathrm{box}}/N_{\mathrm{tube}}$ beyond $r_{\mathrm{h}}$
can naturally be understood as the stellar particles that have interacted with
the SMBH being deposited at larger radii, as described in the analytical
scouring model in section 5. In particular, many tube orbits that interact
with a BH binary are converted to either radial box orbits, or unbound from
the system completely (Frigo et al., 2021). In both cases,
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ of the final merger remnant is raised
compared to the ratio following the first merger.

Turning our attention to the disky mergers without a BH, we see a distinct
difference to the case with a BH. Following the first merger, both D-MA and
D-MI are dominated by box orbits within $r_{\mathrm{h}}$, with
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ in excess of unity. Following the final
merger, the fraction of box orbits within $r_{\mathrm{h}}$ has decreased to
less than 0.1: more than an order of magnitude compared to the first merger.
In both instances, we find that the central stellar densities of the systems
are very high. Comparing the snapshot following the final merger in the D-MA
case with the equivalent D-MA-BH case, we see that the there is an equal
amount of mass within $\sim 10^{-2}\>\mathrm{kpc}$ between the two
simulations. In the D-MA-BH case, the majority of this mass is contributed by
the SMBH, whereas in the D-MA case the mass is due almost exclusively to
stellar particles. Additionally, in the D-MA case, there is more mass than in
the D-MA-BH case for radii $10^{-2}\>\mathrm{kpc}\lesssim r\lesssim
10\>\mathrm{kpc}$ due to violent relaxation. The high central stellar
densities provides the necessary gravitational potential to sustain a large
fraction of tube orbits, thus lowering the overall value of
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ for the D-MA and D-MI systems following
the final merger in their respective sequences.

We now consider the spherical progenitor mergers with SMBHs in the top right
panel of Fig. 14. Immediately apparent is that the S-MI-BH merger sequence
displays only minor variation in $N_{\mathrm{box}}/N_{\mathrm{tube}}$
following the first and final mergers, taking a value of $\sim 0.1$ for $r\leq
r_{\mathrm{h}}$, and increasingly mildly with radius to $\sim 0.2$ at the
outermost radii. There is a slight increase in
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ following the first merger as compared to
the final merger. The case is drastically different for S-MA-FROT-BH and S-MA-
SROT-BH merger sequences. S-MA-FROT-BH displays an evolution in
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ between mergers consistent with the disky
mergers, where the ratio increases by the final merger compared to the first
merger. By $4r_{\mathrm{h}}$, $N_{\mathrm{box}}/N_{\mathrm{tube}}$ has
exceeded unity, indicating that at large radii box orbits are the dominant
orbital family. Conversely, for the S-MA-SROT-BH sequence,
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ decreases from the first merger to the
final merger. Tube orbits always dominate over box orbits at all radii. The
peculiar change in ratio between the first and final mergers is reflected in
the galaxy shape plot of Fig. 13. For D-MA-BH, D-MI-BH, and S-MA-FROT-BH the
triaxiality parameter $T$ increases from the first merger to the final merger,
with a final value $0.1\lesssim T\lesssim 0.4$. For the S-MA-SROT-BH case, the
galaxy following the first merger is relatively prolate ($T\sim 0.7$), and
becomes increasingly oblate, with a final $T$ of $\sim 0.2$. Tube orbits, in
particular $z$-tube orbits with their regular rotation about the minor axis,
are responsible for a more oblate remnant.

We see similar trends for the spherical mergers without SMBHs as for their
SMBH merger counterparts, with the exception of the S-MA-FROT box-tube ratio
never exceeding unity: for all spherical mergers without a BH, tube orbits
dominate over box orbits at all radii. The magnitude of
$N_{\mathrm{box}}/N_{\mathrm{tube}}$ for radii within $r_{\mathrm{h}}$ is also
very similar for the spherical mergers without a SMBH as for the mergers with
a SMBH, with the exception of S-MA-FROT. In the S-MA-FROT-BH sequence, the
magnitude of the box-tube ratio is $\sim 30\%$, compared to $\sim 10\%$ in the
S-MA-FROT case.

###  6.3 Velocity anisotropy profiles

Figure 15: The orbital anisotropy parameter $\beta_{\mathrm{rot}}$ for the
simulated merger remnants after seven minor merger generations. The radial
coordinate $r$ on the x-axis is scaled with the two-dimensional break radii
$R_{\mathrm{b}}$ of the models, except for the sample S-MI-BH (dashed line)
which did not scour a flat central core. For the S-MI-BH sample we use here
$R_{\mathrm{b}}=1$ kpc for visualization purposes. The simulated merger
remnants are tangentially biased ($\beta_{\mathrm{rot}}<0$) within $r\lesssim
2R_{\mathrm{b}}$. The observed orbital anisotropies of cored, local ETGs
(Thomas et al., 2014) in the SINFONI BH survey (Saglia et al., 2016) are shown
as the shaded region.

Next we study the velocity anisotropy structure of our merger remnants after
the minor merger generations. The orbital anisotropy parameter
$\beta_{\mathrm{rot}}$ (Thomas et al., 2014) we use is defined as

$\begin{split}\beta_{\mathrm{rot}}&=1-\frac{\sigma_{\mathrm{rot}}^{2}}{2\sigma_{\mathrm{r}}^{2}}\\\
\sigma_{\mathrm{rot}}^{2}&=\sigma_{\mathrm{\theta}}^{2}+\sigma_{\mathrm{\phi}}^{2}+v_{\mathrm{\phi}}^{2}\end{split}$
(14)

using the $\phi$ component of the rotation velocity $v_{\mathrm{\phi}}$ and
the three components of the velocity dispersion $\sigma_{\mathrm{r}}$,
$\sigma_{\mathrm{\theta}}$ and $\sigma_{\mathrm{\phi}}$ in spherical
coordinates.

We show the orbital anisotropies $\beta_{\mathrm{rot}}$ of the final merger
remnants after seven generations of minor mergers including SMBHs in Fig. 15.
All the merger remnants are tangentially biased ($\beta_{\mathrm{rot}}<0$) at
separations $r\lesssim 2$–$3R_{\mathrm{b}}$ from the center due to the
prevalence of tube orbits at these radii. Perhaps surprisingly, the model
S-MI-BH which did not scour a proper flat core has the lowest central
anisotropy, $\beta_{\mathrm{rot}}\sim-1$. The models S-MA-SROT-BH, S-MA-FROT-
BH, D-MA and D-MI have similar $\beta_{\mathrm{rot}}(r)$ profiles and central
values of $-0.6\lesssim\beta_{\mathrm{rot}}\lesssim-0.3$. We compare these
results to galaxies from the SINFONI black hole survey (Saglia et al., 2016).
The orbital anisotropy profiles of observed galaxies (Thomas et al., 2014) are
very consistent with the profiles of the simulated merger remnants.

###  6.4 Stellar kinematics in two dimensions

Figure 16: Two-dimensional kinematic maps of the simulation sequence S-MA-
FROT-BH. The isophote coinciding with the effective radius $R_{\mathrm{e}}$ of
each merger remnant is highlighted. The vertical bars in the left panels mark
the extent of $3$ kpc. The spatial extent of each panel is six times the
effective radius $R_{\mathrm{e}}$ of the merger remnant. The various features
of the maps are discussed in the main text. Figure 17: Two-dimensional
kinematic maps of the simulation sequence D-MA-BH.

We show the two-dimensional kinematic maps of the simulation sequences S-MA-
FROT-BH and D-MA-BH in Fig. 16 and in Fig. 17. The two simulation samples were
chosen for the figures as they have the strongest rotation signatures among
all the runs of this study. The maps are constructed to mimic how the
simulated merger remnants would be seen if observed with an integral field
unit spectrograph (see e.g. Naab et al. 2014). The simulation particle
velocity data is first binned on a regular grid with a spatial resolution of
$0.03\times R_{\mathrm{e}}$ of each galaxy. The regular grid is then Voronoi
binned into spaxels with a constant signal-to-noise ratio, i.e. the number of
particles using a Voronoi tessellation algorithm (Cappellari & Copin, 2003).
Specifically, we use a version of the Voronoi algorithm in a modified version
of the pygad package (Röttgers et al., 2020).

The four panels of Fig. 16 and Fig. 17 show in each row, spaxel by spaxel, the
line-of-sight velocity $V_{\mathrm{avg}}$, the velocity dispersion $\sigma$
and the parameters $h_{\mathrm{3}}$ and $h_{\mathrm{4}}$ that measure the
skewness and kurtosis of the velocity distribution in the spaxels (Gerhard,
1993; van der Marel & Franx, 1993). In general, the anti-correlation of
$V_{\mathrm{avg}}$ and $h_{\mathrm{3}}$ points towards rotating, disky
structures while negative $h_{\mathrm{4}}$ is a proxy for a tangentially
biased velocity distribution with preferentially circular orbits. Finally, the
panels also show the two-dimensional flux (surface density) contours with
separations of one magnitude.

The kinematic maps of the merger remnants in the simulation sequence S-MA-
FROT-BH are shown in Fig. 16. After the major merger, the remnant is rapidly
rotating with $|V_{\mathrm{avg}}|\sim 240$ km/s. The region of fast rotation
extends beyond $R=2R_{\mathrm{e}}$. The peak $V_{\mathrm{avg}}$ decreases with
increasing number of minor mergers while the region of fast rotation shrinks
inside the increasing $R_{\mathrm{e}}$. After the minor merger generations
$|V_{\mathrm{avg}}|\lesssim 70$–$80$ km/s as the mergers from isotropic
directions bring in more mass but on average no net angular momentum. The
velocity dispersion maps show a single, central peak with $\sigma\sim 400$
km/s after the major merger. As the central stellar density is lowered due to
core scouring, the central velocity dispersion also decreases, reaching
$\sigma\sim 275$ km/s after the minor mergers. The $h_{\mathrm{3}}$ maps
anticorrelate with the corresponding $V_{\mathrm{avg}}$ maps. With increasing
number of minor mergers, the $h_{\mathrm{4}}$ is increasingly more negative in
the core regions and increasingly more positive outside it. Finally, the
isophotal shapes are in general ellipsoidal (but circular at the center) and
turn somewhat more boxy in the minor merger sequences. Overall, the evolution
of the isophotal shapes is relatively minor.

After the disk-disk major merger, the merger remnant D-MA-BH has a maximum
$|V_{\mathrm{avg}}|\sim 120$ km/s peaking at $R\sim 0.5R_{\mathrm{e}}$.
Outside the effective radius there is less rotation. In the minor merger
generations the rotation weakens, reaching $|V_{\mathrm{avg}}|\lesssim 50$
km/s after seven minor mergers. Similarly, the central velocity dispersion
decreases considerably from $\sigma\sim 500$ km/s to $\sigma\sim 300$ km/s as
the initially exponential model transforms into a cored ETG. Again, the
$h_{\mathrm{3}}$ and $V_{\mathrm{avg}}$ anticorrelate, and $h_{\mathrm{4}}$ is
increasingly negative within the core and positive outside it. The negative
central $h_{\mathrm{4}}$ appears earlier in the minor merger generations
compared to the sequence S-MA-FROT-BH. The isophotal shapes are initially
ellipsoidal, also near the center, and become increasingly round after each
minor merger generation.

We note that our simulated merger remnants lack kinematic misalignments
observed in a subset of local, massive ETGs (Franx et al., 1991; Krajnović et
al., 2011). This is most probably due to the fact that the satellite galaxies
infall from isotropic directions in the merger sequences, preventing the
formation of distinct kinematic features.

###  6.5 The slow and fast rotating ETGs in the merger sequences

Figure 18: The spin parameter $\lambda_{\mathrm{R}}$ and the classification of
our simulated merger remnants into slow ($\lambda_{\mathrm{R}}<0.1$) and fast
($\lambda_{\mathrm{R}}>0.1$) rotators. The remnants in the sequences S-MA-
SROT-BH and S-MI-BH are always slowly rotating. In the sequences S-MA-FROT-BH
and D-MA-BH, the remnants are initially fast-rotating, but become (or will
soon become) slow rotators in the minor merger sequences. The remnants in the
sequence D-MI-BH remain fast-rotating even after seven minor mergers.

We use the spin parameter $\lambda_{\mathrm{R}}$ as determined from the two-
dimensional kinematic data to characterize the rotation properties of our
simulated merger remnants. The common classification of ETGs into slow and
fast rotators is based on $\lambda_{\mathrm{R}}$ (Emsellem et al., 2007).
Slowly rotating ETGs have $\lambda_{\mathrm{R}}<0.1$ while for the fast
rotators $\lambda_{\mathrm{R}}>0.1$. Following Emsellem et al. (2007), the
parameter $\lambda_{\mathrm{R}}$ is defined from the two-dimensional kinematic
data as

$\lambda_{\mathrm{R}}=\frac{\Sigma_{\mathrm{i}}F_{\mathrm{i}}R_{\mathrm{i}}|V_{\mathrm{i}}|}{\Sigma_{\mathrm{i}}F_{\mathrm{i}}R_{\mathrm{i}}\sqrt{V_{\mathrm{i}}^{2}+\sigma_{\mathrm{i}}^{2}}}$
(15)

in which $F_{\mathrm{i}}$ is the flux (surface density) of each spaxel,
$R_{\mathrm{i}}$ their distance from the center of the galaxy and
$V_{\mathrm{i}}$ and $\sigma_{\mathrm{i}}$ the line-of-sight velocity and
velocity dispersion in each spaxel. The sum runs over all spaxels within the
effective radius $R_{\mathrm{e}}$ of the galaxy. For rotation-dominated
systems $\lambda_{\mathrm{R}}\sim 1$ and for non-rotating systems
$\lambda_{\mathrm{R}}\sim 0$. Our isolated initial galaxy models have rotation
parameters of $\lambda_{\mathrm{R}}\sim 0.0$ (host-S) and
$\lambda_{\mathrm{R}}\sim 0.7$ (host-D).

We show the evolution of the rotation parameter $\lambda_{\mathrm{R}}$ in the
simulated merger sequences in Fig. 18. The galaxy models in the merger
sequences S-MA-SROT-BH and S-MI-BH are always slowly rotating with
$\lambda_{\mathrm{R}}<0.04$ for each merger remnant. The disky-disky major
merger remnant of the sequence D-MA-BH is marginally fast rotating and
increases its rotation parameter to $\lambda_{\mathrm{R}}\sim 0.18$ in its
first minor merger. Afterwards, the minor mergers decrease the
$\lambda_{\mathrm{R}}$ bringing in stellar mass but no net angular momentum.
After $\sim 6$ minor mergers, the remnant D-MA-BH would be classified as a
slow rotator. The major merger of the sequence S-MA-FROT-BH results in a fast-
rotating remnant due to its merger orbit with $\lambda_{\mathrm{R}}\sim 0.42$.
Each merger generation gradually lowers the $\lambda_{\mathrm{R}}$ and after
seven minor mergers $\lambda_{\mathrm{R}}\sim 0.13$, a marginal fast rotator.
Extrapolating the evolution of $\lambda_{\mathrm{R}}$ in the minor mergers, we
argue that one more minor merger would be required to render the remnant of
the sequence S-MA-FROT-BH a slow rotator. Finally, the galaxy models of the
sequence D-MI-BH remain fast rotating with $\lambda_{\mathrm{R}}\sim 0.34$
after the minor merger generations.

Most importantly, we have shown that besides the $M_{\star}$–$R_{\mathrm{e}}$
size growth and the formation of flat central cores up to $R_{\mathrm{b}}\sim
1$ kpc size, the end products of merger sequences starting from compact early-
type galaxies can resemble $z=0$ massive ETGs also in their rotation
properties. Both initially rapidly rotating compact ETG models and spherical
models gaining rapid rotation due to merger orbit geometry can get rid of
their rotation in $6$–$9$ minor mergers. In our minor merger sequences this is
due to the fact that on average the minor mergers bring in more stellar mass
but no net angular momentum.

##  7 Summary and conclusions

We have studied the evolution of $z=2$ massive compact early-type galaxies
into present-day massive, slowly rotating, cored ETGs in simulated sequences
of major and minor galaxy mergers with and without supermassive black holes
using the KETJU code. Our main result is that $z=2$ compact ETGs (both
spherical and disky) can evolve through mergers into $z=0$ galaxies consistent
with massive ETGs in their effective radius, core size and rotation
properties. The presence of SMBHs in the galaxies is crucial for the models as
without SMBHs no low-density core is scoured. We further summarize our
findings in detail below.

Compact $z=2$ ETGs can reach the $z=0$ mass-size relation after $6$-$9$
generations minor mergers. Minor mergers are more efficient than major merger
in promoting the ETG size growth as suggested by virial arguments (Bezanson et
al., 2009; Naab et al., 2009). Without SMBHs, we find
$1.8\lesssim\alpha\lesssim 2.33$ for the mass growth relation
$R_{\mathrm{e}}\propto M_{\star}^{\alpha}$, consistent with the simulations of
Hilz et al. (2012); Hilz et al. (2013), with $\alpha\gtrsim 2.0$ required to
explain the observed ETG size growth (van Dokkum et al., 2010) since $z=2$. We
find consistently steeper galaxy size growth in simulations including SMBHs
with $2.31\lesssim\alpha_{\bullet}\lesssim 3.10$. Our results highlight the
paramount importance of multiple minor mergers including SMBHs on the
evolution of massive ETGs.

Repeated minor mergers with SMBHs can scour kpc-size cores into initially
compact galaxies, especially if the sequence begins with a major merger. The
break radii reached after seven generation of minor mergers are $0.62$ kpc
$\lesssim R_{\mathrm{b}}\lesssim 1.34$ kpc and are consistent with local, core
ETGs of similar SMBH mass (Rusli et al., 2013; Thomas et al., 2016). The model
S-MI-BH from a spherical progenitor without a major merger does not scour a
core even though its center is less dense than in the corresponding comparison
simulations without SMBHs. Hence major mergers are important to explain the
entire properties of massive ETGs (see also Kluge & Bender 2023). The
eccentric, unequal-mass binaries of the merger sequence S-MI-BH lose their
orbital energy mainly to GW emission and not to the stellar component, and
merge without significantly affecting their host nuclei. Without SMBHs, the
central stellar densities of the galaxies increase, and the core-like almost
flat exponential central parts of the initially disky models are turned into
mild power-law cusps. Thus, SMBHs are crucial for the $z=2$ compact ETGs to
evolve into the most massive ETGs in the local Universe: without SMBHs, there
is no flat central core, even when starting from a core-like exponential
progenitor galaxy.

We find 3D mass deficits up to $M_{\mathrm{def}}\sim
2.0N_{\mathrm{merg}}M_{\bullet}$ in our merger sequences, a value up to $\sim
4$ times higher than in the binary + bulge experiments of Merritt (2006). The
more observationally motivated 2D mass deficits of $2.2\lesssim
M_{\mathrm{def,2D}}/M_{\bullet}\lesssim 3.1$ are in good agreement with the
population of local cored ETGs (Hopkins & Hernquist, 2010; Rusli et al.,
2013). The scaling of the 2D mass deficit $M_{\mathrm{def,2D}}$ as a function
of SMBH mass $M_{\bullet}$ is steeper than the linear relation
$M_{\mathrm{def,2D}}\propto M_{\bullet}$, again consistent with the observed
mass and light deficits in the local core ETG population (Rusli et al., 2013;
Kormendy & Bender, 2009).

The effect of core scouring is not strictly restricted within the influence
radii of the central SMBHs as the SMBH binary can eject any star on an orbit
with low enough angular momentum. Thus, global properties of ETGs such as
their effective radii can be affected as well, connecting the small-scale
post-Newtonian SMBH binary dynamics to the large-scale galactic structure.
Extreme mass deficits of $>20\%$ of the total stellar mass caused by either
overmassive SMBHs or a large number of merger generations can increase the
effective radii $R_{\mathrm{e}}$ by up to a factor of $2$ due to scoured stars
being ejected or re-deposited in the outer parts of the galaxies with low
binding energies. We have presented a simple analytic scouring model which
qualitatively explains the additional galaxy size growth due to SMBH scouring
in the numerical simulations.

The shapes of our simulated merger remnants are mostly close to spherical or
mildly axisymmetric with only one model having an axis ratio less than
$b/a<0.9$. We attribute the shapes of the simulated merger remnants to the
fraction of tube orbits compared to radial box orbits within the half-mass
radius of the remnant. As a consequence of the dominance of tube orbits, the
central regions of the final simulated merger remnants are tangentially biased
($\beta_{\mathrm{rot}}<0$) within $r\lesssim 2$–$3R_{\mathrm{b}}$, consistent
with the observed orbital anisotropies of local core galaxies (Thomas et al.,
2014; Mehrgan et al., 2019; Neureiter et al., 2023; de Nicola et al., 2024).

Our disky galaxy models (D-MA-BH) and spherical setups (S-MA-FROT-BH) that
become fast-rotating due to a major merger orbit gradually evolve towards slow
rotation ($\lambda_{\mathrm{R}}<0.1$) in repeated minor mergers. This is due
to the fact that isotropic minor mergers bring in stellar mass and no net
angular momentum. The models would be classified as slow rotators after $6$
and $8$ minor mergers, respectively. Only the disky models D-MI-BH without a
major merger remain fast rotating after seven minor mergers.

Simulating the evolution of $z\gtrsim 2$ compact ETGs into present-day
massive, cored, slowly rotating ETGs self-consistently in cosmological zoom-in
simulations with KETJU remains an important future goal. In cosmological zoom-
in simulations, the hierarchical galaxy assembly provides the major and minor
mergers naturally, without needing to specify the numbers, times and mass
ratios of mergers by hand. The role of triple SMBH interactions occurring in
cosmological simulations (Mannerkoski et al., 2021) cannot be analyzed in
idealized binary merger simulations, such as this study. In addition, KETJU
has support for relativistic GW recoil kicks which we for simplicity did not
include in this study.

The first KETJU zoom-in simulations performed by Mannerkoski et al. (2021,
2022) did not specifically target compact ETGs and enabled KETJU integration
only at $z\lesssim 0.8$ and thus cannot address the early evolution
($0.8\lesssim z\lesssim 2.0$) of the nuclei of compact ETGs in mergers.
Nevertheless, the stellar mass resolution and the stellar gravitational
softening length of Mannerkoski et al. (2022) are within a factor of $\sim 3$
of the parameter values in this study. Thus, besides increased computational
costs due to initially very compact host galaxies of the targeted SMBHs there
is no fundamental obstacle for studying the evolution of red nuggets in self-
consistently cosmological zoom simulations with KETJU. Such simulations would
help to disentangle the relationships of high-redshift red nuggets, local
relic galaxies and intermediate-redshift relic candidates with varying degrees
of relicness, further scrutinizing the role of minor and major mergers in the
evolution of massive early-type galaxies and the corner-cases of galaxy
evolution.

## Data availability statement

The data relevant to this article will be shared on reasonable request to the
corresponding author.

## Acknowledgments

A. Rantala thanks Stavros Pastras and Natalia Lahén for their custom version
of the velocity map program. A. Rawlings acknowledges the support by the
University of Helsinki Research Foundation and the support of the Academy of
Finland grant 339127\. TN acknowledges support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311 from the DFG Cluster of Excellence
"ORIGINS". PHJ acknowledges the support by the European Research Council via
ERC Consolidator Grant KETJU (no. 818930) and the support of the Academy of
Finland grant 339127. The numerical simulations were performed using
facilities hosted by the Max Planck Computing and Data Facility (MPCDF) in
Garching, Germany.

## References

  * Appel (1985) Appel A. W., 1985, SIAM Journal on Scientific and Statistical Computing, 6, 85
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
  * Begelman et al. (1980) Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307
  * Belli et al. (2014) Belli S., Newman A. B., Ellis R. S., 2014, ApJ, 783, 117
  * Berczik et al. (2006) Berczik P., Merritt D., Spurzem R., Bischof H.-P., 2006, ApJ, 642, L21
  * Bezanson et al. (2009) Bezanson R., van Dokkum P. G., Tal T., Marchesini D., Kriek M., Franx M., Coppi P., 2009, ApJ, 697, 1290
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press, Princeton, NJ USA
  * Blanchet (2014) Blanchet L., 2014, Living Reviews in Relativity, 17, 2
  * Boylan-Kolchin et al. (2005) Boylan-Kolchin M., Ma C.-P., Quataert E., 2005, MNRAS, 362, 184
  * Buitrago et al. (2008) Buitrago F., Trujillo I., Conselice C. J., Bouwens R. J., Dickinson M., Yan H., 2008, ApJ, 687, L61
  * Buitrago et al. (2013) Buitrago F., Trujillo I., Conselice C. J., Häußler B., 2013, MNRAS, 428, 1460
  * Bulirsch & Stoer (1966) Bulirsch R., Stoer J., 1966, Numerische Mathematik, 8, 1
  * Cannarozzo et al. (2023) Cannarozzo C., et al., 2023, MNRAS, 520, 5651
  * Caon et al. (1993) Caon N., Capaccioli M., D’Onofrio M., 1993, MNRAS, 265, 1013
  * Cappellari & Copin (2003) Cappellari M., Copin Y., 2003, MNRAS, 342, 345
  * Cappellari et al. (2013) Cappellari M., et al., 2013, MNRAS, 432, 1709
  * Carrasco et al. (2010) Carrasco E. R., Conselice C. J., Trujillo I., 2010, MNRAS, 405, 2253
  * Cassata et al. (2011) Cassata P., et al., 2011, ApJ, 743, 96
  * Cenarro & Trujillo (2009) Cenarro A. J., Trujillo I., 2009, ApJ, 696, L43
  * Chandrasekhar (1943) Chandrasekhar S., 1943, ApJ, 97, 255
  * Ciotti et al. (2007) Ciotti L., Lanzoni B., Volonteri M., 2007, ApJ, 658, 65
  * Comerón et al. (2023) Comerón S., et al., 2023, A&A, 675, A143
  * Courteau & Dutton (2015) Courteau S., Dutton A. A., 2015, ApJ, 801, L20
  * Daddi et al. (2005) Daddi E., et al., 2005, ApJ, 626, 680
  * Damjanov et al. (2014) Damjanov I., Hwang H. S., Geller M. J., Chilingarian I., 2014, ApJ, 793, 39
  * Dehnen (1993) Dehnen W., 1993, MNRAS, 265, 250
  * Dehnen (2000) Dehnen W., 2000, ApJ, 536, L39
  * Dehnen (2002) Dehnen W., 2002, Journal of Computational Physics, 179, 27
  * Dehnen (2014) Dehnen W., 2014, Computational Astrophysics and Cosmology, 1, 1
  * Deuflhard (1983) Deuflhard P., 1983, Numerische Mathematik, 41, 399–422
  * Dullo (2019) Dullo B. T., 2019, ApJ, 886, 80
  * Ebisuzaki et al. (1991) Ebisuzaki T., Makino J., Okumura S. K., 1991, Nature, 354, 212
  * Emsellem et al. (2007) Emsellem E., et al., 2007, MNRAS, 379, 401
  * Ene et al. (2020) Ene I., Ma C.-P., Walsh J. L., Greene J. E., Thomas J., Blakeslee J. P., 2020, ApJ, 891, 65
  * Faber et al. (1997) Faber S. M., et al., 1997, AJ, 114, 1771
  * Feldmann et al. (2010) Feldmann R., Carollo C. M., Mayer L., Renzini A., Lake G., Quinn T., Stinson G. S., Yepes G., 2010, ApJ, 709, 218
  * Ferrarese & Ford (2005) Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523
  * Ferrarese & Merritt (2000) Ferrarese L., Merritt D., 2000, ApJ, 539, L9
  * Ferrarese et al. (1994) Ferrarese L., van den Bosch F. C., Ford H. C., Jaffe W., O’Connell R. W., 1994, AJ, 108, 1598
  * Ferré-Mateu et al. (2012) Ferré-Mateu A., Vazdekis A., Trujillo I., Sánchez-Blázquez P., Ricciardelli E., de la Rosa I. G., 2012, MNRAS, 423, 632
  * Ferré-Mateu et al. (2017) Ferré-Mateu A., Trujillo I., Martín-Navarro I., Vazdekis A., Mezcua M., Balcells M., Domínguez L., 2017, MNRAS, 467, 1929
  * Franx et al. (1991) Franx M., Illingworth G., de Zeeuw T., 1991, ApJ, 383, 112
  * Franx et al. (2008) Franx M., van Dokkum P. G., Förster Schreiber N. M., Wuyts S., Labbé I., Toft S., 2008, ApJ, 688, 770
  * Frigo et al. (2021) Frigo M., Naab T., Rantala A., Johansson P. H., Neureiter B., Thomas J., Rizzuto F., 2021, MNRAS, 508, 4610
  * Gebhardt et al. (2000) Gebhardt K., et al., 2000, ApJ, 539, L13
  * Gerhard (1993) Gerhard O. E., 1993, MNRAS, 265, 213
  * Gragg (1965) Gragg W. B., 1965, SIAM Journal on Numerical Analysis, 2, 384
  * Graham (2004) Graham A. W., 2004, ApJ, 613, L33
  * Graham et al. (2003) Graham A. W., Erwin P., Trujillo I., Asensio Ramos A., 2003, AJ, 125, 2951
  * Graham et al. (2016) Graham A. W., Durré M., Savorgnan G. A. D., Medling A. M., Batcheldor D., Scott N., Watson B., Marconi A., 2016, ApJ, 819, 43
  * Greengard & Rokhlin (1987) Greengard L., Rokhlin V., 1987, Journal of Computational Physics, 73, 325
  * Gualandris et al. (2017) Gualandris A., Read J. I., Dehnen W., Bortolas E., 2017, MNRAS, 464, 2301
  * Hairer et al. (1993) Hairer E., Nørsett S. P., Wanner G., 1993, Solving Ordinary Differential Equations I: Nonstiff Problems. Springer Berlin, Heidelberg
  * Häring & Rix (2004) Häring N., Rix H.-W., 2004, ApJ, 604, L89
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hilz et al. (2012) Hilz M., Naab T., Ostriker J. P., Thomas J., Burkert A., Jesseit R., 2012, MNRAS, 425, 3119
  * Hilz et al. (2013) Hilz M., Naab T., Ostriker J. P., 2013, MNRAS, 429, 2924
  * Holley-Bockelmann & Khan (2015) Holley-Bockelmann K., Khan F. M., 2015, ApJ, 810, 139
  * Hopkins & Hernquist (2010) Hopkins P. F., Hernquist L., 2010, MNRAS, 407, 447
  * Jaffe (1983) Jaffe W., 1983, MNRAS, 202, 995
  * Jesseit et al. (2005) Jesseit R., Naab T., Burkert A., 2005, MNRAS, 360, 1185
  * Johansson et al. (2009) Johansson P. H., Naab T., Burkert A., 2009, ApJ, 690, 802
  * Johansson et al. (2012) Johansson P. H., Naab T., Ostriker J. P., 2012, ApJ, 754, 115
  * Khan et al. (2011) Khan F. M., Just A., Merritt D., 2011, ApJ, 732, 89
  * Khan et al. (2013) Khan F. M., Holley-Bockelmann K., Berczik P., Just A., 2013, ApJ, 773, 100
  * Khan et al. (2016) Khan F. M., Fiacconi D., Mayer L., Berczik P., Just A., 2016, ApJ, 828, 73
  * Khochfar & Burkert (2006) Khochfar S., Burkert A., 2006, A&A, 445, 403
  * Kluge & Bender (2023) Kluge M., Bender R., 2023, ApJS, 267, 41
  * Kormendy & Bender (1996) Kormendy J., Bender R., 1996, ApJ, 464, L119
  * Kormendy & Bender (2009) Kormendy J., Bender R., 2009, ApJ, 691, L142
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Kormendy & Richstone (1995) Kormendy J., Richstone D., 1995, ARA&A, 33, 581
  * Kormendy et al. (2009) Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 216
  * Krajnović et al. (2011) Krajnović D., et al., 2011, MNRAS, 414, 2923
  * Lauer et al. (2007) Lauer T. R., et al., 2007, ApJ, 662, 808
  * Lynden-Bell (1967) Lynden-Bell D., 1967, MNRAS, 136, 101
  * Mannerkoski et al. (2021) Mannerkoski M., Johansson P. H., Rantala A., Naab T., Liao S., 2021, ApJ, 912, L20
  * Mannerkoski et al. (2022) Mannerkoski M., Johansson P. H., Rantala A., Naab T., Liao S., Rawlings A., 2022, ApJ, 929, 167
  * Mannerkoski et al. (2023) Mannerkoski M., Rawlings A., Johansson P. H., Naab T., Rantala A., Springel V., Irodotou D., Liao S., 2023, MNRAS, 524, 4062
  * McMillan & Aarseth (1993) McMillan S. L. W., Aarseth S. J., 1993, ApJ, 414, 200
  * Mehrgan et al. (2019) Mehrgan K., Thomas J., Saglia R., Mazzalay X., Erwin P., Bender R., Kluge M., Fabricius M., 2019, ApJ, 887, 195
  * Merritt (2006) Merritt D., 2006, ApJ, 648, 976
  * Merritt (2013) Merritt D., 2013, Dynamics and Evolution of Galactic Nuclei. Princeton University Press, Princeton
  * Merritt & Milosavljević (2005) Merritt D., Milosavljević M., 2005, Living Reviews in Relativity, 8, 8
  * Mikkola & Aarseth (2002) Mikkola S., Aarseth S., 2002, Celestial Mechanics and Dynamical Astronomy, 84, 343
  * Mikkola & Merritt (2006) Mikkola S., Merritt D., 2006, MNRAS, 372, 219
  * Mikkola & Merritt (2008) Mikkola S., Merritt D., 2008, AJ, 135, 2398
  * Mikkola & Tanikawa (1999) Mikkola S., Tanikawa K., 1999, MNRAS, 310, 745
  * Milosavljević & Merritt (2001) Milosavljević M., Merritt D., 2001, ApJ, 563, 34
  * Milosavljević & Merritt (2003) Milosavljević M., Merritt D., 2003, ApJ, 596, 860
  * Mirza et al. (2017) Mirza M. A., Tahir A., Khan F. M., Holley-Bockelmann H., Baig A. M., Berczik P., Chishtie F., 2017, MNRAS, 470, 940
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Naab et al. (2006) Naab T., Khochfar S., Burkert A., 2006, ApJ, 636, L81
  * Naab et al. (2009) Naab T., Johansson P. H., Ostriker J. P., 2009, ApJ, 699, L178
  * Naab et al. (2014) Naab T., et al., 2014, MNRAS, 444, 3357
  * Nasim et al. (2021) Nasim I. T., Gualandris A., Read J. I., Antonini F., Dehnen W., Delorme M., 2021, MNRAS, 502, 4794
  * Neureiter et al. (2023) Neureiter B., Thomas J., Rantala A., Naab T., Mehrgan K., Saglia R., de Nicola S., Bender R., 2023, ApJ, 950, 15
  * Nipoti et al. (2009) Nipoti C., Treu T., Bolton A. S., 2009, ApJ, 703, 1531
  * O’Leary et al. (2021) O’Leary J. A., Moster B. P., Naab T., Somerville R. S., 2021, MNRAS, 501, 3215
  * Oser et al. (2010) Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312
  * Oser et al. (2012) Oser L., Naab T., Ostriker J. P., Johansson P. H., 2012, ApJ, 744, 63
  * Partmann et al. (2023) Partmann C., Naab T., Rantala A., Genina A., Mannerkoski M., Johansson P. H., 2023, arXiv e-prints, p. arXiv:2310.08079
  * Peters (1964) Peters P. C., 1964, Phys. Rev., 136, B1224
  * Peters & Mathews (1963) Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  * Poggianti et al. (2013) Poggianti B. M., et al., 2013, ApJ, 762, 77
  * Poisson & Will (2014) Poisson E., Will C. M., 2014, Gravity. Cambridge University Press
  * Preto & Tremaine (1999) Preto M., Tremaine S., 1999, AJ, 118, 2532
  * Preto et al. (2011) Preto M., Berentzen I., Berczik P., Spurzem R., 2011, ApJ, 732, L26
  * Qu et al. (2017) Qu Y., et al., 2017, MNRAS, 464, 1659
  * Quenneville et al. (2024) Quenneville M. E., Blakeslee J. P., Ma C.-P., Greene J. E., Gwyn S. D. J., Ciccone S., Nyiri B., 2024, MNRAS, 527, 249
  * Quilis & Trujillo (2013) Quilis V., Trujillo I., 2013, ApJ, 773, L8
  * Quinlan (1996) Quinlan G. D., 1996, New Astron., 1, 35
  * Quinlan & Hernquist (1997) Quinlan G. D., Hernquist L., 1997, New Astron., 2, 533
  * Rantala et al. (2017) Rantala A., Pihajoki P., Johansson P. H., Naab T., Lahén N., Sawala T., 2017, ApJ, 840, 53
  * Rantala et al. (2018) Rantala A., Johansson P. H., Naab T., Thomas J., Frigo M., 2018, ApJ, 864, 113
  * Rantala et al. (2019) Rantala A., Johansson P. H., Naab T., Thomas J., Frigo M., 2019, ApJ, 872, L17
  * Rantala et al. (2020) Rantala A., Pihajoki P., Mannerkoski M., Johansson P. H., Naab T., 2020, MNRAS, 492, 4131
  * Rawlings et al. (2023) Rawlings A., Mannerkoski M., Johansson P. H., Naab T., 2023, MNRAS, 526, 2688
  * Rodriguez-Gomez et al. (2016) Rodriguez-Gomez V., et al., 2016, MNRAS, 458, 2371
  * Röttgers et al. (2020) Röttgers B., Naab T., Cernetic M., Davé R., Kauffmann G., Borthakur S., Foidl H., 2020, MNRAS, 496, 152
  * Rusli et al. (2013) Rusli S. P., Erwin P., Saglia R. P., Thomas J., Fabricius M., Bender R., Nowak N., 2013, AJ, 146, 160
  * Saglia et al. (2016) Saglia R. P., et al., 2016, ApJ, 818, 47
  * Saulder et al. (2015) Saulder C., van den Bosch R. C. E., Mieske S., 2015, A&A, 578, A134
  * Sérsic (1963) Sérsic J. L., 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Shen et al. (2003) Shen S., Mo H. J., White S. D. M., Blanton M. R., Kauffmann G., Voges W., Brinkmann J., Csabai I., 2003, MNRAS, 343, 978
  * Spiniello et al. (2021) Spiniello C., et al., 2021, A&A, 654, A136
  * Spiniello et al. (2024) Spiniello C., et al., 2024, MNRAS, 527, 8793
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * Szomoru et al. (2011) Szomoru D., Franx M., Bouwens R. J., van Dokkum P. G., Labbé I., Illingworth G. D., Trenti M., 2011, ApJ, 735, L22
  * Szomoru et al. (2012) Szomoru D., Franx M., van Dokkum P. G., 2012, ApJ, 749, 121
  * Szomoru et al. (2013) Szomoru D., Franx M., van Dokkum P. G., Trenti M., Illingworth G. D., Labbé I., Oesch P., 2013, ApJ, 763, 73
  * Taylor et al. (2010) Taylor E. N., Franx M., Glazebrook K., Brinchmann J., van der Wel A., van Dokkum P. G., 2010, ApJ, 720, 723
  * Thomas et al. (2014) Thomas J., Saglia R. P., Bender R., Erwin P., Fabricius M., 2014, ApJ, 782, 39
  * Thomas et al. (2016) Thomas J., Ma C.-P., McConnell N. J., Greene J. E., Blakeslee J. P., Janish R., 2016, Nature, 532, 340