Nonlinear Color-Metallicity Relations of Globular Clusters. XI. Nonlinearity Effect Revealed by NGC 5128 (Centaurus A) and NGC 4594 (Sombrero) Galaxies

Metallicity distributions (MDs) of globular clusters (GCs) provide crucial clues for the assembly and star formation history of their host galaxies. GC colors, when GCs are old, have been used as a proxy of GC metallicities. Bimodal GC color distributions (CDs) observed in most early-type galaxies have been interpreted as bimodal MDs for decades, suggesting the presence of merely two GC subpopulations within single galaxies. However, the conventional view has been challenged by a new theory that nonlinear metallicity-to-color conversion can cause bimodal CDs from unimodal MDs. The unimodal MDs seem natural given that MDs involved many thousand protogalaxies. The new theory has been tested and corroborated by various observational and theoretical studies. Here we examine the nonlinear nature of GC color-metallicity relations (CMRs) using photometric and spectroscopic GC data of NGC 5128 (Centaurus A) and NGC 4594 (Sombrero), in comparison with stellar population simulations. We find that, with a slight offset in color, the overall shapes of observed and modeled CMRs agree well for all available colors. Diverse color-depending morphologies of GC CDs of the two galaxies are well reproduced based on their observed spectroscopic MDs via our CMR models. The results corroborate the nonlinear CMR interpretation of the GC color bimodality, shedding further light on theories of galaxy formation.


INTRODUCTION
Globular clusters (GCs) are among the oldest stellar systems in the Universe. They are thought to have formed during the major star formation events in galactic history and survived relatively intact throughout the galactic evolution. Since they retain the chemical properties of the surrounding environment at the formation epoch, their metallicity distribution (MD) provides a vital clue to the formation and evolution of galaxies (Bastian & Lardo 2018). Due to the limitations of observational instruments, GC MDs are often inferred from photometric data rather than obtained directly through spectroscopy. GC colors have been used as proxies for * Both authors have contributed equally to this paper.
GC metallicity under the assumption that, when GCs are old, colors correlate linearly with metallicity.
One of the most important features observed in GC color distribution (CD) is bimodality, in that the CD of GCs belonging to one galaxy exhibits two peaks (Brodie & Strader 2006). Since this phenomenon does not appear only in certain galaxies, but in most large galaxies, proper interpretation of this phenomenon is crucial to the study of galaxy formation and evolution. With empirical linear GC color-metallicity relations (CMRs), bimodal CDs have generally been translated to bimodal MDs, implying the presence of two populations of GCs within individual galaxies (e.g., Park & Lee 2013;Lee & Jang 2016;Forbes & Remus 2018;El-Badry et al. 2019). However, Yoon et al. (2006, hereafter Paper I) give a plausible alternative explanation for the origin of GC color bimodality, in which the nonlinear nature of GC CMRs can lead to the bimodal CDs even from arXiv:2209.02738v1 [astro-ph.GA] 6 Sep 2022 a single-peaked MD. In this theory, the underlying GC MDs could be quite different from that of GC CDs characterized by bimodality. Although many studies have attempted to examine the new explanation by various means (e.g., Cantiello & Blakeslee 2007;Kundu & Zepf 2007;Strader et al. 2007;Peng et al. 2008;Spitler et al. 2008;Blakeslee et al. 2010;Schuberth et al. 2010;Blakeslee et al. 2012;Brodie et al. 2012;Chies-Santos et al. 2012;Park et al. 2012;Pota et al. 2013;Park & Lee 2013;Vanderbeke et al. 2014;Cantiello et al. 2014;Richtler et al. 2015;Salinas et al. 2015;Cho et al. 2016;Villaume et al. 2019;Fahrion et al. 2020), this issue is still under debate. Elucidating what the nonlinearity is like in GC CMR and how it affects the color-tometallicity transformation is at the heart of this debate.
Since Paper I, a series of papers (Yoon et al. 2011a(Yoon et al. ,b, 2013Chung et al. 2016;Kim & Yoon 2017;Lee et al. 2019Lee et al. , 2020Kim et al. 2021, hereafter Papers II − X) have been scrutinized the nonlinear CMR theory, both observationally and theoretically. Papers II and IV reported strong color-dependent variations in the CD morphology of M87 and of M84 GCs, respectively, which were incomprehensible with linear, straight CMRs but were well reproduced by the nonlinear CMRs. Paper III proposed a highly viable solution based on the CMR nonlinearity for the perplexing conundrum of the MD discrepancy between GCs and field halo stars in given parent galaxies. Papers V and VII paid attention to bimodal distributions of absorption-line indices such as Balmer and Mgb lines for GCs in highly studied M31 and NGC 5128, respectively. The bimodality is readily explained by nonlinear index−metallicity relations, analogous to the nonlinear CMR projection. Paper VI approached theoretically to the MDs derived using the Ca Triplet index, a widely used metallicity indicator. They showed that the Ca Triplet−color−metallicity relation is highly nonlinear, resulting in a sharp difference between color and Ca Triplet distributions of old GCs. Paper VIII proposed a simple, coherent explanation for the observed diversity in the GC CDs of early-type galaxies in the Virgo and Fornax galaxy clusters, where the diversity is due mainly to changes in the mean metallicity of GC systems. Paper IX explained the difference in the radial number density profile between blue and red GCs by combining the radial metallicity gradient of the GC system and the nonlinear nature of metallicity-to-color transformation. Paper X observed GCs in M87 with the Subaru/FOCAS multiobject spectroscopy mode. Their high-quality, homogeneous dataset corroborated the theoretical prediction of inflected CMRs.
As part of the series papers, we investigate the presence and degree of nonlinearity along GC CMRs for var-ious colors by combining multiband photometric data with spectroscopic data of the two nearby massive galaxies, NGC 5128 (Centaurus A) and NGC 4594 (Sombrero). The photometry of NGC 4594 GCs is from our observations and the other data are taken from the literature. The NGC 5128, the central giant elliptical galaxy in the Centaurus group of galaxies, is one of the best targets to study the nature of GC CMRs thanks to its proximity (∼4 Mpc), the richness of GCs (> 3000; Taylor et al. 2017), and relatively low differential reddening of E(B − V )=0.07 (Rejkuba et al. 2002). The NGC 4594 is a lenticular galaxy with an unusually large bulge and a prominent dust lane in the nearly edge-on disk. Since it is also nearby (∼9.55 Mpc) and has a relatively large number of GCs (∼1150; Larsen et al. 2001), this galaxy's GC system is well suited for GC CMR studies.

NGC 5128 (Centaurus A) Globular Clusters
The photometry for GCs in NGC 5128 was obtained from the compilation of Woodley et al. (2007), which combined the U BV RI photometry by Peng et al. (2004) and the CM T 1 photometry by Harris et al. (2004). We used the U BV RI data because CMRs of CM T 1 have too large scatter to examine the CMR nonlinearity. To compare the photometry with theoretical predictions, we applied a correction for Galactic extinction based on the COBE/DIRBE dust maps (Schlegel et al. 1998) with the correction terms by Schlafly & Finkbeiner (2011) and the extinction law by Cardelli et al. (1989).
The spectroscopic data were taken from Beasley et al. (2008). They obtained integrated light spectra for 254 GCs in NGC 5128 using the 2-degree Field instrument on the Anglo-Australian Telescope and provided 21 Lick indices for 146 GCs and empirical spectroscopic metallicities for 243 GCs. The left panel of Figure 1 shows the spatial distribution of photometric data of Peng et al. (2004, black dots) and Woodley et al. (2007, blue dots) and spectroscopic data of Beasley et al. (2008, red dots).
We crossmatched these data using a matching radius of 1 , resulting in a total number of 234 GCs with both the U BV RI photometry and spectroscopic metallicity. The empirical metallicity ([M/H]) provided by Beasley et al. (2008) was transformed to [Fe/H] by subtracting 0.217 (Kim et al. 2002). Table 1 gives astrometry, spectroscopic metallicity and U BV RI photometry of crossmatched GC sample in NGC 5128.

NGC 4594 (Sombrero) Globular Clusters
We carried out U BV I photometry for the GC system in NGC 4594 on the nights of 2008 May 4 and 5 (UT) with the Mosaic II CCD imager mounted on the prime focus of the 4 m Blanco telescope at Cerro Tololo Inter-American Observatory (CTIO). The observing conditions were clear and photometric with an average seeing of around 1 . The Mosaic II consists of eight 2048×4096 pixel CCDs with a pixel scale of 0. 27 providing a total field of view of 36 ×36 . The total exposure times in U , B, V , and I were 18000, 4500, 3000, and 2000 seconds, respectively. We split our observations in each band into five dithered exposures (two sets for U ) to fill in the gaps between the individual CCD chips in the mosaic, and to minimize the impact on our photometry from CCD blemishes such as hot pixels and bad columns.
All images were preprocessed using the MSCRED package (Valdes 1998) in IRAF 1 . We used the CCD-PROC task to correct for cross-talk between the CCD chips, trim the overscan, correct for the bias level, and apply the flat-field correction. The MSCPIXAREA task is used to correct for the variations in pixel scales across the frame produced by geometric distortion in the Mosaic II imager. Each mosaic image was then split into eight individual images using the MSCSPLIT task. We combined the images using the MONTAGE2 in the ALLFRAME package (Stetson 1994) to produce a clean high signal-to-noise coadded image in each filter. Sources were detected from the combined images using the DAOPHOT II/ALLSTAR (Stetson 1987), producing a master catalog of the sources. The master catalog was then input to the ALLFRAME along with all of the images and their point-spread function (PSF) models for final PSF photometry. Aperture corrections were computed through the DAOGROW (Stetson 1990) using bright and isolated point sources, and we applied them to our photometry above. Finally, an error-weighted mean instrumental magnitude for each object was obtained using the DAOMATCH/DAOMASTER (Stetson 1993).
The photometric calibration was achieved using standard stars in the fields of Landolt (1992) and Stetson (2000). We corrected for the foreground Galactic extinction using the reddening maps by Schlegel et al. (1998) with the correction terms by Schlafly & Finkbeiner (2011) and the extinction law by Cardelli et al. (1989). Astrometric solutions were calculated using the USNO-B 1.0 catalog stars (Monet et al. 2003), resulting in the typical r.m.s. level of ∼ 0. 2.
The spectroscopic data were taken from Alves-Brito et al. (2011), who measured three metallicity-sensitive indices (CH,Mgb,& Fe5270) for 247 GC candidates and transformed the measurements into metallicities using the Brodie & Huchra (1990) method. They provided error-weighted mean metallicities for 135 GC candidates for which all three indices are measured. From the sample, we selected GCs using radial velocity and projected galactocentric distance (see Figure 3 in Alves-Brito et al. 2011) and used their mean metallicities for our analysis. The right panel of Figure 1 shows the spatial distribution of photometric GC candidates (see Section 2.4, gray dots) and spectroscopic data (yellow and red dots).
We crossmatched our U BV I catalog with the spectroscopic data using a matching radius of 1 , resulting in a total number of 122 GCs with both U BV I photometry and spectroscopic metallicity. Table 2 gives GC ID, coordinate, spectroscopic metallicity, and U BV I photometry of crossmatched GC sample in NGC 4594.

Theoretical Evolutionary Population Synthesis Model
The theoretical models used in this study are constructed using the Yonsei Evolutionary Population Synthesis (YEPS) code (Chung et al. 2013b). The YEPS model comprehensively deals with the effect of horizontal-branch stars of simple stellar populations, which are the main source causing the nonlinearity of GC CMRs. GC colors, when GCs are old, have been used as a proxy of GC metallicities. We use the YEPS model of the best-matching age (12.5 Gyr) for both NGC 5128 and NGC 4594. We adopt the enhanced αelements ([α/Fe] = 0.3) and solar-scaled abundance ratios for other elements. Readers are referred to Chung et al. (2013bChung et al. ( , 2017 for details of the ingredients and input parameters of the YEPS model. Table 3 provides the YEPS predictions of metallicity and UBVRI colors for the 12.5 Gyr simple stellar population, which is used in our present analysis.

Selection of Old Globular Clusters
In order to properly examine the nonlinearity of CMRs for old GCs, possible young and intermediate-age GCs should be excluded from the sample. Both NGC 5128 and NGC 4594 have strong dust lanes across the central regions, hinting at recent wet merger events and subsequent star formation. The cosmological hydrodynamical simulations by Du et al. (2021) suggested that Sombrero-like galaxies (halo-dominated galaxies with huge bulges) have experienced residual star formation after major merger events during the formation history. It is thus naturally expected that young/intermediate-age GCs are intermingled with the dominant old GCs of the galaxies.
The GC age from integrated light is usually estimated from the comparison of age-sensitive indices such as Balmer lines and metallicity-sensitive indices in spectroscopy (Puzia et al. 2005(Puzia et al. , 2006Worthey 1994), or from two-color diagrams of the combination of age-and metallicity-sensitive colors in photometry (Puzia et al. 2002;Hempel et al. 2003;Yi et al. 2004;de Grijs et al. 2005). We chose the photometric method to select old GCs from the crossmatched data to secure as a large number of GCs as possible for the CMR nonlinearity test. Figures 2 and 3 show the two-color diagrams of GCs in NGC 5128 and in NGC 4594, respectively. The selection region is denoted by the polygon drawing by thick black lines in each panel, and the YEPS models for ages from 1 Gyr (blue) to 13 Gyr (red) at 2 Gyr intervals are overplotted for reference. The model loci are slightly shifted to the direction of observations using the measured offsets, which will be discussed in the following section. The use of U -and B-band data enables us to discern young/intermediate GCs from old GCs, especially for metal-poor GCs, since young/intermediate GCs are bluer in U − B and B − V colors. We determine the selection regions to exclude these blue GCs in U − B and B−V colors and the outliers laid away from the GC clustering in the two-color diagrams. However, we note that, as young GCs with intermediate and high metallicities are indistinguishable from old GCs given the observational uncertainties, they are eliminated only partially, contaminating the sample to some extent.

COLOR-METALLICITY RELATIONS
It has been known that there are some offsets between observations and models in the color-magnitude planes (VandenBerg & Clem 2003;VandenBerg et al. 2010;Worthey & Lee 2011;Angelou et al. 2015) and in the integrated colors (Barmby & Huchra 2000;Lee et al. 2002;Conroy et al. 2009). The discrepancy stems mainly from the incompleteness of models, such as imperfect treatments of binaries, blue stragglers, horizontal-branch stars, and stars in the late evolutionary stages (Chung et al. 2013a(Chung et al. , 2017Jang et al. 2021). To compensate for the defect, we shift the model CMRs in the direction of colors to match the observed CMRs, allowing a more direct comparison of the model predictions with observations. The amount of shifts is determined using the total least squares regression method with sigma clipping. After transforming colors and metallicities into a normalized format, we calculate the inversely error-weighted total least squares between observations and models by mov-ing the model loci by the interval of 0.01 mag in the color direction. We then find the amount of shifts, which minimizes the sum of the orthogonal residuals. Figure 4 shows an example of the process of determining the offset value between observations and models. The dotted line represents the original model prediction of the GC CMR and the red solid line shows the shifted model locus after applying the offset value denoted in the top-left corner of the figure. The orange lines show the orthogonal distances from the data points to the model CMR. The black circles are the final sample used in the total least squares method and the open circles are outliers removed by 5σ clipping. Figure 5 shows the CMRs of NGC 5128 GCs in nine color combinations. The size and opacity of each symbol are set to be inversely proportional to the quadratic sum of errors calculated after normalization. The red dashed and solid lines respectively represent the original 12.5 Gyr YEPS model predictions and the models shifted by the offset value denoted on the top-left corner of each panel. We verify that the overall shape of the CMRs agrees well between observations and models for all the colors. The scatters may result from both observational uncertainties and the spread of intrinsic properties of GCs such as age and α-element abundance. It is a common feature that the slope of GC CMRs is steep in the blue region and decreases as the color becomes redder, which agrees with the features of the YEPS model prediction. The characteristics are also consistent with previously suggested nonlinear empirical CMRs such as broken-line CMRs (e.g., Peng et al. 2006;Usher et al. 2012) and high-order polynomial CMRs (e.g., Blakeslee et al. 2010;Vanderbeke et al. 2014).
The presence of a quasi-inflection point along BV RI CMRs is an essential feature in view of the color bimodality origin debate. The observations, although inconclusive, give hints of inflection along the BV RI CMRs [panels (d)-(i)]. A paucity of GCs occurs at the vicinity of the quasi-inflection point where the color changes faster than the metallicity. The paucity near the inflection point is evident along the color axis, which makes CDs bimodal even from the unimodal MD (Paper I). We will investigate this effect further in the next session. Figure 6 shows the CMRs of NGC 4594 GCs for six color combinations. The symbols and lines are the same as in Figure 5. Although the scatters are rather larger compared to the NGC 5128 GC CMRs, there are similar features in that the CMR slope is steeper in the blue and the slope-changing point appears to be consistent with the model prediction.
Noticeably, there are outliers away from the model CMR predictions, particularly among red GCs. To scrutinize the nature of these GCs, in Figure 7, we present their metallicities as a function of galactocentric distance (left), Mgb vs. [Fe/H] along with model predictions for different α-elements mixture (middle), and the B − V CMR along with the 5 Gyr and 12.5 Gyr model predictions for [α/Fe] = 0.0 and 0.3 (right). For sanity check, the middle panel includes the MILES model predictions (Vazdekis et al. 2015) for [α/Fe] = 0.0 and 0.3. GC systems of massive galaxies usually exhibit a metallicity gradient (e.g., Forbes & Remus 2018;Lee et al. 2020), in the sense that metal-rich GCs are more centrally concentrated than metal-poor GCs. In the left panel, the inner-halo GCs (R < 3 ) 2 show a large metallicity spread. The outer-halo GCs (R > 3 ) are divided into two groups, one of which appears to naturally follow the metallicity gradient while the other group has unexpectedly high metallicity given the radial distance. We mark the outer-halo, metal-rich GCs as red dots. The middle panel indicates that they are less enhanced in α elements than the other GCs with similar metallicities (at [Fe/H] = −0.5 ∼ 0.5). In the right panel, most of the outer-halo, metal-rich GCs are located in the youngerage region of CMRs. As we mentioned in Section 2.4, our sample could be contaminated by young GCs with intermediate/high metallicities. We therefore suspect that these metal-rich outliers are young GCs with a history different from that of most GCs in the galaxy. For example, they might be created from metal-enriched gas in other galaxies through a relatively prolonged star formation and later accreted to the Sombrero. For the metalpoor red outliers, on the other hand, no anomalous properties are found in terms of spatial distribution, radial velocities, metallicity gradient, and α-enhancement.

NONLINEAR METALLICITY-TO-COLOR CONVERSION AND COLOR BIMODALITY
This section focuses on the effect of the nonlinear GC CMR on the metallicity-to-color transformation. For the GC samples of NGC 5128 and NGC 4594 that are crossmatched between photometry and spectroscopy, we present the MD and CDs for each galaxy and analyze the morphologies of the distributions in terms of the degree of bimodality. We then demonstrate that the nonlinear metallicity-to-color conversion, combined with the observed GC MDs of NGC 5128 and NGC 4594, reproduces the observed bimodal CDs.
Figures 8 and 9 present the MD and CDs for NGC 5128 GCs and NGC 4594 GCs, respectively. The MD and CDs in the lower panels of each figure are for the identical crossmatched GC samples to those used for the CMR nonlinearity test performed in Section 3. The upper panels of each figure show, for reference, the CDs of photometrically selected old GC candidates, which have the same properties as those of the crossmatched GC samples in terms of the magnitude range and the radial coverage on the sky. Each histogram is smoothed using a nonparametric kernel density estimate with the bandwidth determined by the mean observational error. To mitigate spiky features in the histogram that arises from the small number statistics, we use the same bandwidth for both photometrically selected GCs (upper panel) and crossmatched GCs (lower panel). We note that, for the MD of NGC 4594 GCs, Alves-Brito et al. (2011) mentioned that their metallicity uncertainties might be underestimated, and the bandwidth is taken as twice of the mean observational error. The smoothed histograms are denoted as red lines for MDs and green lines for CDs. For NGC 5128 GCs, the MD can be described as a broad, unimodal distribution, but the CDs appear clearly bimodal, with the detailed shape being different depending on colors. The MD has a metal-rich peak but the CDs have prominent blue peaks. Given that metal-rich GCs are red, it is interesting that the shape of the MD is deformed dramatically after conversion into the CDs. For NGC 4954 GCs, the MD is unimodal and slightly tilted toward the metal-poor side, and the CDs exhibit more discernible bimodality and more prominent blue peaks and dips compared to the MD. In agreement with NGC 5128 GCs, the degree of CD bimodality of NGC 4594 GCs varies depending on color. This is consistent with the result for the GC systems of M84, M87, and NGC 1399 (Papers II and IV; . To evaluate the degree of bimodality of the MDs and CDs, we calculate the bimodality coefficient (BC; SAS Institute Inc. 1990), defined as where n refers to the number of samples, and m 3 and m 4 are respectively the skewness and excess kurtosis of the distribution corrected for sample bias. The BC exploits the skewness and the kurtosis of a distribution, which measure asymmetry and central clustering of the distribution, respectively. The basic idea is that a bimodal distribution generally has higher skewness and lower kurtosis than a simple normal distribution. The higher BC indicates stronger bimodality, and the lower BC indicates weaker bimodality. Readers are referred to Knapp (2007), Pfister et al. (2013), and references therein for further details. Table 4 presents the BCs for MD and CDs of GCs in NGC 5128 and NGC 4594. The BCs alike indicate that the CDs are more bimodal than the MD for both galaxies. This is consistent with the prediction of the nonlinear CMR theory, in which the inflection point along CMRs evokes or enhances bimodality by creating a dip in the CDs. The BCs also suggest that the degree of CD bimodality varies with color. Consistent with the visual impression above, both indices for U − B CDs of the two galaxies designate weaker bimodality compared to the other CDs, except for the BC value for NGC 4594 GCs 3 . We also carry out two statistical tests; Hartigans' dip test (Hartigan & Hartigan 1985) and the bimodality test using Gaussian Mixture Model (GMM) procedure. The former test is widely used for assessing the unimodality of a distribution by computing the maximum difference between the input distribution and the unimodal (uniform) distribution. The latter test assumes that the input distribution consists of multiple Gaussian components and constructs the GMMs that best describe the input distributions by maximum likelihood using the expectation-maximization algorithm. We use the diptest (Maechler 2021) and mclust (Scrucca et al. 2016) packages implemented in R (R Core Team 2022) to perform the dip test and the GMM test (we refer to it as GMM R ) respectively. In the GMM R test, we select the 'V' model that assumes the Gaussian components can have different variances. We also conduct a revised GMM test provided by Muratov & Gnedin (2010, we refer to it as GMM MG ), that employs three different approaches (the ratio of the likelihood, the peak separation, and the kurtosis) to get a more robust measure of bimodality.
All the tests provide the probability (i.e., p-value) for the null hypothesis that a distribution is unimodal. We note that the two tests occasionally give erroneous results. For example, the dip test hardly distinguishes between multimodal and unimodal distribution when the sample size is small and/or the significance level is set to be low (Kang & Noh 2019). The GMM MG test often falsely indicates bimodal when the true distribution is unimodal but skewed because of its sensitivity to the assumption of Gaussian modes (Muratov & Gnedin 2010). As these issues are beyond the scope of this paper, we focus on the relative comparison of the p-values for the MDs and CDs for each test. Table 5 provides the p-values from the dip test and the GMM tests for the MDs and CDs of NGC 5128 and NGC 4594 GCs. We remind that the smaller the pvalues, the stronger the bimodality. The dip test clearly shows that the MDs are closest to a unimodal distribution, and the CDs are relatively close to bimodal distribution for both galaxies. The GMM MG test gives three probabilities of a unimodal distribution based on maximum likelihood, peak separation, and kurtosis. They show that, with a few exceptions, the GC MDs for both galaxies are closer to unimodal distributions than their GC CDs. The GMM R test shows that, for NGC 5128, the p-value for MD is slightly larger, but a bimodal distribution is preferred over a unimodal distribution for all distributions. The p-values for NGC 4594 indicate that the MD is close to unimodal and the CDs are close to bimodal. Figure 10 presents Bayesian information criterion (BIC; Schwarz (1978)) and integrated complete data likelihood criterion (ICL; Biernacki, Celeux, & Govaert (2000)) values in the GMM R test as a function of the number of Gaussian components K. BIC and ICL are usually used to select the optimal model in a way that K-components GMM with the highest BIC/ICL value is preferred. BIC tends to select the number of mixture components, whereas ICL is more useful for selecting the number of discrete groups in the data (Biernacki, Celeux, & Govaert 2000). Given the properties of our data, it is more appropriate to focus on BIC rather than on ICL. Consistent with the p-values, BIC values are highest at K = 2 for all distributions except for the MD of NGC 4594 GCs. Note that, for the MD of NGC 5128, the difference in BIC values between K = 1 and K = 2 is ∼1.5. The difference between BIC values less than 2 is considered a "barely worth mentioning" difference (Raftery 1995). Thus, the GMM R test indicates that the MDs are relatively unimodal in both galaxies, as opposed to the CDs being bimodal.
The statistical tests corroborate that the distributions are significantly deformed when GC metallicities are converted to GC colors. The large variation in p-values among different colors (e.g., the p-value of the dip test for the NGC 4594 GC system) also indicates a significant change in their degree of bimodality due to the color-dependent shape of the CMR.
In Figures 11 and 12, we perform Monte Carlo simulations aiming at the CDs of NGC 5128 and NGC 4594 GCs to validate the nonlinear CMR effect on the conversion of MDs into CDs. The top row of the figures shows the observed GC CMRs as shown in Figures 5 and 6, overplotted with the YEPS model predictions. The observed MDs are rotated and given on the y-axis. The second rows show the observed color-magnitude diagrams of the crossmatched GCs (black dots) and the photometrically selected GCs (gray dots). Observational uncertainties as a function of M V are denoted by error bars on the left of each panel. The third and the fourth rows show the CDs of the crossmatched GCs and photometrically selected GCs, respectively (same as the Figures 8  and 9). The bottom row of the figures shows the simulated GC CDs. We generate 100,000 model GCs with metallicities derived from the Kernel density estimation of the observed MDs. We then transform the metallicity into colors using the CMRs predicted from the YEPS model, taking into account the observational uncertainties. For each color, we randomly select the same number of model GCs with the observed GCs 10,000 times and calculate the 1σ distribution ranges, which are shown by the shaded bands.
The most conspicuous feature of the observed CDs is the prominent blue peak. Our simulations with nonlinear CMRs can successfully reproduce the observed prominent blue peak that cannot be explained by the conventional linear CMR scheme given the unimodal MD with a metal-rich peak. Metal-poor GCs are projected through the steeper part of the CMR into a narrow color region, producing the sharp blue peak in the CD. By contrast, metal-rich GCs are projected through the shallower part of the CMR into a wider color region, resulting in an attenuated red peak in the CD. This even evokes a reversal of peak intensity for the NGC 5128 GCs.
The simulations show that projections on nonlinear CMRs with inflection points cause bimodal CDs from the observed MDs close to unimodal distributions. For the B − V , B − I, and V − I colors, the inflection point in the CMR causes a dip between blue and red peaks, which makes a discernible bimodal feature in the CDs. On the other hand, for the U − B color, the shape of the CD shows a skewed Gaussian without bimodality since there is no obvious inflection point in the U −B CMR. As a consequence, the simulation results for U − B, B − V , B − I, and V − I are consistent with the corresponding observations.

SUMMARY AND DISCUSSION
In this study, we have examined the nonlinearity of GC CMRs and its effect on the conversion of MD into CDs. We have used multiband photometric data combined with spectroscopic data for NGC 5128 and NGC 4594 GCs. We have demonstrated clear nonlinearity of the old GC CMRs for both galaxies, which is in good agreement with the model predictions.
Our visual inspection and statistical tests confirm the significant changes in the morphology of the CDs compared to the MD and even between the CDs. The most remarkable findings for NGC 5128 are that the metalrich peak of the GC MD is significantly reduced into a weak red bump in the CDs, and that strong blue peaks unexpectedly emerge in the (U -band-free) CDs. Also, a distinct dip that is not obvious in the MD appears in between the blue and red peaks in the CDs. Similarly, for NGC 4594 GCs, the blue peaks of the CDs are sharper compared to the corresponding metal-poor peak in the MD, and the dips are clearly visible in the CDs, making the CDs bimodal. We have demonstrated through Monte Carlo simulations that all of these features arise from the nonlinear conversion from the MDs to the CDs.
The CDs of photometrically selected old GCs are very similar to those of the GC samples crossmatched between photometry and spectroscopy in both galaxies (Figures 8 and 9). This implies that the MDs observed in this study are representative MDs of old GCs that are detectable at the corresponding galactocentric distance range. The shape of the GC MDs in both galaxies is broad and unimodal. This is a natural expectation in the Λ cold dark matter paradigm, where large galaxies such as NGC 5128 and NGC 4954 form through the hierarchical merging of numerous protogalaxies.
It is interesting to note that for NGC 5128, the shape of the GC MD is fairly consistent with that of the halo field stars' MD ( Figure 1 of Rejkuba et al. 2011) and closely matches that of the MD inferred from the inverse transformation of CDs using nonlinear CMRs ( Figure 6 of Yoon et al. 2011a). In the nonlinear CMR scenario, there is no discrepancy in MDs between GCs and halo field stars, and no further explanation for the excess of blue GCs is required. The readers are referred to Yoon et al. (2011a) for an in-depth discussion about solving the apparent discrepancy between GC MDs and halo stellar MDs.
Many studies on GC bimodality begin by dividing GCs into two groups in CDs-blue and red GCs. This is because, in CDs, GCs are distinctively divided into two groups by the dip. However, in our theory, the clear division in CDs do not necessarily indicate the presence of the two groups in an MD in a single galaxy. The crossmatched GC samples for NGC 5128 and NGC 4594 show that the CDs are bimodal but the corresponding MDs are closer to unimodal than bimodal. A crucial point of GC color bimodality is that it is a universal phenomenon commonly observed in massive early-type galaxies (Brodie & Strader 2006). The detailed forma-tion process of these galaxies varies from case to case, but there should be an essential common process to give rise to the phenomenon. In this respect, the nonlinear CMR scenario has its own strength in that the projection effect can simply solve the ubiquity of GC color bimodality. This study provides further direct observational evidence (see also Papers I-X) favoring the projection effect that significantly deforms MDs into various CD morphologies. Therefore, we conclude that the nonlinear CMR effect is the basic cause of the color bimodality phenomenon, and the detailed characteristics of CDs vary depending on the formation history of each galaxy.   Note-This table is available in its entirety in a machine-readable form in the online journal. A portion is shown here for guidance regarding its form and content. a The metallicities are transformed from the empirical metallicity ([M/H]) provided by Beasley et al. (2008).
b The photometric data are from Woodley et al. (2007) and corrected for Galactic extinction (see text). Note-This table is available in its entirety in a machine-readable form in the online journal. A portion is shown here for guidance regarding its form and content. a The metallicities are from Alves-Brito et al. (2011). Figure 1. Spatial distribution of GC data set used in this study for NGC 5128 (left) and NGC 4594 (right). For NGC 5128, U BV RI photometry is taken from Peng et al. (2004, black dots) and Woodley et al. (2007, blue dots), and spectroscopic metallicity is taken from Beasley et al. (2008, red dots). A total of 234 GCs are crossmatched between the photometry and the spectroscopy. For NGC 4594, U BV I photometry is taken from our observations. Gray dots are GC candidates selected using two-color diagrams (see Section 2.4). Spectroscopic metallicity is taken from Alves-Brito et al. (2011), who provide spectroscopic metallicity derived from three metallicity-sensitive indices (CH,Mgb,& Fe5270). We only use the metallicities derived from all three indices (red dots), except for metallicities from one or two indices (orange dots). A total of 108 GCs are crossmatched between the photometry and the spectroscopy.        We suspect that they are young GCs with a different formation and evolutionary history from most GCs in the galaxy, as they are less α-element enhanced than other GCs with similar metallicities in the middle panel and are in the younger-age region of CMR in the right panel. Figure 8. MD and CDs for GCs in NGC 5128. Each histogram is smoothed using a nonparametric kernel density estimate with the bandwidth corresponding to the mean observational error. The smoothed histograms are denoted as red lines for MDs and green lines for CDs. Upper panels: CDs for the GC candidates selected only from the photometric data set. Old GC candidates are selected from the photometric data of Peng et al. (2006) and Woodley et al. (2007) using the same criteria as the two-color diagrams, magnitude range, and radial extent used for the crossmatched GC sample. Lower panels: MD and CDs for the GCs crossmatched between photometry and spectroscopy. Figure 9. MD and CDs for GCs in NGC 4594. Each histogram is smoothed using a nonparametric kernel density estimate with the bandwidth corresponding to the mean observational error. The kernel bandwidth of MD is set to be twice the mean metallicity error since Alves-Brito et al. (2011) mentioned that the metallicity error they provide may be underestimated. The smoothed histograms are denoted as red lines for MDs and green lines for CDs. Upper panels: CDs for the GC candidates selected only from the photometric data set. Old GC candidates are selected from our photometric data using the same criteria as the two-color diagrams, magnitude range, and radial extent used for the crossmatched GC sample. Lower panels: MD and CDs for the GCs crossmatched between photometry and spectroscopy.  The models with the highest BIC/ICL values are preferred. Note that BIC tends to select the number of Gaussian components, whereas ICL is better suited for selecting discrete groups from data. Model (N = 100,000) Figure 11. Monte Carlo simulations for the CDs of NGC 5128 GCs. Top row: same as the CMRs in Figure 5. The error bars indicate the typical errors in [Fe/H] and colors. The observed MD is rotated and given on the y-axes, and the black solid line shows a kernel density estimate with a bandwidth set to the mean observational error. Second row: the observed colormagnitude diagrams for the GCs. Black and gray dots represent the crossmatched and photometrically selected GCs, respectively. Observational uncertainties are shown by error bars as a function of MV . Third row: observed CDs of the crossmatched GCs (see Figure 8). Fourth row: observed CDs of the photometrically selected GCs. Bottom row: simulated CDs of model GCs.
We generate 10 5 model GCs with metallicities derived from the kernel density estimate of the observed MDs and convert the metallicities into colors via the model CMRs presented in the top panel taking into account the observational errors. Given that the number of observed GCs is much smaller, we estimate the CD of 151 model GCs (solid line) and its 1σ distribution range (gray shaded bands) from the 10,000 bootstrap sampling.  Figure 6. The error bars indicate the typical errors in [Fe/H] and colors. We note that Alves-Brito et al. (2011) remark that their metallicity errors might be underestimated. The observed MD is rotated and given on the y-axes, and the black solid line shows a kernel density estimate with a bandwidth set to twice the mean observational error. Second row: the observed color-magnitude diagrams for the GCs. Black and gray dots represent the crossmatched and photometrically selected GCs, respectively. Observational uncertainties are shown by error bars as a function of MV . Third row: observed CDs of the crossmatched GCs (see Figure 9). Fourth row: observed CDs of the photometrically selected GCs. Bottom row: simulated CDs of model GCs. We generate 10 5 model GCs with metallicities derived from the kernel density estimate of the observed MDs and convert the metallicities into colors via the model CMRs presented in the top panel taking into account the observational errors. Given that the number of observed GCs is much smaller, we estimate the CD of 108 model GCs (solid line) and its 1-σ distribution range (gray shaded bands) from the 10,000 bootstrap sampling.