Substructure in the Globular Cluster Populations of the Virgo Cluster Elliptical Galaxies M84 and M86

According to the current paradigm, massive galaxies formed hierarchically, through the continuous merging and accretion of smaller protogalaxies. In this picture, massive structures like giant galaxies started to form in the early universe when small fluctuations in the matter density distribution began to collapse and accrete surrounding matter. These structures continued to merge and build up larger structures over the course of cosmic time. Simulations that involve this type of hierarchical merging and matter accretion within the context of a ΛCDM universe have been fairly successful at recreating the observed overall distribution of matter (e.g., Springel et al. 2005; Vogelsberger et al. 2014; Schaye et al. 2015) and producing galaxies with realistic properties that are generally well matched to the properties observed today (e.g., Springel et al. 2008; Hopkins et al. 2014, 2018).

Much of the recent work on this topic has focused on gathering observational evidence for hierarchical growth through the identification of stellar streams and tidal interactions, which are thought to be the hallmark of low-mass galaxies having been accreted by more massive galaxies. The Sagittarius dwarf galaxy (Ibata et al. 1994) and the stellar streams associated with it (Belokurov et al. 2006) are some of the most well-studied examples of this type of hierarchical accretion in the Milky Way. Recent data from the Gaia mission (Gaia Collaboration et al. 2016) and other large-scale surveys, as well as analyses that combine kinematic data with information about chemical abundances and/or specific stellar populations, have led to discoveries of additional substructures in the Galaxy and brought renewed attention to this topic (e.g., Belokurov et al. 2018; Deason et al. 2018; Helmi et al. 2018; Myeong et al. 2018). The Andromeda galaxy (M31) has been studied in detail as well. Structures formed from spatial overdensities of the red giant branch (RGB) population, e.g., the Giant Stellar Stream, have been found in the outer halo of M31 by using resolved RGB star maps (Ibata et al. 2001; Ferguson et al. 2002). Using the surface density of resolved RGB stars, it is possible to infer the surface brightness distribution of the underlying light; structures revealed by this technique, such as faint streams and clumps of RGB stars, suggest a rich history of accretion events and tidal interactions for M31. Deep observations of galaxies outside the Local Group have also revealed streams, shells, and satellites using resolved RGB stars and unresolved low surface brightness (LSB) features (e.g., Janowiecki et al. 2010; Crnojevic et al. 2016; Mihos et al. 2017). However, in galaxies well beyond the Local Group, RGB stars are unresolved and a different tracer of hierarchical assembly processes and tidal interactions is required. Globular clusters (GCs) have a number of properties that make them well suited for this task.

GCs are compact stellar systems with typical masses ranging from ∼10⁴ to ∼10⁶ M_⊙. They are highly luminous objects and therefore visible up to distances of hundreds of megaparsecs (e.g., Blakeslee 1999; Mieske et al. 2004). They are ubiquitous across all galaxy types, with giant galaxies typically hosting hundreds or thousands of GCs (e.g., Rhode 2012; Young 2016). Furthermore, GC metallicities have been shown to contain information about hierarchical galaxy assembly processes. The Milky Way has two populations of GCs with different mean metallicities, and the abundances and other properties of these two populations have been used to probe the Galaxy's past history (Zinn 1985; Armandroff & Zinn 1988). Most other massive galaxies have GC systems that show multiple peaks in their color distributions (e.g., Brodie & Strader 2006 and references therein). For integrated colors of stellar populations older than ∼2−6 Gyr, these color differences can indicate a difference in metallicity, with bluer colors corresponding to lower metallicities and redder colors corresponding to higher metallicities (e.g., Young 2016). The presence of multiple color peaks in the GC systems of some galaxies is commonly interpreted as evidence that such galaxies have undergone two or more major epochs of star formation, possibly triggered by mergers or accretion events (e.g., Ashman & Zepf 1992; Zepf & Ashman 1993; Li & Gnedin 2014; El-Badry et al. 2019). Moreover, the spatial distributions of the GC subpopulations have been shown to differ in many galaxies, with blue/metal-poor GCs typically covering a larger radial extent than their red/metal-rich counterparts (e.g., Brodie & Strader 2006, and references therein). While this spatial difference may have multiple causes, one hypothesis is that a large proportion of the metal-poor GCs in massive galaxies were deposited through accretion events (Côté 1999; Forbes & Bridges 2010).

GCs have also been shown to trace substructure within and around galaxies. For example, Lim et al. (2017) performed a study of the GCs in the early-type "shell" galaxy NGC 474 and found a significant correlation between the shell and stream structures and the GC spatial positions. Furthermore, previous studies have reported similar results for the M31 GC system (see Mackey et al. 2010; Ferguson & Mackey 2016, and references therein), finding a strong correlation between the spatial position of GCs within the system and stellar streams.

Since not all galaxies have readily visible substructure to compare to GC positions, another approach is to identify deviations in the spatial positions of GCs from the expected symmetric distribution of an unperturbed GC system. Bonfini et al. (2012) found an azimuthal asymmetry along the northeast−southwest direction in the GC system of the elliptical galaxy NGC 4261. D'Abrusco et al. (2013) recovered this same feature in NGC 4261 using the K-Nearest Neighbors algorithm to estimate the surface density at various points in the field. In doing so, D'Abrusco et al. (2013) probed for radial and azimuthal asymmetries within the NGC 4261 GC system and found that the overdensity forms a broken shell or spiral-like pattern. D'Abrusco et al. (2015) probed the cores of the GC systems of the 10 brightest Virgo galaxies for anomalies and found substructures of varying complexity and sizes (∼05 to several arcminutes in length) in all of the observed systems. When studying the Fornax galaxy cluster, D'Abrusco et al. (2016) found an overdensity within the NGC 1399 GC system. This overdensity is dominated by the blue GCs and stretches from east to west, connecting NGC 1399's GC system to the GC systems of neighboring galaxies. Furthermore, the overdensity of blue GCs in Fornax has been shown to be associated with a region of intracluster light within the core of the cluster (Iodice et al. 2017). Durrell et al. (2014) and Powalka et al. (2018) identified GC population features within M87 (NGC 4486) and in the larger Virgo Cluster environment using kernel density estimation (KDE) surface density maps. Durrell et al. (2014) found that the spatial extent of the metal-poor GC populations extends much farther (∼400 kpc from the galaxies) than metal-rich populations around M87 and M84 (NGC 4374), implying that intracluster GC populations are primarily composed of metal-poor GCs. Dividing their sample into three age bins, Powalka et al. (2018) detected a spatial overdensity of young GCs to the south of M87 that is not present in the intermediate-age and oldest GC populations. This overdensity was detected at the 5σ level, and no host galaxy was detected within the overdensity region. Madrid et al. (2018) created a surface density map of the Coma Cluster that highlighted evidence of several ongoing galaxy interactions. Such studies make it clear that GCs can be useful for discovering substructure in galaxy halos and exploring the history and evolution of galaxies and their stellar populations.

In this paper, we make use of existing wide-field CCD imaging data to search for evidence of substructure and galaxy interactions in the GC system of the Virgo Cluster elliptical galaxy M86 (NGC 4406) and the neighboring elliptical galaxy M84 (NGC 4374). Previous efforts have been made to study the substructure present within the central few arcminutes of the M86 GC system (D'Abrusco et al. 2015) and across the Virgo Cluster as a whole (Durrell et al. 2014). The spatial coverage of our data falls in between these two past works; our images include the two massive elliptical galaxies, along with several other low-mass galaxies located at similar distances. Our study is intended to explore whether the GC systems of the galaxies in the M86 field show any evidence for interactions between the galaxies and their cluster populations. Theoretical studies predict that environmental effects due to a galaxy cluster tidal field and galaxy interactions can significantly affect the morphological and kinematic properties of GC systems (e.g., Ramos et al. 2015), strip some GCs from their original host galaxies (e.g Forte et al. 1982; Muzzio et al. 1984; Ramos-Almendares et al. 2018, 2020), and produce intracluster GC populations (e.g., Yahagi & Bekki 2005; Bekki & Yahagi 2006, 2009). The goal here is to carefully search for signatures of these processes in our observational data.

The organization of this paper is as follows. Section 2 describes our data set and the steps we carried out to produce a list of positions and photometric measurements of GC candidates in the M86/M84 field. Section 3 describes both our analysis methods and our results. In this section we report possible substructures within the M86/M84 field, variations in these substructures among the red and blue GC populations, and our findings when studying the kinematics of GCs within the field. The last section of the paper presents a summary of this work.

To investigate substructure in the GC systems of M86 and M84 and the surrounding region, we used the same imaging data that were originally presented in Rhode & Zepf (2004). The goal of the Rhode & Zepf (2004) study was to quantify the global properties—i.e., the total number of GCs, along with the specific frequency, radial distribution, color distribution, and color gradient—of the GC system of M86 and a few other giant early-type galaxies. The images of the M86 field were obtained in 1999 March with the Mosaic imaging camera on the Mayall 4 m telescope at Kitt Peak National Observatory,¹ as part of a wide-field multicolor optical survey of the GC systems of a sample of several spiral, S0, and elliptical galaxies (e.g., Rhode & Zepf 2001, 2003, 2004; Rhode et al. 2007). A ∼36' × 36' area around M86 was imaged for a total integration time of 3900 s in B, 2700 s in V, and 2100 s in R. A detailed description of how the Mosaic images were processed and calibrated is given in the original Rhode & Zepf (2004) paper. The pixel scale of the final science images we used is 026 pixel⁻¹, and the FWHM of the point-spread function (PSF) is ∼098, 11, and 12 in the B, V, and R images, respectively.

In the original study, Rhode & Zepf (2004) found that M86 hosts a population of ∼2900 GCs that extends ∼80 kpc from the galaxy center. The GC system shows at least two peaks in the B − R color distribution, as well as a detectable radial color gradient that arises because the red population of GCs is more centrally concentrated than the blue population. Later, Hargis & Rhode (2014) used the same Mosaic imaging data to examine the spatial distribution of the GC system of M86 and compare it to the galaxy's light distribution. They found that the diffuse galaxy light and the GCs had similarly flattened, elongated distributions. The total GC population, as well as the blue and red subpopulations, have nearly identical azimuthal distributions to the galaxy itself, with ellipticities  ∼ 0.4 and position angles of ∼ −60°.

Here, we are once again making use of the Mosaic images of the M86 field from Rhode & Zepf (2004), but this time our goal is to search for evidence of substructure across the entire field. Because M86 was the primary target of the original study, it appears at the center of the Mosaic images, and M84 is positioned on the western edge (see Figure 1). At the distance of M86 (∼17 Mpc), the 36' × 36' field of view of the Mosaic images corresponds to a physical area of roughly 180 kpc × 180 kpc. Besides M86 and M84, several other less luminous Virgo Cluster galaxies appear in the frames. The basic properties of the galaxies that appear in the Mosaic images are shown in Table 1 and drawn from the NASA/IPAC Extragalactic Database (NED), unless otherwise specified here or in the table note. The table lists the galaxy name (i.e., the NGC number, along with the Messier catalog number and/or the VCC number); the morphological type from RC3 (de Vaucouleurs et al. 1991); the total absolute magnitude in the V band (calculated by combining ${V}_{T}^{0}$ from RC3 with the distance modulus m − M in table column (5)); an estimate of the stellar mass of the galaxy, computed by combining the V-band absolute magnitude with the appropriate mass-to-light ratio for spiral, S0, or elliptical galaxies, drawn from Zepf & Ashman (1993); the galaxy distance modulus and distance in Mpc; and the heliocentric radial velocity of the galaxy. The table also includes (when available) information about the total number of GCs in the system (N_GC), the specific frequency S_N of the GC system (the number of GCs normalized by the V-band absolute magnitude of the galaxy, as defined by Harris & van den Bergh 1981), and T, the total number of clusters normalized by the stellar mass of the galaxy (as defined by Zepf & Ashman 1993). As the information in Table 1 shows, nearly all of the galaxies that appear in our images have distances within the range of ∼16–18 Mpc; one galaxy, NGC 4387, has a somewhat larger distance (21 Mpc), and another, IC 3303, has no redshift-independent distance measurement in NED. For the latter, we have assumed a distance modulus of 31.15 (17 Mpc) because the galaxy is included as a likely Virgo Cluster member in Binggeli et al. (1985).

Zoom In Zoom Out Reset image size

Table 1.  Properties of the Galaxies in the M86/M84 Images

Name	Type	${M}_{{\rm{V}}}^{{\rm{T}}}$	log10(Mass*)	m − M	Distance	Vr	NGC	SN	T
		(mag)	(log(M⊙))	(mag)	(Mpc)	(km s−1)
NGC 4374 (M84, VCC 0763)	E1	−22.3	11.8	31.32 ± 0.11	18.4	1017 ± 5	2800 ± 200	3.5 ± 0.5	4.1 ± 0.6
NGC 4387 (VCC 0828)	E5	−19.6	10.8	31.65 ± 0.73	21.4	565 ± 5	41 ± 9	0.6 ± 0.3	0.6 ± 0.3
NGC 4388 (VCC 0836)	Sb	−21.7	11.4	31.29 ± 0.20	18.1	2524 ± 1	70 ± 10	0.3 ± 0.1	0.6 ± 0.3
NGC 4402 (VCC 0873)	Sb	−20.0	10.7	31.02 ± 0.20	16.0	232 ± 5	110 ± 10	1.5 ± 0.3	2.2 ± 0.7
NGC 4406 (M86, VCC 0881)	E3	−22.3	11.9	31.17 ± 0.14	17.1	−224 ± 5	2900 ± 400	3.5 ± 0.5	4.1 ± 0.6
NGC 4425 (VCC 0984)	SB0	−19.4	10.6	31.13 ± 0.80	16.8	1908 ± 5	90 ± 10	1.6 ± 0.5	2.5 ± 0.8
NGC 4406B (VCC 0882)	dE3,N	−15.2	⋯	31.17 ± 0.14	17.1	1101 ± 55	⋯	⋯	⋯
IC 3303 (VCC 0781)	dS0,N	−17.5	⋯	⋯	⋯	−332 ± 6	⋯	⋯	⋯

Note. For NGC 4402, no V-band magnitude is given in RC3 or in the NED listing. Therefore, to calculate absolute V-band magnitude for this galaxy, we have combined ${B}_{T}^{0}$ from RC3 with an assumed B − V color of 0.65, which is typical for Sb spiral galaxies (Robert & Haynes 1994), and the distance modulus in Column (5). Distance modulus values for NGC 4374, NGC 4837, and NGC 4406 are from Tonry et al. (2001). Distance modulus values for NGC 4388 and NGC 4402 are from Tully et al. (2013), and the value for NGC 4425 is from Tully (1988). The morphological type for NGC 4406B is taken from Binggeli et al. (1985). There is no redshift-independent distance listed in NED for NGC 4406B, so we assume the same distance as that of NGC 4406, because these two galaxies are thought to be interacting (Elmegreen et al. 2000). IC 3303 also has no redshift-independent distance measurement in NED, but because it is identified as a Virgo Cluster member, we have assumed a distance of 17 Mpc (m − M = 31.15). For both NGC 4406B and IC 3303, we have taken ${B}_{T}^{0}$ from Binggeli et al. (1985) and combined it with the median B − V value of 0.77 for early-type dwarf galaxies from van Zee et al. (2004) and the distance modulus to calculate the M_V^T value listed in Column (3). The sources for the heliocentric radial velocities listed in Column (6) are as follows: the values for M84, NGC 4387, M86, and NGC 4425 are from Cappellari et al. (2011); the value for NGC 4388 is from Lu et al. (1993); the value for NGC 4402 is from Binggeli et al. (1985); the velocity for NGC 4406B is from Strauss et al. (1992); and the velocity for IC 3303 is from Albaretti et al. (2017). The GC system properties are drawn from Young (2016) in all cases except NGC 4406, which is taken from Rhode & Zepf (2004), and NGC 4402, which is drawn from K. L. Rhode et al. (2020, in preparation).

Download table as:  ASCIITypeset image

For the analysis of M86 that was published in Rhode & Zepf (2004), the other galaxies and GC populations in the Mosaic field were masked out so that the global properties of M86's GC system (i.e., the total number of GCs, along with the specific frequency, radial distribution, color distribution, and color gradient of the system) could be quantified. For the current analysis, our aim was to look for substructure and environmental effects in the GC populations in all the galaxies in the field, so we utilized the full Mosaic images without masking any large areas. To generate a list of GC candidates across the full field, we carried out the same basic steps that had been used for the Rhode & Zepf (2004) analysis. We began by creating smoothed versions of the Mosaic images using a ring filter with a diameter equal to 6 times the FWHM PSF of the images. We then subtracted the smoothed images from the original images, thereby removing the diffuse starlight associated with the galaxies. A constant background level was added to the galaxy-subtracted images in order to restore the sky background before the source detection and photometry steps. Next, we masked a small central portion of each galaxy (e.g., the central ∼200 × 200 pixels in M86 and M84) where the CCD image was saturated and/or where high noise levels in the galaxy-subtracted images prevented the detection of point-source GC candidates. We then used IRAF DAOFIND to detect the sources in each image, and we matched the source lists to produce a list of objects that appear in all three filters. We followed the same procedures that were described and illustrated in Rhode & Zepf (2004) to remove objects that appear extended in the Mosaic images: we measured the FWHM of the radial profile of each source and discarded objects that had larger-than-typical FWHM values for their instrumental magnitudes. At the distances of the Virgo Cluster galaxies in our field, GCs will be unresolved in our ground-based images, so objects with large measured FWHM are likely to be faint background galaxies rather than genuine GCs.

Again following the same procedures as were used in the original analysis in Rhode & Zepf (2004), we performed aperture photometry on the sources that survived the extended source cut. We measured the light from each source within a small aperture (with a radius equal to the FWHM PSF value a given image) and then applied the appropriate aperture correction and calibration coefficients to calculate calibrated V magnitudes and B − V and V − R for each source. The aperture corrections were calculated by measuring the magnitudes of a set of ∼20−40 bright stars in each image and computing the mean difference between the total magnitude and the magnitude within the small aperture. The aperture corrections (total magnitude minus magnitude within the small aperture) ranged from −0.15 to −0.21 mag with standard deviations of 0.01−0.02 mag. The Mosaic images of the M86+M84 field were acquired on a night with variable sky conditions, so we calibrated them using BVR exposures of the same field acquired under photometric conditions with the WIYN 3.5 m telescope² and Minimosaic camera. We measured magnitudes and colors for a set of bright stars that appeared in both sets of images and used these values to calculate a set of calibration coefficients (color transformation coefficients and zero-point constants) that could be applied to the Mosaic data.

We then used our final, photometrically calibrated point-source photometry to select GC candidates, i.e., unresolved objects with V magnitudes, B − V colors, and V − R colors that are consistent, within the photometric errors of each object, with the magnitudes and colors expected for GCs at the distances of the target galaxies. The BVR magnitude and color selection criteria were originally developed for the GC system survey mentioned earlier in this section (Rhode & Zepf 2001, 2004). For the current analysis of the galaxies that appear in the M86/M84 field, we selected all point sources in the field with V magnitude greater than 20.0 and colors within a specific region of the V − R versus B − V color–color relation defined by Milky Way GCs (Rhode & Zepf 2001). The V magnitude criterion is based on the assumption that the brightest GC in the target galaxies will have M_V ∼ −11 mag. We note that the galaxies in the M86/M84 field have a range of distance moduli (from 31.02 to 31.65), so M_V = −11 corresponds to a slightly different apparent V magnitude for each galaxy. However, because the B − V and V − R color also play a role in the selection, the final list of GC candidates remains the same whether we apply the V > 20.0 magnitude criterion across the entire field or apply a slightly different V threshold for objects around each galaxy that takes into account the varying distance moduli of the galaxies. The B − V and V − R color criteria are explained in detail in Rhode & Zepf (2001); briefly, we select objects that have B − V values between 0.56 and 0.99 (which is the B − V range expected for GCs with [Fe/H] between 0.0 and −2.5, based on the observed colors of Galactic GCs) and V − R values that put them within 3 times the standard deviation of the V − R versus B − V relation for Milky Way GCs. The selection process yielded a final list of 2250 GC candidates across the entire M86/M84 field. Figure 2 shows the GC candidates (blue filled symbols) in the V − R versus B − V color–color plane, along with the other unresolved sources in the field that were not selected (green open symbols). Note that because the photometric errors for each object are taken into account in the GC candidate selection, the selected candidates (the blue symbols in Figure 2) extend over a slightly larger region of the color–color plane than they would if the errors were not considered. The spatial locations of the GC candidates in the M86/M84 images are shown in the bottom panel of Figure 1.

Zoom In Zoom Out Reset image size

Removing extended objects from the source lists and then selecting GCs based on their magnitudes and colors in three filters significantly reduces the contamination in our sample (e.g., Rhode & Zepf 2001, 2004), although some contamination inevitably remains. Our analysis has shown that many of the contaminating objects that coincide with the GC selection box in the V − R versus B − V plane are compact background galaxies that are faint enough that they cannot be distinguished from point sources, and some are foreground Galactic stars (Rhode & Zepf 2001, 2004; Hargis & Rhode 2012).

To assess the amount of contamination in the GC candidate sample and then correct for it in subsequent steps, we created radial surface density profiles for the GC systems of M84 and M86. We assigned the GC candidates around each galaxy to a set of 1'-wide annuli, based on the candidates' projected radial distances from the galaxy center. We calculated the area of the portion of each annulus where GCs could be detected, i.e., the area minus those parts of the annulus that were masked out (because, e.g., they included bright foreground stars or areas with cosmetic defects on the CCD) or extended off the edges of the images. We then computed the surface density (number per unit area on the sky) of GC candidates for each annulus to construct the radial surface density profile of the system. The surface density values were highest near the center of each galaxy and then tapered off to a constant, flat value in the outer regions. We calculated the weighted mean of the surface density in the outer few annuli, where the surface density was constant, and used this as an estimate of the contamination level in the GC candidate sample. We carried out this analysis for M84 and M86 independently (masking a large region around the neighboring galaxy, and then carrying out the steps to construct the radial profile of the unmasked galaxy) and found that these estimated contamination levels matched each other within the uncertainties. We also checked to make sure that the contamination correction was the same on the east and west sides of M86 and on the north and south sides of M84. The estimated contamination level in the GC candidate list for the M84/M86 field is 0.2997+/−0.0268 per arcmin². Given the area of the images (1378 arcmin²), this works out to approximately 413 objects over the entire field, out of the 2250 objects in the GC candidate list.

Because we select GC candidates via their V magnitudes and B − V and V − R colors, any contaminants that are present in the GC candidate list will have similar magnitudes and colors to the true GCs in the sample. Nevertheless, we decided to apply a statistical correction to account for contamination in the GC candidate sample before carrying out our search for spatial features and substructure in the M86/M84 field. Details about the contamination corrections applied for each part of the analysis are given in Section 3.

We carried out a series of artificial star tests in order to quantify the detection limits for point sources in each of the Mosaic images and determine how much of the GC luminosity function (GCLF) we are able to detect in the target galaxies. We carried out these tests on the B, V, and R images separately. For each test, we added artificial point sources to a given image, executed the same set of steps we had used to detect the GC candidates, and then determined how many of the artificial sources were recovered. We added 800 artificial objects to each image, in steps of 0.2 mag, and repeated the process until we had covered a magnitude range of 5–6 mag in each filter. The results of the artificial star tests show that our detection process is 50% complete for point sources with B = 25.0, V = 24.0, and R = 23.4. We can combine the results of the completeness testing in each of the three filters in order to estimate how much of the intrinsic GCLF we have detected. Given our point-source detection limits in B, V, and R and the requirement that cluster candidates are detected in all three filters, and assuming an intrinsic GCLF with a peak magnitude (M_V ∼ −7.4 mag) and dispersion (σ ∼ 1.4 mag) (Brodie & Strader 2006, and references therein), we estimate that we have detected approximately 50% of the intrinsic GCLF (i.e., the brightest 50% of the clusters within the field of view of our images) of the two massive ellipticals M84 and M86.

The numbers above show that the detection completeness varies in each filter—the B image is substantially deeper than the R image. In addition to producing the full list of all 2250 GC candidates in our images, we also wanted to produce a list that we could use to explore how certain features in the spatial distribution of clusters vary when we examine subpopulations of clusters selected by their colors (see Section 3.2). In other words, we need to construct a sample of GCs that is equally complete in all three filters, to be sure that no color selection bias is present in the subsample being used to investigate trends with GC color. To accomplish this, we defined a sample of objects that we refer to as the "90% color sample," because it is at least 90% complete in all three filters. The reddest GC candidate in the full list of 2250 clusters has B − R = 1.87, and in the B-band our detection is 90% complete at 24.77, so a 90%-complete color sample would include all GC candidates with R brighter than 22.9. Applying this criterion to our GC candidate list yields a sample of 1718 objects. The locations of these 1718 objects in the V − R versus B − V color–color plane are shown in Figure 2, plotted in red open symbols on top of the 2250 GC candidates in the full sample.

In the GC system color distributions of elliptical galaxies, the blue, metal-poor peak often occurs around B − R ∼ 1.1, the gap between the two populations is typically around B − R = 1.23, and the second peak appears at B − R ∼1.4 (e.g., Rhode & Zepf 2001, 2004). Therefore, to investigate how the substructures in the M86/M84 field might vary for the different GC subpopulations, we divide our 90% color sample into two subsamples at B − R = 1.23. We use these blue and red subsamples in the analysis presented in Section 3.2. We also explain in that section how we used the estimated contamination level (discussed earlier) to apply a statistical correction for contamination to these subsamples.

We began our analysis by creating surface density maps of the GC candidates in the M84/M86 field. To look for evidence of substructure, we divided the field into a 150 × 150 grid and used KDE (e.g., Silverman 1986; Ivezić et al. 2014 for a detailed description of the technique) to calculate the surface density maps. Isodensity contours were then added to each map to help emphasize surface density features within the field. We used a Gaussian kernel with a fixed bandwidth of 086 (corresponding to ∼200 pixels). This bandwidth was chosen because it appeared to strike the best balance between smoothing away noise within the data and still allowing large-scale substructure to emerge.

3.1. M86 Surface Density Map

Figure 3 shows the surface density map for the full GC candidate sample. As explained in the previous section, some of the objects in the sample will not be true GCs but will instead be compact objects (foreground stars and background galaxies) that have magnitudes and colors like GCs at the distances of the target galaxies, and our sample of 2250 GC candidates should contain approximately 413 of such contaminating objects. To correct for contamination, we artificially cleaned the data by randomly selecting and removing 413 objects from the list of 2250 GC candidates. We then performed the KDE analysis on the cleaned data set. We carried out these steps a total of 100 times, removing a different random set of 413 objects each time. From the 100 realizations, we calculated the average surface density for each cell and used those averaged values to create the final KDE map shown in Figure 3. We find three major features, marked with an A, B, or C in the figure, along with three minor features. We discuss each of these features below.

Zoom In Zoom Out Reset image size

We begin by describing the minor features in the surface density map. First, there appears to be a surface density peak in the M84 GC system located ∼1' to the east of the center of M84. It is important to note that this is an artifact of our method rather than a true surface density anomaly. This offset is a product of edge effects caused by M84's position near the edge of our field. Due to the nature of the KDE technique, surface densities near the edge of the field will be systematically underestimated because of the fact that we cannot account for GCs that are present but fall outside our field of view. We varied the kernel size and found that the position of this surface density peak largely depended on this parameter, with smaller kernel sizes causing the peak to migrate to the west, toward the center of M84. However, decreasing the kernel size to the degree needed to remove the offset of the surface density peak greatly increased the prominence of noise within the image.

We also see a small connection between the GC systems of M86 and NGC 4402 that appears to indicate a tidal interaction. We used the results from our analysis of NGC 4402's GC system (see Table 1; K. L. Rhode et al. 2020, in preparation) to help us investigate whether this connection might be genuine evidence of a tidal interaction or simply the result of the two GC systems overlapping in projection on the sky. We removed the GC candidates around NGC 4402 within 2' (one R_eff for the GC system) of the galaxy center; there were 23 candidates inside 2' in the full sample and 13 blue candidates and 5 red candidates in the 90% color sample. We then repeated the KDE surface density map analysis but with those objects removed. Even with the removal of the GC candidates associated with NGC 4402, the connection still appears in all of the KDE surface density maps, although at a reduced strength.

Lastly, in the southeast corner of the field lies NGC 4425, an SB0 galaxy with a remarkably symmetric GC system. While the galaxy is relatively modest in size and mass (with a total magnitude of M^T_V = −19.4 and an estimated population of 90 ± 10 GCs; see Table 1), it is interesting that it has survived within the Virgo Cluster without, on first inspection, having had any significant tidal interactions to distort its GC system.

3.1.1. A Spatial Offset in the Peak Surface Density of the M86 GC System

Feature A in Figure 3 is a peak in the GC system surface density that is approximately 07 away from the center of M86. This surface density peak is similar to an overdensity reported by D'Abrusco et al. (2015) that lies in roughly the same area. In addition to the unexpected position of the surface density peak, there is a distortion in the central region of the M86 GC system where the surface density contours are elongated along the northeast-to-southwest direction.

Since the masked portion of the galaxy (the central ∼200 × 200 pixels, which are masked because they are saturated in the CCD images) is close to this surface density peak, we first attempted to determine whether this offset might be due to the same kinds of edge effects that created the offset in the M84 surface density peak. Decreasing the kernel size did not cause the peak to migrate toward the center of M86, so we created a 1'-wide annulus just outside and surrounding the masked region and calculated the surface density within this annulus. We then used this surface density along with the area of the masked region to determine the expected number of GC candidates that are likely to be inside the masked region and thus undetected in our images. We randomly distributed the same number of artificial GC candidates within the masked region and used KDE to create a new surface density map. The lopsided morphology of the innermost region of the M86 GC system seen in the original KDE map is still present in this new map.

To quantify the significance of this offset peak, we compared this finding to the expected case of an unperturbed, azimuthally symmetric GC system. To do this, we first needed to know the parameters that describe the M86 GC system. As mentioned in Section 2, Hargis & Rhode (2014) used the Mosaic images from which our sample was derived to study the spatial distribution of the GC system of M86. They found that the morphology of the inner region of the M86 GC system can be characterized by an ellipse with a semimajor axis length of 92, an ellipticity of 0.38, and a position angle of −63° east of north. We created an ellipse with these same characteristics and then divided it into six equal-area wedges with the center of the ellipse positioned to coincide with the center of M86. In an unperturbed GC system, we would expect each wedge to contain approximately the same number of GC candidates, but as can be seen in Figure 4, this is not the case. The left panel in Figure 4 shows our ellipse with the GC candidates and KDE contours plotted; the wedges are numbered, and we plot the number of GC candidates in each wedge in the right panel. The horizontal blue line is the expected number of GC candidates in each wedge, and the blue shaded region marks the ±1σ Poisson error on that number. The northeast wedge, number 6, has an excess of 20 GC candidates, ∼2.5σ higher than the expected number, suggesting that this offset GC candidate surface density peak is modestly significant.

Zoom In Zoom Out Reset image size

It is notable that the dwarf galaxy NGC 4406B (VCC 882) coincides with this overdensity. To estimate the number of GCs NGC 4406B might be contributing, we used the data gathered from a survey of dwarf elliptical galaxies by Miller & Lotz (2007). This survey included 13 nucleated dwarf ellipticals with absolute magnitudes between M_v = −15 and −16. The estimated number of GCs hosted by these 13 dwarf galaxies ranged from ∼2 to 22 clusters. Considering this, it is reasonable to assume that NGC 4406B contributes, at least in part, to the overdensity seen in feature A. We compared the color distribution of GC candidates within this overdensity to that of the GC candidates in the surrounding area. If NGC 4406B is contributing a significant number of clusters to this overdensity, we might expect the color distribution of the GC candidates within the overdensity to be shifted toward slightly bluer colors compared to the GC candidates in the surrounding area. For example, Peng et al. (2006) found that lower-luminosity galaxies have GC systems with bluer mean colors compared to higher-luminosity galaxies, although with plenty of scatter in the relation (see Figure 3 in their paper).

We found 40 GC candidates within 09 of the overdensity and compared their B − R colors to those of the GC candidates in the surrounding regions. We found no evidence to suggest that the GC candidates within this overdensity tend to be bluer than the candidates in the surrounding regions, but with such a small sample size and small expected differences in the color distribution this excess is likely to be difficult to detect. It is also possible that this overdensity is a result of interactions among the local GCs. Elmegreen et al. (2000) found that NGC 4406B is undergoing tidal stripping from M86, so a perturbation in the M86 GC system caused by NGC 4406B would potentially explain the existence of this overdensity.

3.1.2. A High Surface Density Bridge between M86 and M84

To the southwest of M86, in the area labeled with a B in Figure 3, an overdense region characterized by the presence of a few structural anomalies connects the GC systems of M86 and M84. This region incorporates the GC candidates associated with the moderate-luminosity elliptical galaxy NGC 4387, which has a contour around it that directly connects it to M84. Removal of the GC candidates associated with NGC 4387 does not noticeably affect the appearance of this bridge. The region labeled with a B contains several complex features—such as the overdense tail emerging to the east of M84 and the general lopsided morphology on the southwest side of M86—that are not consistent with the expected spatial distributions of the M86 and M84 GC systems, but instead may be due to an interaction between the two populations.

This region also shares qualitative similarities to structures formed in simulations of galaxy interactions in clusters of galaxies (Rudick et al. 2011; Ramos-Almendares et al. 2018), supporting the possibility that it is the result of an interaction between M86 and M84.

3.1.3. Flattened Contours on the Southeast Side of M86

The area of the GC system along the southeast edge of M86, labeled with a C in Figure 3, has flattened isodensity contours. This flattening breaks from the expected elliptical shape of the M86 GC system morphology (see Section 2). We do not see a similar flattening in the contours along the northwest edge of the M86 GC system. This gives an "egg-like" shape to the isodensity contours on the southeast side of M86 and suggests that an accretion event or tidal interaction may have occurred that redistributed the GCs in that region.

We compared the locations of the features we found in the GC populations to results from two studies that searched for LSB substructures in this area of the Virgo Cluster. Janowiecki et al. (2010) used V-band surface photometry to probe the halos of the five brightest ellipticals in the Virgo Cluster for diffuse light features. To observe features nearer the galaxy center, they fitted elliptical isophotes to the galaxy light and subtracted them, leaving behind the diffuse light features. Mihos et al. (2017) searched the intracluster light between and around galaxies within the Virgo Cluster to detect LSB features that could indicate past interactions or hierarchical processes. To achieve the depth required to observe features formed from the intracluster light, Mihos et al. (2017) stacked their images and then removed contamination from Galactic cirrus. Figure 5 is our surface density map, with the yellow regions marking features detected in Janowiecki et al. (2010). From Figure 5 we can see that the diffuse light structures reported by Janowiecki et al. (2010) have little to no obvious correlation with the features in our GC candidate surface density map. The region directly to the north of M84 appears to coincide with a small overdensity in the surface density map, but because of the proximity to the edge of our field, it is possible that this overdensity in the GC system is a spurious result due to edge effects. On the other hand, we do see a similarity between our detected features and a feature that was found in the Mihos et al. (2017) study of the M86/M84 field. Mihos et al. (2017) detected an extended LSB feature along the southeast (SE) edge of the M86 halo that did not match the elliptical isophotes and may, they suggested, be a consequence of a major merger event. This substructure is in the same direction as feature C, and both are characterized by boxier contours, although the feature we find is comparatively closer in to the galaxy center, whereas the "SE Shelf" identified by Mihos et al. (2017) is formed from diffuse light in the outermost regions of the galaxy.

Zoom In Zoom Out Reset image size

3.2. Surface Density Maps of the Red and Blue GC Candidates

Figure 6 shows the surface density maps of the subsamples of GC candidates that are classified as either blue (metal-poor) or red (metal-rich) according to the criteria described in Section 2. To construct this sample and correct it for contamination, we began with the full list of 2250 GC candidates. We randomly removed 413 objects from the sample and then applied the appropriate R magnitude cut (selecting GC candidates brighter than R = 22.9; see Section 2) to produce a contamination-corrected 90% color sample. We then divided the contamination-corrected sample into blue and red subsamples at B − R = 1.23. As we had done when we constructed the main KDE map shown in Figure 3, we repeated these steps (randomly removing GC candidates to account for contamination, and then applying completeness cuts and splitting the samples according to color) 100 times, removing a different selection of GC candidates each time. In each realization, the surface density was estimated at each cell in our grid, and the average surface density for each grid cell was calculated from this collection of estimates.

Zoom In Zoom Out Reset image size

The blue GC candidates, shown in the top panel of Figure 6, have a more extended spatial distribution compared to their red counterparts. In the surface density map of the blue GC candidates, we see a recurrence of the surface density peak to the northeast of the M86 GC system center that was visible in the map of the full GC candidate sample (feature A from Figure 3). The contours in the inner region of the blue component of the M86 GC system are increasingly elongated along the northeast–southwest axis, suggesting that this elongation comes primarily from the blue GC candidates in the sample. Likewise, the bridge connecting the M86 and M84 GC systems (feature B from Figure 3) is prominent, with a high surface density arm connecting M84 to NGC 4397 and another high surface density arm emerging from the M86 GC system along the southwest direction. The southeast shelf is still present but has broken down into patchier, high surface density regions. Furthermore, many of the less luminous galaxies in the field have high-density regions associated with them in the blue GC candidate KDE map, but these high-density regions are much weaker, if they exist at all, in the red GC candidate KDE map. This is not unexpected, because lower-luminosity galaxies tend to have smaller proportions of metal-rich GCs than more luminous galaxies do (e.g., Peng et al. 2006).

In contrast to the blue GC candidates, the red GC candidates (bottom panel of Figure 6) are more centrally concentrated. The bridge linking the M86 and M84 GC systems is much weaker, but still present, with a connection formed just north of NGC 4387. The asymmetric surface density peak to the northeast of the M86 GC system center has almost entirely disappeared. The surface density contours in the central region of the red component of the M86 GC system also do not show the same distortion along the northeast–southwest direction that is seen in the blue component, reinforcing that these two features are dominated by metal-poor GCs. The southeast shelf is largely nonexistent in the surface density map of the red GC candidates. Instead, the central region of the red component of the M86 GC system has boxier isodensity contours than those in the blue GC candidate surface density maps.

3.3. Spatial Positions of GCs with Velocities

Having kinematic data for some of the GC candidates in our sample should provide us with additional insight into the accretion history and ongoing interactions between the GC populations in this well-populated Virgo Cluster field. To further investigate the possibility of substructure in the GC populations of the galaxies, we explored kinematic data published by Park et al. (2012) and Ko et al. (2017). Park et al. (2012) measured the radial velocities of 25 GCs around M86 and reported a mean GC velocity of ${v}_{p}=-{354}_{-79}^{+81}$ km s⁻¹ with a velocity dispersion of σ_p = 292 ± 32 km s⁻¹. Ko et al. (2017) carried out a wide-field spectroscopic survey of GCs within the Virgo Cluster. Of the 201 GC spectra in the Ko et al. (2017) sample, 94 GCs were within our field, with one of these GCs also appearing in the Park et al. (2012) sample. After combining the Park et al. (2012) and Ko et al. (2017) data sets and removing the one duplicate GC, we had a total of 118 GCs with radial velocities in our field to study.

Figure 7 shows a map of the spatial distributions of the collated sample of GCs within our field from the Park et al. (2012) and Ko et al. (2017) surveys plotted on top of the surface density map of our full, contamination-corrected GC candidate sample. The points representing the GCs are color-coded according to their radial velocity in km s⁻¹, and the spatial positions of galaxies within the field are denoted by stars, also color-coded according to their radial velocities as reported in Table 1. From this figure we see that GCs with measured velocities populate all of the major features discussed previously. The GC candidate surface density peak to the northeast of M86's spatial position has three clusters with negative radial velocities, consistent with the velocity of M86 (−224 ± 5 km s⁻¹). The bridge of elevated GC candidate surface density that connects the M86 and M84 GC systems is populated by several GCs with varying radial velocities, suggesting that it is composed of GCs from both galaxies and, possibly, intracluster GCs. The flattened shelf along the southeast side of the M86 GC system is mostly populated by GCs with velocities consistent with the radial velocity of M86, but there are two GCs with high recession velocities more akin to the radial velocities of M84 (1017 ± 5 km s⁻¹) or NGC 4406B (1101 ± 55 km s⁻¹).

Zoom In Zoom Out Reset image size

Most of the less massive galaxies in the field appear to host no, or very few, of the GCs with measured velocities. The barred lenticular galaxy NGC 4425 and the spiral galaxy NGC 4388 each have two clusters in close proximity with recession velocities similar to those of the galaxies, suggesting that the clusters are associated. NGC 4402 also has two clusters with similar velocities nearby, but its close proximity to M86 makes it difficult to assess to which galaxy's GC system these clusters belong.

We plot the radial velocity of these GCs as a function of their projected radial distance from the center of M86 in Figure 8. The clusters have been color-coded based on their proximity to either M86 or M84, with green points representing clusters with a smaller projected radial distance to the center of M86 and orange points representing clusters with a smaller distance to the center of M84. The radial velocities of M86 and M84 (as reported by Cappellari et al. 2011) are represented as the green and orange dashed lines, respectively, with the vertical, black dashed line denoting the projected radial distance between M86 and M84. There is clearly a large gap between ∼300 and ∼800 km s⁻¹ where GCs are lacking; the large kinematic separation could be indicative of a fly-by interaction between M86 and M84. This plot also confirms the existence of several clusters with kinematics that are not consistent with the expected velocity based on their closest massive elliptical neighbor.

Zoom In Zoom Out Reset image size

We attempted to quantify the significance of differences between the radial velocities of 114 GCs in this sample and the radial velocities of each GC's closest massive elliptical galaxy. We removed both sets of clusters nearest NGC 4425 and NGC 4388, as these GCs have radial velocities in close agreement with those of NGC 4425 and NGC 4388 and are unlikely to be hosted by either M86 or M84. To determine whether a GC had a radial velocity that was significantly different from the mean velocity of the nearest massive elliptical galaxy, we needed to estimate the velocity dispersion of the M86 and M84 GC systems. We did this using the biweight method described in Beers et al. (1990) with GCs within one effective radius of M86 and GCs within one effective radius of M84. We limited ourselves to GCs within one effective radius of each galaxy because of the significant overlap between the two systems. With this method we report a velocity dispersion of ${\sigma }_{v,{\rm{M}}86}={308}_{-49}^{+60}$ km s⁻¹ and ${\sigma }_{v,{\rm{M}}84}={265}_{-62}^{+83}$ km s⁻¹ calculated from 46 and 25 GCs, respectively.

In Figure 9, we color-code each GC based on the difference between the cluster's velocity and the average velocity of the closest massive elliptical galaxy (normalized to the velocity dispersion of the closest massive elliptical galaxy). The bridge connecting the GC systems of M86 and M84 hosts 10 GCs. The GCs within this connection form a rough line from the northwest to the southeast, across the feature itself. These GCs have a range of velocities, from −454 to +1081 km s⁻¹, and all of them have velocities consistent with those expected for GCs associated with either M86 or M84. However, as shown in Figure 9, six of these GCs have radial velocities that are 3σ−6σ different from the radial velocity of their closest massive elliptical neighbor. Several of the clusters in the field surrounding the two galaxies, along with two clusters located between NGC 4425 and NGC 4388, have velocities consistent with the GCs of M84. Similarly, there are a few clusters that are located closer to the center of M84 but have velocities consistent with that of M86. Although some of the clusters in question may have been stripped by tidal interactions between the two massive ellipticals, they may also simply represent clusters in the outer regions of the GC systems of M84 and M86. A significantly larger sample of velocities would be needed to construct a complete characterization of the phase space of these GC systems, to explore the possible fingerprints of the tidal interaction between M84 and M86, and to identify a population of intracluster GCs (see, e.g., Romanowsky et al. 2012; Longobardi et al. 2015, 2018).

Zoom In Zoom Out Reset image size

In this paper we presented results from our search for substructure within the GC populations of M86 and M84 using the spatial positions of 2250 GC candidates selected from wide-field, multicolor broadband imaging. We estimated the amount of contamination in the GC candidate sample and then applied statistical corrections to account for it. We then used the KDE technique to create surface density maps of the GC systems within the M86/M84 field and of the red and blue GC candidate subpopulations separately. We combined our KDE analysis results with kinematic data drawn from Park et al. (2012) and Ko et al. (2017) and analyzed the spatial positions and radial velocities of a subset of the GCs to see whether the GC kinematics correlate with the features we detect. Our results can be summarized as follows:

• 1.  
  We identify three main substructures within the field. First, we find that the surface density peak of the M86 GC system is offset from the system's center by approximately 07. This overdensity coincides with the dwarf elliptical galaxy, NGC 4406B. Although the GC system of NGC 4406B has not been studied, dwarf elliptical galaxies with similar luminosities can host enough GCs to account, at least partially, for this overdensity. Moreover, Elmegreen et al. (2000) have found evidence for an interaction between M86 and NGC 4406B, so a perturbation within the M86 GC system caused by NGC 4406B is within the realm of possibility. The second structure we detect is a bridge connecting the M86 and M84 GC systems, perhaps as a result of interactions between the two galaxies. The final structure is a flattened isodensity contour along the southeast edge of the M86 GC system. Mihos et al. (2017) found a similar structure in the LSB starlight around M86 and called it the "SE Shelf," arguing that it may be the result of a major accretion event in the galaxy's past.
• 2.  
  When we divide our GC candidates into red (metal-rich) and blue (metal-poor) subpopulations, we find that the blue GC candidates are the dominant population in the structures mentioned above: the various features that are apparent when we examine the total GC candidate sample (the surface density peak offset from the M86 GC system center, the bridge connecting the M86 and M84 GC systems, the flattened isodensity contour along the southeast of the M86 GC system) are also apparent when we observe the blue GC candidate subsample. In contrast, these features have a marked decrease in strength when analyzing the surface density map of the red GC candidate sample.
• 3.  
  We find several examples of GCs with radial velocities that are inconsistent with the velocity of the nearest giant elliptical galaxy. The majority of such GCs are located within the bridge connecting the M86 and M84 GC systems. It is possible that the GCs within this bridge have been stripped from their host galaxy, but they may instead simply be part of the outer GC population of these two massive elliptical galaxies. Our analysis also shows possible evidence for a stream of intracluster GCs with velocities consistent with that of M84 emerging along the east–west axis of M86 and to the south of the M86 GC system.

The research described in this paper was supported in part by a Daniel Kirkwood Fellowship from the Indiana University Astronomy Department to R.A.L. We are grateful to the anonymous referee for their careful review of the manuscript and valuable suggestions that improved the paper. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research made use of Astropy,³ a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018).