Connecting Core Galaxy Properties to the Massive Black Hole Binary Population

We investigate how the properties of massive black hole binaries influence the observed properties of core galaxies. We compare the observed trend in stellar mass deficit as a function of total stellar mass in the core galaxy with predicted trends in IllustrisTNG. We calculate mass deficits in simulated galaxies by applying sub-grid, post-processing physics based on the results of literature N-body experiments. We find the median value of the posterior distribution for the minimum binary mass ratio capable of creating a core is 0.7. For the gas mass fraction above which a core is erased we find a median value of 0.6. Thus low mass ratio binaries do not contribute to core formation and gas-rich mergers must lead to star formation within the nucleus, ultimately erasing the core. Such constraints have important implications for the overall massive black hole binary population, black hole–galaxy co-evolution, and the origin of the gravitational wave background.


INTRODUCTION
The most massive elliptical galaxies have low-density centers ≲ 1 kpc in size (Lauer et al. 2007b), manifesting observationally as shallow inner surface brightness profiles compared to less massive ellipticals and to spiral galaxy classical bulges.The inner profile is significantly flatter than the steep, inward extrapolation of the outer profile (see Figure 1).Systems for which this is true are called "core" galaxies, while galaxies that have a steep inner profile are called "power-laws".Core galaxies are likely to have significantly different evolution, relative to other ellipticals, as they tend to be more spherical, have boxier isophotes, and rotate more slowly (Lauer 2012;Krajnović et al. 2013).
It is thought that massive black hole binaries (MBHBs) are responsible for the formation of cores.Binary black holes are sources of gravitational waves and the precursor to black hole mergers, which is a channel of black hole growth.The connection between supermassive black holes (SMBHs) and their host galaxies is important for understanding their co-evolution, including BH scaling relations and active galactic nuclei (AGN) feedback.
A number of observational facts about ETGs has lead to a theoretical picture of how galaxy cores are formed through MBHB activity.Most massive early-type galaxies (ETGs) like those in question have undergone at least one major merger since  ∼ 1 (Bell et al. 2006).Between 50-70% of the mass growth of ETGs has occurred in the same time-frame as a result of merger events, with about half of the contribution from major mergers (López-Sanjuan et al. 2012).For galaxies with stellar masses ≳ 10 11  ⊙ major mergers play a dominant role in mass assembly since  ∼ 1, with 30% of those ★ E-mail: cordellh@umich.edu(CJH) Stellar mass deficits in core galaxies may be calculated if the surface brightness profile is known and a model of the extrapolated light distribution is created.The light deficit is the integrated difference between the two profiles, which can be converted into a stellar mass deficit by multiplication of a stellar mass-to-light ratio in the appropriate band.The "break radius" indicates the distance from the galaxy center where the profile transitions from the shallow inner region to the steep outer region, representing the size of the core.Based on the profile of NGC 1016.
mergers being gas-poor (Bundy et al. 2009).Accretion of less massive ETGs (i.e., minor mergers) also contributes to mass assembly (Naab et al. 2006).Concerning rotation, the gas content of remnant galaxies plays a crucial role (Yoon et al. 2022;Walo-Martín et al. 2020), with gas-poor mergers tending to reduce the specific angular momentum of the stellar component by an average of ∼ 30%, while gas-rich mergers are shown to increase angular momentum by ∼ 10% (Lagos et al. 2018).Given these properties, the formation of cores is best explained by stellar scattering in the galactic nucleus by MBHBs (Begelman et al. 1980;Faber et al. 1997) in the aftermath of gas-poor, ellipticalelliptical merger events.Once the progenitor galaxies merge, their respective SMBHs are initially unbound, eventually coupling near the remnant galaxy's nucleus via dynamical friction (DF).Once the black holes form a binary, they scatter neighboring stars, overcome depletion of the loss-cone, and begin to emit gravitational waves before coalescing.Thus, to successfully scour out a core, the secondary black hole must join the primary at the bottom of the remnant potential well.
Prescriptions based on previous analyses (e.g., Chen et al. 2022) estimate the efficiency of DF as a function of the mass of the secondary, the central velocity dispersion of the remnant, the distance from the galactic center, the number of nuclear stars present, and the gas mass fraction.The DF timescale increases for larger merger remnants, since the secondary black hole must travel further to reach the primary.A high stellar velocity dispersion also leads to a longer timescale, as only stars moving slower than the secondary can contribute to friction.The DF timescale decreases for more massive secondaries and for a larger gas content.Once the binary is formed, stellar scattering ensues via three-body interactions.Scattering will tend to eject stars from the nucleus at the expense of the black hole binary's orbital energy, shrinking the separation between the black holes.The result is a low density core which is likely to persist for a long period, unless the gas content in the remnant is high enough to trigger star formation in the galactic center, which has the effect of erasing the core (Faber et al. 1997).
If the core is not erased, then it will appear in observations as a luminosity deficit in the galaxy center, which can be converted into a mass deficit given a mass-to-light ratio.N-body experiments have suggested that the mass in ejected stars is proportional to the mass of the central black hole and the number of major mergers, with between one and three major mergers ejecting a mass in stars consistent with observed mass deficits (Quinlan 1996;Milosavljević & Merritt 2001;Merritt 2006).Additionally, MBHB N-body experiments predict a dearth of radial stellar orbits in the nucleus of core galaxies as such orbits have small pericenters, allowing for strong MBHB interactions that kick the stars out of the center.This prediction is seen in detailed orbital modeling of elliptical galaxies with high-spatial resolution spectroscopy (Thomas et al. 2014).
Despite the observational evidence that galaxy cores are formed by MBHBs, there are few confirmed binaries (e.g., Rodriguez et al. 2006;Eracleous et al. 2012;Graham et al. 2015;Runnoe et al. 2015Runnoe et al. , 2017;;Charisi et al. 2016), though it is thought they should exist in large number.One line of evidence is the population of dual AGN, where widely separated, gravitationally unbound pairs are known (Inada et al. 2012;More et al. 2016;Lemon et al. 2018;De Rosa et al. 2019;Foord et al. 2019;Silverman et al. 2020;Foord et al. 2020Foord et al. , 2021;;Tang et al. 2021), though it is unclear if they will ever merge to form binaries.Recently, the North American Nanohertz Observatory for Gravitational waves (NANOGrav) along with other pulsar timing arrays (PTAs) have determined that the common red noise process seen in measurements is caused by gravitational waves (Agazie et al. 2023a;Reardon et al. 2023;Antoniadis et al. 2023).It has been interpreted as the stochastic gravitational wave background (GWB) caused by MBHBs (Agazie et al. 2023b) or other exotic physics (Afzal et al. 2023).Given that MBHBs are thought to also dominate the formation of galaxy cores, it is possible to combine both electromagnetic and gravitational wave signals to lead to a better understanding of the MBHB population.This multimessanger approach is essential for astrophysical interpretations of the GWB, which has degeneracies that can be broken with EM observations (Agazie et al. 2023b).
In this Paper we investigate how MBHB mass ratio distributions leave observable effects via the stellar cores they leave behind.In section 2 we describe the methods used to calculate mass deficits and stellar masses from an observed sample of ETGs.Section 3 details our post-processing of IllustrisTNG simulations to make predictions of mass deficits utilizing native cosmological merger trees.In section 4 we describe the statistical comparison between synthetic data sets generated from IllustrisTNG and the observed sample to obtain constraints on the MBHB mass ratios capable of creating a core and the maximum gas mass fraction allowed in a remnant before the core is erased.Lastly, we discuss our results in section 5.

METHODS: OBSERVATIONAL DATA
In this section we describe the the observed sample and how we derive the stellar mass of each galaxy and its mass deficit.
The upper limit for the angular size of a galaxy core in the local universe is on the order of a few arcseconds (Faber et al. 1997), requiring high spatial resolution imaging.As such, our parent sample consists of 219 local, ETGs with Hubble Space Telescope (HST) observations.This sample was compiled by Lauer et al. (2007b) using data from Ravindranath et al. (2001); Lauer et al. (2005); Quillen et al. (2000); Laine et al. (2003); Lauer et al. (1995) and Rest et al. (2001).Each galaxy within the sample has reported absolute V-band magnitudes,   , morphological types, and distances, along with derived Nuker-profile parameters.The Sérsic profile (Sérsic 1963) is typically used when information about the entire galaxy, including outer envelopes, is of interest.The Nuker profile (Lauer et al. 1995), on the other hand, was developed with the inner galaxy in mind and best reveals inner cores.See Graham et al. (2003); Trujillo et al. (2004); Hopkins et al. (2009) for further discussion on the pros and cons of each profile.The Nuker profile has the form where  is the semi-major axis on the plane of the sky,   is the break radius,   is the surface brightness at the break radius,  and  are the negative logarithmic surface brightness slopes at large and small radii respectively, and  determines the sharpness of the transition between inner and outer profiles.The sample is split into three categories based on the inner slope  of the surface brightness profile.Galaxies for which  < 0.3 are "cores", and galaxies that have  > 0.5 are "power-laws".Galaxies with inner slope values between 0.3 and 0.5 are referred to as "intermediates".Given a surface brightness profile, and assuming circular symmetry (which is satisfied for the inner portions of most galaxies) the total luminosity is found by integration: (2) In order to obtain a luminosity deficit, one must know the un-cored profile of the galaxy.This is not an observable property, so instead an extrapolated profile must be used.The model we adopt is such that the inner slope  for a given profile is equal to the slope at the break radius.We adopt this local extrapolation over a global model to best capture conditions near the core.The extrapolated profile becomes (3) for regions  ∈ (0,  b ].Regions beyond the break radius continue to be described by the full Nuker profile.We then define the luminosity deficit as the difference between the integrated surface brightness profiles The surface brightness values provided are given in mag arcsec −2 , which we convert to  ⊙ pc −2 via the usual: where  ⊙ is the absolute magnitude of the sun in the V-band. In order to produce a mass deficit from the luminosity deficit a mass-to-light ratio, Υ, for each galaxy is required.We use the Υ-color relation developed by Bell et al. (2003): such that   = −0.628and   = 1.305.For optical Υ the scatter is ∼ 0.1 dex.Calculated values of Υ  are also used to determine the stellar mass of each galaxy from their V-band luminosities.A total of 149 galaxies in our sample had extinction corrected ( − ) indices available in the HyperLeda 1 database (Makarov et al. 2014), which we used to compute Υ  .We calculated luminosity deficits for each galaxy numerically using equation ( 4) and converted them into stellar mass deficits using the mass-to-light ratio (Δ = Υ  Δ  ). Figure 2 shows the results of the stellar mass calculations plotted against the reported V-band absolute magnitude, and Figure 3 shows the calculated mass deficit versus stellar mass.The bimodality of cores and power-laws is readily apparent, with cores being more massive and having higher mass deficits (Lauer et al. 2007b).
If MBHBs are responsible for core souring, it is expected that the mass in ejected stars should positively correlate with the cusp radius   (or the break radius   ) (Lauer et al. 2007a).The cusp radius is related to the break radius by where  ′ is defined as with  0 , in this case, the HST resolution scale.In principle  0 can be any limiting radius.Figure 4 shows the expected relationship between mass deficit and both the cusp and break radii, demonstrating that the numerical prescription used here to calculate mass deficits is consistent with expectations.The mass deficit values are also reasonable, with the highest values belonging to the brightest cluster galaxies (BCGs), which have Δ ∼ 10 10  ⊙ .Further, our mass deficit values show a trend consistent with those calculated from Core-Sérsic models (e.g., Graham 2004;Ferrarese et al. 2006 Here the values are plotted against the reported absolute V-band magnitude.The sample spans over 2 dex in stellar mass, allowing good probes of any trends in  ★ 2014), or else our values tend to be about an order of magnitude lower (e.g., Kormendy et al. 2009;Kormendy & Ho 2013;Rusli et al. 2013), which Figure 5 demonstrates.Disparate values are to be expected due to the sensitivity of observed mass deficits to model choices, including the adopted mass-to-light ratios, surface brightness profile, and how one extrapolates the outer profile (Hopkins & Hernquist 2010;Dosopoulou et al. 2021).Kormendy & Ho (2013), for example, use mass-to-light ratios Υ  ∝  0.36   which end up being a factor of ∼ 3× higher than our values.Additionally, Kormendy & Ho (2013) use a global Sérsic fit and extrapolate inwards, resulting in a sensitivity to larger-scale deviations.In this work we extrapolate using a local Nuker fit, making our model more sensitive to local values.
Of the 149 ETG galaxies for which we have mass deficits, 60 are cores.These are plotted in Figure 6 and colored by morphology.The BCGs, colored in purple, seem to fall below the trend suggested by the elliptical and lenticular galaxies.This could be explained by the fact that BCGs are not statistical, being heavily dominated by past mergers (Lauer et al. 2007a), for which the linear relationship between number of dry mergers and stellar mass deficit breaks (Merritt 2006).Of note, we remove the Sa (non-ETG) galaxy NGC 7213 and NGC 7785 which appears to be an outlier, resulting in 58 core galaxies used in our final analysis.

METHODS: ILLUSTRIS TNG300-3
In this Paper we use IllustrisTNG as an avenue for constraining the minimum mass ratio of MBHBs capable of scouring a core, as well as the maximum gas mass fraction, above which a core will be erased.The IllustrisTNG project is a cosmological magnetohydrodynamical simulation suite for galaxy formation, spanning a range of volume and resolution.Its stated purpose is to elucidate when and how galaxies evolve into the structures observed at  = 0.A major component Δ is calculated as outlined in section 2 and shows a positive correlation with  ★ .The bimodality of cores and power-laws can be seen in the data, with cores tending to be more massive with higher mass deficits.
of the simulation is the SubLink algorithm (Rodriguez-Gomez et al. 2015) used to construct merger trees at the subhalo level, where "subhalo" is interchangeable with "galaxy".The merger trees for each galaxy can be traced back, where information about each progenitor can be obtained, including the mass ratio of the subhalos involved in a merger event as well as their gas mass.We use the results of the publicly available TNG300-3 simulation volume and resolution, as described in Nelson et al. (2019), to make predictions about how MBHB properties influence the creation of galaxy cores.We selected the 300 Mpc volume since it gives the best properties and statistics for the massive galaxies of interest, including the stellar mass (which informs SMBH mass) and merger history of each galaxy.Since mass assembly for galaxies with stellar masses  ★ ≳ 10 11  ⊙ is dominated by mergers, we expect core scouring to be most prominent in this mass regime.We take all 1898 subhalos above this lower limit as our simulation sample.IllustrisTNG does not possess the resolution scale required to resolve core structure.The softening length in TNG300-3 for collisionless particles (dark matter and stars) is  DM,★ ( = 0) = 5.90 kpc (Weinberger et al. 2017;Pillepich et al. 2018), while the range of break radii in our core galaxy sample is 11.8 pc <   < 1870 pc.Accordingly, we perform sub-grid post-processing prescriptions to calculate stellar mass deficits.We use the results of N-body simulations presented by Merritt (2006) to calculate the expected mass deficit for each subhalo.In these simulations a primary SMBH is initially located at the center of a galaxy, and a secondary SMBH is allowed to spiral in.Early in the simulation the orbit of the secondary decays due to DF and at later times the system evolves via both DF and stellar scattering.During the latter phase, the mass in ejected stars is tracked and constitutes the mass deficit.Five binary mass ratios were considered, with values  = (0.5, 0.25, 0.1, 0.05, and 0.025).For each value of  three different density profiles were explored, where the inner density slopes were  =  log / log  = (0.5, 1.0, and 1.5).These models were integrated with two different values of : 1.2 × 10 5 and 2.0 × 10 5 particles.It was found that the mass in ejected stars follows the linear relation where  is a scaling factor, N is the number of mergers contributing to core scouring and  BH is the mass of the remnant black hole.We use this relation along with merger histories provided by the TNG300 merger trees to calculate the mass deficit.
Observations have shown a range of values 0.5 ≲ Δ/ BH ≲ 10 (e.g., Ferrarese et al. 2006;Kormendy & Ho 2013;Rusli et al. 2013;Dullo & Graham 2014).Thus the product  N should also fall in this range.Theoretical estimates of  range between 0.5 for spherical systems (Merritt 2006) and 5 for triaxial systems, though the ejected mass in triaxial systems does not solely come from the innermost regions (Khan et al. 2012), a point we return to in section 5. We determine the number of past mergers N a subhalo has had by following the main (primary) progenitor branch back in time as well as following the branches of the next (secondary) progenitors to obtain the full progenitor history.This is done starting at snapshot 100 (corresponding to  ≈ 0 and continuing back to snapshot 33 (mapping to  ≈ 2).Mergers are counted only when a minimum mass ratio  and a maximum gas mass fraction  gas are satisfied.The values are set by our fitting procedure described in section 4. Establishing a mass ratio cut-off reflects the minimum MBHB mass ratio required to successfully scour out a core, while a gas mass fraction constraint models the condition for core erasure.The mass ratio is defied to be where  ★,2 is the stellar mass of the secondary subhalo within the half mass radius, and similarly for the primary subhalo mass  ★,1 .
The gas mass fraction is where  ,1 and  ,2 are the gas mass within the half-mass radius of the primary and secondary respectively.A merger resulting in a remnant which exceeds a provided gas mass fraction is considered a gas-rich (wet) merger while values less than the limit count as gaspoor (dry) mergers.We assume that once a wet merger has occurred any core formation will be erased by subsequent star formation.Thus only the most recent, consecutive dry mergers count toward calculating mass deficits.Since N is determined by  and  gas ,  is a free parameter.The mass of the merged black hole is determined using the  BH - bulge relation (Kormendy & Ho 2013) ;  = 0.29 dex, (12) assuming for our mass range all subhalos will be bulge dominated.When calculating the mass deficit using equation 9 we introduce the scatter of 0.29 dex to determine  BH for each subhalo.Since equation 12 shows a power-law index of about unity,  is approximately the BH mass ratio.
It is expected that a low mass ratio threshold will result in a larger population of core galaxies.Physically, this is a statement about the efficiency of core scouring.Setting a low  means that smaller secondary black holes do make it to the bottom of the potential well and that the binary contributes to stellar scattering.Similarly, a high gas mass fraction threshold for what counts as a gas-poor  (Lauer et al. 2007a).The figure shows the calculated values of the mass deficit to be consistent with this relationship, with Δ and   showing less scatter.This is expected, as our prescription for mass deficit depends on   merger increases the chances a core will not be erased, according to our prescription.Physically, this means that the gas funneled to the nucleus does not result in star formation.If stellar scattering is more efficient and star formation less efficient, then the average mass deficit across the population of core galaxies will increase.At the same time, inefficient scattering and efficient nuclear star formation must result in a lower mass deficit in the core population.Figures 7 &  8 show different realizations of calculated mass deficits for TNG300 subhalos using combinations of  , , and  gas to illustrate the above statements.

COMPARISON OF SIMULATION TO OBSERVATIONS
Here we compare the parameterized model of cores created by MB-HBs with the observed data to determine which mass ratio binaries are capable of core scouring and what minimum gas mass fraction leads to core erasure.Once a synthetic data set is created using our model we could directly compare it to a volume-limited sample of core galaxies, but our sample is heterogeneous.Thus we require an extra step to calculate the expected mass deficit given a galaxy stellar mass.
The observed data are compared to the TNG300 subhalos using the emcee Python MCMC package (Foreman-Mackey et al. 2013).The likelihood function is calculated for a given data point, , with stellar mass and mass deficit values  ★, , Δ  as Here, ì  = [,  gas ] are our parameters of interest.The function  ( ★, , Δ  ) is a kernel density estimation (KDE) evaluated at the data point.The KDE is calculated with standard routines, taking as an argument a synthetic sample of TNG300 subhalos generated as described in section 3 using the specified combination of  and  gas .
To ensure a smooth KDE a scatter of 0.2 dex is introduced to the stellar mass before the mass deficit is calculated.This method is iterated until the synthetic sample has ≥ 10000 data points.Combinations of  and  gas which result in fewer than 50 core galaxies before iteration are ignored as they are not capable of producing the observed prevalence of core galaxies.The function  is normalized by its integral over all Δ in order to penalize models that predict values of Δ  that are too high or low for a given value of  ★, .This avoids problems with an unnormalized , which would penalize models that do not predict distributions that are as clustered as the heterogeneous sample.Our observational sample is not a volume-limited sample, and selection effects would severely bias direct comparison in the two-dimensional space without instead calculating the expected distribution of Δ for a fixed value of  ★ , as we do here.
Figure 9 shows 100 × 100 grid calculations of the log-likelihood as a function of  ∈ [0.01, 1] and  gas ∈ [0.01, 1] for a range of values  = 0.5, 1, 2, 3, normalized to the maximum value among all calculations.The case where  = 0.5 has the highest log-likelihoods, obvious despite the stochastic nature of the calculations.We therefore take this this to be the constant value of  when conducting MCMC fitting.
The MCMC algorithm is provided the log-likelihood (ln L =  ln L  ) and explores the parameter space using flat priors on both the minimum mass ratio and maximum gas mas fraction.The bounds on these parameters are taken to be [0.01,1.00].Figure 10  Comparison of mass deficits for our sample galaxies which overlap with those observed in other studies.Mass deficits calculated in this work are marked as blue circles while those of Graham (2004) are lavender squares, Ferrarese et al. (2006) are green diamonds, Kormendy & Ho (2013) are yellow triangles, Rusli et al. (2013) are blue-green hexagons, and Dullo & Graham (2014) are purple inverted triangles.Mass deficits of other works were taken from tables or else calculated using methods described in the respective papers.Our Δ- ★ trend is consistent with previous literature, except in the cases of Kormendy & Ho (2013) and Rusli et al. (2013) where our values tend to be an order of magnitude lower.Disparate values are due to differing model choices, including adopted mass-to-light ratios, surface brightness profiles (Core-Sérsic, Sérsic, or Nuker), and outer profile extrapolations.
terior samples and 68% confidence intervals are  = 0.7 +0.2 −0.3 and  gas = 0.6 +0.3  −0.2 .The distribution of  values is unimodal and sharply peaked with a rapid fall-off towards small values.Mass ratios of  < 0.25 are ruled out at the 95% confidence level.The maximum gas mass fraction distribution is broader, but shows a clear peak near 0.45. Figure 11 shows a simulated sample with calculated mass deficits using the median values of the posterior distributions of  and  gas .The observed sample is over-plotted for visual comparison.We find the predicted trend to be in good agreement with observation.

DISCUSSION
The results of the fitting procedure show that the minimum mass ratio for core formation is well constrained in our model.Since  BH ∼   ★ with  approximately unity (see equation 12), this is a direct constraint on the MBHB mass ratio.Low mass ratio binaries, therefore, do not contribute to core scouring.Further, up to around half the mass within the half-mass radius of a remnant ETG can consist of gas without the core being erased.Gas content much larger than this would predict larger mass deficits than we observe (as shown in Figure 7).
Throughout the analysis we have made a number of practical assumptions to make comparisons between observational data and synthetic data.In particular, we have assumed that (1) the extrapolated surface brightness profile of each galaxy is a simple power-law ( 2 .Mass deficit versus stellar mass for core galaxies with  ★ ≳ 10 11  ⊙ colored by morphological type.The BGC galaxies appear to fall below the trend suggested by the elliptical and lenticular galaxies.NGC 7213 is an Sa galaxy, while NGC 7785 appears to be an outlier.Both galaxies are removed from the sample when compared with TNG300 data.a single gas-rich merger is capable of erasing all core structure in a galaxy, (3) all MBHBs coalesce before a subsequent merger (i.e., there are never any triple BH systems), (4) the eccentricity of the MBHB orbit is irrelevant, (5) the mass ejected in stars is always a linear function of N , (6) the  BH −  bulge relation does not evolve with redshift, and (7) DF and stellar scattering timescales are irrelevant.We consider below the effects of the above assumptions on our results.
For any single galaxy, missing-light measurements are sensitive to the selected model (Hopkins & Hernquist 2010;Dosopoulou et al. 2021), though measurements over a population of galaxies should be largely accurate, making our assumption of a power-law, un-cored profile a fair approximation.As mentioned above, the total mass ejected is a simple product of N , the number of mergers, and   BH , an amount of mass scattered per merger that scales with mass of the black hole.For most of our calculations, we assume  = 0.5, but there are simulation-based estimates for  as large as 5 (Khan et al. 2012).The value for N comes from the TNG300-3 simulations, whose merger tree is consistent with observed pair-fractions and merger rates (Bell et al. 2004;Man et al. 2012;Rodriguez-Gomez et al. 2015).We also find the best agreement between theory and observations for  = 0.5, though our preferred values of  and  gas do not change (Figure 9).
It is worth understanding why the more realistic triaxial simulations due to Khan et al. (2012) result in values of  > 1 whereas the spherical simulations due to Merritt (2006) predict  = 0.5 in best agreement with our data.The triaxial systems include orbits that come from large radii to refill the loss cone, thus enabling the binary to merge within a Hubble time (Yu 2002;Holley-Bockelmann & Sigurdsson 2006;Khan et al. 2012;Holley-Bockelmann & Khan 2015;Mannerkoski et al. 2023).Because stars on these orbits are responsible for the removal of a substantial amount of energy from the binary, they constitute a significant amount of the total mass deficit. .The grid of plots shows mass deficit versus stellar mass for both the observed sample (large blue circles) and synthetic data sets (small purple circles) generated from combinations  and  gas , determining N, and the factor  = 0.5 in equation 9.The top row holds the maximum gas mass fraction required to produce a core fixed while the minimum mass ratio is allowed to vary, while in the bottom row the reverse is true.A low mass ratio threshold results in more core galaxies compared to high thresholds when keeping  gas fixed.A high gas mass fraction threshold, for fixed , results in more core galaxies.Additionally, for low  and high  gas , mass deficits are higher on average than for other combinations.
These stars, however, likely do not contribute a substantial amount of light at the center of the galaxy.Thus while the total amount of mass ejected is larger (as is necessary to solve the final parsec problem), the amount of observable mass deficit is much smaller and likely to be well approximated by the amount of mass ejected in a spherical galaxy, since this will only remove stars from the central part of the galaxy where light deficits are observed.
Assuming a single wet merger can completely erase core structure will tend to underestimate the final size of the core, and therefore the stellar mass deficit.It may be an oversimplification that all gas-rich mergers result in full core erasure, since not all gas will be cold, star-forming gas and the central SMBH of a remnant can act as a guardian, protecting the core (Faber et al. 1997).
N-body simulations would be needed to understand the effect of a triple BH system on stellar scattering, including the possibility of the ejection of the primary and secondary from the system (e.g., Hoffman & Loeb 2007;Gualandris & Merritt 2008;Mannerkoski et al. 2023).A tertiary black hole can possibly decrease the mass deficit by stalling stellar scattering, extracting energy from the primary and secondary, hardening the orbit and "skipping" a large portion of the scattering phase (e.g., Blaes et al. 2002).A highly eccentric MBHB can also decrease the mass deficit, as binaries with large eccentricities coalesce on a shorter timescale (Peters & Mathews 1963;Peters 1964).Thus the size of cores also depends on the distribution of MBHB eccentricities.
Assuming a linear relation between number of mergers and mass deficit may be an oversimplification, as shown by Merritt (2006).In N-body experiments the the mass deficit as a function of N is no longer linear after ≳ 3 dry merger events.However, Figure 12 shows there are no galaxies in the synthetic sample whose number of recent, consecutive gas-poor mergers exceeds three.This is consistent with the results described in Merritt (2006) and with hierarchical structure formation (e.g., Bell et al. 2004;Man et al. 2012).
If the  BH - bulge relation has a larger normalization at higher redshifts, then we might expect larger cores.MBHB binaries formed at higher redshifts would have larger masses compared to their host galaxy's total stellar mass, leading to more stellar scattering.Such an evolution would imply larger values of  and/or smaller values of  gas .
The minimum mass ratio of MBHBs capable of carving out a core and the maximum gas content are affected by the assumption DF and stellar scattering timescales can be neglected.We have not ) q = 0.75, f gas = 0.75 considered in our treatment how the effective radius, central velocity dispersion, or the nuclear stellar density of the progenitor galaxies contribute to whether the secondary BH forms a binary with the primary.In future work these physical considerations can be added to the post-processing phase when determining predicted mass deficits and thus lead to further constraints on MBHB properties.An additional constraining factor is the number of core galaxies observed within some volume.The current sample is not volume-limited and is therefore not a representative population of core galaxies, though a statistical comparison between the number of core galaxies observed and those predicted by IllustrisTNG would provide an important constraint.This can be achieved in the future by creating a volumelimited sample using archival HST data or new observations.Even when making these assumptions, we are able to see a clear picture emerge in which (i) we can rule out low-mass mergers as relevant for core scouring and (ii) we find gas-rich mergers must be capable of erasing cores or else we would find larger mass deficits than are seen in the data.

SUMMARY
In this Paper we presented a general method for computing stellar mass deficits, which involved parameterization of the light distribution of a galaxy in question and the adoption of a model for the original central light profile.Once we derived a surface brightness profile, and selected a model, the light deficit was computed by taking the integrated difference between them (see Figure 1).The missing light was then translated into a stellar mass deficit by the multiplication of a mass-to-light ratio.We followed this prescription for galaxies with fitted Nuker profiles in section 2.   Mass deficit versus stellar mass for for the set of TNG300-3 subhalos designated as core galaxies when the minimum mass ratio for core scouring and the maximum gas mass fraction are set to the median values obtained from MCMC fitting.Each subhalo is colored according to the number of most recent, consecutive gas-poor mergers it has had.There are no subhalos with N > 3 in this simulated sample, which is consistent with the results described in Merritt (2006) and with hierarchical structure formation (e.g., Bell et al. 2004;Man et al. 2012).
Although the model of core scouring has been well developed, the connection between the properties of core galaxies and the massive black hole binary population is poorly constrained.A major reason for this is that the merger history of any individual galaxy is largely unknown.We made progress in solving this problem by using TNG300-3 to obtain subhalo stellar masses and and merger histories via SubLink merger trees.We applied sub-grid, post-processing physics (described in section 3) to create synthetic Δ- ★ trends to compare with the observed trend using MCMC methods (section 4).
We found the median values of the posterior distribution for our parameters to be  = 0.7 +0.2 −0.3 and  gas = 0.6 +0.3 −0.2 .Low mass ratio binaries do not contribute to core formation, and gas-rich mergers must be capable of erasing cores.These constraints add to our understanding of the distribution of MBHB mass ratios and the expected rates of these binaries.Further, the distribution of mass ratios for MBHBs contributing to core scouring is consistent with constraints on the MBHB determined using NANOGrav's 15-year PTA data set (Agazie et al. 2023b).These results and future endeavours to increase understanding of the MBHB demographic will affect how the origins of the GWB are approached along with future science done by PTAs like NANOGrav and in the course of the LISA mission (Ellis et al. 2023;Amaro-Seoane et al. 2023).
Figure1.Stellar mass deficits in core galaxies may be calculated if the surface brightness profile is known and a model of the extrapolated light distribution is created.The light deficit is the integrated difference between the two profiles, which can be converted into a stellar mass deficit by multiplication of a stellar mass-to-light ratio in the appropriate band.The "break radius" indicates the distance from the galaxy center where the profile transitions from the shallow inner region to the steep outer region, representing the size of the core.Based on the profile of NGC 1016.

Figure 2 .
Figure 2. The stellar mass of each galaxy was calculated with equation (6).Here the values are plotted against the reported absolute V-band magnitude.The sample spans over 2 dex in stellar mass, allowing good probes of any trends in  ★

Figure 3 .
Figure3.Relationship between the stellar mass deficit and stellar mass.Δ is calculated as outlined in section 2 and shows a positive correlation with  ★ .The bimodality of cores and power-laws can be seen in the data, with cores tending to be more massive with higher mass deficits.

Figure 4 .
Figure4.The stellar mass deficit Δ is expected to be a function of the cusp radius   or the break radius  (Lauer et al. 2007a).The figure shows the calculated values of the mass deficit to be consistent with this relationship, with Δ and   showing less scatter.This is expected, as our prescription for mass deficit depends on Figure5.Comparison of mass deficits for our sample galaxies which overlap with those observed in other studies.Mass deficits calculated in this work are marked as blue circles while those ofGraham (2004) are lavender squares,Ferrarese et al. (2006) are green diamonds,Kormendy & Ho (2013) are yellow triangles,Rusli et al. (2013) are blue-green hexagons, andDullo & Graham (2014) are purple inverted triangles.Mass deficits of other works were taken from tables or else calculated using methods described in the respective papers.Our Δ- ★ trend is consistent with previous literature, except in the cases ofKormendy & Ho (2013) andRusli et al. (2013) where our values tend to be an order of magnitude lower.Disparate values are due to differing model choices, including adopted mass-to-light ratios, surface brightness profiles (Core-Sérsic, Sérsic, or Nuker), and outer profile extrapolations.
Figure6.Mass deficit versus stellar mass for core galaxies with  ★ ≳ 10 11  ⊙ colored by morphological type.The BGC galaxies appear to fall below the trend suggested by the elliptical and lenticular galaxies.NGC 7213 is an Sa galaxy, while NGC 7785 appears to be an outlier.Both galaxies are removed from the sample when compared with TNG300 data.
Figure7.The grid of plots shows mass deficit versus stellar mass for both the observed sample (large blue circles) and synthetic data sets (small purple circles) generated from combinations  and  gas , determining N, and the factor  = 0.5 in equation 9.The top row holds the maximum gas mass fraction required to produce a core fixed while the minimum mass ratio is allowed to vary, while in the bottom row the reverse is true.A low mass ratio threshold results in more core galaxies compared to high thresholds when keeping  gas fixed.A high gas mass fraction threshold, for fixed , results in more core galaxies.Additionally, for low  and high  gas , mass deficits are higher on average than for other combinations.

Figure 8 .
Figure8.Mass deficit versus stellar mass for both the observed sample (large blue circles) and synthetic data sets generated from  =  gas = 0.75, determining the value of N, and the factor  in equation 9 Synthetic data sets are marked according to the value of  chosen, with  = 0.5 corresponding to purple circles,  = 1 to green diamonds,  = 2 to yellow triangles, and  = 3 to light blue squares.Increasing  can be seen to simply increase the average mass deficit at fixed  ★ .

Figure 9 .
Figure9.Log-likelihood calculations with a 100 × 100 grid of  and  gas values with a range of [0.01, 1].The four panels represent different calculations using  = 0.5, 1, 2, 3, starting from the top panel.The log-likelihood is normalized to the maximum value found among all calculations.The scenario where  = 0.5 has the highest log-likelihoods, despite the stochastic nature of the calculations.

Figure 10 .
Figure10.Corner plot for the best fit minimum mass ratio and maximum gas mass fraction for the observed data.The histogram for  is unimodal, with a sharp peak and rapid fall-off towards low mass ratios.The model rules out mass ratios  < 0.25 at 95% confidence.On the other hand, the histogram for  gas shows a peak near 0.45, though the distribution is broad.

Figure 11 .
Figure11.A comparison between the observed mass deficit as a function of stellar mass and a sample synthesised from TNG300-3 using the median values for the  and  gas posterior samples.The main trend in the observed data is recovered.Note that the absolute number of core galaxies is not a constraint, as our sample is heterogeneous.It is only the main trend and overall values of Δ for a given  ★ we are able to reproduce. ,