Composite Bulges. III. A Study of Nuclear Star Clusters in Nearby Spiral Galaxies

We present photometric and morphological analyses of nuclear star clusters (NSCs)—very dense, massive star clusters present in the central regions of most galaxies—in a sample of 33 massive disk galaxies within 20 Mpc, part of the “Composite Bulges Survey.” We use data from the Hubble Space Telescope including optical (F475W and F814W) and near-IR (F160W) images from the Wide Field Camera 3. We fit the images in 2D to take into account the full complexity of the inner regions of these galaxies (including the contributions of nuclear disks and bars), isolating the NSC and bulge components. We derive NSC radii and magnitudes in all three bands, which we then use to estimate NSC masses. Our sample significantly expands the sample of massive late-type galaxies with measured NSC properties. We clearly identify NSCs in nearly 80% of our galaxies, putting a lower limit on the nucleation fraction in these galaxies that is higher than previous estimates. We find that the NSCs in our massive disk galaxies are consistent with previous NSC mass–NSC radius and galaxy mass–NSC mass relations. However, we also find a large spread in NSC masses, with a handful of galaxies hosting very low-mass, compact clusters. Our NSCs are aligned in PA with their host galaxy disks but are less flattened. They show no correlations with bar or bulge properties. Finally, we find the ratio of NSC to BH mass in our massive disk galaxy sample spans a factor of ∼300.


INTRODUCTION
Most galaxies today have a central massive object that is indicative of current or past extreme activity at their centers.These central massive objects can be a supermassive black hole (SMBH), a nuclear star cluster (NSC), or a combination of both.A review of NSC properties has recently been compiled by Neumayer et al. (2020).NSCs are extremely luminous objects that are present in the centers of a majority of all types of galaxies (e.g., Phillips et al. 1996;Carollo et al. 1998;Böker et al. 2002;Côté et al. 2006;Georgiev et al. 2009).They are compact, massive star clusters with an effective radius ranging from 3 -20 pc (e.g., Böker et al. 2004;Côté et al. 2006;Walcher et al. 2006;Georgiev & Böker 2014;Georgiev et al. 2016) and dynamical and stellar populaton based masses from 10 5 -10 8 M ⊙ (e.g., Walcher et al. 2005;Rossa et al. 2006;Erwin & Gadotti 2012;Spengler et al. 2017;Nguyen et al. 2018).NSCs are known to be more massive and denser than GCs (Walcher et al. 2005;Hopkins & Quataert 2010).
NSCs are mostly frequently found in intermediatemass galaxies (> 10 8−10 M ⊙ ) with a nucleation fraction (the fraction of galaxies studied that host an NSC) of ∼ 70-90% in the early-type galaxies (Côté et al. 2006;Turner et al. 2012;den Brok et al. 2014;Sánchez-Janssen et al. 2019;Hoyer et al. 2021) and > 75% in the latetypes (Böker et al. 2002;Seth et al. 2006;Georgiev & Böker 2014;Neumayer et al. 2020;Hoyer et al. 2021) in this same mass range.More recent studies of earlytype galaxies have shown a strong mass dependence and milder environmental dependence on the nucleation fraction.Sánchez-Janssen et al. (2019), Zanatta et al. (2021), and Carlsten et al. (2022) found the nucleation fraction to be as high as 90% for galaxy stellar masses of ∼10 9 M ⊙ , decreasing towards both lower and higher mass galaxies.Neumayer et al. (2020) and Hoyer et al. (2021) found a similar trend for late-type galaxies at lower masses, but at higher masses (> 10 10 M ⊙ ) a lack of data means we do not know if the nucleation fraction of late-types decreases in the same way as early-types.
NSCs are located at the dynamical centers of their galaxies (Neumayer & Walcher 2012).
NSCs and SMBHs are known to co-exist, including in the Milky Way (e.g., Seth et al. 2008a;Graham & Spitler 2009;Seth et al. 2010;Georgiev et al. 2016;Nguyen et al. 2019;Neumayer et al. 2020).The relative masses of NSCs and SMBHs appears to span a wide range, with NSCs being more massive than SMBHs in many lower mass galaxies, while SMBHs are the dominant component in higher mass galaxies (Graham & Spitler 2009;Neumayer et al. 2020).The relationship between NSCs and SMBHs is complicated and unclear.NSCs may provide a seeding mechanism to create and/or grow massive black holes at the centers of galaxies (e.g.Stone et al. 2017;Inayoshi et al. 2020).A possible consequence of NSCs and SMBHs co-existing is the presence of Tidal Disruption Events (TDEs).This occurs when tidal forces pull apart a star as it approaches the SMBH.TDEs are observed mostly in low-mass galaxies with some intermediate-mass galaxies (galaxy stellar masses between 10 9 -10 10 ; Wevers et al. 2019).
The formation history of NSCs can be directly probed through their stellar populations.Spectroscopic analyses of NSCs show they have have multiple stellar populations, extended star formation histories, and strong rotation (e.g., Walcher et al. 2006;Rossa et al. 2006;Seth et al. 2006;Kacharov et al. 2018;Pinna et al. 2021;Fahrion et al. 2021;Hannah et al. 2021;Fahrion et al. 2022b).Another way of studying the formation mechanisms of NSCs is through their correlation with the SMBHs at their centers and their surrounding host galaxies.Mass measurements of NSCs and host galaxies have shown that these quantities are strongly correlated; initially these were thought to form a relation sismilar to SMBHs (e.g.Ferrarese et al. 2006;Wehner & Harris 2006;Balcells et al. 2007;Graham 2012;Scott et al. 2013); however, more recent work shows the NSC and SMBH scaling relations are quite different (e.g.Erwin & Gadotti 2012;Georgiev et al. 2016;Sánchez-Janssen et al. 2019).
An in-depth analyses of NSC scaling relations in different galaxy environments can help us understand what physical mechanisms play an important role in NSC formation and how NSCs impact the overall evolution of the host galaxy.Evidence from previous studies (e.g., Hartmann et al. 2011;Antonini et al. 2012;Neumayer et al. 2020;Fahrion et al. 2021Fahrion et al. , 2022a) ) show that there are two primary formation mechanisms that drive the growth of NSCs, (1) star cluster mergers (Tremaine et al. 1975;Gnedin et al. 2014), or (2) in situ formation (Bekki 2007;Antonini et al. 2015).The relative importance of these mechanisms seems to depend on the galaxy stellar mass.Specifically, recent work (Neumayer et al. 2020;Fahrion et al. 2021Fahrion et al. , 2022a,b) ,b) finds that the NSCs of low mass galaxies primarily grow from cluster mergers, while in higher mass galaxies the NSCs grow from in situ formation.A reflection of this transition seems to be present in the scaling relations as well, with the NSC masses scaling with the square root of their host galaxy masses at lower masses but steepening to a more linear relationship in higher mass galaxies (den Brok et al. 2014;Sánchez-Janssen et al. 2019).
The changing properties of NSCs as a function of galaxy morphology and the resulting implications for their formation are not yet fully understood.Overall, there have been more studies of early-type than latetype galaxies, especially towards higher masses where dust and bulge contributions make studying late-type NSCs more challenging.Therefore, most observations of NSCs in late-type galaxies have focused on lowermass galaxies; this includes both HST imaging to quantify structure (Georgiev et al. 2009;Böker et al. 2002Böker et al. , 2004;;Georgiev & Böker 2014;Carson et al. 2015;Hoyer et al. 2023a), and spectroscopic observations focused on kinematics and stellar populations (Walcher et al. 2005;Kacharov et al. 2018;Pinna et al. 2021;Fahrion et al. 2022b).These observations are broadly consistent with the NSC mass trends discussed in the previous paragraph.However, some differences have been suggested in NSCs in early vs. late-type galaxies.For instance, the compilation of data by Georgiev et al. (2016) found that the sizes of NSCs in massive early-type galaxy are ∼2× larger than those in late-type galaxies, although Neumayer et al. (2020) suggested this difference may only exist at the highest masses.Late-type galaxies also seem to show stronger rotation on average than earlytype galaxies (Pinna et al. 2021).
Some studies exist of higher-mass, more bulgedominated late-type galaxies as well, most notably the studies Carollo et al. (1998Carollo et al. ( , 2002)), which found that a a majority of massive spiral galaxies do host NSCs.However, the information available on these NSCs is quite heterogeneous (i.e.photometric bands, sizes) making them challenging to interpret together.Stellar population measurements of a subset of these galaxies by Rossa et al. (2006) found that these galaxies clusters tended to be older than those in lower mass late-type galaxies.The recent paper by Hoyer et al. (2023b) shows the promise of JWST for studying NSCs in massive galaxies due to its high angular resolution, ability to penetrate dust, and the broad spectral energy distribution measurements it can obtain.Despite these studies, we still know little about the populations of NSCs in massive late-type galaxies (like the Milky Way; MW hereafter), which leads to a lack of knowledge about the NSCs in these galaxies and their co-existence with the ubiquitous SMBHs in these galaxies.In this paper, we focus on NSCs in high-mass (log(M ⋆ /M ⊙ ) > 10) late-type galaxies, presenting a study of NSCs in 33 galaxies with high resolution, uniform imaging from the Composite Bulges Survey (CBS; Erwin et al. in prep).
Section 2 describes the galaxy sample selected from CBS, the high resolution HST data, and the other data used in this work.In Section 3, we describe in detail the morphological modelling process for the galaxies and NSCs, including derivation of morphological NSC parameters.We also discuss the quality of the fits and error estimates.In Section 4, we focus on the properties of the NSCs, including the nucleation fraction of our sample, NSC magnitudes and masses, and correlations between NSC internal properties.In section 5, we present scaling relations of the NSCs with their host galaxies.We also discuss briefly the presence of SMBHs in our galaxy sample.Section 6 summarizes our work in this paper.

DATA SAMPLE
The work shown in this paper is part of the Composite Bulges Survey (CBS; Erwin et al. in prep).This survey is aimed at a detailed analysis of the stellar morphology and populations in the inner 1-2 kpc of massive disk galaxies, using both HST optical and near-IR imaging and VLT-MUSE IFU data.
CBS is based on a mass-and volume-limited sample of disk galaxies.These galaxies were selected from the de Vaucouleurs et al. (1995, RC3) catalog, restricted to galaxies with distances ≤ 20 Mpc, stellar masses ≥ 10 10 M ⊙ , S0-Sbc morphologies, inclinations between 35 • and 60 • , declinations δ ≤ +20 • , and galactic latitudes |b| > 20 • .Based on these selections, a total of 53 galaxies were selected; these include a large number of Virgo Cluster galaxies.

NSC and MW-like sample selection
For the work shown in this paper we present NSC fits for a sub-sample of 33 galaxies from the parent CBS sample (hereafter the "NSC sample").This sub-sample of the CBS survey was chosen to include a complete set of 20 galaxies form a complete Milky-Way like sample (hereafter the "MW-like" subsample).Our MW is a special galaxy which hosts the nearest NSC that can be studied in unparalleled detail.Understanding how typical or atypical the NSC in the MW is with the NSCs in similar galaxies is important.Hence we create this sample of MW-like galaxies which represents all galaxies having spiral morphologies (t ≥ 1) and stellar masses from 10 10.4 to 10 11.1 M ⊙ from CBS.We present NSC fits for the complete sample of CBS galaxies that have fit these criteria.In addition to these galaxies, we fit the NSCs of 13 additional galaxies -these galaxies are a random subsample of galaxies for which we had available larger scale models; the full sample of CBS galaxies morphological fits will be presented in Erwin et al. in prep.The circles in Figure 1 show all the 33 of the NSC sample galaxies; of these the MW-like subsample are shown as squares and the rest of the CBS sample as triangles.As can be seen, a majority of the CBS galaxies that are not included in our NSC sample are S0 type galaxies, while the NSC sample represents a nearly complete set of the Sa-Sbc galaxies.The properties of the NSC sample galaxies are summarized in Table 1.

HST data from CBS
For modeling the central morphology of each galaxy, we use a consistent set of high resolution HST data obtained using the Wide Field Camera 3 (WFC3) in the UVIS and IR modes (Cycle 25, Proposal ID 15133).We use the full-field WFC3/IR image in the F160W filter with a total integration time of 600s.The optical images were restricted to the C1K1C aperture for efficient readout and data transmission.This aperture is a 1024 X 1024-pixel subarray of the full-field WFC3/UVIS images centered on the galaxy nuclei.We use the F814W filter with a total integration time of 500s and the F475W filter with a total integration time of 700s.All of our images are divided into four dithered exposures to provide sub-pixel sampling.We note that the galaxies more than fill both detectors in many cases.
All the individual exposures from each band are combined using the Python-based DrizzlePac code.We set the output image scale for the UVIS images to 0.03 ′′ /pixel and for the IR images to 0.06 ′′ /pixel with a pixfrac = 0.7 in all cases.The sky subtraction during the processing of these images is turned off; we estimate the sky background using larger-scale ground-based or Spitzer IRAC1 images as explained in detail in Section 3.1.

Other data used in this work
For the large scale galaxy fitting and to determine the F160W image sky levels, we use Spitzer IRAC1 (3.6µm) images.The images for the galaxies come from either the S 4 G survey with a final mosaic pixel size of 0.75 ′′ /pixel (Sheth et al. 2010) or archive-generated mosaic images (Watkins et al. 2022) for galaxies not in S 4 G, with a default archive mosaic pixel of 0.6 ′′ /pixel.The program IDs for the IRAC1 images are listed in Table 1.
For determining the background sky levels in the UVIS images (Section 3.1), we use either SDSS g and i images or (for galaxies without SDSS images) B, V , and I images from the Carnegie-Irvine Galaxy Survey (Ho et al. 2011), kindly provided by Luis Ho.

DISSECTING THE GALAXIES' STELLAR COMPONENTS
In this section we explain in detail the process of modeling each of the galaxy components for the NSC sample galaxies.We describe how we identify and create accurate models for the inner regions of the galaxy, especially the NSC component.The fits for each of the galaxies is done as a three-step process.First, we use the skysubtracted IRAC1 images to fit the large-scale components (mainly the bar and disc related components) as explained in section 3.2.Using this fit, we then use the HST images (first the wider field-of-view (FOV) F160W images and then the UVIS images) to fit for the inner components (e.g., bulge, nuclear disc).Once we have a good model for the overall structure of the galaxy, we then identify, refine and constrain the NSC component using mostly the HST UVIS images as explained in section 3.3.
We fit our galaxies using the fast, multi-component image fitting program imfit1 version 1.8 (Erwin 2015).This code creates 2-D models for each galaxy component using user-defined input parameters, adds them together, and then convolves the summed image with a user-supplied PSF image.These are fitted to the images with the default χ 2 minimization using a Levenberg-Marquardt (L-M) minimization algorithm.In some cases where we run into local minima issues while modeling the images, we also use the Nelder-Mead (N-M) minimization algorithm to explore a wider range of solutions.The per-pixel uncertainties are estimated from the data values using the Gaussian approximation to Poisson statistics.When fitting models to the IRAC1 images, we convert the pixel values in the latter to ADUs with an assumed A/D gain of 3.7; the units of the HST images are in electrons and require no scaling.
We use a wide range of 2-D functions to accurately fit our galaxies disk, bulge and bar components as well as their NSCs.The most commonly used components are listed in detail in Appendix B.

Sky subtraction
While the primary goal of this paper is the nuclear morphology of galaxies, our modeling still requires accurate sky background estimates so that we can correctly model the large scale components of the glaaxies.The galaxies in our sample are massive and nearby, and thus they are almost always larger in angular size than the FOV of our HST images.(This is nearly always true for the 160 × 160 ′′ FOV of WFC3-IR, and always true for our 41 × 41 ′′ WFC3-UVIS C1K1C images.)It is thus not possible to accurately estimate the sky back- ground from the HST images themselves.To determine reasonable estimates for the sky backgrounds, we match surface-brightness profiles from ellipse-fitting of the HST images to those from ground-or spaced-based images at similar wavelengths with larger FOV (i.e., large enough to determine the overall background outside the galaxy).
The surface brightnesses were measured with the iraf ellipse task using ellipses with position angles and ellipticities matching the galaxy main-disk orientation; on the non-HST images, we used masks reproducing the orientation and FOV of the HST images.The resulting HST profile was then scaled to match the corresponding larger-FOV profile, using measurements outside the region strongly affected by differences in the PSFs, including an additive component representing the unknown HST background.
For the UVIS F475W and F814W images, the reference images are SDSS g and i-band images, respectively, where we apply a mask mimicking the orientation and FOV of the HST images.In the case of galaxies lacking SDSS images, we make use of images from the Carnegie-Irvine Galaxy Survey (Ho et al. 2011), matching their I-band images to the F814W images and averaging the results from matching the B-and V -band images to the F475W images.For a small number of galaxies where the matching with ground-based profiles failed we use the average sky values obtained from the rest of the galaxies in the NSC sample: 9.14 electrons for F814W images and 13.56 electrons for the F475W images.For the HST F160W images we generally used Spitzer IRAC1 images as the reference images.(Exceptions were NGC 3351, where we used an H-band image from the 2MASS Large Galaxy Atlas (Jarrett et al. 2003), and NGC 1300 and NGC 4321, where we made use of H-band images from Grosbøl & Dottori (2012).)Just as with the optical images, we match the surface brightness profiles within the larger HST F160W FOV to that of the IRAC1 images.(See Erwin et al., in prep, for more details.) The sky background is incorporated into our imfit modeling using a constant (fixed value) FlatSky function.

Large scale fitting
We initially use the large scale IRAC1 images to fit the larger, outer components of the galaxy.For all of our galaxies, we started with the combination of a single Sérsic component and a single Exponential component, and then added components to best fit the galaxy.These models were fit using an appropriate PSF; we have used the official in-flight Pixel-Response-Function (PRF) images2 that are down-sampled to the appropriate pixel scale from the original scale of 0.24 ′′ /pixel at column, row = 129,129; the approximate central location of most of our galaxies.We also generate a bad-pixel mask to remove bright foreground stars and background galaxies; we do this by running SExtractor (Bertin & Arnouts 1996) on the image, and then scaling detected objects using circles (for stars) or ellipses (for background galaxies).Notable image defects are masked by hand.After running our initial fit, we then refine the fit by adding components based on (1) visible patterns in the residual images, (2) surface brightness profiles, and (3) changes in the ellipticity and position angle of the galaxy isophotes.For some galaxies, we find significantly better fits if we replace the initial Exponential with a Broken-Exponential component.(Figures 5 and 8 in Erwin et al. (2021) show examples of this for the galaxies NGC 4608 and NGC 4643.)We also incorporate information from the galaxy morphology (e.g.known rings or bars) from previous studies.We continue to add additional components to the model until the galaxy is well represented; specifically, when there are no clearly visible systematic residual patterns that could be fit with an additional component.Typically, this includes an exponential or broken exponential disk component, a bar component, and at least one bulge component.Occasionally, we also fit ring and spiral arm components, especially when this is critical for understanding the bar and bulge structures.These components are always an addition to the initial Sérsic + Exponential model, except for cases where we have replaced the initial Exponential with a BrokenExponential.
An example of this fitting procedure can be seen in the upper panel of Figure 2 which shows the IRAC1 image, best-fitting model and the residual ratio image for NGC 4689.

HST image fitting
Once we have the best-fit IRAC1 model, we then model the galaxy's inner components in more detail using the HST images, beginning with the F160W image.This includes components such as a boxy-peanut bulge, classical bulge, star forming disk, etc.. PSFs -For each model we fit to the HST data (in all three band images), we provide an appropriate PSF image for convolution.For the HST images, we generated the PSF images using the grizli software3 .An "empirical PSF" image from Anderson (2016) 4 is inserted into the each of the four individual exposures at the location of the galaxy center; these are then run through the same drizzling process used to prepare the final HST image (as explained in section 2.2).The final PSF image is then extracted from the combined, drizzled image.
Masks -Foreground objects and dust can prevent us from getting an accurate model fit to the data.We initially identify the foreground stars, background galaxies and other image defects using SExtractor (Bertin & Arnouts 1996).We then create a mask to flag all those pixels in the image.We also mask other regions of the galaxies that are not well fit by our models, these include spiral arms in many cases, as well as other nonaxisymmetric features that are challenging to model.
Dust extinction near the centers of our galaxies can also significantly impact the best-fit nuclear models.To mask dust features, we use UVIS F475W-F814W color maps to find reddened pixels.We mask pixels using the distribution of pixel values in unreddened regions, choosing a slightly different color threshold for each galaxy.We then translate this mask to the F160W image as well.In addition to masking reddened regions, we also mask star-forming regions in some cases using a blue color threshold.
Fitting the F160W images -We start the modeling of the F160W image by using the best-fitting IRAC1 model as an initial guess (this involves translating size and angular parameters to account for differences in the pixel scales and image orientation, and estimating the difference in intensity values for intensity parameters).
Since some galaxy components (i.e. the disk) from the IRAC1 model can extend well outside the F160W FOV, we model these components by holding the translated best-fit shape, orientation and size parameters fixed, and fit only for the intensity of these largest components.For components that are mostly or entirely within the F160W [or WFC3-IR] FOV, we use the IRAC1 values as initial guesses, but leave all parameters free.We then iterate by running imfit to determine the best-fit parameters, followed by inspection of the residuals, adding additional components to the model if needed and rerunning the fit.
In addition to examining the residuals, we also apply a more quantitative approach to adding components.Specifically, we add components when (1) we see residuals in the radial profile that are > 10% of the data value and when (2) the addition of a component improves (i.e., reduces) the Akaike Information Criterion (∆AIC; Akaike 1974) of the fit by 1000 or more relative to the original fit (as in Erwin et al. 2021).Often at this stage we include an initial NSC component that fits for the central excess light; we describe our final NSC determinations below.
An example of this fitting procedure can be seen in the middle panel of figure 2 which shows the F160W image, best-fitting F160W model image and the residual ratio image for NGC 4689.We obtain a decent F160W model for this galaxy with a preliminary NSC component that fits for the central excess light.
Translation to the F814W images -The F814W images are read out in a sub-array, and thus they span a smaller region (the inner 41 ′′ × 41 ′′ ) of the galaxy than the F160W images.They are also higher resolution than the F160W images, with an original pixel scale of 0.04 ′′ and a final processed image scale of 0.03 ′′ pixels (critically sampling the ∼0.′′ 07 FWHM PSF) and thus provide better constraints on the nuclear structure than the F160W images as long as the galaxies are not too dusty.We translate the best-fit F160W model to create the initial guess for the F814W model by (1) scaling the intensity parameters by the ratio of fluxes in a fixed angular area in the two images, and (2) scaling the sizes of components to the higher resolution pixel scale.For components that extend well beyond the edge of our chip, we hold their shape parameters fixed and fit just for their intensity parameters.As with the F160W images, we add components as necessary in the higher resolution F814W images; again this often includes an initial NSC fit to the central light excess.
Final NSC fitting -Once we have a good fit to the model that accurately represents the whole F814W band image, we focus on modeling the NSC component.We choose to use a Sérsic function to describe the NSC component, following previous work by several authors (e.g.Graham & Spitler 2009;Carson et al. 2015;Nguyen et al. 2017).
As a first step, we check for the presence of a nuclear excess in the galaxy.For this, we exclude any NSC component from the F814W model and use the remaining components as initial conditions for a new fit -this allows us to check whether any nuclear light can be fit through adjusting the parameters of the larger components.We inspect the residuals of these fits and the surface brightness profiles in the central 10-15 pixels for the presence of the excess light.In all galaxies we find there is central excess light; in many cases this is due to an NSC.However, we expect central light excess can also be due to an AGN component, in which case we expect the light to be unresolved.Therefore, as a next step, we add a PointSource component to the model in an attempt to fit the nuclear excess.
We inspect the residuals of the centers for any visible pattern that is not being modeled.If the central light excess appears to be resolved, we then replace the PointSource function with a Sérsic function.We consider this central component to be unresolved if the (i) ∆AIC between the PointSource and Sérsic function is < 3000 and (ii) the difference in the surface brightness profile residuals between the PointSource and Sérsic function fits is < 5%.We consider resolved sources to be NSCs.Table 2 indicates sources that are unresolved as well as the PointSource vs. Sérsic ∆AIC values.For some galaxies with known strong AGN, we use an additional PointSource component along with the NSC Sérsic component.The evidence for AGN and details on the fits are given in Appendix C. Out of the 33 galaxies we fit, we find four of them to have unresolved nuclear components (using the ∆AIC < 3000 constraint discussed above).We discuss these four sources in more detail in Section 4.1 below.
An example of this fitting procedure can be seen in the lower panel of Figure 2, which shows the HST F814W image, the best-fitting F814W model image and the residual ratio images for three models: the first with no central (NSC) component, the second with a central PointSource component and the the third with the best-fitting NSC Sersıc component.From these three residual ratio images, we can visually judge the presence of a NSC component.The superiority of the third model, with the absence of any residual nuclear excess is evident.The AIC values for each of the 3 models are given below their respective residual ratio images; the NSC-Sérsic model has AIC ∼ −32, 000.
In most cases, we use the F814W image as the primary filter for modeling the NSCs.Once we have a good fitting model in the F814W filter, we then translate all the components to the F475W image as well as the F160W image -we fit only for the intensity of these components, and keep the shape parameters fixed to scaled best-fit F814W values.This ensures that we measure accurate colors for the NSCs (and other components).Allowing the NSC shape parameters to be fit independently in each filter can result in unphysical colors when e.g. the NSC component is much larger in one filter than another.
For galaxies that are very dusty in their centers (10 out of the 33), we instead use the F160W image as the primary band to fit our NSCs, and perform the methods outlined above on the F160W image directly, without translating to F814W.The lower extinction in F160W helps us better fit the nuclear structure despite the lower spatial resolution.In these cases, we still fit the F475W and F814W images with shape parameters fixed to the best-fit F160W values in order to derive colors for the NSCs in all three filters.Once we have the best-fitting model in all three filters, we then calculate the magnitudes for all of the components -explained in detail in section 4.2.Table 2 includes the primary filter used for modelling the NSCs.
The top panel in figure 3 shows an example of the 1-D surface brightness profiles we create to inspect the quality of the fit.These profiles were derived using circular aperture photometry with python's photutils library.We ignore masked pixels in deriving the surface brightness of both the original image (green points) as well as the best-fitting model (blue solid line).The NSC Sérsic component (purple solid line) can be identified as the visible bump seen at the smallest radii.The bottom panel shows the residuals (model -data) in magnitudes; radii where the data is brighter than the model have positive values in this panel.We also show the residuals for models with no central component (shown in red) and a central PointSource component (shown in orange).Here, we can observe changes in the residuals due to improving the NSC model fits.The residual based on the best-fitting model has residual below 0.05 magnitudes all the way out to 10 ′′ .

Fit Quality, AGN, and Exceptions
Here we discuss the process of evaluating the goodness of the fits to the galaxy data.As mentioned in section 3.2 and section 3.3, our quality check and criteria for adding new components for the each fit to the galaxy are based on (1) improvements in the residuals and (2) improvements in the AIC.
Some of the galaxies in our NSC sample (especially in the MW-like subsample) have literature data that indicates they have an AGN.Detailed descriptions of these AGN components and associated X-ray sources are provided in Appendix A).For galaxies with a known AGN, we test whether an additional point source component needs to be included in the model as well, but find no cases where both a point-source and resolved NSC component provide a significantly improved fit.In one case, NGC 1566, the AGN is bright enough (and in fact satu- Dust can also prevent us from obtaining good fits to the NSCs.In two galaxies (NGC 289 and NGC 7177) we are unable to determine accurate NSC morphologies due to dust obscuring the centers of all three bands.
We rate the overall quality of our fits using a single number.The quality values are as follows: Quality 0: These galaxies are those for which we do not obtain good fits for the NSCs -this includes the galaxies discussed above, NGC 289, NGC 1566 and NGC 7177.No fits to these galaxies are shown in Table 3. Quality 1: Unresolved NSC components in four galaxies.When comparing a point source component to a Sérsic component fit the ∆AIC between these fits is ≲ 1000, giving minimal evidence that the nuclear component is resolved.The nature of these sources is discussed in Section 4.1.Quality 2: dusty galaxies where dust absorption is evident all the way to the center.However, we are able to estimate NSC properties using the F160W images in these galaxies.Quality 3: galaxies with complex central regions with poorly modeled structures (e.g., spiral arms) that result in large radial surface-brightness profile residuals (>0.1 mag).Quality 4: galaxies with some dust in the nuclear regions (but not crossing the center).Fits have surface brightness profile residuals <0.1 magnitudes.Quality 5: galaxies with no dust within the central 2 ′′ .Fits have surface brightness profile residuals <0.1 magnitudes.
There are 27 galaxies with quality ≥2.For the plots in the paper, we exclude all NSCs with Quality < 2, while Quality 2 fits are shown with open symbols.

Error Estimation
We estimate the errors on our imfit models using bootstrap resampling.This method allows us to capture asymmetric errors and the covariance of NSC parameters with the other fitted components.It also provides more accurate errors than those estimated from the Leavenberg-Marquardt algorithm.We performed 200 iterations of bootstrapping for each galaxy model.For the NSC models, the effective radii have median errors of ∼9%, while the median error on Sérsic indices are  ∼14%.For ellipticities above 0.05, the errors on the ellipticity are just ∼7%.

Nucleation Fraction
Of the 33 galaxies in our sample, we find unambiguous, resolved NSCs with radii <50 pc in 26.Previously measured NSCs have typical effective radii of ∼3 pc with a small tail towards larger sizes and a cutoff suggested at ∼50 pc by Neumayer et al. (2020).In three galaxies, the presence of a nuclear cluster can't be constrained due to dust opaque enough to obscure the nucleus even in F160W (NGC 289 and NGC 7177), or the presence of a very bright AGN (NGC 1566) as discussed in Section 3.4.In another four galaxies, the potential NSC components are unresolved: NGC 613 and NGC 1300 have a low luminosity AGN, and thus may not be NSCs.On the other hand in NGC 2775 and NGC 1440 the nuclear light appears to be stellar; specifically, the nuclear sources in both have very similar colors to the surrounding galaxy light, and the nuclear spectra (Ho et al. 1995) shows very little emission.Therefore it is likely these two galaxies host compact NSCs that we just cannot resolve with HST.The best-fit Sérsic r e are used here as upper limits; for NGC 2775 and NGC 1440 these are 3.87 pc and 4.65 pc or 0. ′′ 042 and 0. ′′ 048.Our ability to resolve NSCs in our CBS galaxies is complicated due to varying galaxy backgrounds, however, it does extend to more compact sources than these in some galaxies -the most compact clearly resolved NSC is in NGC 4377 with an r e of 1.91 pc or 0. ′′ 22; this is similar to the limit on resolving NSCs found for galaxies at similar distances by Côté et al. (2006).
Using only the unambiguous resolved NSCs, we get a nucleation fraction of 78.8 +6.2 −7.9 % (26/33) for the full sample, with errors calculated using the Wilson interval.The nucleation fraction is 68.4 +9.5 −11.4 % (13/19) for the MW-like subsample.On the other hand, we cannot exclude the presence of NSCs in any of our galaxies, thus the nucleation fraction in both samples could be as high as 100%.We note that our galaxy sample are all high mass (log M ⋆ ) > 10.1) and mostly late-type galaxies; if we take all late-types in our sample, we get 76.0 +7.4 −9.4 % (19/25) with NSCs.
Two previous measurements exist for the nucleation fraction of massive late-type galaxies.We took the data from Neumayer et al. (2020) and Hoyer et al. (2021) to find a comparable nucleation fraction for galaxies with (log M ⋆ ) > 10.1); from the Neumayer et al. (2020)   cleation measurements are nucleated.Thus we find a higher nucleation fraction than either study, with our results being consistent with the measurement in Hoyer et al. (2021).Our values are much higher than the nucleation fraction of ∼30% seen in early-type galaxies with similar mass (Neumayer et al. 2020;Hoyer et al. 2021).

Magnitude and color
For all the NSCs, we estimate the magnitude using the imfit makeimage program, which can calculate fluxes and magnitudes for each component in the model.To determine the magnitudes of the NSC in the 3 HST bands, we divide the total counts by the total exposure time; 600s for F160W band, 500s for F814W band and 700s for F475W band.We then use the following zeropoints: 24.6949 in F160W, 24.684 in F814W, and 25.801 in F475W5 .These are Vega based zeropoints, and thus all magnitudes listed here are in the Vega system.We correct for foreground extinction using the A F 814W values and conversions to the other two bands in Table 1 and the notes to that table.The extinction corrected NSC magnitudes in each filter band are presented in Table 3.
Figure 4 shows the color-magnitude and color-color (UVIS-IR) diagrams of the NSCs.Padova PARSEC 1.2S single-stellar population models (Bressan et al. 2012) with ages from 1-13 Gyr and at two metallicities are overploted.A majority of the galaxies are consistent with these models with modest extinction up to (A V ≲ 2).Of these, only two require populations younger than ∼10 Gyr -due to the age-extinction degeneracy it is not possible to separately constrain the ages and extinctions using our colors.Almost half of the galaxies fall redwards of these SSP models in the F814W−F160W color in a way that is inconsistent with the reddening vector (which assume R V = 3.1), in some cases by >1 magnitude.This offset in F814W−F160W color could be due to (1) a mismatch between the data and SSP models, or (2) issues with calculating the NSC magnitudes.To investigate this last issue, we measured aperture photometry within the center 0. ′′ 5, and compared these aperture colors to the model NSC colors.For clusters with bright, prominent NSCs, the aperture and integrated NSC values agreed.However, for fainter NSCs that make up a smaller light fraction of the galaxy, we see a blueward offset of up to 0.5 mags for the aperture magnitudes relative to the NSC model magnitudes.This blueward offset can be explained by a combination of the lower encircled energy in F160W and bluer surrounding bulge components.Overall, this test suggests we are accurately measuring the NSC colors with our model magnitudes.While some NSCs have clear evidence for significant dust absorption (i.e. the open circles), others appear dust-free, suggesting that the PAR-SEC models may under-predict the F814W−F160W colors of NSCs.

Mass
We determine the NSC masses using the F814W magnitude and color-M/L relations from Roediger & Courteau (2015).Specifically, we first convert our F814W magnitude to the Sloan i band and then convert our HST F475W−F814W color to g − i using relations derived from PARSEC 1.2S stellar population models (Bressan et al. 2012): i − F814W = 0.099 × (F475W − F814W) + 0.404 We then use the g − i color vs i band magnitude relation of Table 2 in (Roediger & Courteau 2015) to determine the i-band M/L ratio.The resulting NSC masses are presented in Table 3.For our NSCs the log(M NSC /M ⊙ ) values range from 5.7 to 8.74 with a median of 7.16 for all our NSCs and 7.26 for the NSCs in the MW-like subsample.We determine the errors for the derived NSC mass using the bootstrap sampling (see Section 3.5), recalculating the luminosities, colors, and derived M/Ls for each sample.

NSC size, mass relations
We use the derived color and mass to understand our NSCs and compare them with the available literature from both early-and late-type galaxies.In Figure 5 we plot the radius of the NSCs (in pc) versus their derived masses.The solid circles denote our NSC sample galaxies, with surrounding squares indicating the MW-like subsample.All galaxies are color coded into early-(red) and late-(blue) type.The dashed line show the relationships for early-and late-type galaxies from Georgiev et al. (2016).As has been found in many previous studies (e.g.Hopkins & Quataert 2010;Norris et al. 2014;Georgiev et al. 2016;Neumayer et al. 2020), the NSC radii correlate with their masses, i.e., massive NSCs have a larger radii.Our NSC sample improves the available literature, especially for the massive late-type galaxies (> 10 10 M ⊙ ).Overall, the masses and radii of these clusters agree well with the overall trend seen in previously published data (shown as small stars in Fig. 5).Georgiev et al. (2016) fit mass-radii relations and find that the NSCs in the late-type galaxies (blue-dashed line) are roughly two times smaller than the NSCs in early-type galaxies (red dashed line).Our data does not seem to support this difference; in particular most of the latetype galaxies in our sample fall above the blue dashed line, while all of the early-type galaxies fall below the red-dashed line.This weakens the previous literature findings that there is a difference in the mass-radius relationship in early and late-type galaxies.
The biggest NSCs (> 25 pc) in our sample are found in IC 2051, NGC 4380, NGC 4501, NGC 4699, NGC 5248 and NGC 7513.These largest objects may be ambiguous in their classification as NSCs, however, the continuity of the mass-radius relationship suggests these are related components.For NGC 7513, Carollo et al. (2002) found the NSC to be a very compact source unlike in our model.The difference in the NSC model is discussed in detail in section 4.6.
Four objects (NGC 613, NGC 1300, NGC 1440 and NGC 2775) are unresolved in our sample but nonetheless appear to be NSCs; these are shown as upper limits on the mass radius diagram.These galaxies fall on the compact side of the locus of previous NSC measurements and thus are significantly denser than typical NSCs.

NSC mass vs. ellipticity and Sérsic index
In Figure 6, we compare the NSC mass with the NSC Sérsic index (top panel) and ellipticity (bottom panel).
We see no correlation between NSC masses and Sérsic indices, indicating that the NSCs have a wide range of concentrations (see also Hoyer et al. 2023a).One of the possible reasons for a high Sérsic index (> 5) is the presence of multiple components within the NSC, as seen in previous studies of very nearby NSCs where these components can be resolved (e.g., Seth et al. 2006Seth et al. , 2010;;Nguyen et al. 2018).The Sérsic indexes of our NSCs range from 0.1-5.9 with a median value of 2.7.The median value of the Sérsic indices of the MW-like subsample is 2.2.We have no evidence for multiple component NSCs in our sample except for in NGC 4612, in which Note-(1) Galaxy name, (2) NSC position angle (PA), (3) ellipticity, (4) Sérsic index, (5) effective radius, (6) F475W magnitude, (7) F814W magnitude, (8) F160W magnitude, (9) logarithmic stellar mass estimated using the M/L ratio (see Section 4.3, (10) Photometric bulge to total ratio of the galaxy (see Section 5.2.1), (11) bar luminosity determined using the bar component magnitudes.All magnitudes are corrected for Galactic extinction.
we obtain a two component fit.We only plot the average ellipticity obtained from the two component NSC in Figure 6 and do not plot the Sérsic index of this NSC.More information about the NSC fit for this galaxy is provided in the Appendix C. A trend of higher ellipticities in higher mass NSCs was seen in early-type galaxies by Spengler et al. (2017), however, we find no correlation between the NSC mass and ellipticity (bottom panel in figure 6).The ellipticities of all our NSCs range from 0-0.5 with a median value of 0.14 (0.16 for the MW-like subsample).

NSC literature comparison
The NSCs of five of the galaxies from our sample -NGC 289, NGC 1566, NGC 4237, NGC 4612 and NGC 7513 -have been studied previously in the literature.Two of these galaxies (NGC 289 and NGC 1566) have NSCs presented in Carollo et al. (2002); due to the bright AGN component (in NGC 1566) and dust (in NGC 289), we are unable to obtain reliable, unambiguous NSC fits for them, and so we focus on the other three objects, below.
NGC 4237 -This unbarred Virgo Cluster spiral galaxy has been previously studied in Georgiev & Böker (2014).They determine the effective radius of the NSC to be 0. ′′ 07 with an F814W magnitude of 17.74.They model the NSC using multiple images from HST includ- ing F606W and F814W band.From our work, the NSC in this galaxy has an effective radius of 0. ′′ 08, with an F814W magnitude of 18.37.So while the effective radii are similar, our fit to the NSC is considerably fainter than Georgiev & Böker (2014).This is likely due to our more careful modeling of the galaxy background.We note we also are fitting higher resolution images than the wide field chip WFPC2 images fitted in Georgiev & Böker (2014).
NGC 4612 -This barred S0 Virgo Cluster galaxy was previously studied in Côté et al. (2006) and Spengler et al. (2017) as part of the ACS Virgo Cluster Survey.They determine the NSC in this galaxy to be unresolved with an effective radius of 0. ′′ 024.In our model, we find that the NSC is best fit by a two component model, a compact Sérsic surrounded by a larger exponential with a combined effective radius of 0.14 ′′ .This model is preferred to a point source model with a ∆ AIC of 5810.We can directly compare the F475W magnitudes of the sources; they find 17.73 (after conversion to Vega magnitudes), while our combined NSC has a magnitude of 17.66, thus these agree quite well.We note that our approaches differ significantly; Côté et al. (2006) fit 1-D profiles with a single Sérsic background galaxy models, while we fit a more sophisticated galaxy model and fit in 2D.
NGC 7513 -This barred galaxy was previously studied in Carollo et al. (2002).They modeled the NSC using NICMOS2 data in the F110W and F160W filters and determine the NSC in this galaxy to have an F160W magnitude of 18.3 and an effective radsius of 0. ′′ 06 (0.97 pc), slightly smaller than the 0. ′′ 075 pixels.We find a much bigger and brighter NSC component in this galaxy, with an effective radius of 0.45 ′′ (43.54 pc) and an F160W magnitude of 15.66; an unresolved point source or compact component provides a much worse fit to the nuclear regions.The Carollo et al. (2002) fits did not model the galaxy background at all, and we suspect that this methodological difference may be responsible for this discrepancy.

NSC-GALAXY RELATIONS
In this section, we discuss in detail the relation of the NSCs to the properties of their host galaxies.We also briefly discuss the relations between the NSC mass and supermassive black hole mass from Saglia et al. (2016) for 4 galaxies in our sample in Section 5.3.

Structural Parameters
In Figure 7 we compare the ellipticity (ε = b/a) and position angle (PA) of our NSCs relative to their host galaxies.For the host galaxies, we use the PA and ε from Table 1.It is important to note that the CBS selection criteria selects galaxies with inclinations between 35 • and 60 • .In the left panel of the figure, we see the NSCs have ellipticies equal to or smaller than the ellipticites of their host galaxies, suggesting that NSCs are typically rounder than their host galaxy disks.
The right panel of the figure shows the difference in PA between the NSCs and the galaxies.Here we plot the difference for all NSCs with ellipticities > 0.05 where we can robustly estimate the NSC PA.We see that the distribution does not appear to be uniform as would be expected if there was no correlation between the galaxy disks and NSCs, but instead most of the NSCs have PAs within 25 • of their host galaxies.Using Kolmogorov-Smirnov and Anderson-Darling tests, we can reject the relative PAs being drawn from a uniform distribution at high significance (p-values of 0.0072 and 0.00075).(δPA).We plot the difference for all those NSCs whose ϵ > 0.05.There is a clear preference for near-alignment between NSCs and their host galaxies' disks.
This result is similar to what is seen in our Milky Way and in edge-on galaxies, where the NSC and galaxy PAs are typically aligned (Seth et al. 2006(Seth et al. , 2008b;;Feldmeier et al. 2014).However, this correlation of PAs is not seen by Georgiev & Böker (2014).Their sample included a much wider range of inclinations, and less correlation between NSC and galaxy PAs would be expected to be visible in more face-on galaxies.
The correlation of NSC and galaxy PAs suggests that NSCs are flattened and aligned with their large-scale galaxy disks.This favors NSC formation from material in the disk -either from gas accretion followed by in situ star formation, or by formation and inspiral of young star clusters (e.g.Seth et al. 2006;Agarwal & Milosavljević 2011;Tsatsi et al. 2017), and is consistent with strong rotation seen in many NSCs (Pinna et al. 2021).It disfavors NSC formation from inspiral of a more spherical distribution of globular clusters (e.g.Tremaine et al. 1975;Hartmann et al. 2011).This result is in agreement with previous work that suggests NSC formation is dominated by in situ star formation in more massive galaxies (log(M ⋆ /M ⊙ )≳9) like those in the CBS sample (Neumayer et al. 2020;Fahrion et al. 2021Fahrion et al. , 2022a)).

Correlations of NSC Mass with Host Galaxy Properties
Previous studies have found scaling relations between the mass of the NSC and its host galaxy properties.This includes scaling relations with bulge luminosity, velocity dispersion, and total stellar mass (e.g., Ferrarese et al. 2006;Wehner & Harris 2006;Rossa et al. 2006).Initially, the NSC scaling relations were found to be similar to the BH scaling relations, but recent studies with more data over wider ranges of host galaxy properties have shown that, unlike BH scaling relations, the NSC mass correlates better with galaxy mass than bulge mass or stellar velocity dispersionf, and does not follow the same scaling relations (e.g., Erwin & Gadotti 2012;Georgiev et al. 2016;Sánchez-Janssen et al. 2019).
In Figure 8, we plot the NSCs mass against their host galaxy stellar mass.We observe that the bulk of our galaxies fall among the highest NSC masses as expected given their high galaxy stellar mass.All the NSCs in our sample have masses > 10 6 M ⊙ , with a median mass of 4.2 × 10 7 M ⊙ .The typical masses fall along the relation for late-type galaxies from Georgiev et al. (2016), with a tight cluster of points around the median sample mass.However there is also a very large, >2 order of magnitude spread in the NSC mass, with several significant low outliers including NGC 3351 (log M ⋆ = 10.56 and log M NSC = 6.16) and NGC 6744 (log M ⋆ = 11.07 and log M NSC = 7.01) and the very compact NSC in NGC 2775 (log M ⋆ = 10.98 and log M NSC = 6.98).This broad range of masses suggests a wide range of formation and evolutionary processes in the NSCs in these massive (mostly late-type) galaxies.As noted in Section 4.1, the nucleation fraction in our NSC sample is much higher than the nucleation fraction of early-type galaxies with similar mass (Côté et al. 2006;Neumayer et al. 2020), where binary BH mergers might have destroyed NSCs (e.g.Milosavljević & Merritt 2001).With the ongoing formation of NSCs in late-type galaxies (e.g.Walcher et al. 2006;Rossa et al. 2006), the absence of an NSC after a binary BH merger would likely be shortlived and the existence of low mass outliers may trace galaxies where NSCs have been destroyed in the relatively recent past.The compact radii (< 5 pc) of two of these clusters may be due to the reformation resulting in more compact NSCs; recent star formation is seen to be centrally concentrated in nearby NSCs including the Milky Way (Feldmeier-Krause et al. 2015), M31 (Lauer et al. 2012), and other nearby late-type NSCs (Carson et al. 2015).

Bulge and Bar Relations
Given that NSCs are located in the centers of galaxies, it is interesting to understand how they relate with the bulge and bar components in their hosts.Understanding this relationship might provide insights into their formation mechanism.In this section, we discuss in detail the correlations of the NSC mass with the fraction of light in the photometric bulge and the luminosity of the bar components hosted by our sample of galaxies.In Figure 9, we plot the NSC mass against the photometric bulge to total (B/T ) light ratio (top panel) and the bar luminosities (bottom panel).We define the photometric bulge to consist of all the components (including the NSCs) except the bar and the disk.We determine the luminosity of the bar components in the galaxy, integrating the flux of the bar and the boxy-peanut bulge components (we note that 8/26 of the galaxies with NSCs lack bar components and are not included).Also, since we do not model the disk in NGC 1300 and NGC 5248, we exclude these galaxies from the top panel in the figure.
We find no correlation between the NSC mass and the photometric B/T ratio in our sample.This indicates that the bulges in galaxies appear to be uncorrelated with the formation of the NSCs.Similarly, we do not see any correlation between the NSC mass and the bar luminosities.

NSC-Black Hole relations
NSCs and BH are found to co-exist in many massive galaxies (> 10 9 M ⊙ , e.g.Seth et al. 2008a;González Delgado et al. 2008).Early studies of scaling relations and the relative masses of these quantities suggested there may be a transition between NSCs dominating the nuclear mass in low mass galaxies and BHs dominating in higher mass galaxies (Ferrarese et al. 2006;Wehner & Harris 2006;Graham & Spitler 2009).The recent compilation of co-existing BH and NSC mass measurements in Neumayer et al. (2020) shows a clear trend where NSCs are typically more massive than BHs in lowermass galaxies, while the opposite is true in higher mass galaxies.NSCs become less common in massive early- type galaxies (> 10 10 M ⊙ ).This could be due to the dynamical impact of binary BH mergers (e.g.Milosavljević & Merritt 2001;Antonini et al. 2015).However, the trend in NSC/BH mass is not a simple one -This large scatter can be seen in the very different relative masses of the BHs and NSCs in the Milky Way (where the NSC is ∼10× the BH mass) and M31 (where the opposite is true).Unfortunately, a lack of NSC mass measurements in massive spiral galaxies has limited our ability to make this comparison more widely for MWlike galaxies.
From our galaxies with good NSC measurements, 4 galaxies have black hole mass measurements available in the Saglia et al. (2016) compilation (NGC 3368, NGC 3412, NGC 4501, and NGC 4699).In Figure 10 we show the ratio of BH and NSC mass of these galaxies against their host galaxy stellar mass (left panel) and the NSC mass (right panel) added to the data from Neumayer et al. (2020).The solid line in the figure represents equal NSC and BH mass.The objects above this line (including our measurement for the very massive latetype galaxy NGC 4699) have a more massive BH than NSC.Our measurement confirm that there are a wide range of BH to NSC mass ratios in massive galaxies.

CONCLUSIONS
In this paper, we have presented the photometric and morphological analyses of nuclear star clusters in 33 nearby galaxies from the Composite Bulges Survey.This includes a subsample of 20 MW-like galaxies with spiral morphologies (T = 1-4) and stellar mass from 10 10.4 M ⊙ to 10 11.1 M ⊙ .This MW-like subsample is a complete volume-limited sample of galaxies similar to the Milky Way that also meet the distance (< 20 Mpc), inclination (35 to 60 • ), and Galactic latitude (|b| > 20 • ) criteria of the complete CBS sample.
Using imfit, we obtain accurate models for the nuclear regions of the galaxies.We model the NSCs using Sérsic profiles in three HST filters and derive their sizes, colors, and masses.We present the Sérsic profile fit parameters of the NSCs in Table 2 and their derived properties in Table 3.
Our main results are: • We clearly identify NSCs in 78.8 +6.2 −7.9 % of our 33 galaxies, and 68.4 +9.5  −11.4 % for the MW-like subsample.NSCs may be present in other galaxies, but are missed due to dust or AGN, thus these nucleation fractions are lower limits.This work significantly expands the number of nucleated galaxies known in higher mass, late-type galaxies.The nucleation fractions are higher than, but consis- • We calculated the mass of our NSCs using color-M/L relations, and find a median mass of NSCs in our galaxies of log(M NSC /M ⊙ ) = 7.16.
• • We find a large scatter in NSC mass over a small range of galaxy mass, with two prominent lowmass outliers.These outliers also have small radii, suggesting a possible difference in NSC formation mechanism or evolution.
• Our NSCs are preferentially aligned with but are less flattened than their host galaxy disks.This alignment suggests these NSCs are forming either from gas accretion or star clusters inspiraling from the disk due to dynamical friction.
• Our NSCs do not show any correlation with the bar luminosity or photometric B/T ratios.
• We add four more galaxies to the small number of galaxies with known NSC and BH masses.One has a BH ∼ 10 times the mass of the NSC, while the others have NSCs that greatly exceed their BH masses.This confirms that massive galaxies have a wide range of NSC-to-BH mass ratios.2013,2018), DrizzlePac (Hack et al. 2013), Imfit (Erwin 2015, v1.8), grizli (Brammer 2019)

Figure 1 .
Figure 1.The stellar mass, distance, and Hubble type of the 33 NSC sample galaxies (circles) selected from the 53 galaxy CBS survey (triangles).The MW-like subsample of 20 galaxies are shown as squares.

Figure 2 .
Figure 2. 2D image modeling of an example galaxy, NGC 4689, using imfit.Each row describes fits on a different scale; the top row shows the largest scale fits to the Spitzer IRAC1 image, the middle row, the HST/WFC3-IR F160W image, while the bottom row shows the smallest scale and highest resolution fits to the HST/WFC3-UVIS F814W image.The left column shows the data, the next column rigth shows the best-fit model, and the right most columns show the residuals.In the bottom row, three residual images are shown, the left one is from a fit with no NSC, the middle with a point source NSC component, and the right with the best-fit model image shown; the stretch of all three of these residual images is identical and clearly shows the need for a resolved NSC component in this galaxy.

Figure 3 .
Figure 3. Top panel: An example 1-D surface brightness profile of NGC 4689 from our F814W imaging and modeling.Note this is the same galaxy shown in Fig. 2. The data is shown as green crosses, the best-fit model is shown as the blue solid line, and the best-fit NSC component is shown as the purple solid line.Lower panel: The residuals (model -data) in magnitudes.Three different residuals are shown corresponding to the right three panels in Fig. 2. The best-fit model with a resolved NSC, point source NSC model, and no NSC model are shown as blue, orange, and red lines respectively.

Figure 4 .
Figure 4. Top panel: CMD with the NSC UVIS color VS F814W band absolute magnitude.The circles denote the NSCs in our sample.The MW-like subsample is denoted by the squares.The open circles represent the dusty galaxies with a NSC quality fit = 2. Bottom panel: NSC color-color plot (UVIS vs. IR).Colored pluses and squares show PARSEC 1.2S SSP models(Bressan et al. 2012) at two metallicities (-1.0 and +0.3) and with ages from 1-13 Gyr (indicated by the color).An extinction vector of corresponding to AV = 1 is also shown.All plotted magnitudes in both panels are foreground extinction corrected.

Figure 6 .
Figure 6.Top panel: NSC mass vs. NSC Sérsic index.Bottom panel: NSC mass vs. NSC ellipticity.The solid circles denote our NSC sample galaxies, with squares indicating the MW-like subsample.All dusty galaxies (i.e., quality flag = 2) are denoted by open circles.The points are colored red or blue based on their Hubble type (Early-and Latetype respectively).

Figure 7 .
Figure 7. Left panel: Galaxy ellipticity vs NSC ellipticity.The population of NSCs is typically rounder than their host galaxy disks, suggesting a less flattened distribution.Circles denote the NSC sample with squares showing the MW-like subsample.The galaxies are colored red and blue based on their Hubble type (early and late type respectively).The NSCs with quality flag = 2, i.e., dusty centers are denoted by open circles.Right panel: Difference between NSC and Galaxy position angles(δPA).We plot the difference for all those NSCs whose ϵ > 0.05.There is a clear preference for near-alignment between NSCs and their host galaxies' disks.

Figure 8 .
Figure 8.The NSC mass is plotted against the host galaxy mass.The markings are the same as in Figure 5.The literature sample contains both early-and late-type galaxies from Côté et al. (2006); Turner et al. (2012); Georgiev et al. (2016); Sánchez-Janssen et al. (2019).

Figure 10 .
Figure 10.The ratio of black hole mass to NSC mass is plotted against the host galaxy mass (left panel) and against the NSC mass (right panel) based on plots from Neumayer et al. (2020) but adding in data with new NSC masses from our sample.The circles denote our NSC sample galaxies, with surrounding squares indicating the MW-like subsample.Open circles indicate uncertain measurements due to dust (i.e., quality flag = 2 in Table 3).The points are colored red or blue based on their Hubble type (Early-and Late-type respectively).Literature sample from Early-(red stars) and Late-(Blue stars) type galaxies are obtained from Seth et al. (2008a); Graham & Spitler (2009); Neumayer & Walcher (2012); Georgiev et al. (2016); Nguyen et al. (2018).Literature results with upper and lower limits for the BH and NSC masses are represented by the arrows.The solid horizontal line indicates equal NSC and BH masses in the galaxies.tent within the errors of the determination in Neumayer et al. (2020).
Our NSCs are consistent with the mass-radius relationship of literature NSCs (Figure.5).They also follow the galaxy stellar mass-NSC mass relation for late-type galaxies from Georgiev et al. (2016), significantly expanding the sample of NSCs at the high mass end (Figure.8).
Note-(1) Galaxy name, (2) NSC sample or MW-like subsample, (3) primary band used for modelling the different components in the galaxy, (4) total number of components we fit for in the galaxy (including the NSC component), (5) NSC component PointSource vs. Sérsic ∆AIC values (explained in detail in section 3.3), (6) quality of the NSC fits, see Section 3.4, (7) brief notes on the NSC fits.rationartifacts affect the r ≲ 0.3 ′′ central region in the optical images) that no clear NSC component is visible.Accordingly, we do not present NGC 1566 in our NSC parameters table.

Table 3 .
NSC Sérsic Parameters and derived properties.
that program.J.M.A. acknowledges the support of the Viera y Clavijo Senior program funded by ACIISI and ULL.This work is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.It has also made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Labora-tory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.We also made use of data from the Sloan Digital Sky Survey.Astropy (Astropy Collaboration et al. Software: