LeMMINGs. V. Nuclear activity and bulge properties: a detailed multi-component decomposition of $e$-MERLIN Palomar galaxies with $HST$

[Abridged] We use high-resolution $HST$ imaging and $e$-MERLIN 1.5-GHz observations of galaxy cores from the LeMMINGs survey to investigate the relation between optical structural properties and nuclear radio emission for a large sample of galaxies. We perform accurate, multi-component decompositions of new surface brightness profiles extracted from $HST$ images for 163 LeMMINGs galaxies and fit up to six galaxy components (e.g., bulges, discs, AGN, bars, rings, spiral arms, and nuclear star clusters) simultaneously with S\'ersic and/or core-S\'ersic models. By adding such decomposition data for 10 LeMMINGs galaxies from our past work, the final sample of 173 nearby galaxies (102 Ss, 42 S0s, 23 Es plus 6 Irr) with bulge stellar mass (typically) M_*, bulge ~ 10^6-10^12.5 M_sun, encompasses all optical spectral classes (LINER, Seyfert, ALG and H II). We show that the bulge mass can be significantly overestimated in many galaxies when components such as bars, rings and spirals are not included in the fits. We additionally implement a Monte Carlo method to determine errors on bulge, disc and other fitted structural parameters. Moving (in the opposite direction) across the Hubble sequence, i.e., from the irregular to elliptical galaxies, we confirm that bulges become larger, more prominent and round. Such bulge dominance is associated with a brighter radio core luminosity. We also find that the radio detection fraction increases with bulge mass. At M_*,bulge>10^11 M_sun, the radio detection fraction is 77%, declining to 24% for M_bulge<10^10 M_sun. Furthermore, we observe core-S\'ersic bulges tend to be systematically round and to possess high radio core luminosities and boxy-distorted or pure elliptical isophotes.


Introduction
In the standard structure formation model, elliptical galaxies and massive bulges 1 of lenticular galaxies (S0s) seen today Dutton & van den Bosch 2009;Hopkins & Quataert 2011;Anglés-Alcázar et al. 2021) gives rise to star formation, the growth of supermassive black holes (SMBHs), and active galactic nuclei (AGN) activity. Feedback from an accreting SMBH, which is often evident in the form of radio jets and outflows, particularly for the more common, lower luminosity AGN and low-ionisation nuclear emission-line regions (LINERs) which dominate the local Universe, is believed to inject energy and momentum that expel gas and regulate (the mass of the SMBH itself, e.g. Debuhr et al. 2010) and the star formation histories of the host galaxy (e.g. Silk & Rees 1998;Murray et al. 2005;Fabian et al. 2006;Robertson et al. 2006;Croton et al. 2006;Bower et al. 2006;Vogelsberger et al. 2014;Weinberger et al. 2017). How such a process impacts the structural properties of the host galaxy (e.g. stellar light distribution, stellar mass, concentration, and disc-to-total light ratio), remains a key unanswered question. The empirical scaling relations between the SMBH mass and the properties of its host bulge are understood as an indication of co-evolution of SMBHs and their host bulges (e.g. Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Graham et al. 2001;McLure & Dunlop 2002;Marconi & Hunt 2003;Häring & Rix 2004;Gaspari et al. 2019;Dullo et al. 2020).
The effect of AGN feedback has been invoked to explain the presence of two types of structurally distinct elliptical galaxies: (i) core-Sérsic ellipticals with depleted cores whose (massive, n 4) bulges likely experienced 'dry' (gas-poor) merger events involving binary SMBHs, and (ii) Sérsic ellipticals which are coreless, with less massive (n ∼ 3 ± 1) bulges thought to be built in gas-rich processes (e.g. Faber et al. 1997;Somerville et al. 2008). Kormendy et al. (2009) argued that X-ray emitting, hot gas in core-Sérsic ellipticals is prevented from cooling and forming stars by radio-mode AGN feedback, explaining the dry formation scenario. In contrast, Sérsic ellipticals undergo nuclear starburst events since the AGN feedback is weaker and the galaxies' shallow potential well is incapable of retaining hot gas (e.g. Cattaneo et al. 2009). While the influence of the AGN activity on the host galaxy structures is not yet clear, a growing body of observational evidence points to heating by (episodic) radio AGN from central galaxies to prevent the intracluster medium from cooling in galaxy groups and clusters (Sijacki et al. 2007;Somerville et al. 2008;McNamara et al. 2009McNamara et al. , 2011Fabian 2012). High-resolution, multi-wavelength data are crucial to investigate the central structures of galaxies, AGN activities, and jet structures.
The Legacy e-MERLIN Multi-band Imaging of Nearby Galaxies Survey (LeMMINGs; Beswick et al. 2014;Baldi et al. 2018Baldi et al. , 2021a is designed to reveal the physical origin of the observed relationships between different wavebands by combining high-resolution observations from radio (e-MERLIN), through optical -Hubble Space Telescope (HST) -to X-ray (Chandra). As part of LeMMINGs, we have recently completed e-MERLIN 1.5 GHz observations of all 280 galaxies above declination δ > +20 • from the Palomar bright spectroscopic sample of nearby galaxies (Filippenko & Sargent 1985;Ho et al. 1995Ho et al. , 1997a. Using our e-MERLIN 1.5 GHz data in Baldi et al. (2021b), we reported different radio production mechanisms for the different optical classes and SMBH masses: above M BH ∼ 10 6.5 M AGN are the main driver of the nuclear radio emission in galaxies, whereas below this characteristic mass stellar processes power the majority of the radio emission. Furthermore, our new 5 GHz e-MERLIN data for the full LeMMINGs sample of galaxies have also been calibrated and are currently in the process of being analysed (Williams et al., in prep.). To allow for a proper understanding of the relationships between the radio properties revealed by e-MERLIN and the structural properties of the host galaxies as revealed by optical observations, we require optical imaging with similar spatial resolution to the e-MERLIN observations. Only HST and the James Webb Space Telescope (JWST) can provide such observations. The main aim of this paper is to provide high-quality host galaxy optical structural parameters and then to investigate the relationship between those parameters and the radio parameters.
To understand the true nature of these correlations between optical galaxy structural parameters and radio core properties, it is necessary first to identify and model the distinct photometric components of the galaxies such as bulges, discs, bars, rings, spiral arms, AGN, nuclear star clusters (NSCs), depleted cores, and haloes. Clues as to how bulges form and evolve are thought to have been imprinted in the shapes of the galaxies' surface brightness profiles, which exhibit a degree of diversity. Simulations have shown that multiple major mergers would naturally produce a spheroidal structure with a large Sérsic index n 4.0 (e.g. Nipoti et al. 2006;Naab et al. 2006;Hopkins et al. 2009a,b). This is contrary to galaxies that underwent dissipative, unequal mass mergers or have been shaped by secular processes, as they possess a bulge with a low Sérsic index (n 3.0). To measure robust bulge structural properties, and investigate their connection with sub-kiloparsec radio core properties, the bulge and the remaining galaxy components should be carefully separated by performing physically motivated, multi-component decomposition (e.g. Laurikainen et al. 2005Laurikainen et al. , 2010Gadotti 2008Gadotti , 2009Salo et al. 2015;Saglia et al. 2016;Méndez-Abreu et al. 2017;Sahu et al. 2019;Davis et al. 2019;Dullo & Graham 2014;Erwin et al. 2015;Dullo et al. 2016Dullo et al. , 2017Dullo et al. , 2019Kim et al. 2017;Gao & Ho 2017;Dullo 2019;Su et al. 2021). When fits are restricted to account only for the two main galaxy components, that is a bulge-disc profile, as is commonly the case for fits performed on a large sample of galaxies in an automated fashion (e.g. Allen et al. 2006;Simard et al. 2011), the flux contribution from prominent other components can wrongly brighten the bulge luminosity and give rise to erroneous central and global structural parameters for the galaxy (e.g. Laurikainen et al. 2010;Salo et al. 2015;Dullo et al. 2016Dullo et al. , 2019Su et al. 2021). This can lead to wrong conclusions about key correlations such as those between bulges and SMBHs and radio core properties.
Poor spatial resolution has also been shown as a key limitation in deriving galaxy structural properties. Structural studies that rely on ground-based data (e.g. Laurikainen et al. 2005Laurikainen et al. , 2010Méndez-Abreu et al. 2017) and Spitzer Space Telescope data (e.g. Salo et al. 2015;Sahu et al. 2019;Davis et al. 2019) do not have the appropriate resolution to study sub-kiloparsec structures, such as depleted cores, AGN, NSCs, nuclear discs, and inner bars, which are commonly detected at HST resolution.
To date, HST studies of nearby galaxies have mainly focussed on early-type galaxies (Lauer et al. 1995(Lauer et al. , 2007aFaber et al. 1997;Carollo et al. 1997a;Rest et al. 2001;Ravindranath et al. 2001;Laine et al. 2003;Trujillo et al. 2004;Ferrarese et al. 2006;Kormendy et al. 2009;Richings et al. 2011;Turner et al. 2012;Dullo & Graham 2012, 2013, 2014Dullo et al. 2021;Krajnović et al. 2013Krajnović et al. , 2020Rusli et al. 2013;Dullo 2019). Multi-component structural analysis of HST images for a large sample of late-type galaxies is lacking. Carollo et al. (1997bCarollo et al. ( , 1998 studied 40 spiral galaxies using HST WFPC2 and NICMOS images. These authors used the inner (radially limited) 10 HST light profiles and restricted their fits to describe only the bulge light distribution using R 1/4 , R 1/4 +exponential, exponential or double-exponential models. Fisher & Drory (2008) modelled the HST light profiles of 64 spiral galaxies using a two-component Sérsic bulge + exponential disc model. The robustness of these fits is limited as innermost data points (R 1 ), bars, rings, and spiral arms were excluded from the fits after being subjectively identified through visual inspections.
Over the past two decades a few studies measured galaxy structural properties with HST and explored the relation with radio-quiet and radio-loud dichotomy (e.g. Capetti & Balmaverde 2005, 2007de Ruiter et al. 2005;Baldi et al. 2010;Richings et al. 2011). Capetti & Balmaverde applied the Nuker model, which does not provide an accurate description of galaxy surface brightness profiles Ferrarese et al. 2006;Dullo & Graham 2012), to radially limited 10 HST light profiles of early-type galaxies, resulting in an approximate bulge structural classification ('core' versus 'power-law'). Subsequent work by Richings et al. (2011; see also Baldi et al. 2010) fit a Sérsic, core-Sérsic, and double-Sérsic model to the HST surface brightness profiles of 110 galaxies with Hubble type earlier than Sbc, selected to avoid disc dominated systems. They reported that their fits were robust for 62/110 galaxies. However, there are limitations. Most of their fits (>75%) rely on radially limited (R 20 ) NICMOS light profiles. Also the fits do not account for large-scale galaxy components such as bars and rings which potentially contribute significant light to the inner domain of the galaxy surface brightness profiles.
In this paper we perform accurate multi-component decompositions of the one-dimensional HST surface brightness profiles of a representative sample of 173 active and inactive galaxies from the full sample of 280 LeMMINGs galaxies (see also Dullo et al. 2016Dullo et al. , 2018Dullo et al. , 2019Dullo 2019). Uniquely, our sample contains a large number of 108 late-type galaxies. We fit up to six galaxy components (e.g. bulges, discs, depleted core, AGN, nuclear star clusters, bars, spiral arms, rings, and stellar haloes) simultaneously to the galaxy surface brightness profiles, which extend out to large radii of R 80−100 , adequate to accurately quantify the shapes of the bulge and disc profiles along with other large-scale galaxy components (e.g. bars, spiral arms, rings, and stelar haloes). Our work, which represents the largest, most detailed structural analysis of nearby galaxies with HST to date, offers the possibility to make a plausible connection between the optical structural properties of galaxies and their nuclear optical and radio activities in a homogenous manner and at a much higher resolution than has hitherto been possible. Dullo et al. (2023) investigate the connection between the radio core luminosity and various properties of the host bulge.
The outline of this paper is as follows. In Sect. 2, we discuss the LeMMINGs sample, radio and optical emission line data, our HST imaging data and data reduction details. Section 3 presents an overview of the analytical models used to describe galaxy components and then discusses the 1D (one-dimensional) and 2D (two-dimensional) multi-component decompositions of the galaxy stellar light distributions. This section also describes our colour calibration equations we derived to transform the magnitudes obtained in various filters into V-band and the stellar massto-light ratios employed to calculate the stellar masses from the luminosities. In Sect. 4, we study the connection between bulge structural properties, optical emission-line classes and nuclear radio activity. Finally, in Sect. 5 we summarise our main results and conclude. Notes on our decompositions and comparison with past fits in the literature for 54 selected individual galaxies are given in Appendix A. The data tables, individual galaxy HST images, 1D and 2D decompositions and our results are presented in Appendices B-E, respectively. The Monte Carlo (MC) method which we employ to estimate realistic errors on fitted structural parameters and to test the robustness the 1D multi-component decomposition of the galaxy light profiles, is described in Appendix F.

The LeMMINGs sample and data
The LeMMINGs (Beswick et al. 2014;Dullo et al. 2018;Baldi et al. 2018Baldi et al. , 2021a is the deepest, high-resolution radio survey of the local Universe and constitutes radio continuum observations of 280 galaxies for a total of 810 h with e-MERLIN at 1.25−1.75 GHz (L-band) and C-band, centred at 5 GHz, with a 512 MHz band width, complemented by synergies from highresolution, multi-band (X-ray, optical and near-IR) data. The survey sample is a subset of the magnitude-limited (B T ≤ 12.5 mag and declinations δ > 0 • ) Palomar spectroscopic sample of 486 bright, nearby galaxies (Ho et al. 1995(Ho et al. , 1997a, which in turn were drawn from the Revised Shapley-Ames Catalog of Bright Galaxies (Sandage & Tammann 1981) and the Second Reference Catalogue of Bright Galaxies (de Vaucouleurs et al. 1976). In order to discriminate between star-formation and AGN activities in Palomar galaxies, Ho et al. (1995Ho et al. ( , 1997a utilised optical emission-line ratios and classified the emission-line nuclei into four spectral classes: Seyferts, LINERs, H ii, and Transition galaxies. For LeMMINGs galaxies, Baldi et al. (2018Baldi et al. ( , 2021a used the diagnostic line ratios mostly from Filippenko & Sargent (1985), Ho et al. (1995Ho et al. ( , 1997a plus new spectroscopic classifications from the literature and revised the classification for the entire sample employing the emission line diagnostic diagrams by Kewley et al. (2006) and Buttiglione et al. (2010). In this revised classification, used in this paper, each object with emission lines was categorised as H ii, Seyfert, or LINER after Transition galaxies (Ho et al. 1995) were identified either as LINER or H ii. There are also galaxies in the LeMMINGs sample which lack emission lines and were classified as Absorption Line Galaxies (ALGs), Baldi et al. 2018Baldi et al. , 2021a To allow for better radio visibility coverage for the e-MERLIN array, all LeMMINGs galaxies were selected to have δ > +20 • (see Fig. 1). An extensive description of this legacy survey goals, radio data reduction technique, radio detection, maps, and flux measurements is presented in Baldi et al. (2018Baldi et al. ( , 2021a. With an angular resolution of ∼0 . 15 and sub-mJy sensitivity afforded by e-MERLIN at 1.5 GHz, we were able to measure radio core luminosities L R,core ∼ 10 34 −10 40 erg s −1 . Furthermore, the analysis of Chandra X-ray observations for 213 LeMMINGs galaxies is presented in Williams et al. (2022).
In this work we perform accurate multi-component analysis of the central and global structures for LeMMINGs/Palomar galaxies with HST imaging data (Tables 1 and 2). Our sample comprises 173 galaxies (∼62%) from the full sample of 280 LeMMINGs galaxies (Tables 1 and 3). As the best balance between image availability and immunity to the obscuring effect of dust, preference was given to HST images taken with the redder optical band filter F814W. We did not use HST images taken in narrow-band or mid-UV (F300W) filters. Of the 107 LeMMINGs galaxies excluded in this work, 31 have peculiar shown in red are the 173 LeMMINGs galaxies with decent HST imaging studied in this paper, the remaining 107 (shown in blue) are galaxies either (a) with centres obscured by dust, (b) that belong to a pair of merging/interacting galaxies, (c) with peculiar and/or edge-on morphology, (d) with HST images that miss the galaxy centre, or that were taken in narrow-band or mid-UV (F300W) filters only and (e) with no HST observations. Right panel: a plot of the radio core luminosity (L Radio,core ) as a function of total galaxy B-band absolute magnitude from Hyperleda (M B,glxy ). The sample studied in this work is representative of the full LeMMINGs sample which is constructed by selecting all the galaxies in Palomar spectroscopic sample (Ho et al. 1995(Ho et al. , 1997a with δ > 20 • . Notes. The LeMMINGs galaxies that are omitted are those (a) with centres obscured by dust, (b) that belong to a pair of merging/interacting galaxies, (c) with peculiar and/or edge-on morphology that render the extraction and modelling of the galaxy light profiles difficult and thus deferred to a future work, (d) with HST images that miss the galaxy centre, or that were taken in narrow-band or mid-UV (F300W) filters only and (e) with no HST observations. and/or edge-on morphologies that make the extraction and modelling of the galaxy light profiles difficult and thus deferred to a future work, five exhibit centres markedly obscured by dust (See Fig. 2), two belong to a pair of merging/interacting galaxies, 16 have HST images which miss the galaxy centre or were taken in narrow-band or mid-UV (F300W) filters and 53 lack HST observations (see Table 1). We rely on newly extracted HST surface brightness and associated geometrical profiles for 163/173 LeMMINGs galaxies, augmented by the HST structural data from our previous work for the remaining 9 galaxies (Dullo & Graham 2014;Dullo et al. 2016, 2018, see Fig. 3 andAppendix B).
As one of the goals of this paper is to find out how radio luminosity is related to optical properties, it is important to investigate the parameter space coverage for our sample of 173 LeMMINGs galaxies and ensure that the sample is representative of the full LeMMINGs sample. As noted above the par- Notes. The FOV for the SDSS images used in this work is ≥9 .0 × 9 .0. The science images retrieved from the HLA are (re)sampled and aligned north up using the DrizzlePac (Gonzaga et al. 2012) andMultiDrizzle (Koekemoer et al. 2003) software packages. For the 104 sample galaxies, we used the full WFPC2 mosaic images. Notes. The sample galaxies are first separated based on the galaxy morphological and optical spectral classes (Cols. 1-2) and then further divided based on their radio non-detection and core-Sérsic type central structure (Cols. 3-4). The term 'full sample' refers to the total LeMMINGs sample of 280 galaxies, whereas the term 'our sample' refers to the sub-sample of 173 LeMMINGs galaxies studied in this paper.

Fig. 2.
Central 25 × 25 regions of three LeMMINGs galaxies, taken with HST, which have their centres copiously obscured by dust (i.e. dark areas). Consequently, the galaxies were excluded from our analysis as the extraction of surface brightness profiles was not possible.
ent Palomar sample is statistically complete and so is the full LeMMINGs sample as it represents all the Palomar galaxies with δ > 20 • . Figure 1 shows the distribution for our 173 galaxies and the remaining 107 LeMMINGs galaxies. The twodimensional Kolmogorov-Smirnov (KS) test (Peacock 1983;Fasano & Franceschini 1987) on these two sub-samples ( Fig. 1 Fig. 3. Histograms of distance (mostly from NED), absolute B-band galaxy magnitude and from HyperLeda (M B,glxy ) and morphological Ttype shown for the subset of 173 LeMMINGs galaxies in this work and the full, statistically complete sample of 280 LeMMINGs galaxies. The Hubble classes based on the numerical T-type and absolute B-band galaxy magnitudes are from HyperLEDA. left) yields a significance level P ∼ 0.4, thus we cannot reject the null hypothesis that the sub-sample are drawn from the same distribution. In Fig. 1 (right) we also show a plot of the radio core luminosity (L R,core ) as a function of total galaxy B-band absolute magnitude (M B,glxy ) from Hyperleda (Makarov et al. 2014). It is evident that our sample provides good coverage of the entire range of radio core luminosity and galaxy luminosity afforded by the full LeMMINGs sample. We also employ the 2D KS test on the (L R,core , M B,glxy ) dataset for the sub-sample used in this work and the full LeMMINGs sample, finding a significance level of P ∼ 0.3 which indicates that the null hypothesis, that is the (L R,core , M B,glxy ) data for the two samples are drawn from the same distribution, cannot be rejected and that our sample is representative (albeit being statistically incomplete). Furthermore, Fig. 3 also shows that our sample covers a wide range in distance (D ∼ 0.7−100 Mpc) and morphological type between E and Im (Hubble 1926;de Vaucouleurs 1959). These important galaxy property considerations make the sample well suited for robust characterisation of galaxy structural scaling relations. The mean distance for our sample is ∼22 Mpc, which translates to spatial scales by HST of 5−10 pc. Details of our sample are summarised in Tables B.1−B.5.

Archival HST and ground-based imaging data
High-resolution HST images for the 163 target galaxies were retrieved from the Hubble Legacy Archive (HLA 2 ). These images were processed using the HLA pipeline which performs the standard data reduction procedures including bias subtraction, geometric distortion correction, dark current subtraction, and flat fielding. The galaxy images were observed with the HST ACS, WFPC2, NICMOS, or WFC3 cameras and taken in the following filters: F547M (equivalent to the Johnson-Cousin Vband), F555W (∼V-band), F606W (broad ∼V-band), F702W (Cousins ∼R-band), F814W (∼I-band), F850LP (∼SDSS-z band), F110W (∼J-band), and F160W (∼ H-band).
The fitted radial extent is of critical importance when decomposing a galaxy's surface brightness profile. The recovery of accurate model parameters in particular for large-scale (∼3.5 ± 2.0 kpc) galaxy components such as discs, rings, and spiral arms is prone to how well the curvature of the light profile 3 that represents the outer most galaxy components is sampled by the data. 2 HLA is the site for high-level HST archival products at https:// hla.stsci.edu 3 The light profile of a galaxy is the radial distribution of its stellar light.
A related concern is the uncertainty in the shapes of the 1D light profiles of the outer galaxy components due to sky background errors. HST images do not always probe a sufficiently large radial domain to perform a firm assessment of the background level and to sample the outer parts of the surface brightness profiles for nearby galaxies. For 81/173 (46.8%) of the sample galaxies their large-scale structures (e.g. discs, rings, and spiral arms or low surface brightness stellar envelopes) extend beyond the FOV of the high-resolution HST images. Measuring accurate structural parameters for such galaxies may necessitate the high-resolution HST data at small radii to be matched and combined with radially extended ground-based data at large radii (see e.g. Dullo et al. 2016Dullo et al. , 2018Dullo 2019). To extract 1D light profiles we therefore use ground-based SDSS-r, i, z, and DSS (103aE, i.e. optical red) band images 4 for the galaxies, which are downloaded from the SDSS (York et al. 2000) Data Release 13 (DR13 5 , Albareti et al. 2017) and NED 6 , respectively (see Table 2 and Sect. 2.4). A summary of the HST imaging, instrument, and filter is presented in Tables B.2-B.4, respectively.
The extraction of accurate surface photometry for the sample galaxies proceeds in several steps: sky background subtraction (Sect. 2.2), creation of mask files (Sect. 2.3) and then fitting elliptical isophotes to the galaxy images using the IRAF task ellipse and for some sample galaxies the creation of composite surface brightness, ellipticity, PA, and B 4 profiles by matching HST and large FOV ground-based data (Sect. 2.4).

Sky background subtraction
As noted above, 81/173 sample galaxies extend beyond the FOV of the high-resolution HST images, as such background oversubtraction by the HLA pipeline is a potential concern. Indeed, we find that, because of sky over-subtraction by the HLA, for most angularly extended galaxies in the sample the HST-derived brightness profiles tend to depart downwards with respect to the SDSS profiles at radii of R 50-60 (or R 10 for NICMOS images. Therefore, for such galaxies we extract light profiles from the large FOV SDSS/DSS images and matched them with our inner HST data (as detailed in Sect. 2.4) to construct composite (i.e. inner HST plus outer ground-based) light profiles. The SDSS images used for the 81 sample galaxies have ≥9 .0 × 9 .0 FOV, considerably larger than the median major-axis diameter for our galaxies which is 2 .6 ± 2 .0. This allowed us to characterise the background level for our 1D and 2D decompositions more accurately.
For the remaining 92 sample galaxies, the quality of the HLA pipeline sky subtraction is examined by computing manually the median sky levels from several (>5) 20 × 20 pixel boxes at the edges of the HST CCDs. We find that the galaxy flux in the outermost parts is 10% of the HST sky level, that is ∼21.42 ± 0.33 mag arcsec −2 I-band, which converted here from a V-band value (Lauer et al. 2005) using V − I = 1.08 (Fukugita et al. 1995). We ensure that the average of the medians is near zero (∼0.004 ± 0.004 e s −1 ) when the galaxies are within the field of view of the HST CCDs.

Masking of images
We inspected the galaxy images visually for dust structures, gaps between individual HST CCD detectors, image defects, bright foreground stars and background galaxies, which were then masked by taking extra care as discussed here. For each galaxy, an initial mask region file was first created using SExtractor (Bertin & Arnouts 1996). To achieve this, we used a threshold background value, detection threshold of (typically) 2σ above the background and convolved the galaxy's image with a Gaussian filter gauss_2.0_3×3 to generate catalogue sources with their image coordinates. This was first converted to a DS9 7 region file and then into a mask file using iraf. The initial mask was then added to a manual mask. In some cases we first removed a model galaxy image, defined in terms of the isophotal parameters derived for the science image, created by the iraf task bmodel. This was crucial to better identify dust-obscured regions, contaminating light sources and thus improve the initial mask. Finally, we run the iraf mskregions task to create mask files with '.pl' extension, automatically recognisable by the iraf task ellipse (see Sect. 2.4).

Surface brightness profiles
The extraction of major-axis surface brightness profiles for our galaxies was carried out following data reduction steps similar to those described in Dullo et al. (2016Dullo et al. ( , 2017Dullo et al. ( , 2019, Dullo (2019).
The iraf task ellipse was used to fit elliptical isophotes to the sky-subtracted galaxy images along logarithmically increasing semi-major axis (Jedrzejewski 1987; see also Kent 1983), hence giving the most weight to the innermost portion of the galaxy. We used 3σ clipping for flagging deviant pixels from contamination by cosmic rays and image defects. ellipse samples the 2D intensity distribution in the images as a function of an azimuthal angle φ starting from a first guess elliptical isophote defined by initial values of the isophote centre, position angle (PA), and ellipticity, = 1 − b/a, where a and b are the semi-major and semi-minor axes of the isophote, respectively. While and PA were always set free to vary, the isophote centres were held fixed for most galaxies. When concentric isophotes are assumed for a galaxy that is actually well modelled by a nested isophotes and has a centre that wobble by few pixels it only slightly affects the galaxy surface brightness profile and the associated best-fitting structural parameters for the galaxy remain unchanged. For eight galaxies (NGC 2770, NGC 2782, NGC 3034, NGC 3448, NGC 4605, andNGC 5474) with a low signal-to-noise ratio and complex morphology, the extraction of isophotal parameters and surface brightness profiles with ellipse was difficult. To overcome the problematic nature of such galaxies we ran ellipse multiple times using several sets of conditions: different initial semi-major axis lengths, occasionally fixing the ellipticity and PA, increasing the logarithmic spacing between successive ellipses and in a few cases a residual image was created after a model image was built from ellipse table ('.tab') file. The surface brightness profiles were then extracted when convergence is achieved with ellipse and the residual images do not show strong galaxy light. We classify the fits to these eight galaxies' surface brightness profiles as Quality 2 (see Sect. 3.2). 7 https://ds9.si.edu, SAOImageDS9.
A Fourier expansion of the intensity distribution with an average intensity I 0 can be written as (1) Higher order coefficients of the Fourier series ( j ≥ 3) carry crucial information as they quantify any deviations of the isophotes from perfect ellipses. Of particular importance is the coefficient of the cos(4φ) term namely the fourth-order moment B 4 : a positive B 4 value indicates the isophote is 'discy', whereas negative B 4 values signify that the isophotes are 'boxy'. For edge-on disc galaxies, Ciambur (2015) implemented modifications to iraf/ellipse task and introduced a new fitting routine dubbed isofit to better represent the galaxy isophote shapes. While this latter task was not employed in this work to extract the light profiles for the 13 (7%) sample galaxies with high inclination as determined from the galaxy axial ratio i 75 • , our fit quality flag, which quantifies the reliability of the 1D light profile decompositions, reflects the effect of inclination (see Table B.2). Nonetheless, we find a small dispersion on the best-fitting parameters, comparing, for each of these 13 galaxies, the decompositions of the original 1D light profile and its (≥100) simulations (see Sect. 3.2 and Appendix F). Also, we perform 2D decompositions for six of the 13 high inclination galaxies, finding good agreement between the 1D and 2D modelling (see Sect. 3.2.5 and Appendices D and E). As such, we decided to include all the 13 galaxies in our analyses. In the future, we plan to fit a 2D, edge-on exponential disc model using imfit (Erwin 2015) and remodel some of the high inclination sample galaxies to further explore the robustness of the galaxies' structural parameters. Likewise, for the 1D decompositions, isofit will be used to extract new light profiles for the galaxies.
In constructing composite surface brightness, ellipticity, PA and B 4 profiles for the 81 sample galaxies for which we used ground-based imaging to extend the HST FOV, see Sect. 2.1, we manually shifted the ground-based profiles to match HST data by applying constant zero-point offsets. With the exception of two galaxies (NGC 3198 and NGC 3665), we could ensure that the HST and ground-based images from SDSS and DSS for the individual galaxies are taken with filters that are similar or close (e.g. HST F814W and SDSS-i; HST F606W and SDSS-r), thus tracing similar stellar populations with no obvious colour gradients. For NGC 3198 and NGC 3665, data from two different filters (NICMOS F160W and SDSS-z) are combined. The matching procedure resulted in excellent overlap between the two sets of data over a large radial range (i.e. typically R ∼ 2 −50 ), but the range is smaller when HST NICMOS and ground-based images are used (R ∼ 2 −10 ). We find poor sky subtraction by the HLA pipeline has negligible effect on our HST-derived light profiles at R < 50 (R < 10 for the NICMOS data). We note that the ground-based SDSS and DSS profiles depart downwards with respect to the HST profiles inside R ∼ 1 where the blurring from seeing in the ground-based images becomes important.
Overall, the galaxy light profiles extend out to semi-major axis R 80−100 , covering 2R e,bulge for 97% of the sample. The light profiles provide large ranges in surface brightness down to 26.0 mag arcsec −2 in I-band with a 1σ error of 0.14 mag arcsec −2 . Furthermore, we follow similar steps as in Pohlen & Trujillo (2006) and determine from SDSS-i images the limiting surface brightness (µ lim ) as the surface brightness in which the galaxy profiles perturbed by 3σ sky deviate by 0.2 mag from the actual galaxy light profiles, where σ sky is the uncertainty on the sky level measured using 20, 100 × 100 pixel boxes. This yields µ lim,I ∼ 24.0 mag arcsec −2 . For the sample galaxies with SDSS-i data, the median surface brightness value of the outermost data points in the light profiles is µ med,I 24.0 ± 0.5 mag arcsec −2 . This suggests that the main bodies of the large-scale galaxy components such as bars, discs, rings, and spiral arms in the light profiles are brighter than the surface brightness limit. That is, owing to the large radial extent of the light profiles, we are able to characterise and quantify the properties of bulges and constrain the shape of the stellar light distribution that describes large scale galaxy components. We quote all magnitudes in the VEGA magnitude system.

Light profile models
The analytic description of the stellar light distribution of galaxies provides a means to quantify their photometric structures and to examine how these structural properties change across galaxies, wavelength and environment. A number of analytic functions have been used to describe the stellar light distributions in early-and late-type galaxies. The three-parameter Sérsic (1963Sérsic ( , 1968) R 1/n model, a generalisation of the two-parameter de Vaucouleurs (1948) R 1/4 has been shown to represent the stellar light profiles of the main body of elliptical galaxies and the bulges of disc galaxies very well over a large luminosity and radial range (e.g. Caon et al. 1993;D'Onofrio et al. 1994;Young & Currie 1994;Andredakis et al. 1995). This model can be defined as where the Sérsic index n (a good proxy for galaxy light concentration; Graham et al. 2001) describes the shape (curvature) of the radial brightness profile and I e denotes the intensity at the half-light radius R e . The variable b n , defined to ensure that the half-light radius encloses half of the total luminosity, is determined by solving numerically Γ(2n) = 2γ(2n, b n ), where Γ(n) and γ(n, x) are the complete and incomplete gamma functions, respectively. For 1 n 10, b n ≈ 2n − 1 3 (Caon et al. 1993). When n = 0.5 the Sérsic model reduces to a two-parameter Gaussian function, whereas for n = 1 it yields an exponential function. We use the Sérsic model to describe the light profile of the bar component in disc galaxies.
The luminosity for a Sérsic model component within any radius R can be determined as follows: where γ(2n, x) is the incomplete gamma function and x = b n (R/R e ) 1/n (Graham & Driver 2005). In general, the radial intensity distribution of the discs of spiral and lenticular galaxies can be well modelled with an exponential function (Freeman 1970), given as where I 0,d and h are the central intensity and scale length of the disc, respectively. We adopt the three-parameter Gaussian radial profile to describe the radial intensity distribution of the outer ring and spiral-arm components of the sample galaxies, given by where I 0,r denotes the highest intensity value of the Gaussian ring with a semi-major axis R r and width σ. Setting R r = 0 yields a two-parameter Gaussian function, which is used in this paper to model additional nuclear light components (i.e. AGN and nuclear star clusters). The bulge component of luminous early-type (core-Sérsic) galaxies have long been known to contain a depleted core (i.e. a luminosity deficit), exhibited as a flattening of the inner stellar light distribution relative to the inward extrapolation of the outer Sérsic profile (e.g. King 1978;Young et al. 1978;Lauer 1985;Trujillo et al. 2004;Ferrarese et al. 2006;Dullo & Graham 2012;Dullo 2019). In contrast, Sérsic (elliptical, S0, spiral, and irregular) galaxies do not have such a luminosity deficit. Throughout this work, we regard bulgeless spiral and irregular galaxies to be Sérsic galaxies. For core-Sérsic galaxies, the radial intensity distributions are well described by the core-Sérsic model which is a blend of an inner power-law and an outer Sérsic model with a transition region . This model is defined as with where I b is the intensity at the core break radius R b , γ is the slope of the inner power-law profile, and α moderates the sharpness of the transition between the inner power-law and the outer Sérsic profile. The half-light radius of the outer Sérsic model represented by R e and the quantity b has the same meaning as in the Sérsic model (Eq. (2)). The total luminosity for the core-Sérsic model (see e.g. Trujillo et al. 2004, their Eq. (A19)) can be written as For R b → 0, this expression becomes the Sérsic expression given by Eq. (3).

Fitting methodology and multi-component decomposition
Constraining how well galaxy optical properties such as mass, size, and stellar concentration correlate with each other, and also with the radio, X-ray, and optical emission line properties ultimately depends on an accurate characterisation of the bulge and disc structural properties. High-resolution near-infrared and optical imaging and IFU spectroscopy have revealed distinct morphological and/or kinematical galaxy structures including bulges, discs, nuclear star clusters, nuclear discs and rings, bars, and spirals (e.g. Faber et al. 1997;Carollo et al. 1998;Graham & Guzmán 2003;Ferrarese et al. 2006;Krajnović et al. 2013Krajnović et al. , 2020Erwin et al. 2015;Catalán-Torrecilla et al. 2017;Dullo & Graham 2014;Dullo et al. 2016Dullo et al. , 2019Bittner et al. 2020;Gadotti et al. 2020;Johnston et al. 2020Johnston et al. , 2021. As noted in the Introduction, if these structures are not modelled separately due to their coupling the derived bulge and disc structures would be biased. Multi-component decompositions are A105, page 7 of 73 Fig. 4. Distribution of inner, sub-kiloparsec stellar structures (point sources, nuclear star clusters, discs and rings, bars, and depleted stellar cores) detected in our photometric decompositions. The size of the sub-kiloparsec structure is defined to be the break radius for galaxies with a depleted core (i.e. core-Sérsic galaxies), whereas for others the size determined by eye from the fits indicates the radius where the subkiloparsec structure contributes significantly to the galaxy light. While a core-Sérsic galaxy can be nucleated, here we only consider its depleted core as a sub-kiloparsec structure. The overwhelming majority (95%) of our sample have sub-kiloparsec stellar structures that are detected in our analyses of the HST imaging with a high spatial resolution ≤0 . 1 (blue vertical line). With high sensitivity e-MERLIN L-band (1.5 GHz) observations, we detected nuclear radio emission, at high resolution of ≤0 . 2, as such close comparison can be made with the analogous optical subkiloparsec structure in the HST data. The red and magenta vertical lines delineate the spatial resolutions for SDSS and Spitzer images. While we show the pixel scale of SDSS, the actual resolution is set by the seeing, that is typically larger than ∼1 . 5 and varying across the sky. Accounting for the effect of the PSF, only 15% (22%) of the kiloparsec-scale structures from HST can be well resolved with Spitzer (SDSS) images. Furthermore, we see no clear trend between the physical size of the inner, sub-kiloparsec structures and total galaxy stellar mass (M * ,glxy ) or T-type. The total galaxy mass (M * ,glxy ) is based on our light profile decompositions ( needed if we are to understand properly the different mechanisms that build up these very distinct photometric components. To this aim, we perform 1D (one-dimensional) and 2D (two-dimensional) multi-component decompositions. This section discusses our 1D and 2D decompositions and the reasons why we regard the results from the 1D (multi-component) light profile modelling to be more appropriate for this study and thus base the analyses throughout this paper on them.

1D fitting procedure
The 1D multi-component fitting procedure utilised to decompose the one-dimensional galaxy light profiles with our personal Fortran programming code was discussed in Dullo et al. (2016Dullo et al. ( , 2017, 2019), Dullo (2019). The fitting code not only accounts for the PSF (see Sect. 3.2.2) but is also sophisticated, capable of fitting simultaneously up to six galaxy components, which are summed up to a full model with (up to) 16 free parameters. The cornerstone of our fitting strategy concerns remedying two commonly existing limitations of galaxy light profile parameterisation in the literature. These are summarised and addressed here.
The first limitation is excluding innermost structures. Spacebased NIR imaging data from Spitzer and most ground-based optical/NIR imaging do not have the necessary resolution to study reliably inner, sub-kiloparsec stellar structures such as NSCs, discs, rings, bars, depleted stellar cores, and non-thermal sources (AGN). Excluding the innermost portions of galaxies' light profiles has been a preferred strategy to overcome such resolution limitations. In Fig. 4 we plot the distribution of subkiloparsec optical/NIR structures detected in our analyses of the HST imaging with a high spatial resolution ≤0 . 1. For core-Sérsic galaxies, the sub-kiloparsec structure is defined to be the depleted core with size R b ( 0.5 kpc, see Table B.2). For Sérsic galaxies, we inspect the results of the decompositions and measure the size of the sub-kiloparsec structure by eye as the radius within which the innermost component contributes significantly to the galaxy light. An overwhelming majority (95%) of our sample have sub-kiloparsec stellar structures detected with HST. The high sensitivity e-MERLIN L band (1.5 GHz) observations allow for analogous radio continuum detections at high resolution of ≤0 . 2. Accounting for the effect of the PSF, we find only 15% (22%) of the sub-kiloparsec, optical/NIR structures from HST can be well detected with Spitzer (SDSS) imaging data.
The second limitation is restricting fits to bulge+disc profiles. It is common, particularly for automated fits, to model elliptical galaxies with a single Sérsic model, and to consider Sersic+exponential bulge-disc fits adequate for disc galaxies. The LeMMINGs sample covers a morphology range that encompasses all Hubble types between E and Im. Dissecting the galaxies into bulges, discs, bars, rings, spiral arms, nuclear sources, cores, and haloes, we derive the fraction of galaxy light that is not ascribed to the bulge and disc components ( f other ) to highlight the need for fits beyond the two main galaxy components, that is bulge-disc profile (see Appendix D). Excluding irregular galaxies, we find that the fraction of LeMMINGs galaxies with f other > 20% is ∼26% (43/167), increasing to ∼40% (67/167) for f other > 5% (Fig. 5). Due to the assignment of f other to the bulge, one plausible outcome of a simple bulge+disc fit to disc galaxies which contain bars, rings, and spiral arms is an unreasonably high B/T ratio. In Fig. 5 we show the dependence of f other on the bulge mass and luminosity. We also plot green and blue tracks which trace the median trend for the binned spiral, S0 and elliptical galaxy data. For spirals, in general, the fraction f other tends to increase as the bulge mass and luminosity increases, whereas S0s typically show f other ∼ 1.3%. For most (71%) elliptical galaxies f other ∼ 0, however at the most massive end dominated by BCGs and cD galaxies the fractional contribution from light components such stellar haloes and ICL increases f other typically to 5%. In summary, the prevalence of galaxy components beyond bulge+disc and the unprecedented level of detail required to decompose our galaxies suggest that Sérsic bulge (exponential disc) fits, especially when forced in an automated fashion are incapable to reproduce our work and, in general, to derive accurate structural parameters of the galaxies' inner regions.
The best-fit model parameters which provide an optimal match to the data are calculated iteratively using the The fraction of galaxy light that is not ascribed to the bulge and disc components ( f other ) as a function of the bulge stellar mass and absolute V-band bulge magnitude corrected for Galactic and internal extinction. S0 and spiral galaxies are shown as green boxes and blue crosses, whereas the elliptical galaxies from this work are denoted by black circles. We also add 13 massive elliptical galaxies (stars) from Dullo & Graham (2014), Dullo (2019). Excluding irregular galaxies, the fraction of LeMMINGs galaxies having f other ≥ 20% is ∼26% (43/167), this figure increases to ∼40% (67/167) for f other ≥ 5% (i.e. above the dashed horizontal line). The blue, green, and black tracks trace the median trend for the binned spiral, S0 and elliptical galaxy data. Most (71%) elliptical galaxies have f other ∼ 0. A track is not shown for our irregular galaxies which are excluded due to their small number (6 objects). Representative error bars are shown at the bottom.
Levenberg-Marquardt minimisation algorithm (see Dullo et al. 2017Dullo et al. , 2019Dullo 2019). We quantify the global goodness of the fits using the root-mean-square (rms) residuals (∆). We let all fitting parameters free to vary. Our identification of a component in a galaxy is with guidance from physical considerations such as the galaxy's morphological type, velocity dispersion, and structural features visible in the optical and NIR images. We visually inspect the high-resolution HST images as well as the deeper ground-based (SDSS/DSS) data (when available) before and after the fitting procedure. We further check for any systematic tendencies in radial , PA, and B 4 profiles 8 ; galaxy components such bars and prominent rings are typically accompanied by local minima/maxima in , B 4 , and a twist in PA For each galaxy, the fits were carefully examined to determine the one that provides the best description. Fits that we regard optimal are those with: (i) a minimum number of galaxy components, (ii) low rms residual values and (iii) physical meaningful best-fitting model parameters that describe the full complexity of the galaxies' stellar distributions well over the entire radial range. Ideally, acceptable fits do not exhibit systematic ('snake-like') residual patterns which commonly signify a mismatch between the data and fitted model, admittedly however a few sample galaxies are subject to modest residual patterns due to the presence of dust, noise in the data or simply because of difficulties in modelling the galaxy light profiles. In Sect. 3.2.3, we measure the best-fitting model component parameters, implementing an MC approach, which are found to agree very well with the adopted parametrisation from the multi-component decompositions discussed here.
In addition to the rms residuals, each galaxy's goodness of fit is given a 'quality flag' designated by numbers ('1' or '2') and listed in Tables B.2-B.4 along with the best-fitting model 8 A caveat is the presence of strong dust structure in a galaxy can also cause patterns in the in and B 4 profiles. parameters form the decompositions. Quality flag '1' indicates fits with good or higher quality levels that meet all our aforementioned criteria for an optimal fit, whereas those labelled '2' are deemed questionable. The latter are mainly associated with highly inclined (i 75 • ) and morphologically complex faint galaxies (see Sect. 2.4) whose surface brightness profiles could not be accurately extracted and modelled, and in few cases the galaxy centre was affected by dust that results in unreliable decomposition (e.g. NGC 3665). Those with quality flag '2' make up ∼13% of the sample.
As an illustration of our fitting procedure, in Fig. 6 (left) we show a core-Sérsic light profile and a six-Sérsic-component light profile: a model fit to the light profile of the giant elliptical NGC 3348 and a six-component decomposition of a doubledbarred lenticular LeMMINGs galaxy NGC 2859 into a nucleus, a bulge, an inner bar, plus three outer galaxy components (bar+ring+disc). The bulge and bars are each described with a Sérsic model, while the disc is modelled with an exponential function. We fit a two-parameter Gaussian function to the nuclear source tentatively identified here as AGN (see Baldi et al. 2018Baldi et al. , 2021a, whereas the outer disc was well reproduced by a three-parameter Gaussian ring model. Each galaxy model was convolved with a Gaussian PSF during every fitting iteration. Inspection of the HST and SDSS images (Fig. 6, right) strongly corroborates the results of the decompositions.
Figures D.1-D.3 display the decompositions of the majoraxis surface brightness profiles of the galaxies listed in Tables B.1−B.4 and the residual profiles from the fits. The rootmean-square (rms) residuals ∆ are also shown. The median value of ∆ for our fits is ∼0.065 mag arcsec −2 .

PSF treatment
When decomposing a galaxy light profile, at each fitting iteration, the individual model components were convolved with the The stellar light distribution of the nuclear component (which we tentatively identify as AGN, see Baldi et al. 2018Baldi et al. , 2021a is modelled with a two-parameter Gaussian function and dominated by the Sérsic (n ∼ 3.6) bulge. The inner and outer bars are described using two Sérsic models each with n ∼ 0.2. The outer, extended exponential disc is the dotted blue curve. The three-parameter Gaussian ring model (cyan dot-dashed curve) describes the outer ring. The two bars and the large-scale ring are accompanied by three local maxima in and B 4 and twists in PA. Each model component is convolved with a Gaussian PSF. The residual profile along with the small rms residual ∆ ∼ 0.047 mag arcsec −2 . We fit six model components which are summed up to a full model with 16 free parameters. Right-hand panel: SDSS image of NGC 2859. The top and bottom insets show the surface brightness contours of the galaxy's SDSS and HST ACS images, respectively. North is up and east is to the left. point-spread function (PSF) in 2D. The PSF implementation is discussed Dullo et al. (2017), the mathematical expressions are presented in Trujillo et al. (2001a, their Eq. (2)). For more reference see Pritchet & Kline (1981), Saglia et al. (1993). The PSF of a telescope + instrument, which scatters the stellar light from the innermost concentrated regions of the imaged galaxy to more extended regions, gives a measure of the image's spatial resolution. Over the years, several analytic functions have been used to describe the PSF. Trujillo et al. (2001b) remarked that the Moffat function is numerically more well behaved in modelling HST PSFs than a Gaussian function. While the design consideration of our multi-component decomposition code is to perform PSF convolution using either a Gaussian or Moffat function, we find that modelling the HST PSF using the latter function is computationally intensive, a striking contrast to the markedly faster convolution with the former one. Motivated by this difference we devoted efforts to investigate the effect of our treatment of the PSF on the best-fitting structural parameters.
As our intention is to illuminate severe effects of the PSF, we select a dozen LeMMINGs galaxies with a steep inner light profile slope and perform decompositions with multiple Sérsic models, PSF-convolved using a Moffat and Gaussian function (Appendix D). The FWHMs of the Gaussian PSFs as well as FWHMs and β values for the Moffat PSFs were determined from Gaussian and Moffat fits to the radial flux profiles of several unsaturated stars present in the galaxy HST images using the IRAF task imexamine. In Fig. 7, we compare the Sérsic indices, effective radii, effective surface brightnesses of the bulges and point source magnitudes. The good agreement between the structural parameters from the two PSF convolution approaches sug-gests that no matter which of the two modelling functions are selected to describe the PSF in the HST data the light profile decompositions remain largely unaffected. We therefore adopt the Gaussian function to convolve the HST PSFs throughout this work.

1D fitting analysis and galaxy components
For luminous early-type galaxies with cores (i.e. core-Sérsic galaxies) we show that their bulges can be very well described with the core-Sérsic model, whereas for coreless (i.e. Sérsic) galaxies the underlying bulge light distributions are well fitted by the Sérsic model (Appendix D). We identify 20 (5 S0s and 15 Es) core-Sérsic galaxies out of the total sample of 173 galaxies (∼11.6%) which are defined to have a deficit of light at the centre with respect to the outer Sérsic profiles. HST's resolution implies we can detect depleted cores, as measured by R b , as small as ∼10 pc. However, we caution that the two small core galaxies NGC 5631 and NGC 6482 have R b ∼ 0 . 04 (∼10 pc) and as these cores are from HST imaging with a scale of 0 . 05/px and defined by few innermost data points, their authenticity can be questioned. While we tentatively classify two lenticular galaxies, NGC 1167 and NGC 3665, as core-Sérsic, they behave as intermediate objects, ergo we also decomposed their light profiles fitting multiple Sérsic functions (see Appendix D). For NGC 1167, both the Sérsic and core-Sérsic approaches give good description to the galaxy light profile. Despite the uncertainty which arises from the nuclear dust in the galaxy, the core-Sérsic fit of NGC 3665 is superior in quality than the Sérsic one. For the 5 core-Sérsic S0s (NGC 0507, NGC 1167, NGC 2300, To explore the marked effects of the HST PSFs in the HST data, the galaxies were selected to have a steep inner light profile slope. The solid lines represent one-to-one relations. The good agreement between the best-fitting structural parameters suggests that the light profile decomposition does not depend appreciably on selecting either of the two HST PSF modelling functions employed. NGC 3665, and NGC 5631) the light profiles are modelled as the sum of core-Sérsic bulge and an exponential disc. Of the 15 core-Sérsic elliptical galaxies in the full sample, six (NGC 0315, NGC 0410, NGC 2832, NGC 3193, NGC 5322, and NGC 6482) have outer stellar halo components modelled with an exponential function (Dullo 2019;den Brok et al. 2021). We regard the outer bump (i.e. light excess) over the bulge profile of the core-Sérsic elliptical (E6) galaxy (NGC 3613) to be a disc light that is well described using a Sérsic n ∼ 0.81 profile. Most (80%) core-Sérsic bulges in the sample have n > 4. On the other hand, an overwhelming majority (94%) of Sérsic bulges have n < 4. Furthermore, the bulges of 25 disc galaxies, 80% (68%) of which are late-type (H ii) galaxies, in the sample follow a nearexponential n ∼ 1.12 ± 0.20 light profile and have median values of M * ,bulge ∼ 9.2, ∼ 0.37, B 4 ∼ 0.003, and B/T ∼ 0.17. Given the exponential nature coupled with the median properties these bulges are suspected to be rotationally supported pseudo-bulge candidates. We identify 10 bulgeless galaxies (IC 2574, NGC 0672, NGC 0784, NGC 1156, NGC 2537, which are all latetypes galaxies typically with low stellar masses (M * ,glxy < 10 11 M ). Furthermore, for three low-mass late-type galaxies (NGC 2541, NGC 3031, and NGC 4151) in our sample, the bulges are very compact (R e ∼ 241−268 pc) and with high Sérsic indices (n ∼ 3.9−5.1), which appear to be tell-tale signs of dynamically hot classical bulges akin to those reported in Erwin et al. (2015).
In general we find that the outer disc light distribution is well fitted with an exponential n = 1 model, for 17 late-type galaxies however their outer disc light are better fitted with a nearexponential, n ∼ 0.85±0.50 Sérsic model (see also Dutton 2009;Head et al. 2015). Of the 150 disc galaxies (S, S0, and Irr) in the sample, our structural analyses have identified bar structures in 47 (30 S, 15 S0, and 2 Irr). The light profiles of bars are modelled with a low-n Sérsic profile. Detailed analysis of bar structures and the dependence on host galaxy properties is deferred to a future work.
With spiral arms and rings in a galaxy disc typically resulting in excess flux ('a bump') at large radii relative to the exponential disc profile, they are generally well modelled using the threeparameter Gaussian ring function. In Appendix A, we present notes on our decompositions and literature comparison for 54 selected individual galaxies.

2D fitting
Past studies (e.g. Ciambur 2016; Savorgnan & Graham 2016) highlighted the advantage of 2D decompositions of the full galaxy light distributions over 1D decompositions of brightness profiles obtained from major-axis cuts (see Sect. 3.2.5). This is particularly true for spiral galaxies which are known to have multiple components such as nuclei, bars, rings, and spiral arms. In order to test the robustness of our 1D modelling, we perform 2D decompositions of the HST images of 65 sample galaxies using imfit 9 v1.8 (Erwin 2015). The 65 galaxies with 2D decompositions, make up ∼40% of the 163 newly analysed galaxies in the paper. We note that our 1D brightness profiles for the sample galaxies are not derived from major axis cuts (Sect. 2.4). Nonetheless, in our 2D analysis we focus on spiral galaxies and perform 2D decompositions for nearly half of the spiral galaxies in the paper (i.e. 49/102). These 49 spiral galaxies are representative of the full spiral galaxy sample in terms of the morphological T-type, galaxy mass, and AGN power as measured by the O iii λ5007 Å line luminosity; we find that the KS tests on the datasets from the sub-sample and the full sample give significance levels P ∼ 0.52, 0.39, and 0.19 for the three quantities, respectively.
Our 2D fitting contains the same type and number of galaxy structural components as the corresponding 1D fitting. The only exception is NGC 2859: while we fit a six component model to the galaxy's 1D, major-axis light profile (see Fig. 6), the 2D decomposition of the HST image reveal a seventh component, that is a faint inner ring which encloses the inner bar. In order to convolve the 2D model images, we created a Moffat PSF using the imfit task makeimage. The FWHM and β values of the PSF were measured using stars in the HST images of the galaxies. The 2D mask images were generated by converting the '.pl' mask files from the ellipse run into images using the iraf task wfits. Appendix E presents the best-fitting model parameters from the 2D decompositions, whereas Appendix A provides a discussion by comparing our 2D decompositions with those in literature for 23 galaxies. Appendix E shows the residual images for a representative sample of 28/65 galaxies created after subtracting the imfit model images from the HST images of the galaxies.  Table B.1).

Comparison of 1D and 2D decompositions
We now compare our 1D and 2D decompositions and explain why the results and conclusions of this work are based on the 1D decompositions. Figure 8 plots the 1D and 2D bulge properties including the Sérsic indices, n (a), effective radii, R e (b), surface brightnesses at R e , µ e , (c) and absolute V-band magnitudes, M V , (d) for the 65 sample galaxies. We note that the 2D R e values were converted into major-axis values using the corresponding bulge ellipticities (Fig. 8b). The strong agreement between the 1D and 2D measurements of n, R e , µ e , and M V is such that 85, 88, 92, and 86% of the 1D versus 2D dataset pertaining to n, R e , µ e , and M V , respectively, are within 1σ (see Table B.1 and Sect. 3.4) of perfect agreement. The 1D and 2D n, R e , µ e , and M V values are within the 2σ error ranges for 91, 95, 97, and 97% of the cases. Similarly, Fig. 9 reveals good agreement between the 1D and 2D bulge-to-total ratios (B/T ). Overall, we find consistency between the 1D and 2D best-fitting parameters when obtained carefully through multi-component decompositions. Mild discrepancies that we observe between the 1D and 2D measurements are expected to arise from intrinsic differences between the two methods (see Figs. 8 and 9). An advantage of performing a 2D decomposition over a 1D decomposition is that the former makes use of all available information in a 2D galaxy image to determine the best fit model image. As noted above, the 2D method of galaxy decomposition is therefore considered to provide a better description to non-axisymmetric stellar structures, which are common in spiral galaxies. However, we note that non-axisymmetric structures are also captured by azimuthally averaged 1D light profiles (e.g. Ciambur 2016). 2D decompositions have limitations; each 2D model component, generated using for example galfit (Peng et al. 2010) or imfit (Erwin 2015), is represented by a single ellipticity, PA, and boxy/discyness parameter constant over all scales. In contrast, distinct galaxy components can have radial gradients in ellipticity and PA, as is the case for triaxial systems which are seen in projection (e.g. Mihalas & Binney 1981). Determining whether residual structures from a 2D fitting are caused by radial variations in ellipticity and PA or are real galaxy features that were unaccounted for in the 2D modelling can therefore be difficult (e.g. Peng et al. 2002). This poses a problem when determining the appropriate number of 2D model components that are needed to fit a galaxy image.
However, the radial tendencies of ellipticity, PA, and boxy/discyness profiles are well represented in 1D brightness profiles determined from the azimuthal averaging of the data in the galaxy image. As such 1D residual profiles are usually infor- mative about distinct galaxy components, as long as they are not based on major-axis cut fits. This in part has motivated several authors to advocate the 1D method (Gilhuly & Courteau 2018; see also Ciambur 2016;. These 1D methods can also be more beneficial than the 2D counterparts for sky background determination since 1D brightness profiles can be easily constructed by combining high-resolution space-based data and deep ground-based data. To summarise, the strong agreement between our 1D and 2D decompositions, the capabilities of our 1D decompositions to capture radial tendencies of ellipticity, PA, and boxy/discyness profiles of galaxies and the informative nature of the 1D residual profiles are the reasons why we adopte the 1D modelling throughout this paper.

Ubiquity of nucleation
The bulk of the sample galaxies (124/173) have central light excesses (sub-kiloparsec nuclear optical components), with respect to the inward extrapolation of the Sérsic or core-Sérsic model fitted to the bulge, detected in our analyses of the high-(spatial resolution) HST imaging. This yields a high nucleation frequency of ∼72%, but we note that nuclei are less common (∼10−20%) in core-Sérsic galaxies, which are among the most massive galaxies (see also Côté et al. 2006;Dullo & Graham 2012; A105, page 12 of 73

Magnitudes, M/L, and stellar masses
Apparent magnitudes of the galaxy components from the 1D fitting are determined using the best-fitting major-axis, Sérsic and core-Sérsic structural parameters along with the ellipticity values, which are integrated to R = ∞ (e.g. Dullo & Graham 2014;Dullo 2019). For each fitted galaxy component, we use our ellipticity profiles (Appendix D) to measure an average ellipticity value over the region where the component contributes significant light to the galaxy light distribution. For the 2D galaxy components, the apparent magnitudes are measured by integrating the best-fitting structural parameters. Foreground Galactic extinction corrections of the magnitudes are based on the reddening values from Schlafly & Finkbeiner (2011). For disc (S0, spiral, and irregular) galaxies, we additionally correct for the inclination (i) dependent, internal dust attenuation using Driver et al. (2008, their Eqs. (1) and (2) and Table 1). These corrections differ between the bulge and disc components and can be expressed as: where m corr and m obs denote the observed and corrected magnitudes, respectively, and the values of the coefficients b 1 , b 2 , b 3 , d 1 , d 2 , and d 3 vary with the filter type as tabulated in Driver et al. (2008 , Table 1). To compute cos(i) = b/a, the minor and major diameters of the galaxies were obtained from NED. The same amount of internal dust correction applied to the bulge was applied to the nuclear and intermediate galaxy components ,   Table 4. Transformation into V-band magnitude.
that is including nuclear discs. Analogously, the disc and other outer galaxy components were treated as the same for the internal dust correction. With component apparent magnitudes corrected for dust, they are then converted into stellar masses in a three-step procedure. First, we computed luminosities in solar units after transferring the apparent magnitudes into absolute magnitudes. To achieve this the galaxy distances from NED (3K CMB) and other sources (Karachentsev et al. 2004;Bonanos et al. 2006) and absolute magnitude for the Sun from Willmer (2018) were used. Next, we calculated the stellar mass-to-light ratios (M/L) for the galaxies using the waveband dependent (B − V colour)-(M/L) λ relations from Zibetti et al. (2009, their Table B1) appropriate for elliptical galaxies and those from Into & Portinari (2013, their Table 4) for the disc galaxies assuming a Salpeter stellar initial mass function IMF 10 (Salpeter 1955). The total B − V colours are from HyperLEDA and for each sample galaxy all the fitted components are assumed to have the same M/L. Finally, the luminosities coupled with the estimates of stellar mass-to-light ratios in the corresponding bands yield the stellar masses. To calculate the total stellar mass of a galaxy we add the stellar masses from all its components (see Table B.5).
In order make a direct comparison with previous work and among the sample galaxies, as the data are measured in various HST photometric bands (see Table B.2), the component absolute magnitudes are transformed into V-band. To do so, the galaxies' total V magnitudes and B − V colours were taken from HyperLEDA. The bulk (61%) of the sample galaxies have their magnitudes measured in HST F814W (∼I-band), which we calibrated into V-band by comparing their total V − F814W and B − V colours. Figure 10 shows the correlation between the two colours and the resulting ordinary least squares (OLS) bisector fit, which is the adopted transformation equation (see Table 4). Applying the same prescription, transformation equations were also derived for the HST F702W and F850LP filters. For the remaining HST filters when necessary the calibration was performed using Fukugita et al. (1995, their Table 3) and Lauer et al. (2007b), see Table 4. The absolute magnitudes, stellar masses, M/L for the sample galaxies are reported in Table B.5. Figure 11a plots the total stellar mass of a galaxy (M * ,glxy ) as a function of its bulge stellar mass (M * ,bulge ). We also show the correlations between the absolute V-band galaxy magnitude (M V,glxy ) and absolute V-band bulge magnitude (M V,bulge ) for the sample galaxies (Fig. 11b). We find strong correlations between M * ,glxy and M * ,bulge for our early-type, late-type and full sample of galaxies (Pearson correlation coefficients r ∼ 0.77−0.94, see Table 5 Table B.5). Bulgeless galaxies have been excluded from our analyses. The OLS bisector fits to the early-type data (solid line), late-type data (dashed line) and full galaxy sample (dotted line).
Notes. Galaxy stellar mass (M * ,glxy ) and absolute V-band galaxy magnitude (M V,glxy ) as a function of bulge stellar mass (M * ,bulge ) and absolute V-band galaxy magnitude (M V,bulge ) for early-type galaxies (ET), late-type galaxies (LT) and combination of both. The different columns represent the OLS bisector fits to the data and the Pearson correlation coefficients (r p ).
(r ∼ 0.72−0.94, Table 5). After excluding bulgeless galaxies and applying the ordinary least squares (OLS) bisector regressions (Feigelson & Babu 1992), our fits to the (M * ,glxy , M * ,bulge ) and (M V,glxy , M V,bulge ) data are listed in Table 5. For late-type galaxies, we find that the scatter in the M * ,glxy − M * ,bulge and M V,glxy − M V,bulge relations increases as the galaxy/bulge stellar mass and luminosity decrease. This suggests that SMBH scaling relations which are based on M * ,glxy and L glxy do not have the same predictive power as those constructed using M * ,bulge and L bulge for a galaxy sample containing late-type galaxies (e.g. . Given that the masses of SMBHs are known to correlate well with the properties of the bulges (e.g. Kormendy & Ho 2013; Graham 2016), Fig. 11 underscores the importance of separating the bulge component from the rest of the galaxy through decomposition, especially for late-type galaxies.

Error analysis
We determine realistic errors for the 1D best-fitting structural parameters, the bulge magnitudes and stellar masses after performing multi-component decomposition of simulated galaxy light profiles. We generated over 100 realisations of the surface brightness profile for each galaxy by running a series of MC sim- ulations. This prescription, discussed in detail in Appendices F, was also used to test the robustness of the multi-component decomposition of the galaxy light profiles. As an illustration of this exercise of measurements of the errors, in Fig. 12 we show a corner plot for the bulge of the LeMMINGs galaxy NGC 959 fitting each of the 300 realisations of the galaxy with a four-component (bulge+disc+spiral-arm+nucleus) model, as done for the original galaxy light profile (see Appendix D).

Bulge structural properties and connection with emission-line class and radio core emission
In this section we study sets of bulge structural properties over large stellar mass and morphology ranges obtained from detailed modelling of HST surface brightness profiles for a sample of 173 LeMMINGs galaxies. We also examine how the bulge properties vary as a function of host optical emission class and radio morphological structures from e-MERLIN (Baldi et al. 2018(Baldi et al. , 2021a. In the analysis of the galaxies properties, we divided the sample into two optical morphological classes, early-type galaxies (Es and S0s) and late-type galaxies (Ss and Irrs), Hubble 1926;Sandage 1961. While the study of the dependency of bulge structural properties on galaxy properties such as stellar mass and morphology is not new (e.g. Kormendy & Bender 2012; Conselice 2014), robust characterisation of galaxy structures using homogeneously measured high-resolution optical and radio data for a large sample of galaxies has not been possible to date. In particular, our analysis focuses on the bulge component as it is known to correlate better with the mass of the SMBH and other central galaxy properties than the galaxy disc at large radius.

Inner light profile slope of the bulge
Inner slopes of the surface brightness profiles of early-type galaxies have been extensively studied as a diagnostic of galaxy formation mechanisms. The advent of high-resolution data afforded by HST prompted the detection of a bimodal distribution of the negative, logarithmic, inner profile slopes (γ) for early-type galaxies (e.g. Ferrarese et al. 1994;Lauer et al. 1995;Gebhardt et al. 1996;Byun et al. 1996). With their bulges preferentially pressure-supported and built amid major dry mergers, bright galaxies with depleted cores display shallow inner profiles with power-law indices γ 0.3 (e.g. Ferrarese et al. 2006;Dullo & Graham 2012;Dullo 2019). In contrast, early studies have shown that low-and intermediateluminosity early-type bulges-which are likely a consequence of gas-rich processes-are consistent with being rotationally supported and possess profiles that are steep (γ 0.5) all the way into the centre. The evidence for the bimodal distribution of γ however remains a matter of hot debate. Some studies identified a few galaxies having intermediate slopes 0.3 γ 0.5 and cautioned that the bimodality is weak (Rest et al. 2001;Ravindranath et al. 2001) while others later argue the bimodality is entirely nonexistent after finding more objects with intermediate slopes. Another key unknown is the distribution of γ for late-type bulges and how it compares with those of early-type bulges. As discussed extensively in Dullo & Graham (2012), see also Sect. 3.1, our core identification is not based on γ but instead a 'core-Sérsic' galaxy is defined to have a central stellar light deficit relative to the inward extrapolation of the bulge's outer Sérsic profile. In contrast, a Sérsic galaxy has no such central light deficit (Sects. 3.1 and 3.2). Nonetheless it is beneficial to explore any systematic trend of γ with bulge stellar mass for our Sérsic and core-Sérsic galaxies.
In Fig. 13 we plot the negative, logarithmic, inner profile slope of the bulge, unaffected by PSF convolution and not biased by the nuclear light component, against its stellar mass (M * ,bulge ) for 211 (66 core-Sérsic plus 145 Sérsic) galaxies. Of the 211 galaxies, 173 are from this work, whereas the remaining 38 are non-LeMMINGs core-Sérsic galaxies from Dullo & Graham (2014), Dullo (2019), Dullo et al. (2021). For a core-Sérsic bulge, γ is the slope of the inner power-law profile of the fitted core-Sérsic model given by Eq. (6), while for a Sérsic bulge γ determined from the fitted Sérsic model at any radius R can be written as follows: We use local slopes measured at a representative HST resolution limit R = 0 . 1, that is γ Ser (R = 0 . 1) for Sérsic bulges. As Ferrarese et al. (2006), Côté et al. (2007), Glass et al. (2011), we find no evidence for a bimodal γ distribution for early-type galaxies in the γ − M * ,bulge plane, instead the distribution appears continuous. All core-Sérsic galaxies have M * ,bulge 8 × 10 10 M and low values of γ 0.35, except for NGC 584 with γ ∼ 0.47, whereas Sérsic early-type galaxies (M * ,bulge ∼ 4.0 × 10 8 −3.6 × 10 10 M ) are spread over a large γ range covering from 0.01 to 0.70 and only 11/46 (24%) of them have γ > 0.5. For early-type galaxies, a Sérsic to core-Sérsic regime transition occurs across a mass range M * ,bulge ∼ 6 × 10 10 −3 × 10 11 M where 13 early-type Sérsic galaxies (10/13 with 0.3 γ 0.5) and 27 core-Sérsic galaxies coexist. For the first time, we examine the inner logarithmic slope distribution for the bulges of late-type galaxies as a function of stellar mass. Akin to early-type Sérsic galaxies, late-type Sérsic galaxies exhibit a wide range of γ, from shallow slope values (γ ∼ 0) to very steep ones (γ ∼ 1.71). For late-types, we also find a mild tendency for the slopes to be systematically steepened as the bulge masses decreased.
A caveat for our inner slope analysis of Sérsic galaxies is the effect of distance. The local Sérsic slope γ(R = 0 . 1) is calculated at an angular radius of R = 0 . 1 instead of using a fixed physical radius, thus identical Sérsic galaxies at different distances will have different γ(R = 0 . 1) values. Concerned by this, we measure local slopes for the Sérsic bulges at constant R/R e values of 0.002 and 0.005. Figure 14 shows local slopes γ(R/R e = 0.002), γ(R/R e = 0.005) and γ(R = 0 . 1) as a function of bulge mass for our Sérsic galaxies. The plot reveals that the γ distributions for Sérsic galaxies obtained from the three methods are largely consistent. That is there is no strong impact of the galaxy distance on the γ(R = 0 . 1) distribution discussed above. This is expected given the bulk (60%) of our galaxies have distances ∼10−30 Mpc. To summarise, our results demonstrate that γ alone cannot uniquely separate core-Sérsic and Sérsic galaxies.   Fig. 15. Hubble type as a function of (a) the size of the bulge as measured by its effective (half-light) radius R e , (b) the bulge-to-total light ratio (B/T) and (c) the average ellipticity of the galaxy within R e calculated excluding the most-PSF affected data points inside R ∼ 0 . 1−0 . 2. The mean trends together with the standard deviation within each bin 1σ are shown. Going across the Hubble sequence from the irregular types to elliptical galaxies, the bulge becomes larger, dominant, and round.
for the sample of 173 galaxies, consisting of 23 Es, 42 S0s, 102 Ss, and 6 Irrs. When calculating the average galaxy ellipticities, we excluded the most PSF-affected data. The mean trends for the different correlations together with the standard deviation within each bin are shown. A clear trend emerges as one moves across the Hubble sequence from the irregular types to elliptical galaxies, the bulges are generally larger, more prominent, and round. Our results confirm previous work ( based on detailed decompositions of high-resolution HST data for a large sample of galaxies. For spiral, S0, and elliptical galaxies, we find median values of B/T ∼ 0.21, 0.46, and 0.86, respectively. While most Sab galaxies in the sample have flatter ellipticities than S0s and Es, four (IC 520, NGC 278, NGC 3344, and NGC 7217) have roughly the same round ellipticities as S0s. Most massive elliptical galaxies posses a single spheroidal component, therefore their B/T ∼ 1. All the 10 bulgeless (i.e. B/T = 0) galaxies in our sample, which are not shown in Fig. 15, have morphological type later than Sbc (i.e. Hubble stage T ∼ 5.2−9.8).

The relations between optical emission-line class and bulge properties
In Fig. 16, we plot the relation between the bulge Sérsic index (n) and nuclear emission-type colour-coded on the basis of the bulge stellar mass (M * ,bulge ). Seyferts, ALGs, LINERs, and H ii galaxies constitute 5.8%, 13.3%, 41.0%, 39.9% (6.4%, 10.0%, 33.6%, and 50.0%) of our (the full LeMMINGs) sample, respectively (Table 3). Regarding morphological type nearly half (52.0%) of the elliptical galaxies in the sample are ALGs and the remaining are LINERs (48.0%). We find that 65.0% of elliptical galaxies have bulge mass M * ,bulge 10 11 M and that all those with M * ,bulge < 10 11 M are classified as ALGs. For S0s the emission-line type breakdown is 7.1% Seyfert, 57.2% LINER, 26.2% ALG, and 9.5% H ii galaxies. An overwhelming majority (94.0%) of H ii galaxies are spiral galaxies, the remaining 6.0% are S0s. Furthermore, all irregular type galaxies in the sample are H ii galaxies.
We find that n is closely correlated with M * ,bulge , which is expected as high-n galaxies are shown to be brighter (see Sect. 3.2). Galaxies with large values of n are massive, except for the low mass tiny classical bulges having large values of n (see also Sect. 3.2.3). We find that most (∼74.0%) disc galaxies have a bulge Sérsic index n 3. LINERs, ALGs, and H ii span a broad range in n, unlike Seyferts which have a narrow n distribution: n ALG ∼ 1.3−6.9, n LINER ∼ 0.  (Table B.5) on the Sérsic index (n) and stellar mass of the bulge (M * ,bulge ), Table B.2. Filled circles show the galaxies in our sample that are radio detected with e-MERLIN at 1.5 GHz, whereas open circles are for undetected galaxies. Core-Sérsic galaxies are enclosed in blue boxes and 7/20 (35%) of them are undetected with e-MERLIN at 1.5 GHz. For comparison, the radio detection rate with e-MERLIN at 1.5 GHz for the full LeMMINGs sample of 280 galaxies is 44%. Radio-loud galaxies are enclosed in orange diamonds. A representative error bar is shown at the bottom. and n Seyfert ∼ 1.8−3.7. As can be seen, LINERs and ALGs exhibit similar distributions in the n − M * ,bulge diagram, they are among the most massive galaxies with high n values in the sample. Conversely, H ii galaxies are almost exclusively associated with low-mass galaxies typically with n 2.5. At fixed n, LINERs, ALGs, and Seyferts have higher bulge masses than H ii galaxies. Of the 69 H ii galaxies in the sample 61 (88%) have While not shown in Fig. 16, all the 10 bulgeless galaxies in our sample are H ii galaxies. This points to the low incidence of AGN for bulgeless galaxies, some of which are known to house a SMBH but rarely an AGN (see Kauffmann et al. 2003;Greene et al. 2008;Simmons et al. 2013). At M * ,bulge 2 × 10 11 M , core-Sérsic galaxies (enclosed in boxes, Fig. 16) make up 61% of the galaxy population, although a core-Sérsic bulge can possess a stellar mass as low as M * ,bulge ∼ 8 × 10 10 M (Fig. 16) 4.2.3. The relations between radio detection and bulge properties Figure 17 shows the n − M * ,bulge relation colour-coded by galaxy radio detection. Of the 173 sample galaxies, we detect radio emission 0.2 mJy from 83 with e-MERLIN at 1.5 GHz (filled circles). This gives a detection rate of 48%, in fair agreement with that from the full sample (44%). Of the 83 radio-detected galaxies, five (NGC 3034, NGC 3838, NGC 4242, NGC 5273, and NGC 5907) are radio-detected but 'core-unidentified 11 ' (see Table B.5); they have low Sérsic indices (n ≤ 3), faint bulge magnitudes (M V,bulge −18.9 mag), and low stellar masses (M * ,bulge 10 10 M ). For core-Sérsic galaxies the detection rate is 65%. The remaining 90 undetected sources are denoted by open circles. The radio detection fraction increases with bulge mass, n, and B/T . At M * ,bulge 10 11 M (n ∼ 4.00 ± 1.50, B/T ∼ 0.78 ± 0.24), this fraction is 77%, but it plummets to 24% for M * ,bulge < 10 10 M (n ∼ 1.55 ± 1.06, B/T ∼ 0.10 ± 0.10). Large M * ,bulge and high n values, however, are not strictly associated with radio detection at e-MERLIN sensitivity. Most massive undetected sources are ALGs, the second most common being LINERs. Baldi et al. (2021b, see Sect. 3) used L core /L [OIII] together with SMBH masses for discriminating 'radio-loud' galaxies from 'radio-quiet' ones in the LeMMINGs sample. They classified as radio-loud those galaxies having log (L core /L [OIII] ) > −2 and log (M BH /M ) > 7.7. Using the Baldi et al. (2021b) classification, we find 18 radio loud galaxies out of the sample of 173 galaxies, 17 of which are LIN-ERs and one is a (jetted) H ii galaxy (see Table B.5).

Correlations between radio luminosity and bulge size, B/T, and ellipticity
Here we compare the global optical properties of our galaxies with the radio core luminosity. Figure 18 shows a plot of L Radio,core as a function of the bulge effective (half-light) radius (R e ), dust corrected bulge-to-total light ratio (B/T ) and average galaxy ellipticity value within R e (i.e. ). L Radio,core very nicely correlates with R e (Pearson coefficients r p ∼ 0.53 and Pearson probability P ∼ 6 × 10 −14 ) and B/T (r p ∼ 0.33, P ∼ 8 × 10 −6 ). L Radio,core also correlates reasonably well with (r p ∼ −0.24, P ∼ 1 × 10 −3 ), but the scatter in this relation is large  Fig. 18. The radio core luminosity (L Radio,core ) against (a) the size of the bulge as measured by its effective (half-light) radius R e , (b) the bulge-tototal light ratio (B/T ) and (c) the average ellipticity within R e calculated excluding the most-PSF affected data points inside R ∼ 0 . 1−0 . 2. The dotted curves trace the mean trend. Apparent is the tendency for galaxies with brighter radio core luminosities to be bulge prominent, round, and with large bulge sizes. A representative error bar is shown at the bottom of each panel.
(see Fig. 18c). The mean trends observed here are reminiscent of those in Fig. 15: moving towards brighter radio core luminosities the bulge becomes progressively dominant, round, and larger. Although not shown in Fig. 18, we find that the galaxies' radio core luminosities correlate well with their Hubble types (r p ∼ −0.46, P ∼ 3 × 10 −10 , see Fig. 19).
By separating the sample into Sérsic and core-Sérsic galaxies, we can identify any trends between the radio core luminosity (L R,core ) and the median of the isophote shape parameter of the bulge inside R e (B 4 ) and average bulge ellipticity within R e (Fig. 19). For B 4 , we opt for the median as the mean B 4 values can easily be biased negative or positive (away from the central tendency) by a few extremely high and low B 4 values in the dataset. We witness a trend for core-Sérsic galaxies (enclosed in boxes in the figure) to be systematically round. They also show a tendency to possess high radio core luminosities and boxy-distorted or pure elliptical (i.e. neutral) isophotes, in broad agreement with past work (e.g. Hummel et al. 1983;Faber et al. 1997;Disney et al. 1984;Wrobel & Heeschen 1991;Sandage 1965;Sadler et al. 1995;Capetti & Balmaverde 2005;Richings et al. 2011). We consider |B 4 | values < 0.001 as neutral. In general, Sérsic galaxies tend to have discy isophotes, high ellipticities, and lower radio core luminosities. However, we did not find evidence for the previously alleged strong tendency of the central structures to correlate with radio core luminosity, radio-loudness, isophote shape, and bulge ellipticity (e.g. Lauer 2012;Kormendy et al. 2009;Kormendy & Ho 2013) as revealed by the fact that roughly half (11/20) of the core-Sérsic galaxies show boxy/neutral isophotes, while the remaining half exhibit discy isophotes.
Furthermore, all core-Sérsic galaxies fall inside the region defined by L R,core 10 35 erg s −1 and 0.3 but they are cospatial with early-and late-type Sérsic galaxies which cover a wide range of L R,core , B 4 , and . For core-Sérsic, Sérsic, earlytype and late-type galaxies in our sample, we measure median B 4 values ∼(−0.007 ± 0.08) × 10 −2 , (0.202 ± 0.90) × 10 −2 , (0.174 ± 0.80)×10 −2 , and (0.162 ± 0.90) × 10 −2 and average values ∼ 0.21 ± 0.08, 0.30 ± 0.15, 0.24 ± 0.10, and 0.32 ± 0.16. For the sample disc galaxies, we find that the bulge ellipticity differs (by more than 15%) from the average galaxy ellipticity exterior to R e in 77% of the cases, indicating the ellipticities in these systems are not due to projection. As for the radio-loudness of our core-Sérsic galaxies, we show that only 7/20 (35%) are radio-loud. However, as in the case of core-Sérsic galaxies the hosts of most radio-loud sources are massive bulges with low ellipticities and boxy-distorted isophotes. The fraction of radioloud sources with M * ,bulge > 10 11 M , B 4 < 0, and < 0.32 are 72%, 61%, and 67%, respectively. We determine the typical uncertainties associated with and B 4 to be ∼20% and ∼35%. In summary, Sérsic and core-Sérsic bulges cannot be distinctively distinguished by their radio-loudness, L R,core , B 4 , or , in clear departure from past conclusions.
In their IFU stellar kinematic study of the ATLAS 3D galaxy sample of 260 early-type galaxies, Emsellem et al. (2007Emsellem et al. ( , 2011, Krajnović et al. (2013) revealed that most slow rotators (SR) which are relatively round systems with <0.3 are core-Sérsic galaxies, further reinforcing the dry merger scenario, while most fast rotators (FR) which span a large range in ∼ 0.05−0.6 are coreless galaxies (see also Naab et al. 2014). Krajnović et al. (2013) also went on to remark that the correspondence between these galaxies' central structure and kinematics is not one-toone. There are 30 early-type galaxies in common between our sample and the ATLAS 3D galaxy sample (Emsellem et al. 2011) which encompass 4 SRs and 26 FRs. We confirm the tendency for SRs (FRs) to be core-Sérsic (Sérsic) galaxies: 3/4 SRs are core-Sérsic galaxies, while 22/26 FRs are Sérsic galaxies. While the radio detection fraction that we find for FRs (12/26) Fig. 19. The radio core luminosity against the median of the isophote shape parameter of the bulge inside R e (B 4 ) and average bulge ellipticity within R e both derived after excluding the most-PSF affected data points inside R ∼ 0 . 1−0 . 2. Downward arrows denote upper limits for undetected sources. Boxy and discy isophotes have B 4 < 0 and B 4 > 0, respectively. Core-Sérsic galaxies (enclosed in boxes) show a tendency to possess brighter radio core luminosities, they are also systematically rounder (see Sect. 4.4). However, only 7/20 (35%) of our core-Sérsic galaxies are radio-loud. Radio-loud galaxies are enclosed in orange diamonds. In the left panel, the red-dashed line separates boxy and discy galaxies, whereas in the right panel such lines are shown to demarcate the location of core-Sérsic galaxies. Nuclear radio emission is more prevalent in round, boxy sources than in flat, discy ones. A representative error bar is shown at the bottom of each panel.
different from that of the SRs (1/4), conclusions, about the connection between the bulge kinematics and nuclear radio emission, cannot be drawn due to the very low number of SRs in the sub-sample.

Conclusions
We present an accurate structural analysis of high-resolution HST imaging for a representative sub-sample of 173 (23 Es, 43 S0s, 102 Ss, and 6 Irrs) galaxies drawn from the full sample of 280 nearby galaxies in the e-MERLIN legacy survey (LeMMINGs, Beswick et al. 2014). The aim is to investigate the nuclear activity, optical and radio properties at sub-arcsec resolution using HST and 1.5 GHz e-MERLIN radio observations. This work focuses on the results from our HST imaging analyses, which is coupled in more detail with our 1.5 GHz e-MERLIN radio data in Dullo et al. (2023). Using HST (ACS, WFPC2, WFC3, and NICMOS) images, we have extracted new, 1D major-axis surface brightness, B 4 , PA, and ellipticity profiles for a sample 163 LeMMINGs galaxies and these are combined with data for an additional 10 LeMMINGs galaxies from our previous work. We perform accurate multi-component decompositions of the surface brightness profiles, which extend out to R 80−100 and cover 2R e,bulge for 97% of the sample, fitting simultaneously up to six galaxy components (i.e. bulge, disc, depleted core, AGN, NSC, bar, spiral arm, ring, and stellar halo), which are summed up to a full model with (up to) 16 free parameters. The median rms residual for our fits is ∆ ∼ 0.065 mag arcsec −2 . Galaxy components were carefully identified before models were fitted to the galaxy profiles. We also perform 2D decompositions of the HST images 65 sample galaxies, includ-ing nearly half of the spiral galaxies in the paper (i.e. 49/102). We find that, regardless of the galaxy morphology, careful 1D and 2D galaxy decompositions result in strong agreements. This strong agreement, the capabilities of our 1D decompositions (unlike 2D decompositions) to capture radial tendencies of ellipticity, PA, and boxy/discyness profiles of galaxies and the informative nature of 1D residual profiles motivate us to adopt the results from the 1D decompositions throughout the this work. The LeMMINGs HST sample encompasses all morphological types and spans over six orders of magnitude in stellar mass (6 log M * ,bulge 12.5). Our work represents the largest, most detailed structural analysis of nearby galaxies with HST to date, providing accurate structural parameters for 173 galaxies, a major improvement over past studies especially given the large number of 108 late-type galaxies (Tables B.2, B.3, and B.5). Having carefully isolated galaxy components, we derived luminosities and stellar masses for the bulges, nuclear components, and host galaxies (Table B.5). We have also implemented an innovative method to estimate realistic uncertainties on the fit parameters after creating over 100 realisations of each galaxy's light profile that were later decomposed akin to the original galaxy profiles.
The main results are as follows: (1) We have highlighted the need for performing fits beyond the two main galaxy components (i.e. bulge-disc profile), and calculated, from the detailed decompositions, the fraction of galaxy light outside the bulge+disc component ( f other ): f other 20% for 43/167 of sample galaxies and f other > 5% for 67/167. To do this, we excluded irregular galaxies. Fitting a bulgedisc model to light profiles of disc galaxies which have components such as strong bars, rings, and spiral arms could thus A105, page 19 of 73 overestimate the actual bulge/disc mass by 25% for ∼26% of the cases. For spirals, f other tend to increase as the bulge mass and luminosity increases, whereas S0s typically show f other ∼ 1.3%. For most (71%) elliptical galaxies f other ∼ 0, although for the most massive galaxies (such as BCGs and cD galaxies) f other is typically 5%. For nearby galaxies, we suggest using imaging data with a resolution of 0 . 2, particularly in the inner regions, for accurate multi-component decompositions and to reproduce the level of detail achieved in our galaxy decomposition analysis.
(2) An overwhelming majority (95%) of our sample galaxies have sub-kiloparsec stellar structures (including depleted cores) detected with the analysis of the HST imaging. Correspondingly, we find that our high-sensitivity e-MERLIN L-band observations permit for analogous radio continuum detections at a comparable resolution to that of the HST, however, only ∼15% (∼22%) of the sub-kiloparsec structures can be well studied with Spitzer or ground-based imaging data (e.g. SDSS).
(3) Our analysis of the high-(spatial resolution) HST imaging reveals that nuclei, presented as central light excesses, are ubiquitous in nearby galaxies. We find a nucleation frequency of ∼72% (124/173).
(5) We find no evidence for a previously reported bimodal distribution of the inner logarithmic light profile slope γ (e.g. Lauer et al. 2007b) in the γ − M * ,bulge plane due to Sérsic and core-Sérsic galaxies, that could be interpreted as a diagnostic of galaxy formation mechanisms. Instead, the distribution appears continuous, in good agreement with for example Côté et al. (2007). Overall, core-Sérsic bulges have γ 0.35, whereas Sérsic early-type bulges exhibit a large γ range from 0.01 to 0.70. For the first time, we show that late-type Sérsic galaxies exhibit a wide range of slopes (γ ∼ 0.00−1.71), analogous to early-type Sérsic galaxies.
(6) There are strong trends between galaxy Hubble type and bulge size (R e ), bulge-to-total light ratio (B/T ), and ellipticity ( ). Moving from the irregular types to elliptical galaxies, bulges on average are larger, more prominent, and round, confirming past work. Similarly, tight correlations between the radio core luminosity (L R,core ) and R e , B/T , and are such that bulges with brighter L R,core are dominant, round, and larger. Unsurprisingly, all the 10 bulgeless galaxies (5.8%) in our sample have a morphological type later than Sbc. We find median values of B/T ∼ 0.21, 0.46, and 0.86 for spiral, S0, and elliptical galaxies, respectively. The median values of bulge Sérsic indices (n) and B/T for Seyfert, ALG, LINER, and H ii nuclei are n ∼ 2.69, 3.01, 2.84, and 1.48 and B/T ∼ 0.29, 0.57, 0.44, and 0.11.
(7) We find that the fraction of galaxies harbouring emissionline AGN is a strong function of M * ,bulge and M V,bulge . The majority of AGN (83%) hosts have M bulge 10 10 M (M glxy 10 10.5 M ). For low mass galaxies (M * ,glxy 8 × 10 9 M ), we find an AGN fraction of 6.2%. All the 10 bulgeless galaxies in our sample are H ii galaxies, confirming the low incidence of AGN for bulgeless galaxies, some of which are known to house a SMBH but rarely an AGN (see Kauffmann et al. 2003;Greene et al. 2008;Simmons et al. 2013).
(8) We find that the radio detection fraction increases with bulge mass M * ,bulge , n, and B/T . At M * ,bulge 10 11 M , the radio detection fraction is 77%, declining to 24% for M * ,bulge < 10 10 M . Furthermore, we report a tendency for core-Sérsic galaxies to be systematically round and to possess high radiocore luminosities and boxy-distorted/pure elliptical isophotes but there is no evidence for the previously alleged strong correlation of the central structures (i.e. a sharp Sérsic/core-Sérsic dichotomy) with radio-loudness, B 4 , L R,core , and (e.g. Kormendy et al. 2009;Lauer 2012). Of the 20 core-Sérsic galaxies in the sample, only 7/20 (35%) are radio-loud. Also, all core-Sérsic galaxies are confined to the region defined by L R,core 10 35 erg s −1 and 0.3 but they are cospatial with Sérsic galaxies, the latter cover a large range of L R,core , B 4 , and . Nonetheless, our results, are overall in accordance with cosmological models which predict that the most massive early-type galaxies are more round, and have boxy isophotes, compatible with their formation and evolution scenario as a more evolved object which have undergone several major (dry) mergers. This formation process then results in a dominant bulge housing a massive black hole, a crucial precondition for supporting the launch of jets and outflows in the radio band (e.g. Heckman & Best 2014;Baldi et al. 2021b).
We will investigate the relations between the radio core luminosity and the host bulge properties in an upcoming paper. A multi-wavelength view of the origin and formation mechanisms of nuclei (NSCs/AGN) and the AGN triggering processes and their relation with host galaxy environments will be presented in future papers. Further observations and analysis are in progress to exploit synergies from a large sample of multiwavelength (HST optical/near-infrared, e-MERLIN radio and Chandra X-ray) data.
Acknowledgements. We thank the anonymous referee for their careful reading of the manuscript and for their suggestions. We would like to thank Alex Rosenthal for the help with the 2D imfit decompositions of the HST images. The structure and evolution of galaxies and their central regions' with reference PID2019-105602GB-I00/10.13039/501100011033, from the ACIISI, Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant with reference PROID2021010044, and from IAC project P/300724, financed by the Ministry of Science and Innovation, through the State Budget and by the Canary Islands Department of Economy, Knowledge and Employment, through the Regional Budget of the Autonomous Community. JSG thanks the University of Wisconsin-Madison for partial support of this research. CGM acknowledges financial support from Jim and Hiroko Sherwin. Based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA). We would like to acknowledge the support from the e-MERLIN Legacy project 'LeMMINGs', upon which this study is based. e-MERLIN, and formerly, MERLIN, is a National Facility operated by the University of Manchester at Jodrell Bank Observatory on behalf of the STFC. We acknowledge Jodrell Bank Centre for Astrophysics, which is funded by the STFC. This work has made use of numpy (van der Walt et al. 2011), matplotlib (Hunter 2007) and corner (Foreman-Mackey 2016) and astropy, a community-developed core python package for Astronomy (Astropy Collaboration 2013, 2018), and of topcat (i.e. 'Tool for Operations on Catalogues And Tables', Taylor 2005). The data underlying this article are available in the article and in its online supplementary material.

Appendix A: Notes on Selected Individual Galaxies
Direct comparison of our multi-component structural decomposition with past fits in the literature is not straightforward as detailed decompositions of HST data for a large sample of (nearby) late-type galaxies are not available. Here, we comment on our fits for 54/173 galaxies and compare them, primarily, with those in Salo et al. (2015) and Davis et al. (2019) and, in a few cases, with those from other studies in the literature. Salo et al. (2015) provided two-dimensional, multi-component decompositions of 3.6 µm images for 2352 nearby galaxies, while Davis et al. (2019) presented one-dimensional, multi-component decompositions for 43 nearby spiral galaxies with measured SMBHs relying mostly on 3.6 µm Spitzer data. We include all the sample galaxies in common with Davis et al. (2019). Of the 101 galaxies in common between us and Salo et al. (2015), here we discuss 38. In selecting these 38 galaxies we try to create a representative sub-sample (in terms of the morphology and number of fitted galaxy components) of the 101 overlapping galaxies and to include all galaxies when there are notable disagreements with Salo et al. (2015).
IC 2574. The 1D HST surface brightness profile is well described by the exponential disc + Sérsic bar + two-parameter Gaussian nucleus model ( Fig. D.1). Salo et al. (2015) fitted a 2D exponential disc model to their 3.6 µm Spitzer data.
NGC 2273. We decompose the 1D HST surface brightness profile into five components (bulge, disc, bar, ring, and nucleus), which are best fitted with four Sérsic models and a Gaussian function. The HST surface brightness profile decomposition of this galaxy by Davis et al. (2019) shows a bulge, a disc, a bar, plus six, 3-parameter Gaussian components. They did not detect the nuclear component.  2012) showed the presence of multiple outer shell structures in the galaxy, which may indicate of recent gravitational interaction or merger. We therefore fit 1D and 2D Sérsic bulge + outer exponential shell + Gaussian nucleus model to the galaxy data.
NGC 2655. We fit the HST data with a Sérsic bulge + an outer exponential disc + a Gaussian nucleus model. Salo et al. (2015) fitted a 3-component bulge+bar+disc 2D model to their 3.6 µm Spitzer image. The best-fitting Sérsic indices from our 1D and 2D Sérsic bulge models are n ∼ 2.7, smaller than that from Salo et al. (2015, n ∼ 4.4). The latter value is high for a Sérsic galaxy.
NGC 2748. The 2D and 1D HST data is well fitted by a (n ∼ 2.2) Sérsic bulge, an outer exponential disc, plus a Gaussian nucleus model. Davis et al. (2019) fitted their major-axis HST surface brightness profile for the galaxy with a (n ∼ 1.59) Sérsic bulge+ Edge-on disc model and two, faint, Gaussian components. Salo et al. (2015) performed 2D decomposition of the 3.6 µm Spitzer galaxy image into an exponential disc plus a Gaussian nucleus.
NGC 2768. This galaxy is classified as an E6 in the RC3 (de Vaucouleurs et al. 1991), but it was classified as an S0 by Sandage & Tammann (1981). We fit a Sérsic bulge, an outer exponential disc, plus a Gaussian nucleus model to our composite HST+ground-based surface brightness profile. Salo et al. (2015) fit a 2D Sérsic bulge + outer exponential disc model to their 3.6 µm Spitzer image.
NGC 2770. The galaxy is classified as an SA(s)c in the RC3 (de Vaucouleurs et al. 1991). We therefore fit a Sérsic bulge, an outer exponential disc, and a Gaussian nucleus model, whereas Salo et al. (2015) reported the galaxy has a bar. They fitted the bulge component with an exponential model, in good agreement with our Sérsic (n ∼ 0.9) bulge model. NGC 2787. We decompose our 1D HST surface brightness profile in to bulge, a disc, a bar, an inner disc, and a nucleus, whereas the 2D bulge+disc+bar model fitted to the Spitzer image by Salo et al. (2015) did not include a nucleus and an inner disc. Our best-fitting Sérsic bulge model (n ∼ 1.3), compare to that from Salo et al. (2015, n ∼ 2.9).
NGC 2859. The 1D HST+ground-based light profile is well decomposed into six components: a bulge, a disc, two bars, a ring, and a nucleus (Figs. 6 and D.2). Our 2D decomposition of the galaxy's HST image identified a seventh component, that is an inner ring which encloses the inner bar. In their 2D decomposition of the galaxy's Spitzer image, Salo et al. (2015) fitted a bulge, an outer exponential disc, a secondary disc fitted by a Ferrers function and a bar.
NGC 2964. Our 1D and 2D decompositions show that the galaxy has a bulge, a disc, a spiral-arm component, and a nucleus. We fit the spiral-arm component with a 3-parameter Gaussian function. Salo et al. (2015) fitted a 2D bulge+disc model to their Spitzer data. The bulge+disc model was also used by Ganda et al. (2009) to fit the 1D HST+ground based data for the galaxy.
NGC 3031. We fit the 1D HST+SDSS light profile with a (n ∼ 4) Sérsic bulge, an exponential disc, a weak Sérsic Bar and a Gaussian nucleus. The 2D decomposition of the Spitzer 3.6 µm stellar distribution of the galaxy by Salo et al. (2015) show a (n ∼ 3.6) Sérsic bulge plus an outer, exponential disc components. In contrast, Davis et al. (2019) modelled the 1D Spitzer 3.6 µm major-axis surface brightness profiles of NGC 3031 with a core-Sérsic bulge + exponential disc model, they also added four, relatively faint Gaussian components.
NGC 3073. We decompose the 1D and 2D HST light distributions into four components (bulge-disc-bar-nucleus), fitting three Sérsic models and a Gaussian function. Salo et al. (2015) modelled the 2D Spitzer 3.6 µm stellar distribution of the galaxy with a Sérsic bulge and an outer exponential disc model. NGC 3077. This galaxy is classified as an I0 in the RC3 (de Vaucouleurs et al. 1991). We classify NGC 3077 as 'bulgeless' and fit its 1D and 2D HST(+SDSS) data with an outer (n ∼ 1.8) Sérsic disc, a Sérsic inner disc and a Gaussian nucleus model. Salo et al. (2015) modelled their 2D Spitzer 3.6 µm stellar distribution with a single, (n ∼ 3.2) Sérsic bulge model. NGC 3079. We decompose the HST+SDSS surface brightness profile into a (n ∼ 2.3 Sérsic bulge)-disc-bar-(spiral-arm)nucleus model profile. The 1D decomposition of the Spitzer 3.6 µm major-axis data by Davis et al. (2019) was performed using a (n ∼ 0.5) Sérsic bulge, an edge-on disc, a Ferrers bar, plus three, additional, Gaussian components. Salo et al. (2015) modelled the 2D Spitzer 3.6 µm stellar distribution of the galaxy with a 2D (n ∼ 1.7) Sérsic disc plus an edge-on exponential disc model. 2D Spitzer 3.6 µm stellar distribution with a (n ∼ 1.4) Sérsic, an outer disc, and a nuclear component.
NGC 3198. This galaxy is a barred spiral galaxy (de Vaucouleurs et al. 1991). We find the galaxy's 1D HST+SDSS light profile is well fitted with a five-component bulge (n ∼ 3.4)+disc+bar+spiral-arm+nucleus model. In contrast, Salo et al. (2015) modelled the 2D Spitzer 3.6 µm stellar distribution with a bulge+disc 2D model and the resulting best-fitting bulge Sérsic index ( n ∼ 10.0) is unusually large for a Sérsic galaxy.
NGC 3486. This galaxy is well modelled by a fourcomponent bulge+disc+ spiral-arm+nucleus decomposition model. Salo et al. (2015) fitted the 2D Spitzer 3.6 µm stellar distribution of the galaxy with a bulge and two exponential disc components. The Sérsic index of the bulge from our 1D and 2D decompositions is n ∼ 2.4, in close agreement with that from Salo et al. (2015, n ∼ 2.0).
NGC 3504. The surface brightness profile of this galaxy is well described by the bulge+disc+ring+bar+spiral-arm+nucleus decomposition model. Salo et al. (2015) fitted a 2D, threecomponent bulge+disc+bar decomposition model to their Spitzer 3.6 µm data.
NGC 3613. This galaxy is classified as an E6 in the RC3 (de Vaucouleurs et al. 1991). However our decomposition of the 1D HST +SDSS surface brightness profile for the galaxy reveals a core-Sérsic bulge and an outer (n ∼ 0.8) Sérsic disc. Trujillo et al. (2004) also modelled their 1D HST brightness profile and identified a 'possible core' in the galaxy.
NGC 3665. While the HST NICMOS+SDSS surface brightness profile is well fitted by a (n ∼ 3.9) core-Sérsic bulge and an outer exponential disc model, the nuclear dust in the galaxy implies our core identification is tentative. Laurikainen et al. (2010) performed a 2D decomposition of to their ground-based data for the galaxy into a (n = 2.7) Sérsic bulge plus an outer exponential disc model. NGC 3718. Our 1D and 2D five-component (bulge + outer disc + inner disc + spiral-arm + nucleus) decompositions of the HST data yield a bulge Sérsic index of n ∼ 2. In contrast, the 2D, three-component (Sérsic bulge + exponential disc + Ferrer disc) decomposition of the Spitzer 3.6 µm data by Salo et al. (2015) yielded a bulge Sérsic index of n ∼ 9.0, which is unusually high for a Sérsic galaxy.
NGC 3884. This galaxy is classified as an SA(r)0/a in the RC3 (de Vaucouleurs et al. 1991), we therefore fit a (n ∼ 2.7) bulge + disc + ring + nucleus model to the 1D HST light profile and 2D HST image. Because Dong & De Robertis (2006) modelled their 2MASS data in 2D with a single, Sérsic bulge model and did not account for the disc and ring components, the resulting best-fitting bulge n is biased high (i.e. ∼ 3.8).
NGC 3898. Decomposing the 1D HST+SDSS light profile of NGC 3898 into a bulge, a disc, a spiral-arm component, and a nucleus, we find n ∼ 3.4 for the bulge. In contrast, Salo et al. (2015) fitted a, 2D, (n ∼ 5.2) Sérsic plus an exponential disc model to the Spitzer 3.6 µm data.
NGC 3949. We describe the 1D HST surface brightness profile and the 2D HST image for the galaxy as the sum of a (n ∼ 1.0 − 1.3) Sérsic bulge, an exponential disc, a Gaussian spiral-arm, and a Gaussian nucleus. Our decomposition can be compared to Ganda et al. (2009) who fitted a (n ∼ 1.1) bulge plus an exponential disc model to their 1D HST+ground-based light profile.
NGC 3982. This galaxy has a weak bar component. In our (bulge + disc + bar + spiral-arm + nucleus) decomposition, the bar was fitted with a (n ∼ 0.2) Sérsic model. In contrast, Salo et al. (2015) fitted a two-component (exponential disc + Gaussian nucleus) decomposition model to their Spitzer 3.6 µm data.
NGC 4026. This galaxy has an inner stellar disc (Gültekin et al. 2009). In our decomposition, we fitted the HST surface brightness profile as the sum of a (n ∼ 1.3) Sérsic bulge, an outer exponential disc, two (n ∼ 0.3) Sérsic discy components, and a Gaussian nucleus.
NGC 4036. We find that the three-component bulge-discnucleus light profile of NGC 4036 is well fitted by the (n ∼ 1.4) Sérsic bulge + exponential disc + (n ∼ 0.4) Sérsic nucleus decomposition model. Similarly, in their 2D decomposition of the galaxy's i-band SDSS data, Beifiori et al. (2012) fitted a (n ∼ 1.7) Sérsic bulge plus an outer exponential disc model. NGC 4041. Our 1D and 2D decompositions reveal that this galaxy has a (n ∼ 0.9) Sérsic bulge, an outer exponential disc, and a (n ∼ 0.8) Sérsic nucleus. Our structural analysis can be compared to the three-component, 2D fit by Salo et al. (2015): a (n ∼ 0.9) Sérsic disc, an outer exponential disc, and a Gaussian nucleus.
NGC 4062. We describe the 1D HST light profile and 2D HST stellar light distribution of NGC 4062 as the sum of a (n ∼ 0.4 − 0.7) Sérsic bulge, an exponential disc and a Gaussian nucleus. Salo et al. (2015) modelled the galaxy's 2D, 3.6 µm galaxy stellar light distribution with a (n ∼ 0.7) Sérsic bulge plus an outer exponential disc.
NGC 4096. This galaxy is classified as an S0 in the RC3. We did not identify a bar, thus we fit a 1D/2D bulge-disc-(spiral-arm) model to our HST light profile/image. In contrast, Salo et al. (2015) fitted the 3.6 µm galaxy light distribution with a outer disc, a bar, and a nucleus.
NGC 4144. We fit this galaxy's HST surface brightness profile with a (n ∼ 0.9) Sérsic bulge and an exponential disc decomposition model. Salo et al. (2015) identified our bulge component as a bar and modelled the 2D, 3.6 µm stellar light distribution for the galaxy with a Ferrers bar plus an outer exponential disc.
NGC 4151. The galaxy has a bright point like Seyfert nucleus (Williams et al. 2017(Williams et al. , 2020) and a large-scale oval bar, an inner disc-like feature seen in molecular gas and optical dust maps (Dumas et al. 2010). The inner disc-like component is likely due to the gaseous response to the primary bar Asif et al. 2005). Inside this inner gaseous feature is an outflowing ionised gas from the AGN (Kishimoto et al. 2022). Wang et al. (2010) reported that the galaxy might have experienced an Eddington-rate outburst ∼ 10,000 yrs ago.
We decompose the 1D HST+SDSS brightness profile of NGC 4151 into a (n ∼ 3.9) Sérsic bulge, an outer exponential disc, plus a bar, an inner disc, a spiral-arm component, and a nucleus. In their five-component decomposition of the 1D Spitzer 3.6 µm light profile, Davis et al. (2019) fitted a (n ∼ 2.2) bulge, a disc, a bar, plus two Gaussian components. On the other hand, Salo et al. (2015) decomposed the Spitzer 3.6 µm galaxy light distribution in 2D with a (n ∼ 1.8) bulge and an outer exponential disc.
NGC 4203. This galaxy is classified as an SAB0 in the RC3. While Salo et al. (2015) fitted a bar component in their 2D decomposition, it was too weak to be included in our decomposition. We therefore describe the HST brightness profile as the sum of a Sérsic bulge, an exponential disc, and a Gaussian nucleus, while Salo et al. (2015) fitted a 2D (Sérsic bulge)-(exponential disc)-(Ferrers bar) model to their Spitzer 3.6 µm data. Our bestfitting bulge Sérsic index (n ∼ 2.5) is in fair agreement (within the quoted 12% error, Table B.2) with that from Salo et al. (2015, n ∼ 3.0).
NGC 4245. The five-component light profile for this galaxy is well described by the (n ∼ 1.7) Sérsic bulge + exponential disc + Sérsic bar + Gaussian ring + Sérsic nucleus decomposition model. Our decomposition can be compared to the 2D decomposition by Salo et al. (2015) who fit a three-component decomposition model: (n ∼ 1.1) Sérsic bulge + exponential disc + Ferrers bar.
NGC 4258. We describe our 1D HST brightness profile for this galaxy as the sum of a (n ∼ 3.6) Sérsic bulge, an exponential disc, plus a bar (Siopis et al. 2009), a spiral-arm, and a nuclear component. Similarly, the decomposition of the galaxy's Spitzer 3.6 µm data by Davis et al. (2019) shows a (n ∼ 3.2) Sérsic bulge, an exponential disc, plus a Ferrers bar, and Gaussian nucleus, but they did not account for the spiral-arm component.
NGC 4448. In our 1D and 2D five-component (bulge+disc+bar+spiral-arm+nucleus) decomposition of this galaxy's stellar light distributions, the bulge is described by a (n ∼ 2) Sérsic model, while in the bulge+disc+bar decomposition by Salo et al. (2015) the bulge was fitted with a (n ∼ 0.9) Sérsic model. NGC 4449. Our 1D and 2D bulge+disc+bar+nucleus decomposition of the HST data for this galaxy reveals a small (R e ∼ 0.13 − 0.16 kpc) bulge which is well fitted by a (n ∼ 2.1 − 3.2) Sérsic model. This is contrary to the fit in Salo et al. (2015) who modelled the 2D Spitzer 3.6 µm stellar distribution of the galaxy with the exponential disc + Ferrers bar + Gaussian nucleus decomposition model.
NGC 4565. Our multi-component decomposition of the HST data for this galaxy shows a bulge, an outer disc, an inner disc, a boxy bulge, and a nucleus. Our fit agrees well with the 1D HST surface brightness decomposition in Kormendy & Barentine (2010), although they did not account for the faint, inner disc, and nucleus (see also Laurikainen & Salo 2016).
NGC 4605. We model the 1D and 2D light distributions of the galaxy with a (n ∼ 0.3 − 0.5) Sérsic bulge, a nucleus, and an outer exponential disc model. Salo et al. (2015) fitted a (n ∼ 0.7) Sérsic disc plus an outer exponential disc decomposition model to their 3.6 µm Spitzer data.
NGC 4648. This galaxy is classified as an E3 in the RC3 (de Vaucouleurs et al. 1991), therefore we describe the 1D and 2D HST data as a sum of a (n ∼ 2.6 − 2.8) Sérsic bulge, an outer exponential halo, a Sérsic ring, and a Gaussian nucleus. Baldassare et al. (2014) modelled the two-dimensional HST light distribution of the galaxy with an inner Sérsic component plus an outer exponential disc.
NGC 4800. This galaxy show a near-exponential (i.e. n ∼ 1.3) bulge, an outer exponential disc, a circumnuclear ring (Comerón et al. 2008), a spiral-arm component, and a Gaussian nucleus. Salo et al. (2015) fitted a two exponential discs plus a Gaussian nucleus.
NGC 4826. The 1D HST+SDSS light profile is well described as a sum of a Sérsic bulge, an exponential disc, a Gaussian ring, and a Gaussian nucleus. The Sérsic index of the bulge is n ∼2.3. Salo et al. (2015) modelled their 3.6 µm Spitzer data for the galaxy with a (n ∼ 4.2) Sérsic bulge and two exponential discs, whereas Davis et al. (2019) modelled their 1D (majoraxis) 3.6 µm Spitzer light profile with a (n ∼ 0.7) Sérsic bulge, an outer exponential disc, plus five Gaussian components.

5005.
While this SAB(rs)bc galaxy (de Vaucouleurs et al. 1991) is reported to have a bar by Salo et al. (2015), it is not visible in our structural analysis. Also, the multi-wavelength structural analysis by Richards et al. (2015) did not detect a bar component in the galaxy. We fitted the 1D HST+SDSS light profile as the sum of a bulge, an outer disc, an inner clumpy disc, a spiral arm, and a nucleus.
NGC 5055. In our decomposition of this galaxy's 1D and 2D HST+SDSS data, we fit a (n ∼ 0.8 − 1.3) Sérsic bulge, an outer exponential disc, plus a spiral-arm component, and a nucleus which are described by two Gaussian functions. In their 2D decomposition of the galaxy's Spitzer 3.6 µm data, Salo et al. (2015) fitted two exponential discs plus a Gaussian nucleus. On the other hand, Davis et al. (2019) modelled the 1D major-axis Spitzer 3.6 µm data with a (n ∼ 2.2) Sérsic bulge, an exponential disc, and a Gaussian nucleus.
NGC 5308. Seifert & Scorza (1996) reported the presence of a rapidly rotating inner disc component in this galaxy. In addition, the kinematic maps by Emsellem et al. (2004) reveal that the central dispersion structure in the galaxy is box-shaped. Our decomposition shows a near-exponential (i.e. n ∼ 1.4 Sérsic) bulge and a (n ∼ 0.9) Sérsic nuclear disc. We also fit an intermediate-scale component which we refer to as a 'boxy bulge'.
NGC 5353. The 1D HST+SDSS surface brightness profile is well described as the sum of a (n ∼ 1.3) Sérsic bulge, an exponential disc, and an X-shaped Sérsic bar (e.g. Laurikainen et al. 2011).
NGC 5377. Our six-component decomposition of this galaxy's 1D HST+SDSS light profile reveals a (n ∼ 2.2) Sérsic plus outer disc, an inner disc, a bar, a spiral-arm component, and a nucleus. Decomposing the galaxy's 2D Spitzer 3.6 µm stellar distribution, Salo et al. (2015) reported a (n ∼ 2.2) Sérsic bulge, an exponential disc, and a Ferrers bar.
NGC 7741. This galaxy is classified as a barred S(s)cd in RC3. The bulge component in our 1D and 2D three-component (bulge-disc-spiral-arm) decomposition may be a bar. Salo et al. (2015) modelled their Spitzer 3.6 µm stellar light distribution of the galaxy in 2D with a Ferrers bar and exponential disc decomposition model. Tables   Table B.1 lists the tables, which present the best-fitting structural parameters for our sample of 163 newly analysed LeM-MINGs galaxies and the corresponding figures which show the decompositions. Tables B.2−B.4 list the best-fitting structural parameters from the surface brightness profile decompositions. Table B.5 provides global and central properties of sample galaxies including distance, morphological classification, velocity dispersion, bulge and galaxy stellar masses, optical and radio luminosities, ellipticity, isophote shape parameter, and inner logarithmic slope of the galaxies' inner light profiles.
(2) indicates the tables, which list the galaxy structural parameters for our sample of 163 newly analysed LeMMINGs galaxies and the associated figures where the fits are shown. The full table is available at the CDS. Notes.-Best-fitting model parameters for LeMMINGs galaxies from our 1D (major-axis) and 2D decompositions (see Figs. D.1 and E.1). Unless labelled otherwise, we present 1D best-fitting parameters. Shown are the galaxy name, the HST instrument and filter, the goodness of the fit, the number and nature of the fitted components, and the ground-based data at large radii for galaxies with composite profiles, that is inner HST plus outer ground-based data. Bulgeless galaxies are indicated by the superscript ' †'. Quality flag '1' indicates fits with good or higher quality levels whereas those labelled '2' are questionable primarily as a result of the difficulty in modelling highly inclined (i 75 • ) galaxies, see the text for further detail. For each 2D galaxy component, we present the best-fitting imfit ellipticity ( ) and position angle (P.A.). The imfit 'boxy' and 'discy' parameter c 0 is shown for 2D fits that were performed using generalised ellipses (Erwin 2015). We do not show the errors on the best-fitting 2D parameters from imfit as they are unrealistically low, see Erwin (2015). The fitted models (see  Notes.-Similar to Table B.2 but here showing Sérsic galaxies with a large-scale Gaussian component (Figs. D.2 and E.1). The fitted 3-parameter Gaussian ring model parameters (see Section 3.1) are µ r , σ, and R r , where µ r denotes the brightest value of the Gaussian ring with a semi-major axis R r and a full width at half maximum of FWHM = 2 √ 2ln2 σ. The errors associated with the best-fitting parameters for µ r , σ, and R r are 0.28, 15%, and 12%, respectively. The full table is available at the CDS. Notes.-As Table B.2 but here for core-Sérsic LeMMINGs galaxies (see Figs. D.3 and E.1). A '?' shows that we tentatively classify the lenticular galaxies NGC 1167 and NGC 3665 as core-Sérsic. The fitted core-Sérsic model parameters (see Section 3.1) are , R e [kiloparsec], γ, and α. The errors associated with the best-fitting core-Sérsic parameters are 0.15, 5%, 5%, 10%, 10%, 10%, and 10%, respectively. The full table is available at the CDS.

Appendix C: HST images
For each galaxy newly studied in this paper, we show a 15 ×15 HST image (Fig. C.1).    Table B.2; galaxies with a large-scale Gaussian ring component are shown in Fig. D.2, whereas core-Sérsic galaxies are in Fig. D.3). The residual profiles along with the rms residual ∆ for each fit are shown. The ellipticity, position angle (P.A., measured in degrees from north to east), and isophote shape parameter (B 4 ) profiles (red boxes) are given in the lower panels. The HST filter associated with each galaxy profile is listed in Table B.2. The surface brightnesses (µ) are in units of mag arcsec −2 . The dashed red curves represent the bulges, while the dotted blue curve shows the large-scale discs or stellar halo light which we model with an exponential function. Nuclear components (i.e. AGN, NSCs, inner discs, and bars) were typically fit with the two-parameter Gaussian model (dash-dot-dot-dot green curve). 'SNucleus' ('SNuc') (dot-dashed green curve) and 'SDisc' (dot-dashed blue curve) denote nuclei and large-scale disc components that were fitted with Sérsic models, respectively. Galaxy components such as bars, small-scale discs, rings, spirals, and lenses are described by Sérsic models (i.e. orange and magenta dot-dashed curves). For two galaxies in our sample (IC 239 and IC 520) their spiral arm components were modelled with a Sérsic model (i.e. 'SSarm'). The complete fits (solid orange curves) fit the observed galaxy profiles with a median rms residual

Appendix E: 2D decomposition
Fig. E.1 shows the residual images from our 2D fitting which were generated after subtracting the 2D imfit model images from the original HST images for 28 sample galaxies, which are representative of the 65 sample galaxies with 2D decompositions, listed in Tables B.1-B.4. We note that the 2D multicomponent modelling of the full HST images for the 65 galaxies are performed using the same type and number of galaxy structural components as the corresponding 1D modelling. For the 2D fits, the 2D models generated using imfit have the same type and number of galaxy structural components as the corresponding 1D . This exercise has a twofold merit: to estimates of the uncertainties on the best-fitting structural parameters and to perform an MC-based estimates of the structural parameters for the galaxy components . Panels show the histograms of the best-fitting parameters for the bulge, disc, and second Sérsic component from fits to the realisations for 50 LeMMINGs galaxies (see the text for details). Gaussian model fits to the best-fitting parameter distributions (red curves) together with the associated mean and standard deviation values are shown. These mean values are MC-based measurements of the best-fitting structural parameters for the galaxy components, which overall are in good agreement with those adopted in this work indicated by dashed vertical line (see Tables B.2  This section describes our considerable endeavour dedicated to explore the robustness of the decomposition of the surface brightness profiles and the accuracy of the measured best-fitting structural parameters using Monte Carlo (MC) simulated profiles. Doing so has a twofold merit: it provides realistic error estimates for the best-fitting structural parameters and allows us to derive a MC-based galaxy component structural parameters. More than 100 realisations of each galaxy's light profile were generated by running a series of MC simulations. To achieve this, for each galaxy, we perturb the data points of the galaxy light profile by sampling from a Gaussian with sigma, where the sigma value is the mean of all the residuals obtained from the actual galaxy light profile fitting (see Appendices B and D) We then add a correlated noise (i.e. a single, constant value) to the entire profile to account for potential errors from inaccurate sky subtraction. The standard deviation around the mean of the median sky background values measured in Section 2.2 were consulted to roughly simulate potential sky subtraction errors. Each realisation was then decomposed following the exact same fitting procedure as the modelling of the original galaxy light profile, including PSF convolution.
In Fig. F.1, we display the distributions of best-fitting parameters for the bulge, second galaxy component, and the outer disc for 50 LeMMINGs galaxies together with Gaussian models fitted to the parameter distributions. The adopted errors for the best-fitting parameters (Tables B.2−B.4) are computed using the standard deviations (i.e. 1σ values) from the fitted Gaussians to the 100 MC best-fit values (i.e. assuming a normal distribution) together with the difference between the mean values of the Gaussians and the best-fitting parameters adopted in this work. We note that the median 1σ values (see Fig. F.1) associated with the bulge's µ e , R e , and n for the LeMMINGs galaxies are ∼ 0.20, 10%, and 10%, respectively. We also remind the reader that the errors we derive here are consistent within the framework of the Sérsic models. However the total errors could modestly exceed those estimated from the model fits as the fitting functions are not perfect and cannot possibly reproduce all structural features in nearby objects such as those dealt with in this work. We note that the mean values of the Gaussians (see Fig. F.1) are the MCbased best-fitting structural parameters for the galaxy components, which overall are in good agreement with those adopted in this work. The application of the galaxy decomposition methods on the simulated light profiles reveal that they are fairly robust.
To quantify the errors on the bulge magnitudes and stellar masses (Table B.5), the uncertainties on the associated bestfitting parameters were propagated. Uncertainties in M/L are also considered in the stellar mass error budget (see Section 3.3).