THE TWO-DIMENSIONAL PROJECTED SPATIAL DISTRIBUTION OF GLOBULAR CLUSTERS. I. METHOD AND APPLICATION TO NGC 4261

A considerable body of work has been accumulated in the past two decades on the properties of globular cluster (GC) populations in elliptical galaxies (see review by Brodie & Strader 2006). The studies available in the literature that use both Hubble Space Telescope (HST) and larger ground-based telescopes, have established the existence of red (metal rich) and blue (metal poor) subpopulations of GCs in most galaxies. The radial distributions of these populations differ, with the red/metal rich GCs being more centrally concentrated than blue GCs. These results suggest different formation scenarios for the two GC subpopulations, which may reflect the formation history of their parent galaxy. Brodie & Strader (2006), in particular, argue that blue GCs may have formed in low-mass dark matter halos in the early universe, while the red GCs were built in subsequent dissipative buildup of the parent galaxy. Other studies have investigated possible formation mechanisms for the families of low-metallicity and high-metallicity GCs using different approaches. These studies include semi-analytical models that use galaxy assembly history from cosmological simulations and observed scaling relations to estimate the amount and metallicity of cold gas available for star formation (Muratov & Gnedin 2010), physical models of the collapse leading to the formation of GCs based on high mass resolution simulations (Griffen et al. 2010), and models based on the observed galaxy mass–metallicity relation, the galaxy stellar mass function and theoretical merger rates (Tonini 2013). Other authors have explained the formation of low-metallicity GCs as starburst remnants of old dwarf galaxies that could have entered the halos of spiral galaxies (Elmegreen et al. 2012). While these works have explored multiple possible explanation of the peculiar observational traits of the GC populations, no agreement has yet been reached in the literature.

While the spectral and photometric properties of GC populations, as well as their radial distributions, have been explored in depth, relatively little work has addressed their two-dimensional (2D) distributions (e.g., NGC 4471 in Rhode & Zepf 2001, NGC 1399 in Dirsch et al. 2003; Bassino et al. 2006, NGC 3379, NGC 4406 and NGC 4594 in Rhode & Zepf 2004, NGC 4636 in Dirsch et al. 2005, multiple galaxies in Hargis & Rhode 2012, NGC 3585 and NGC 5812 in Lane et al. 2013). This type of work requires us to extract the entire GC population with full or near-full coverage of the parent galaxy. Augmented by kinematics of the GC systems and subsystems (Strader et al. 2011 for M87; Blom et al. 2012 for NGC 4365), these studies are bringing forth a picture of complex and diverse GC populations in elliptical galaxies, consistent with a history of gravitational interactions and merging with neighboring galaxies.

Here we study the 2D GC distribution of NGC 4261, for which deep observations of the innermost region of the galaxy, obtained with HST, are reported by Bonfini et al. (2012, hereinafter B+12). NGC 4261 first came to attention because of the inhomogeneous spatial distribution of its low-mass X-ray binaries (LMXBs), which suggested a possible minor merger event (Zezas et al. 2003). Since a number of these LMXBs were associated with GCs (Giordano et al. 2005), the entire population of 718 GCs was extracted and studied with HST, revealing an azimuthal asymmetry in its distribution, which was attributed to past minor merging or interaction (B+12). Ferrarese et al. (1996) also reported boxy isophotes which support the hypothesis of recent gravitational interaction experienced by NGC 4261.

Although the results of B+12 are convincing, our more advanced analysis methods provide a clearer picture of the asymmetry and its statistical significance. Here we report the re-analysis of the 2D distribution of this GC population aimed at quantifying the reality and shape of the spatial features in the 2D GC distribution. We have used the K-Nearest Neighbor method (KNN; Dressler 1980) to identify the 2D features, supplemented by Monte-Carlo simulations to test their statistical significance. In Section 2 the data used in this paper are described. The method and the results of the its application to NGC 4261 are discussed in Section 3. Our findings are discussed in Section 4 and summarized in Section 5.

We have used the B+12 catalog of GC positions and properties, which lists 718 GCs, within the D₂₅ ellipse of NGC 4261 (de Vaucouleurs et al. 1991). From this sample, we have extracted color and magnitude based subsamples. Figure 1 (left) shows the V − I histogram where blue (V − I < 1.15) and red (V − I ⩾ 1.15) GCs are separated following B+12. These authors found that V − I = 1.15 is the color corresponding to the 50% probability of the GC to belong to either the red or the blue subpopulations, assuming a two-Gaussian model of the color distribution. However, as discussed in B+12, the GC color distribution does not show the clear bimodality typical of the GC color distribution in other early-type galaxies (Brodie & Strader 2006; Peng et al. 2006). Figure 1 (right) shows the histogram of I magnitudes; we arbitrarily define high luminosity (high L) GCs those with I < 23 mag, and low luminosity (low L) GCs those with I ⩾ 23 mag. The I < 23 value used to separate low-L from high-L GCs was set to obtain equipopulated classes of sources (see Table 1). The density and residual maps obtained for the luminosity classes defined using I < 23 are described in Section 3.1. Figure 2 shows the spatial distributions of GC positions in the plane of the sky, where the azimuthal asymmetry in the 2D projected GC distribution is evident. Table 1 lists the number of GCs in each of the main samples used in the following analysis. Changing the color and luminosity boundaries (shaded regions in the plots in Figure 1) does not affect our results. Details can be found in Section 3.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have excluded from our analysis the central circular region with r < 042 where incompleteness in the detection of the GCs is substantial (B+12). We have also excluded the undersampled regions outside the D₂₅ isophote (de Vaucouleurs et al. 1991).

For each of the GCs samples listed in Table 1, we have determined the 2D spatial distributions by applying the KNN density estimator (Dressler 1980). This density is based on the local distribution of GCs, i.e., on the distances of the closest GCs. For each knot of a regular grid covering the region where the density is to be determined, we measured the distance of the Kth nearest neighbor GC (D_K) from the position of the knot. We used this approach in order to have density estimates even in the low density regions of the GC distribution. The point-density is estimated as:

$$\begin{equation} D_{K}=\frac{K}{V_{D}(D_{K})} \end{equation} \tag{ 1 }$$

where K is the index of the nearest neighbor used to calculate the density; for example, for K = 5 only the five GCs nearest to the grid knot are used. V_D is the volume of the region within the distance D_K of the Kth nearest neighbor from the point where the density is evaluated. In the case of 2D spatial density, the volume $V_{D}=\pi \cdot d_{K}^2$ is equivalent to the area of the circle with radius equal to the distance of the Kth nearest neighbor. The assumption behind this method is that the density is locally constant. The uncertainty on the KNN density scales with the square root of K, so that the relative fractional error is

$$\begin{equation} \frac{\sigma _{D(K)}}{D(K)}=\frac{1}{\sqrt{K}}. \end{equation} \tag{ 2 }$$

The fractional accuracy of the method increases with increasing K at the expense of the spatial resolution.

For each subsample, we have determined the 2D GC surface density by using values of K ranging from 2 to 10. The density values have been calculated on a regular grid with spacing ∼00018 (∼65) and ∼00015 (54) along the R.A. and declination respectively, in order to have the same number of knots covering the whole region of the sky occupied by the GC distribution along the two coordinate axes. Different spacings along the R.A. and declination axes yield similar results to those described in the following. The density in the pixels overlapping the boundaries of the area covered by the observations has been weighted according to the fraction the area of the pixel located within the observed region.

To assess the significance of features suggested by the KNN density maps (hereafter "observed" density maps), we performed Monte Carlo experiments, by creating random samples of GCs, each conforming to the observed radial density distribution of the relevant observed samples. The simulated random samples contain the same number of GCs as the corresponding observed samples. The simulated radial positions were randomly drawn from the histogram of the observed radial distribution of GCs integrated over elliptical annuli with the same eccentricity and position angle of the D₂₅ isophote (e = 0.45, de Vaucouleurs et al. 1991). Independently, the azimuthal distributions were extracted from a uniformly random distribution between 0° and 360°. With this prescription, we do not need to determine a 2D model of the observed GC distribution in NGC 4261. We introduced a geometrical correction on the expected number of GCs to take into account the eccentricity of the galaxy.

We also tried using a histogram generated from the best fit power law of the observed radial density profile to seed our simulations. This approach has the advantage of smoothing local density fluctuations and therefore, in principle, makes us more sensitive to small-scale annular features. However, it introduces a model which may not represent the real radial distribution of GCs in NGC 4261. Both approaches gave consistent results. In the following we will describe only the results obtained using the observed histogram to seed the simulations.

We performed 5000 simulations for each of the samples in Table 1 and constructed from each simulation the 2D density maps using the same KNN method applied to the observed data, for values of K ranging from 2 to 10. From these maps we derived average pixel-by-pixel simulated density maps for each sample under consideration (in the remainder of the paper, we will use "pixel" to indicate each cell in the R.A. versus declination grid used to evaluate the 2D density associated to the GCs spatial distribution). Figure 3 compares the observed and simulated mean radial density profiles for the color-based samples.

Zoom In Zoom Out Reset image size

Since our purpose is to characterize the deviations from the average radial trends of the GC distributions, our next step was to calculate the maps of residuals R by subtracting pixel-by-pixel the mean simulated density map 〈S〉 from the corresponding observed map O. These residual maps were also normalized by the value of the density in the average simulated map 〈S〉. The residual R_i of the ith pixel of the map is thus defined as

$$\begin{equation} R_{i}=\frac{(O_{i}-\langle S\rangle _{i})}{\langle S\rangle _{i}}. \end{equation} \tag{ 3 }$$

We note that the pixel-by-pixel distributions of the simulated KNN densities are well approximated by Gaussians, simplifying the calculation of the statistical significance of the observed over densities. A typical example of the pixel density distributions is shown in Figure 4.

Zoom In Zoom Out Reset image size

To evaluate statistical significance of the observed residuals, for each set of simulated density maps we calculated the fraction of pixels with values above the 90th percentile of the densities in the observed maps (the "extreme" pixels).

Since the extreme pixels in the observed density and residual maps tend to be clearly spatially correlated, we also evaluated the fraction of simulations with at least one group of contiguous extreme pixels as large as the observed. For each simulation, we counted the number of groups of extreme contiguous pixels with area (measured in pixels) equal to the area of the region within the density contours corresponding to the 90th percentile observed density threshold. Since we did not impose any specific geometry to the groups of contiguous extreme pixels in the simulated density maps, these fractions represent upper limits to the fraction of extreme contiguous pixels expected for a given spatial distribution of residuals.

The results, compiled in Table 2, show that for small values of K, simulated density maps have a total fraction of extreme pixels larger than the observed. However, this fraction decreases with increasing K, and it becomes negligible for K ⩾ 5. The percentage of contiguous extreme pixels (in parenthesis in Table 2) does not exceed 20% even for K = 2 and rapidly decreases to zero at K = 5. Therefore, the probability of a random distribution of high-density pixels with the observed spatial distribution is very low even for low K and becomes null with increasing K. Based on Table 2, in the following we will only discuss our results for K ⩾ 5. Figure 5 shows the histogram of the simulated number of extreme pixels compared to the observed for K = 9 for the entire sample, and the corresponding histograms for spatially correlated extreme pixels. Together with Table 2, these comparisons shows that the degree of correlated features observed (discussed in detail below) in the data cannot be randomly produced.

Zoom In Zoom Out Reset image size

3.1. Density and Residual Maps

Figure 6 shows the 2D KNN density maps for the spatial distribution of the entire GC sample for K ranging from 5 to 10. In all these maps two high-density structures emerge in the SW and NE quadrants, on opposite sides of the galaxy major axis. The SW density enhancement resembles a spiral arm and is composed of two distinct sections: a short elongated structure aligned along the radial direction and a second longer feature spanning azimuthally more than 45°. The NE high-density region has an elongated shape and occupies azimuthally ∼100°. These over-densities roughly correspond to the local enhancements of the GC distribution shown in Figure 9 of the B+12 paper, which was obtained by applying adaptive kernel smoothing and Voronoi Percolation–Tessellation. This figure shows that maps with K = {7, 8, 9} are similar, while the K = 10 maps smooth over some of the spatial features. In the remainder of the paper, we will discuss only the density and residual maps obtained for K = 9 because these maps are the most statistically significant (as discussed in Section 3 and shown in Table 2).

Zoom In Zoom Out Reset image size

We show in Figure 7 the K = 9 residual map, which highlights the existence of the two main over-densities in the SE and NW quadrants already observed in the density maps and discussed above, together with similarly significant under-densities in the SE region and the western region of the galaxy. The latter is partially due to the incomplete HST coverage, see Figure 1 in B+12 for details. Figure 8 maps the single pixels with significance larger than 1σ, 2σ and 3σ. In particular the SW features stand out, since multiple spatially correlated nearby pixels with significant positive residuals are located in this region. Although the NE over-density is weaker, the spatial correlation of the pixels with positive residuals in this case cannot be due to chance, as shown in Table 2. The number of GCs located within the pixels associated to positive residuals is 175 (∼24%), 86 (∼12%), and 60 (∼8%) for 1σ, 2σ, and 3σ respectively. Table 3 shows the fractions of GCs located within the statistically significant over-densities for all GCs classes used in this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Number of GCs Located in Over-densities Regions with Significance Larger than 1σ, 2σ and 3σ Respectively for all Classes of GCs Used in this Paper

	All GCs			Red GCs			Blue GCs
	1σ	2σ	3σ	1σ	2σ	3σ	1σ	2σ	3σ
K = 5	241(33.5%)	137(19.1%)	78(11.0%)	203(66.3%)	109(35.6%)	79(25.8%)	191(46.4%)	100(24.3%)	67(16.3%)
K = 6	241(33.5%)	133(18.5%)	80(11.1%)	217(70.9%)	109(35.6%)	79(25.8%)	219(53.2%)	117(28.4%)	60(14.6%)
K = 7	247(34.4%)	125(17.4%)	88(12.3%)	210(68.6%)	122(38.9%)	80(26.1%)	228(55.3%)	125(30.3%)	70(17.0%)
K = 8	247(34.4%)	133(18.5%)	79(11.0%)	218(71.2%)	141(46.1%)	93(30.4%)	225(54.6%)	112(27.2%)	76(18.4%)
K = 9	252(35.1%)	137(19.1%)	90(12.5%)	211(70.0%)	128(41.8%)	91(29.7%)	230(55.8%)	126(30.6%)	63(15.3%)
K = 10	274(38.1%)	140(19.5%)	93(13.0%)	220(71.9%)	135(44.1%)	78(25.5%)	224(54.4%)	143(34.7%)	77(18.7%)
	High-L GCs			Low-L GCs
	1σ	2σ	3σ	1σ	2σ	3σ
K = 5	194(61.4%)	104(33.0%)	58(18.3%)	200(62.9%)	103(32.4%)	73(23.0%)
K = 6	210(66.5%)	104(33.0%)	63(19.9%)	208(65.4%)	123(38.7%)	86(27.0%)
K = 7	216(68.4%)	107(33.9%)	55(17.4%)	221(69.5%)	138(43.4%)	96(30.2%)
K = 8	232(73.4%)	113(35.8%)	64(20.3%)	239(75.1%)	146(45.9%)	80(25.2%)
K = 9	231(73.1%)	112(35.4%)	67(21.2%)	262(82.4%)	142(44.7%)	98(30.8%)
K = 10	220(69.6%)	127(40.2%)	66(20.9%)	245(77.1%)	138(43.4%)	93(29.2%)

Note. In parenthesis, the percentage relative to the total number of sources in each class, as shown in Table 1.

Download table as:  ASCIITypeset image

Tal et al. (2009) report the presence of two low surface brightness structures, a NW tidal arm and a faint SE "fan" in NGC 4261. These structures are located within the HST observations used by B+12 to extract the catalog of GCs, but within the central avoidance region not used in this study (see Section 2). We have compared the spatial distribution of this faint optical features with the large scale shape of the GCs distribution in the region of the galaxy closest to the avoidance region, but we have not noticed any resemblance or clear correlation between the two different structures.

The K = 9 maps of red and blue GC classes (Section 2), show significant differences in the over-densities, although the overall residuals follow those of the entire sample (Figure 9). In particular, a stronger concentration of ⩾3σ over-dense pixels can be seen in the north of the blue GC residual map. These results do not depend on the value of the B+12 color threshold used to define blue and red subpopulations, as long as the threshold is chosen within the interval V − I = [1, 1.3]. The density and residual maps of the two color classes using ten regularly spaced threshold values in the above interval are self consistent. For values outside this interval, the significance of the residual map for one of the classes degrades rapidly because of the small number of GCs.

Zoom In Zoom Out Reset image size

The observed radial density profiles of the GC distributions in wedges aligned along the axes of the galaxy (Figure 10) indicate that GCs have different radial distributions along the directions of the major and minor axes of the galaxy. The density profile of the entire sample along the major axis is flatter than the average density profile for radii r ⩾ 12, while the minor axis density profile is significantly steeper (Figure 10, left), after correcting for geometrical effects. Both red and blue GC profiles are more extended along the major axis, compared with the minor axis distributions. This difference is more striking for the blue GCs, which have a flat major axis profile while the minor axis profile first increases at small radii, and then plummets toward D₂₅.

Zoom In Zoom Out Reset image size

The K = 9 density maps for the high and low luminosity subsamples (Table 1) are shown in Figure 12. The density distributions show that the high-L GCs are more centrally concentrated than the low-L ones. This is also reflected by the azimuthally averaged radial profiles (Figure 11). This effect could be at least in part due to incompleteness that will be more pronounced in the inner radii where the galaxy stellar light is more intense. Although there is a small radial dependence, this effect is minimized by only using GCs with I ⩽ 24 for which, according to B+12, a completeness of 75% is achieved. However, there are azimuthal differences at the same radius that cannot be due to incompleteness effects. Figure 13 (lower panels) shows the position of the pixels with residuals larger than 1σ obtained using K = 9 for high-L and low-L GCs (the same plots for blue and red GCs are shown in the upper panels of Figure 13 for reference). These plots also show significant differences: there is an excess of high-L GCs to the W of the major axis in the N, while the low-L GCs present an excess to the opposite SE side. The same conclusions about the spatial distributions of low and high luminosity classes of GCs can be drawn using different values of the threshold magnitude within the interval I = [22, 23.5], using 10 regularly spaced magnitude thresholds. The high-L GCs are, in each case, more centrally concentrated than low-L GCs, and high-L residuals show a unique significant positive density enhancement in the NW region of the galaxy. Other values of the I magnitude have not been considered because they would generate luminosity classes too unbalanced to correctly estimate the significance of the results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Color versus Spatial Distribution of GCs

Blom et al. (2012) found significant effects in the color versus galactocentric distance distribution of the GC system of the galaxy NGC 4365. We have investigated if similar effects are present in NGC 4261. We have produced the density maps of the distribution of GCs in the radial distance versus V − I color with the KNN method for K = {5, 6, 7, 8, 9, 10}. The map for K = 9 is shown in Figure 14. The clumpy distribution of over-densities does not show clear separation between different classes of color except for a tendency of red GCs to be more centrally concentrated than blue GCs, in agreement with the density and residual maps of the two separate color classes of GCs (Figure 9). In particular, at large radial distances, the only interesting density enhancement of Figure 14 is located in r ∼ 175 and V − I ∼ 1.12. This feature corresponds to the sum of the two significant density enhancements in the spatial density map of the blue GCs class located along the major axis of the galaxy in the NW quadrant and in the SW quadrant (see upper left plot in Figure 9). A more detailed representation of the GCs color distribution in NGC 4261 as a function of the position in the galaxy can be seen in Figure 15, where each pixel is colored according to the average color of the GCs placed within the pixel. The size of the pixels is proportional to the density evaluated with the KNN method for K = 9 and the density contours are derived from the same density map. Figure 15 confirms that red and blue GCs tends to follow the same overall spatial patterns, except for the different small differences visible in the upper panels in Figure 9, obtained using the V − I = 1.18 color threshold value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our study of the 2D projected distribution of GCs in the elliptical galaxy NGC 4261 by means of KNN density maps has confirmed at high statistical significance the presence of the 2D large-scale anisotropy reported by B+12. The morphology of this anisotropy suggests a large-scale spiral-like over-density which originates along the major axis and follows the D₂₅ elliptical isophote in the periphery of the galaxy. Alternatively, there may be a system of broken shells enclosing the major axis in both directions roughly along the D₂₅ ellipse, with local density enhancements especially along the northern segment of the galaxy major axis. We find the same overall "grand pattern" when the sample is subdivided by either color or luminosity, although there are specific differences between the subsamples. In particular, while the southern feature is seen in both red and blue GCs, the northern anisotropy shows some segregation of red and blue GCs. Similarly the morphology of the residuals differs in the NW (see Section 3.1) suggesting a lack of high-L GCs in the inner major axis. Overall there are more red GCs at smaller radial distances (see Figure 3). However, we note that these effects may be correlated, since the high-L sample has a slightly larger fraction (∼56%) of blue GCs (see Table 1).

At small radial distances, the density enhancements of the GC distribution are aligned along the minor axis of the galaxy. The overall distribution is driven by high-L GCs which are located along the minor axis direction while low-L GCs do follow the broken-shell/spiral arm features more closely. These results differ significantly from recent results from Wang et al. (2013), which show that GCs of the early-type galaxies in the ACS Virgo Cluster survey tend to be aligned along the major axis in galaxies with visible elongation and intermediate to high luminosity. The GC system in NGC 4261 shows that the geometry of the GC distribution observed in the Wang et al. (2013) sample is not universal, even though the different conclusions based on the results provided by our method may depend on the rich sample of NGC 4261 GCs used, which, nonetheless, does not cover the whole GC system.

As already discussed in B+12, these large-scale features suggest global shaping events, however their nature is not clearly understood. Tidal tails and shell-like ripples of the stellar surface brightness have been reported in several elliptical galaxies, and have been modeled in terms of mergers and interactions with satellites (e.g., Schweizer 1980; Quinn 1984; Hernquist & Quinn 1988; Hernquist & Quinn 1989). However, these simulations did not address GC systems.

In NGC 4261, no significant sign of perturbation of the stellar surface brightness of the galaxy (except for two faint central structures noticed by Tal et al. 2009), is observed, indicating that a major merging event has not likely occurred in the last ∼1 Gyr. However, Ferrarese et al. (1996) observed boxy isophotes which are consistent with past merger. As discussed in B+12, a major merger scenario may be an alternative for the global 2D anisotropy of the GC distribution only if the relaxation timescale for the GCs is larger than that of the stars. Another possibility in the major merger scenario is that the GC over-densities may be related with GC formed during the merger in shocked regions and tidal tails, such as the evolved remnant of the young massive cluster populations detected in the Antennae (Whitmore et al. 2010) and NGC 7252 (Bastian et al. 2013).

The peculiar spatial distribution of GCs in NGC 4261 could also partially be influenced the GC systems of the smaller galaxies visible in the NGC 4261 field and observed by Tal et al. (2009). In particular, two dwarf galaxies are visible in the HST images of NGC 4261 GCs. The first, with B = 16.2 mag (Binggeli et al. 1985), is located in the NW corner of the HST mosaic (R.A. = 12: 19: 17.6, Decl. = +5: 52: 39) and is too far from NGC 4261 to contribute to the NGC 4261 GCs system. The second (B = 15.3) is close to the SW spiral-arm feature visible in the density map of the NGC 4261 GC system (R.A. = 12: 19: 23, Decl. = +5: 47: 51). Similar galaxies have compact GC systems (see Dirsch et al. 2005). Thus, if this galaxy has an GC system overlapping with the NGC 4261 GC population, this would only influence a small region of the spiral-like over-density associated to the NGC 4261 GC spatial distribution.

A third, brighter, galaxy (B = 13.8) (R.A. = 12: 19: 35, Decl. = +05: 50: 47), could in principle have a GC population spatially overlapping the NGC 4261 GC system. Nonetheless, we notice that this galaxy is located on the E side of NGC 4261 and no significant density enhancements is visible in the E outskirt of NGC 4261 (Figure 6). A feature resembling a broken shell is apparent in the E quadrants of the NGC 4261 closer to the center of the galaxy. Assuming a power-law radial density profile for the GCs system of the companion galaxy, the result of its contribution to the NGC 4261 GC density map would be an enhancement of the density decreasing from the outskirts to the center of NGC 4261. For this reason, the contamination from the GC system of the neighboring galaxy could only reduce the significance of the observed over-density in the E side of NGC 4261, and would not change qualitatively the results of our analysis.

We cannot exclude that present and past gravitational interactions with these neighboring galaxies may have affected the GC system of NGC 4261, but the molding of these effects is beyond the scope of this paper.

Recent smoothed particle hydrodynamics simulations of the formation of GCs in merging and interactive galaxies suggest that newly formed metal-rich GCs tend to concentrate at the center of the merger remnant elliptical galaxy while metal-poor GCs are distributed in the outer parts due to strong angular momentum transfer (see Bekki et al. 2002). However, our residual maps, while showing a few differences in the over-densities of red and blue GCs, demonstrate that the overall patterns are similar, as already noted by B+12. This result argues against a "rejuvenation" of the GC populations that would result in an increase of the higher metallicity red population. B+12 also noted that a minor merger scenario cannot be excluded, but the incoming galaxy should be particularly GC-rich, and that a displacement of the GC population by a galaxy fly-by may be a possible explanation.

The characterization of the distribution of GCs in the radial distance versus V − I space (Section 3.2) demonstrates the lack of a clear bimodality in the color distribution of GCs over the whole interval of radial distances considered in our analysis. Nonetheless, the expected differences of the spatial distributions of red and blue GCs are evident in our analysis, regardless of the specific color value used to separate the two classes (Section 3.1).

While our paper cannot solve the nature of the perturbations in the spatial distribution of GCs in NGC 4261, we have contributed detailed 2D maps of the morphology of the asymmetry at different spatial scales. The morphology of the perturbations suggests a stream or an orbiting GC system, but clearly this could also be a projection effect. While difficult to obtain, a large kinematic sample of GCs could be used to compare the observed spatial distribution of GCs with the results from simulations of the formation of elliptical galaxies, providing a way to constrain the dynamical evolution of the host galaxy (see, for example, Bekki et al. 2005).

We have developed a new approach to the study of the 2D distribution of GCs in elliptical galaxies. The method is an implementation of the KNN method presented in Dressler (1980), supplemented by Monte Carlo simulations to establish the statistical significance of the results. We have applied this method to NGC 4261, a "test galaxy" where significant 2D anisotropy in the GC distribution has been reported (B+12). We confirm that the 2D distribution of GC is not azimuthally isotropic. Moreover, we demonstrate that the 2D distribution departures from the average GC radial distribution results in highly significant spiral-like or broken shell features. While the southern feature is seen in both red and blue GC subsamples, the northern anisotropy shows some evidences of segregation of red and blue GCs. The good characterization of the 2D distribution of GCs that we have achieved provides a new benchmark for future simulations of galaxy merging evolution.

R.D'.A. gratefully acknowledges the financial support of the U.S. Virtual Astronomical Observatory, which is sponsored by the National Science Foundation and the National Aeronautics and Space Administration. T.F. acknowledges support from the CfA and the ITC prize fellowship programs. This work was partially supported by the Chandra X-Ray Center (CXC), which is operated by the Smithsonian Astrophysical Observatory (SAO) under NASA contract NAS8-03060.