The Formation of Extremely Diffuse Galaxy Cores by Merging Supermassive Black Holes

Our simulations are run using the recently developed KETJU code, which is built on the widely used galaxy simulation code GADGET-3 (Springel 2005). Here we briefly review the main features of the code and refer the reader to R17 for a complete description. The central idea of the KETJU code is to include special regions around every black hole particle, in which the dynamics of SMBHs and stellar particles is modeled using a nonsoftened algorithmic chain regularization technique (Mikkola & Merritt 2006, 2008) that also includes post-Newtonian (PN) corrections (e.g., Will 2006). The remaining particles far from the SMBHs are evolved using the standard GADGET-3 leapfrog integrator, with the softened gravitational force calculated using either a pure tree or a hybrid tree–mesh TreePM algorithm (Barnes & Hut 1986; Springel et al. 2001).

The KETJU code operates by dividing the simulation particles into three categories (see Figure 2 in R17 for an illustration). All the SMBH particles and stellar particles, which lie within a user-defined chain radius (r_chain), are marked as chain subsystem particles. Particles that lie just outside the chain radius but induce a strong tidal perturbation on the chain system are marked as perturber particles. Finally, the remaining particles that are far from any SMBHs are treated as ordinary GADGET-3 particles with respect to the force calculation. The chain membership of the simulation particles is evaluated at the beginning of each timestep. KETJU also allows for both multiple simultaneous chain subsystems and several SMBHs in a single subsystem.

The chain radius should be chosen carefully to ensure that the regularized regions are large enough to accurately simulate the dynamics around SMBHs but at the same time small enough in order to not excessively degrade the performance of the code. The minimum chain radius is set by the Plummer-equivalent gravitational softening length in GADGET-3, r_chain > 2.8 × , which ensures that the star–SMBH and SMBH–SMBH interactions always remain nonsoftened in KETJU. The maximum chain radius is limited by the number of particles in the chain subsystem. The scaling of the AR-CHAIN integrator with particle number is of the order of ${ \mathcal O }({N}^{2.13\mbox{--}2.20})$ (R17). In the current implementation, the chain region is limited to roughly N_c ∼ 500 particles for stellar density profiles ρ_⋆ ∝ r⁻¹, beyond which the simulations becomes unfeasible to run. Thus, for steeper density profiles ρ_⋆ ∝ r^−γ, with γ > 1, and higher numerical resolution, the chain radius must be reduced in size. For the typical applications presented in this paper, the chain radius, r_chain, is set at a few tens of pc, resulting in a few hundred particles within a single chain region.

The regularized subsystems are also tidally perturbed by their surroundings, imposing a perturbing acceleration ${{\boldsymbol{f}}}_{i}$ on every particle i in the subsystem. Particles in the immediate vicinity of a chain subsystem are predicted forward in time during the AR-CHAIN external perturbations calculation, while the tidal field from the more distant particles remains constant during the AR-CHAIN integration. The distinction between the nearby perturbing particles and the particles in the far field is set by the perturber radius,

$$\begin{eqnarray}&&{r}_{\mathrm{pert}}={\gamma }_{\mathrm{pert}}\times {\left(\displaystyle \frac{m}{{M}_{\bullet }}\right)}^{1/3}\times {r}_{\mathrm{chain}},\end{eqnarray} \tag{ 1 }$$

in which m is the mass of the perturbing simulation particle and M_• is the mass of the SMBH in the chain subsystem. The constant γ_pert is typically chosen so that r_pert = 2 × r_chain.

During each global GADGET-3 timestep, the particles in the chain subsystems are propagated using a novel reimplementation of the AR-CHAIN algorithm (Mikkola & Merritt 2008) by R17. The AR-CHAIN algorithm has the following three main aspects. First, the equations of motion of the simulation particles are time-transformed using a new time coordinate in the regularized regions (e.g., Mikkola & Aarseth 2002), with the integration proceeding using a leapfrog algorithm. In essence, algorithmic regularization works by transforming the equations of motion by introducing a fictitious time variable such that integration by the common leapfrog method yields exact orbits for a Newtonian two-body problem including two-body collisions. Second, numerical round-off errors are greatly reduced by the introduction of a coordinate system based on chained interparticle vectors (e.g., Mikkola & Aarseth 1993). Finally, the Gragg–Bulirsch–Stoer (GBS; Gragg 1965; Bulirsch & Stoer 1966) extrapolation method yields high numerical accuracy in orbit integrations at a preset user-given error tolerance level (η_GBS).

The regularized integrator has multiple advantages compared to the standard GADGET-3 leapfrog algorithm when integrating the equations of motion of the simulation particles near SMBHs. In addition to high numerical accuracy and sophisticated error control, gravitational softening is not applied in the regularized regions. Thus, the N-body dynamics is accurately resolved even at very small particle separations. Furthermore, the extension of phase space with an auxiliary velocity variable in the algorithm (Hellström & Mikkola 2010; Pihajoki 2015) allows for the efficient implementation of the velocity-dependent PN corrections in the motion of the simulation particles (e.g., Will 2006). In this study, we include PN terms for the SMBH–SMBH interaction only, since the star–SMBH terms would be unphysically large, given that the stellar particles have very large masses, m_⋆ ≫ M_⊙. In the included PN terms, we go up to order PN2.5, the highest of which includes the radiation reaction term due to losses caused by the emission of GWs. Higher-order PN corrections, spin terms, and their cross terms are not used in this study.

To summarize, KETJU is able to follow SMBH dynamics accurately from global galactic scales to particle separations of a few tens of Schwarzschild radii of the SMBHs where the PN approach breaks down. Before this happens, the SMBH particles are merged using a merger criterion that is based on the analytic Peters & Mathews (1963) estimate for the GW-driven decay of the orbital binary semimajor axis (see Equation (19)). For additional details concerning the AR-CHAIN integrator, see Appendices A and B in R17.

Here we present code updates compared to the previous code version presented in R17. The far-field acceleration part of the tidal acceleration ${{\boldsymbol{f}}}_{i}$ is now calculated individually for every particle in the regularized subsystem. In the previous version, the far-field acceleration was imposed on the center of mass of the chain subsystem only.

In addition, the regularized subsystems are now temporarily deconstructed for the GADGET-3 tree gravity computation. The gravitational forces between particles in the same subsystem are ignored in the tree calculation, and the center of mass of the regularized subsystem is propagated as an ordinary GADGET-3 SMBH particle in the simulation. This procedure eliminates the demand for the distinct force correction used in Karl et al. (2015) and the previous version of the KETJU code. The old force correction procedure included a step in which the tree force on the center of mass of a subsystem was canceled using a directly summed term. This step produced low-level spurious noise to the total gravitational force as the directly summed component, and the tree component did not always cancel exactly. Overall, the code updates have only a minor effect on the simulation results but are added for the sake of code clarity and to ensure a consistent formulation of Newtonian gravity that also results in a slight performance upgrade.

Most importantly, the interface between GADGET-3 and the AR-CHAIN integrator is rewritten for an increased performance for simulations with high particle numbers, in excess of N ≳ 5 × 10⁶. In addition, the regularized subsystems are now built and deconstructed every smallest GADGET-3 timestep. Now the tree rebuild on every smallest GADGET-3 timestep is no longer necessary. Nevertheless, the tree domain update frequency is kept smaller than the typical value ∼0.05 in GADGET-3 (Springel 2005). In this study, we use the value of 0.01. This code interface formulation anticipates the future use of KETJU for hydrodynamical simulations that also include GADGET-3 subresolution feedback modeling (Hu et al. 2014, 2016). However, in the present study, we still restrict ourselves to simulating collisionless gas-poor (dry) galaxy mergers that include central SMBHs.

In setting up the galaxies, we use the Dehnen density-potential (Dehnen 1993) models, which are commonly used for modeling the stellar components of early-type galaxies and bulges. The spherically symmetric one-component Dehnen model is defined by three parameters: the total mass M, scale radius a, and central slope of the density profile γ. The allowed values for the central slopes γ span 0 ≤ γ < 3. The most commonly used Dehnen models are the Hernquist profile with γ = 1 (Hernquist 1990), the Jaffe profile with γ = 2 (Jaffe 1983), and the γ = 3/2 profile, which closely resembles the de Vaucouleurs profile (de Vaucouleurs 1948; Dehnen 1993) when projected.

The Dehnen density-potential pair family is defined as

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{(3-\gamma )M}{4\pi }\displaystyle \frac{a}{{r}^{\gamma }{(r+a)}^{4-\gamma }},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\phi (r)=\displaystyle \frac{{GM}}{a}\times \left\{\begin{array}{ll}-\tfrac{1}{2-\gamma }\left[1-{\left(\tfrac{r}{r+a}\right)}^{2-\gamma }\right], & \gamma \ne 2\\ \mathrm{ln}\tfrac{r}{r+a}, & \gamma =2.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

From the density profile, we can solve the cumulative mass profile M(r) and the three-dimensional half-mass–radius r_1/2,

$$\begin{eqnarray}&&M(r)=4\pi {\int }_{0}^{r}\rho (r){r}^{2}{dr}=M{\left(\displaystyle \frac{r}{r+a}\right)}^{3-\gamma },\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{r}_{1/2}=a{({2}^{1/(3-\gamma )}-1)}^{-1},\end{eqnarray} \tag{ 5 }$$

whereas the projected half-mass (effective) radius R_e can be well approximated by R_e ≈ 3/4 r_1/2.

We construct spherically symmetric, isotropic multicomponent ICs consisting of a stellar bulge, DM halo, and central SMBH. The ICs are generated using the distribution function method (e.g., Merritt 1985; Ciotti & Pellegrini 1992) following the approach of Hilz et al. (2012). The distribution functions f_i for the different components of the galaxy are obtained from the corresponding density profile ρ_i and the total gravitational potential Φ_T using Eddington's formula (Binney & Tremaine 2008),

$$\begin{eqnarray}&&{f}_{i}({ \mathcal E })=\displaystyle \frac{1}{\sqrt{8}{\pi }^{2}}{\int }_{{{\rm{\Phi }}}_{{\rm{T}}}=0}^{{{\rm{\Phi }}}_{{\rm{T}}}={ \mathcal E }}\displaystyle \frac{{d}^{2}{\rho }_{i}}{d{{\rm{\Phi }}}_{{\rm{T}}}^{2}}\displaystyle \frac{d{{\rm{\Phi }}}_{{\rm{T}}}}{\sqrt{{ \mathcal E }-{{\rm{\Phi }}}_{{\rm{T}}}}},\end{eqnarray} \tag{ 6 }$$

in which ${ \mathcal E }=-\tfrac{1}{2}{v}^{2}-{{\rm{\Phi }}}_{{\rm{T}}}+{{\rm{\Phi }}}_{0}$ is the (positive) energy relative to the chosen zero point of the potential Φ₀. In general, the zero point is chosen so that f_i > 0 for ${ \mathcal E }\gt 0$ and f_i = 0 for ${ \mathcal E }\leqslant 0$. For an isolated system extending to infinity, such as the Dehnen profiles in our case, we can set Φ₀ = 0. The particle positions of the components are drawn from the cumulative mass profile (Equation (4)), after which the velocities of the particles are sampled from precomputed and tabulated distribution functions (R17).

The SMBH is represented by a single point mass M_• at rest at the origin. The stellar component is parameterized by the total stellar mass M_⋆, the effective radius of the stellar component R_e, and the central density profile slope γ. The DM halo is modeled as a Hernquist sphere (γ = 1) with a total mass M_DM (Springel et al. 2005). The scale radius of the DM component is computed assuming a DM fraction f_DM inside the three-dimensional stellar half-mass–radius r_1/2. Defining the DM fraction inside a radius r as

$$\begin{eqnarray}&&{f}_{\mathrm{DM}}(r)=\displaystyle \frac{{M}_{\mathrm{DM}}(r)}{{M}_{\star }(r)+{M}_{\mathrm{DM}}(r)},\end{eqnarray} \tag{ 7 }$$

the scale radius of the DM component can be derived by inserting Equation (4) describing the cumulative mass profile into Equation (7) and using the definition of the half-mass–radius from Equation (5), resulting in

$$\begin{eqnarray}&&{a}_{\mathrm{DM}}\approx \displaystyle \frac{4}{3}\left[\sqrt{\displaystyle \frac{2{M}_{\mathrm{DM}}}{{M}_{\star }}\left(\displaystyle \frac{1}{{f}_{\mathrm{DM}}({r}_{1/2})}-1\right)}-1\right]{R}_{{\rm{e}}}.\end{eqnarray} \tag{ 8 }$$

We study a simplified scenario in which the current structural properties of massive core elliptical galaxies were shaped in a final dry major merger of two early-type galaxies. The progenitor galaxies are identical to each other in each simulation run. The structural properties of the progenitors are motivated by the observations and dynamical modeling of the core elliptical galaxy NGC 1600 (T16). In Table 1, we list the physical parameters of the progenitor models (M_⋆, R_e, M_DM, f_DM), which remain the same for every IC in this study. We assume that all of the stellar mass was already present in the progenitor galaxies at the time of this final major merger, yielding progenitor stellar and DM halo masses of M_⋆ = 4.15 × 10¹¹ M_⊙ and M_DM = 7.5 × 10¹³ M_⊙ (T16).

Using virial arguments, Naab et al. (2009) showed that the effective radius of a galaxy grows in a dry equal-mass major by roughly a factor of ∼2. We set the effective radius of our progenitor galaxies to R_e = 7 kpc, as observations have shown that the effective radius of NGC 1600 is R_e ∼ 14–16 kpc (e.g., Matthias & Gerhard 1999; Trager et al. 2000; Fukazawa et al. 2006). The DM fraction was set to f_DM = 0.25, which is a rather conservative estimate, as the observed central DM fractions inside R_e of local early-type galaxies are typically found to be relatively low (e.g., Cappellari et al. 2013; Courteau & Dutton 2015).

We construct a sample of 14 ICs for major-merger simulations. In half of the sample, the stellar component of the progenitor galaxies is modeled using a Hernquist profile (γ = 1), and in the other half of ICs, we use a steeper γ = 3/2 profile. All progenitor models are listed in Table 2.

For 12 ICs, we include an SMBH with initial masses ranging from M_• = 8.5 × 10⁸ to 8.5 × 10⁹ M_⊙. For comparison, we also include one set of models without a central SMBH. The largest SMBH masses are consistent with the final SMBH mass obtained using the dynamical modeling of NGC 1600 (M_• = 1.7 × 10¹⁰ M_⊙), while the ICs with lower-mass SMBHs lie closer to the observed M_•–σ relation (T16). Note that even if the stellar and DM density profiles are identical in the ICs with different SMBH masses, the circular velocity and velocity dispersion profiles differ due to the different gravitational potential of the SMBHs, especially near the central regions of the galaxies.

Similarly to R17, the progenitor galaxies are set on a nearly parabolic merger orbit with an initial separation of d = 30 kpc. From this orbit, the approach of the galaxies is swift, and the central stellar cusps merge before t ∼ 300 Myr.

The values of the GADGET-3 and AR-CHAIN numerical accuracy parameters are based on the simulations of R17. We set the GADGET-3 integrator error tolerance to η = 0.002 and the force accuracy to α = 0.005 using the standard cell opening criterion (Springel 2005). For the AR-CHAIN integrator, we choose a GBS tolerance of η_GBS = 10⁻⁶.

The chain radius is set to r_chain = 30 pc in the runs with an initial slope for the stellar density profile of γ = 1 and r_chain = 10 pc for the runs with a steeper γ = 3/2 slope. These choices for the chain radii result in a few hundred particles in the regularized regions at the start of the simulations, making the runs numerically feasible, as discussed in Section 2.1. The gravitational softening lengths for the GADGET-3 stellar particles are selected accordingly: _⋆ = 10 pc in the γ = 1 runs and _⋆ = 3.5 pc in the γ = 3/2 simulations, chosen to fulfill the criterion r_chain > 2.8 × . The softening length for the DM component is _DM = 100 pc in all simulation runs. The tidal perturber parameter γ_pert is set following the condition described in Section 2.1, resulting in r_pert = 2 × r_chain for each simulation run.

The individual progenitor galaxies are sampled using N_⋆ = 4.15 × 10⁶ stellar particles and N_DM = 1.0 × 10⁷ DM particles for each galaxy, yielding particle masses of m_⋆ = 1.0 × 10⁵ M_⊙ and m_DM = 7.5 × 10⁶ M_⊙, respectively. Consequently, the ratio between the SMBH mass and the stellar particle mass in the simulation sample is in the range 8500 ≤ M_•/m_⋆ ≤ 85,000, which is sufficient to realistically study the evolution of SMBH binaries in the field of lighter particles (Mikkola & Valtonen 1992). With a sufficiently high SMBH–to–stellar particle mass ratio, the stochastic effects for the SMBH binary evolution are minimized.

When the central regions of the merging galaxies coalesce, the SMBHs are rapidly driven by the force of dynamical friction to within the influence radius r_h of each other. We define r_h as

$$\begin{eqnarray}&&{r}_{{\rm{h}}}=\displaystyle \frac{{{GM}}_{\bullet }}{{\sigma }_{\star }^{2}},\end{eqnarray} \tag{ 9 }$$

where σ_⋆ is the LOS stellar velocity dispersion of the host galaxy, measured inside the effective radius R_e. Almost simultaneously when arriving within r_h, the SMBHs also form a gravitationally bound binary. In this section, we focus on the simulation sample with γ = 3/2, and for these simulations, the binaries become bound between t = 248 and 275 Myr, with the earlier times referring to the initially more massive SMBHs.

The subsequent evolution of the binary toward smaller semimajor axis (a) values is very swift. Merritt (2013) provided an analytical estimate for the timescale on which the binary transfers energy to the surrounding stellar population, and a similar computation of the energy transfer rate can also be found in Milosavljević & Merritt (2001). The orbital energy of the binary, defined as

$$\begin{eqnarray}&&E=-\displaystyle \frac{{{GM}}_{1}{M}_{2}}{2a},\end{eqnarray} \tag{ 10 }$$

is transferred into the stellar component on the timescale

$$\begin{eqnarray}&&{T}_{{\rm{E}}}={\left|\displaystyle \frac{1}{E}\displaystyle \frac{{dE}}{{dt}}\right|}^{-1}\approx \displaystyle \frac{{\sigma }_{\star }^{3}}{{{CG}}^{2}{M}_{2}{\rho }_{\star }}.\end{eqnarray} \tag{ 11 }$$

Here σ_⋆ is now the 3D stellar velocity dispersion, ρ_⋆ is the stellar density, and M₁ and M₂ are the masses of the SMBHs, with M₂ being the less massive SMBH. The constant C depends on the primary energy transfer mechanism, which is either dynamical friction or three-body slingshots. However, the value of C is roughly of the order of C ∼ 10 for both physical processes. Note that when r_h > a but the binary is still very wide, the energy transfer mechanism is a mixture of the two processes. Computed within the central r < 0.19 kpc of the SMBHs (see T16), the energy transfer timescales are T_E ∼ 4.3–6.9 Myr for the binaries when a = r_h. This timescale is small compared to the crossing time of the merger remnant, which is t_cross = R_e/σ_⋆ ∼ 35 Myr.

The shrinking SMBH binaries rapidly become hard. We adopt the definition of Merritt (2013) for the hardness of the binary,

$$\begin{eqnarray}&&{a}_{{\rm{h}}}=\displaystyle \frac{G\mu }{4{\sigma }_{\star }^{2}}=\displaystyle \frac{q}{1+q}\displaystyle \frac{{r}_{{\rm{h}}}}{4},\end{eqnarray} \tag{ 12 }$$

in which μ = M₁M₂/(M₁ + M₂) is the reduced mass and q = M₂/M₁ < 1 is the mass ratio of the binary. For an equal-mass SMBH binary (M₁ = M₂, i.e., q = 1), the hard semimajor axis becomes

$$\begin{eqnarray}&&{a}_{{\rm{h}}}=\displaystyle \frac{{r}_{{\rm{h}}}}{8}.\end{eqnarray} \tag{ 13 }$$

For our simulation sample, the influence radii of the SMBHs are in the range 0.09 kpc < r_h < 0.87 kpc, and the hard semimajor axes are in the range 0.01 kpc < a_h < 0.11 kpc. Note that these separations are quite large when compared to values commonly found in the literature, the reason being that most of the SMBH masses in our study are above the values expected from the scaling relations between M_• and the host galaxy properties. The SMBH binaries shrink from a = r_h to a = a_h in 13–20 Myr, which is within a factor of 2–3 of both the binary energy transfer timescale T_E ∼ 4.3–6.9 Myr and the crossing time of the host galaxy, t_cross ∼ 35 Myr. However, using another popular definition,

$$\begin{eqnarray}&&{a}_{{\rm{h}}}=\displaystyle \frac{1}{{(1+q)}^{2}}\displaystyle \frac{{r}_{{\rm{h}}}}{4},\end{eqnarray} \tag{ 14 }$$

for the hard separation (e.g., Merritt 2006), the definition differing from Equation (12) by a factor of 2 for an equal-mass binary, we obtain shrinking timescales of 44–55 Myr. These timescales are somewhat longer than the crossing time of the host galaxy. For further discussion of the numerical factors in the definition of a_h, see Merritt (2013). In Figure 1, we plot the evolution of the binary semimajor axes around a = a_h for the γ = 3/2 simulation sample.

Zoom In Zoom Out Reset image size

The rapid shrinking of the SMBH binaries is accompanied by a sudden decline in the stellar density around the binary, transforming the central stellar cusp of the host galaxy into a flat low-density core (e.g., Milosavljević & Merritt 2001). The evolution of the stellar density profiles in simulations γ-1.5-BH-0 to γ-1.5-BH-6 during the core scouring process is presented in Figure 2. At the time when the SMBHs enter the influence radius r_h of each other, the central stellar cusps of the progenitor galaxies have not yet merged. In order to study the disruption of the stellar cusps, we compute the radial stellar density profiles centered at only one of the SMBHs, with this choice remaining fixed throughout the simulations. Calculating the profile with respect to the center of mass of the SMBHs would yield unphysical results for the density of the central region when the stellar cusps are still separated by several kpc and the center of mass of the SMBHs lies between the high-density cusps. In the simulation without SMBHs, the center point for the density profile computation is found using the shrinking sphere algorithm, taking into account only stars from one of the progenitor galaxies.

Zoom In Zoom Out Reset image size

The single panel at the top of Figure 2 shows that in the simulations without SMBHs, the stellar density profiles remain practically unchanged after the coalescence of the stellar cusps of the progenitor galaxies. This is consistent with earlier simulation results that showed that the flattening of the density profiles is expected to be mild in collisionless mergers of cuspy progenitors (e.g., Boylan-Kolchin & Ma 2004). The other panels in Figure 2 show the evolution of the density profiles in the runs with SMBHs, where the separation of the SMBHs is marked with an open circle for the cases where the separation is larger than the plotting range of 10 pc.

The shrinking binaries displace stellar mass from the central region, with the final mass deficit scaling with the mass of the SMBH binary (Merritt 2006). This is clearly seen especially in the bottom panels with the more massive progenitor SMBHs, for which a flattened central section appears in the density profile when the semimajor axis of the binary shrinks from a = r_h to a = a_h. The extent of the flat core and the amount of displaced stellar mass increases systematically as a function of increasing initial SMBH mass. Most of the stellar mass to be displaced from the central region (area between the blue and red lines in Figure 2) is ejected before a = a_h. The density profiles of the run γ-1.5-BH-0 without SMBHs are presented between t = 280 and 363 Myr, which coincides approximately with the moments of the SMBH binary becoming bound and subsequently hard in the simulations containing SMBHs.

The decline of the central stellar density can also be characterized by the increasing size of the sphere of influence r_SOI of the SMBHs. Here we define r_SOI as the radius enclosing a stellar mass equaling the SMBH mass, M_•. The evolution of the sphere of influence of an individual SMBH during the core scouring process is presented in Figure 3. The color-coding and symbols describing the different simulations are the same as in Figure 1. Before the SMBHs form a gravitationally bound binary, the extent of the sphere of influence fluctuates as a response to the changing central stellar density caused by the close passages of the nuclear stellar cusps. After a = r_h (open diamonds), the sphere of influence r_SOI rapidly increases as the stellar density around the shrinking SMBH binary strongly declines. This is especially clearly seen in the simulations with more massive SMBHs, which result in larger mass deficits. The phase of rapid evolution of the r_SOI lasts only until the binary becomes hard (open circles), after which the rate of evolution of r_SOI is much slower.

Zoom In Zoom Out Reset image size

The evolution of the two-dimensional surface densities of the merger remnants is presented in Figure 4 around the time of coalescence of the central stellar cusps. The four rows of Figure 4 show the stellar surface density in the simulation run γ-1.5-BH-0 without SMBHs and in the run γ-1.5-BH-6 with the most massive SMBHs. The spatial extent of the panels is 20 kpc in the two top rows and 2 kpc in the two bottom rows. In the run without SMBHs, the central stellar surface density remains high after cusp coalescence, as seen in the first and third rows of Figure 4, with some stars ejected in shells during the pericenter passages of the cusps. In the simulations that include SMBHs, the evolution is qualitatively different, with the central stellar surface density strongly declining due to the shrinking binary (fourth row in Figure 4). The stars ejected from the central region end up in the outer parts of the merger remnant and are visible as a diffuse stellar component, as can be seen in the rightmost panel of the second row in Figure 4.

Zoom In Zoom Out Reset image size

Next, we study the evolution of the stellar velocity dispersion and velocity anisotropy in the center of the merger remnant. Following T16, we define the central region as r < 0.19 kpc. In spherical coordinates, the components of the velocity dispersion can be expressed as σ_r, σ_θ, and σ_ϕ. The two angular dispersions can be combined into the tangential velocity dispersion component σ_t, defined as

$$\begin{eqnarray}&&{\sigma }_{{\rm{t}}}=\sqrt{\displaystyle \frac{{\sigma }_{\theta }^{2}+{\sigma }_{\phi }^{2}}{2}}.\end{eqnarray} \tag{ 15 }$$

The most commonly used parameter to describe the velocity anisotropy structure of a stellar system is the velocity anisotropy β (Binney & Tremaine 2008), defined as

$$\begin{eqnarray}&&\beta =1-\displaystyle \frac{{\sigma }_{\theta }^{2}+{\sigma }_{\phi }^{2}}{2{\sigma }_{{\rm{r}}}^{2}}=1-\displaystyle \frac{{\sigma }_{{\rm{t}}}^{2}}{{\sigma }_{{\rm{r}}}^{2}}.\end{eqnarray} \tag{ 16 }$$

The velocity anisotropy β is connected to the orbits of the stellar population. If all the stellar orbits are radial, σ_t = 0 and β = 1. On the other hand, for a stellar population on purely circular orbits, σ_r = 0 and $\beta =-\infty$. The SMBH binaries can interact strongly with stars that have small pericenter distances, corresponding to low angular momenta. Consequently, SMBH binaries are expected to primarily eject stars on radial orbits and produce a tangentially biased (β < 0) stellar population in the core region. The main aspects of this model have been verified by both numerical simulations (Quinlan & Hernquist 1997; Milosavljević & Merritt 2001) and observations (e.g., Thomas et al. 2014).

The evolution of the central radial and tangential velocity dispersions (σ_r, σ_t), as well as the velocity anisotropy β, is presented for the γ = 3/2 simulation sample in Figure 5. Initially, the radial and tangential velocity dispersions are equal in the vicinity of the SMBHs, resulting in β ∼ 0 before the SMBHs form a bound binary. As the binary forms and shrinks from the influence radius a = r_h toward the hard separation a = a_h, both the central radial and tangential velocity dispersions increase as more mass is brought into the center of the galaxy. The velocity anisotropy remains β ∼ 0 in this evolutionary phase.

Zoom In Zoom Out Reset image size

The radial velocity dispersion σ_r reaches its maximum value almost coincidentally with the SMBH binary becoming hard. After this, the tangential velocity dispersion remains roughly constant while the radial velocity dispersion begins to decline. The decline is stronger for galaxies with higher initial SMBH masses. Consequently, a tangentially biased region with β < 0 is formed in the central region of the galaxy. We performed an additional test simulation to confirm that the SMBH binary scouring indeed causes the decline of the central β, and that it is not due to numerical effects. The SMBHs of simulation γ-1.5-BH-6 were merged by hand when β = −0.4 and the run was restarted. The evolution of the velocity anisotropy in this run is presented with a dotted line in the bottom panel of Figure 5. The decline of β ceases immediately if the SMBHs are merged by hand, confirming that the hard, shrinking SMBH binaries cause the decline of the central velocity anisotropy.

A star that experiences a strong interaction with the SMBH binary is typically ejected with a velocity v_⋆ comparable to the circular orbital velocity of the binary,

$$\begin{eqnarray}&&{v}_{\star }\sim {V}_{\mathrm{bin}}={\left(\displaystyle \frac{2{{GM}}_{\bullet }}{a}\right)}^{1/2},\end{eqnarray} \tag{ 17 }$$

even though the distribution f(v_⋆) of the ejection velocities is very broad (Valtonen & Karttunen 2006; Merritt 2013). The three-body slingshots become important when the typical ejected star is able to completely escape the host galaxy. A commonly used criterion for escape from a stellar system in N-body studies is ${v}_{\mathrm{esc}}=2\sqrt{3}{\sigma }_{\star }$, simplified from $\langle {v}_{\mathrm{esc}}^{2}\rangle =12\langle {\sigma }_{\star }^{2}\rangle$, which assumes only the virial theorem and an isotropic velocity dispersion tensor (Spitzer 1987). For the merger remnants in our simulations, this criterion yields v_esc ∼ 1000 km s⁻¹. In order to study when this occurs in our simulations, we insert the definitions of the influence radius and hard separation from Equations (9) and (13) into the definition of the circular binary velocity (Equation (17)). We obtain the typical ejection velocities

$$\begin{eqnarray}{v}_{\star }=\left\{\begin{array}{ll}\sqrt{2}{\sigma }_{\star }\sim 410\,\mathrm{km}\,{{\rm{s}}}^{-1}, & a={r}_{{\rm{h}}}\\ 4{\sigma }_{\star }\sim 1160\,\mathrm{km}\,{{\rm{s}}}^{-1}, & a={a}_{{\rm{h}}}.\end{array}\right.\end{eqnarray} \tag{ 18 }$$

At a = r_h, the mean stellar ejection velocity is still low, and the slingshot mechanism is ineffective. The hard binaries eject stars at typical velocities of v_⋆ = 4σ_⋆ > v_esc; thus, the slingshot mechanism is very effective. When the typical ejection velocity equals the escape velocity, v_⋆ = v_esc, the semimajor axis of the binary is a = r_h/6 ≲ a_h. This is consistent with the radial velocity dispersion σ_r reaching its maximum value a few Myr before the binaries become hard, as can be seen in the top panel of Figure 5. After this point, the binary efficiently destroys radial orbits. As stars on more circular orbits cannot interact strongly with the SMBH binary (pericenter distance p ≫ a), the tangential velocity dispersions of the merger remnants remain relatively unaffected.

The results in Sections 4.1 and 4.2 indicate that the formation of the low-density core region and the development of the tangentially biased velocity dispersion occur mainly at different times during the evolution of the merger remnant. Most of the stellar mass to be displaced from the core region is removed before the SMBH binary becomes hard at a = a_h, whereas most of the central negative velocity anisotropy β is built up after this stage. In addition, the timescales of the two processes are very different. The central stellar density declines on a timescale of the order of tens of Myr, which is below the crossing time of the merger remnant. On the contrary, the buildup of the tangentially biased velocity distribution occurs on a timescale of ∼a few hundred Myr, which is several times the crossing time of the merger remnant.

Using KETJU, we include in this study PN corrections up to order 2.5 in the equations of motion of the SMBHs. The PN corrections, and thus the GW emission effects, are still negligibly small when a ∼ a_h, excluding the possible but very rare head-on collisions of SMBHs. Eventually, depending on the SMBH binary mass and eccentricity, the PN radiative loss terms in the equations of motion become important. This occurs at the latest when a ∼ a_GW ∼ 0.01 × a_h (Quinlan 1996), yielding a_GW ∼ 0.1–1.1 pc for the binaries in our sample. The binary emits the rest of the orbital angular momentum and energy away as GWs, resulting in an SMBH merger and the formation of a single SMBH.

Averaged over the orbital period, the rates of change of the binary semimajor axis and eccentricity are (Peters & Mathews 1963)

$$\begin{eqnarray}{\left\langle \displaystyle \frac{{da}}{{dt}}\right\rangle }_{\mathrm{GW}} & = & -\displaystyle \frac{64}{5}\displaystyle \frac{{G}^{3}{M}_{1}{M}_{2}({M}_{1}+{M}_{2})}{{c}^{5}{a}^{3}}\displaystyle \frac{1+\tfrac{73}{24}{e}^{2}+\tfrac{37}{96}{e}^{4}}{{(1-{e}^{2})}^{7/2}}\\ {\left\langle \displaystyle \frac{{de}}{{dt}}\right\rangle }_{\mathrm{GW}} & = & -\displaystyle \frac{304}{15}e\displaystyle \frac{{G}^{3}{M}_{1}{M}_{2}({M}_{1}+{M}_{2})}{{c}^{5}{a}^{4}}\displaystyle \frac{1+\tfrac{121}{304}{e}^{2}}{{(1-{e}^{2})}^{5/2}} & ,\end{eqnarray} \tag{ 19 }$$

assuming that the binary evolution is driven purely by GW emission in the PN order PN2.5. Here M₁ and M₂ are again the masses of the SMBHs, a is the semimajor axis, and e is the eccentricity of the binary. An important property of Peters' formulae (Equation (19)) is the extreme dependence of the GW-driven binary evolution on the binary eccentricity. Both $\langle \dot{a}\rangle$ and $\langle \dot{e}\rangle$ diverge in the limit e → 1.

Most of the simulations are run until the SMBHs merge. However, as our main priority lies in studying the core scouring process in the simulated galaxies, we do not continue simulation runs with low-eccentricity binaries beyond t = 2 Gyr. The maximum possible SMBH coalescence timescale in our simulation sample is t_coal ∼ 6 Gyr, assuming the lowest hardening rate among the simulations (run γ-1.0-BH-6) and a circular binary with e = 0 (R17).

The final surface brightness profiles μ(r) of the merger remnants are presented in Figure 6. As in T16, we assume a constant stellar mass-to-light ratio of M_⋆/L = 4.0 in order to compare the simulated surface density profiles with actual observations. The profiles are azimuthally averaged over 100 random viewing angles. The outer parts of the simulated merger remnants appear very similar, as expected from the Dehnen profile ICs and identical merger orbits. The surface brightness of the simulated galaxies falls below the observed profile of NGC 1600 at r ∼ 30 kpc.

Zoom In Zoom Out Reset image size

Studying the central parts of the simulated merger remnants, the surface brightness profiles in the runs that include SMBHs turn almost flat at 0.1 kpc ≲ r_b ≲ 1 kpc. The profiles of the merger remnants without SMBHs remain cuspy until the very center. We also immediately see that the central surface brightness is lower for the γ = 1 ICs compared to the γ = 3/2 runs with the same SMBH mass. In addition, the central surface brightness decreases and the core size (extent of the flat section) increases as the initial SMBH mass increases in both progenitor galaxy samples. This is consistent with the mass deficit arguments in the literature, which state that the mass deficit is proportional to the total mass of the SMBH binary (e.g., Merritt 2006).

The agreement with the observations is best in the simulations with a final SMBH mass similar to M_• = 1.7 × 10¹⁰ M_⊙ inferred for NGC 1600. However, there is some ambiguity here because of the assumed stellar mass-to-light ratio M_⋆/L = 4.0 (van Dokkum et al. 2017). With larger values for the mass-to-light ratio M_⋆/L > 4.0, simulated galaxies with lower M_• or steeper initial stellar density profiles could also match the central surface brightness of NGC 1600. However, with M_⋆/L > 4.0, the outer parts of the simulated galaxies would end up too dim compared to the observations. As NGC 1600 is in a group environment, the galaxy has most likely experienced several minor mergers that have puffed up and brightened the outer parts of the galaxy, as opposed to our isolated merger models (e.g., Hilz et al. 2012, 2013). Based solely on surface brightness data, we cannot exactly state which one of the simulated galaxies most closely resembles the observed NGC 1600, and additional velocity data are required for this.

The final velocity anisotropy profiles β(r) of the simulated galaxies are presented in Figure 7. The velocity anisotropy is computed in radial bins, with the bin width Δr = 0.19 kpc motivated by the spatial resolution of the spectral data in T16. For comparison, the isotropic β profiles of the progenitor galaxies are also shown as horizontal lines. The top panels compare the β(r) of the simulated galaxies to the dynamical model of NGC 1600 (T16). Without the inclusion of SMBHs, the velocity structure of the merger remnants of both the γ = 1 and the γ = 3/2 samples is close to isotropic in the center of the galaxy.

Zoom In Zoom Out Reset image size

In the outer parts, the stellar population is dominated by radial orbits (β ∼ 0.40–0.45 at r = 10 kpc), which is typical for collisionless mergers (Röttgers et al. 2014). The outer parts of the merger remnants with SMBHs are almost identical to the remnants without SMBHs. The dynamical model of NGC 1600 is more radially biased at radii beyond r ≳ 3 kpc, with β ∼ 0.5–0.6 at r = 10 kpc. The smaller amount of stellar material on radial orbits in the outer parts can be explained by the lack of minor mergers in our isolated major-merger simulation, as already discussed in the previous section. However, it should be noted that the velocity structure of the stellar population in the outer parts of the simulated merger remnants still shows qualitatively similar trends to the dynamical model of NGC 1600, at least when compared to the isotropic progenitor galaxies.

The central parts (r ≲ 1 kpc) of the simulated merger remnants have a tangentially biased stellar population (β < 0). There is a clear trend of increasingly more negative β with increasing initial SMBH mass. In addition, more cuspy initial stellar profiles produce steeper β profiles in the core region, consistent with the results of Bortolas et al. (2018). The dynamical model of NGC 1600 has an even steeper β profile near the center when compared to the simulations. This suggests that the initial stellar density profiles for the progenitors of NGC 1600 were even steeper than γ = 3/2. We note that the final central values of β depend on the chosen isotropic ICs. A second generation of mergers would already have β < 0 initially, leading to an increasingly tangential central stellar population. In addition, other physical processes, such as the adiabatic growth of the single central SMBH, may decrease β as well (e.g., Goodman & Binney 1984; Thomas et al. 2014).

The bottom panels of Figure 7 present the same simulated anisotropy profiles as the top panels, but the radial coordinate is now scaled with the break radius r_b of the surface brightness profile. The determination of the break radius r_b is discussed in Section 5.2. The β profiles of the dynamical models of six observed core galaxies from the SINFONI black hole survey (Thomas et al. 2014; Saglia et al. 2016) are presented for comparison. Observed massive elliptical galaxies with cores have remarkably similar velocity anisotropy profiles when scaled with the core radius r_b. In the outer parts (r/r_b > 2–3), the β profiles of the simulated merger remnants and the galaxies from the SINFONI black hole survey now agree very well. At even larger radii (r/r_b > 10), the simulated profiles agree with the lower end of the distribution of observed profiles. In the core region (r/r_b < 1), the profiles of the simulated merger remnants coincide with the upper end of the observed anisotropies in the SINFONI black hole survey. Increasing the initial SMBH mass and selecting a steeper stellar density profile slope results in merger remnants in better agreement with the dynamical models of the galaxies from the SINFONI black hole survey. Again, we note that the assumption of a single generation of galaxy mergers with isotropic ICs affects the results discussed here. An additional generation of mergers would start with β < 0 in the center and β > 0 in the outer parts of the galaxies, potentially leading to a better agreement with the observed galaxies.

In summary, the results of our surface brightness profile and velocity anisotropy analysis support the SMBH mass M_• = 1.7 × 10¹⁰ M_⊙ obtained by the dynamical modeling of T16 and suggest steep initial stellar profiles with γ ≥ 3/2 for the progenitor galaxies of NGC 1600.

We show in Figures 8 and 9 the two-dimensional kinematics maps for the γ = 1 and 3/2 simulation samples, respectively. The maps are constructed similarly to what an observer would see if the merger remnants were observed with an integral field unit spectrograph (see Naab et al. 2014). Following the same procedures used in analyzing observational data, each panel is divided into "spaxels" with the same signal-to-noise, using the Voronoi tessellation algorithm (Cappellari & Copin 2003). In our case, the signal-to-noise is represented by the number of particles in each bin. For each row in the two figures, the four panels show, from left to right, the average LOS velocity V_avg, the LOS velocity dispersion σ, and the h₃ and h₄ parameters, which represent the skewness and kurtosis of the LOS velocity distribution (van der Marel & Franx 1993; Gerhard 1993). These parameters are calculated by fitting the histogram of LOS velocities in every spaxel with the following modified Gaussian function:

$$\begin{eqnarray}&&f(v)={I}_{0}\,{e}^{-{y}^{2}/2}(1+{h}_{3}\,{H}_{3}(y)+{h}_{4}\,{H}_{4}(y)),\end{eqnarray} \tag{ 20 }$$

where y = (v − V_avg)/σ and I₀ is a normalization constant. Here H₃ and H₄ are the Hermite polynomials of third and fourth order:

$$\begin{eqnarray}&&{H}_{3}(y)=(2\,\sqrt{2}\,{y}^{3}-3\,\sqrt{2}\,y)/\sqrt{6},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{H}_{4}(y)=(4\,{y}^{4}-12\,{y}^{2}+3)/\sqrt{24}.\end{eqnarray} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before performing this analysis, the galaxies are oriented using the reduced inertia tensor so that their major axis is on the x axis of the figure, and the LOS velocities are centered so that the spaxels within the central kpc in projection have zero average velocity. We also plot luminosity contours with a spacing of 1 mag on top of each velocity map.

In the merger remnants without SMBHs, the central regions show a positive h₄ peak. This is typical for nearly isotropic galaxies with SMBHs and not too steep central stellar cusps (Baes et al. 2005). However, in simulations including SMBHs, we find systematically more negative h₄ values as a function of increasing SMBH mass (see bottom rows of Figures 8 and 9). The size of the negative h₄ region also grows with increasing SMBH mass and size of the sphere of influence. The lower values of h₄ likely result from the missing weakly bound stars on low angular momentum orbits.

In the velocity maps, a central kpc-scale, decoupled rotating region becomes increasingly more visible as the SMBH mass increases. There are even indications of more complex kinematic subsystems in the merger remnants γ-1.5-BH-4 to γ-1.5-BH-6. As can be seen in the bottom panels of Figure 9, a large (∼5 kpc) region counterrotating to the central rotating region becomes apparent in the simulations containing the more massive SMBHs. The h₃ map shows an anticorrelation to the LOS velocity map, which is commonly observed in rotating systems (e.g., Bender et al. 1994; Krajnović et al. 2011; Veale et al. 2017). Finally, the stellar velocity dispersion also increases because of the increased mass in the central region for larger SMBH masses.

In conclusion, Figures 8 and 9 clearly demonstrate that the effects of black hole scouring are indeed observable in the kinematic maps of quiescent merger remnants.

In T16, tight scaling relations are presented between the observed galaxy core sizes (or break radii) r_b, SMBH masses M_•, and spheres of influence of the SMBHs r_SOI. We perform a similar analysis for our simulated merger remnants. In addition, we use a simple core scouring model to interpret the results of our merger simulations.

It should be noted that while our merger simulation sample consists of single-mass remnants with ${M}_{\star }={M}_{\star }^{\mathrm{NGC}\,1600}\,=8.3\times {10}^{11}\,{M}_{\odot }$, the scaling relations of T16 are obtained from a broader core galaxy sample. The core galaxies in the T16 sample have total stellar masses in the range from M_⋆ = 5.0 × 10¹⁰ to 2.0 × 10¹² M_⊙, with a mean stellar mass of $\langle {M}_{\star }\rangle =5.8\times {10}^{11}\,{M}_{\odot }$ (Rusli et al. 2013a, 2013b; Saglia et al. 2016). The mean effective radius is approximately $\langle {R}_{{\rm{e}}}\rangle =8.5\,\mathrm{kpc}$. The velocity dispersions of the observed sample are in the range 168 km s⁻¹ < σ_⋆ < 380 km s⁻¹ with a mean of $\langle {\sigma }_{\star }\rangle =289\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This agrees reasonably well with the velocity dispersions of our simulated merger remnants, for which σ_⋆ ∼ 290 km s⁻¹. In the T16 sample, the SMBH masses span from 4.2 × 10⁸ to 2.1 × 10¹⁰ M_⊙ with an average of $\langle {M}_{\bullet }\rangle =3.5\times {10}^{9}\,{M}_{\odot }$. This is in a similar range to the final SMBH masses in our merger simulations: 1.7 × 10⁹ M_⊙ < M_•< 1.7 × 10¹⁰ M_⊙.

Even though the stellar mass and velocity dispersion in our simulated galaxies lie close to the mean values of the T16 observational sample, we stress that some care should be taken in the interpretation of our simulated core scaling relations. The interpretation boils down to the central question of how local the core scouring process is. Thus, the question is whether the properties of the host galaxy beyond r_SOI strongly affect the core formation process or whether the core formation process is mainly dictated by the mass of the SMBH. The results in Section 4.1 support the latter scenario, as the central core forms rapidly compared to the crossing time of the merger remnant, when the SMBH binary shrinks from a = r_h to a = a_h. However, on the other hand, if the large-scale properties of the host galaxies strongly affect the core scouring process, we should not be able to reproduce the observed scaling relations using only single-mass progenitor galaxy models.

There are a number of widely used methods to determine the size of the cored region in the surface brightness profile μ(r). One popular model is the six-parameter core-Sérsic profile (Graham et al. 2003; Trujillo et al. 2004), which is defined as

$$\begin{eqnarray}&&\mu (r)=\mu ^{\prime} {\left[1+{\left(\displaystyle \frac{{r}_{{\rm{b}}}}{r}\right)}^{\alpha }\right]}^{\gamma /\alpha }{e}^{-b{[({r}^{\alpha }+{r}_{{\rm{b}}}^{\alpha })/{R}_{{\rm{e}}}^{\alpha }]}^{1/(\alpha n)}}.\end{eqnarray} \tag{ 23 }$$

In the outer parts, the profile follows the ordinary three-parameter Sérsic profile (Sérsic 1963; Caon et al. 1993) defined by the Sérsic index n, effective radius R_e, and overall normalization μ'. Toward the central regions, the core-Sérsic profile breaks away from the Sérsic profile at the break radius (i.e., the core size) r_b to allow for a shallow central power law with a slope γ. The parameter α determines how sudden this break is. Finally, the constant b is determined by requiring that the effective radius R_e encloses half of the total light of the core-Sérsic profile.

We perform the profile fitting using the CMPFIT library.⁴ The radial fitting range for all of the surface brightness profiles shown in Figure 6 is 0.01 kpc < r < 62 kpc in this section, similar to T16. For all of these fits, the core-Sérsic index is fixed to n = 4. The reason for this procedure is to avoid any possible degeneracy between the Sérsic index and the break radii in the fitting procedure (Dullo & Graham 2012) and to compare the results from subsamples γ = 1 and 3/2 with each other in a more straightforward and robust manner. If we allow the Sérsic index to freely vary, we typically get n ≤ 3 for the γ = 1 simulation sample, which is on the low side for massive early-type galaxies; see, e.g., Appendix A of Naab & Trujillo (2006) for details. The initial surface brightness profiles of the γ = 3/2 sample are best fit with the Sérsic indices close to n = 4 (Dehnen 1993), making these ICs more in tune what would be expected for massive early-type galaxies.

Another option is to fit the so-called Nuker profile (Lauer et al. 1995) to the surface brightness profile, defined as

$$\begin{eqnarray}&&\mu (r)={\mu }_{{\rm{b}}}{2}^{(\beta -\gamma )/\alpha }{\left(\displaystyle \frac{r}{{r}_{{\rm{b}}}}\right)}^{-\gamma }{\left[1+{\left(\displaystyle \frac{r}{{r}_{{\rm{b}}}}\right)}^{\alpha }\right]}^{(\gamma -\beta )/\alpha },\end{eqnarray} \tag{ 24 }$$

where r_b is the break radius, γ is the logarithmic slope of the profile inside, and β is the logarithmic slope outside the break radius. The parameter α again measures the sharpness of the transition from the inner to the outer power law at the break radius, r = r_b. The Nuker profile was originally used to fit the Hubble Space Telescope observations of early-type galaxies within the central 10''. At the distance of NGC 1600, this corresponds to ∼3.1 kpc. Later, the Nuker profile was also used to fit the surface brightness profiles in the entire radial range of observed galaxies. However, it is well known that the parameters of the Nuker profile depend on the radial fitting range (Graham et al. 2003).

We performed the Nuker profile fits first in the same radial range as for the core-Sérsic fits (0.01 kpc < r < 62 kpc). In this case, the best-fit Nuker profiles for all the galaxies in our simulated galaxy sample have α < 1. In general, when α ≲ 1, the presence of two power laws may not emerge, and instead, one continuous curving arc describes the profile (Graham et al. 2003). Consequently, the break between the inner and outer profiles is not well defined. We tested whether stacking of the surface brightness profiles from different viewing angles or using elliptical rather than circular bins in Figure 6 would affect the results, but we still found α < 1 for all the Nuker fits. However, we found that when the outer edge of the radial fitting range is chosen to be below 2 × R_e, the sharpness parameter α becomes larger than unity. The core size r_b thus also depends on the radial fitting range, which is consistent with findings in the literature (Graham et al. 2003). Next, we performed the Nuker fits only in the central parts of the surface brightness profile in the range 0.01 kpc < r < 3.1 kpc. Now we found α > 1 for all the simulated galaxies except γ-1.5-BH-2. The core sizes obtained from the core-Sérsic and central Nuker fits are compared in Figure 10. We find that the Nuker break radii are systematically larger than the core-Sérsic core radii, especially in the simulation set with γ = 1.

Zoom In Zoom Out Reset image size

Finally, a simple nonparametric method can also be used to estimate the core size by studying the logarithmic slope of the surface brightness profile. The "cusp radius" r_γ can then be defined as

$$\begin{eqnarray}&&{\left[\displaystyle \frac{d\mathrm{log}\mu (r)}{d\mathrm{log}r}\right]}_{r={r}_{\gamma }}=-\displaystyle \frac{1}{2},\end{eqnarray} \tag{ 25 }$$

i.e., as the point where the logarithmic slope of the surface brightness profile, μ(r), equals −1/2 (Carollo et al. 1997; Lauer et al. 2007). We find that the cusp radii r_γ of our simulated galaxies are consistent with the core radii derived from the core-Sérsic fits (see Figure 10).

In this study, we use the core-Sérsic break radius r_b to measure the sizes of the cores of our simulated galaxies. This selection makes the comparison with T16 straightforward. The Nuker fits are not used in the following sections because of their dependence on the radial fitting range and the α < 1.0 problem encountered while fitting our simulated surface brightness profiles. The cusp radii r_γ could also be used in the following analysis, as the core sizes using the method are consistent with the results obtained using the core-Sérsic fit, which was also found by T16.

We use a simple core scouring model to interpret the results of our simulations. Starting from the stellar density profile ρ(r) in our ICs, we first remove a mass deficit equaling the mass of the SMBH, i.e., M_def = M_•. This choice is motivated by both observations and N-body simulations (see, e.g., Merritt 2013). The stellar mass is removed from the density profile in such a way that the central part of the scoured density profile ρ_s(r) becomes flat up to the deficit radius r_def:

$$\begin{eqnarray}{\rho }_{{\rm{s}}}(r)=\left\{\begin{array}{ll}\rho (r), & r\geqslant {r}_{\mathrm{def}}\\ \rho ({r}_{\mathrm{def}}), & r\lt {r}_{\mathrm{def}}.\end{array}\right.\end{eqnarray} \tag{ 26 }$$

The mass deficit is therefore

$$\begin{eqnarray}&&{M}_{\mathrm{def}}=4\pi {\int }_{0}^{{r}_{\mathrm{def}}}[\rho (r)-\rho ({r}_{\mathrm{def}})]{r}^{2}{dr},\end{eqnarray} \tag{ 27 }$$

from which the mass deficit radius can be solved numerically. The deficit radius r_def is sometimes also referred to as the core radius (e.g., Volonteri et al. 2003), but we reserve this label for the break radius r_b obtained from the core-Sérsic fit (Equation (23)) of the surface brightness profile.

Next, we measure the extent of the sphere of influence of the SMBH from the flattened density profile ρ_s(r). Instead of using the definition of r_h from Equation (9), we follow T16 and define the sphere of influence r_SOI using the enclosed mass as

$$\begin{eqnarray}&&{M}_{\star }({r}_{\mathrm{SOI}})=4\pi {\int }_{0}^{{r}_{\mathrm{SOI}}}{\rho }_{{\rm{s}}}(r){r}^{2}{dr}={M}_{\bullet }.\end{eqnarray} \tag{ 28 }$$

The sphere of influence r_SOI is then again most conveniently solved numerically.

The two-dimensional surface brightness profile μ(r) is obtained by projecting the analytically scoured density profile. Labeling the constant stellar mass-to-light ratio as M_⋆/L = ϒ_⋆, the surface brightness profile can be computed from

$$\begin{eqnarray}&&\mu (r)={{\rm{\Upsilon }}}_{\star }{\rm{\Sigma }}(r)=2{{\rm{\Upsilon }}}_{\star }{\int }_{r}^{\infty }\displaystyle \frac{{\rho }_{{\rm{s}}}(r^{\prime} )r^{\prime} }{\sqrt{{r}^{{\prime} 2}-{r}^{2}}}{dr}^{\prime} .\end{eqnarray} \tag{ 29 }$$

Finally, the break radii r_b of the generated surface brightness profiles are obtained from the core-Sérsic fitting procedure, as described in the previous section.

We derive the break radii r_b of the surface brightness profiles of our simulated merger remnants using the fitting procedure described in Section 5.2. The residuals of the fits are of the same order as the residuals of the observed core-Sérsic of NGC 1600 presented in T16, that is, below Δμ < 0.1 mag. The extents of the spheres of influence r_SOI are obtained directly from the simulation data using the enclosed stellar mass definition (see Equation (28)). These parameter values comprise the core data points (r_b, r_SOI), (r_b, M_•) of the merger simulations.

We produce another set of (r_b, r_SOI), (r_b, M_•) data points using our simple scouring model and the 12 ICs of the merger simulations containing central SMBHs. The actual scaling relations for both the simulation samples and the analytical scouring models are obtained by fitting the free parameters a₁, b₁, a₂, b₂ in the fitting functions

$$\begin{eqnarray}{\mathrm{log}}_{10}({r}_{\mathrm{SOI}}/\mathrm{kpc}) & = & {a}_{1}+{b}_{1}{\mathrm{log}}_{10}({r}_{{\rm{b}}}/\mathrm{kpc})\\ {\mathrm{log}}_{10}({M}_{\bullet }/{M}_{\odot }) & = & {a}_{2}+{b}_{2}{\mathrm{log}}_{10}({r}_{{\rm{b}}}/\mathrm{kpc}).\end{eqnarray} \tag{ 30 }$$

The values of the best-fit parameters are listed in Table 3, where the error estimates have been computed using a simple bootstrap method. The results from both the merger runs and the simple scouring model are presented in Figure 11. The simulation data points and simple scouring model points are plotted as red (γ = 1) and blue (γ = 1.5) filled circles, respectively. The scaling relations for both of the γ = 1 and 3/2 subsamples are plotted as dotted lines of the corresponding color. The observed galaxy sample of T16 is shown as black squares, with the solid lines representing the observed scaling relations log₁₀(r_SOI/kpc) = (−0.01 ± 0.29) + (0.95 ± 0.08)log₁₀(r_b/kpc) and log₁₀(M_•/M_⊙) = (10.27 ± 0.51)+ (1.17 ± 0.14)log₁₀(r_b/kpc).

Zoom In Zoom Out Reset image size

Table 3.  The Parameter Values of the Core Scaling Relations (r_b vs. r_SOI) and (r_b vs. M_•) Presented in Figure 11

${\mathrm{log}}_{10}({r}_{\mathrm{SOI}}/\mathrm{kpc})={a}_{1}+{b}_{1}{\mathrm{log}}_{10}({r}_{{\rm{b}}}/\mathrm{kpc})$
	a1	b1
Thomas et al. (2016)	−0.01 ± 0.29	0.95 ± 0.08
Simple model γ = 1	−0.05 ± 0.01	2.10 ± 0.08
Simulation γ = 1	0.25 ± 0.03	2.55 ± 0.22
Simple model γ = 3/2	0.06 ± 0.01	1.06 ± 0.01
Simulation γ = 3/2	0.33 ± 0.07	1.35 ± 0.17

${\mathrm{log}}_{10}({M}_{\bullet }/{M}_{\odot })={a}_{2}+{b}_{2}{\mathrm{log}}_{10}({r}_{{\rm{b}}}/\mathrm{kpc})$
	a2	b2
Thomas et al. (2016)	10.27 ± 0.51	1.17 ± 0.14
Simple model γ = 1	10.11 ± 0.01	3.78 ± 0.19
Simulation γ = 1	10.55 ± 0.09	4.35 ± 0.58
Simple model γ = 3/2	10.47 ± 0.01	1.52 ± 0.01
Simulation γ = 3/2	10.72 ± 0.11	1.93 ± 0.28

Notes. The constant term a from both the simple scouring model and the merger simulations is within the uncertainty of the observed relation of T16. We also see a trend for both the analytic and simulated slopes b toward the observed slope values of b ∼ 1 for steeper initial stellar density profiles, with the steeper γ = 1.5 profiles clearly providing a better match with the observations.

Download table as:  ASCIITypeset image

Studying Figure 11, we see that our simple scouring model produces qualitatively similar results to the major-merger simulations. We also see a trend toward the observed relations with increasing initial stellar density profile slope γ, as the scaling relations flatten for steeper initial stellar profiles. This fact points toward very cuspy initial stellar density profiles of the progenitor galaxies, with γ in excess of 3/2, in the framework of a single generation of galaxy mergers. Mergers of galaxies with relatively flat Hernquist-like stellar bulges (γ = 1) result in too-steep core scaling relations. Finally, we also note that our core scaling relations are consistent with the results of Lauer et al. (2007), who obtained the M_•–r_γ relation from a sample of 11 observed core galaxies using directly determined SMBH masses.

The surprisingly good agreement between the scaling relations derived from our cuspier (γ = 1.5) ICs and the observed relations requires an explanation. All the progenitor galaxies in our merger simulation sample had identical stellar density profiles and, in practice, identical velocity profiles outside of r > r_SOI, while the observed T16 galaxy sample had a broad range of stellar mass and velocity dispersions, as described in Section 5.1.

Both observations and simulations have shown that the total density profiles of power-law galaxies are well represented with simple density profiles ρ ∝ r^−γ over a wide range in radius, with γ ∼ 2.0–2.2 (e.g., Remus et al. 2013, 2017; Cappellari et al. 2015 and references therein). At r ≪ r_1/2, the stellar component dominates the total density profile, so the stellar density profile has a central slope γ ∼ 2. In order to estimate the mean central stellar densities ρ of the power-law galaxies, we model the stellar component using the Dehnen profile as in Section 3.1 with γ = 2. We sample a large number of progenitor galaxy models with parameters M_•, M_⋆, and r_1/2 and their uncertainties from T16 and Saglia et al. (2016). Next, we compute the mean central (r < 1 pc) stellar densities for these idealized galaxy models. Although the stellar mass increases by a factor of 40 from $\min ({M}_{\star }^{{\rm{T}}16})$ to $\max ({M}_{\star }^{{\rm{T}}16})$, the central stellar density increases simultaneously only by a factor of 6.2 ± 3.8 from the most diffuse galaxy to the galaxy with the highest central density. Thus, the central stellar density increases slowly compared to the total stellar mass. The population of power-law galaxies is very homogeneous in its central stellar density properties. At the same time, the mass of the SMBH, M_• in the T16 sample, increases much more significantly, by a factor of ∼50.

We argue that this homology of the central properties of the power-law ellipticals is responsible for the good correspondence between the core scaling relations from our single-M_⋆ merger simulation sample and the scaling relations of the observed core galaxy population. Our results support the scenario in which the progenitors of the most massive core ellipticals are cuspy power-law galaxies (γ > 3/2) with central SMBHs. The central properties of the stellar population are fairly homogeneous over a large range in M_⋆, while the SMBH masses steadily increase as M_⋆ grows. Thus, the SMBH mass M_• stands out as the main parameter in the core scouring process and is also responsible for determining the core scaling relations.