The (Black Hole Mass)-(Spheroid Stellar Density) Relations: $M_{\rm BH}$--$\mu$ (and $M_{\rm BH}$--$\Sigma$) and $M_{\rm BH}$--$\rho$

This paper is the fourth in a series presenting (galaxy morphology, and thus galaxy formation)-dependent black hole mass, $M_{\rm BH}$, scaling relations. We have used a sample of 119 galaxies with directly-measured $M_{\rm BH}$ and host spheroid parameters obtained from multi-component decomposition of, primarily, $3.6\,\mu$m Spitzer images. Here, we investigate the correlations between $M_{\rm BH}$ and the projected luminosity density $\mu$, the projected stellar mass density $\Sigma$, and the deprojected (internal) stellar mass density $\rho$, for various spheroid radii. We discover the predicted $M_{\rm BH}$--$\mu_{\rm 0,sph}$ relation and present the first $M_{\rm BH}$--$\mu_{\rm e, sph}$ and $M_{\rm BH}$--$\rho_{\rm e,int, sph}$ diagrams displaying slightly different (possibly curved) trends for early- and late-type galaxies (ETGs and LTGs) and an offset between ETGs with (fast-rotators, ES/S0) and without (slow-rotators, E) a disk. The scatter about various $M_{\rm BH}$--$\langle\Sigma\rangle_{\rm R,sph}$ (and $\langle\rho\rangle_{\rm r,sph}$) relations is shown to systematically decrease as the enclosing aperture (and volume) increases, dropping from 0.69~dex when using the spheroid \enquote{compactness}, $\langle\Sigma\rangle_{\rm 1kpc,sph}$, to 0.59~dex when using $\langle\Sigma\rangle_{\rm 5kpc,sph}$. We also reveal that $M_{\rm BH}$ correlates with the internal density, $\rho_{\rm soi,sph}$, at the BH's sphere-of-influence radius, such that core-S\'ersic (high S\'ersic index, $n$) and (low-$n$) S\'ersic galaxies define different relations with total rms scatters 0.21~dex and 0.77~dex, respectively.(Abridged)


INTRODUCTION
The number of galaxies with directly-measured black hole masses, i.e., where observations could resolve the black hole's gravitational sphere-of-influence, has grown to about 145 galaxies (Sahu et al. 2019a). Using state-ofthe-art two-dimensional modeling (Ciambur 2015) and multi-component decompositions (Ciambur 2016), we have modeled the surface brightness profiles of 123 of these galaxies 1 and their components. We have discovered morphology-dependent correlations between the black hole mass (M BH ) and various host galaxy properties, such as the galaxy stellar mass (M * ,gal ), the spheroid stellar mass (M * ,sph ), the spheroid central light concentration or Sérsic index (n sph ), the spheroid effective half-light radius (R e,sph ), and the central stellar velocity dispersion (Graham 2012;Scott et al. 2013;Savorgnan et al. 2016;Davis et al. 2018;Davis et al. 2019;Sahu et al. 2019bSahu et al. ,a, 2020. These have improved as the quality, and the quantity of data has grown. The simple (galaxy morphology)independent black hole scaling relations 2 (e.g., Dressler & Richstone 1988;Magorrian et al. 1998;Häring & Rix 2004;Gültekin et al. 2009;Kormendy & Ho 2013;Mc-Connell & Ma 2013) are, in fact, too simple to accurately trace the coevolution of the different types of galaxies and their black holes. For example, the small and massive bulges of spiral and lenticular galaxies follow M BH -M * ,sph relations different from that of elliptical galaxies (Sahu et al. 2019b). It is hoped that the advances with morphology-dependent correlations will help identify which correlation is more fundamental, i.e., primary versus secondary. However, one of the potential candidates remains to be explored; it involves stellar density.
Most of the morphology-dependent black hole scaling relations are significantly different to the familiar but now superseded "single relations" obtained when all galaxy types are combined. Crucially, diagrams with differing numbers of different galaxy types can yield "single relations" with different slopes and intercepts. As a result, many of the past "single relations" (built by grouping galaxies of different morphological types) are physically meaningless, merely representative of the ratio of the different galaxy types in that sample.
The above realization is fundamental if we are to adequately understand the co-evolution of galaxies and their central massive black holes. This is because the black hole mass is, in a sense, aware of the different formation history and physics which went into building its host galaxy. What is important is not simply the amount of mass in stars and perhaps dark matter, but how that mass was assembled (and moves) to create a galaxy's substructure/morphology.
Using a sample of 27 galaxies, Graham & Driver (2007) observed a strong correlation between M BH and the bulge central concentration, which is quantified by the shape parameter of spheroid's surface brightness profile the Sérsic index (n sph , Trujillo et al. 2001). Graham & Driver (2007) found a comparable level of (intrinsic) scatter about the M BH -n sph relation as seen in the M BH -(central stellar velocity dispersion: σ) relations observed at that time, which was about 0.3 dex. The observed stellar velocity dispersion traces the underlying mass distribution and radial concentration of light ). Thus, Graham & Driver (2007) suggested that a combination of the central stellar density and the central light concentration of a spheroid may be governing the black hole-host spheroid connection.
Some studies (e.g., Graham & Guzmán 2003;Merritt 2006) presented a correlation between central concentration and central stellar density, suggesting that one of these quantities can be written in terms of the other; although, it is still not clear which quantity is more fundamental. Graham & Driver (2007) combined this relation with their linear and curved 3 M BH -n sph relations, to predict both a linear and a curved M BH -(spheroid central surface brightness or central projected density, µ 0,sph ) relation (Graham & Driver 2007, their equations 9 and 10). Moreover, they suggested that an even better correlation might exist between M BH and the (threedimensional) deprojected density (ρ, aka the internal or spatial density) at the center of the spheroid. For the first time, here we explore these predicted M BH -(stellar density) relations and others.
We present new correlations between M BH and the spheroid surface brightness (projected/column luminosity density), the projected (or column) stellar mass density (Σ), and the deprojected stellar density at various spheroid radii. Our sample of 123 galaxies is described in the following Section 2. That section also describes the linear regression applied and the parameter uncertainty used. The calculation of the deprojected density is detailed in the Appendix A, where we also compare our numerically calculated internal density with an approximation from the model of Prugniel & Simien (1997).
Section 3 presents the correlation between M BH and the spheroid projected (luminosity and stellar mass) density at various radii (center, 1 kpc, 5 kpc, and halflight radius). In Section 4, we reveal additional new correlations obtained between M BH and the bulge internal mass density at various (inner and larger) radii, including the sphere-of-influence radius of the black hole. We also provide the projected and deprojected density profiles, µ(R) and ρ(r), to help explain various trends obtained between M BH and µ, and between M BH and ρ at different spheroid radii.
In all of these diagrams, we also investigate possible dependence on galaxy morphology, e.g., early-type galaxies (ETGs: elliptical E, ellicular ES 4 , and lenticular S0) versus late-type galaxies (LTGs: spirals S), centrally-fast (ES, S0, S) versus slow (E) rotators, core-Sérsic 5 versus Sérsic 6 galaxies, and barred versus nonbarred galaxies. We compare our findings with the morphology-dependent substructures seen in our recently published correlations (Sahu et al. 2019b(Sahu et al. ,a, 2020. In Section 5 we discuss our results and some of the more notable implications. Finally, we summarize the main results of this work in Section 6. We have used the terms spatial density, and internal density interchangeably for the (3D) deprojected density throughout this paper. All uncertainties are quoted at the ±1 σ (≈ 68%) confidence interval.

DATA
We have used the spheroid's structural parameters from Savorgnan & Graham (2016b), Davis et al. (2019), and Sahu et al. (2019b), which were obtained from multicomponent decompositions of 123 galaxies with directlymeasured central black hole masses reported in the literature. The direct methods for black hole mass measurement include stellar dynamical modeling, gas dynamical modeling, megamaser kinematics, proper motions (for Sgr A * ), and the latest direct imaging (for M87*). The majority of the galaxy images (81%) were in the 3.6 µm-band taken by the infrared array camera (IRAC, Fazio et al. 2004, resolution ∼ 2 ) onboard the Spitzer Space Telescope. The remaining images came from the archives of the Hubble Space Telescope (HST, 11%), the Sloan Digital Sky Survey (SDSS, 3%), and the Two Micron All Sky Survey (2MASS, 6%). For full details of the image analysis, we refer readers to the aforementioned three studies.
Briefly, we performed 2D modeling of the galaxy images using our in-house 7 software isofit and cmodel (Ciambur 2015), which were built into the image reduction and analysis facility (iraf, Tody 1986Tody , 1993. isofit fits quasi-elliptical isophotes at each galactic radii. It uses an elliptical coordinate system, thereby improves upon the spherical coordinate system imple-5 The core-Sérsic galaxies are generally the most massive galaxies, likely formed through major gas-poor mergers. The eventual coalescence of their central massive black holes scours out the stars from the central "loss cone" (through the transfer of the binary black hole's orbital angular momenta) and creates a deficit of light at the center, referred to as a "core". The bulge surface brightness profile for a core-Sérsic galaxy is described by a core-Sérsic function , which consists of a shallow inner power-law followed by a Sérsic function (Sérsic 1968(Sérsic , 1963 at larger radii. Such cores were first noted by King & Minkowski (1966). 6 Sérsic galaxies do not have a deficit of light at their center. 7 The software isofit, cmodel, and profiler are publicly available at the GitHub platform (see Ciambur 2015Ciambur , 2016 mented in ellipse (Jedrzejewski 1987a,b). The angular coordinate known as the "eccentric anomaly" is used for uniform sampling of the quasi-elliptical isophotes, and the code employs Fourier harmonics to capture the isophotal deviations from a pure ellipse (Carter 1978;Kent 1984;Michard & Simien 1988). Thus, isofit generates an (azimuthally-averaged) one-dimensional surface brightness profile along any galaxy axis, together with the radial variations of the isophotal ellipticity (e), position angle, and Fourier coefficients. These parameters are used to create a 2D galaxy model via cmodel.
The model captures all symmetric features about the major-axis (mirror symmetry) and leaves behind disturbances and star clusters which can be explored in the "residual image".
We disassemble the galaxy model into components with the help of various functions inbuilt in the software profiler (Ciambur 2016). A galaxy can have a bulge, intermediate-or large-scale disk, bar, ansae, rings, depleted core, and nuclear components (e.g., star cluster, nuclear bar, disk, or ring). The presence of disks and bars in our decompositions were verified, whenever possible, through recourse to the literature, including kinematic evidence for disk rotation. We perform this multicomponent decomposition using the surface brightness profile along the galaxy's major-axis as well as the socalled "equivalent-axis", which represents a radial-axis equivalent to a circularised form of the galaxy's quasielliptical isophotes, such that the total enclosed luminosity remains conserved 8 .
The multi-component decomposition process provides us with the surface brightness profiles of individual galaxy components and the detailed galaxy morphology, indicating the presence of a rotating disk, a depleted core, a bar, etc. One of the most noted of all galactic components is the bulge, whose surface brightness distribution is described using the Sérsic (1963) function (Appendix Equation A1), which is parameterized by the Sérsic index (n sph ), effective half-light radius (R e,sph ), and the surface brightness at the halflight radius (µ e,sph = −2.5 log I e,sph ). The equivalentaxis spheroid surface brightness profiles for our 3.6 µm (Spitzer) sample are shown in Figure 1.
To obtain the spheroid's internal (deprojected) stellar mass density distribution, ρ(r), we performed an inverse 8 The equivalent-axis radii, Req, for an isophote is the geometricmean of the isophote's major-and minor-axis radii (R maj and R min , respectively), i.e., Req = R maj × R min or Req = R maj 1 − e maj (see Ciambur 2015Ciambur , 2016, for more details on the isophotal galaxy modeling, multi-component decomposition, and the circularised equivalent-axis).
Abel transformation (Abel 1826) of the (circularly symmetric) equivalent-axis spheroid surface brightness profiles. The numerical calculation of the internal density profiles and a comparison with the Prugniel & Simien (1997) density model (an approximation to the exact deprojection of the Sérsic profile) is presented in the Appendix Section A.
The spheroid parameters required to calculate the projected and the internal stellar mass densities, e.g., the bulge surface brightness parameters (n sph , R e,sph , µ e,sph ), along with the galaxy morphology, distances, physical (arcsec-to-kpc) scale 9 , stellar massto-light ratio, and the image band information for all 123 galaxies are available in Sahu et al. (2020, their Appendix Table A1). Sahu et al. (2020) also tabulates the directly-measured central black hole masses and the bulge stellar masses (M * ,sph ) of these 123 galaxies.
Here, we excluded the galaxies NGC 404, NGC 4342, NGC 4486B, and the Milky Way throughout our investigation, unless expressly stated otherwise. NGC 404 is the only galaxy with a black hole mass below 10 6 M (Nguyen et al. 2017) and it may skew/bias the results. Its published black hole mass has a sphere of influence five times smaller than the seeing (∼ 0.1 ) under which it was measured. NGC 4342 and NGC 4486B have been heavily stripped of their mass due to the gravitational pull of their massive companion galaxies (see Batcheldor et al. 2010;Blom et al. 2014). For the Milky Way, the available surface brightness profile (Kent et al. 1991;Graham & Driver 2007) was not flux-calibrated to obtain a calibrated density profile. The exclusion of these galaxies leaves us with a reduced sample of 119 galaxies. All the galaxies excluded from the linear regressions (performed to obtain the scaling relations presented here) are shown with a different symbol in the ensuing diagrams.
We use the bivariate correlated errors and intrinsic scatter (bces) regression (Akritas & Bershady 1996) to obtain our black hole scaling relations. bces is a modification of the ordinary least squares regression. It considers measurement errors in both variables (and their possible correlation) and allows for intrinsic scatter in the distribution. We prefer to use the bces(bisector) 10 line obtained by symmetrically bisecting the bces(Y |X) line (which minimizes the error-weighted root mean square, rms, vertical offsets about the fitted line) and the bces(X|Y ) line (which minimizes the error-weighted rms horizontal offsets about this different fitted line). We do this partly because it is unknown whether the host spheroid density is an independent variable and the central black hole mass is a dependent variable, or viceversa, or if there is an interplay. We also check our best-fit parameters using a symmetric application (Novak et al. 2006) of the intrinsically non-symmetric (modified fitexy) known as mpfitexy 11 regression (Markwardt 2009;Williams et al. 2010).
Uncertainties in M BH , and spheroid profile parameters n sph (±0.09 dex), R e,sph (±0.13 dex), and µ e,sph (±0.58 mag arcsec −2 or ±0.23 dex in µ e,sph /2.5) are taken from Sahu et al. (2020, see their section 2 for more details). The uncertainty in the internal density (ρ e ) at an internal radius equal to the projected half-light radius (R e,sph )-obtained by propagating errors in the spheroid parameters through the analytical expression (Equation A7)-are ∼ ±0.30 dex. For densities at other radii, the error propagation (assuming independent parameters) through the internal density expression (Equation A6) provides even higher uncertainties due to multiple occurrences of n sph and R e,sph , in addition to ρ e . Such uncertainties are likely to be overestimated and can affect the best-fit lines. Therefore, we used a constant uncertainty of ±0.23 dex on the projected mass densities (Σ) and, similarly, a constant uncertainty of ±0.30 dex on the internal densities for all the correlations, unless stated otherwise. Additionally, we test the stability of our correlations (their slopes and intercepts) using a range of (zero to ±0.38 dex) uncertainties for the projected and internal densities. To study the correlation between black hole mass and the host spheroid's central surface brightness (µ 0,sph ), which is dependent on the image wavelength band, we used our 3.6 µm-sample comprised of 97 galaxies from the reduced sample of 119 galaxies (see Section 2). This includes 72 Sérsic galaxies, i.e., galaxies with a Sérsic spheroid surface brightness profile, and 25 core-Sérsic galaxies, i.e., galaxies with a depleted central core whose spheroid profile is described by a shallow central powerlaw followed by a Sérsic function at larger radii (see Graham et al. 2003).
Using the (equivalent-axis) surface brightness parameters (n sph , R e,sph , and µ e,sph ) for the spheroids, we cal-  (Sérsic) surface brightness profiles for our 3.6 µm-sample. Right-hand panel: The horizontal axis is normalized at the (projected) half-light radii of each spheroid. The color sequence blue-white-red traces the increasing black hole mass and helps with the understanding of the positive/negative trends observed between MBH and the spheroid surface brightness (and the projected stellar mass density) presented in Section 3. Figure 2. Left-hand panel: Black hole mass versus the central surface brightness (in the AB magnitude system) of the spheroids using our 3.6 µm sample. Right-hand panel: Black hole mass versus the projected central stellar mass density, including the 22 non-Spitzer galaxies. The dark green line represents the best-fit obtained from the bces(bisector) regression. The dark green shaded region around the best-fit line delineates the ±1σ uncertainty on the slope and intercept, and the light green shaded area outlines the ±1σ rms scatter in the data. The same description follows for all other correlations presented in this paper. For both Sérsic and core-Sérsic galaxies, µ 0,3.6µm,sph and Σ 0,sph have been obtained through the inward extrapolation of the Sérsic portion of their spheroid profiles. Core-Sérsic galaxies have a deficit of light at their core and, hence, their µ 0,3.6µm,sph values depicted here are brighter than the actual value. The galaxies excluded from the regression are marked with a black star. culated µ 0,sph via µ 0 = µ e − 2.5 log e b (Equation A1 at R=0), i.e., an inward extrapolation of the Sérsic fit to the spheroid's surface brightness profile. It is important to note that for our core-Sérsic galaxies, µ 0 has been obtained through the inward extrapolation of the Sérsic part of their spheroid profile (as in the L gal -µ 0 diagram of Jerjen et al. 2000). This is because the size of the depleted core is generally much smaller than the ∼ 2 spatial resolution of IRAC images (see Dullo & Graham 2014), and, as such, the (3.6 µm-band) parameters for the central power-law of our core-Sérsic spheroids are not accurate 12 . Thus, the µ 0 used here for the cored galaxies represents their central surface brightness before the damaging effect of binary black holes, which will cause a departure of cored galaxies from an initial M BH -µ 0 trend line. This is the case with cored ETGs in the M * ,gal -µ 0 diagram shown in Graham & Guzmán (2003, their figure 9), which accounted for the central mass/light deficit in the cored galaxies.
The high-n sph galaxies M 59, NGC 1399 (cored), and NGC 3377 are marked by black stars and have the brightest µ 0,3.6µm,sph (∼ 3 mag arcsec −2 ) in Figure 2. They reside beyond the 2σ scatter of the remaining dataset and have significant leverage on the best-fit line, such that including these three galaxies in the regression changes the slope by 1σ. Therefore, these three galaxies were excluded from the regression (in addition to the four exclusions mentioned in Section 2) to obtain the M BH -µ 0,3.6µm,sph and the M BH -Σ 0,sph relations reported here.
The M BH -µ 0,3.6µm,sph relation 13 plotted in the lefthand panel of Figure 2 was obtained using 94 (Sérsic +core-Sérsic) galaxies with 3.6 µm imaging data, and can be expressed as, log M BH M = (−0.41 ± 0.04) µ 0,3.6µm,sph − 13 mag arcsec 2 + (7.97 ± 0.10). (1) The total (measurement error and intrinsic scatter) rms scatter (∆ rms|BH ) is 1.03 dex in the log(M BH )-direction. This correlation quantifies how the (Sérsic) spheroids hosting more massive black holes have a brighter central surface brightness, qualitatively consistent with the linear prediction of the log(M BH )-µ 0,sph relation in Graham & Driver (2007, their equation 9 based on B-band data). This trend can also be inferred from the spheroid profiles plotted in the left-hand panel of Figure 1, where for R tending to zero, µ becomes brighter when moving from low-M BH (blue profiles) to high-M BH (red profiles). In our M BH -µ 0,3.6µm,sph diagram (Figure 2), the core-Sérsic galaxies are represented with the µ 0,3.6µm,sph value that they presumably would originally have if their cores did not undergo a depletion 14 of light due to co-13 The uncertainty we assigned to µ 0,3.6µm,sph is 0.58 mag arcsec −2 (the same as assigned to µ e,3.6µm,sph ); however, consistent relations are obtained upon using up to 2 mag arcsec −2 uncertainty.
14 The deficit of light at the center of core-Sérsic galaxies is generally only a small fraction (5%, average value from alescing BH binaries in dry major-mergers (Begelman et al. 1980). Hence, one should not use the above relation to estimate M BH using the actual (depleted) central surface brightness (µ 0,core ) for cored galaxies, instead, the µ 0 extrapolated from the Sérsic potion of their spheroid profile can be used. The actual µ 0,core for the core-Sérsic galaxies dims with increasing M * ,gal (Graham & Guzmán 2003).
To include our remaining (non-Spitzer) sample of 22 galaxies, we mapped the central surface brightness, µ 0,sph , values to the central surface stellar mass density (Σ 0,sph ) with the units of solar mass per square parsec (M pc −2 ) using Equation A5. We obtained a positive log-linear M BH -Σ 0,sph relation, which is represented in the right-hand panel of Figure 2. The best-fit relation is provided in Table 1, along with the (intrinsic 15 and total) rms scatter, Pearson correlation coefficient, and Spearman rank-order correlation coefficient. The M BHµ 0,3.6µm,sph and M BH -Σ 0,sph relations obtained using only Sérsic galaxies are consistent with the relations obtained when including the core-Sérsic galaxies.

Projected Mass Density within 1 kpc: The
Spheroid Compactness Σ 1kpc,sph The projected stellar mass density ( Σ 1kpc ) within the inner 1 kpc of a galaxy has been used as a measure of galaxy "compactness" (Barro et al. 2017;Ni et al. 2020), and to identify compact star forming galaxies (Suess et al. 2021). Interestingly, for star-forming galaxies Σ 1kpc has been found to correlate with the black hole growth 16 , and it has been suggested that this correlation is stronger than the connection between the black hole growth and host galaxy stellar mass (see Ni et al. 2020). Additionally, it has been suggested that Σ 1kpc is a better indicator of black hole growth than the projected stellar mass density within the galaxy half-light radius and the projected density within other (smaller or larger) constant (e.g., 0.1 kpc, 10 kpc) radii (Ni et al. 2019(Ni et al. , 2020. Here we investigate a possible correlation between the black hole mass and the average projected stel-15 It should be noted that this depends on the adopted parameter uncertainties. 16 The connection between black hole growth and the host galaxy's stellar mass density is explained by the assumption of a linear correlation between stellar density and gas density for star-forming galaxies (see Lin et al. 2019, and references therein). Thus, a high value of Σ 1kpc for a star-forming galaxy infers a high gas density, and the abundance of gas at the inner galactic regions is known to boost the black hole growth (Dekel et al. 2019;Habouzit et al. 2019).
Central Surface Brightness ( Projected Density within R e,sph ( Figure  a Regression performed after excluding three outliers, see Section 3.1 for more details. lar mass density ( Σ 1kpc,sph ) within the inner 1 kpc of the host spheroid. Thus, we refer to Σ 1kpc,sph as the spheroid compactness. Most of our sample with directly-measured M BH are quiescent, with Σ 1kpc,sph greater than the critical/threshold value ( Σ 1kpc = 3 × 10 3 M pc −2 , Cheung et al. 2012) used to identify the quiescent galaxies (Hopkins et al. 2021). Additionally, we explored how the correlation between M BH and spheroid compactness compares against the correlation between M BH and the spheroid densities at/within other radii. We find a tight M BH -Σ 1kpc,sph correlation (lefthand panel in Figure 3), where all galaxy types (ETGs+LTGs) seem to follow a single positive rela- with ∆ rms|BH = 0.69 dex (see Table 1 for correlation coefficients). Similarly, we also see a positive trend between M BH and the column (stellar mass) density within other projected spheroid radii (e.g., 0.01 kpc, 0.1 kpc, 5 kpc, 10 kpc). This positive trend is evident from the  distribution of spheroid profiles in the left-hand panel of Figure 1, where, at all fixed radii, the redder profiles with higher M BH are brighter than the bluer profiles with lower M BH . The correlation between M BH and densities within fixed physical radii smaller than 1 kpc ( Σ 0.01kpc,sph and Σ 0.1kpc,sph ) are not as tight as the above relation (Equation 2). However, we find better M BH -Σ R,sph correlations for R > 1 kpc, with a gradually shallower slope and smaller scatter than the M BH -Σ 1kpc,sph relation. For example, the relation between M BH and Σ 5kpc,sph shown in the right-hand panel of Figure 3.
It can be expressed as, with the total rms scatter ∆ rms|BH = 0.59 dex. A plot of the rms scatter about the M BH -Σ R,sph relation as a function of R is shown in Figure 4. The scatter asymptotes to ∼ 0.58 ± 0.01 dex beyond 5 kpc. The low (0.69 dex) scatter about the M BH -Σ 1kpc,sph relation relative to the 0.95 dex scatter about the M BH -Σ 0,sph relation (Table 1) and the (soon to be discussed) M BH -Σ e,sph relations suggests that Σ 1kpc,sph is a better predictor of M BH than the latter projected mass densities. However, the M BH -Σ R,sph relations for R > 1 kpc is stronger, reflective of the separation of the µ (and Σ ) profiles at large radii.

Surface Brightness and Projected Density at the
Half-Light Radius: µ e,3.6µm,sph & Σ e,sph Using our 3.6 µm-sample, we see a positive trend between M BH and the surface brightness (µ e,sph ) at the projected half-light radius of spheroids ( Figure 5). A higher magnitude of µ e,sph corresponds to a lower luminosity density; thus, we find a declining relation between M BH and the effective luminosity density. We observe that ETGs and LTGs in our sample define two different M BH -µ e,3.6µm,sph relations, which are represented in the left-hand panel of Figure 5. ETGs define the following relation, log M BH M = (0.47 ± 0.04) µ e,3.6µm,sph − 19 mag arcsec 2 + (7.95 ± 0.11), (4) Figure 5. Black hole mass versus the bulge surface brightness µ e,3.6µm,sph at R e,sph (left-hand panel, Equations 4 and 5), the projected stellar mass density Σ e,sph at R e,sph (middle panel, Table 1), and the average stellar mass density Σ e,sph within R e,sph (right-hand panel, Table 1). ETGs and LTGs seem to define different relations in these diagrams. However, the complete picture of these relations is curved, as shown by the expected black and green curves for ETGs and LTGs, respectively. Note that the horizontal-axes in the middle and right-hand panels are inverted, and in the first panel, the horizontal-axis presents µ e,3.6µm,sph in a dimming order from left to right.
In order to include our full sample, we mapped µ e,sph (mag arcsec −2 ) to Σ e,sph (M pc −2 ) using Equation A5, and recover two trends defined by ETGs and LTGs in the M BH -Σ e,sph diagram. Similar trends due to ETGs and LTGs are observed in the M BH -( Σ e,sph , average projected density within R e,sph ) diagram. The M BH -Σ e,sph and M BH -Σ e,sph relations are depicted, respectively, in the middle and the right-hand panel of Figure  5. The fit parameters and the correlation coefficients for these distributions are provided in Table 1.
For early-type galaxies, the galaxy luminosity (or mass) has a curved relation with the galaxy surface brightness at/within any scale radius, R z,sph , enclosing (a non-zero) z% of the galaxy's total light (Graham 2019b). This includes the relation between galaxy luminosity and galaxy surface brightness at/within the half-light (z = 50%) radius, i.e., M * ,gal -µ e,gal or M * ,galµ e,gal (see Graham 2019b, their figure 3). Similarly for spheroids, the M * ,sph -µ e,sph and M BH -µ e,sph (also Σ e,sph , and Σ e,sph ) relations are expected to be curved, as shown in Figure 5. These curved relations for ETGs and LTGs are predicted using the M BH -n sph (Sahu et al. 2019b) and M BH -µ 0,3.6µm,sph relations defined by the two morphological classes for the M BH range of our sample, and further applying the equation µ e,sph = µ 0,sph + 2.5 b n /ln(10) (or Σ e,sph = Σ 0,sph − b n /ln (10)) for a Sérsic distribution. Here, Σ e,sph , µ e,sph , and µ e,sph can be related through equation 9 in (Graham & Driver 2005) and equation 11 in (Graham et al. 2006).
Thus, the slopes of the fitted M BH -µ e,sph lines (e.g., Equations 4 and 5) obtained here depend on the M BH and µ e,sph range of the fitted sample. Moreover, as the ETGs and LTGs seem to follow different M BHµ e,3.6µm,sph trends, the foreseeable M BH -µ e,3.6µm,sph (or Σ e,sph , or Σ e,sph ) curves shall also be different for the two galaxy types.
3.3.1. Offset between ETGs with and without a disk Sahu et al. (2019b) observed an offset of 1.12±0.20 dex in the M BH -direction, between ETGs with a disk (ESand S0-types) and ETGs without a disk (E-type), which defined almost parallel relations in the M BH -M * ,sph diagram. The calculation of M * ,sph is based on the spheroid Sérsic profile, quantified by the parameters n sph , R e,sph , and µ e,3.6µm,sph , and thus the offset between ES/S0-and E-types must have propagated to M * ,sph from these parameters. Further, Sahu et al. (2020) re-observed this offset (1.38 ± 0.28 dex in the M BH -direction) between ETGs with and without a disk in the M BH -R e,sph diagram, where again these categories defined almost parallel relations. (Sahu et al. 2020) did not report any significant offset between ETG subsamples in the M BHn sph diagram. Here, we next investigated if there is any such offset between ES/S0-and E-types in the M BHµ e,3.6µm,sph diagram.
Upon separating the ETGs with and without a disk, we do see the two groups offset from each other in the M BH -µ e,3.6µm,sph , M BH -Σ e,sph , and M BH -Σ e,sph diagrams ( Figure 6). However, the quality of fit for the two samples is poor (see Table 1); thus, it is difficult to quantify the offset accurately. Moreover, as discussed before, the complete M BH -µ e,3.6µm,sph rela-  Table 1. The dashed and dot-dashed curves represent the expected relations for E-type and ES/S0-type galaxies, respectively.
tions are curved. The expected M BH -µ e,3.6µm,sph curves for E-and ES/S0-types are shown in Figure 6. These curves are also calculated by using the M BH -n sph and M BH -µ 0,sph lines for the two populations combined with µ e,sph = µ 0,sph + 2.5 b n /ln(10). Additionally, as these curves are not parallel, the offset between the relations for E-type and ES/S0-type may not be constant throughout.
The declining M BH -(effective surface brightness) relation can also be inferred from the distribution of the spheroid surface brightness profiles for our sample shown in the right-hand panel of Figure 1. At the half-light radii (R/R e = 1) of spheroids, the surface brightness dims when going from low-M BH (blue) to high-M BH (red) profiles. Also, given the radially declining projected density profile, the ES/S0-types with a smaller R e,sph than the E-types, have a brighter µ e,sph (higher Σ e,sph and Σ e,sph ) than E-types hosting similar M BH . This is why the direction of the offset between E-and ES/S0-types in these diagrams ( Figure 6) is opposite to the offset seen in the M BH -R e,sph and M BH -M * ,sph diagrams, where ES/S0-types have a smaller R e,sph and M * ,sph than the E-types hosting a similar M BH .

BLACK HOLE MASS VERSUS SPHEROID SPATIAL DENSITY
We deprojected the (equivalent-axis) Sérsic surface brightness profiles of the spheroids to obtain their spatial (i.e., internal) mass density profiles, as described in the Appendix A. These internal density profiles 18 are displayed in Figure 7. We used a sequential blue-whitered color map to represent the central black hole masses in increasing order from low-mass (blue) to high-mass (red). The density profiles in the three panels of Figure 7 will help one understand the upcoming correlations observed between the black hole mass and the host spheroid's internal density at various radii.
Similar to the projected surface brightness profiles, the (deprojected) internal density profiles, ρ(r), are monotonically declining and can be characterized using the Sérsic surface brightness profile parameters (n, R e , µ e = −2.5 log I e , see Equation A3). Smaller and less massive spheroids, generally quantified by smaller Sérsic parameters (n and R e ), have a shallow inner density profile that descends quickly at outer radii (see the bluer profiles in the left-hand panel of Figure 7). On the contrary, more massive spheroids, generally indicated by higher Sérsic parameters (n and R e ), have a steeper inner density profile with a higher density and a shallower decline at large radii (see the red profiles in the left-hand panel of Figure 7).
The horizontal-axes in the middle and the righthand panels of Figure 7 are scaled using the sphere-ofinfluence radius (r soi ) of the black holes and the internal (or spatial) half-mass radius (r e,sph ) of the spheroids, respectively. This accounts for some of the different size scales used and will help with the understanding of the Figure 7. Spheroid internal stellar mass density profiles. The sequential color map blue-white-red depicts an increasing order of MBH. In the left-hand panel, the black and green colored stars mark the black hole's influence radius for the core-Sérsic and Sérsic galaxies, respectively. The horizontal axes are scaled with respect to the sphere of influence radius (rsoi) of black holes and the spheroid's spatial half-light radius re for the profiles shown in the middle and right-hand panels, respectively.
observed (M BH )-(spheroid internal density) relations revealed in the following sub-sections.

Spatial Density at the Black Hole's
Sphere-of-Influence: ρ soi,sph Based on the exact deprojection of the Sérsic model (Equations A3) the internal density near the spheroid center, ρ(r → 0), tends to infinity 19 for n > 1. Hence, as a measure of the central internal density, we chose the internal density at another central radius, where the gravitational potential of the black hole is in dynamical equilibrium with that of the host galaxy, known as the sphere-of-influence radius (r soi ) of the black hole. We denote the spheroid spatial density at r soi by ρ soi,sph .
We first calculated r soi using the following standard definition (Peebles 1972;Frank & Rees 1976;Merritt 2004;Ferrarese & Ford 2005), where σ is the host galaxy's central (projected) stellar velocity dispersion, which is likely to be dominated by the spheroid component of our galaxies. The stellar velocity dispersions of our galaxies are primarily taken from the HyperLeda (Makarov et al. 2014) database 20 , and are listed in Sahu et al. (2019a, their  at r = r soi . We also include the core-Sérsic galaxies in the M BH -ρ soi,sph diagram ( Figure 8), for whom ρ soi,sph is based on the de-projection of the (inwardly extrapolated) Sérsic component of their core-Sérsic surface brightness profile.
NGC 404, which is the only galaxy with an intermediate-mass black hole (IMBH) in our sample, is a genuine outlier, and it is possible that IMBHs may not follow the log-linear M BH -(stellar density) scaling relations defined by SMBH hosts. It has already been excluded from our correlations (as mentioned in Section 2). NGC 5055 has an unusually small central stellar velocity dispersion relative to its black hole mass 21 (see the M BH -σ diagram in Sahu et al. 2019a, their figure 2), resulting in a large r soi and thus a small ρ soi,sph .
Galaxies IC 2560, NGC 3079, NGC 4388, NGC 4826, and NGC 6323 have spheroid Sérsic indices between 0.58 and 1.15; thus, they have a shallow inner density profile and a small ρ soi,sph . The Sérsic galaxy NGC 0821, on the other hand, has a Sérsic index of 6.1 and hence, a steep inner density profile and a high ρ soi,sph . It also contains a faint edge-on intermediate-scale disk (Savorgnan & Graham 2016b), suggestive of an accretion event.
Including the above eight galaxies significantly biases the best-fit relation defined by most of the sam- Figure 8. Black hole mass versus internal stellar mass density at rsoi (top and middle panels) and black hole mass versus (averaged) internal stellar mass density within rsoi (bottom panel). The top panel shows a single regression (Footnote 22), where core-Sérsic galaxies (black stars) are distributed in a manner that suggests a different MBH-ρ soi,sph trend for these high-n sph systems. All the data points in this panel are color-coded according to their Sérsic indices. The middle panel shows the two different MBH-ρ soi,sph relations defined by Sérsic (blue) and core-Sérsic galaxies (red) galaxies. Similar substructure due to Sérsic (low n sph ) and core-Sérsic (high n sph ) galaxies are observed in the MBH-ρ soi,sph diagram (bottom panel). The excluded galaxies are named. See the text in Section 4.1 for details. Note that the horizontal axis is inverted such that the density decreases when going from left to right. ple; hence, we have excluded these galaxies (plus the Milky way) from our M BH -ρ soi,sph relations. Here, we do not exclude the stripped galaxies NGC 4342 and NGC 4486B (described in Section 2) because we are dealing with the spheroid spatial density at a central radius, which may not be affected by the outer mass stripping of these galaxies.
Initially, we performed a single regression between M BH and ρ soi,sph using our (Sérsic + core-Sérsic) sample 22 , as shown in the top panel of Figure 8. We noticed that the distribution of core-Sérsic galaxies traces a substructure systematically offset from the best-fit line for the ensemble of galaxies, suggesting a different trend for this sub-sample. Therefore, we further performed different regressions for the core-Sérsic and Sérsic galaxies, presented in the middle panel of Figure 8. We observed a tight, shallower relation for the core-Sérsic galaxies, given by log M BH M = (−0.68 ± 0.06) log ρ soi,sph 10 2.5 M pc −3 + (9.06 ± 0.05), with ∆ rms|BH = 0.21 dex. Curiously, this relation has the lowest total rms scatter of all the black hole scaling relations 23 . For Sérsic galaxies (with n 1), we found a relatively steeper relation, log M BH M = (−1.18 ± 0.10) log ρ soi,sph 10 2.5 M pc −3 + (8.39 ± 0.10), with ∆ rms|BH = 0.77 dex. The correlation coefficients for the above two relations are presented in Table 2.
Here, we needed to know M BH in advance to measure r soi and thus ρ soi,sph , voiding Equations 7 and 8 as black hole mass predictor tools but leaving them as constraints for simulations and to predict ρ soi,sph for a given M BH when the host spheroid surface brightness parameters are not known (as done in Biava et al. 2019, using other black hole scaling relations). The scatter in the above relations is smaller than that about the M BH -µ 0,sph (and M BH -Σ 0,sph ) relations, indicating that M BH has a better relation with ρ soi,sph , supporting the prediction in Graham & Driver (2007).

22
The single regression provides the relation log M BH /M = (−1.27 ± 0.07) log ρ soi,sph /10 2.5 M pc −3 + (8.55 ± 0.07), with ∆ rms|BH = 0.77 dex. 23 It is noted that r soi is derived from M BH (Equation 6). For a roughly similar σ among core-Sérsic galaxies, those with bigger M BH have a larger r soi (see the left-hand panel of Figure 7). The slope in Equation 7 tracks the average slope across 20-1000 pc of the high-n (red) profiles in Figure 7.
On their own, the core-Sérsic galaxies appear to have no correlation in the M BH -µ 0,sph diagram ( Figure 2). However, the overlapping nature of (the Sérsic component of) their density profiles in the left-hand panel of Figure 7 (also see Dullo & Graham 2012, their figure 18), coupled with Equation 6, supports the tight trend for cored galaxies seen in Figure 8. The smaller scatter observed for the core-Sérsic relation can be understood from the tight distribution of black points marking r soi and ρ soi,sph on the density profiles of the core-Sérsic spheroids in the left-hand panel of Figure 7. The green points, marking r soi and ρ soi,sph on the density profiles of the Sérsic spheroids are more scattered, explaining the higher rms scatter about the M BH -ρ soi,sph relation for the Sérsic galaxies. Figure 7 (left-hand panel) also explains why there will be a correlation between black hole mass and the isophotal or isodensity radius measured at faint/low densities. It is easy to see that the use of ever-lower densities will result in an ever greater separation of the curves. A result due to the different Sérsic indices (n) and the trend between M BH and n (e.g., Graham & Driver 2007;Sahu et al. 2020).
As with the M BH -ρ soi,sph , a similar apparent separation between core-Sérsic and Sérsic galaxies is recovered in the M BH -ρ soi,sph diagram involving the average spatial density within r soi , as shown in the bottom panel of Figure 8 (see Table 2 for fit parameters). However, the scatter is a bit higher than about the M BH -ρ soi,sph relations. Again, it should be noted that the values of ρ soi,sph (and ρ soi,sph ) for the core-Sérsic spheroids are higher than the actual values because these are based on the de-projection of the (inwardly extrapolated) Sérsic portion of their surface brightness profiles, which intentionally do not account for the deficit of light (see footnote 14) in the core, r R b .
For the core-Sérsic galaxies, the M BH ∝ (stellar mass deficit: M 0.27 * ,def ) relation (Dullo & Graham 2014, their equation 18) suggests that galaxies with high M BH have a higher mass deficit. Upon accounting for the mass deficit to obtain the actual ρ soi,sph,core (and ρ soi,sph,core ), all the core-Sérsic galaxies will move towards a lower ρ soi,sph (and ρ soi,sph ), i.e., towards rightside in Figure 8 (where the horizontal axes are inverted). However, galaxies with higher M BH shall shift more than the galaxies with lower M BH , generating a slightly shallower (negative/declining) slope than the slope of the relation presented here (Equation 7), but still preserving the apparent core-Sérsic versus Sérsic substructuring.
The negative correlations between M BH and ρ soi,sph (and ρ soi,sph ) can be visualized from the vertical ordering of blue-to-red shades (i.e., low-to-high M BH ) of the spheroid density profiles, shown in the middle panel of Figure 7, with the radial-axis normalized at r soi . Broadly speaking, at the influence radius (and any fixed multiple of this radius substantially beyond r/r soi = 1), the stellar density increases while going from the high-M BH (reddish profiles) to low-M BH (bluer profiles). The general (negative) M BH -ρ soi,sph trend for our sample arises from massive black holes having larger spheresof-influence, relative to low-mass black holes, combined with the spheroid's radially declining density profiles. However, the resultant relations for the core-Sérsic and Sérsic galaxies are dependent on the sample selection and, thus, the range of Sérsic profiles included in each subsample, as discussed in the following subsection.

Investigating the core-Sérsic versus Sérsic substructure
One may wonder if the substructures in the M BHρ soi,sph (and M BH -ρ soi,sph ) diagrams seen between core-Sérsic and Sérsic galaxies may be related to a similar division observed in the L-σ and M BH -σ diagrams (see Davies et al. 1983;Held & Mould 1994;Matković & Guzmán 2005;Bogdán et al. 2018;Sahu et al. 2019a). This may be because some of the division seen in the M BH -ρ soi,sph and M BH -ρ soi,sph diagram may be influenced by the use of the central stellar velocity dispersion while calculating r soi . Or conversely, the substructures observed in the M BH -σ diagram may partly be a reflection of the M BH -ρ soi,sph (or ρ soi,sph ) relations, if ρ soi,sph influences σ.
To test this connection, we tried an alternative estimation of the black hole's influence radius denoted by r soi,2BH . The radius r soi,2BH marks the sphere within which the stellar mass is equivalent to twice the central black hole's mass (Merritt 2004). Upon using the internal density (ρ soi,2BH,sph ) calculated at r soi,2BH , we recover the substructure between core-Sérsic and Sérsic galaxies in the M BH -ρ soi,2BH,sph diagram 24 (not shown), albeit with an increased scatter. This test demonstrated that the substructuring seen in Figure 8 is not due to the propagation of σ via Equation 6.
To investigate another scenario underlying the apparent substructures in the M BH -ρ soi,sph (and M BHρ soi,sph ) diagrams, we color-coded the data points in the top panel of Figure 8 according to their Sérsic indices. This Sérsic index color map divides the data in 24 The core-Sérsic galaxies follow the relation log M BH /M = (−0.65 ± 0.07) log ρ soi,2BH,sph /10 2.5 M pc −3 + (8.73 ± 0.09), and Sérsic galaxies follow log M BH /M = (−0.98 ± 0.07) log ρ soi,2BH,sph /10 2.5 M pc −3 + (7.76 ± 0.10), with ∆ rms|BH = 0.29 dex and 0.87 dex, respectively. the M BH -ρ soi,sph diagram in different diagonal zones, in a sequential order of n sph , such that one can obtain a set of M BH -ρ soi,sph relations applicable for different ranges of n sph . For example, roughly, we can point out three zones in the top panel of Figure 8: the excluded data points near the bottom right of the plot with the smallest Sérsic indices (n sph 1.5); the blue-purple-magenta points with 1.5 n sph 5 in the middle, and the redorange-yellow points with n sph 5 in the upper-left part of the diagram. Most of our core-Sérsic galaxies fall in the third zone, which is why we observe them defining a different M BH -ρ soi,sph relation than the majority of the Sérsic galaxies which fall in the second zone.
The distribution of data-points in the top panel of Figure 8 can be better represented on an M BH -ρ soi,sphn sph plane. This plane will be investigated in our future exploration of a black hole fundamental plane. We note that our calculation of ρ soi,sph depends on n sph and M BH , and thus these terms are not independently measured quantities. As noted, a high M BH , associated with a large n sph (see Sahu et al. 2020, for the M BH -n sph relation), will generate a large r soi and thus lower ρ soi,sph .

Spatial Mass Density within 1 kpc: The Spheroid
Spatial Compactness ρ 1kpc,sph The internal mass density is a better measure of the inner density than the projected column density. Hence, we introduce ρ 1kpc,sph , the spatial version of the projected spheroid compactness Σ 1kpc,sph (Section 3.2), defined as the mean internal stellar mass density within the inner 1 kpc of the spheroids.
We find a positive correlation between the black hole mass and the spheroid spatial compactness without any detectable substructuring due to the morphological classes of galaxies. The single-regression M BHρ 1kpc,sph relation 25 , shown in the left-hand panel of Figure 9, can be expressed as log M BH M = (2.96 ± 0.21) log ρ 1kpc,sph 10 0.5 M pc −3 + (8.47 ± 0.07), and has ∆ rms|BH = 0.75 dex. The M BH -ρ 1kpc,sph relation is marginally steeper than the M BH -Σ 1kpc,sph relation (Equation 2), and has a slightly higher vertical scatter. However, the orthogonal (perpendicular to the best-fit line) scatter in both the diagrams is comparable (∼0.24 dex).
We find positive trends between M BH and the internal spheroid density within other constant radii (e.g., 0.1 kpc, 5 kpc, 10 kpc) as well. The left-hand panel in Figure 7 shows that, in general, the high-M BH profiles reside above the low-M BH profiles at all radii; thus, the galactic spheroids with higher M BH are relatively denser than the spheroids with lower M BH , when compared at a fixed physical radius. This partly explains the positive trends obtained for the correlations of black hole mass with the spatial compactness, ρ 1kpc,sph , and the internal density at/within any fixed spatial radii. However, there is a varying scatter in the relations that decreases with larger radii.
A plot of the vertical ∆ rms|BH versus r sph for the M BHρ r,sph relations is shown in Figure 10. For r sph < 1 kpc, the M BH -ρ r,sph relations have a higher scatter than Equation 9, whereas, for r sph > 1 kpc the M BH -ρ r,sph relations are relatively stronger and have a gradually decreasing scatter with increasing r sph , analogous to the M BH -Σ R,sph relations (Section 3.2). This can be readily understood by again looking at the left-hand panel of Figure 7, even though it shows the density profiles, ρ, rather than the somewhat similar mean density profiles, ρ . There, one can see a cleaner separation of profiles of different M BH (and Sérsic index, n) when moving to larger radii, which is due to the increasingly longer tails of the high-n light profiles. For a comparison, the M BH -ρ 5kpc,sph relation (see the right-hand panel of Figure 9), which has ∆ rms|BH = 0.61 dex, can be expressed as, log M BH M = (1.99 ± 0.11) log ρ 5kpc,sph 10 −1.5 M pc −3 + (7.85 ± 0.06).
The smaller scatter in the above relation when compared to the M BH -ρ 1kpc,sph relation, and the quasisaturation of ∆ rms|BH for r sph 5 kpc (Figure 10), suggests that ρ 5kpc,sph can be preferred over ρ 1kpc,sph to predict M BH . Overall, the M BH -ρ r,sph relations are steeper than the M BH -Σ R,sph relations for any fixed spheroid radius (r = R), with a marginally higher vertical scatter and similar orthogonal scatter. Hence, potentially both properties ( ρ r,sph and Σ R,sph ) of a spheroid are equally good predictors of the central black hole's mass.

Internal Density at and within the Spatial
Half-Light Radius: ρ e,int,sph & ρ e,int,sph Using the spheroid internal density profiles, we calculated the spheroid spatial half-mass radius, r e,sph , which represents a sphere enclosing 50% of the total spheroid mass (or luminosity, for a constant mass-to-light ratio).  The ratio r e,sph /R e,sph is approximately 1.33 (Ciotti 1991).
We find that ETGs and LTGs define different (negative) trends between M BH and the internal stellar mass density (ρ e,int,sph ) at r = r e,sph , as shown in panel-a of Figure 11. The M BH -ρ e,int,sph relation followed by ETGs can be expressed as log M BH M = (−0.64 ± 0.04) log ρ e,int,sph M pc −3 + (7.81 ± 0.10), with ∆ rms|BH = 0.73 dex. The steeper relation followed by LTGs, with ∆ rms|BH = 0.69 dex, is given by log M BH M = (−1.02 ± 0.13) log ρ e,int,sph M pc −3 + (7.20 ± 0.11). (12) These two relations have a smaller scatter than the M BH -Σ e,sph relations for ETGs and LTGs (Table 1). The relatively smaller scatter and smaller uncertainties on the fit parameters suggests that ρ e,int,sph can be a better predictor of M BH than Σ e,sph (see Table 1). As we have repeatedly found, the shallower slope for the ETGs is physically meaningless. Its value reflects the sample selection and thus the relative number of ETGs with and without a disk. Further analysis of the M BHρ e,int,sph diagram reveals an offset between the ETGs with a rotating stellar disk (ES, S0) and ETGs without a rotating stellar disk (E), as shown in panel-c of Figure 11. The parameters for the M BH -ρ e,int,sph relations obtained for the two ETGs sub-populations are presented in Table 2. Notably, these two sub-categories of ETGs follow steeper M BH -ρ e,int,sph relations than Equation 11, almost parallel to each other but offset from each other by more than an order of magnitude in the M BH -direction. This offset is analogous to the offset found in the M BH -M * ,sph (Sahu et al. 2019b), M BH -R e,sph (Sahu et al. 2020), and M BH -Σ e,sph diagrams (Section 3.3). This offset originates from the smaller effective sizes (R e,sph ) and higher Σ e,sph of the ES/S0type galaxies relative to that of E-type galaxies possibly built from major mergers.
Similar trends and morphological substructures are found between M BH and the average internal density, Figure 11. Black hole mass versus internal stellar mass density at r = r e,sph (left-hand panels) and within r e,sph (right-hand panels). Top panels show that ETGs and LTGs follow two different MBH-ρ e,int,sph relations (panel a, Equations 11 and 12) and MBH-ρ e,sph relations (panel b, Table 2). The bottom panels present only ETGs, where ETGs with a disk (ES-and S0-types) and ETGs without a disk (E-type) are found to follow almost parallel MBH-ρ e,int,sph (panel c) and MBH-ρ e,sph (panel d) relations, offset in the vertical direction by more than an order of magnitude (see Table 2 for best-fit parameters). Note that the horizontal axes of all the panels are inverted, such that the internal density decreases when going from left to right. ρ e,int,sph , within r e,sph (see panels b and d of Figure  11). The parameters for the M BH -ρ e,int,sph relations are provided in Table 2. The right-hand panel in Figure  7 presents the spheroid spatial density profiles for our sample with the radial-axis normalized at r e,sph . At the spatial half-light radius, where log(r/r e,sph ) = 0, the increasing spatial density when going from high-M BH to low-M BH profiles is quite clear. This explains the negative M BH -ρ e,int,sph (and ρ e,int,sph ) correlations.
Table 2 also provides the morphology-dependent relations obtained between M BH and the spatial density ρ e,sph within the projected half-light radius (r = R e,sph ), which are analogous to the substructures in the M BH -ρ e,int,sph diagram. The M BH -ρ e,sph relation defined by our ETGs is consistent with that of Saglia et al. (2016). However, they do not report any of the vital substructures in this diagram due to ETGs (E, ES/S0) and LTGs. Without this awareness of the host galaxy morphology, the slope and intercept of the M BH -ρ e,sph relation is meaningless because it is biased by the randomness of one's sample selection. Indeed, this is why our ETGs relation has a slope of -0.64 rather than roughly ∼ −1.1, as followed by the E-type galaxies, the ES/S0type galaxies, and the spiral galaxies (see Table 2).
Finally, we again note here that similar to the M BH -Σ e,sph (and Σ e,sph ) relations (see Figures 5 and 6 in Section 3.3), the complete picture of the M BH -ρ e,int,sph (and ρ e,int,sph ) distributions are curved, which may be revealed in future using a larger sample. The slopes of the linear relations presented here are dependent on the mass range of our sample. Internal Density at r soi (Figure 8 Note-Column names are same as in Table 1.
a After excluding eight significant outliers, marked in Figure 8, which can significantly affect the best-fit line for the ensemble of Sérsic galaxies.

Prediction of M BH
We have shown how and explained why the BH mass correlates with a range of projected and internal stellar densities of the host spheroid. Plotting the density profiles of the (123 − Milky Way=) 122 spheroids together in the same figure reveals that the spheroids with larger BH masses reside in profiles with larger halflight radii, higher Sérsic indices, and longer tails to the light profile (Figures 1 and 7). At larger radii, the separation of spheroids with low-M BH and low-n profiles from those with high-M BH and high-n profiles becomes cleaner. Consequently, and counter-intuitively, the use of densities calculated at larger radii yields less scatter in the M BH -density diagram (Section 3.2 and 4.2).
The M BH -Σ 5kpc,sph relation (Equation 3) and the M BH -ρ 5kpc,sph relation (Equation 10) have similar rms scatters (0.59 dex and 0.61 dex), and are applicable to all galaxy types. The scatter in these diagrams is comparable to the morphology-dependent M BH -M * ,sph relations (cf., 0.50 dex, 0.57 dex, and 0.64 dex for E-, ES/S0-, and S-types, respectively) and the M BH -R e,sph relations (cf., 0.59 dex, 0.61 dex, and 0.60 dex for E-, ES/S0-, and S-types, respectively), and smaller than the morphology-dependent M BH -n sph relation (cf., 0.73 dex and 0.68 dex for ETGs and LTGs, respectively). Thus, Σ 5kpc,sph and ρ 5kpc,sph can predict M BH as good as predicted using M * ,sph and R e,sph , and better than M BH predicted using n sph . However, the density at 5 kpc may be very low for spheroids with R e less than half a kpc, and these relations become more of a reflection of the M BH -n relations, and to a lesser degree the M BH -R e relations (Sahu et al. 2020). The 3.6 µm M BH -µ 0,sph (Equation 1) and the M BHµ e,sph relations (see Table 1) offer an alternative way to predict M BH using µ 0,sph or µ e,sph , just from a calibrated (3.6 µm) spheroid surface brightness profile, without requiring either galaxy distance (for local galaxies where the cosmological corrections are very small) or a stellar mass-to-light ratio which can be complicated to choose. However, due to a higher scatter about these relations, the error bars on the predicted M BH will be higher than obtained using the M BH -n sph and M BH -R e,sph relations (see Sahu et al. 2020). The values n sph and R e,sph can also be obtained from an uncalibrated surface brightness profile. Plausibly, the high scatter in the M BH -µ 0,sph diagram is due to the use of a column density, and the high scatter in the M BH -µ e,sph diagram arises from a curved distribution of points.
For comparison, the M BH -φ relation for spiral galaxies (Seigar et al. 2008;Berrier et al. 2013) also has a small scatter of 0.43 dex (Davis et al. 2017), where φ is the pitch angle, i.e., the winding angle of the spiral arms. This relation can provide good estimates of M BH for spiral galaxies. Including all galaxy types, the M BHσ relation has a scatter of 0.53 dex; however, the M BH -σ diagram has different relations for core-Sérsic (cf., 0.46 dex) and Sérsic (cf., 0.55 dex) galaxies, which can provide a better estimate of M BH than the single relation, if the core-Sérsic or Sérsic morphology is known. Another, preferred relation to predict M BH may be the morphology-dependent M BH -M * ,gal relation (cf., 0.58 dex and 0.79 dex for ETGs and LTGs, respectively Sahu et al. 2019b), where, one does not need to go through the multi-component decomposition process to obtain the galaxy stellar mass, M * ,gal .
The tight M BH -ρ soi,sph relation for the core-Sérsic galaxies has the least total scatter (0.21 dex, see Table  2) among all the black hole scaling relations; whereas the M BH -ρ soi,sph relation obtained for the Sérsic galaxies has a higher scatter (0.77 dex). The relation for core-Sérsic galaxies only captures the upper envelope of high-n sph spheroids in the M BH -ρ soi,sph diagram, while the relation for Sérsic galaxies describes the average relation for spheroids with a medium value of n (between ∼ 1.5 to ∼ 5). Overall, the M BH -ρ soi,sph diagram suggests that the inclusion of n sph as a third parameter will lead to a black hole plane with a considerably re-duced scatter. However, if it was to turn out that the mass of the black hole is better connected to the stellar density within its sphere of influence and the stellar concentration (quantified by n), it is not useful for predicting M BH , because ρ soi,sph requires knowledge of r soi and thus M BH .

Dependence of the Black Hole Scaling Relations on the Galaxy Morphology
Sahu et al. (2020) did not report on the offset between the ETG subpopulations (E vs ES/S0-types) in the M BH -µ e,sph (or µ e,sph , or Σ e,sph , or Σ e,sph ) diagrams, that we reinvestigated here. Our investigation here has revealed an offset between the E-and ES/S0type galaxy samples ( Figure 6). However, the M BHµ e,sph correlations obtained for the E-and ES/S0-types are weak, and their slopes and the offset are not established. This is plausibly because they follow a curved relation with varying slopes, and we have sampled the bend points of the curves (see Figure 6). Consequently, there is not a strong correlation between M BH and the various effective densities for our sample (Section 3.3).
Morphology-dependent divisions in the M BH -n sph (ETG vs LTG), M BH -R e,sph (E vs ES/S0 vs LTG), and, as seen here, the M BH -µ e,sph (E vs ES/S0 vs LTG) diagrams, propagate into the M BH -M * ,sph (E vs ES/S0 vs LTG, Sahu et al. 2019b) diagrams. Similarly, these morphological substructures are also propagated to the M BH -ρ e,int,sph (and ρ e,int,sph ) diagrams presented here ( Figure 11). Although the ETGs and LTGs seem to define distinct tight relations, there is an order of magnitude offset 26 in the M BH -direction between ETGs without a disk (E-type or slow-rotators) and ETGs with a disk (ES/S0-types or fast-rotators). The offset between E-and ES/S0-type galaxies is a combined effect of a smaller bulge size (R e,sph ) and brighter µ e,sph (higher Σ e,sph and Σ e,sph ) of the ES/S0-type galaxies compared to that of E-type galaxies hosting a similar black hole mass (Section 3.3).
As discussed in Section 4.1, the Sérsic versus core-Sérsic division in the M BH -ρ soi,sph diagram ( Figure 8) remains independent of whether or not r soi is calculated using the central stellar velocity dispersion. Hence, the Sérsic versus core-Sérsic substructures observed in the M BH -σ diagram (Sahu et al. 2019a) and the M BHρ soi,sph (or ρ soi,sph ) diagrams are not directly related. Nonetheless, the M BH -σ and M BH -ρ soi,sph relations are respectively aware of the galaxy morphology, and thus the galaxies' evolutionary tracks and their central light/mass concentration, i.e., Sérsic index (see the top panel in Figure 8 and the description in Section 4.1).

Fundamental Black Hole Scaling Relation
Many studies have suggested that the M BH -σ relation may be the most fundamental/universal relation (e.g., Ferrarese & Merritt 2000;Gebhardt et al. 2000;Ferrarese & Ford 2005;de Nicola et al. 2019;Marsden et al. 2020) between a black hole and the host galaxy due to its obvious link with the galaxy's gravitational potential and the appearance of M BH ∝ σ 4−5 relations in theories trying to explain black hole feedback (Silk & Rees 1998;Fabian 1999). These claims are based on past observations (van den Bosch 2016; Saglia et al. 2016) which reported a single M BH -σ relation for all galaxy types (including bulge-less galaxies), and also a smaller scatter seen in the M BH -σ diagram relative to the M BH -M * ,sph relation. However, over the years, increments in the scatter about the M BH -σ relation to ∼ 0.5 dex with growing sample size (see the introduction in Sahu et al. 2019a), plus the revelation of a Sérsic (M BH ∝ σ ∼5 ) versus core-Sérsic (M BH ∝ σ ∼8 ) division in the M BHσ diagram (e.g., Bogdán et al. 2018;Sahu et al. 2019a;Dullo et al. 2020), undermine the perceived superiority of σ.
Importantly, if the relation with the least scatter should be the primary criteria for deciding the fundamental black hole scaling relation, recent studies further confound the situation. For example: the M BHρ soi,sph relation (Equation 7) for core-Sérsic galaxies has a total rms scatter of 0.21 dex; the M BH -(R b : break radius) relation for core-Sérsic galaxies has ∆ rms|BH = 0.29 dex (Dullo et al. 2020); the M BH -(pitch angle) relation for spiral galaxies has ∆ rms|BH = 0.43 dex (Davis et al. 2017); and the M BH -M * ,sph relation for ETGs has ∆ rms|BH = 0.52 dex (Sahu et al. 2019b). Moreover, the substructure in the M BH -ρ soi,sph diagram ( Figure 8) due to different ranges of n sph values suggest the existence of a possibly stronger M BH -ρ soi,sph -n sph plane, which shall be investigated in future work. Of course, ρ soi,sph is calculated using M BH , so some care will be required in such an exploration.

Super Massive Black Hole Binary Merger Timescale
The stellar density around a super massive black hole binary (SMBHB) plays an essential role in accelerating the merger of the black holes through dynamical friction (Chandrasekhar 1943;Begelman et al. 1980;Arca-Sedda & Capuzzo-Dolcetta 2014). During a galaxy merger, dynamical friction pushes the black holes towards the core of the galaxy merger remnant, forming a binary at parsec scales. The SMBHB goes through a hydrodynamical interaction with the surrounding stars (and dust/gas), entering a hardening phase, i.e. when the binding energy of the binary exceeds the average kinetic energy of stars around it (Holley-Bockelmann 2016). The binary then transitions from the hardening to the gravitational wave (GW) emission phase, which eventually drives the binary to merge (Celoria et al. 2018). The major part of a binary lifetime is spent in this transition phase/separation (Sesana & Khan 2015), the orbital frequency at this transition separation is known as the transition frequency. This time period (≈ binary lifetime) can be estimated using the average stellar density ( ρ soi ), and stellar velocity dispersion (σ soi ) at the sphere-of-influence of the binary and the binary's orbital eccentricity (e.g., Sesana & Khan 2015, their equation 7). The transition frequency, which is a part of GW strain model (discussed next), is also estimated using ρ soi , σ soi , and eccentricity (Chen et al. 2017, their equation 21).
Recently, Biava et al. (2019) estimated the SMBHB lifetime, as discussed above, using Sérsic parameters of a remnant-bulge hosting a given (binary) black hole mass. They used the M BH -M * ,sph relation (Savorgnan et al. 2016), the M * ,sph -R e,sph relation (Dabringhausen et al. 2008), and the M BH -n sph relation (Davis et al. 2019) to obtain the Sérsic parameters of bulges hosting 10 5 − 10 8 M binary black holes, with the assumption that the merger remnants follow these relations. Using these bulge parameters, they applied the Prugniel & Simien (1997) density model to obtain ρ soi to estimate the binary lifetime using the model from Sesana & Khan (2015, their equation 7). Now, using our M BH -ρ soi relations obtained here and the M BH -σ relations (e.g., Sahu et al. 2019a), one can directly obtain the ρ soi values and the central σ, respectively, for a given M BH , and using the central σ as a proxy for σ soi , one can estimate the typical binary lifetime more directly. One can also apply the expression of mean aperture correction for stellar velocity dispersion (from e.g., Jorgensen et al. 1995;Cappellari et al. 2006) to drive σ soi using the central σ (normalized at aperture size of 0.595 kpc) obtained from our M BH -σ relation and r soi . Similarly, using the ρ soi and σ values for a given M BH (and some binary eccentricity), the estimation of the transition frequency can be more straightforward (see Chen et al. 2017, their equation 21). This way, one would not need to go through various black hole scaling relations for the bulge parameters to obtain ρ soi , using an approximation for σ, and choosing an approximate density model, e.g., as suggested in Sesana & Khan (2015) and followed in Biava et al. (2019).
However, one should note that for galaxies with either a nuclear disk or nuclear star cluster, the ρ soi will be higher than estimated using the M BH -ρ soi relations for just spheroids. Whereas, for core-Sérsic galaxies, the ρ soi will be lower than estimated using the M BH -ρ soi relation.

Predicting the Gravitational Wave Strain
The long-wavelength gravitational waves (GWs: mHz -nHz), emitted during the SMBHB merger, fall in the detection band of pulsar timing arrays (PTAs: µHz -nHz), laser interferometer space antenna (LISA: 0.1 Hz to 0.1 mHz, Amaro-Seoane et al. 2017), and other planned space interferometers, such as TianQuin (Luo et al. 2016). These detectors aim to detect the stochastic GW background (GWB) and individual GWs, which are challenging to predict (Sesana et al. 2009;Mingarelli et al. 2017). The detectable amplitude (per unit logarithmic frequency) of perturbations due to the GWB is quantified by the characteristic strain (h c ), a typical estimate of which is required for different detectors sensitive to different wavelength ranges of GWs (e.g., see the sensitivity curves for various detectors in Moore et al. 2015).
The GWB characteristic strain can be modeled by integrating the SMBHB merger rate across redshift for a range of chirp-mass 27 (see the model described in Chen et al. 2019). The estimation of SMBHB merger rate is dependent on the observed galaxy mass function, galaxy pair fraction, SMBHB merger time scale (galaxy merger time scale + binary lifetime), and the (black hole)-galaxy scaling relations (Sesana 2013). The (black hole)-galaxy scaling relations convert the galaxy mass function and the galaxy pair fraction into the black hole mass function (BHMF) and the black hole pair fraction (BHPF).
Often, a constant M * ,sph /M * ,gal ratio has been combined with the old linear M BH -M * ,sph relation to obtain an M BH /M * ,gal ratio, which is used to convert the galaxy mass function into the BHMF (e.g., Shannon et al. 2015;Chen et al. 2019). This causes a bias in the estimated GWB characteristic strain (e.g., Mapelli et al. 2012, show that a quadratic M BH -M * ,sph relation, instead of a linear relation, changes the predicted extreme 27 Chirp mass of a binary comprising of objects with masses M 1 and M 2 is given by M = (M 1 M 2 ) 3/5 /(M 1 + M 2 ) 1/5 (e.g., see Cutler & Flanagan 1994). It influences the orbital evolution of the binary, e.g., the orbital frequency which governs the emitted GW frequency.
mass-ratio inspiral event rate by an order of magnitude). The use of our new morphology-dependent M BH -M * ,gal relations (Sahu et al. 2019b) will provide a direct way to obtain a better BHMF and BHPF. Coupled with these, the better estimates of the binary lifetime (Section 5.4) will improve the SMBHB merger rate, which will ultimately improve the predictions for the detectable GWB strain for PTAs and GW space missions.

Tidal Disruption Event Rate
The M BH -ρ soi relation obtained here can also help model the rate of tidal disruption events (TDEs, Hills 1975). This is important because apart from probing the black hole population and their environments (especially for BHs in inactive galaxies), TDEs are used to estimate the black hole mass (Mockler et al. 2019;Zhou et al. 2021), and electromagnetic counterparts of the extreme mass-ratio inspirals (EMRIs).
TDEs are expected to occur more frequently in galaxies with an elevated central stellar density or a nuclear star cluster (Frank & Rees 1976). The TDEs also require M BH 10 8 M because the weaker tidal forces at, and beyond, the Schwarzchild-Droste radii (Schwarzschild 1916;Droste 1917) of more massive black holes are insufficient to tear open stars and produce a TDE (Rees 1988;Komossa 2015). The TDE rate (Γ TDE ) versus ρ soi relation in Pfister et al. (2020, their equation 8) provides a lower limit of Γ TDE for a given ρ soi . Combining their Γ TDE -ρ soi relation with our M BH -ρ soi relation for Sérsic galaxies, we can obtain a relation between M BH and TDE rate as which can be used to obtain a typical estimate of the TDE rate for a given M BH . This can be refined further through the use of a set of M BH -ρ soi relations, applicable for different ranges of Sérsic index, or the creation of an M BH -ρ soi -n sph plane. We shall leave this for future work. It is worth noting that exact estimates of Γ TDE can vary depending on the presence of a nuclear star cluster.

CONCLUSION
We used the largest-to-date sample of galaxies which have a careful multi-component decomposition of their projected surface brightness profile (Savorgnan & Graham 2016b;Davis et al. 2019;Sahu et al. 2019b) and a directly-measured central black hole mass present in the literature (Section 2). We build upon our recent (published) work, where we revealed morphology-dependent M BH -(M * ,sph and M * ,gal ) relations (Davis et al. 2018;Davis et al. 2019;Sahu et al. 2019b), M BH -σ relations (Sahu et al. 2019a), and M BH -(n sph and R e,sph ) relations (Sahu et al. 2020).
Here, we investigated the connection between the black hole mass and the host spheroid's projected and internal stellar mass densities (Sections 3 and 4, respectively). More specifically, we presented the scaling relations of M BH with the spheroid projected luminosity density (µ, mag arcsec −2 ) and projected stellar mass density (Σ and Σ , M pc −2 ) at and within various spheroid radii (e.g., R = 0, 1 kpc, 5 kpc, and R e,sph ).
Importantly, we explored the correlation of M BH with the internal stellar mass density ρ (M pc −3 ), which is a better measure of density than the projected column density. We deprojected the (Sérsic) surface brightness profiles of our galactic spheroids using the inverse Abel transformation (Appendix A) and numerically calculated the internal densities at various internal radii, including the black hole's sphere-of-influence radius (r soi ), fixed physical internal radii (e.g., 1 kpc, 5 kpc), and the spatial half-mass radius r e,sph . We investigated possible correlations between M BH and the internal stellar mass density at and within these spheroid radii. We also presented the density profiles (Figure 7), which help in understanding the various observed M BH -ρ correlations ( Table 2).
In all these cases, we explored the dependence of the black hole scaling relations on the host galaxy morphology, i.e., possible division/substructure in the scaling diagrams due to ETGs versus LTGs, Sérsic versus core-Sérsic spheroids, barred versus non-barred galaxies, and galaxies with and without a stellar disk. The main results are summarized below.
• Spheroids with higher M BH have a brighter central surface brightness µ 0,sph (Equation 1) or higher central projected stellar mass density Σ 0,sph (Figure 2). This is true for Sérsic spheroids without depleted cores. This is qualitatively consistent with the linear M BH -µ 0,sph relation predicted in Graham & Driver (2007). However, the total rms scatter in the M BH -µ 0,sph (and Σ 0,sph ) diagrams are notably high (∼ 1 dex, see the fit parameters in Table 1).
• M BH defines a positive correlation with the average projected density, Σ 1kpc,sph , within the inner 1 kpc of the host spheroid (aka the spheroid compactness). The relation has ∆ rms|BH = 0.69 dex (Equation 2), and is followed by all galaxy types (see the left-hand panel of Figure 3, and Section 3.2).
• M BH has a stronger correlation with Σ R,sph for R > 1 kpc, than with the spheroid compactness Σ 1kpc,sph , such that the slope of the relation and the scatter decreases with increasing R. The total scatter starts saturating at ∼ 0.59 dex beyond ∼ 5 kpc (see Figure 4).
• In the M BH -µ e,sph and M BH -Σ e,sph (and Σ e,sph ) diagrams, ETGs and LTGs (S-types) follow different negative relations ( Figure 5, Table 1). The negative trend is because spheroids with higher M BH have a larger half-light radius with a lower density at/within these radii relative to that of spheroids with lower M BH . Further investigation reveals an offset between the E-and ES/S0-type galaxies in these diagrams, with suggestively similar slopes as that of LTGs (see Figure 6). However, the correlation coefficients are very poor, and the high scatter across these relations makes it difficult to quantify this offset correctly. Moreover, the actual distributions for the E-, ES/S0-, and S-types are expected to be curved; the predicted curves are also presented in Figures 5 and 6 (Section 3.3).
• M BH correlates with the internal density at and within the corresponding sphere-of-influence radius (ρ soi,sph and ρ soi,sph , Figure 8). The Sérsic and core-Sérsic galaxies seem to define two different relations with a negative slope. The core-Sérsic galaxies define a shallower M BH -ρ soi,sph relation with ∆ rms|BH = 0.21 dex, whereas, the Sérsic galaxies with n 1 follow a steeper relation with ∆ rms|BH = 0.77 dex (see Table 2). This substructuring is primarily due to the range of high Sérsic index profiles for the core-Sérsic spheroids (see the top panel of Figure 8 and Section 4.1). The data suggests an M BH -ρ soi,sph -n sph plane, which will be the subject of future work.
• Analogous to the (projected) spheroid compactness Σ 1kpc,sph , we introduced the spheroid spatial compactness, ρ 1kpc,sph , which is a measure of density within a sphere of 1 kpc radius. The quantity ρ 1kpc,sph defines a positive correlation with the M BH , which has ∆ rms|BH = 0.75 dex (see Equation 9 and the left-hand panel in Figure 9). As with Σ 1kpc,sph , we do not find a morphological dependence in the M BH -ρ 1kpc,sph diagram.
• Analogous to the M BH -Σ R,sph diagram, we find stronger correlations between M BH and ρ r,sph for r > 1 kpc. The slope of the M BH -ρ r,sph relation and the total scatter decreases with increasing internal radius r, where ∆ rms|BH asymptotes at ∼ 0.6 dex for r 5kpc (Figure 10). The M BH -ρ 5kpc,sph relation (Equation 10) is shown in Figure 9. Given the comparable scatter in the M BH -Σ R,sph and M BH -ρ r,sph diagrams, both the relations seem equally good predictors of M BH , where the density within 5 kpc is preferred over the density within 1 kpc (Section 4.2).
• In the M BH -ρ e,int,sph and M BH -ρ e,int,sph diagrams, ETGs and LTGs appear to define two different relations with a negative slope (top panels in Figure 11). Further analysis reveals that ETGs with a disk (E) and ETGs without a disk (ES/S0) appear to follow two different almost parallel relations, offset by more than an order of magnitude in the M BH -direction (bottom panels in Figure 11). They roughly have the same slope (∼ −1) as the relation for LTGs (Table 2). However, the relation may be curved, in which case the observed slope is a function of our sample's range of black hole mass. This morphology-dependent pattern has also been seen in the M BH -M * ,sph (Sahu et al. 2019b), M BH -R e,sph (Sahu et al. 2019a), and M BH -Σ e,sph diagrams ( Figures 5 and 6).
The revelation of morphology-dependent substructure in diagrams of black hole mass with various host spheroid/galaxy properties makes it more complex to conclude which relation may be the best to predict M BH or the most fundamental relation. It also rewrites the notion of the coevolution of galaxies and their black holes. The black holes appear to be aware of the galaxy morphology and thus the formation physics of the galaxy.
The central densities (µ 0,sph , Σ 0,sph , ρ soi,sph , and ρ soi,sph ) are based on the inward extrapolation of the Sérsic component of the spheroid's surface brightness model; however, additional nuclear star clusters or partially depleted cores will modify these densities. In future work, we hope to use high-resolution HST images to measure the depleted cores of the core-Sérsic galaxies and extract the nuclear star clusters from the host galaxy profile. This will enable us to revisit the M BHcentral density relations.
The M BH -density relations revealed in this paper have a wide range of applications (Section 5). For example: an alternative way to estimate the black hole mass in other galaxies; forming tests for realistic simulated galaxies with a central black hole; estimating the SMBH binary merger time scales; constraining the orbital frequency of the SMBHB during the transition from binary hardening to the GW emission phase; modeling the tidal disruption event rates (e.g., Equation 13); estimating/modeling the SMBH binary merger rate; and modifying the characteristic strains for the detection of long-wavelength gravitational waves for pulsar timing arrays and space interferometers.
This research was conducted with the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. This project was supported under the Australian Research Council's funding scheme DP17012923. This material is based upon work supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CAP 3 . I additionally thank the astrophysics group at the University of Queensland for hosting me and providing an office space for one year during the covid-19 pandemic. The surface brightness (projected/column luminosity density) profile of a galactic spheroid or an elliptical galaxy is very well described using the (Sérsic 1963(Sérsic , 1968) function, which can be expressed as It is parametrized by the Sérsic index (n), the scale radius (R e ), and the intensity (I e ) at R e . The term b is a function of n, defined such that the scale radius R e encloses 50% of the total spheroid light; therefore, R e is known as the (projected) effective half-light radius 28 . As noted by (Ciotti 1991), the exact value of b can be obtained using Γ(2n) = 2γ(2n, b) or it can be approximated as b = 1.9992 n − 0.327 for the value of n between 0.5 to 10 (Capaccioli 1989). The parameter n is the profile shape parameter and quantifies the central light concentration of the spheroid . The intensity, I e , is related to the surface brightness (µ in mag/arcsec 2 ) at R e , via µ e ≡ −2.5 log(I e ). As mentioned in Section 2, the bulge Sérsic profile parameters used here are taken from Sahu et al. (2020, their table A1). Who provide both major-axis bulge surface brightness parameters (n maj , R e,maj , µ e,maj ) and the equivalent-axis bulge surface brightness parameters (n eq , R e,eq , µ e,eq ) obtained by independent multi-component decompositions of galaxy surface brightness profiles along the major-axis and geometric mean-axis (equivalent to a circularised axis), respectively, (see Sahu et al. 2019b, for more details on the decomposition process). As described below, we used the equivalent-axis bulge parameters to utilize their circular symmetry while calculating the (deprojected) internal bulge density profile.
Using the inverse Abel (integral) transformation (Abel 1826;Anderssen & de Hoog 1990), the spatial luminosity density, j(r), for a spherical system can be expressed in terms of derivative (d I(R)/d R) of the projected luminosity density profile (see Binney & Tremaine 1987) as, where R represents a projected radius and r represents a 3D spatial (or internal) radius. Using the Sérsic profile (Equation A1) for I(R) and a stellar mass-to-light ratio, Υ λ (suitable for the corresponding image wavelength λ), to convert the luminosity into stellar mass, the spatial mass density profile, ρ(r) ≡ Υ λ j(r), can be expressed in a simplified integral form (Ciotti 1991;Graham & Colless 1997) as, ρ(r) = Υ λ I e be b πr r R e 1/n 1 0 e (−b(r/Re) 1/n /t) t 2 √ t −2n − 1 dt.
The above transformation (Equation A2) is applicable for a spherical system, and we can use the equivalent-axis bulge Sérsic parameters (n eq , R e,eq , and I e,eq ) to obtain ρ(r).
Using R e in parsec (pc), I e in solar luminosity per unit area (L /pc 2 ), and Υ λ in the units of solar mass per solar luminosity (M /L ), we obtain ρ(r) in the units of M /pc 3 from the above integral (Equation A3). The surface brightness µ e (in mag/arcsec 2 ) at the half-light radius can be mapped into I e (L /pc 2 ) using the following equation taken from Graham et al. (2006), where DM = 25 + 5 log[Distance (Mpc)] is the distance modulus, M ,λ is the absolute magnitude of the Sun in the corresponding wavelength-band λ, and s is the physical size scale for a galaxy in pc arcsec −1 . The projected (or surface) mass density Σ (M pc −2 ) at any projected radius R is calculated via, −2.5 log(Σ R [M pc −2 ]) = µ R − DM − M ,λ − 2.5 log(1/s 2 ) − 2.5 log(Υ λ ).
The solution of the above integral (Equation A3) can be expressed with the Meijer-G function 29 (Meijer 1936;Bateman & Erdélyi 1953;Mazure & Capelato 2002), which we numerically calculated using a Mathematica script to obtain the internal densities at various spheroid radii, used in Section 4. Prugniel & Simien (1997, PS hereafter) provides a remarkably simple, and one of the closest, approximation to the deprojected Sérsic profile (Equation A3), which can be expressed as, Here, ρ e is the spatial mass density at r = R e , and p is a function of n, obtained by maximizing the agreement between this approximated ρ(r) profile (Equation A6) and the exactly deprojected (Sérsic) density profile (Equation A3 ). The value of p is given by p = 1.0 − 0.6097/n + 0.05563/n 2 for a radial range of 10 −2 r/R e 10 3 and index Figure 12. The spatial density at Re calculated using the PS model (ρe,approx) plotted against the numerically calculated (Equation A3) spatial density (ρe,exact).
range of 0.6 n 10 (Lima Neto et al. 1999;Márquez et al. 2000). On equating the total mass obtained from the projected Sérsic profile (Equation A1) with the total mass calculated using the PS spatial density profile (Equation A6), considering a constant mass-to-light ratio, one has ρ e = Υ I e 2R e b n(p−1) Γ(2n) Γ(n(3 − p)) .
Owing to its simple analytical form and the model parameters common to the Sérsic luminosity profile, the PS model makes it easy to estimate the internal density profile of elliptical galaxies and the spheroids of multi-component (i.e., ES-, S0-, and Spiral-type) galaxies. Thus, Equations A3 and A6, both, are applicable for a galaxy/component whose surface brightness profile can be described using a Sérsic function; however, Equation A3 can provide the most accurate value.
For core-Sérsic galaxies, i.e., galaxies with power-law + Sérsic spheroid surface brightness profiles, Terzić & Graham (2005, their equation 5) modified the PS model and presented an expression for the deprojected core-Sérsic spheroid density profile. However, as we do not have precise parameters for the power-law core of our core-Sérsic galaxies, we use the Sérsic part of their surface brightness profile and deproject its inward extrapolation to obtain the central/inner rho for core-Sérsic galaxies.
The approximation of the deprojected density profiles can be imprecise (< 10% diferrence at 0.01 < R/R e < 100 R e for n > 2) to emulate the actual density profiles at the central radii, especially for low Sérsic index spheroids. See the comparisons in Terzić & Graham (2005, their figure 4), Emsellem & van de Ven (2008), and Vitral & Mamon (2020). Therefore, for our black hole-internal density correlations, we prefer to use the numerically calculated internal densities from Equation A3.
In Figure 12, we have compared ρ e,approx at r = R e calculated using the PS model (Equation A7) against ρ e,exact , numerically calculated using Equation A3. Here, we see an almost one-to-one match between the two values, except for galaxies M 59, NGC 1399, and NGC 3377, the three offset galaxies in Figure 2 with n sph,eq 8.8. The two offset points shown in Figure 12 are M 59 and NGC 1399, whereas, for NGC 3377, the exact integral (Equation A3) did not converge to provide an appropriate value of ρ e,exact .
Given the agreement between the exact and approximate internal densities at r = R e for the majority of the sample, for NGC 3377, we have used ρ e obtained from the PS model. Similarly, for some instances, where the exact ρ(r) integral (Equation A3) did not converge or provide a valid density value, we used the internal densities obtained using the PS model (Equation A7). This does not have a significant effect on the best-fit relations presented here. The extended density profiles in Figure 7 are obtained using the PS model, as it can still explain the qualitative nature of the trends observed in the M BH -ρ diagrams.