The Next Generation Virgo Cluster Survey (NGVS). XXVII. The Size and Structure of Globular Cluster Systems and Their Connection to Dark Matter Halos

We study the size and structure of globular cluster (GC) systems of 118 early-type galaxies from the NGVS, MATLAS, and ACSVCS surveys. Fitting Sérsic profiles, we investigate the relationship between effective radii of GC systems (R e,gc) and galaxy properties. GC systems are 2–4 times more extended than host galaxies across the entire stellar mass range of our sample (108.3 M ⊙ < M * < 1011.6 M ⊙). The relationship between R e,gc and galaxy stellar mass exhibits a characteristic “knee” at a stellar mass of M p ≃ 1010.8, similar to the galaxy R e –stellar mass relationship. We present a new characterization of the traditional blue and red GC color subpopulations, describing them with respect to host galaxy (g′−i′) color (Δgi): GCs with similar colors to their hosts have a “red” Δgi, and those significantly bluer GCs have a “blue” Δgi. The GC populations with red Δgi, even in dwarf galaxies, are twice as extended as the stars, suggesting that formation or survival mechanisms favor the outer regions. We find a tight correlation between R e,gc and the total number of GCs, with intrinsic scatter ≲0.1 dex spanning two and three orders of magnitude in size and number, respectively. This holds for both red and blue subpopulations, albeit with different slopes. Assuming that N GC,Total correlates with M 200, we find that the red GC systems have effective radii of roughly 1%–5% R 200, while the blue GC systems in massive galaxies can have sizes as large as ∼10% R 200. Environmental dependence on R e,gc is also found, with lower-density environments exhibiting more extended GC systems at fixed mass.


INTRODUCTION
The size of a galaxy, i.e., the spatial extent of its baryons, is the result of a number of important physical processes that govern its evolution: angular momentum (Mo et al. 1998), dissipation, and merging.Dissipational star formation, possibly the result of wet mergers or interactions, can reduce the size of a galaxy by driving gas to the center and preferentially forming stars there (Tacchella et al. 2016).By contrast, dissipationless interactions like dry mergers and tidal heating, can increase the sizes of galaxies without adding very much stellar mass (Naab et al. 2009).As a result the effort to measure the sizes of galaxies as a function of their stellar mass, halo mass, and other fundamental properties as a function of redshift, as well as reproduce those sizes in simulations (e.g., Furlong et al. 2017;Genel et al. 2018), has been central to our efforts to understand galaxy evolution.
While it may be possible to have a single number for a galaxy's "size" (usually the effective or half-light radius), studies of nearby galaxies have long shown that galaxies in fact have multiple distinct structural components (disks, bulges, bars, halos), each of which has its own characteristic size scale, and which reflects a different phase of evolution.Galaxy stellar halos are the most obvious example of how different structural components in a galaxy can have quite different spatial distributions and evolutionary histories.Stellar halos are the most extended and diffuse structures in galaxies, and yet also contain a stellar population that represents some of the most intense star formation in a galaxy's evolution: globular clusters.This apparent contradiction, as well as many other properties, makes globular clusters, and their apparently large spatial extent relative to their host galaxies, an intriguing target for study.
Globular clusters (GCs) -compact, dense associations of stars that surround galaxies spanning wide ranges in mass, type and environment -are valuable probes of galaxy formation processes.Similar to the stars that make up the bulk of their host galaxies, the ubiquity of GCs requires any complete theory of galaxy formation to explain not just the existence of GC systems but also their detailed properties, including the observed trends with host galaxies properties.Compared to the stars that make up galaxies, GCs offer some distinct advantages as galaxy formation probes.From a purely observational perspective, their compact nature (⟨r h ⟩ ≃ 3 pc) and high luminosities (⟨L⟩ ≃ 10 5 L ⊙ ) allow them to be readily recognized against the diffuse stellar light of the underlying galaxy.As Nature's closest approximations to simple stellar populations (SSPs), it is possible from imaging and spectroscopy to deduce ages and abundances for individual clusters, and thereby investigate the past history of star formation and chemical enrichment in their host galaxies (e.g., Brodie & Strader 2006, and references therein).If radial velocity measurements are also available, then GCs can be used as dynamical probes for studying the amount and spatial distribution of gravitating mass within galaxies -from their centers to their outermost regions, where dark matter dominates (e.g., Alabi et al. 2017;Toloba et al. 2018;Fensch et al. 2019;Müller et al. 2020).
Over the past two decades, our understanding of the detailed properties of GC systems increased dramatically, thanks mainly improvements in observing facilities.For example, the deployment of the ACS and WFC3 cameras on the Hubble Space Telescope (HST) made it possible to obtain high-resolution imaging across the ultraviolet-optical-infrared region for GCs belonging to many nearby galaxies, including both dwarfs and giants (e.g., Côté et al. 2004;Jordán et al. 2007).On the ground, wide-field imaging performed with 2.5 to 4m-class telescopes, and multi-object spectroscopy with 8-10-m-class telescopes, have also led to significant improvements in our understanding of the properties of GC systems, especially those associated with high-and intermediate-luminosity galaxies: i.e., their overall numbers, formation efficiencies, luminosity and mass functions, color distributions, color-magnitude relations, and dynamical properties (see, e.g., the review of Brodie & Strader 2006).
On the other hand, our knowledge of spatial distributions for GC systems remains far more limited.This issue has taken on a renewed importance in recent years, for several reasons.First, cosmological simulations have reached a level of sophistication where it is now possible to make quantitative predictions for the spatial distributions of GC systems and their dependence on galaxy mass and assembly history (Kruijssen 2015;Kruijssen et al. 2019Kruijssen et al. , 2020;;Reina-Campos et al. 2019, 2022;Chen & Gnedin 2022).Such predictions must be carefully tested against observations for GCs in real galaxies, ideally selected to span wide ranges in mass, type and environment.Meanwhile, an increasing number of studies have used data from the literature to explore the connection between the total numbers of GCs in galaxies and the underlying mass in baryons and dark matter -with the latter quantity usually inferred from an assumed stellarto-halo mass relation (Peng et al. 2008;Hudson et al. 2014;Harris et al. 2015Harris et al. , 2017;;Forbes et al. 2016;Hudson & Robison 2018;Forbes et al. 2018).These studies suggest that GC systems contribute a nearly constant fraction of the total dark matter mass, supporting earlier claims based on observations of GCs in massive cluster galaxies (e.g., Blakeslee et al. 1997;Blakeslee 1999;McLaughlin 1999).
From a practical perspective, measurements of the size and the richness of GC systems are linked in a fundamental way: i.e., a reliable determination of the total number of GCs associated with a galaxy generally requires a measurement of its GC density profile, which can be integrated to give the total number of GCs.In principle, this is straightforward, but our knowledge of GC density distributions, and their variation along the mass spectrum of galaxies, remains surprisingly limited.Most investigations into the structure of GC systems either combined photographic and CCD data for a few bright galaxies (Harris 1986;McLaughlin et al. 1994), performed GC counts using mosaic CCD cameras on ground-based telescopes for small samples of intermediate-mass galaxies based on GC counts (e.g.Rhode & Zepf 2001, 2004;Rhode et al. 2007;Hargis & Rhode 2012, 2014;Ko et al. 2019).Several studies have examined the spatial distribution of GCs using larger sample sizes.(Kartha et al. 2014) proposed a linear relationship between the effective radius of galaxies and the effective radius of GC systems based on six galaxies from the SLUGGS survey (Brodie et al. 2014).Zaritsky et al. (2015) studied the GC properties of 97 early-type galaxies using S 4 G survey images (Sheth et al. 2010), and they estimated the total number of GCs by fitting the GC radial number densities with a power-law function with a fixed power value.Hudson & Robison (2018) and Forbes (2017) studied GC spatial distributions using 20 − −30 galaxies, and suggested the linear relations between galaxy stellar mass and GC system effective radii, but slopes of the relations in two studies are significantly different.Caso et al. (2019) andDe Bórtoli et al. (2022) showed different slopes for the relationships between galaxy stellar mass and effective radius of GC systems, distinguishing between high and low mass galaxies.A power-law relationship between the effective radius of GC systems and the total number of GCs has been suggested.However, these studies are limited by the small number of galaxies considered; or by the use of GC counts from narrow-field HST images, which often require substantial, sometimes uncertain, corrections to account for missed GCs at larger radii; or by the use of more assumptions about the GC radial distribution, which provide limited information about the GC spatial distribution.
In this paper, we build on this latter approach by presenting a comprehensive study of the structure and sizes of the GC systems associated with more than one hundred galaxies in the local Universe.Our analysis relies on imaging for 87 early-type galaxies in the Virgo Cluster acquired as part of the Next Generation Virgo Cluster Survey (NGVS, Ferrarese et al. 2012Ferrarese et al. , 2020) ) -a multi-band imaging survey undertaken with the Canada-France-Hawaii Telescope (CFHT).The depth, areal coverage and multi-band nature of the NGVS 1 makes it an ideal resource for studying GCs in nearby galaxies.Previous papers in the NGVS series have examined: the distribution of GCs within galaxies and across the entirety of the cluster (Muñoz et al. 2014;Durrell et al. 2014;Sánchez-Janssen et al. 2019;Lim et al. 2020); the stellar populations within GCs (Powalka et al. 2016a(Powalka et al. ,b, 2017(Powalka et al. , 2018;;Ko et al. 2021); and galaxy kinematics and dynamics across a range of spatial scales using GC radial velocities (Zhu et al. 2014;Zhang et al. 2015;Toloba et al. 2016Toloba et al. , 2018;;Longobardi et al. 2018;Li et al. 2020;Taylor et al. 2021).For 75 of the galaxies in our sample, HST imaging is also available from the ACS Virgo Cluster Survey2 (Côté et al. 2004), which allows us to accurately measure GC density profiles into the core region.To supplement the NGVS sample, we include 31 additional galaxies from the Mass Assembly of early-Type GaLAxies with their fine Structures (MAT-LAS, Duc et al. 2015;Duc 2020) survey.Like the NGVS, the multi-band MATLAS3 survey was carried out with CFHT and has the benefit of depth, wide field coverage and sub-arcsecond image quality.
This paper is structured as follows.In §2, we describe the data analyzed in this paper and the methodology used to study the spatial distribution of GCs in our sample galaxies.In §3, we present our determinations of the effective radii of our GC systems and show how the measured sizes correlate with other galaxy parameters, including stellar masses, stellar radii and halo masses.In §4, we discuss the implications of our results within the context of current models of galaxy formation.In §5, we summarize our findings and outline some directions for future work.
Our analysis is based on data acquired in two complementary imaging surveys undertaken with MegaCam on the Canada-France-Hawaii Telescope and supplemented by imaging from the ACS/WFC instrument on HST.Table 1 summarizes the data used in this paper, including survey names and filters, numbers of galaxies, and the average values and full ranges spanned by our respective samples in magnitude, color and stellar mass.

Surveys and Sample Definition
The primary survey used in this analysis is the NGVS, a u * g ′ i ′ z ′ imaging survey of the Virgo cluster carried out with Megacam on CFHT (Boulade et al. 2003).The survey covers an area of 104 deg 2 -from the cluster core out to the virial radii of its two main sub-clusters, Virgo A and B. The survey is comprised of 117 distinct pointings, with a uniform (point-source) limiting magnitude of 25.9 mag, and a median seeing of 0. ′′ 80 in the g-band.Full details on the survey, including observing and data processing strategies, are available in Ferrarese et al. (2012Ferrarese et al. ( , 2020)).
We also rely on data from MATLAS (Duc et al. 2015;Bilek et al. 2022), a deep imaging survey of ATLAS 3D galaxies (Cappellari et al. 2011) that is also based on imaging from CFHT Megacam.The ATLAS 3D survey targeted 260 massive early-type galaxies with distances d < 42 Mpc and magnitudes M K < −21.5 mag (Cappellari et al. 2011).There are 58 ATLAS 3D galaxies within the NGVS footprint, and the MATLAS survey uses NGVS data for these objects.Outside of the NGVS region, the MATLAS survey was designed to be carried out with the u * g ′ r ′ i ′ filters but, in practice, most galaxies were observed in a subset of these filters (g ′ and r ′ being most common).For this study, we consider only the subset of MATLAS galaxies with MegaCam photometry in at least three passbands, usually g ′ r ′ i ′ or u * g ′ r ′ .The seeing varies with targets and filters but, on average, image quality in the g-band is ≃ 0. ′′ 84.The details for dithering strategies and stacking images are described in Duc et al. (2015).The MegaPipe data reduction and processing pipeline (Gwyn 2008) provides the astrometric and photometric calibration of GC candidates.The surveys have a similar limiting surface brightness of µ g ≃ 29.0 mag arcsec −2 , but limiting magnitudes can vary significantly depending on the seeing for specific targets.
Among the NGVS and MATLAS samples, we have chosen galaxies for our analysis in the following way.First, we selected all MATLAS galaxies (including the 58 in the NGVS region) observed in three or more filters.Because GCs at distances of ∼ 10-45 Mpc appear mostly as point sources in ground-based images, we dis- tinguish GCs from foreground stars using their colors.Although we could select GC candidates using magnitudes and a single color index, we have used color-color diagrams to decrease contamination by foreground stars and compact, background galaxies.Second, we limited our targets to distances of 25 Mpc or less in order to detect enough GCs for an analysis of their spatial distribution.Finally, we included a sample of dwarf galaxies in the NGVS that were also observed in ACS Virgo Cluster Survey (Côté et al. 2004), where images taken with the ACS instrument (Ford et al. 1998) in WFC mode can be used to select GCs with high confidence (see, e.g., Jordán et al. 2004Jordán et al. , 2009)).We have 139 galaxies with the above criteria as a parent sample.As described in §2.3, not all of these galaxies yielded reliable GC radial number density profile fits, and so from a parent sample of 139 galaxies, the final sample used in our analysis consists of 118 galaxies.Of these, 87 galaxies are from the NGVS, with 75 of them having ACSVCS data.Additionally, 31 galaxies are from the MATLAS.Given that the NGVS and MATLAS surveys have data in common, to reduce confusion we refer to "MATLAS galaxies" as those MATLAS galaxies outside of Virgo, i.e., not in the NGVS footprint.Note-a ACSVCS data is supplementary for the NGVS, so these galaxies are included in the NGVS.Full list of galaxies is available in Table A.1

Sample characteristics
Figure 2. The distribution of distances for galaxies in this study.Red and blue histograms show results for the NGVS and MATLAS samples, respectively.The black solid histgram shows results for total samples in this study.The MATLAS galaxies in this study have distances d ≲ 25 Mpc.For ACSVCS galaxies, distances are based on HST surface brightness fluctuation measurements (Mei et al. 2007;Blakeslee et al. 2009;Cantiello et al. 2018), while distances for MATLAS and NGVS (not in ACSVCS) galaxies are taken from Cappellari et al. (2011).These are stacked histograms, so the top edge line shows the distribution of distances for the entire sample.
Figure 1 presents a color-magnitude diagram for our program galaxies.All magnitudes and colors of galaxies are estimated on NGVS and MATLAS images.As expected, our targets roughly follow the red sequence of galaxies in the central region of Virgo cluster (Roediger et al. 2017).The galaxies in the MATLAS survey also lie on the red sequence but seem to have a little larger scatter than those in the NGVS survey.Figure 2 are the distribution of galaxy distances.All galaxies are more distant than 10 Mpc, with the great majority of our targets belonging to the Virgo Cluster at d ≃ 16.5 Mpc (Mei et al. 2007).Our sample includes a handful of galaxies beyond 20 Mpc -some from the MATLAS program, and some from the NGVS, which includes a few members of the W ′ Cloud at d ∼ 23 Mpc (Mei et al. 2007;Blakeslee et al. 2009).
Stellar masses for our program galaxies were calculated by modelling of their spectral energy distributions (SEDs) in all available bands (i.e., as many of u * g ′ r ′ i ′ z ′ as were available).We used the Flexible Stellar Population Synthesis models of Conroy et al. (2009) and assumed exponentially declining star formation histories and a Chabrier initial mass function.The SEDs were fitted to a grid of 50,000 synthetic models with metallicities in the 0.01 ≤ Z/Z ⊙ ≤ 1.6 range, star formation timescales 0.5 ≤ τ ≤ 100 Gyr −1 , and luminosityweighted ages between 5 and 13 Gyr.The mass errors are derived from the 16th and 84th percentiles of the marginalized posterior for the mass: we obtain samples of the posterior from our MCMC analysis, find the 16th and 84th percentiles, and set the error to 0.5× (P84 -P16).Mean mass error of our sample is about 0.04 in log with maximum and minimum of 0.10 and 0.01, respectively.All details of stellar mass estimation are in Roediger et al. (2023, in prep.). Figure 3 shows the distribution of these stellar masses.For comparison, we also show in this figure stellar mass distributions for galaxies targeted in two recent studies of GC systems (Alabi et al. 2017;Hudson & Robison 2018).Compared to these previous studies, our program galaxies span a significantly wider range in stellar mass (a factor of roughly ∼2000 with the faintest dwarfs in our program having stellar masses of M * ≃ 2 × 10 8 M ⊙ .

Photometry, Source Catalogs and Globular Cluster Selection
Our methods for source detection and photometry are similar to those used in Durrell et al. (2014) and Liu et al. (2015).Briefly, we used Source Extractor (Bertin & Arnouts 1996) in dual-image mode to detect sources and perform photometric measurements.For NGVS target galaxies, we used the g ′ -band image as a detection image.Fluxes were measured within circular apertures for a variety of radii, and magnitudes corrected to 16-pixel (3 ′′ ) diameter aperture magnitudes.These magnitudes were then transformed onto the standard AB system by using PSF magnitudes of bright SDSS stars (Albareti et al. 2017) within the field.A similar methodology was adopted for the MATLAS galaxies, although in this case, r ′ -band images were used for source detection when the g ′ -band images suffered from significantly poorer seeing (≳ 1. ′′ 5).
We subtracted the diffuse component of each galaxy from the original image to ensure reliable photometry for GCs located in the central regions of galaxies.For NGVS galaxies, we used customized two-dimensional models for diffuse light subtraction, as described in Ferrarese et al. (2020).Galaxy models were fitted and subtracted independently in each band.For the MATLAS galaxies, we used ring median models generated with the RMEDIAN task in IRAF (Tody 1986(Tody , 1993) ) to subtract the underlying galaxy light.Galaxy subtraction aside, all photometric measurements were carried out using the same methodology as for the NGVS targets.
Completeness tests were carried out using artificial stars generated with DAOPHOT (Stetson 1987).We in-jected 5000 artificial stars into model subtracted images in each observing field, and derived the recovery rate as a function of both magnitude and underlying surface brightness.Most regions have recovery rates in excess of 90% for g ≲ 24.5 mag, apart from the expected sharp decrease in the core regions of high surface brightnesses galaxies.
GC candidates were then selected on the basis of their sizes (i.e., compactness) and colors.Most GCs appear as point sources in our survey, so we selected point sources using a measure of how extended an object is relative to a point source (i.e., "inverse concentrations"); these indices are the magnitude differences between 4-pixel (0. ′′ 75) and 8-pixel (1.′′ 5) diameter aperture magnitudes, ∆m 4−8 , normalized so that point sources have a mean value of zero (Durrell et al. 2014).We chose sources with −0.08 < ∆m 4−8 < 0.08 as point sources.This index is actually measured in two filters (g ′ and i ′ for NGVS and top two best seeing filters for MATLAS), where the final ∆m 4−8 is the combination in quadrature of the measurement of each object in the two filters, i.e., ∆m 4−8 = (∆g ′ 4−8 ) 2 + (∆i ′ 4−8 ) 2 for NGVS.Because GCs in nearby galaxies (d ≲ 20 Mpc) can be slightly extended in images taken during conditions of good ground-based seeing (≲ 0.8), we used a slightly expanded range, (−0.08 < ∆m 4−8 < 0.16, Figure 4).We also required GC candidates to have m g < 24.5 mag because ∆m 4−8 values for point sources and extended sources begin to merge below g ′ ∼ 24.5 mag.
From these "point-like" sources, we selected GC candidates based on their colors.Using spectroscopically confirmed GCs in M87 as a guide, u * g ′ i ′ or g ′ r ′ i ′ colorcolor diagrams were used to select likely GCs (see Figure 5).For M87, the polygons shown in the u * g ′ i ′ and g ′ r ′ i ′ color-color diagrams include all spectroscopically confirmed GCs.
Finally, we used ACSVCS photometry (Côté et al. 2004) in the central regions to supplement the NGVS data for galaxies overlapped with ACSVCS galaxies in our NGVS sample.GC catalogs were taken from Jordán et al. (2009), and transformed from the HST to CFHT MegaCam systems using relations fitted with photometry of GCs in M87.Followings are the fitted relations: GC catalogs generated from these deep, high-resolution images are crucial for measuring the density profiles of GC systems in galaxy cores. .Density map for GC candidates in the vicinity of two of our program galaxies, NGC4649 and NGC4621.The image measures 100 ′ × 100 ′ , with a pixel size of 1 ′ × 1 ′ , and has been smoothed using a 2-pixel Gaussian filter.North is up and East is to the left.Inner and outer circles are drawn at radii of Re,gc and 3Re,gc centered on each galaxy.The two GC systems have been fitted simultaneously in our analysis, while two other galaxies that are not members of our sample (NGC4660 and NGC4647) have been masked.

Globular Cluster Density Profiles
Our goal is to measure radial number density profiles for GC candidates in our program galaxies and use these profiles to explore the spatial distribution of the GC systems as a function of galaxy properties.An example of our methodology is presented in Figure 6 which shows a Gaussian-smoothed density map for GC candidates in the vicinity of NGC4649 and NGC4621 -two program galaxies that happen to fall within this single 100 ′ ×100 ′ field.As described below, we mask interloping galaxies and extract a radial number density profile for GC candidates associated with each program galaxy.Note that, unlike the case shown in Figure 6, the vast majority of our program galaxies are relatively isolated and can be fitted individually.
Figure 7 shows a radial number density profile of GC candidates in NGC4564 (VCC1664), a typical galaxy for which both ACSVCS and NGVS imaging is available.For display purposes, we plot number densities computed in circular, logarithmically-spaced bins, although we rely on the locations of individual GC candidates in the profile fitting.The density profiles are then parameterized with a modified Sérsic function where Here θ is the position angle measured in the customary sense: i.e., from North to East.The fitted Sérsic function provides us with an estimate for both the size (R e ) and concentration (n) of the GC system.Because the GC catalogs for each galaxy inevitably include some residual contamination (i.e., from foreground stars, background galaxies, and interloping GCs from nearby galaxies), we include a constant background term, Σ b , in our fits.This function was fitted to the data with a Markov Chain Monte Carlo (MCMC) method using the emcee code in python (Foreman-Mackey et al. 2013).Flat priors were set on several parameters: Σ e > 0.001 arcmin −2 , 0.05 < n < 8.0, and 0.05 < R e < 30 arcmin.For the ellipticity and position angle, we imposed flat priors of 0 ≤ ϵ < 0.1 and −10 • < θ < 10 • due to difficulties in determining these values in the presence of GC contamination.For a handful of clearly elongated galaxies, we adopted a flat prior of 0 ≤ ϵ < 0.4 and no prior for θ.
For the background level, we adopted a Gaussian prior with mean and standard deviation estimated from the data.The mean of the Gaussian was estimated using the mean GC number density in a large annulus around the galaxy.For most fields, we used a ring with an inner radius of 30 ′ and a width of 1 ′ .For the largest galaxies, an inner ring radius of 60 ′ was used.The standard deviation of the Gaussian was estimated in the same annulus by dividing it into nine angular sub-sections (i.e., "annular arcs"), and calculating the standard deviation of the mean backgrounds in these nine sub-annular regions.Using these values, we could adopt a reasonable Gaussian prior on the background level.
We employed a likelihood function as follows: where is the probability of finding the datum i at radius R i , given the three Sérsic parameters, background level, and estimated completeness.To calculate the probability, we employed an areal integration of the modified Sérsic function.While the standard Sérsic function has a theoretical integrated functional form, this is not true of our modified Sérsic function, so we used a numerical integration method to generate the likelihood function.The integration range was normally from the radius corresponding to the 50% completeness limit to 30 ′ .For NGVS galaxies, where ACSVCS imaging is available, the integration was started at the galaxy center.As with the largest galaxies, the outer limit was taken to be 60 ′ .Figure 8 shows representative corner plots from MCMC fitting of the GC system in NGC4564.
Aside from the familiar correlation between R e and Σ e , there are generally no strong correlations between any of the other parameters except for Σ e and Σ b .In a few cases, there are also weak correlations between Sérsic index, n, and effective radius, R e .In general, the fitted models provide good matches to the measured profiles, as shown in Figure 7.
As a final check on our results, we visually inspected the fitted profiles for each our program galaxies and excluded galaxies when either: (1) The GC system is strongly affected by that of a neighbouring galaxy.For example, the GC system of NGC4486A is difficult to separate from that of M87.
(2) The fitted parameters were unreliable due to a small numbers of GCs.Specifically, results were rejected when the error on the fitted effective radius exceeded 100%.
Of our parent sample of 139 galaxies, we find 118 GC systems to have well measured GC density profiles with reliable model fits.

Analysis of Globular Cluster Colors
Up to this point, the analysis has been carried out for the full GC system for each galaxy: i.e., with no restriction on GC colors.However, there is considerable interest in understanding how the size and structure of GC systems depends on both age and metallicity, so the fitting process was repeated after dividing each GC system into blue and red sub-components on the basis of their (g ′ − i ′ ) 0 colors.
To divide the sample, we fitted each GC color distribution using a Gaussian Mixture Modelling (GMM, Muratov & Gnedin 2010) code that identifies the dividing point in color between the blue and red GC systems.If the D value returned by GMM (a statistic that measures the separation of the means relative to their widths) was found to be smaller than 2.0, including errors, then we considered the GC system to be unimodal based on the manual of the GMM code.A more detailed study of the GC color distributions will be presented in future paper; for this analysis, we used these GC subsamples to explore of the sizes of the chemically distinct GC systems, fitting the separate radial number density profiles with Sérsic profiles in the same way as for the Lim et al.
full GC systems.Altogether we find 51 galaxies have statistically-significant bimodal GC color distributions.
Although separating the blue and red GC subsystems is relatively straightforward in high-mass galaxies, it can be problematic in intermediate-and low-mass systems -a regime in which some galaxies have unimodal GC color distributions and, even for bimodal systems, the separation between the GC peaks can become subtle (Peng et al. 2006).Moreover, the interpretation of these GC colors becomes problematic across a wide range in galaxy mass.In high mass galaxies, the often-seen blue and red GC sub-populations are respectively interpreted to be accreted and in situ populations of GCs.However, in dwarf galaxies, their often unimodal populations are almost always blue in color because of their low metallicity (e.g., Larsen et al. 2022).While this provides a link between the accreted population of GCs in more massive galaxies and their likely low-mass, metal-poor progenitors (e.g., Côté et al. 1998), it obscures the origins of the GCs in dwarf galaxies themselves.Are the blue GCs in low-mass galaxies accreted as well?This is unlikely given the merger rate for dwarf galaxies, and the inefficiency of star cluster formation in the halos that accrete onto them.Given the similarity in color between dwarfs and their GCs, it is likely that the star formation episodes that formed the GCs also formed a significant fraction of the galaxy, and that they are thus associated.
For this reason, we introduce a new parameter to describe the color of GC sub-populations, ∆ gi , that measures the difference between the color of the host galaxy and that of a (sub-)population of GCs.In the case of bimodal GC color distributions, a galaxy can have two ∆ gi values, one for each sub-population peak color, and for unimodal GC populations, it is simply the difference between the galaxy color and the median GC color.
Figure 9 illustrates the definition of this parameter for two program galaxies.For galaxies with bimodal GC color distributions (as is the case for NGC4621, shown in the left panel), we denote the blue and red GC peak colors by B and C. Relative to the underlying galaxy, which has a color4 A, the GC subpopulations have values of ∆ gi,BGC = (B − A) for the blue GCs, and ∆ gi,RGC = (C − A) for the red GCs.For galaxies with unimodal GC color distributions (such as NGC4612, shown on the right), we define a single parameter, ∆ gi = (B − A), as the difference between the median GC and galaxy color.We will return to the use of this parameter when discuss the color dependence of GC scaling relations.

RESULTS
Our program galaxies span a factor of ∼2000 in stellar mass and, as we show below, an equivalently wide range in GC system size, with R e,GC ∼ 1 kpc for dwarfs and ∼ 100 kpc for the brightest giants.In general, the blue GC systems have effective radii similar to those for the full GC systems while the largest red GC systems have effective radii of ∼ 30 kpc -much smaller than the blue or composite GC systems.Although the scatter is large, the measured Sérsic indices mostly fall in the range 1 ≲ n ≲ 4 with a median value of n ∼ 2 (albeit with fairly large measurement errors in most cases).The GC systems of some massive galaxies have large values of n ∼ 4, which is comparable to the galaxies themselves (Ferrarese et al. 2006).Errors on the measured effective radii are generally much smaller, with uncertainties of ∼ 25% being typical.R e,gc will thus be our primary tracer of the spatial extent of the GC systems.In this section, we examine how GC system effective radius scales with the properties of the host galaxy and dark matter halo.Since stellar mass is a fundamental parameter that dictates many other galaxy properties, we begin our analysis by considering the the dependency of GC system effective radius on the stellar mass of the host galaxy.The measurements that form the basis of our analysis are recorded in Table 2 of Lim et al. (2024 in preparation).
Figure 10 shows the effective radius of the full GC system plotted against the stellar mass of the host galaxy.
Results from this study are plotted as black symbols.For comparison, the gray curve shows the best-fit polynomial for galaxies in the Virgo core region from Côté et al. (2023, in prep.).This relation exhibits the familiar steepening of the galaxy size-mass relation above ∼10 10.5 M ⊙ .We can see from Figure 10 that this behavior is mirrored by the GC systems of galaxies in this same environment (Virgo), with an unmistakable "break" in the size-mass relation at intermediate masses.In this log-log representation, the trend is clearly nonlinear; to capture this behaviour, we fitted the data using a broken power-law function (Mowla et al. 2019) that has been used previously to analyze the size-mass scaling relation of galaxies: Here R p is the "pivot" radius of the GC system which marks the change in slope, M p is the pivot stellar mass, α is the slope at the low-mass end, β is the slope at the high-mass end, and δ is a smoothing factor.After some experimentation, we opted to fix the smoothing parameter to δ = 6 rather than allowing it to vary as a free-parameter (see also Mowla et al. 2019).
Our fitting results are summarized in Table 2.The slope we measure at the high-mass end, β = 1.30+0.11  −0.17 , is consistent with values found in previous studies: 1.30 ± Figure 11.The effective radii of blue and red GC subsystems plotted against the stellar masses of their host galaxies.The data points show the effective radii for blue and red GC systems, or for the total GC system of galaxies with unimodal GC color histograms.Data points are color-coded by the difference in color between the GC system and host galaxy as shown by the color bar at the right.GC systems with well fitted effective radii are only shown.The blue dashed curve shows the best-fit broken power-law relation for blue GCs.The red dotted curve shows the best-fit result for red GCs.In both cases, the light curves show 500 random MCMC samples of based on these best-fit results.The gray curve shows the relation between effective radius and stellar mass for NGVS galaxies (Côté et al. 2023, in prep.).
0.14 (Hudson & Robison 2018) and 0.97 ± 0.4 (Forbes 2017).The slope at the low-mass end -a regime that is largely unexplored -is significantly shallower than this, with α = 0.34 +0.04 −0.04 .We find the break in the size-mass relation to occur at a stellar mass of M p = 6.5 +1.9 −1.5 × 10 10 M ⊙ .Caso et al. (2019) reported a similar estimate for this transition (i.e., ∼ 4 × 10 10 M ⊙ ) although the functional form they adopted does not fit the data well at low and high masses (see §4.2).
Figure 11 shows these same relations for the blue and red GC systems.In this plot, symbols for individual GC systems -either their separate blue and red components when bimodal, or the total color distribution, if unimodal -have been color-coded by the ∆ gi parameter.This parameter, which has been discussed in §2.5, measures the color offset relative to the underlying galaxy.In this plot, we see that the size-mass relation separates into two branches at high masses, having distinct "blue" and "red" ∆ gi values.Below stellar masses of ∼ 10 10 M ⊙ , where unimodal or weakly bimodal GC systems are common, ∆ gi values are intermediate in color, with many "green" values.In this regime, the GC systems appear to define a single, shallow relation.Best-fit broken power-law parameters for the blue and red GC systems are recorded in Table 2.These were estimated using a weighted fit, where each point was weighted by both its uncertainty and its ∆ gi parameter.Therefore, there is no strict cut that includes or excludes data points from each fit.For example, the power law fit for the blue GC systems includes all points, but gives greater weight to the ones that have bluer ∆ gi .For galaxies with bimodal GC color distributions, this effectively gives the red modes nearly zero weight, whereas the unimodal GC color distributions are often given an intermediate weight.
As a rule, nearly all blue GC systems are larger in size than the red GC systems within the same host galaxy.The trends for the blue and red GC systems show differences from that of the composite GC systems.Although the slope of blue GC systems at high-mass, and the break point of red GC systems, are consistent with those of the composite GC systems, the slope of low-mass blue GC systems, the break point of blue GC systems, and both slopes of the red GC systems are different with those of composite GC systems.Interestingly, the slope at low-mass for composite GC systems is consistent with a mean value of the slopes at low-mass for blue and red GC systems.The above results suggest that the break observed in composite GC systems is not only driven by a transition to having more blue GCs at higher mass, but also by an increase in size of red GC systems (mirroring the stars).We will discuss the implications of this finding in §4.

Re,gc vs. Galaxy Size
We now turn our attention to the variation in the relative sizes of the GC systems and underlying galaxies.Figure 12 shows a comparison between the effective radii of the GC systems, R e,gc , and those of their host galaxies, R e, * .In the left panel, the black symbols show results based on the size of the full GC system; on the right, we show results for the blue and red GC systems separately, with the symbols color-coded as in the right panel of , the GC systems are observed to be significantly more extended the the galaxies themselves -conclusion that holds whether one considers the composite GC system, or separate the blue and red subcomponents.We also note that the trend for red GCs appears to show less scatter than that for the blue GCs (rms = 0.03 vs. 0.06).A comparison of Figure 10 and 12 shows the correlation between GC size and stellar mass to be significantly tighter than the correlation with host galaxy size.Of course, much of the scatter in Figure 12 may be attributable to measurement errors in R e, * , which can be significant (see, e.g., Chen et al. 2010).Finally, we point out that the green dashed line in Figure 12 shows the relation R e,gc = 1.5R e, * , which has sometimes been adopted as an estimator of GC system size, most notably in the study of ultra-diffuse galaxies (e.g., van Dokkum et al. 2016).As this figure shows, this relation significantly underestimates the size of the GC system.If we fix the slope of fitted relation at unity in log-log spaces, we obtain a best-fit scaling relation of R e,gc = (4.5 ± 0.1)R e, * for all GCs.If we consider the blue and red GCs separately, then we find factors of 7.9 ± 0.1 and 2.2 ± 0.1, respectively.
Here we show that even the red GC population-i.e., the ones with colors most closely matched to those of the host galaxy population and are most likely to have formed in situ-have an extent that is twice as large as the underlying stars (on average), and that this is true over a wide range of mass.Given that intense star formation episodes in which GCs are most likely to form are often centrally concentrated, this result requires that even GCs that form in situ encounter evolutionary mechanisms (whether formation or destruction) that bias their spatial distribution toward a galaxy's outer regions.We note that we all of our measurements are for roughly the brighter half of the GC luminosity function, and that if there are mass-dependent processes (such as dynamical friction), then deeper studies might reveal a different radial distribution for fainter GCs.

Scaling Relations: II. Dark Matter Parameters
Having established in the previous section that the sizes of GC systems -including those of their blue and red subcomponents -are correlated with the stellar masses of their host galaxies, we now turn our attention to possible connections with the total numbers of GCs and, by extension, the masses of their host dark matter halos.Our discussion of these parameters is linked by the assumption that the total number of GCs belonging to a galaxy is directly proportional to the total mass of the galaxy, which is in turn dominated by the dark matter halo.This possibility, which dates back to Blakeslee et al. (1997) and McLaughlin (1999), has gained momentum in recent years, with several studies providing support for the claim (e.g., Blakeslee 1999;Hudson et al. 2014;Harris et al. 2015;Forbes et al. 2016;Harris et al. 2017).
Figure 13 shows the dependence of the GC system effective radius on the total number of GCs in each galaxy, N GC,T otal , using data from Table 2 of Lim et al. (2023 in prep.).We emphasize that both of these parameters are measured directly from our GC density profile fits: i.e., the effective radius for each GC system is measured as a free parameter, while the total number of GCs is calculated by numerically integrating the backgroundcorrected density profile (Equation 1).Following this integration, we correct the GC numbers for our magnitude limit in each galaxy using the dependence of the GC luminosity on host galaxy luminosity (Jordán et al. 2007;Villegas et al. 2010).Figure 13 presents the result of this exercise in the same way as Figure 12: i.e., the left panel shows results for the full GC system while results for the individual blue and red GC systems are shown on the right.
It is apparent that galaxies show a tight correlation between the effective radii of their GC systems and the total number of GCs.This is not surprising given that the two parameters are related to each other via the Sérsic law used in our analysis (although the concentration of the GC system, as measured by Sérsic index, and the density at the effective radius also contribute to the total system richness).However, we stress that, in Figure 13, the two-dimensional error ellipses for individual points have randomly distributed ellipticities and position angles.This indicates that correlated errors are certainly not responsible for the observed trend, which spans two full decades in effective radius and more than three decades in N GC,T otal .
We fitted the relations in Figure 13 with a linear function in log-log space: log Re,gc = (0.65 ± 0.02) log N GC,T otal − (0.77 ± 0.06) all (12) = (0.69 ± 0.042 log N GC,T otal − (0.92 ± 0.06) blue (13) = (0.47 ± 0.03) log N GC,T otal − (0.48 ± 0.08) red( 14) The tightness of this relation was previously noted by Caso et al. (2019).However, they found the trend to flatten at low galaxy masses (corresponding to N GC,T otal ≲ 100) and opted to fit the relation using a second order polynomial.In this study, we find that a single linear relation (shown as the dashed line in the left panel) provides an excellent representation of the data from N GC,T otal ∼ 10 to ∼ 10000.
The right panel of this figure shows this same relation but now replaced with the effective radii of the blue and red GC subcomponents.At the high-mass end, there is a clear separation between the blue and red GC systems.However, at low masses, the separate branches disappear entirely and only a single relation is apparent below N GC,T otal ∼ 300, with the best-fit relations formally crossing at N GC,T otal ∼ 100.In other words, the blue and red GC systems belonging to these low-mass galaxies no longer show distinct spatial distributions.We will return to this issue in §4, since it may support the view that blue GCs in these low-mass hosts represent bonafide "in-situ" populations from which the GC systems of higher-mass galaxies have been assembled.
Finally, the horizontal axes along the top of both panels in Figure 13 indicate the halo mass, M dm , corresponding to the total number of GCs measured from our density profile fits and plotted along the lower horizontal axes.To convert from N GC,T otal to M dm , we used the relation of Harris et al. (2017) η M ≡ ⟨m gc ⟩N GC,T otal /M dm = (2.9 ± 0.2) × 10 −5 (15) and adopted the same average GC mass of ⟨m gc ⟩ = 2.8 × 10 5 M ⊙ used in that study.Harris et al. (2017) established this relation based on the halo mass derived from K-band galaxy luminosity, calibrated through gravitational lensing, and assuming a mean mass-tolight ratio of GCs with ⟨M (GCS)/L V ⟩ ∼ 1.3.They used the total number of GCs from compiling literature data.If we choose to fit the log-log relations shown in Figure 13 with log M gc as the independent variable, we find With halo masses in hand, we now turn our attention to an exploration of how the GC system properties correlate with various halo parameters, and what those correlations imply for the formation of galaxies and their GC systems.

A Check on Halo Masses
In this section, we use our measurements to examine the connection between the GC systems and dark matter halos associated with our program galaxies.Of course, we do not probe the halos directly, so this discussion, like many before it, hinges on the ansatz that the total number of GCs belonging to a galaxy scales in proportion to the halo mass.Before proceeding with the discussion, we therefore pause to perform a simple check on this assumption.
Because our sample includes a number of bright, nearby, well studied galaxies, we can test the key assumption from §3.2 that the observed number of GCs, N GC,T otal is linearly related to the halo mass (i.e., Equation 15). Figure 14 shows a comparison, for 18 of our program galaxies, between the halo mass derived in this way (plotted on the abscissa) and the dynamically measured halo mass (plotted on the ordinate).6Dynamical masses are taken from the SLUGGS survey (Brodie et al. 2014) and are based on GC radial velocity measurements.These masses were derived by using the tracer mass estimator of Watkins et al. (2010) to derive homogeneous masses within 5R e .These masses were then corrected to M 200 (i.e., Table 3

of Alabi et al. 2017).
The dotted line in Figure 14 shows the one-to-one relation, which should reflect the trend if our assumption is correct that N GC,T otal is an accurate tracer of halo mass.This does indeed appear to be the case, with unweighted best-fit relation showing good agreement with the expected relation (with an rms scatter of 0.23 dex).Although this comparison is limited in that the plotted galaxies span only ∼ two orders of magnitude in halo mass, as opposed to the full sample which nearly four orders of magnitude (see Figure 13), we conclude that the available evidence suggests Equation 15 is a faithful predictor of halo mass for our program galaxies.

Connection to Dark Matter Halo Parameters
By fitting the GC density profiles for our program galaxies, we determine two fundamental parameters for each GC system: (1) the total number of GCs, N GC,T otal , which is used to infer the dark matter halo mass ( §3.2); and (2) the effective radius, R e,gc , which allows us to compare the spatial extent of the GC system to that of the underlying halo.For this comparison, we use three different radial measures for the halo, which is assumed to have an NFW profile in all galaxies.
First, we use our GC-based estimate for the mass of the halo, M 200 , to calculate R 200 (e.g., Equation 7 where to find the radius, R h , containing half the halo mass.This provides us with three radial scales for the assumed NFW profile: R s , R h and R 200 . Figure 15 shows relations between R 200 and the GCS effective radius, R e,gc .Plots for the full GC systems are Comparison with various models and other observations.All y-axes represent the effective radius of the GC system, Re,gc.Panels in the left column show relations between the virial radius, R200, and Re,gc, while panels in the right column show relations between host galaxy stellar mass, M * , and Re,gc.We calculated R200 directly from the halo mass (which itself was based on N gc,T otal ), so corresponding halo masses are noted on the top for reference.All panels' solid black, blue, and dashed red lines show the best-fit relations for the total, blue, and red GC systems.The grey dotted lines in each of the left column panels are drawn at fixed fractions of logR200 (i.e., 1%, 2%, 5%, 10%, 20% or 50%).Data points in the bottom panels are color-coded, as shown in the right color bar. Figure 16.Variation in globular cluster system size as a function of host galaxy stellar mass.The panels on the left show results for the full globular cluster systems (black symbols).Results for the separate blue and red globular clusters are shown in the right panels.For both columns, the globular cluster effective radius is normalized to three different measures of the size of the dark matter halo: i.e., the virial radius, R200, the radius containing half the virial mass, R h,dm , and the NFW scale radius, R S,dm (top to bottom, respectively).
Figure 17.Variation in globular cluster system size as a function of inferred halo mass.The panels on the left show results for the full globular cluster systems (black symbols).Results for the separate blue and red globular clusters are shown in the right panels.For both columns, the globular cluster effective radius is normalized to three different measures of the size of the dark matter halo: i.e., the virial radius, R200, the radius containing half the virial mass, R h,dm , and the NFW scale radius, R s,dm (top to bottom, respectively).
shown in the top panels, while the bottom panels show the results of fitting the blue and red GC systems separately.As described above, the symbols in the bottom panel have been color-coded by the ∆ gi parameter from §2.For reference, the dot-dash green lines in these panels show the relations obtained by Reina-Campos et al. ( 2022) based on their simulated GC systems from E-MOSAICS.The shaded green band in each panel correspond to the median and 25-75th percentiles from the E-MOSAICS simulations.In the two panels at the left column, magenta dotted line and cyan dotted line show the relations from simulated GC systems presented by Chen & Gnedin (2022).They also provide the relations for ex situ and in situ GC systems which are shown with cyan dotted line and magenta dotted line, respectively.The shaded cyan and magenta bands in these panel represent 1σ confidence levels.The effective radii of GC systems plotted against the total number of GCs.This is the same as in Figure 13 but with galaxies color-coded by local density, log n10, and MATLAS galaxies highlighted.(Right Panel) Residuals from the best-fit relation in the preceding panel plotted as a function of local density, log n10 but with galaxies color-coded by the total number of GCs.Red solid and blue dashed lines show the best-fit relation for the low mass (N GC,T otal < 100) galaxies and the high mass (N GC,T otal > 500) galaxies, respectively.The blue dotted line shows the relation found when fitting the low mass NGVS galaxies alone.MATLAS galaxies are also highlighted with open magenta squares.
We begin by noting that conclusions drawn from these comparisons should be viewed with some caution for several reasons.First, our halo masses are, of course, not measured directly, but are based on an empirical relationship between number of GCs and halo mass.7 Second, the halo radii plotted in this figure were calculated under the assumption that the halos have NFW profiles, which allows us to link the halo radii to the virial masses; although this is a standard assumption, it is an assumption nonetheless.Finally, the Reina-Campos et al. ( 2022) simulations are predictions for the GC systems of central galaxies, whereas most of the galaxies in our sample are satellites.These caveats aside, we may draw a few conclusions from this comparison.
Beginning with the full GC systems, we see that the scaling relations are continuous over a range of 3.7 decades in stellar mass and have a roughly linear behaviour with surprisingly small scatter -a direct consequence of the tight relationship between R e,gc and N GC,T otal noted in §3.1 (see Figure 13).The rms scatter around the best-fit scaling relations involving R 200 , R h and R s falls in the range 0.13 to 0.20 dex.This is far smaller than the scatter in the simulations, making the tightness of the observed relations remarkable in light of the stochastic and diverse formation paths inherent in the simulations.[Note that the green bands plotted in Figure 15 represent the 25-75th percentiles from the simulations, and magenta band in Figure 15 shows 1σ confidence level].Considering the overall trends, the simulations and observations match reasonably well around R e,gc ∼ 10 kpc -a radius corresponding to stellar and halo masses of ∼ 10 11 M ⊙ and ∼ 10 12.5 M ⊙ , respectively.At the high-and low-mass ends, however, the agreement is less good.At these extremes, both simulations over-predict the GC system sizes in the highest mass galaxies, and Reina-Campos et al. (2022) underpredicts the sizes in the lowest mass galaxies.On the other hand, for the highest-mass galaxies, the simulations are in reasonable agreement with the measurements for the blue GC systems which, having an accretion origin in these galaxies, are more spatially extended than their red counterparts.
As already noted from Figure 10, we find the GC effective radii to vary with galaxy stellar mass in roughly the same way as the galaxies' effective radii -apart from a roughly constant offset of a factor of ∼ 2.5×, in the sense that the GC systems have more extended distributions.This similarity includes a prominent inflection at M p ≃ 6.5 × 10 10 M ⊙ .A standard interpretation for this bend in the galaxy relation is that (mostly) "dry" mergers are responsible for the high masses and steeply increasing radii of galaxies above this mass.As we have shown here (see Figures 10 and 11), this characteristic mass also roughly divides galaxies that have strongly bimodal GC color distributions (at higher mass) from those that do not (at lower masses).Consistent with the basic picture developed to explain the size-mass relationship for galaxies, we interpret this trend as evidence that galaxies at higher mass have grown their GC systems largely through dissipationless accretion of lower mass GC systems, while galaxies below this characteristic mass have largely formed their GCs in situregardless of whether the GCs are metal-poor or metalrich.In other words, at low masses, where two GC subpopulations have similar spatial extents and colors, their origins may not be that different.
The size-M 200 relation of the red GC systems aligns well with the in situ GC systems found in Chen & Gnedin (2022).Chen & Gnedin (2022) also presents this relation for ex situ GC systems, which coincides with the upper limit of the blue GC systems.Interestingly, the size relation of the ex situ simulated GC systems corresponds closely to that of the blue GC systems in massive galaxies.This may suggest that the simulation traces well the pure accreted GCs, while a significant fraction of GCs in low mass galaxies may constitute an in situ population, even though they appear relatively bluer than the galaxy's color as mentioned in Choksi & Gnedin (2019).
We also present a comparison of our results with other simulation results and previous observations in the right column of Figure 15.Consistent with our findings mentioned in Section 3.1, previous observational results align with the outcomes of our study within the ranges where they were constrained by data.Notably, the simulation results of Doppel et al. (2023) exhibit a good match with our findings in the high mass range, but they overpredict the size of GC systems at low masses, contrary to other simulation results.The simulation result by Creasey et al. (2019) traces a totally different GC formation mechanism, focusing solely on accreted GC populations formed in halos that collapse before reionization, and their results are also displayed in the bottom-right of Figure 15.Their simulations produce GC systems that are much larger than the observed blue GCs.While this does not rule out such early-forming GCs, it indicates that the blue GCs need at least one other formation mechanism.
To explore this issue in more detail, the comparisons between the effective radii of GC systems and various DM parameters are shown in Figures 16 and 17.These figures compare ratios between R e,gc and the halo radii, with the panels on the left showing results for the full GC systems, and the panels on the right illustrating the trends after separation into blue and red components.The figures differ in the choice of independent variable, with Figure 16 and Figure 17 showing the trends as a function of galaxy stellar mass and halo mass, respectively.Relative to their dark halos, the GC systems thus do not have a constant fractional extent across the range of stellar and halo masses explored here, pointing to a somewhat more complicated formation and assembly picture.On the other hand, if we exclude the highest mass galaxies (i.e., those with log M * ≳ 10 11 or log M dm ≳ 10 12.5 ), then the apparent variation in fractional size of the GC system is significantly reduced, with an average value relative to the halo scale radius of R e,gc /R s,dm ≈ 0.18 ± 0.03.The fractional size of the red GC systems is more nearly constant across stellar mass, at R e,gc /R s,dm ≈ 0.17 ± 0.01.

Trends with Environment
To conclude this section, we return to the tight correlation reported in §3 between GC system effective radius and total number of GCs.Although the scatter around this relation (Figure 13) is small -with an observed rms scatter of 0.21 dex which, given our measurement errors, implies an intrinsic scatter of just ∼ 0.1 dexit is interesting to note that about ten galaxies lie far above the line and about half of them are MATLAS galaxies.These galaxies are similar in many respects to the Virgo/NGVS galaxies (i.e., they have early-type morphologies, lie along the red sequence, and overlap in stellar mass), so it is natural to ask if this apparent offset is statistically significant and, if so, whether it could be related to a difference in environment.
The left panel of Figure 18 shows the log R e,gclog N GC,T otal relation once again but now with the MATLAS galaxies highlighted.The dashed line in this figure is same to the one in Figure 13.We have colorcoded galaxies in this plot by local density, where the parameter n 10 is the density (in galaxies per arcmin 2 ) measured for all Virgo galaxies by computing the distance to the 10th nearest neighbour (and applying a geometric correction for those galaxies close to the survey boundary).For MATLAS galaxies, we estimate n 10 from the ATLAS 3D local density parameter, Σ (Cappellari et al. 2011), by defining an empirical relation between these two parameters using ATLAS 3D galaxies that fall within the NGVS footprint.Although this bootstrapping does introduce some uncertainty, the color coding in Figure 18 confirms expectations that the MATLAS galaxies occupy less dense environments than do galaxies in the Virgo Cluster -one of the richest environments in the local universe.
The situation depicted in the figure is complex.While most galaxies positioned well above the fitted line are in low-density environments, there are also galaxies in low-density environments that are close to the line (dark symbols in the left panel of Figure 18).These factors make it challenging to discern the environmental effects on the N gc,T otal − R e,gc relation.To simplify the analysis, we refocus on the MATLAS sample (symbols with open squares in the left panel of Figure 18), which represents the lowest-density environments in our dataset.Within the MATLAS galaxies, we observe instances of galaxies both above and below the line, with above-theline galaxies generally being more distant from the fitted line than below-the-line galaxies.Interestingly, most above-the-line MATLAS galaxies have small numbers of GCs with logN GC,T otal ≲ 100.This implies that environmental effects may vary depending on the total mass.
The right panel of Figure 18 shows residuals (∆ log R e,gc ) about the fitted relation as a function of log n 10 , but we divided samples for fitting into subgroups based on the total number of GCs.The solid red line shows the best-fit relation between log n 10 and residuals of massive galaxies with N GC,T otal > 500, and the blue dashed line shows the fitted relation for low-mass galaxies with N GC,T otal < 100.These two datasets reveal clearly distinct trends.Effective radii of GC systems in massive galaxies do not depend on the environments.In contrast, low-mass galaxies show a clear trend for GC systems in lower-density environments to have larger effective radii at fixed N GC,T otal .Given that the log n 10 estimates are more uncertain for MATLAS galaxies, we also show the best-fit relation for NGVS low-mass galaxies alone (the blue dotted line).In this case, the trend is even steeper, although some caution is probably appropriate given the more limited span in density.Regardless of which sample is used, though, the slope of low-mass galaxies is significant at the 2-3σ level.In addition, this trend persists even when the low-mass cut is varied up to N GC,total ≈ 300.This trend, if confirmed by future observations, may indicate a tendency for GC systems of low-mass galaxies in low-density environments to have formed with a more extended spatial distribution.Alternatively, it may be evidence for the dynamical evolution of GC systems in the low-mass halos, with those systems in dense environments having been tidally truncated through interactions.On the other hand, GC systems in massive halos exhibit little dependence on their environments for spatial distributions, implying that massive halos may prevent tidal truncation on the dynamical evolution of GC systems in dense environments.

SUMMARY
The size of a galaxy is a fundamental property that encodes the initial conditions of galaxy formation, as well as a Hubble time of dissipation and dynamics.The relative spatial distributions of a galaxy's gas, stars, and dark matter reflect the dominant physical processes for each component and their correspondingly different evolutionary histories.GCs, which trace the oldest stellar populations and often have the large spatial extent characteristic of stellar halos, have the potential to provide new insights into galaxy assembly histories, although spatial distributions remain one of the most poorly understand properties of GC systems.
We have combined ground-and space-based imaging for early-type galaxies in the nearby universe in order to measure the size and structure of their GC systems, and to explore correlations with the properties of their host galaxies.We have targeted early-type galaxies from the NGVS and MATLAS surveys -two CFHT large programs undertaken with the MegaCam mosaic camera.For most of our sample galaxies, high-resolution imaging is also available from HST/ACS.Our program thus boasts several advantages over previous studies of this sort: i.e., wide-field coverage, high angular resolution in the crowded inner regions, multi-band coverage, depth, uniformity and completeness.The sample itself is both large (N=118) and spans a wide range (a factor of ∼2000) in stellar mass, from dwarf to giant galaxies.
Using a homogeneous method to measure and fit the two-dimensional density distributions of GCs in our program galaxies, we find the following results: 1.The relationship between GC system size and galaxy stellar mass (log R e,gc -log M * ) for earlytype galaxies is best characterized by a broken power law, with an rms scatter of ∼ 0.2 dex, over the range M * ≈ 10 8.5 to 10 11.5 M ⊙ , and a "break" at M p = 6.5 +1.8 −1.3 × 10 10 M ⊙ .The relation is significantly steeper above this mass than below, with power-law slopes of β = 1.30+0.22  −0.17 and α = 0.34 ± 0.04 in the giant and dwarf regimes, respectively.Above a mass of log M * ∼ 10, the size-mass relation of GCs separates into two distinct branches, with the effective radii of the blue

Lim et al.
GCs exceeding those of the red GCs by a factor of ∼ 2.1.At low masses, the GC systems define a single, shallow relation with enhanced scatter compared to high masses.
2. The relation between GC system size and galaxy size (log R e,gc -log R e ) has a roughly linear form over the full range in mass, albeit with significant scatter (σ = 0.3 dex).Although the slope of the best-fit relation is steeper than unity (1.22 ± 0.01), the GC systems have effective radii that are, on average, ∼ 4.0 times larger those of the galaxies.When considering the blue and red GCs separately, these factors are ∼ 5.1 and ∼ 2.2, respectively.
3. We find a remarkably tight relation between the total number of GCs (N GC,T otal ) and the GC system size (R e,gc ), with an rms scatter of just σ ≃ 0.23 dex about the best-fit linear relation.We estimate the intrinsic scatter of this relation to be just ∼ 0.1 dex.Similarly tight linear relationships apply to the blue and red GC systems.The relation defined by the blue GCs is found to be steeper than that of the red GCs -as expected given that the blue components dominate the GC systems of high-mass galaxies.
4. We compare halo masses for our program galaxies -estimated from N GC,T otal using the empirical relation of  2022), the GC systems in our program galaxies exhibit a relatively small scatter about the relations between log R e,gc and the various halo radii (i.e., log R 200 , log R h and log R s ).Typically, the rms scatter about these relations is just ∼ 0.1 dex, which would appear to present a challenge to the current generation of GC system formation models.In addition, although the predicted and observed relations "intersect" at intermediate sizes (R e,gc ∼ 10 kpc) and masses (log M * ∼ 10 11 , log M dm ∼ 10 12.5 ), the agreement is rather poor at the high-and lowmass ends, where the simulations over-predict the GC system sizes in the highest mass galaxies, and show a wide range of possibilities for the lowest mass galaxies.
6. Apart from the highest-mass galaxies (i.e., those with log M * ≳ 10 11 or log M dm ≳ 10 12.5 ), we find the GC systems to have a roughly constant size relative to the halo scale radius, albeit with some scatter: i.e., ≈ 0.23 ± 0.05%.The fractional size of the red GC systems are found to be nearly constant across all stellar masses, at 0.14 ± 0.01%.
7. Our sample galaxies span a significant range in environment -from the densest regions of Virgo's sub-clusters, to the cluster periphery, to the small groups occupied by many MATLAS galaxies.We find some evidence that the deviations from the best-fit size-number (log R e,gc -log N GC,T otal ) relation are correlated with environment in the sense that GC systems of low-mass galaxies in lowdensity regions are, at fixed N GC,T otal , more extended than those in high-density regions.This may indicate that GC systems of low-mass galaxies in low-density environments have more extended distributions at the time of formation or, alternatively, that GC systems of low-mass galaxies in dense environments been truncated by tidal interactions.Otherwise, massive halos may reduce the effect of tidal interactions to spatial distributions of GC systems in dense environments.
There are some obvious extensions to this work.Observations for additional galaxies, chosen to occupy an even wider range in local density, would be useful for confirming and/or characterizing possible trends with environment.High-quality measurements for additional low-mass galaxies would be helpful in understanding the size and structure of the GC systems belonging to these presumed "building blocks" of high-mass galaxies.Such measurements will be challenging given the small numbers of GCs associated with any single dwarf galaxy, so high-resolution, multi-band imaging will be needed to assemble GC samples that have both high completeness and low contamination.We have a follow-up study for the MATLAS dwarfs with HST observation (Marleau et al. in prep.).Dynamical mass measurements for an expanded sample of galaxies, including dwarfs, would make it possible to connect the GC system parameters to those of halo directly, without resorting to empirical scaling relations.On the theoretical front, the next generation of cosmological simulations will need to match these observed scaling relations, including the apparently tight correlations between GC effective radii and host galaxy properties.Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.The specific observations analyzed can be accessed via DOI: 10.17909/yc2t-rr81.

Figure 1 .
Figure 1.Color-magnitude diagram for our program galaxies.Red filled circles, black open circles, and blue crosses show galaxies from NGVS, ACSVCS and MATLAS, respectively.The dashed curve shows the red sequence defined by ∼400 galaxies in the Virgo cluster core (Roediger et al. 2017).

Figure 3 .
Figure 3. Distribution of stellar masses, M * , for our program galaxies (gray histogram).For comparison, the green and blue hatched histograms show the distribution of stellar masses for the samples of Alabi et al. (2017) and Hudson & Robison (2018).There are 32 and 35 galaxies in Alabi et al. (2017) and Hudson & Robison (2018), respectively.

Figure 4 .
Figure 4. Magnitude versus inverse concentration index, ∆m4−8, for sources in the field of NGC4472 (M49, VCC1226), one of our NGVS program galaxies.The ∆m4−8 index represents the magnitude difference between apertures with diameters of 4 and 8 pixels.Globular cluster candidates are selected to fall within the dotted region.

Figure 6
Figure6.Density map for GC candidates in the vicinity of two of our program galaxies, NGC4649 and NGC4621.The image measures 100 ′ × 100 ′ , with a pixel size of 1 ′ × 1 ′ , and has been smoothed using a 2-pixel Gaussian filter.North is up and East is to the left.Inner and outer circles are drawn at radii of Re,gc and 3Re,gc centered on each galaxy.The two GC systems have been fitted simultaneously in our analysis, while two other galaxies that are not members of our sample (NGC4660 and NGC4647) have been masked.

Figure 7 .
Figure 7. Number density radial profile for GC candidates surrounding NGC4564.The filled circles show radial number densities in logarithmically spaced radial bins centered on the galaxy core.The red solid curve shows our best-fit Sérsic (plus background) profile.The orange curves show 100 random resampling results that illustrate the range of fitting errors.The vertical dashed line and horizontal dotted line show the fitted effective radius of the GC system, Re,gc, and the underlying background level, Σ b .

Figure 8 .
Figure 8. Two-dimensional and marginalized posterior probability density functions for the effective radius (Re), Sersic index (n), number density at the effective radius (Σe) and background level (Σ b ) for GC candidates in NGC4564.

Figure 9 .
Figure9.GC color distributions for two of our sample galaxies which illustrate the approach used in this paper to analyze the color dependencies of GC system scaling relations.(Left panel).The GC color distribution for NGC4621, which is classified as bimodal by GMM(Muratov & Gnedin 2010).In this case, the peak colors of the blue (B) and red (C) components are combined with the color of the underlying galaxy (A) to define a color index for each of the blue and red GC subsystems: ∆gi,BGC = (B − A) and ∆gi,RGC = (C − A), respectively, and shown pointing to a corresponding position on the color bar at right.This color will be used in §3.(Right panel) The GC color distribution for NGC4612, which is classified by GMM as unimodal.In this case, we define a single index (B − A) as the difference between the color of the galaxy (A) and the median color (B) of its GC system.

Figure 10 .
Figure10.The effective radii of GC systems plotted against the stellar masses of their host galaxies.Filled black circles show GC systems from this study; well-fitted systems (error of Re,gc < 50%) are shown with large heavy symbols.Stellar masses for Virgo galaxies are taken fromRoediger et al. (2023 in prep.).Stellar masses for MATLAS galaxies are taken from ATLAS 3D , but adjusted to match the NGVS results computed using the subsample of galaxies in common to ATLAS 3D and NGVS (logM * ,NGV S = logM * ,ATLAS 3D − 0.11).The heavy orange line shows the best-fit broken power-law function to the well-fitted systems; light orange lines show 500 random MCMC samples based on this fit.The relation between effective radius and stellar mass for our sample galaxies are shown by the green symbols (binned data) and gray solid curve (model,Côté et al. 2023,  in prep.), respectively.

Figure 12 .
Figure 12. (Left panel) The effective radii of GC systems plotted against the effective radii of their host galaxies.Measurements for the full GC system are shown as black symbols.The heavy orange line shows the best-fit linear relation in log-log space.The light orange lines represent 100 random resampling results based on this fitted relation.The gray dashed line shows one-to-one relation.The green dashed line shows the relation Re,gc = 1.5Re, * which has sometimes adopted in studies of the GC systems belonging to ultra-diffuse galaxies (van Dokkum et al. 2016; Lim et al. 2018).The cyan dashed line shows the relation from Caso et al. (2019).The typical error bar of galaxy effective radii is noted at the bottom right.(Right panel) Same as in the previous panel, except the GC systems have been divided by color.Symbols colors are same as in Figure 11.The blue dotted and red dashed curves show the best-fit linear relations in log-log space for blue and red GCs, respectively.The light blue and light red curves shown 100 random re-sampling results based on these fits.

Figure 11 .
Figure 13.(Left panel) The effective radii of GC systems plotted against the total number of GCs.Measurements for the full GC systems are shown as black symbols; the large, heavy symbols show those galaxies with well-fitted GC density profiles.The best-fit relation is shown by the black line.The upper x-axis shows the dark matter halo mass (M200) corresponding to the computed number of GCs in each galaxy.(Right panel) Same as in the previous panel, except symbols have been color-coded by the color difference between GC systems and their host galaxy, as indicated by the color bar on the right.The blue and red lines show the best-fit fitted linear functions for the blue and red GC subsystems.See §2.4 for details on the separation of blue and red GCs.The cyan line in the left panel shows the relation of Caso et al. (2019) 5 .Comparing our measurements to the oneto-one relation (which is shown as the dotted line in each panel), the GC systems are observed to be significantly more extended the the galaxies themselves -conclusion that holds whether one considers the composite GC system, or separate the blue and red subcomponents.We also note that the trend for red GCs appears to show less scatter than that for the blue GCs (rms = 0.03 vs. 0.06).A comparison of Figure10and 12 shows the correlation between GC size and stellar mass to be significantly tighter than the correlation with host galaxy size.Of course, much of the scatter in Figure12may be attributable to measurement errors in R e, * , which can be significant (see, e.g.,Chen et al. 2010).Finally, we point out that the green dashed line in Figure12shows the relation R e,gc = 1.5R e, * , which has sometimes been adopted as an estimator of GC system size, most notably in the study of ultra-diffuse galaxies (e.g.,van Dokkum et al. 2016).As this figure shows,

Figure 14 .
Figure 14.Comparison of halo masses for the subset of our sample galaxies with dynamical halo mass measurements in the literature.The abscissa shows halo masses inferred from the observed number of GCs, N GC,T otal , as described in §3.2.The ordinate shows dynamically measured halo masses for a sample of 18 galaxies from the SLUGGS survey (Alabi et al. 2017).The dotted line shows the one-to-one relation, while the dashed line shows the best-fit linear relation.

Figure
Figure15.Comparison with various models and other observations.All y-axes represent the effective radius of the GC system, Re,gc.Panels in the left column show relations between the virial radius, R200, and Re,gc, while panels in the right column show relations between host galaxy stellar mass, M * , and Re,gc.We calculated R200 directly from the halo mass (which itself was based on N gc,T otal ), so corresponding halo masses are noted on the top for reference.All panels' solid black, blue, and dashed red lines show the best-fit relations for the total, blue, and red GC systems.The grey dotted lines in each of the left column panels are drawn at fixed fractions of logR200 (i.e., 1%, 2%, 5%, 10%, 20% or 50%).Data points in the bottom panels are color-coded, as shown in the right color bar.(Top left) The dot-dash green line shows a relation based on the simulated GC systems from E-MOSAICS (Reina-Campos et al. 2022) while the shaded green band in each panel corresponds to the median and 25-75th percentiles in each parameter from the simulations.The magenta dotted line and shaded regions shows a relation from simulated GC systems based on the Illustris TNG50-1 (Chen & Gnedin 2022).(Bottom left) Cyan and magenta dotted lines with shaded regions represent relations for in situ and ex situ simulated GC systems from Chen & Gnedin (2022), respectively.(Top right) Dashed cyan (Caso et al. 2019), dotted blue (Forbes 2017), and dashed-dot magenta lines (Hudson & Robison 2018) represent previous observational results, whereas the light green solid line shows a simulation result from Doppel et al. (2023) based on Illustris TNG50.(Right bottom) The green dashed line shows the relation from the simulation by Creasey et al. (2019).
Figure15.Comparison with various models and other observations.All y-axes represent the effective radius of the GC system, Re,gc.Panels in the left column show relations between the virial radius, R200, and Re,gc, while panels in the right column show relations between host galaxy stellar mass, M * , and Re,gc.We calculated R200 directly from the halo mass (which itself was based on N gc,T otal ), so corresponding halo masses are noted on the top for reference.All panels' solid black, blue, and dashed red lines show the best-fit relations for the total, blue, and red GC systems.The grey dotted lines in each of the left column panels are drawn at fixed fractions of logR200 (i.e., 1%, 2%, 5%, 10%, 20% or 50%).Data points in the bottom panels are color-coded, as shown in the right color bar.(Top left) The dot-dash green line shows a relation based on the simulated GC systems from E-MOSAICS (Reina-Campos et al. 2022) while the shaded green band in each panel corresponds to the median and 25-75th percentiles in each parameter from the simulations.The magenta dotted line and shaded regions shows a relation from simulated GC systems based on the Illustris TNG50-1 (Chen & Gnedin 2022).(Bottom left) Cyan and magenta dotted lines with shaded regions represent relations for in situ and ex situ simulated GC systems from Chen & Gnedin (2022), respectively.(Top right) Dashed cyan (Caso et al. 2019), dotted blue (Forbes 2017), and dashed-dot magenta lines (Hudson & Robison 2018) represent previous observational results, whereas the light green solid line shows a simulation result from Doppel et al. (2023) based on Illustris TNG50.(Right bottom) The green dashed line shows the relation from the simulation by Creasey et al. (2019).

Figure 18 .
Figure18.(Left Panel) The effective radii of GC systems plotted against the total number of GCs.This is the same as in Figure13but with galaxies color-coded by local density, log n10, and MATLAS galaxies highlighted.(Right Panel) Residuals from the best-fit relation in the preceding panel plotted as a function of local density, log n10 but with galaxies color-coded by the total number of GCs.Red solid and blue dashed lines show the best-fit relation for the low mass (N GC,T otal < 100) galaxies and the high mass (N GC,T otal > 500) galaxies, respectively.The blue dotted line shows the relation found when fitting the low mass NGVS galaxies alone.MATLAS galaxies are also highlighted with open magenta squares.
ACKNOWLEDGMENTSS.L. acknowledges the support from the Sejong Science Fellowship Program by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No.NRF-2021R1C1C2006790).CL acknowledges support from the National Natural Science Foundation of China(NSFC, Grant No. 12173025,  11833005, 11933003), 111 project (No.B20019), and Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education.C.S. acknowledges support from ANID/CONICYT through FONDE-CYT Postdoctoral Fellowship Project No. 3200959.O.M. is grateful to the Swiss National Science Foundation for financial support under the grant number PZ00P2_202104.IRAF was distributed by the National Optical Astronomy Observatory, which was managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii.The observations at the Canada-France-Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site.Funding for SDSS-III has been pro-vided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science.The SDSS-III web site is http://www.sdss3.org/.SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.These observations are associated with program GO-9401.

Table 1 .
Sample galaxies and their properties.

Table 2 .
GC System Size versus Stellar Mass RelationAll GCs Blue GCs Red GCs