The Spatial Distribution of Globular Cluster Systems in Early Type Galaxies: Estimation Procedure and Catalog of Properties for Globular Cluster Systems Observed with Deep Imaging Surveys

The data used in this study were obtained from two large optical imaging surveys: NGVS and MATLAS. While imaging for numerous galaxies is available from these surveys, we limited our targets to ensure reliable GC studies, mainly focusing on early-type galaxies (ETGs). First, we targeted nearby ($\leq 25$ Mpc) MATLAS galaxies (including galaxies inside the NGVS footprint) that were observed in at least three filters. Second, we also targeted ACSVCS galaxies inside the NGVS footprint to include low-mass (early-type dwarf) galaxies in our sample. Among the galaxies that satisfy the above categories, several systems were too close to neighboring, giant galaxies with their own rich GC systems, making it almost impossible to detect and study GCs in these galaxies; such objects were excluded for our analysis. Ultimately, we targeted 118 galaxies in this study. Table 1 lists our targets and their properties.

The NGVS is a deep, multi-band imaging survey of the Virgo cluster carried out with MegaCam (Boulade et al., 2003) on the Canada-France-Hawaii Telescope (CFHT) from 2008 to 2013. The survey covers an area of 104 deg² (with 117 pointings) within the virial radii of both the Virgo A and Virgo B subclusters. Full survey details, including observing strategy and data processing, are described in Ferrarese et al. (2012). Additional details on the data reduction and analysis procedures are also available in Ferrarese et al. (2020).

The MATLAS survey is a second deep imaging survey using MegaCam on CFHT. The targets of the MATLAS survey are galaxies from the ATLAS^3D sample (Cappellari et al., 2011). The sample contains 260 nearby (within 42 Mpc) bright ($M_{K}<-21.5$) early-type galaxies. Full survey details are available in Duc et al. (2015); Duc (2020); Bílek et al. (2020), and the data reduction process is described in Gwyn (2008).

We used the NGVS aperture photometry catalog for the NGVS samples. The catalog details are fully described in Liu et al. (2015), so we present here only a brief description of the catalog. Source Extractor (Bertin & Arnouts, 1996) was run on the processed images to obtain aperture magnitudes of sources with dual-image mode. We used the g’-band images as detection images and adopted circular apertures with a series of diameters between 2 and 16 pixels to measure the source fluxes, which were then corrected to 16-pixel-diameter aperture magnitudes. Instrumental magnitudes were then calibrated to standard AB magnitudes through a comparison to SDSS PSF magnitudes after conversion to MegaCam filter magnitudes.

The central regions of some galaxies have high surface brightness, making it difficult to estimate the background and detect sources with general photometry programs such as Source Extractor. Since all NGVS galaxies have been modeled with customized two-dimensional isophote (ISO) fitting models (Ferrarese et al., 2020), we subtracted diffuse light from galaxies using ISO fit models to enhance source detection and background estimation. The model subtraction is performed with a cutout images having a $10\arcmin\times 10\arcmin$ field of view. We then ran Source Extractor on these galaxy-subtracted images in the same way as for the NGVS source catalog. The magnitudes measured on the model-subtracted images are then matched with the original NGVS catalog by comparing magnitudes of sources in the outer $2\arcmin$ width area of model-subtracted images. We replaced the NGVS aperture photometry catalogs of the central $8\arcmin\times 8\arcmin$ regions of our target galaxies with the photometric catalogs on the model-subtracted images.

We also generated aperture photometry catalogs for MATLAS galaxies using the same methodology as the NGVS aperture photometry catalog. However, the detection images and galaxy model subtraction for MATLAS galaxies differ slightly from those for NGVS galaxies. For the detection image, we chose the best seeing filter image for MATLAS. We used ring median filtered galaxy models instead of two-dimensional ISO fit models for the model subtraction. The ring median filtering method can produce diffuse images by setting the inner and outer ring sizes. We set the radii of the inner and outer rings to $15$ and $20$ pixels, respectively. We subtracted these ring-median-filtered model galaxies from the MATLAS cutout images with a field of view of $10\arcmin\times 10\arcmin$, similar to the NGVS data. Additionally, we replaced the MATLAS aperture photometry catalogs for the central $8\arcmin\times 8\arcmin$ regions of our target galaxies with photometric catalogs based on the model-subtracted images.

We selected GC candidates based on a combination of size information and colors. Since most GCs at the distances of our target galaxies appear as point-like, or slightly extended, sources in the images, we chose point-like sources first based on the inverse concentration index, $\Delta m^{\prime}_{4-8}$. This is the difference in magnitudes between apertures of 4-pixel and 8-pixel diameters. These aperture magnitudes are corrected for missing point source fluxes, so the $\Delta m^{\prime}_{4-8}$ value of point sources is defined to be zero. We measured $\Delta m^{\prime}_{4-8}$ values using the g’- and i’-band images for the NGVS targets and combined them with error-weighting. As for the MATLAS targets, we calculated $\Delta m_{4-8}$ values on the best seeing filter image and the second best seeing images. These $\Delta m_{4-8}$ values are also combined with error-weighting. We used these error-weighted mean $\Delta m_{4-8}$ values for the point source selection. The $\Delta m_{4-8}$ values of point sources show scatter with a mean of zero due to photometric errors. Therefore, we chose point-like sources with a range of $\Delta m_{4-8}$ values, $-0.08\leq\Delta m_{4-8}\leq 0.08$. We limited the selection to sources brighter than $g^{\prime}=24.5$ mag to mitigate the effect of large photometric errors. GCs at smaller galactocentric distances can be partially resolved in high-quality images, so we expanded the range of IC values as $-0.08\leq\Delta m_{4-8}\leq 0.16$ for galaxies within $20$ Mpc for the NGVS targets and for galaxies within $20$ Mpc in MATLAS having high image quality (Seeing $\leq 1.0\arcsec$). After selecting point-like sources, we used color information to choose GCs. We chose GC candidates using polygons in the $(u^{*}-g^{\prime})-(g^{\prime}-i^{\prime})$ color-color diagrams (Figure 1) when u-band data is available. Otherwise, we used $(g^{\prime}-r^{\prime})-(g^{\prime}-i^{\prime})$ color-color diagrams (Figure  2). All GC selection polygons are defined based on the M87 spectroscopically confirmed GCs (see Lim et al., 2017).

Although model-subtracted images allow us to detect additional sources in the central regions of galaxies, the GC samples will still be incomplete in the cores of bright galaxies. Because we have targeted galaxies in the ACSVCS, we also have HST/ACS photometric catalogues that are much more complete than is possible with ground-based imaging. Our analysis thus uses the GC catalog from ACSVCS (Jordán et al., 2007b) which provides $g_{\rm ACS}$ and $z_{\rm ACS}$ magnitudes, and GC probability. We transformed these HST magnitudes to CFHT $g$ and $i$ magnitudes using the following linear equations derived from matched sources in the M87 region.

where $g^{\prime}_{\rm NGVS}$, $i^{\prime}_{\rm NGVS}$, $g_{\rm ACS}$, and $z_{\rm ACS}$ represent magnitudes of CFHT $g^{\prime}$, CFHT $i^{\prime}$, HST/ACS $F475W$, and HST/ACS $F850LP$, respectively.

For the ACSVCS sources, we define GCs with ${\rm p}_{GC}\geq 0.5$, and $g^{\prime}_{\rm NGVS}\leq 24.5$ mag, where ${\rm p}_{GC}$ is the probability that the object is a GC. We matched these GC candidates to those from the NGVS catalog, and those without a match are added to the total GC samples for further analysis.

To estimate the limit of our ground-based photometry, we performed completeness tests by injecting artificial stars into the $10\arcmin\times 10\arcmin$ cutout images. Artificial stars were constructed using Point-Spread Functions (PSFs) empirically obtained for each observing field. PSFs were generated with DAOphot, and PSF stars were selected using SExtractor, with inverse concentration indices applied for point source selection. Detailed information on PSF generation for the NGVS is available on the NGVS webpage, and the process for generating PSFs for MATLAS data closely follows that of the NGVS. When adding these artificial stars, we used a power-law luminosity function with a magnitude range of $20<g^{\prime}<25$. Our GC survey is limited to $g^{\prime}_{0}=24.5$, so we set the faint magnitude limit to $g^{\prime}=25$. When performing these tests, the number of added artificial stars did not exceed $10\%$ of the number of sources with the same magnitude range. We repeated this experiment over a thousand times, resulting in more than $150,000$ artificial stars. On average, we added about $200,000$ artificial stars to each target image. As expected, the completeness varied with both the magnitude of sources and the background brightness, which generally varies radially from the center of galaxies for ETGs. To account for this variation, we divided the artificial star test into subgroups based on the radial distance from galaxy center and fitted the results with a step function as follows:

where $\alpha$ is the slope of the decreasing recovery rate; $m_{50}$ is the magnitude for the $50\%$ recovery rate; and $m$ is the magnitude of sources. By applying this method, we obtained the recovery rates that vary with source magnitudes and distances from the galaxy center.

We found that the ground-based photometry suffers from incompleteness near center of galaxies for the sources with $g^{\prime}\sim 24.5$ mag, whereas the ACSVCS catalog is $100\%$ complete for the sources with $g^{\prime}_{0}\leq 24.5$ mag in the entire fields.

We used an analytic function to investigate the spatial distributions of GCs, especially the GC number density profiles. We fitted these profiles using a two-dimensional Sérsic function plus a constant background, given by:

where

Here $\theta$ is the position angle of the GC candidate, measured from north to east, while $X_{0}$ and $Y_{0}$ are the coordinates of the center of the host galaxy. $X$ and $Y$ are the coordinates of the GC candidates, and $\epsilon$ is the ellipticity. In the Sérsic function, $R_{e}$ is the effective radius, $\Sigma_{e}$ is the GC number density at the effective radius, $\Sigma_{b}$ is the background GC number density, $n$ is the Sérsic index, and $b_{n}$ is a constant that depends on $n$. We used an approximation of $b_{n}$ (Ciotti & Bertin, 1999).

We fitted this function to the data using a Markov Chain Monte Carlo (MCMC) method with the emcee code in python (Foreman-Mackey et al., 2013). We set flat prior distributions on several parameters, including $\Sigma_{e}>0.001$ arcmin^-2, $0.25<n<8.0$, and $0.05<R_{e}<30$ arcmin. Due to difficulties in constraining both ellipticity and position angle in the presence of GC contamination, we imposed priors of $0\leq\epsilon<0.1$ and $-10^{\circ}<\theta<10^{\circ}$ for most cases. For a handful of elongated galaxies, we also imposed priors of $0\leq\epsilon<0.4$ and $-180^{\circ}<\theta<180^{\circ}$. We adopted a prior from a Gaussian function with pre-estimated values of the background mean and background standard deviation for the background level.

To obtain the logarithmic probability, we used the following likelihood function:

Here, $\ell_{i}(R_{i}|\Sigma_{e},n,R_{e},\Sigma_{n})$ is the probability of finding datum $i$ at radius $R_{i}$ given the Sérsic parameters. We applied a completeness correction to each data point $R_{i}$ for galaxies without HST data. This completeness correction involved convolving the source detection probability function for CFHT/MegaCam data with the intrinsic luminosity function, which allowed the density probability function to be corrected for each data point.

To construct the probability function, we need to integrate the modified Sérsic function. As our modified Sérsic function does not have a theoretical integrated function form, we used a numerical integration method to generate the probability function. The integration range is normally from the radius corresponding to the 50% completeness limit to 30^′. If ACSVCS imaging is available, the integration starts at the galaxy center. For the largest galaxies, the outer limit was taken to be 60^′.

We also masked unfavourable areas for GC density profile fitting: i.e., those regions affected by saturated stars, nearby galaxies, or outside of the observation fields, etc. In the case of target galaxies with neighbors, we fit two 2D Sérsic functions, simultaneously, to include the contribution from the the neighboring galaxy. However, if there are more than two neighboring target galaxies, then we masked them except for the most dominant neighbor.

The colors of GCs in early-type galaxies are often used as a metallicity indicator, and they typically show unimodal or bimodal distributions. We used Gaussian Mixture Modelling (GMM, Muratov & Gnedin, 2010) to test for color-bimodality among the GC systems. To create input catalogs for the GMM test, we selected GCs within a radius of $2.5R_{e,gc}$, and used their $(g-i)$ colors. However, these GCs can be contaminated by various sources, including intra-cluster GCs, foreground stars, and background galaxies. To deal with background contamination, we followed several steps: (1) We defined a background area outside the galaxy region and calculated the areal fraction between the background field and the galaxy field (within $2.5R_{e,gc}$); we then randomly chose background objects based on this areal fraction; (2) We subtracted GCs in the galaxy region that had the closest colors to the GC candidates in the background region; (3) Using this background-subtracted catalog, we ran the GMM code and obtained the results; (4) We repeated this process 30 times to account for errors in background subtraction. After these steps, we calculated the mean and standard deviation of D-values, representing the distances between the peaks of the fitted Gaussian functions. If the D-value was greater than 2, it indicates the data was better fit by two Gaussian functions, and we considered the distribution bimodal when the D-values were greater than 2 standard deviations. For GCs with a D-value greater than 2, we divided them into blue and red GCs based on the crossover values of the two fitted Gaussian functions.

The total number of GCs is a fundamental parameter in the study of GC systems. However, estimating this number requires overcoming the observational limits of magnitudes and spatial coverage. In this study, we used the 2D-Sérsic fitting of the GC distribution to estimate the total number of GCs above the magnitude limit ($g^{\prime}_{0}=24.5$ mag). We integrated the fitted 2D-Sérsic function to obtain the magnitude-limited total number of GCs.

It is well known that GC luminosity functions (GCLFs) have a roughly universal Gaussian form with a peak luminosity of $M_{V}\sim-7.5$ (e.g., Harris, 2001; Lee, 2003). More precisely, peaks and widths of GCLFs change slightly depending on the host galaxy luminosity (e.g., Jordán et al., 2007b; Villegas et al., 2010). We made a similar assumption by adopting different peak luminosities and widths of the Gaussian function based on the host galaxy luminosity, using the relation from Villegas et al. (2010). Villegas et al. (2010) also provides mean and sigma values for the GCLF of many target galaxies in this study, so we used these values directly for these common targets. With these GCLFs, we could correct the observation’s magnitude limit and obtain the total number of GCs. We were also able to estimate the errors of the GC total numbers based on the errors of the Sérsic fitting, but we did not include the errors from GCLF. Note that the total number of GCs for ACS-VCS galaxies has been calculated based on the HST data (Peng et al., 2008); we will compare both numbers and discuss them in the next section.

Table 2 presents the estimated properties of the GC systems in our target galaxies. The table provides information on the effective radii and total numbers of GC systems with their $1\sigma$ errors. Additionally, we include the Sérsic indices, $n$, for all GC systems. The table also notes the bimodality of GC colors, indicating the effective radii and peak colors of the blue and red GCs for those galaxies that exhibit a bimodal GC color distribution. The table also includes the median GC colors for those galaxies exhibiting an unimodal GC color distribution. Furthermore, the GC specific frequency is also noted (see the bottom of this section).