Disc cloaking: Establishing a lower limit to the number density of local compact massive spheroids/bulges and the potential fate of some high-z red nuggets

The near-absence of compact massive quiescent galaxies in the local Universe implies a size evolution since $z\sim2.5$. It is often theorised that such `red nuggets' have evolved into today's elliptical (E) galaxies via an E-to-E transformation. We examine an alternative scenario in which a red nugget develops a rotational disc through mergers and accretion, say, at $1\lesssim z\lesssim2$, thereby cloaking the nugget as the extant bulge/spheroid component of a larger, now old, galaxy. We have performed detailed, physically-motivated, multi-component decompositions of a volume-limited sample of 103 massive ($M_*/\rm M_{\odot} \gtrsim 1\times 10^{11}$) galaxies within 110\,Mpc. Among our 28 galaxies with existing elliptical classifications, we found that 18 have large-scale discs, and two have intermediate-scale discs, and are reclassified here as lenticulars (S0) and elliculars (ES). The local spheroid stellar mass function, size-mass diagram and bulge-to-total ($B/T$) flux ratio are presented. We report lower-limits for the volume number density of compact massive spheroids, $n_\mathrm{c,Sph}\sim (0.17$-$1.2) \times 10^{-4}\,\rm Mpc^{-3}$, based on different definitions of `red nuggets' in the literature. Similar number densities of local compact massive bulges were reported by de la Rosa et al. using automated two-component decompositions and their existence is now abundantly clear with our multi-component decompositions. We find disc-cloaking to be a salient alternative for galaxy evolution. In particular, instead of an E-to-E process, disc growth is the dominant evolutionary pathway for at least low-mass ($1\times10^{10}<M_*/\rm M_{\odot} \lessapprox 4 \times 10^{10}$) red nuggets, while our current lower-limits are within an alluring factor of a few of the peak abundance of high-mass red nuggets at $1\lesssim z\lesssim2$.

While 'relic' galaxies, believed to be untouched red nuggets from the local Universe, have been found, they are certainly not common. In addition, relics show complex dynamical structures (Trujillo et al. 2014;Ferré-Mateu et al. 2017) that warrant being treated as more than a single component system. Could the growth of two-dimensional discs, rather than three-dimensional envelopes, be a prelevant mechanism transforming these galaxies? Such a two-phase disc-building scenario, rather than an envelope-building process, was suggested by Graham et al. (2011), and Nelson et al. (2012) subsequently reported on the rapid build up of discs at ∼ 1, which could yield S0 galaxies with old stellar populations by today. At stake is the evolutionary pathway for some/many of the Universe's old and massive galaxies, and perhaps also for some of the less massive, high-'red nuggets', which we term 'red pebbles', which may now be the bulges of some spiral galaxies.
The distinct size difference between the passive galaxy populations at high and low redshifts sparked deliberations as to the dominant evolution mechanism of massive early-type galaxies (ETGs). While regular size ETG are not rare at low and high-, one cannot say the same for the compact massive red nuggets. The red nuggets are rather abundant at high redshifts. For instance, Saracco et al. (2010) reported a comoving number density of compact galaxies at 0.9 < < 1.92 of 2.3 × 10 −5 Mpc −3 for * /M > 3 × 10 10 . Barro et al. (2013) plotted the evolution of the number density across different redshifts, finding that the number density, , of red nuggets peaked at a lookback time of 8 Gyr ( = 1), with peak ≈ 1.5 × 10 −4 Mpc −3 . However, regardless of the assorted size and mass criteria used to define red nuggets, it has repeatedly been claimed that these objects are virtually nonexistent at ∼ 0 (Trujillo et al. 2009;Taylor et al. 2010). The sharp decrease in the compact massive galaxy population indicates a drastic change in galaxy structure over this period of time. However, disagreement among observations (e.g., Trujillo et al. 2007;Valentinuzzi et al. 2010;Poggianti et al. 2013a) led to debates over whether there is indeed a drop in the number of compact massive quiescent systems across time.
The structure of an ETG is often overlooked and oversimplified. Many are misclassified as purely elliptical (E) galaxies while in reality, most ETGs contain a disc (Carter 1987;Capaccioli et al. 1987;Nieto et al. 1988;Capaccioli et al. 1990a;D'Onofrio et al. 1995;Graham et al. 1998;Emsellem et al. 2011;Scott et al. 2015); either a large-scale disc making it a lenticular (S0) galaxy, or an intermediatescale disc making it an ellicular 2 (ES) galaxy (Liller 1966;Graham 2019). Given that different structures, for example, triaxial spheroids or relatively flat discs, form via different physical mechanisms, advances may be made by not simply treating galaxies as if they are single-component systems but rather examining their bulge/disc nature and the role of S0 and ES galaxies in the grand scheme of galaxy evolution. 3 From a random incomplete sample, Graham et al. (2015) reported a lower-limit to the volume number density for local compact massive spheroids in fair agreement with some of the number densities reported for compact massive galaxies at high-(see also de la Rosa et al. 2016;Costantin et al. 2020). Here, we investigate the issue more thoroughly, using a mass and a volume-limited sample of nearby galaxies combined with a careful, homogeneous, physicallymotivated 4 , multi-component decomposition of their surface brightness profile. We do this to detect and quantify the primary spheroidal component, a.k.a. bulge.
While we do not advocate any particular origin for the high-red nuggets, we do note that some bulges may form in other ways. For example, star-forming clumps may form in turbulent disks at high-2 The concatenated name (elliptical + lenticular) 'ellicular' was introduced for the ES galaxy type in Graham et al. (2017). 3 Those unfamiliar with the historical discovery of, and continuity between, the different morphological sub-structures in large, high surface brightness galaxies may like to refer to the galaxy classification grid presented in Graham (2019). It builds upon past works by better recognizing the range of disc sizes in ETGs, which can also contain bars. 4 Rather than randomly/automatically fitting two or three Sérsic components, we fit physical structures such as bars, ansae (e.g., Martinez-Valpuesta et al. 2007), discs, etc., as detected through imaging and kinematic data. redshift before migrating inwards to form the bulge (e.g., Bournaud et al. 2007;Elmegreen et al. 2008). This process does not preclude subsequent disk accretion and growth but it does start with a disk rather than the monolithic collapse of a spheroid. In addition, disks (with and without a stabilising dark matter halo) can experience instabilities leading to bars, which can in turn experience instabilities of their own, leading to the formation of buckled bars often with distinct (peanut shell)-shaped structures known as pseudobulges (e.g., Bardeen 1975;Hohl 1975;Combes & Sanders 1981;Combes et al. 1990). While both low-mass classical bulges and pseudobulges can rotate and have a Sérsic index for their light profile around 1 (Graham 2015), the bar+pseudobulge has a tendency to result in elongated isophotes having a maximum 6 Fourier harmonic term at around half the length of the bar (Ciambur 2016, and references therein) 5 . Classical bulges, bars and pseudobulges, and disks are known to coexist (e.g., Norman et al. 1996;Erwin et al. 2003;Athanassoula 2005). In our multi-component decompositions, we focus on recovering what would be considered the "classical bulge". We will present the decompositions of the galaxy light in a forth-coming paper which will compare the sizes and masses obtained from fitting a single Sérsic function, a Sérsic plus an exponential function, and our multi-component analysis used here to obtain the spheroid size and mass.
In Section 2, we describe the data selection of local galaxies that potentially contain a compact massive spheroid. Section 3 goes through the process of galaxy model fitting and some of the nuances of multi-component decomposition. In Section 4, the spheroid mass function, the spheroid size-mass relation, and the number density of local compact massive spheroids are presented. A comparison with local spheroid and high-samples are shown in Section 5. Based on the newly-found local spheroids, we discuss the implications for galaxy evolution as a whole in Section 6. We summarise our findings in Section 6.5. Finally, in a set of appendices, we provide an accounting of distance corrections (Appendix A), assumptions governing our adopted mass-to-light colour relations (Appendix B), host galaxy sample (Appendix C), and structural parameters for the local spheroids (Appendix D).
Throughout this paper, we assume the following cosmographic parameters: 0 = 68 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7. The absolute magnitude of the sun ( ) in −band is set to 4.53 (in AB mag). Unless otherwise stated, all bulge parameters ( e , e and ) correspond to the geometric-mean axis ( eq = √︁ maj min , aka the equivalent axis) which is equivalent to a circularised form of the isophotes.

DATA
In Section 2.1, we explain the parent data selection process. Section 2.2 briefly describe the images sources. Section 2.3 explores the various stellar mass-to-light ratio implemented in this study. Section 2.4 describes our process of finding the appropriate distance for each galaxy. Section 2.5 is a general summary of our volume-limited sample.

Parent sample
To calculate the number of (compact, massive) spheroids per cubic Mpc, hereafter the number density, as well as study the statistical properties of these local spheroids, a well-defined volume-limited sample of galaxies is required. The galaxy images should also be of sufficiently high spatial resolution for us to properly separate the substructures such as discs and bars from the spheroidal components. This limits our sample's maximum co-moving distance to approximately 100 Mpc, because, beyond this point, the multi-component decomposition becomes challenging with ∼1 seeing.
We elected to use Sloan Digital Sky Survey (York et al. 2000) photometric data because of its wide sky coverage of ∼14,555 deg 2 (Aihara et al. 2011), and its completeness of galaxies at the bright end ( < 23 mag) of the galaxy magnitude distribution. Our sample selection is based on the NASA-Sloan ATLAS catalogue 6 , a data set that consists of information from the SDSS DR8 (Aihara et al. 2011), NASA Extragalactic Database (NED, Helou et al. 1991) 7 , Six-degree Field Galaxy Redshift Survey (Jones et al. 2004), Two-degree Field Galaxy Redshift Survey (Colless 1999), ZCAT (Huchra et al. 1983), and ALFALFA (Giovanelli et al. 2005). The catalogue provides a list of unique objects as they remove all the redundant data and stars that were previously misclassified as galaxies. It contains 145,155 galaxies within a redshift ∼ 0.055.
From the parent sample, we limit the angular boundaries to We define these arbitrary boundaries in the attempt to cover the majority of the northern sky in the SDSS field. It contains 55,370 galaxies and captures 8.5% of the sky, or a solid angle of 1.07 sr. Figure 1 shows the distribution of the angular selection. The area contains notable supergalactic structures such as the Virgo and Coma clusters.

The imaging data
Our imaging data come from SDSS −band images. We decided to use -band photometry for its smaller galactic extinction. The initial flux of the objects was previously measured with a single Sérsic component fit by the SDSS pipeline. The NASA-Sloan ATLAS catalogue provided the flux in the bands in units of nanomaggies. The SDSS magnitude of a galaxy in the AB magnitude system (Oke 1974) is: where is the flux of the object in nanomaggies and [22.5 mag] is the photometric zero-point magnitude for all bands in the SDSS. The seeing for each galaxy was estimated by using the Image Reduction and Analysis Facility (IRAF; Tody 1986Tody , 1993 function imexam. We typically measured 5 or 6 stars randomly distributed in each frame. The imexam task was used to fit a Gaussian function to each star. The median value of the full width at half maximum (FWHM) was used to represent the seeing of the frame. The measurements are listed in Column 8 of Tables C1-C3.

Stellar masses
Despite slight differences in definition between studies, compact massive galaxies generally have a stellar-mass larger than * /M ∼ (1-7) × 10 10 . We will focus on local host galaxies with a total stellar mass of * /M 1 × 10 11 , in order to find compact massive spheroidal components. The determination of the stellar mass depends heavily on different factors, including the initial mass function (IMF); stellar population synthesis (SPS) models; the treatment of dust attenuation; and the fine-tuning of the age, metallicity, and reddening parameters. While it is a common practice to obtain the stellar mass by fitting a stellar model to the spectral energy distribution (SED) of the galaxy light, it is not quite appropriate for our purpose. Having the galaxy SED does not provide us with the mass-to-light ratio ( * / ) of the spheroid, which would require a structural decomposition in many bands to provide the SED of the spheroid.
While bulge+disc fits in many bands are increasingly produced on scale (e.g. Kennedy et al. 2016), the proximity of our sample, in which bars and rings are resolved, mandates the implementation of multi-component decompositions. To date, these are performed manually rather than automatically, which allows for greater care and recourse to additional information, such as kinematic maps in the literature. Given time constraints, we elected to focus on the decomposition of −band images (described in detail in Section 3.3). It has also be shown that the -band mass-to-light ratio * / and ( − ) relation has a low scatter intrinsically (Taylor et al. 2011). We use the galaxy colour, rather than the full SED, as a proxy for the spheroid colour when estimating the spheroid * / ratio.
sume an exponentially declining star formation history (SFH ∝ / ) in their models. The following equations can be used to calculate the host galaxy total stellar mass, and later the bulge/spheroid stellar mass: Given the conditions on their respective assumptions, each MLCR self-reported to have the following intrinsic dispersions, Z09: 0.10-0.15 dex, T11: 0.1 dex, IP13: 0.13 dex, and RC15: 0.14 dex. For our sample, the galaxy ( − ) colours exhibit a familiar bimodal distribution (Tully & Verheĳen 1997;Strateva et al. 2003;Blanton et al. 2003;Kauffmann et al. 2003;Baldry et al. 2004;Brinchmann et al. 2004;Balogh et al. 2004). The most massive galaxies occupy the 'red' distribution, centered at ( − ) ∼ 1.2 mag. We demonstrate the log 10 ( * / )-to-( − ) relations for the four equations in Figure 2. Within our sample's colour range (∼ 1.2 ± 0.04 mag), the MLCRs yield decreasing * / ratio in the following order: Given the goal of this work is comparative in nature, it is important to know the stellar mass relative to the high-galaxies. While the underlying assumptions in calculating the stellar masses of highgalaxies have varied from paper to paper, we attempt to make the comparison using the MLCR best matches their assumptions 8 . Throughout this work, we present the result mainly using RC15 stellar mass for the sake of consistency. It represents an intermediate estimation between the lowest (T11) and the highest estimation (IP13 9 ). 8 The theoretical background for each MLCR is discussed in detail in Appendix B. 9 Considering there are many S0 and S galaxies in our sample, we apply the disc-model equation from IP13 (see their Table 5). Although, the two equations in IP13 (in their Table 3 and 5) for ( − ) colour would return a similar result. We provide several estimates of the stellar mass, with the MLCRs from IP13 yielding the highest masses here. It is noted that this Moreover, RC15 have the same priors as some prominent red nuggets studies (Barro et al. 2013;van der Wel et al. 2014;van Dokkum et al. 2015, , etc.), namely the Chabrier (2003) IMF and the Bruzual & Charlot (2003) SPS. That being said, the effect of using a different * / ratio prescription should be considered thoroughly. Having several stellar mass estimations allow us to examine the potential bias inherent in the choice of stellar mass estimation measurement. We shall explore how the number density of the embedded compact massive spheroids may change for each MLCRs, later in Section 5.4.

Distances
We try to account for the peculiar velocity influence on distance measurements. There are several empirical methods to approximate the distance of nearby galaxies: the Tully-Fisher (TF) relation (Tully & Fisher 1977), Fundamental Plane (Djorgovski & Davis 1987), and velocity field reconstruction (Peebles 1989(Peebles , 2001Phelps et al. 2006;Shaya & Tully 2013). To construct our volume-limited sample, the procedure involved several steps.
(i) Redshifts. We commenced this project by first calculating the distances using the corrected redshift (ZDIST) available in the NASA-Sloan ATLAS catalogue. The redshift values are based on the Mark III Catalog of Galaxy Peculiar Velocities (Willick et al. 1997, obtained using the TF or D n -relation. This catalog provides a Malmquist bias-corrected distance based on the velocity field provided by the IRAS 1.2 Jy redshift survey (Fisher et al. 1995) density field. The initial distances and stellar mass estimation are shown in the upper panel of Figure 3 as grey crosses.
(ii) A broad mass-distance sample selection boundary. We select our initial galaxy sample for further examination using the broad mass-distance selection criteria (blue line) in the upper left-hand corner of Figure 3, satisfying the arbitrary condition log 10 ( * /M ) > 10 + (0.014 × Dist.), and Dist. < 115 Mpc. This is intended to produce a manageable sample size of galaxies for us to further check if redshift-independent distances are available on NED. There are 708 data points (green crosses) remaining after imposing the mass-distance sample selection boundaries.
(iii) Updated redshift-independent distances. Each of these green points was checked in NED to see if there were redshift-independent measurements. We primarily focussed on obtaining distances related to stellar phenomena, such as surface brightness fluctuation (SBF), Cepheids, Supernovae Ia (SNIa), the tip of the red giant branch (TRGB), etc. When such measurements were not available for a galaxy, we resorted to using the velocity correction from NED, based on the linear infall velocity model in Mould et al. (2000. This model corrects for the influence imposed by major gravitational bodies, namely the Virgo Cluster, the Great Attractor, and the Shapley supercluster. The revised distances are shown as the coloured dots and purple crosses (labelled 'New Distance', see Column (4) in Table C1-C3) in the lower panel of Figure 3 10 . One can see changes from the upper panel. For example, one can see how the Virgo Cluster at ∼17 Mpc is not an extreme mass estimate, with MLCRs from Bell et al. (2003) and Schombert et al. (2019) yielding yet higher masses. 10 Note that, because we rely on redshift independent distances measured using stellar information, our distance is luminosity distance. Although, in this redshift range (0 < < 0.025), luminosity distance and comoving distance are similar. . The data selection for potential host galaxies containing compact massive spheroids. Top: The grey crosses (×) are the galaxies satisfying the angular sample selection boundaries in Equation 1. The green crosses (×) are the galaxies that satisfy Equation 4. All distances are from W1997. Bottom: Bins 1, 2, and 3 are shaded in green for the right, middle, and left boxes, respectively. All distances have been updated (see Appendix A). The red, blue, and black points correspond to the samples listed in Tables C1-C3, respectively. Violet crosses (×) either do not satisfy the Equation 5 conditions or do not consistently reside in the same bin if alternative distance estimates are adopted (see Figure A1). After careful examination, we excluded the false sample that are stars or galaxies with wrong velocities. The crosses not overlapping with the points in the lower panels are rejected from the analysis. materializes from the data, and quite a number of galaxies escape the blue-line selection boundary. One can similarly expect that some galaxies initially outside of the blue-line selection boundary (and ignored by us) would be shifted inside. We mitigated against this by limiting the final selection to galaxies in three regions, which contained the galaxies previously outside these bins based on the older distances. Within the broad mass-distance selection boundary sample of 708 galaxies, we defined three volumes, aka bins, to acquire both a manageable 11 and useful number of galaxies.
We chose these arbitrary bins to include the massive galaxies that potentially host a massive bulge ( * /M > 10 10 ). The data points of each bin are highlighted in different colours in the lower panel of Figure 3.
(iv) Further refinements. In this step, we further checked the appropriate distance of each galaxy in the three bins. In addition to the redshift-independent measurements, as well as the W1997 and M2000 velocities, we collected another independent set of distances using the newly available Cosmicflow-3 distance calculators (Kourkchi et al. 2020). The Cosmicflow project uses a collection of diverse distance measurements (Tully et al. 2008 to derive the velocity and density field fluctuations in the local Universe. Their two online calculators provide a distinct approach for galaxies with Dist. < 38 Mpc and distances greater than this and up to 200 Mpc. In the former case, the velocity field is based on the non-linear numerical action orbit reconstruction in Shaya et al. (2017) and the velocity model of Tully et al. (2014). For the latter case, the distance is inferred from the three-dimensional linear velocity field model in Graziani et al. (2019), with an emphasis on mitigating inherent biases such as Malmquist bias and the lognormal distribution of peculiar velocities. The main advantage of Cosmicflow is the large library (17,647 individual galaxies) of distances from which the velocity field is constructed.
Armed with 3 to 4 distinct distance estimations for each galaxy in our set of three bins, we found that no matter which distance is chosen, all of our sample lies within the three bins selection volume (within 110 Mpc). The sample provides a good representation of galaxies within these boundaries.
The distance we used for the calculation is chosen on a case-by-case basis. The -independent distance is prioritized if they are available. We have 20 and 1 -independent distances in Bin 3 and Bin 1, respectively. Unless the measurement has a significant error or disagrees significantly with the other measurements, -independent distance is used (in total 18). When -independent distances are not available, we resort to using the Cosmicflow-3 distances (in total 81). There are a few of Cosmicflow-3 distances that fall within the "triple-value" region, as the Virgo Cluster bends the distance-velocity relation into a cubic function and returns three distinct distance solutions. In such cases (in total, only 4 galaxies), we choose the W1997 distances instead.
The final distance we used in the calculation is shown in Column (2) in Tables D1-D3, and the details of specific cases are discussed in the Appendix A. After eliminating all the false samples (stars or galaxies with obviously wrong velocities), this leaves us with 115 galaxies spread over three bins.

Volume-limited galaxy sample
We defined the bins in a manner that would produce a manageable sample size for careful multi-component decomposition work. The volume of each bin and the total survey volume are shown in Table 1. It is labour intensive to create a detailed decomposition. As we will discuss in Section 3.3, the process involves manual attention and cross-referencing evidences from the literature. While it would be best to conduct large-scale studies on all the local massive host galaxies, the current limitation in automatic decomposition programmes prevent us from doing so. We mitigate the challenge by limiting a local volume with respect to the relevant mass range. In the parent sample selection (Equation 5), we have proceeded by intentionally using the highest mass estimations, which come from IP13, in order not to miss potential host galaxies of compact massive spheroids. This ensures the inclusion of even the least likely (the least massive) candidates within our stellar mass selection criteria. For easy comparison across different MLCRs, we present a table of conversion, the lower mass limit for each bin across the four MLCRs in Table 2. As one applied another MLCRs instead of IP13, the galaxies stellar mass within the three bins decreased by roughly a fixed amount ( * /M ) 12 . The lower mass limits for other MLCRs are obtained by reducing the IP13 lower mass limits by * /M . Ultimately, the choice of MLCR would not affect the analysis. As we move forward to discuss in RC15 stellar mass terms, the lower mass boundaries of the host galaxies are now 3.4 × 10 11 M , 1.3 × 10 11 M , and 6.7 × 10 10 M for Bin 1, 2, and 3, respectively.
Given each bin's differences in the lower galaxy mass that was included, throughout this paper, the result of each bin is presented separately. Tables C1-C3 show the galaxies' basic information from the three bins.

METHODOLOGY
This section showcase our data processing methods, including the reduction process (Section 3.1), modelling of the galaxy image (Section 3.2), and the multi-component decompositions (Section 3.3). 12 While * /M is different for each galaxy, the deviation is insignificant. The deviation between the maximum and minimum * /M is equal to 0.087 dex for Bin 1, 0.01 dex for Bin 2, and 0.002 dex for Bin 3.
Based on the result of the decomposition, we explore several important aspect about the embedded spheroids such as: the spheroid colour (Section 3.5), mass (Section 3.6) and size (Section 3.7). Additionally, we reclassify the morphology for each galaxy based on our multi-component decompositions (Section 3.4).

Data reduction
We obtained cutout images of the galaxies from the SDSS DR12 Science Archive Server (SAS) 13 . Each frame spans 2048 pixels × 1489 pixels (819. 2 × 595. 6). The flux originates from multiple sources, including the Poissonian background noise plus the signal coming from foreground stars and our target galaxy. Naturally, the distribution of the pixel intensity values in frames with a sufficient field-of-view to capture the sky background will be a combination of a symmetrical Poissonian bell curve (random noise) centred on a low photon count, plus pixels with target-galaxy light stretching to brighter intensities. We define our sky value as the median of the Poissonian distribution and subsequently subtract it from the frame (Almoznino et al. 1993;Davis et al. 2019;Sahu et al. 2019).
The use of wide-frame cutouts resolves the reported SDSS sky subtraction problem in the SDSS DR8 (West 2005;Blanton et al. 2005;Bernardi et al. 2007;Lauer et al. 2007;Hyde & Bernardi 2009;West et al. 2010;Blanton et al. 2011). The SDSS standard photometric pipeline (Lupton et al. 2001) estimates the background from a 100 × 100 frame and subtracts it from the image. In the case of large extended local galaxies, such an approach will result in over-subtraction because the frame is not wide enough to reach beyond the galaxy. Since our frame is significantly bigger than the target galaxy, the sky subtraction issue is under control.
After the sky subtraction, the images proceed to a masking procedure. It is important to mask out the excess light from foreground stars and nearby galaxies because they can prevent the fitting algorithm (Section 3.2) from creating an accurate model for the target galaxy. We divided the masking process into two stages.
(i) First, we exclude the majority of the bright sources using the SExtractor (Bertin & Arnouts 1996) automatic process. The programme identified the contamination by threshold sigma detection. While it is good at tracking isolated sources, SExtractor is not as effective in deblenidng overlapping objects, in particular, foreground stars in front of a galaxy along the line-of-sight. We handle this issue through the second stage: manual masking.
(ii) Additional foreground objects are identified by human eyes and concealed using the polygon tool from DS9 (Joye & Mandel 2003a). We masked the bright sources, as well as the dust lanes that influence the stellar light.
We construct the mask using the IRAF function mskregion. Combining both the automatic and manual masking, all the irrelevant objects are effectively removed from the frame.

Isophotal fitting
The galaxy images are modelled with the isophotal fitting programme ISOFIT (Ciambur 2015). The method was first described by Carter (1978Carter ( , 1987 and implemented by Jedrzejewski (1987). It has since been a staple technique in modelling early-type galaxies, most notably through the IRAF package ellipse. The algorithm fits a series of nested isophotes centred around the peak brightness region. It performs a least-square fit on each isophote with intensity, ( ), expressed by the Fourier series: where is the average intensity, and are the harmonic coefficients of th order, and is the azimuthal angle. ISOFIT made two significant changes to the Jedrzejewski (1987) algorithm. It replaces the azimuthal angle ( ) with the eccentric anomaly ( ), with the following relation: where is the varying ellipticity ≡ 1 − ( / ) of each nested isophote, and and are the lengths of the semi-major and semiminor axes, respectively.
The incorporation of the ellipticity information in the eccentric anomaly gives us a better description of the shape of each isophote, in particular those that have high . It substantially improves the modelling of edge-on disc galaxies and those with elongated bars. Given the dramatic improvement, ISOFIT additionally allows more higher harmonic coefficients, with amplitudes and typically asymptoting to zero by 12. The high-order Fourier coefficients capture the radial perturbations of each isophote. Each integer represents a specific type of deformation from a perfect ellipse. For example, 4 and 4 describe the amplitude of a four-cornered deformation of an ellipse, in turn showing the 'boxyness' and 'discyness' of the isophote. As illustrated in Ciambur (2015, Figure 4), the high-order ( 6) coefficient reveals additional information about a galaxy that is rarely studied, i.e., the 6 profile captures (peanut shell)-shaped structures nestled around bars (Ciambur & Graham 2016;Ciambur et al. 2017;Saha et al. 2018).
The 2D intensity map described by the ISOFIT task is comprised of a series of 1-D radial distributions, capturing the surface brightness, ellipticity, position angle, and Fourier harmonic terms. The ratio of the model error err = / √ to intensity ( ) at each isophote is shown in Figure 4, where is the standard deviation of the pixel values (data) about the two-dimensional reconstruction of the image created with the IRAF task Cmodel (Ciambur 2015) and is the number of pixels in each isophote. max is the cutoff radius to which we modelled each galaxy 14 . While the model error-to-intensity ratio increases with radius, they do not deviate much from err / = 0, with most points having | err |/ < 0.02. The error from the isophotal fitting is small and the surface brightness information is well captured by the fitting process.
It is common in the community to use two-dimensional (2-D) fitting algorithms such as GIM2D (Simard 1998), GALFIT (Peng et al. 2002(Peng et al. , 2010, BUDDA (de Souza et al. 2004) and IMFIT (Erwin 2015), to fit 2-dimensional analytic functions directly onto the image. While there are pros and cons concerning whether 1-D profile or 2-D image analysis is more adequate (e.g. Ciambur 2016), we favour the use of the above-mentioned 1-D radial brightness profiles. As shown in van der Marel et al. (1990) and Läsker et al. (2014), some triaxial spheroids exhibit twisting in their position angle ( ) and ellipticity profile. The aforementioned 2-D codes can not account for the shape changes because they assume a fixed and for each galaxy component. Furthermore,  conducted an empirical comparison between 1-D and 2-D methods on the same 14 For further information on the galaxy size, see Appendix C. , where is the standard deviation and is the number of pixels comprising the isophotal ring at a given radius. The normalized radius ( / max ) is defined by the maximum radius max to which we modelled each galaxy. We display the average ratio with black points. galaxies and remarked on some key advantages of the 1-D approach. With the 1-D radial profiles, weak features such as embedded discs and small bars were more visible. The changes in the , , and harmonic coefficient profiles can signify the presence of these components. They also found that the 2-D routine has a substantially harder time converging to a meaningful solution. We, therefore, adopted the more manually-intensive, but also more reliable, 1-D approach over the automated 2-D approach.
A few of the galaxies from our three bins had to be excluded from the analysis due to various reasons: (i) edge-on galaxies with substantial dust covering the centre of the galaxy (see Section 3.2); (ii) galaxy pairs where the gravitational interaction significantly distorts the structure of both galaxies; or (iii) unique cases for which the fitting programme failed to produce a convincing 2-D model.
Ultimately, we have reliable data for 27, 24, and 52 galaxies in Bins 1, 2, and 3, highlighted as red, blue, and black in Figure 3, respectively. 12 galaxies are rejected in our analysis.

Multi-component decompositions
In this section, we present the justification and describe the decomposition process. The projected radial intensity profile ( ) of a spheroid or bulge can be described by the Sérsic (1968) function, expressed in Graham & Driver (2005) as: where e is the effective 'half-light' radius, is the Sérsic index, e is the surface brightness at e , and can be obtained by solving Γ(2 ) = 2 (2 , ), involving the complete and incomplete Gamma functions. The surface brightness profile is calculated by substituting ( ) into the following equation: where is the zero-point magnitude and is the pixel angular size.
Additional substructures, such as ansae and rings, are common. They have distinct features and stellar orbits from a bulge or a disc. Furthermore, spheroids and discs can have a more complex structure than the Sérsic and exponential function permit. For instance, massive galaxies often have a depleted core (King 1978;Trujillo et al. 2004), where the light in the centre is fainter than what a Sérsic function predicts.
Discs also often exhibit bending in their outer region (van der Kruit 1987; Pohlen et al. 2004;Erwin et al. 2005), classified by the different behaviour as either Type I (no bend), II (a downward-bend), and III (a temporary upward-bend). Moreover, substructures such as nuclear discs, intermediate-scale discs, and nuclear stellar clusters also exist. They can contribute an indelible amount of stellar light to the galaxy.   (1877); Sellwood & Wilkinson (1993) For years, these extra features have motivated the use of multicomponent decompositions for well-resolved galaxies (e.g., Martin 1995;Prieto et al. 1997;Aguerri et al. 1998;Graham et al. 2003b;Laurikainen et al. 2005Laurikainen et al. , 2010Vika et al. 2012;Läsker et al. 2014;Davis et al. 2019;Sahu et al. 2019). Indeed, Stone et al. (2021) concluded 'that two-component fits (e.g., Sérsic plus exponential) are insufficient to describe late-type galaxies with high fidelity'.
In addition, one need to be cautious while modelling with Sérsic function. A part of its usefulness comes from the function's flexibility. It is able to fit a range of light profile shapes. The flexible nature, however, also presents a problem. If one ignores a prominent substructure, such as a bar or an outer ring, the fitting code will attempt to compensate by bending the Sérsic function, intended to describe the bulge, upward at the tail end, leading to a higher Sérsic index and radius e than is correct for the bulge. Therefore, in a counter-intuitive sense, the most accurate way to model a spheroid is to not solely focus on the Sérsic component but also account for potentially biasing substructures. With that in mind, we proceed to model galaxies' structures. . From top to bottom panel: -band surface brightness ( ); (model − data) residual (Δ ); ellpticity ( ); Position Angle ( ); and 4 th -order Fourier harmonic coefficient ( 4 ). The galaxy components are the Bulge (dashed-dotted red line), nuclear lens (dashed red line) , Bar (orange), Type III truncated disc (blue), plus a ansae, a faint inner, and outer spiral arm (cyan).
We use the programme Profiler (Ciambur 2016) to perform the decomposition. The analytic functions in Profiler that were used to depict the different galactic structures are listed in Table 3. We start our fitting using ISOFIT's major axis surface brightness profile, where the galaxy components often appear more prominent than in comparison with the geometric-mean axis 15 . The key to the decomposition is to identify deviations in the various 1-D profiles for each galaxy, such as their surface brightness, position angle, and ellipticity profiles. Each component that we include in the fitting process has a corresponding physical substructure that is visible and identified in either the 2-D image or (more readily) the 1-D profiles. We additionally spent some time scouring the literature for additional evidence to justify the inclusion of components. This typically included confirmation of bars already detected by others or rotating discs seen in kinematic data.
To illustrate the process, here we present a case study using the late-type galaxy NGC 4045. In Figure 5, we present the ISOFIT modelling result. In the bottom row, from left to right, they are the 15 The geometric mean of each isophotes' major-and minor-axis, equivalent to a circularised version of the galaxy. image of NGC 4045, the model produced by the ISOFIT model, and the residual by subtracting the image by the model. The top row illustrates a 1-D representation of the image, model and residual, by averaging the pixel value along the vertical direction. One can see the ISOFIT model has done a marvellous job capturing the surface brightness profile of the galaxy. Figure 6 shows the decomposition for NGC 4045 along its geometric-mean axis, a.k.a. equivalent axis, eq . The galaxy has a small bulge (dashed-dotted red curve), a large extended type III truncated disc (blue curve), and a bar (orange curve, see Erwin 2005;Erwin & Debattista 2013;Díaz-García et al. 2016;Font et al. 2017Font et al. , 2019. Along the and profiles, there are two 'bumps' at eq ∼ 3 and 7 , corresponding to the inner spiral arm-like structures. Comerón et al. (2014) classified NGC 4045 as a (R 1 L)SAB(rs, nl)ab galaxy containing two outer rings and a 'nuclear lens' 16 . Upon closer inspection (see Figure 7), there are two dust lanes obscuring the light at ∼ 7-8 (see also Erwin & Sparke 2002), resulting in a dip in , making the isophote appear more circular. We modelled the nuclear lens with a Sérsic function (dashed red line) with a very low Sérsic index ( ∼ 0.01). The strong position angle ( ) twist at eq ∼ 27 coincides with the small bump in . This feature associated with the bar and elevation in is an 'ansae', where the bar transitions into the disc rotation. Finally, the outermost plateau of ∼ 0.4 beyond eq ∼ 40 is the projection of the inclined, rotationally-supported disc (see Erwin et al. 2005Erwin et al. , 2008. We use an type III truncated exponential model for this disc, plus a Gaussian function to account for the brightening in at eq ∼ 43 , which thereby captures the weak spiral arm within this disc. Failing to account for the elongated bar structure in NGC 4045 would result in the fitted Sérsic function substantially over-estimating the mass and size of the spheroidal component of this galaxy. As shown in Column (9) of Table C1, our sample has a range of morphologies.
Our decomposition experience also reveals that the outer components can significantly affect the estimation of a bulge. That is, if a prominent structure in the outer region of a galaxy, such as a strong spiral arm or disc truncation, is not considered, then the parameters of the disc will be distorted. Like a domino effect, an incorrect disc model impacts upon the Sérsic model at the centre, changing its parameters to compensate for the ill-modelled disc.
We analysed 103 galaxies with the same individual level of attention to detail, following that done by Davis et al. (2019) and Sahu et al. (2019) in their decomposition of late-and early-type galaxies with directly measured supermassive black hole masses. Tables D1-D3 provide the best-fitting, equivalent-axis, Sérsic parameters for the 103 spheroids. This axis provides circularised component sizes which enable Profiler to integrate each component's model surface brightness profile to obtain the associated luminosity. The resulting spheroid apparent magnitudes are given in Column (7) of the tables. From there, we calculated the stellar mass based on the four MLCRs from equation array 3. These are listed in Columns (9)-(12).

Morphological reclassifications
We are able to reclassify the morphology of the host galaxies based on our decomposition of their light. The previous morphologies were determined by visual inspection of photographic images (Dressler 1980;Binggeli et al. 1987;Sandage et al. 1985;de Vaucouleurs et al. 1991), which can drown out features and limit the accuracy of the  Figure 5. One can see a four armed 'pinwheel' shaped pattern, aka the nuclear lens shown as dashed red line in Figure 6. The two bright fringes of dust lanes spiraling towards the centre obscure some of the light and decreases the ellipticity at ∼ 7 , making the isophote appears more circular than it should be. past classifications. For instance, a faint disc in an early-type galaxy can be hard or impossible to see in images at certain fixed contrasts. In the case of a face-on disc, from a photograph, it is very difficult to distinguish a barless lenticular galaxy from an elliptical galaxy. For this reason, in our search for compact massive spheroids, we have not restricted the initial galaxy sample to galaxies classified as either lenticular or spiral.
Many, and no doubt most, galaxies classified as elliptical are either lenticular or ellicular galaxies with discs and bulges. Analysing the surface brightness profiles is much more reliable for identifying host galaxy substructures due to the bumps and bends in the profiles. Our decompositions, therefore, provide a more complete picture of the host galaxies than past visual inspections. The new information will be particularly valuable to understand the morphological evolution of galaxies.
We present one such 'elliptical galaxy', NGC 4008, to illustrate why we reclassified it as an S0 galaxy instead of an E galaxy. NGC 4008 was classified as an 'E5' in RC3. Figure 8 shows the galaxy, our model, and our residual image of NGC 4008. The galaxy has a central stellar velocity dispersion of 0 ∼ 208-240 km s −1 (Faber et al. 1989;Simien & Prugniel 2002). In its ellipticity profile ( , middle panel in Figure 9), one can see that the ellipticity reaches ∼ 0.4, and the surface brightness profile ( , top panel in Figure 9) has a smooth exponential decline from eq ∼ 20 . Although the and 4 seem to be featureless, the galaxy has clear signs of an extended disc. Simien & Prugniel (2002) found the stellar velocity to be = 123 km s −1 at = 10 (in major axis) and appear to keep rising beyond this point. We first tried to use a Sérsic bulge plus an exponential disc to model the galaxy, but we found that it can be better described with a bulge (top panel in Figure 9, red line) and a Type III exponential disc with a slight slope difference separated at eq = 35 (top panel in Figure 9, blue line). The final fit is an excellent agreement with the isophotal model and has a root mean square residual Δ rms = 0.0516 mag arcsec −2 (see the middle panel in Figure 8). This decomposition also aligns with Simien & Prugniel (2002), who show that at eq > 10 , NGC 4008 transitions from bulge-dominant to disc-dominant. The presence of an extended disc means NGC 4008 is not an 'E5' but an 'SA0' galaxy. There are 17 other galaxies similar to NGC 4008 where we reclassified them as S0.
Two other galaxies are worthy of note. NGC 4459 was the only galaxy originally labelled as an S0 prior to our changing it to an E2 designation. It contains a nuclear disc in the centre, as reported by Peterson (1978) and Krajnović et al. (2011). Gutiérrez et al. (2011) has suggested that it contains a Type III anti-truncated disc.  broke down the galaxy into a Sérsic bulge, a Gaussian nucleus (for the nuclear disc), and an exponential disc. However, the neighbouring disc galaxy is likely the main contribution to the light in the 'extended disc' at large radii. As we limit the fitting range to within 77 along the major axis, a nuclear exponential disc and a Sérsic bulge are found to suffice. Additionally, we sustain the E classification for NGC 2872, but it is noteworthy that it contains a dusty nuclear disc (as shown in Tran et al. 2001).
This result is comparable to the ATLAS 3D kinematic analysis of 260 nearby early-type galaxies (Cappellari et al. 2011a;Krajnović et al. 2011;Emsellem et al. 2011;Cappellari et al. 2011b;Krajnović et al. 2013a). Studying a volume-limited sample within <40 Mpc, Emsellem et al. (2011) reveal only ∼14% (36/260) of their early type galaxies are 'slow-rotators', defined by a luminosity-weighted specific angular momentum R 0.1 (see Emsellem et al. 2007). They found that one-third (∼34%) of the previously-classified elliptical galaxies are slow rotators, and ∼47% (17/36) of them contain a kinematically-distinct core (KDC). We share a mutual sentiment that reclassification is needed for many early-type galaxies. Indeed, Cappellari et al. (2011b) advocated the kinematic classification scheme from Bender (1988a) and Capaccioli & Caon (1992) to correct for the misunderstanding caused by the old visual-based classification. Although, as illustrated in Bellstedt et al. (2017), this kinematic scheme fails to accommodate the ES galaxies which are either fast or slow rotators depending on the radial extent of one's observations, it does help identify which early-type galaxies have discs and therefore require a decomposition to measure the size and mass of the spheroidal component. 17 Among the true ellipticals, there are six cD galaxies: NGC 2832, NGC 4073, NGC 4874, NGC 4889, NGC 4914, and NGC 5444. NGC 4914 is relabelled as a SAB0 galaxy, for its weak bar  and an extended disc with ∼ 0.4. Informed by the three components model by Dullo (2019), NGC 4874 is now labelled as SA0 (see later in Section 4.1).

Spheroid colours
Local early-type galaxies typically, although not universally, exhibit a negative colour gradient such that they are redder in their centre than their outskirts (e.g. Franx et al. 1989a;Peletier et al. 1990;Gargiulo et al. 2012;Peletier et al. 2012;Marian et al. 2018). In the absence of substantial dust and gas, a negative gradient suggests younger (more blue) stellar population in the outer regions, which can be a sign of inside-out disc growth. However, from galaxy colour profiles, it can be problematic to readily establish the colours of the individual overlapping structural components. For instance, a young nuclear disc can make a galaxy core colour blue and result in a more positive galaxy gradient than compared to the surrounding red quiescent spheroid (e.g. Graham et al. 2017).
While we would like individual − colours for the bulges, or better yet, the spectral energy distributions of the bulges for use in software like MAGPHYS (da Cunha et al. 2012) to obtain the stellar mass, we make use of the galaxy colours, recognising that there is not a big difference between bulge and disc colours of early-type galaxies and early-type disk galaxies (e.g., Peletier & Balcells 1996a). Our galaxies have a mean − colour of ∼1.2, in accord with the value reported in Fukugita et al. (1996) for what they thought were elliptical galaxies, but undoubtedly included many S0 galaxies.
In general, the bulges of galaxies are redder than their discs (Bothun & Gregg 1990;Peletier et al. 1999). For bright red galaxies, Kennedy et al. (2016, see their Figures 10-13) find − bulge colours ∼0.5 mag redder than the disc colours, although they remind the reader that they performed a bulge-disc decomposition of every galaxy regardless of whether there was a physical need for two components. As such, their colour difference will also capture metallicty gradients in disk-less elliptical galaxies. However, for some of our galaxies with bulges, their galaxy − colour might be bluer (a smaller value) than the bulge colour.
Our use of the galaxy colour in the MLCRs to derive the stellar mass of the bulge might, therefore, underestimate the stellar mass of the bulges. Although, Vika et al. (2014, see their Figure 10 and Table 2) report that E-S0-Sa galaxies tend to have bulges and discs with a similar ( − ) colour (see also Peletier & Balcells 1996b). Using morphologically-classified galaxies, they find it is the Sb and later spiral galaxies that have notably bluer discs than their bulges, with average ( − ) disc colours that are 0.3 mag bluer than their bulges. Throughout this work, we use the galaxy colour as a proxy for the bulge/spheroid colour. The effect of underestimating how red any bulge is will be to underestimate the * / ratio assigned to that bulge and thus underestimate the stellar mass of that bulge. That is, the number density of local, massive, compact spheroids will not be erroneously overestimated due to this issue.
In an effort to check on this matter, in Figure 10, we compare the galaxy ( − ) gal colour with the SDSS point-spread function (PSF) ( − ) PSF colour. Figure 10 reveals how some galaxies with a global ( − ) gal ∼ 1.2 ± 0.1 mag have relatively blue cores. The PSF magnitude is obtained by fitting a point-spread function to the light source, stars, and galaxies alike (see in Stoughton et al. 2002). A PSF is a good description for point sources like stars, but not so for extended sources like galaxies. The PSF fit on the galaxies should mainly sample the spheroidal part of the galaxies. PSF size variation is biased to the red. We highlighted those with unique central components in Figure 10. The sample with a nuclear disc is marked with blue crosses. Other prominent central components, such as nuclear rings, secondary bar, and AGN are marked with cyan crosses.
All of this simply provides a rough guide for what the error in the spheroid colour may be. Of the 103 galaxies in our sample, 69 have PSF magnitude available in the NED database. The median difference between ( − ) gal and ( − ) PSF is 0.09 mag and the standard deviation ±0.26 mag. Most galaxies reside at 1.0 < ( − ) gal < 1.5 along the dashed centre line. Some of the major outliers can be explained by the influence of prominent central components, as evidenced by the four galaxies at ( − ) PSF < 0.9 mag and the one with ( − ) PSF > 1.8 mag. One outlier deviates from the centre line significantly at ( − ) gal ∼ 1.8 mag and ( − ) PSF ∼ 0.8 mag. This galaxy is UGC 8736 (marked with a red cross in Figure 10). Its image has severe foreground contamination. The galaxy is situated next to a bright star. We suspect its PSF colour might not be accurate due to such an influence. The galaxy ( − ) gal colour provides a reasonable estimation of the spheroid ( − ) for our purpose.
Given that galaxies do not have − colours redder than ∼1.5 (e.g., T11), we recognise the five galaxies 18 in our sample with a − colour of ∼1.8 as erroneous, perhaps due to dust or a poor single Sérsic fit to obtain the − and − 0band magnitudes. We reset these galaxies − colour to 1.5 so as to not overestimate the stellar mass-to-light ratio used and thus not overestimate the galaxies' stellar mass. We move on to calculate the stellar masses of the spheroids accordingly.

Spheroid masses
Multiple high-red nugget studies (Barro et al. 2013;van der Wel et al. 2014;van Dokkum et al. 2015) use the SED fitting code FAST ) to determine the stellar masses of their highgalaxies. These works assume the Chabrier (2003) IMF and Bruzual & Charlot (2003, hereafter BC03) SPS, with a variety of dust extinction treatments. Among the available MLCRs, only the T11 and RC15 relations are constructed with these priors. For the sake of simplicity, the result will be presented in RC15 stellar mass. However, to ensure the result is not biased by the choice of stellar mass measurements, the effect of using the different MLCRs will be later discussed in Section 6.1.2.
Our final uncertainty on the stellar masses have been calculated 18 The five outliers are: UGC 8736, NGC 4527 NGC 3718, NGC 2968, and NGC 2894. All of which are spiral galaxies from Bin 3. We marked their spheroids in Figure 15 with black star signs (★). Three of them contains a low mass spheroids ( * /M < 1 × 10 10 ). Only two lie within a 'compact massive' region. The outliers does not alter the number of compact massive spheroids, nor the shape of the spheroids distribution in any significant way. The galaxies that contain a nuclear disc are marked with blue crosses. Those containing other prominent nuclear components, such as a nuclear ring, secondary bar, or AGN, are marked collectively with cyan crosses, labelled as 'nuclear cpt'. The particular outlier UGC 8737 is marked with a red cross.
using the following error propagation equation: where is our assumed uncertainty in the distance, is our uncertainty in the bulge apparent magnitude, and Υ * is our uncertainty in the * / ratio. The uncertainty in the magnitude, , is determined by the component modelling process. The unknown degeneracy amongst components is the primary source of error. As explored in , adding or subtracting a galaxy component can change the spheroid parameters. Based on our extensive practical experience with the fitting routine, and informed by the results of previous studies Davis et al. 2019;Sahu et al. 2019), we have assigned a probable magnitude error of = 0.3 mag for all of our spheroids. This level of uncertainty encapsulates the average error commonly encountered throughout the combined efforts of this work and its predecessors. The uncertainty in the mass-to-light ratio is taken to be the intrinsic scatter of each MLCR, as recorded in Section 2.3. The error in distance is taken directly from their respective measurement sources (see Appendix A).

Bulge sizes
We present the bulge/spheroid size (in equivalent axis, e,Sph ) distribution in Figure 11. The spheroids from Bin 1, 2, and 3 are depicted with red, blue, and black filled histograms, respectively. The median size e,Sph for Bin 1, 2, and 3 are 1.7 kpc, 1.2 kpc, and 0.8 kpc, respectively. The majority of our spheroids (96/103, 93%) are smaller than e,Sph < 7 kpc.
Over time, as the stars in galaxies evolve and eject their stellar winds, the amount of dust builds up, at least in late-type galaxies where dust sputtering does not operate efficiently. While the presence of dust is known to reduce the apparent bulge and disc magnitudes, and increasingly so with a more edge-on disc inclination ( ), the impact on the effective half-light radius of the bulge is not well known. An investigation of how e changes with , at different wavelengths, Figure 11. The size ( e,Sph , in equivalent axis) distribution of our bulges/spheroids. The spheroids from Bins 1, 2, and 3 are depicted by the red, blue, and black histograms, respectively. The thick black line depicts the size distribution of the entire sample set. The figure only shows the 96 spheroids with e,Sph < 10 kpc of the complete sample, see the spheroids' size-mass distribution later in Section 4.3 for further reference.
would be welcome. Ideally, this information would be known for different galaxy morphological types.
The dust is known to be more centrally concentrated in galaxies than in their outskirts. Qualitatively, this acts to reduce the inner flux relative to the outer flux, which has the effect of increasing the halflight radius of the galaxy, and the disc component of the galaxy (e.g., Graham & Worley 2008, see their Figure 2). It is not known how much the apparent half-light radius of the bulge increases, however, the effect is such that we may be underestimating the true number of local, compact spheroids.
Combined with potentially underestimating the masses of the bulges, we move forward in the knowledge that we are presenting a conservative estimate of the number density of compact massive spheroids. That is, these potential biases do not act to overestimate the number density.

RESULTS
We present our result in the following subsections: Section 4.1 showcase the bulge-to-total flux ratio ( / ) newly obtained from multicomponent decomposition; Section 4.2 shows the bulge/spheroid mass function; and Section 4.3 present the size-mass relation of the local spheroids and the number of compact massive spheroids under different definitions.

Bulge-to-total ( / ) flux ratios
We present the -band bulge-to-galaxy flux ratios in Figure 12, plotted against by their morphologies. Our flux ratios are the observed (no internal dust correction) ratios taken from the equivalent-axis decompositions. For simplicity, we have applied the same * / ratio to all galaxy components; this poses a slight limitation and is an area where future improvements can be made. However, we note that the (dust-free) colours of bulges and discs are, in general, not too dissimilar from each other in massive disc galaxies (Peletier & Balcells 1996a). As such, the * / ratio will typically, although not The spiral (SA-SB) galaxies presented here are dominated by the so-called 'early-type' spirals (Sa-Sb). The boxes denote the lower (Q1) to upper (Q3) quartiles of the data, and the whiskers extend from the minimum to the maximum data points (excluding any outliers). The red lines represent the median / ratios for each morphology. The '+' marker denotes an outlier (NGC 5350) that lies beyond the Q1 − 1.5(Q3 − Q1) limit (Tukey 1977). necessarily, be similar. The grey points indicate the bulge-to-total ( / ) flux ratios for individual galaxies while the red line is the median ratio in each subgroup. We obtained a disc-dominated median / ±1 ≈ 0.33±0.15 for the S0 galaxies (49 in sample). Our spiral galaxies exhibited a notably lower median ratio of ≈ 0.13 ± 0.12, (41 in sample). In both cases, the barred galaxies consistently have a lower / ratio than the unbarred galaxies.
In Table 4, we present the median / ± 1 value for each morphological type from each bin. The two unbarred ES galaxies have an average / ± 1 ≈ 0.86 ± 0.02. While this simply reflects the range of masses and galaxy types, it is important for understanding the transition from galaxy to spheroid mass at the high-mass end of the galaxy stellar-mass function.
Our inclusion of minor components, such as nuclear discs, spiral arms, lenses, and secondary bars can change the / ratio. This is not so much because of the flux held in these components but because of how their presence can skew the galaxy decomposition if they are not properly accounted. Although one may expect this to decrease the / ratio, it is more complicated than that. For example, the more sophisticated models, which can involve truncated-exponential discs rather than a single exponential disc model, can result in a fainter central disc surface brightness and thus a greater assignment of flux to the bulge if the single-disc model was biased by the outer disc profile. We have 32 truncated discs in our sample: 14 Type II (downwardbending) and 18 Type III (upward-bending) discs. Echoing the same point made by Kim et al. (2014b), which claimed that ignoring Type II disc breaks will result in an ≈10% decrease in the / ratio, we agree that broken, and inclined, exponential disc models play an important role in correctly modelling a galaxy. Contrarily, the application of these more sophisticated disc models to Type III discs can be expected to decrease the / ratio, compared to the result using a single exponential model, if the single-disc model is biased by the outer disc slope.
The lack of a strong correlation between / ratio and galaxy stellar mass is reminiscent to the result from Méndez-Abreu et al. (2017), where they performed detailed 2D multi-component decomposition (including nuclear, bar, and broken disc components) to 404 galaxies from the Calar Alto Legacy Integral Field Area (CALIFA, Sánchez et al. (2016)) data. They presented the mean / ± 1 value for S0, Sa, and Sb to be ∼ 0.32 ± 0.17, 0.28 ± 0.17, 0.12 ± 0.11 respectively. The mass of their S0, Sa, and Sb galaxies occupy a similar range of log 10 ( * /M ) ± 1 ∼ 10.79 ± 0.6, 10.80 ± 0.4, and 10.48 ± 0.4, respectively. It has been shown that, if the bar component is not considered, the bulge's Sérsic index and the B/T ratio will be overestimated (Gadotti 2009;Salo et al. 2015). The rather low median / ratio ( / ∼ 0.1-0.3) in both massive S0 and S galaxies reflect a reality that stellar discs and their induced structures take up the majority of the mass budget (see also, Laurikainen et al. 2005).

Spheroid stellar mass function
We present the high-mass end of the local bulge/spheroid stellar mass function in Figure 14. It includes bulges in the mass range 2 × 10 9 * /M (RC15) 1.5 × 10 12 . To our knowledge, these position in Dullo (2019), it contains three components, an inner core-Sérsic bulge, an intermediate component, and an outer exponential halo. The intermediate component is dominant between ∼ 10 and ∼ 100 . In their 606 − 814 colour map (their figure 2), the colour becomes gradually bluer towards larger radii. We found that an exponential component can describe the intermediate component well. The outer envelope comes into play beyond our cut-off radius of 90 in the major axis and was therefore not required in our decomposition 20 Some (13/103) galaxies are below our lower mass selection limit in ( * /M (RC15) = 6.7 × 10 10 ). This is somewhat expected because SDSS −band total galaxy magnitude is measured using a single Sérsic function, while ours is measured by a multi-components model (Bulge+Disc+others). For galaxies with prominent discs (S0 & S), a single Sérsic model can overestimate the size and luminosity of the galaxy (see later in Figure C1). The implication of which, however, shall be better explored in future works. results are the first of their kind derived from highly-detailed, multicomponent decompositions on a volume-and mass-limited sample. For each Bin (1 to 3), the mass function is constructed by counting the number of spheroids per 0.3 dex range in mass, then dividing that number by 0.3 and by the volume of the associated Bin (1 to 3).
Based on the / flux ratios from Section 4.1, we use the individual galaxy / ratios in each bin to estimate the median / ± 1 range for each bin. They are: (i) Bin 1 ( / ∼ 0.33 ± 0.32): E (5), EAS (1), SA0 (13), SB0 (2), SA (1), SAB (3), and SB (2) While we have a well-defined stellar mass 'selection limit' for the galaxies, it is unclear to what completeness the spheroid mass function is covered below this limit. Spheroids/bulges more massive than the 'selection limit' cutoff are of course accounted for, but since we did not sample galaxies below the limit, the spheroid mass function is only 'partially complete' at lower masses. For instance, a massive galaxy of * ,gal ∼ 2 × 10 11 M with / = 0.1 will have a bulge of * ,Sph ∼ 2 × 10 10 M , below our selection limit. Our sample provides a good representation for embedded spheroids coming from galaxies more massive than 6.7 × 10 10 M (RC15), but increasingly less so for less massive spheroids. For example, a less massive galaxy ( * ,gal ∼ 4 × 10 10 M , outside of the limit) with a bulge-to-total ratio / ∼ 0.5 will also have a spheroid of the same stellar mass as above ( * ,Sph ∼ 2 × 10 10 M ), but it will not be in our sample. We refer to the region below the 'selection limit' as the 'partially-complete region'. The number of spheroids found in this mass range represents a lower limit to the true number.
We estimate a rough lower bound to this 'partially-complete region' by multiplying the galaxy mass limit (see Table 2) by the average / ratio for the galaxies with masses greater than this limit. This provides the following rough lower bound to the spheroid stellar masses (when using the / prescription from RC15) for each bin: 1.12 × 10 11 M (Bin 1), 3.77 × 10 10 M (Bin 2), and 1.54 × 10 10 M (Bin 3). Between the galaxy mass limit and these lower values, our spheroid mass function is only partially complete. Above the galaxy mass limit, our numbers are complete. The galaxy mass 'selection limit' has been depicted by the colour-coded solid bars, and the 'partially-complete region' by the semi-transparent colour-coded bars, in the lower right of Figure 14. The data points above the 'selection limit' are shown with solid circles, while those in the 'partially-complete region' are shown using transparent halffilled circles, while those less massive than that are depicted using open circles.
We present the overall mass function across the three bins using a dashed cyan line in Figure 14, obtained by summing up the number of spheroids in each mass interval from all three bins and dividing them by the total survey volume (4.76 × 10 5 Mpc 3 ) and then by 0.3 dex. The purple dashed line depicts the overall mass function of small spheroids with e < 2 kpc (equivalent axis). The overall mass function serves as a lower limit of spheroids in each mass interval. At stellar masses above * ,Sph /M (RC15) ∼ 1.5 × 10 11 , there are no spheroids smaller than 2 kpc in our sample. However, there is no strict definition in the literature for 'compact massive' galaxies. To better explore the situation requires us to next look at the size-mass diagram.

Size-mass relation of the local spheroids
In the top left panel of Figure 15, we show the distribution of local spheroids in terms of their (equivalent-axis) effective half-light radius e,eq and their stellar masses derived from the RC15 MLCR, assuming a BC03 SPS. The spheroids are colour-coded according to the volume-bins (1, 2, and 3) shown in red, blue, and black, respectively. We plot our spheroids against the opaque grey cloud of the general population of bright SDSS galaxies, consisting of both early-and late-type galaxies. It is clear that, for a given mass, the local spheroids are considerably more compact than most galaxies. As expected, the spheroids extracted from Bin 3, with their lower galaxy mass-boundary, occupy the lower to the medium mass range, and the giant spheroidal elliptical galaxies reside in the (high mass)-(large radius) region of the diagram.
In the other panels, we report the number density ( (Mpc −3 )) based on different selection criteria. The arbitrary selection boundary often contains a diagonal line separating the sample into 'compact' and 'non-compact' systems, with a lower-mass limit.
The compact massive quiescent galaxies in the CANDELs survey (Grogin et al. 2011;Koekemoer et al. 2011) within this boundary roughly follow the Newman et al. (2012) size-mass relation. The slope of the selection boundary is motivated by the apparent trend of the quiescent galaxies (with specific star-formation rate log 10 ( UV+IR ) > −0.5 [Gyr −1 ]) at > 1.5.
Note, that unlike most studies, van der Wel et al. (2014) selected their compact sample by the major axis effective radius e,major instead of in the equivalent axis e,eq , where the two quantities follow the relation e,eq ≡ e,major √︁ / , in which / is the semiminor to semi-major axes ratio of a galaxy's apparent elliptical shape. Therefore, the size-mass plot in the lower left of Figure 4.3 is plotted in major axis instead.
• (13) Figure 14. The local bulge/spheroid 'stellar mass' function, with the stellar mass defined using the RC15 MLCR. Red, blue, and black points represent the spheroid sample from the three 'galaxy mass bins' 1, 2, and 3, respectively (see Figure 3). The spheroid stellar masses are divided into 0.3 dex mass bins. The error bars represent the Poisson error. The cyan dashed line, labelled ' Sph ', is the lower limit of the number density for each mass interval, obtained by summing up the total number of spheroids in each interval and dividing it by the total volume of our selection space. The purple dashed line, labelled '< 2 kpc', is the same as the cyan line but including only the spheroids with e < 2 kpc (in equivalent axis). The three colour-coded bars in the bottom right represent the mass range of confidence, in which the solid bars depict the mass 'selection limit' from the host galaxies, and the transparent bars depict a 'partially-complete' region. The points in the 'selection limit' are plotted in solid circles, the ones in the 'partially-complete' region are in half-filled circles, and the ones below these regions are in open circles. The spheroid mass function is complete down to each host galaxies' 'selection limit' where all the hidden spheroids within this range are accounted. In the 'partially-complete' region, our sample includes all spheroids coming from high-mass galaxies ( * /M 0.7-3.0 × 10 11 ) with low / ratios, but not the spheroids coming from low-mass galaxies ( * /M ∼ 0.18-1.5 × 10 11 ) with high / ratios. Since Bin 1 has a volume of 3.25 × 10 5 Mpc 3 , 1 object per 0.3 dex range of spheroid mass will give a number density of Φ ≈ 1 × 10 −5 Mpc −3 dex −1 . See Table 1 for the volume of each bin.
Unlike the high-studies, Damjanov et al. (2014) applied this selection criteria on a sample of intermediate redshift galaxies from the COSMOS field at 0.2 < < 0.8. While having discussed the different choice in the mass limit, we set it to the lowest limit for easier comparison with Barro et al. (2013).
• In the middle right-hand panel, we select for the van Dokkum et al.
This selection criteria are constructed based on the argument that the Barro et al. (2013) boundary is not restrictive enough because it will also include 60% of the star-forming galaxies, which have a factor of two bigger effective radii than the quiescent galaxies. It could enhance the 'progenitor bias', where the quenching of starforming compact galaxies may inflate the number density of the passively evolving quiescent galaxies. It does not, however, concern our sample because the stellar population of local bulges has exist for a long time (MacArthur et al. 2009). It is unlikely the quenching process between 0 < < 1.0 to produce such system.
Building on Graham (2013), Graham et al. (2015) furthered the proof of concept for this study as it supplied generic selection criteria on small, but massive bulges in the local Universe from different studies (Seigar & James 1998;Graham 2001;Möllenhoff & Heidt 2001;Balcells et al. 2007;Laurikainen et al. 2010;van den Bosch et al. 2012;Dullo & Graham 2013;. They found 21 small, but massive spheroids within 90 Mpc and gave an initial (lower) estimate of the volume number density ∼ 6.9×10 −6 Mpc −3 (or per unit dex stellar mass, it is 3.5 × 10 −5 Mpc −3 dex −1 ) as a lower limit for such systems.
While we have focused on galaxies more massive than 10 11 M , given that the presence of a disc results in a bulge-to-total ratio less than 1, we have uncovered some spheroids less massive than 10 11 M . Although our spheroid sample is incomplete below this mass, we are able to obtain lower-limits to the actual number density of compact spheroids defined by a range of selection criteria in the literature which have encompassed masses down to 1, 4, 5, and 7×10 10 M . It must, therefore, be kept in mind that the number densities, from our three bins, shown in Figure 15 are lower-limits.
The majority of the spheroids (66/103) lie within the size-mass range 0.4 kpc < e < 6.0 kpc and 1.0 × 10 10 < * /M (RC15) < 2.5×10 11 . The lower bound simply reflects the mass boundary of our host galaxy sample and reveals that some galaxies have / ratios less than 0.1. The scatter along either the size or mass axis is rather small in this region. With (geometric mean)-axis half light radii less than ≈1 kpc and stellar masses * ,sph < 2 × 10 10 M (RC15), the scatter is much more prominent. This variation of scatter is less apparent when using the major-axis size, as shown in the lower-left panel of Figure 15. We include this because van der Wel et al. (2014) used the major-axis galaxy sizes for their investigation.
Considerable effort has been made searching for the ultra-compact massive galaxies (UCMGs) in the local Universe, defined by as those with * ,gal /M > 8 × 10 10 and e < 1.5 kpc (Trujillo et al. 2009;Ferré-Mateu et al. 2017;Tortora et al. 2018Tortora et al. , 2020Scognamiglio et al. 2020). These objects are extremely rare compared to the regular red nuggets, with the number density of UCMG ∼ 1-9 × 10 −6 Mpc −3 at < 0.5 (Tortora et al. 2018(Tortora et al. , 2020. In our sample, there are two spheroids that satisfy the UGMG definition, giving a lower limit to the number density of ultra-compact massive spheroids of c,Sph ∼ 4.2 × 10 −6 Mpc −3 . We provide an assortment of information in Tables D1-D3 regarding our spheroid sample. The parameters of the Sérsic functions are listed in Columns (4)-(6). Column (7) gives the spheroid apparent magnitudes as measured by Profiler. Column (8) lists the absolute magnitudes of the spheroids, based on the distances shown in Appendix A. The magnitudes are corrected for galactic extinction from Schlegel et al. (1997). Columns (9)-(12) tabulate the spheroid stellar masses, obtained using each of the four MLCRs (Equation Array 3).

COSMIC EVOLUTION OF COMPACT MASSIVE SPHEROIDS
In this section, we compare our local spheroids with the red nuggets in low-to-intermediate-(Section 5.3) and high-(Section 5.2) from the literature, as well as the embedded local bulges obtained by other multi-decomposition works (Section 5.1). Subsequently, we are able to plot out the potential evolution of compact massive spheroids across cosmic time. The effect of adopting a different * / ratio is explored in Section 5.4.

Low redshifts ( < 0.03)
We have compared our local spheroids with the bulges from other studies of ≈ 0 galaxies. We have a preference for decompositions with a similar degree of detail. Relevant studies of the galaxy size evolution (e.g., Trujillo et al. 2006Trujillo et al. , 2007Taylor et al. 2010;Newman et al. 2012;McLure et al. 2013;Damjanov et al. 2014;Fang et al. 2015) have measured the effective radii and stellar masses of local early-type galaxies, but not their bulges nor the bulges of any massive early-type spiral galaxies. This may have prevented them from establishing the evolutionary pathway taken by the high-compact massive galaxies. The underlying assumption of those works was the preservation of the morphological types throughout their evolutionary history. Due to our different approach, a comparison with the galaxy sizes and masses from those studies would be of limited value to test our theory. We can, however, explore if our bulge sizes and masses agree with those from multi-component decompositions of local galaxies.
Several studies of supermassive black hole scaling relations (e.g., Davis et al. 2019;Sahu et al. 2019) have performed a similar style of careful galaxy analysis using IRAC 3.6 µm images of galaxies in the local Universe. Figure 16 shows the size-mass distribution of our spheroids (red points) and the local bulges and giant elliptical galaxies from these studies.  obtained their surface brightness profiles via IRAF task ellipse from 66 early-type host galaxies and measured the sizes and masses by using a bespoke fitting code. Davis et al. (2019) and Sahu et al. (2019) extracted their surface brightness profiles with ISOFIT and performed decompositions via Profiler. The Davis et al. (2019) spheroids are from 44 spiral host galaxies, while the Sahu et al. (2019) spheroids are from 40 early-type galaxies. These studies examine each galaxy on a case-by-case basis and assign physical components (e.g., bars and ansae rather than random Sérsic components) when modelling each surface brightness profile (see the appendices of the respective papers for detailed discussions on each galaxy).
Despite the decompositions being conducted by different personnel, the resulting distribution of bulge sizes and masses appear remarkably consistent. The high-mass end of the near-linear distribution in Figure 16 resembles the bright arm of the spheroid distribution seen in Graham (2013, see Figure 1). Here, we mainly want to illustrate the consistency across the different studies; that is, there is no unusual offset in our sample, nor sign of dichotomy. The implications of this bulge/spheroid size-mass diagram shall be explored in future work.
In passing, we briefly note that we may be seeing evidence for a steepening of the relation at high masses. This is expected for mergerbuilt elliptical galaxies for which the mass scales with 2 eq / yet barely increases over the value in the progenitor galaxies (e.g., Volonteri & Ciotti 2013). Also at play is the influence of the intracluster light (ICL) surrounding BCGs, and we may have over-estimated the galaxy sizes at the top-end. Driver et al. (2007b) produced a -band bulge luminosity function from two-component Sérsic+exponential fits (Allen et al. 2006b) at 0.013 < < 0.18. While that work contained no dust corrections, Figure 2 from that work reveals, for their brightest bulges, a declining number density at = −20 to −20.75 mag (AB magnitude system) -where S0 galaxies likely dominate and there is no substantial dust -of around 3 × 10 −4 to 0.5 × 10 −4 Mpc −3 dex −1 . For , = +5.44 mag, and * / = 8, this roughly corresponds to 1.2-2.4×10 11 M . Their number density, therefore, overlaps well with our result at a similar mass in Figure 14. We do, however, find the same number density at 5×10 11 M in our data, while Driver et al. (2007b) has no bulges at these higher masses.
Our number density drops to ∼0.2×10 −4 Mpc −3 dex −1 at * ,Sph ∼ 10 12 M . The absence of the highest mass bulges in (Driver et al. 2007b) may be due to the 'logic filter' employed by Allen et al. (2006b) for identifying the more secure bulge+disc decompositions. It meant that the ES galaxies, and importantly perhaps some misfit S0 galaxies, remained in the single-component E galaxy bin at these highest masses. Resolving this is, however, beyond the scope of the current investigation. However, it is important to note that in this mass range, there are only two spheroids in our sample. It reflects the reality that spheroids with mass * ,gal /M ∼ 10 12 are rare indeed.

High redshifts (1 < < 3)
The different selection criteria for compact massive systems have been used to illustrate the migration pattern of galaxy distributions, over cosmic time, in the size-mass diagram. As shown by Barro et al. . Additionally, the 13 E+ES class galaxies are highlighted in squares ( ), where E galaxies are in green and ES galaxies are in yellow. The outliers with galaxy colour ( − ) > 1.5 are marked with black star signs (★, see Section 3.5). In our sample, there are 55 (53%), 10 (9.7%), 60 (58.3%), 9 (8.7%), and 7 (6.7%) compact massive spheroids according to Barro et al. (2013) 2015), with decreasing redshift, less and less galaxies occupy the 'compactness' selection boundaries. Damjanov et al. (2015a) demonstrated that while the actual number densities are different -depending on the selection criteria for compact massive systems -the overall trend/reduction with time remains.
We have compared the number density of our local, compact massive spheroids with that of compact massive galaxies in the intermediate and high redshift Universe. Figure 17 shows the evolutionary trend of number density for compact systems across time, based on four different compactness criteria, from top to bottom: Equation 11 (top), 12 (upper-middle), 14 (lower-middle), and 13 (bottom). To represent the high-(1.0 < < 3.3) red nuggets, we have included the number density trends from three works: Barro et al. (2013)  In Figure 17, the three colour-coded circles (•) at = 0 are our lower limits on the number density of compact massive spheroids, as obtained using the different selection criteria applied to our three bins. In all the panels, Bin 3 tends to have the highest number density, and Bin 1 the lowest. This is a result of the bin design. As Bin 3 is complete down to the lowest host galaxy stellar mass ( * /M (IP13) ∼ 1 × 10 11 ) compared to Bin 2 ( * /M (IP13) ∼ 2 × 10 11 ) and Bin 1 ( * /M (IP13) ∼ 5 × 10 11 ), we are able to capture more compact massive spheroids in Bin 3. Note that in the lower-middle panel (van Dokkum et al. (2015) criterion) of Figure 17, the number densities of all three bins are closer than the other panels. This is because van Dokkum et al. (2015)'s criterion selects only for high-mass spheroids and red nuggets with * /M (RC15) 4 × 10 10 . If relatively lowmass (1 × 10 10 * /M (RC15) 4 × 10 10 ) systems are also included, such as when using the selection criteria from Barro et al. (2013) and Damjanov et al. (2015a), an abundance of low-mass compact massive spheroids can be found in Bin 3. These 'red pebbles' are more likely to live inside of galaxies within the mass range of 6.7 × 10 10 < * ,gal /M (RC15) < 1.3 × 10 11 . As a result, there is roughly a factor of 10 more compact massive spheroids in Bin 3 in the upper-and lower-middle panels than the top and bottom panels in Figure 17.
For the sake of discussion, we categorise the number density of our spheroids into c,Sph and E+ES . In Figure 17, across all four panels, the black horizontal line (' c,Sph ') depicts the total number of compact spheroids from all three bins divided by the volume = 4.
These two values depict lower limits to the number densities. For c,Sph , it is the lower limit of compact massive spheroids based on a volume-limited sample of galaxies with * > 6.7×10 10 M (RC15) and within a distance of 110 Mpc. The value for E+ES represents the lower limit to the true number density of ellipitcal (E) and ellicular (ES) galaxies that can be found above * /M ( 15) ∼ 6.7 × 10 10 . Note that E+ES,Bin3 is likely to be close to the true value because true elliptical galaxies tend to be massive; that is, ETGs less massive than 6.7 × 10 10 M tend to be S0 galaxies 21 (Emsellem et al. 2011;Krajnović et al. 2013a;Cappellari et al. 2013;Cappellari 2016). For Bins 1 and 2, we may have missed some E galaxies due to the higher stellar-mass selection criteria required to keep the number of galaxies requiring careful multi-component decomposition at a manageable level.
From Figure 17, it is apparent that the definition of 'compact and massive' makes a difference to the reported number density (see also, Charbonnier et al. 2017). For the less restrictive size-mass sample selection criteria (Barro et al. 2013;Damjanov et al. 2014), the number density of compact massive spheroids is consistently higher than the number density of ellipitcal galaxies ( c,Sph > E+ES ), while  figure 11). The red, blue, and black circles (•), at ∼ 0, represent the number densities from Bins 1, 2, and 3, correspondingly. The error bars of our compact massive spheroids are assumed to be Poissonian. The horizontal black line, labelled ' c,Sph ', is the lower limit of the number densities across the three bins, and the red line, labelled ' E+ES ', is the lower limit for the true elliptical galaxies and ellicular galaxies in the three bins (given by Equation 16). In addition, the Damjanov  applying the more restrictive criteria (van der Wel et al. 2014;van Dokkum et al. 2015), the two numbers are comparable ( c,sph E+ES ). This is, in part, because E+ES does not depend on any size or mass selection criteria, unlike c,sph (see Figure 15).
The discrepancy between the restrictive and the less-restrictive size-mass sample selection criteria in use in the literature is informative. It means that among the local compact systems, the elliptical galaxies and high-mass spheroids (>(4-7)×10 10 M ), formed via an E-to-E process or disc-cloaking, respectively, may have occurred with similar frequency. The disc-cloaking process is likely the dominant mechanism in shaping galaxy evolution within the stellar mass range 1 × 10 10 < * /M (RC15) < 4 × 10 10 .

Low-to-intermediate redshifts (0.03 < < 1)
As a bridge between the local spheroids and the high-galaxies, the information from intermediate redshifts provides insight into the transitioning process. In Figure 17, we compare our results with Damjanov et al. (2015a). They studied the number density of quiescent, compact massive galaxies at 0.2 < < 0.8 from the COSMOS field (Leauthaud et al. 2007). While they acknowledged that most of their sample exhibited a Bulge+Disc structure, they concluded that the single Sérsic profile captures the size of the spheroids well in this redshift range, with perhaps just a slight overestimation in size.
This sample's number density shows consistency with the Barro et al. (2013) and van der Wel et al. (2014) high-sample's peak abundance (see the top and upper middle panels in Figure 17). In the upper panel of Figure 17, the data from Damjanov et al. (2015a) straddles the region between our Bin 2 and Bin 3 number densities. It is also more abundant in comparison to the red nuggets' peaks number densities at = 1.2, and = 1.7. The lower panel shows the trend in number density, with redshift, using the size-mass selection criteria from Damjanov et al. (2014). Unfortunately with this selection, the complete set of data is not available, and Damjanov et al. (2015a) only show a partial result that they selected from the COSMOS quiescent galaxies classified as point sources in the photometric SDSS database. Their Figure 11 depicts only a lower limit to the number density at intermediate-redshift. Nonetheless, our number of compact spheroids is vastly more numerous than the trend shown in the bottom panel. The matching number densities with the high-data led them to conclude that the number of compact objects does not depend strongly on redshift.
Both Valentinuzzi et al. (2010) and Poggianti et al. (2013b) have argued for a constant compact quiescent galaxy number density (see also Saracco et al. 2010). Valentinuzzi et al. (2010) found a substantial number of compact quiescent galaxies in the WIdefield Nearby Galaxy-cluster Survey (WINGS; Fasano et al. 2006) at 0.04 < < 0.07. They reported a minimum number density of 1.3 × 10 −5 Mpc −3 and possibly up to 1.3 × 10 −2 Mpc −3 in cluster environments. Crucially, the majority (67%) of their compact galaxies are S0. In the Padova Millennium Galaxy and Group Catalogue (PM2GC; Calvi et al. 2011), Poggianti et al. (2013b found ∼ 4.23 × 10 −4 ℎ 3 Mpc −3 at 0.03 < < 0.11. We have plotted their number density in the top panel of Figure 17 (marked as ) for comparison. Both studies acknowledge that more compact quiescent galaxies were found in dense cluster environments (see also, Tortora et al. 2020). Indeed, Damjanov et al. (2015b) highlighted the environmental dependency of compact galaxies in the COSMOS field. At intermediate-, the compact quiescent galaxies align with the highdensity regions (see their Figure 8). Although, it is important to note that this might simply be a result of more massive quiescent galaxies (of all sizes) living in cluster environments compared to the field (Tortora et al. 2020).
The depth and field-of-view in our parent sample selection (see Section 2, and Equation 1) is wide enough to cover both clusters and fields. It captures the small Virgo Cluster overdensity at a distance of 17 Mpc and a slight void at 40-60 Mpc. We explore the influence of the Virgo Cluster has on our analysis, later in Section 6.4.3.

Variation arising from different MLCRs
One of the uncertainties on our, and everyone's, spheroid number density comes from the adopted MLCRs and thus the stellar mass measurement. For Figures 14 to 17, we presented the result using the RC15 MLCR. Here, we will explore the consequences of applying different MLCRs. Figure 18 shows the number density of the local compact spheroids (marked with a •) and the ellicular and elliptical galaxies (marked with a ) based on the four different MLCRs presented in Section 2.3 and the four different selection criteria presented in Section 4.3. We have also highlighted the maximum number density of red nuggets from various works with colour-coded dashed lines.
Using the size-mass selection criteria from Barro et al. (2013) and Damjanov et al. (2015a), the number densities of compact massive spheroids are higher than elliptical+ellicular (E+ES) galaxies, for all four MLCRs (see the upper green panel and the lower red panel in Figure 18). Applying the more restrictive size-mass selection criteria from van der Wel et al. (2014) and van Dokkum et al. (2015) reduces the number of compact massive spheroids included in the sample selection boundaries (see the upper-middle panel shaded blue and the lower-middle panel shaded yellow), lowering it to roughly less than or equal to the fixed number (density) of E+ES galaxies among our sample of 103 galaxies. When using the T11 stellar masses, no compact massive spheroid was found with the restrictive size-mass selection criteria. However, when using the higher IP13 stellar masses, the observed number densities of compact spheroids increased to greater than the number of E+ES galaxies, in all criteria.
The comparison showcases how the number density varies under the different stellar mass measurement and compactness selection criteria. Importantly, the MLCRs from T11 and RC15 assume the same IMF and SPS as used by those studying the high-red nuggets. The result from the less restrictive boundaries in the size-mass plane is consistent when using these two MLCRs. With the more restrictive boundaries in the size-mass plane, and since no compact system was found using the MLCR from T11, we can conclude that the number density of local 'compact massive' spheroids is subject to the particular * / ratio employed. Resolving this, beyond having reported its impact, is clearly outside the scope of the current investigation.
Our results are, however, robust when using the more inclusive size-mass selection criteria from Barro et al. (2013) and Damjanov et al. (2014). For all stellar mass estimates, we found a consistent trend of c,Sph > E+ES . Considering that c,Sph is just a lower limit for the local compact massive spheroids, while E+ES is close to the true value, our result calls for a reevaluation of the nature of the sizeevolution experienced by red nuggets. At a minimum, disc growth has to play an essential role in shaping the high-galaxies in the lower mass range (1 × 10 10 < * /M (RC15) < 4 × 10 10 ).

DISCUSSION
We have established that only ∼11% of the galaxies in our sample of massive galaxies are true elliptical galaxies; the majority are multicomponent systems with bulges, discs, bars, and other lesser components. Our trend between the bulge, aka spheroid size, and mass agrees well with data from comparable, multi-component decomposition analyses of nearby galaxies. That is, we do not detect any systematic bias in our measurements. We derive a lower limit to the number density per unit volume of (hidden) compact massive spheroids using four different size/mass selection criteria, giving c,Sph ∼ (0.17-1.20)×10 −4 Mpc −3 . We additionally examined how c,Sph changes when adopting the four different MLCRs. We found that if only high-mass spheroids ( * /M (RC15) > 4 × 10 10 ) are considered, c,Sph can be either lower, comparable or higher than E+ES . The uncertainty in stellar mass makes estimating the correct number density challenging. In the case of low-mass spheroids (1 × 10 10 < * /M (RC15) < 4 × 10 10 ), such an issue does not exist.
In what follows, we further compare local compact massive spheroids with distant red nuggets, and we address proposed sizeevolution hypotheses based on the results from our findings.

Alternative evolutionary processes
Here, we provide a brief overview of different conjectures regarding the size-evolution pathway taken by high-compact massive quiescent galaxies. Trujillo et al. (2007) plotted the half-light galaxy radii of local, luminous galaxies with Sérsic indices > 2.5. These galaxies are a factor of four larger than the compact spheroidal-like systems at ∼ 1.5. From this, they concluded that the ∼ 1.5 compact massive galaxies evolved through dry mergers with themselves, i.e., equalmass major mergers without star formation. van  reported a typical growth in galaxy half-light radii by a factor of 5-6 since ∼ 2.3. They found that only 10% of the massive spheroidlike galaxies at ∼ 2.3 had galaxy sizes equal to that of elliptical galaxies of comparable mass at ∼ 0. They, therefore, argued against 'monolithic collapse' as the formation channel because they found only 10% are fully assembled by ∼ 2.3. They noted an array of processes to possibly explain the size growth, including dry mergers, satellite accretion, expansion after mass loss from stellar winds, and the arrival of new, larger galaxies since ∼ 2.3 onto the red sequence.

(i) Compact massive quiescent galaxies experience rapid size evolution, an E-to-E scenario
The notion of major mergers driving size evolution has been heavily challenged. The observed number of major mergers (1:1 to 1:4 mass ratios) is too low to explain the change in number density of massive (log 10 ( * /M ) 11.1) galaxies (Man et al. 2016). Taylor et al. (2010) found no galaxies in the SDSS DR7 with comparable sizes to the van  and Damjanov et al. (2009) red nuggets at higher redshifts. They also argued against stochastic major mergers as the critical driving force for the size evolution in earlytype galaxies because it will result in some galaxies (those which experienced many mergers) having masses larger than observed in the local Universe. Carrasco et al. (2010), who studied eight massive galaxies at 1 < < 2 from the POWIR survey , reported extreme compactness within < 2 kpc (higher than the stellar densities of local galaxies). With similar reasoning as Taylor et al. (2010), they advocated that minor mergers and accretion are more likely the cause of the size growth.
While some kind of E-to-E scenario is widely accepted in the community (e.g., Buitrago et al. 2008;Bezanson et al. 2009;Hop-kins et al. 2009;Szomoru et al. 2011;Newman et al. 2012), the mechanism which supports a dramatic E-to-E size growth remains a subject of debate. Fan et al. (2008) had proposed that through mass loss (stellar winds, AGN jets), galaxies experience adiabatic expansion as the new gravitational equilibrium relaxes. However, adiabatic expansion has been called into question by numerical simulations, as Ragone-Figueroa & Granato (2011) found galaxies 'puffing-up' due to stellar winds occurs when the stellar population is younger than the formation of the red nuggets (by 0.5 Gyr). Hopkins et al. (2009) and Bezanson et al. (2009) advocated that, through minor mergers and accretion, red nuggets may build up a low-density three-dimensional envelope which hides the compact galaxies as the core of a larger elliptical galaxy (see also Hopkins et al. 2010;Szomoru et al. 2011;López-Sanjuan et al. 2012;Oogi & Habe 2013).
(ii) Compact massive quiescent galaxies experience only a mild, or no, evolution and now reside in local galaxy clusters: an E-to-E/S0 scenario.
(Galaxy cluster)-focused studies at intermediate-have challenged the universality of the significant size-growth scenario. Valentinuzzi et al. (2010) and Saracco et al. (2010) report on the abundance of 'compact' galaxies in cluster environments and postulate that red nuggets live in dense environments nowadays. Similarly, Poggianti et al. (2013b) found that the fraction of local compact massive galaxies are three times higher in clusters than in the field, arguing for the scenario of mild evolution (by a factor of ∼1.6 in size) for high-red nuggets. Damjanov et al. (2015b) echo the same point as they subsequently found compact massive quiescent galaxies at 0.1 < < 0.4 in the COSMOS field (Scoville et al. 2007) resided in denser environments than other quiescent galaxies. Carollo et al. (2013) argued that the addition of more large newly quenched early-type galaxies, as time progresses, is responsible for the apparent size evolution (a.k.a., the 'progenitor bias ', van Dokkum & Franx 1996;Saglia et al. 2010). In this picture, the size growth of quiescent galaxies since ∼ 2 is of secondary importance.
(iii) Compact massive quiescent galaxies experience size evolution via disc growth: a disc-cloaking (E-to-S0/S) scenario Graham (2013) and Graham et al. (2015) expanded the possible evolutionary channels, where instead of only considering elliptical (E) galaxies as the end product of red nugget evolution, lenticular (S0) and early-type spiral (Sa) galaxies may arise by cloaking the compact spheroid with a large-scale disc built through minor mergers and gas accretion. Upon entering a cluster of galaxies with a hot X-ray gas halo, any future disc growth via gas accretion and star formation would naturally be curtailed. The hypothesis of disccloaking was supported by a reported lower-limit to the number density of ∼ 6.9 × 10 −6 Mpc −3 , based on 21 local compact massive spheroids within 90 Mpc . It was a lower-limit because the volume was not sampled completely; the 21 systems simply represented spheroids that had been reported in the literature as having small sizes and high masses and were known to the authors. Their finding was further supported by de la Rosa et al. (2016) using the automatic Bulge+Disc decompositions of local galaxies imaged by SDSS (Mendel et al. 2014). From a sample of SDSS DR7 galaxies and a survey area of 8, 032 deg 2 at 0.025 < < 0.15, they found a number density of local compact bulges ( ∼ 0.28-3.12 × 10 −4 Mpc −3 ) comparable with the number density of the distant red nuggets across multiple size-mass sample-selection criteria (Barro et al. 2013;van der Wel et al. 2014;van Dokkum et al. 2015). Because galaxies can have (truncated) non-exponential discs and more than two components, such as bars and rings, etc., automated Bulge+Disc decompositions often fail to provide a reliable measurement of the bulge size and luminosity, which is why we undertook our careful investigation.

Differing selection criteria and assumptions
The criteria of selecting for the parent sample, measurement of the galaxies' properties, and the definition of compactness, vary from one study to another. It will be a monumental task to homogenise all the data presented in the literature. Instead, we provide a brief description of various approaches and assumptions commonly made (see Table 5), and we highlight a few potential causes of discrepancy.

Parent sample
The majority of studies select their parent galaxy set based on the galaxies' stellar masses, redshifts, and star formation (or its proxies, e.g., colour-colour selection). Sometimes the sample is further restricted by morphology or brightness concentration (Sérsic index). Column (2) in Table 5 summaries such unique conditions for several studies. Trujillo et al. (2007) divided their sample of luminous galaxies into 'spheroidal' (Sérsic index > 2.5) and 'disc-like' (Sérsic index < 2.5) to illustrate the stronger size evolution in spheroids, and the lack thereof in disc galaxies. The Bezanson et al. (2009) local sample is morphology-dependent, using the apparent 'E' galaxies from Tal et al. (2009). Some studies did not construct the parent sample based on star formation (Trujillo et al. 2007(Trujillo et al. , 2009Taylor et al. 2010), where both quiescent and star-forming galaxies are considered. The sample from Carollo et al. (2013) is based on the ZEST+ morphology classification (an upgraded version from Scarlata et al. 2007), from which they select both massive early-type ('E/S0') and late-type ('Sa-cd') galaxies. Studies that have a parent sample inclusive of all morphologies would not be biased by an E-to-E scenario presumption. Although, depending on the lower-mass threshold in their 'compactness' criteria, individual works might be probing a different population of red nuggets. Indeed, for instance, applying a higher mass selection from Trujillo et al. (2009)

Size measurements
The method of size measurement varies slightly between studies. The most commonly used method is fitting a single Sérsic 1/n model to the surface brightness profile of the galaxies (e.g., Daddi et al. 2005;Trujillo et al. 2007Trujillo et al. , 2009Bezanson et al. 2009;Valentinuzzi et al. 2010;Cassata et al. 2011;Mancini et al. 2010;Saracco et al. 2010;Poggianti et al. 2013a;Barro et al. 2013;van der Wel et al. 2014;van Dokkum et al. 2015). From this, the effective radius e (essentially the half-light radius 1/2 ) is generally used to represent the size of a galaxy in discussions. Some studies, however, fit the de Vaucouleurs (1948) 1/4 model instead. For example, Taylor et al. (2010) rely on the parametric measurements from Strauss et al. (2002), who used an 1/4 model to describe the galaxy surface brightness profiles (see also section 3.2 in Taylor et al. 2010, for a relevant discussion). For red nuggets, Damjanov et al. (2009) fit the 1/4 model to their MUNICS (Drory et al. 2001) data and both 1/4 and 1/ models to their GDDS (Abraham et al. 2004) data. If one assumes an 1/4 model, the effective radius could be overestimated if the galaxy has a Sérsic index < 4. In Column (3) of Table 5, we list the measurement method.
The colour of an image also affects its size measurement Vulcani et al. 2014;Kennedy et al. 2016). In Column (4) of Table 5, we show the band in which the study was conducted. Importantly, the size of early-type galaxies is smaller in the red band as compared to the blue band. As seen in La Barbera et al. (2010), -band Salpeter (1955) 3  Baldry & Glazebrook (2003) and Salpeter (1955) IMFs, respectively. there is a size decrease by 30% from -through -band, and a 25% decrease from the 160 to the 625 filter (Marian et al. 2018).

Stellar mass
The galaxy stellar mass can be estimated via broadband SED fitting (e.g., Daddi et al. 2005;Longhetti et al. 2007;Kriek et al. 2008;Barro et al. 2013) 22 . The underlying assumptions for such an operation differ between studies and, as we have shown, the stellar mass calculation affects the reported number density. One such factor is the assumed stellar population model. Most works in Table 5 use the Bruzual & Charlot (2003)  Another factor is the choice of the initial mass function (IMF). Here, we only show, in Column (5) of Table 5, the IMF used in the respective studies for reference purposes. The Chabrier (2003) IMF is the most widely-used and also the most bottom-light compared to the other contemporaries. The CANDELS-based (Grogin et al. 2011;Koekemoer et al. 2011) size evolution studies (Barro et al. 2013;van der Wel et al. 2014;van Dokkum et al. 2015) all assume the Chabrier (2003) IMF. However, some low-to-intermediate redshift studies Valentinuzzi et al. 2010;Poggianti et al. 2013b;Carollo et al. 2013) instead assume the more bottomheavy Salpeter (1955) IMF, or the pseudo-Kroupa IMF (Kroupa 2001) rather than the Kroupa (2002) IMF. Using a Chabrier (2003) IMF will yield 0.24 dex less mass than the Salpeter (1955) IMF (Salimbeni et al. 2009). A bottom-heavy IMF could result in higher stellar mass, and therefore, more spheroids will be included by the 'compact massive' selection criteria.

Compactness criteria
As discussed in Section 4.3, the compactness criteria are usually somewhat arbitrary or it borders the distribution of the general galaxy population's size-mass relation. The lower-mass limit, in particular, can alter the number density of 'compact' systems in one's sample. We have listed the lower-mass limit for each study in Column (6) of Table 5. If the study did not apply a clear size-mass selection boundary, the lowest mass limit (corresponding to the magnitude limit) of the parent sample is provided. These lower-mass limits vary from * /M = 10 10 to 10 11 . If a study assumes a bottom-heavy IMF and a lower-mass limit, one can expect a higher reported number density.

The possibility of E-to-E evolution
Our work is fundamentally different from the previously mentioned works. By not treating galaxies as single-component systems, our analysis provides additional information on the galaxies' structural features. Based on the RC3 morphology, the frequency of alleged 'E' galaxies in our local sample is abundant. However, after performing analysis and decomposition of the galaxy light, the number of elliptical galaxies is revised downward, already weakening the argument that the high-red nuggets transformed into today's elliptical galaxies, simply because the actual number of elliptical galaxies is sparse. Only 13/28 originally labelled 'E' galaxies are true ellipticals (11 E + 2 ES, see Section 3.4). This reflects revelations from a quarter of a century ago, whereby the true number of pressure-supported elliptical galaxies was found to be notably lower than previously thought (e.g., D'Onofrio et al. 1995;Graham et al. 1998).

The morphology of host galaxies with an embedded compact massive spheroid
The literature abounds with examples of what were formerly elliptical galaxies that are reclassified as disc galaxies. Early examples are NGC 4111 (E7 → Sa: Adams & Seahes 1937, their p.31) and NGC 3115 (E7 → S0: Burbidge et al. 1961).
Our results cast doubt on the evolutionary scenario of red nuggets transforming into elliptical galaxies. In Figure 19, we plot the geometric-mean, i.e., circularised, half-light radius and total mass of the host galaxies, shown by the green plus signs ('+'), and their respective embedded spheroids (shown by the solid red circles). The three prominent outliers are marked with red hollow circles 23 . The sample shown here are the 55 spheroids that satisfy the compactmassive selection criteria of Barro et al. (2013), the most inclusive criteria. We divided the sample into two categories, they are: 23 For the sake of mimicking how most studies measure galaxy size, we use a single Sérsic function to measure the host galaxy's effective radius, although we know it is not the best measurement for some galaxies (see also Appendix C). Case in point, the three outliers are NGC 3158 (S0), IC 983 (S), and NGC 3646 (S). All have a disc component with a rather shallow profile and large scale length ℎ > 20 kpc. Combined with a prominent bulge, this creates a curvature in the surface brightness profile that requires a single Sérsic function with an exceptionally high index ( > 5). This highlights a problematic issue with the conventional method of measuring galaxy size and shall be pursued further in future work.
(i) spheroids embedded in lenticular galaxies 24 ('S0 hosts', see the upper panel in Figure 19); and (ii) spheroids embedded in spiral galaxies ('S hosts', see the lower panel in Figure 19).
None of the compact massive spheroids, as identified using the selection criteria of Barro et al. (2013), are pure elliptical galaxies. A significant number of compact massive spheroids with stellar masses above (4-5)×10 10 M come from S0 host galaxies, while those with stellar masses lower than (4-5)×10 10 M are slightly more likely to come from a spiral host galaxy. Additionally, none of the 11 E + 2 ES galaxies (when using the T11, Z09, or RC15 mass prescriptions) satisfy the 'compact-massive' selection criteria of any paper 25 . What this means is that all of our compact massive spheroids are embedded inside either S0 or early-type S galaxies; none are compact massive elliptical galaxies, although a few such galaxies do exist locally. Hopkins et al. (2009) proposed an evolutionary pathway for red nuggets, where, through dry minor mergers, the high-compact galaxies develop a pressure-supported 3-D envelope to become the cores of today's elliptical galaxies. This is a different evolutionary path to acquiring a rotating 2-D disc to become the bulges of today's lenticular galaxies. The answer to differentiating between these two evolutionary paths depends on understanding the surface density profiles of massive local galaxies as a function of radius. One can also resort to kinematic information to ascertain the presence of a rotating disc. Hopkins et al. (2009) investigated 180 early-type galaxies from Lauer et al. (2007) and Kormendy et al. (2009), and found similarities between the inner part, ∼(1-5) kpc, of the local early-type galaxy light profiles and the distant compact massive quiescent galaxies in the sample of van . Earlier in the same year, Bezanson et al. (2009) had reached a similar conclusion, stating that the inner ( < 1 kpc) stellar densities of local elliptical galaxies are similar to the high-red nuggets. Possibly due to the difficulty distinguishing E from S0 among ETGs, papers advocating for an E-to-E evolution kept emerging (e.g. Oldham et al. 2017;Tortora et al. 2020;Zhu et al. 2021), despite the fact that the similarity of the profile's inner density did not actually prove the proposed E-to-E evolution over an E-to-S0 evolution. Arnold et al. (2011) suggested minor mergers were building a pressure-supported halo around NGC 3115. De et al. (2014) reported structural similarity in the inner part at radius < 1 kpc of local early-type galaxies with the red nuggets at 2 < < 7, and virtually no correlation at > 10 kpc. However, information on the ∼ 0 number density of S0 galaxies with compact massive cores is needed. This required identifying which early-type galaxies are S0 or E, or ES.

3-D envelopes
While the studies mentioned above provide a legitimate argument for the inside-out growth scenario, our result provides additional insight into the nature of this evolution. Those studies generally did not attempt to separate the early-type galaxies into a disc (lenticular) and no-disc (true elliptical) systems. By and large, they assumed 24 The morphologies here refers to the reclassification we introduced in Section 3.4. 25 Only when applying the highest stellar mass estimation (IP13), did we find 3, 3, and 4 'compact massive' galaxies when using the Barro et al. (2013), Damjanov et al. (2014), andvan Dokkum et al. (2015) selection criteria, respectively. However, because IP13 assumes a different IMF compared to these studies, it is not exactly a fair comparison. evolution into larger pressure-supported systems rather than entertaining the notion of substantial disc growth to create a different type of galaxy with significant ordered motion. Barbosa et al. (2021b) conducted a detailed kinematic and population analysis on a bona fide elliptical, NGC 3311, and found a compact core, a high-'relic' galaxy, is hidden inside < 2 kpc. Similar studies would be much welcomed to confirm an E-to-E transition. Our bulge-to-disc and disc-to-total ratios reveal that the discs, and their bars and ansae, are nowadays the dominant mass component in many local massive galaxies. This has considerable implications on many fronts, such as the inferred dark matter halos of massive early-type galaxies implied from their central velocity dispersion, , via the expression ∼ 2 e,gal (e.g., Poincare & Vergne 1911;Zwicky 1937;Poveda 1958Poveda , 1961. Obviously, the larger e,gal radii, which are dominated by the size of the disc, is not applicable for use in this expression for a virialised, pressure-supported galaxy. Hopkins et al. (2009) used data from Kormendy et al. (2009), a selected set of 27 elliptical (E) galaxies from the Virgo Cluster and with very few massive galaxies ( * > 3 × 10 11 M ). They modelled all the E galaxies with a single Sérsic fit. However, upon inspection, some of these galaxies exhibit features which may be a signature of a rotational component.  Hopkins et al. (2009) also used data from Lauer et al. (2007) based on 219 early-type galaxies taken from different sources (Lauer et al. 1995;Faber et al. 1997;Quillen et al. 2000;Ravindranath et al. 2001;Rest et al. 2001;Laine et al. 2002;Lauer et al. 2005). In terms of morphologies, roughly half of the sample was thought to be 'E', a quarter of them 'S0', and the other quarter 'BCG' according to the RC3 classifications. Most BCGs are ellipticals, yet without a multi-component analysis of kinematic data, the nature of the alleged 'E' galaxies is not confirmed.
The local galaxy sample used by Bezanson et al. (2009) came from Tal et al. (2009), with well-defined distances (15 Mpc < Dist. < 50 Mpc) and absolute magnitude ( < −20 mag). They pointed out half of their sample exhibits morphological disturbance features, such as shells and tidal tails, implying heavy interactions with their environments. De et al. (2014) supports the 3-D envelope scenario based on the three-Sérsic-component decompositions by Huang et al. (2013) on the Carnegie-Irvine Galaxy Survey data (Ho et al. 2011). The 94 galaxies they investigated used to model listed as 'E' in the RC3 classifications. Huang et al. (2013) used three Sérsic functions to describe the inner, middle, and outer region of the galaxy. They found a high correlation between the inner region of the local galaxies and the high-red nuggets but not the outer region, hence concluding the red nuggets went through a multi-stage buildup in an insideout manner. This, however, does not rule out disc growth. Indeed, as Huang et al. (2013) pointed out in their Table 1, the following galaxies have an edge-on disc: NGC 3585,and NGC 7029,and these: IC 4797,NGC 584,NGC 3904,NGC 4033,NGC 4697,NGC 6673,NGC 7145,NGC 7192,NGC 7507 are possibly S0 galaxies. Instead of a multi-layered elliptical galaxy, there could be many S0 galaxies in the mix, and the inner region of these galaxies is, in fact, the bulge.
As we have established in Section 3.4, without exploring the galactic substructure, the morphology of early-type galaxies is often mistaken. Instead of comparing the high-quiescent galaxies to true elliptical galaxies, the comparisons are essentially being made using many local lenticular galaxies. It is likely that many of the alleged low-density 3D envelopes are actually 2D rotating discs seen somewhat face-on. Rather than evidence of 3D envelope growth, many of the galaxies used in the above studies instead support the 2-D disc cloaking process.

Evidence of disc accretion
Lyman-gas clouds are observed within 500 kpc of high-galaxies of all types (Wakker & Savage 2009;Prochaska et al. 2011;Thom et al. 2011;Stocke et al. 2013;Tumlinson et al. 2013). For instance, Bouché et al. (2013) used a background quasar to trace the nearby gas (damped Lyman absorber 26 kpc away) of a star-forming galaxy at = 2.3 with a typical rotational disc, according to its kinematics (Förster Schreiber et al. 2009). The kinematics of the gas shows a signature of cold accretion, where it appears to be low-metallicity, coplanar, and co-rotational compared to its disc. In the local Universe, Coccato et al. (2013) examined the two counter-rotating discs in both NGC 3593 (S0/a) and NGC 4550 (S0), where the secondary discs differ from the primary discs in metallicity and stellar populations. They favour a scenario of gas accretion forming the second disc ∼2-7 Gyr ago.
High-red nuggets are reported to have reached their peak abundance at ∼ 1.2-1.6, and dropped in number ever since (Barro et al. 2013;van der Wel et al. 2014;van Dokkum et al. 2015). The disappearance of these galaxies through the growth of discs may involve the creation of star-forming galaxies (should much of the disc stars be built from accreted gas) or not (when the bulk of the disc is built from accreted galaxies with little gas).
There is an abundance of evidence from different approaches that confirm the prominence of discs in local early-type galaxies. Kaviraj et al. (2007) investigated the UV and optical photometry of ∼2100 SDSS early-type galaxies at 0 < < 0.11 and found at least ∼30% exhibit recent (∼1 Gyr ago) star formation activity. Fabricius et al. (2014) performed a detailed analysis of NGC 7217, a spheroiddominated galaxy (classified as an early-type spiral in Buta et al. 1995) with two distinct rotational components. The active star formation of the rotational structures suggested an ongoing (re)growth of the stellar disc. While this process is truncated in cluster environments, it can still occur in the field (e.g. Graham et al. 2017).
While 'cosmic noon' has passed, and there are less substantial cold gas clouds nearby for galaxies to accrete, it is an ongoing, downsizing phenomenon. The kinematics of the gas shows a signature of cold accretion, where it appears to be low-metallicity, coplanar, and corotational compared to its disc (see also, Stewart et al. 2013;Kacprzak et al. 2010;Martin et al. 2019;Zabl et al. 2019;Kacprzak et al. 2015;Nielsen et al. 2017). Davis et al. (2011) examined the gaseous content of local early-type galaxies with CO and H interferometric observations. The kinematic misalignment between the gas and stellar components implied that 42% of the gas in field galaxies (fast rotators) is ex-situ in origin, namely from cold accretion and minor mergers 26 . Alatalo et al. (2013) investigated the morphology of molecular gas through CO imaging by the CARMA ATLAS 3D survey, finding half of the 40 CO-rich early-type galaxies exhibit a CO disc. They estimated (with a correction factor from Kaviraj et al. 2012) the lower-mass limit to the total accreted gas to be 2.48 × 10 10 M across 15 galaxies, a value consistent with the mass brought in by minor mergers (Lotz et al. 2008) and thus the accretion of not just gas but also stars from smaller galaxies.
As mentioned before, the extensive kinematic observations from the ATLAS 3D survey revealed that the majority of the local earlytype galaxies are 'fast rotators' (Emsellem et al. 2011). Krajnović et al. (2011) shows that most early-type galaxies (∼82%) are 'regular rotators' (defined in Kinemetry, Krajnović et al. 2006Krajnović et al. , 2008, in the sense that the velocity maps are dominated by ordered rotation. Krajnović et al. (2013b) performed Bulge+Disc decompositions on a set of 180 of the non-barred galaxies. They found that galaxies with high angular momentum ( e ) also have a large disc-to-total flux ratio. The combination of kinematics and photometric results revealed that 83% of their early-type galaxies have a rotational disc that accounts for around ∼40% the stellar mass of a galaxy.
Our result highlighted the abundance of discs in massive local galaxies, providing strong support for the disc growth (cloaking) scenario. The cyclical nature of galaxy metamorphosis has been recognised widely in simulations (White & Rees 1978;White & Frenk 1991;Navarro & Benz 1991;Steinmetz & Navarro 2002;Bournaud & Combes 2002). If disc cloaking were to occur to our spheroids, and, subsequently, build up an extended disc and other substructures in an inside-out manner, the systems will experience growth in both size and mass. They would migrate outward (see the blue arrows, in Figure 19) from the compact massive region (the green shaded area) in the size-mass diagram and become the massive ( * /M > 10 11 ) local S0 and S galaxies. Depending on the nature of the minor mergers (wet or dry) and the amount of cold stream gas accretion, the galaxies may also migrate outside of the 'quiescent' region in the colour-colour diagram to enter the 'star-forming' blue region because of the newly built rotational discs. With time, star formation in these discs turns off, and a quenched, lenticular galaxy now moves back into the 'quiescent' region with a larger overall size.

The role of IMF in stellar mass estimation
It is questionable whether there is a universal IMF among galaxies Bastian et al. 2010;van Dokkum & Conroy 2010;Cappellari et al. 2012;Ferreras et al. 2013;La Barbera et al. 2013;Smith 2014a;Spiniello et al. 2014). The slope of the low-mass end of the IMF is thought to be directly correlated with the velocity dispersion of galaxies Domínguez Sánchez et al. 2019) and, consequently, the colour of galaxies (Dutton et al. 2012;Pforr et al. 2012;Ricciardelli et al. 2012;Vazdekis et al. 2012). For example, Martín-Navarro et al. (2015) analysed the 'relic' galaxy NGC 1277 ) that is compact (1-2 kpc) and has an exceptionally high velocity dispersion ( > 400 km s −1 , see also Ferré-Mateu et al. (2017) for its detailed features). However, the IMF correlation with velocity dispersion has been called into doubt. Importantly, from one BCG, NGC 3311, Barbosa et al. (2021a) presented an important case showing the IMF-correlation to be invalid, a mere coincidence that old stars are found in the area of a galaxy with high dispersion. They found the tightest relations to be between stellar age-to-IMF and radius-to-IMF. They found the IMF to be bottom-heavy and requiring a high Υ * . Ferré-Mateu et al. (2013) has shown the non-universality, whereby adjusting the slope of the IMFs according to the central velocity dispersion yields a more comparable star formation history among the early-type galaxies. While Smith (2014b) pointed out the discrepancy between methods for probing the IMF of galaxies, based on the data from ATLAS 3D and Conroy & van Dokkum (2012b), there was a broad agreement that elliptical galaxies have a more bottom-heavy IMF.
If a bottom-heavy IMF is indeed better suited for high dispersion galaxies, then among the MLCRs in Section 2.3, IP13 perhaps portrays a more accurate stellar mass since it assumes the Kroupa (1998) IMF -yielding log( * / ) ratios 0.225 dex lower than those obtained with the Salpeter (1955) IMF, see Flynn et al. (2006, their Figure 12) -instead of the more bottom-light Kroupa (2002) IMF or the Chabrier (2003 IMF (see also Ferreras et al. 2013;Kroupa et al. 2013). Although, there is plenty of evidence Cenarro et al. 2003;Falcón-Barroso et al. 2003;van Dokkum & Conroy 2010;Conroy & van Dokkum 2012a;Spiniello et al. 2012;Conroy & van Dokkum 2012b;van Dokkum & Conroy 2012) indicating there are more low-mass stars in elliptical galaxies than even the most bottom-heavy IMF (Salpeter 1955) predicts. This implies that even our highest mass estimation (the IP13 masses) is an underestimation of the real stellar mass. Throughout this work, however, the relative stellar mass compared to the high-measurements is more important to us. Because the studies of high-red nuggets assume a bottom-light Chabrier (2003) IMF, for comparison purposes, we need to apply the same assumption as we estimate the spheroids' stellar mass (T11 and RC15 MLCRs), even if the real stellar mass is larger than what was presented in Figure 15.

The similarity between local spheroids and distant red nuggets
In order to fully prove that high-red nuggets and local spheroids are drawn from the same population, we need matching sizes, stellar masses (and thus stellar densities), internal dynamic features and stellar populations. Such a full comparison for our sample is beyond the scope of this paper but is perhaps not required given the extensive evidence in the literature. The bulges of lenticular galaxies and high-quiescent galaxies both have high velocity dispersions. For example, Martinez-Manso et al. (2011) reported on a range of dispersions, ∼ 156-236 km s −1 , for four red nuggets at ∼ 1 taken from Trujillo et al. (2007). van de Sande et al. (2013) measured five galaxies at ∼ 2 to have ∼ 290-450 km s −1 . These highest values appear as outliers after Belli et al. (2014) analysed 56 quiescent galaxies at 0.9 < < 1.6 and obtained each galaxy's mean velocity dispersion within e,gal . They found an average e value of 219 km s −1 . For comparison, in the local Universe, the MASSIVE survey (Ma et al. 2014) reports on the central velocity dispersions c for 90 early-type galaxies within 104 Mpc. They are in the range of c ∼ 200-350 km s −1 (Veale et al. 2017). In the redshift range 0.15 0.5, Scognamiglio et al. (2020) detected 19 UCMGs ( * /M > 8 × 10 10 and e < 1.5kpc) spanning a range of 200 /km s −1 400 (see also, Tortora et al. 2016Tortora et al. , 2018Tortora et al. , 2020Barbosa et al. 2021b).
Local spheroids are also known to contain old stellar populations that existed at high-. Saracco et al. (2009) showed 32 early-type galaxies at 1 < < 2 which exhibit two distinct stellar populations: the old (∼3.5 Gyr) and young (∼1 Gyr) stars, implying the early-type galaxies exist at least since > 3-3.5. Peletier et al. (2007) also reached the same conclusion from the perspective of stellar velocity dispersion dips. MacArthur et al. (2009, see their table 4) found that in local spiral galaxy bulges, old stellar populations (>10 Gyr) make up the majority of the mass budget. Regardless of whether the bulges of local galaxies formed via outside-in (i.e., secular evolution) or inside-out (i.e., mergers and/or accretion) scenarios, the structural similarity to high-red nuggets is now well established. Although, evidence from stellar population studies (Proctor & Sansom 2002;Moorthy & Holtzman 2006;Thomas & Davies 2006;Jablonka et al. 2007;MacArthur et al. 2008;Saracco et al. 2009) favour the insideout process in early-type galaxies.
The stellar densities are also comparable between local bulges and high-red nuggets (Graham 2013, his figure 1), which have halflight densities higher than local early-type galaxies of the same mass. This is because the discs of local lenticular galaxies -which tend to dominate these galaxies' stellar mass budget -increase the galaxies' half-light radii but not in an economical space-filling fashion. The Bulge+Disc decompositions from de la Rosa et al. (2016) also show the size-mass relation overlaps (see their figure 3); the red nuggets at 0.18 < < 1 share almost the same relation and galaxies at 1 < < 1.8 have a similar slope, but are slightly more compact. Graham (2013) further showed that local dwarf 'compact elliptical' (cE) galaxies have similar sizes, stellar masses, and densities as local low-mass bulges, suggesting that they are the (largely) discfree counterparts of lower-mass bulges, just as red nuggets are the (largely) disc-free counterparts of high-mass bulges. Graham et al. (2015) also illustrated that the overlap in the (Sérsic index, )-(effective radius, e ) diagram (their figure 3, primarily from Bulge+Disc decompositions of local galaxies) between bulges and high-red nuggets from Damjanov et al. (2011). Given the structural similarity, they, therefore, concluded that not all compact massive spheroids have gone through significant size evolution since < 2.5.

Contributions from the Virgo Cluster
Bin 3 is of particular interest. Bin 3 captured galaxies within Dist. < 45 Mpc, a space containing the Virgo Cluster. Not only is it complete down to the lowest galaxy stellar mass ( * (IP13) = 10 11 M ), it also yields the highest number density for compact massive spheroids.
In Figure 20, we show the location of our host galaxies in the sky. Small black dots are all the galaxies in the SDSS field bounded by our angular selection (Equation 1). One can clearly see the largescale filament structures in the local Universe manifesting. The blue points are the 103 host galaxies in our sample. Since we only select for massive galaxies, it is natural that most of them live in dense environments. One can see the blue points largely coincide with the over-density regions. The red points are the Virgo Cluster host galaxies. There are 22 out of 52 galaxies in Bin 3 that are in the Extended Virgo Cluster Catalog (EVCC, Kim et al. 2014a).
Using the less restrictive Barro et al. (2013) and Damjanov et al. (2014) sample selection criteria, Bin 3 is more abundant with compact massive spheroids than the peak of the red nugget population ( c,Sph,Bin3 > max(n RN )) across all stellar masses. Note that the c,Sph,Bin3 here corresponds to the respective 'compact massive' selection criteria as the high-nuggets' peak abundance (max(n RN )). When using the more restrictive selection criteria from van der Wel et al. (2014) andvan Dokkum et al. (2015), our result varies depending on the stellar mass estimates (Figure 18). With T11 and Z09 stellar mass, no compact massive spheroid was found in Bin 3 with the restrictive van der Wel et al. (2014) andvan Dokkum et al. (2015) criteria. With RC15 stellar mass, Bin 3 has a lower number of compact massive spheroids to the red nuggets, c,Sph,Bin3 < max(n RN ) (around 0.66-0.78 difference). Finally, in IP13 stellar mass, c,Sph,Bin3 is slightly higher than max(n RN ), by 38% in van der Wel et al. (2014) criteria and52% in van Dokkum et al. (2015) criteria.
Having roughly half of the galaxies in Bin 3 being Virgo Cluster members, we examine if the cluster affects the estimation of the number density of the compact massive spheroids. We do so by removing the Virgo Cluster galaxies from our sample and recalculating the number density of compact massive spheroids in Bin 3. Table 6 shows the resulting number densities using the four criteria for defining compact massive spheroids (Barro et al. 2013;van der Wel et al. 2014;Damjanov et al. 2014;van Dokkum et al. 2015) and using the four stellar mass estimates (T11, Z09, IP13, and RC15).
With the restrictive selection criteria (van der Wel et al. 2014;van Dokkum et al. 2015), we found a similar number of compact massive spheroids as before. With the less restrictive selection criteria (Barro et al. 2013;Damjanov et al. 2014), the numbers of compact massive spheroids are roughly halved. Because the reduction (∼58%-61%) is comparable to the number of removed Virgo members (∼42%), we, therefore, conclude the Virgo Cluster over-density does not increase Bin 3's number density in any significant way. Our wide-sky coverage prevented a potential bias as the Virgo Cluster takes up less than one-eighth of the area in the R.A. and Dec selection window (see Figure 20). Moreover, with the reduced number densities ( c,Sph,Bin3 ≈ (0.3-5.0) × 10 −4 Mpc −3 ), Bin 3 still contains a comparable number of compact massive spheroids to the peak red nugget density (max(n RN ) ≈ (1.8-2.4) × 10 −4 Mpc −3 ). Thus, even without the influence of the Virgo Cluster, there remains a sizable amount of compact massive spheroids in Bin 3, sufficient to account for the missing red nuggets.

Some distant red nuggets may also have discs
It is unclear whether the original red nuggets (from Daddi et al. 2005;van Dokkum et al. 2008;Damjanov et al. 2009) are devoid of any disc growth. From visual inspection, some galaxies exhibit discy, high-ellipticity isophotes. For example, galaxy numbers: '4950' ( = 1.55) from Daddi et al. (2005) The existence of (nascent) discs in some high-galaxies is supported by van der Wel et al. (2011). They found a significant portion of the compact massive quiescent galaxies ( * /M > 10 10.8 , and e < 2 kpc) are flattened in projection, thus concluding that 40%-65% of red nuggets at ∼ 2 are disc-dominated. Ferré-Mateu et al. (2012) also found elongated morphologies in red nuggets at > 1 and classified them as fast rotators. One may therefore question whether those red nuggets are genuinely single-component systems. While spheroids can rotate, it is certainly the case that the low-resolution of images at high-will blend the bulges and discs together to create the appearance of a single-component system. Indeed, even at low , astronomers have struggled for a century to separate S0 galaxies from E galaxies. As such, the disc inclination angle and dust content will further complicate the situation. The shrouding dust heated by star formation will redden the galaxies' light (Pérez-González et al. 2008). Interestingly, a few notable local 'relic' galaxies (NGC 1277, PCG 032873, and Mrk 1216), that are believed to be untouched  Barro et al. (2013) 3.1 × 10 −4 (10) 3.1 × 10 −4 (12) 4.3 × 10 −4 (14) 6.4 × 10 −4 (21) van der Wel et al. by ex-situ processes, have strong stellar rotation. NGC 1277 has a rotational velocity ∼ 300 km s −1 in the outer region and a central velocity dispersion 0 > 300 km s −1 (Trujillo et al. 2014). In (Ferré-Mateu et al. 2017, see their Figure 2), both PCG 032873 and Mrk 1216 exhibit a rotational velocity of r ∼ 200 km s −1 and central dispersion of 0 ∼ 320 km s −1 , implying, indeed, they are not single-component systems. These systems may have accreted and grown their discs. Until they transition through the ES to S0 phase, they may remain compact massive galaxies. This may even conceivably explain some of the youthfulness seen in some ultra-compact massive galaxies (Spiniello et al. 2021), although we note this is just speculation on our part.
It will be interesting to perform structural decompositions on intermediate-and high-galaxies when high-quality images become available from facilities such as the European Southern Observatory's Very Large Telescope interferometer, equipped with PIONIER (1.6 µm, Le Bouquin et al. 2011), GRAVITY (2.0-2.4 µm imager; 4 mas seeing using the four 8-m Unit Telescopes, Gillessen et al. 2010;Eisenhauer et al. 2011), andMAVIS (optical, under development, McDermid et al. 2020;Rigaut et al. 2020). Subsequently, cosmological simulations involving disc cloaking can be tested if the observed high-mass end of the spheroid mass function, and the spheroid size-mass relation can be reproduced, reported here for the local Universe. Careful analysis will enable one to measure the growth of discs in early-type galaxies, shifting the focus away from spheroids -which appear to be relatively inert.

Summary and Conclusions
We defined volume-and mass-limited samples of local galaxies from the SDSS that are nondiscriminatory against galaxy morphology (Section 2). The galaxies have masses greater than 10 11 M and reside within 110 Mpc. We have discussed how the different massto-light ratios (Section 2.3) and distance measurements (Section 2.4) would affect the estimation of the stellar masses and sizes of the galaxies and their spheroids. A detailed, case-by-case study of 103 galaxies was undertaken. We performed multi-component decompositions (Section 3.3) on the surface brightness profiles of each galaxy. Rather than using automated routines which blindly fit several Sérsic components, we carefully examined each image and identified the morphological features that are present. This was further supported by recourse to the literature, in particular, studies including kinematic data or high-resolution HST data. This enabled a considerable improvement upon the conventional Bulge+Disc decomposition. Our analysis encompasses the less prominent yet important substructures (i.e., bars, ansae, spiral arms, nuclear discs, intermediate-scale discs, etc.), which can skew the bulge parameters if not taken into account. From our investigation, we have made the following observations.
• New morphological classifications were required, and made, for the galaxies. Previously overlooked substructures, such as large-and intermediate-scale discs, are now incorporated into the description of the morphology (Section 3.4). This is important because the morphology reflects the physical processes which shaped the galaxies' evolution.
• Implications from the low number of 'true elliptical' galaxies in the local Universe rule out the dominance of the popular E-to-E evolutionary path whereby quasi-spherical 3-D envelopes accumulate around red nuggets for relatively low-mass red nuggets ( * /M < 4 × 10 10 ). The possibility of such an E-to-E scenario is discussed in Section 6.2.
• Morphology-dependent / flux ratios are presented for the massive galaxies in Figure 12 (Section 4.1).
• The size-mass distribution of the local spheroids is presented (Section 4.3). The number density of 'compact' and 'massive' spheroids is obtained based on five different selection criteria.
• We compared the number density of local compact massive spheroids with intermediate-and high-redshift galaxies. The evolutionary trend for the number density of compact massive quiescent systems is presented in Figure 17, based on four different sample selection criteria (Sections 5.2 and 5.3).
Two important insights can be drawn from our results: (i) The local elliptical galaxies are not as abundant as previously believed. Among the 28 supposed elliptical galaxies (E), 20 are actually lenticular galaxies (S0). There are only ∼11% (11/103) true elliptical galaxies, plus 2 ES galaxies, in our sample.
(ii) We found compact massive spheroids hidden in the local lenticular and early-type spiral galaxies. Working with our host galaxies' selection limit ( * ,gal /M (RC15) > 6.7 × 10 10 ), we present a range of lower limits of the number density to the compact massive spheroid of c,Sph ∼ (0.17-1.20) × 10 −4 Mpc −3 , as defined by different definitions of 'compact massive' systems. If relatively low-mass ( * /M ∼ (1-4) × 10 10 ) spheroids are included, then there are sufficient numbers of compact massive spheroids ( c,Sph ∼ 1.20 × 10 −4 Mpc −3 ) to match the peak in the number density of red nuggets when defined using the same mass and size criteria. If we select for only high-mass spheroids ( * /M > 4×10 10 ), the number densities for compact massive spheroids and true ellipticals are comparable ( c,Sph E+ES ). However, in the high-mass range, the uncertainty in stellar mass estimation makes it difficult to conclude which evolution dominates by number density alone (Section 5.4). Interestingly, for * ,gal /M (RC15) > 6.7 × 10 10 and Dist. < 45 Mpc (Bin 3), we found more compact massive spheroids than the red nuggets' peak abundance: c,Sph,Bin3 > max(n RN ), inspite of the influence from Virgo Cluster galaxies.
Collectively, this calls the frequency of the E-to-E scenario into doubt, while the disc-cloaking scenario is now more salient than ever. An E-to-E evolution and disc-cloaking create distinctly different end products. If we are to truly understand galaxies, and the cosmological evolution of galaxies, this distinction is of key importance. Figure 18. The number densities of local compact spheroids (•) and ellip-ticals+elliculars ( ), using four different MLCRs, based on four size-mass selection criteria. The selection criteria applied from top to bottom are given by Equations 11, 12, 14, and 13, respectively. From left to right, the MLCRs are from T11, Z09, RC15, and IP13, respectively. The red, blue, and black colours represent the samples from Bins 1, 2, and 3, respectively. The black and red horizontal lines are c,Sph and E+ES , as depicted in Figure 17. Here, the number density of 'E+ES' includes 11 true E galaxies and 2 ES galaxies. None of these E and ES galaxies are not 'compact massive', that is, these 13 galaxies did not satisfy any of the four 'compact massive' selection criteria. The colour-coded dashed lines are the maximum number density of the red nuggets (labelled 'max(n RN )') from the respective high-publications which defined the size-mass selection criteria. Figure 19. The size-mass (in RC15 mass) distribution of 55 compact massive spheroids according to the size and mass selection criteria shown in Equation 11 (green dashed line and the shaded area). The host galaxies' morphologies are based on their reclassifications in Section 3.4. The sample is separated into 'S0 hosts' (upper panel) and 'S hosts' (lower panel). The red solid circles are the compact massive spheroids hidden within the host galaxies (labelled 'Hidden c,Sph') and the green plus signs are the host galaxies (labelled 'Host galaxies'). The size (circularised radius e,eq ) of the host galaxy is obtained via a single Sérsic fit to the surface brightness profile. The blue arrows indicate the direction of growth if the high-compact massive system were to assemble mass via the inside-out, disc cloaking process. We highlighted the three prominent outliers with red hollow circles. Figure 20. The spatial distribution of the galaxies. The parent sample selected by the angular boundaries from Equation 1 is marked in black dots. The 103 galaxies selected by the three bins (Equation 5) are marked as the blue circles. Among these 103 galaxies, the ones that are the members of the Virgo Cluster are marked as red circles.

APPENDIX A: DISTANCE CORRECTION
The variation among the four different distance measurements that we explored turns out not to be significant in regard to the number density of compact massive spheroids. We show these distance measurements in the forest plot of Figure A1. The W1997, M2000, Cosmicflow-3 (Kourkchi et al. 2020), and -independent distances are labelled in green, blue, purple, and orange points, respectively. The mean differences between each set of distances span 3.16 to 7.54 Mpc, while the standard deviations vary from 2.38 to 4.95 Mpc. The values and sources of the various -independent distances are listed in Table A1.
Cosmicflow-3 has the most sophisticated model to date, and it was built from the largest observational data set among the three model-based choices. The main issue of Cosmicflow-3 is the triplevalue region, where the gravity of the Virgo Cluster is significant enough to bend the velocity-distance relation into a cubic equation. Galaxies residing within 7-21 Mpc, therefore, have two or three solutions (possible distances) for the same recessional velocity. In Figure A1, we highlighted such objects with extended error bars (NGC 4374,NGC 4429,NGC 4459,NGC 4636,NGC 4643,NGC 4845,NGC 4649,and NGC 4772), in which the upper and lower limits represent the high and low estimation of the distance, respectively.
We used the redshift-independent distance as the primary distance to calculate the luminosity of the spheroids and galaxies in our sample. The only exceptions are NGC 4874 (an SBF distance), NGC 5440 (an SNIa distance), and NGC 4233 (an SBF distance), where all of them exhibit a large error (>10 Mpc). In these cases, we use the Cosmicflow-3 value instead. In total, 19 galaxies' total and spheroidal stellar mass are obtained via -independent distances. Other than these 19 galaxies, we primarily used the Cosmicflow-3 distance, in total, for 81 galaxies. There are an additional three galaxies that are within the triple-value region but which do not haveindependent measurements: NGC 4845, NGC 4643, and NGC 4429. In all three cases, the W1997 distances lie within the triple-value region, with a reasonable agreement with the median value of the Cosmicflow-3 distance. For these galaxies, we elect to use the W1997 distances, derived from empirical TF-relation, instead. The adopted distances for all 103 galaxies are provided in Table C1.

APPENDIX B: THE MASS-TO-LIGHT COLOUR RELATIONS
Here, we provide a supplementary description of possible error sources and the underlying assumptions that went into building the mass-to-light colour relations (MLCRs) presented in this work. This speaks directly to the stellar masses that we have derived.
Several key factors contribute to the error budget of the estimated stellar mass.
• The innate stellar population model degeneracy, between age and metallicity, is reported to constitute a small error of 0.1-0.2 dex from idealised mock galaxy studies (Gallazzi & Bell 2009).
• A more prominent error comes from the systematic treatments of the 'priors'. The IMF provides the normalisation factor for the MLCRs, and its uncertainty affects the * / ratio. For example, there can be a change of up to 0.3 dex in mass if one assumes a Chabrier (2003) over a Salpeter (1955) IMF (RC15).
• Another obvious factor is the stellar population model. The most widely-used prescription is the BC03 Stellar Population Synthesis (SPS) model. It provides a comprehensive library of the stellar population that was validated by agreement with early data from the SDSS. Over the years, the community has been continuously updating the model by including previously unconsidered stellar isochrones. The most notable issue plaguing stellar population models is the influence of the thermally-pulsating asymptotic giant branch (TP-AGB) stars. Multiple studies (e.g., Maraston 1998Maraston , 2005 have pointed out that the excess light from the young (∼1 Gyr) TP-AGB star's emission can lead to an underestimation of the * / ratio. This effect is particularly relevant for near-infrared (NIR) studies, as the inclusion of TP-AGB stars (Bruzual 2007a,b)  • Treatment of the interstellar medium (ISM) and dust attenuation also changes the perceived light. IP13 show the effect of omitting dust attenuation in their Figure 13. In a dusty scenario, the mass-to-light relation can change by up to 0.5 dex in the case of an inclined disc (see also Driver et al. 2007aDriver et al. , 2008.
Given the need to appreciate how the stellar masses of the highgalaxies were obtained by different authors, we identify and briefly describe each MLCR they used.
• Z09 assumes the Chabrier (2003) IMF, the CB07 SPS, and the dust attenuation from Charlot & Fall (2000). The * / ratio is produced using the fiducial mass reconstruction method by Monte Carlo sampling from the CB07 library, with parameters advised by Figure A1. Distances using different models. The pink zone is the area outside of the selection volume. The purple points are the Cosmicflow-3 (Kourkchi et al. 2020) distances, the blue points are the M2000 distances, and the green points are the W1997 distances. The orange points are a blanket label for the redshift-independent measurements from various sources; see Table A1. Galaxies are listed in descending (westward) order, according to right ascension. Note that several galaxies (NGC 5490, NGC 5406, NGC 4914, NGC 2824, NGC 5390, NGC 5354, NGC 5326, and NGC 5311) have their distances lie in different Bins. For those cases, we present the result in the Bins according to their Cosmicflow-3 distance as in both Table C1-C3 and Table D1-D3. Kauffmann et al. (2003). The * / ratio is then mapped onto the colour space of ( − ) and ( − ) to build their MLCR shown in Equation 3a. They verified this by comparing it to nine local galaxies within 30 Mpc.
• Both T11 and RC15 used the Chabrier (2003) IMF and the BC03 SPS, which, as we have seen, has been the most commonly used assumption in the literature when calculating the stellar masses of high-galaxies. Given the assumed dust attenuation law from Calzetti et al. (2000), T11 calculated the stellar mass by fitting the SED from early Galaxy And Mass Assembly (GAMA DR1) broadband photometry data via Bayesian probability minimization. Among the above-mentioned MLCRs, T11 is the only one that is built empirically from the ground up. The GAMA data covers the intermediate-redshift range ( < 0.65; median = 0.2) and magnitudes petro < 19.4-19.8 mag. However, claim has been made that T11's MLCR underestimated the / ratio for early-type spirals and ellipticals (Schombert et al. 2022).
• RC15 constructed their MLCR via mock galaxies using the MAGPHYS library (da Cunha et al. 2008). The mock galaxy SEDs undergo Bayesian fitting on bands to recover the stellar masses. The MLCR yields excellent agreement with the observational SHIVir survey ('Spectroscopy and -band imaging of Virgo Cluster galaxies'; McDonald et al. 2011) data, having mean residuals of −0.01 dex and 0.05 dex (see their Figure 6).
With the same IMF and SPS as T11, RC15 is different in several ways. First, RC15 assume a different dust treatment (Charlot & Fall 2000) with two-component models accounting for the dust in molecular clouds and the ambient ISM. Second, RC15 suggests that T11 did not include bursts in their model SFHs, while MAGPHYS contains random burst models superimposed in the SFH. Both works are validated by their respective observational data: one of which focuses on the nearby Virgo Cluster galaxies and the other on galaxies at intermediate redshifts. The considerable discrepancy between them ( Figure 2) might arise from the use of different data sets.
• Finally, IP13 constructed their MLCR uniquely by choosing a Kroupa (1998) IMF and the Marigo et al. (2008) isochrone model (the same used in CB07). Once again, the choice of Marigo et al. (2008) isochrones stems from the concern over the TP-AGB influence. Their equation is constructed assuming a dust-free scenario. While the effect of dust is significant in some colour, they concluded that ( − ) remains a good tracer for the galaxies' stellar masses. The consequence of using such a combination is the slower decline of stellar mass as a function of age (see their Figure 1) in comparison to the more popular BC03 SPS.

APPENDIX C: GALAXY DATA
Here, we present the basic information for the galaxy sample identified in Section 2.4. Tables C1-C3 also provide the new morphological classifications (Column 9) based on our multi-component decomposition, as described in Section 3.3.
In Tables C1-C3, Column (8) lists the past morphology assignment while Column (9) provides our new classifications. The labels for the new classes are inspired by the 'morphology grid' from Graham (2019) for high surface brightness galaxies. It represents the progression from the traditional Hubble-Jeans sequence (Jeans 1919;Lundmark 1925;Hubble 1926aHubble ,b, 1927Hubble , 1936, with an emphasis on the commonly-overlooked features in early-type galaxies, such as bars and the continua of disc sizes (nuclear, intermediate, and large-scale). The grid includes the ellicular (ES) class of galaxies (Liller 1966), denoting the spheroids with an embedded intermediate-scale disc, Figure C1. The distribution of galaxy circularized effective radius ( e )maximum radius ( max ) ratio. We obtained e via a single Sérsic fit and the max is our cutoff radius for the fit. The sample is separated by their morphological type and colour coded with red, orange, and blue for E, S0, and S galaxies, respectively. Here, we only show the distribution between 0 < e / max < 3.0, where 96 galaxies are present. and it removes the redundant E4-E7 classes because they are essentially all lenticular galaxies (Liller 1966;Gorbachev 1970;Michard 1984;Capaccioli et al. 1990b;van den Bergh 1990). For simplicity's sake, we omit the subdivision of spiral arm shapes (a, b, c, d, and m) which are not as relevant to our investigation. We label the galaxies as ellipticals without a disc (E0-E3), elliculars (EAS or EBS), lenticulars (SA0, SAB0, or SB0), and spirals (SA, SAB, or SB), in which A and B are the distinctions between barless and barred galaxies, and AB is designated for weak bars, defined here as those with lower central surface brightness ( 0 ) than the disc.
In accordance with the common practice in the literature, we performed a single Sérsic function fit to determine the effective radius ( e ) for each galaxy. In Figure C1, we show the distribution of e / max for our galaxies, separated by their morphology. We choose the cutoff radius, max , to fitting the light profile at the radius where the / ratio is low and no longer useful for decomposition. The majority (83/103) of our galaxies have e / max < 1. There are 20 galaxies which have effective radii exceeding their cutoff radii ( e / max > 1). They are exclusively S0 and S galaxies. The outliers illustrate the problem of using a single Sérsic fit to measure the size of a galaxy. A shallow disc profile combined with a prominent bulge can result in a long-tailed surface brightness profile. The parameterised effective size of such a galaxy, from a single Sérsic fit, will be an overestimate of the actual half-light radius. The effect of inflated galaxy radii in S0 and S galaxies is not widely appreciated and shall be addressed in a separate work.

APPENDIX D: LOCAL SPHEROIDS' PARAMETERS
We present the structural parameters and the masses of the bulges/spheroids in our host galaxies in Tables D1-D3. Column (2) is the distance we used to calculate the absolute magnitude and stellar masses (see Appendix A and (iv) in Section 2.4). Column (3) showcases the function we used to model the spheroids, either a Sérsic (S) or a core-Sérsic (cS) function (see Graham et al. (2003a). Columns (4)-(6) provide the parameters of the Sérsic and core-Sérsic functions. In column (4), it represents the surface brightness (i.e. e ) at effective radius e (column (5)) in equivalent axis. Column (6) shows the Sérsic index of the spheroids. Columns (7) and (8) give their apparent and absolute magnitudes, respectively. Columns (9)- In descending (westward) order, according to right ascension. Columns: (1) the galaxy name; (2) the right ascension angle in J2000 coordinates (in degrees); (3) the declination angle in J2000 coordinates (in degrees); (4) the adopted (luminosity) distance described in part (iii) of Section 2.4; (5) the -band apparent magnitude and (6) the -band apparent magnitude from the NASA-Sloan ATLAS catalogue (http://www.nsatlas.org/) after correcting both for Galactic dust extinction and the K(z)-correction; (7) the seeing (in arcsec) measured by imexam; (8) the RC3 Hubble morphology classification; (9) our new classification described in Section 3.4; (10) the absolute -band magnitude calculated using columns (4) and (6), with Galactic dust and K-correction from the NASA-Sloan ATLAS catalogue applied; and (11) the total galaxy stellar mass using the IP13 ( − )-dependent * / ratio applied to the SDSS-derived -band magnitude. (12) the cutoff radius (in arcsec) in circularised equivalent-axis.
(12) are the spheroids' stellar masses based on the four * / ratios discussed in Section 2.3. This paper has been typeset from a T E X/L A T E X file prepared by the author. Columns: See Table C1.