Abstract
We present a study of the stellar populations of globular clusters (GCs) in the Virgo Cluster core with a homogeneous spectroscopic catalog of 692 GCs within a major-axis distance Rmaj = 840 kpc from M87. We investigate radial and azimuthal variations in the mean age, total metallicity, [Fe/H], and α-element abundance of blue (metal-poor) and red (metal-rich) GCs using their co-added spectra. We find that the blue GCs have a steep radial gradient in [Z/H] within Rmaj = 165 kpc, with roughly equal contributions from [Fe/H] and [α/Fe], and flat gradients beyond. By contrast, the red GCs show a much shallower gradient in [Z/H], which is entirely driven by [Fe/H]. We use GC-tagged Illustris simulations to demonstrate an accretion scenario where more massive satellites (with more metal- and α-rich GCs) sink further into the central galaxy than less massive ones, and where the gradient flattening occurs because of the low GC occupation fraction of low-mass dwarfs disrupted at larger distances. The dense environment around M87 may also cause the steep [α/Fe] gradient of the blue GCs, mirroring what is seen in the dwarf galaxy population. The progenitors of red GCs have a narrower mass range than those of blue GCs, which makes their gradients shallower. We also explore spatial inhomogeneity in GC abundances, finding that the red GCs to the northwest of M87 are slightly more metal-rich. Future observations of GC stellar population gradients will be useful diagnostics of halo merger histories.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
The assembly of the most massive galaxies in the universe is a story of continual merging and accretion of gas, stars, and dark matter over a Hubble time. These behemoths have long ago stopped forming significant numbers of their own stars in situ, but they are still adding stars to their halos through the accretion of lower-mass galaxies (Oser et al. 2010; Johansson et al. 2012). When massive galaxies are situated at the centers of galaxy clusters, their stellar halos ("cD envelopes") join seamlessly into the larger accreted population sometimes referred to as "intracluster light" (Murante et al. 2004, 2007; Purcell et al. 2007; Contini et al. 2014). While the relatively low surface brightnesses of these halo stars make them difficult to study, they contain a wealth of information on the progenitor components that assembled to form the parent galaxy. In the most massive galaxies, the accreted halo stars may actually be the dominant stellar population (Purcell et al. 2007; Oser et al. 2010; Cooper et al. 2013).
The characterization of extragalactic halo stellar populations is made challenging by their low surface brightnesses, usually well below that of the sky. Deep imaging surveys of massive galaxies have revealed extended stellar halos in cluster galaxies (Bernstein et al. 1995; Gonzalez et al. 2000; Feldmeier et al. 2004; Mihos et al. 2005, 2017; Krick & Bernstein 2007; Rudick et al. 2010), some showing significant fine structure expected from merging (Duc et al. 2015; Bílek et al. 2020), and some have been resolved into individual stars (Mouhcine et al. 2005a, 2005b, 2007; Williams et al. 2012). Obtaining further information on the nature of these stellar halos is more difficult, as even colors are subject to uncertainties with sky subtraction (e.g., Mihos et al. 2013). Recent spectroscopic surveys have been able to derive the integrated ages, metallicities, and α-element abundances of the brighter regions (R < 3 effective radii) of the stellar halos of some massive early-type galaxies (e.g., Spolaor et al. 2010; Greene et al. 2015; Gu et al. 2018, 2018; Iodice et al. 2019; Fensch et al. 2020).
Numerous studies have used discrete tracers like planetary nebulae (PNe; e.g., Hui et al. 1993; Arnaboldi et al. 1996, 1998; Ciardullo et al. 2004; Buzzoni et al. 2006; Ciardullo 2010; Cortesi et al. 2013; Longobardi et al. 2013, 2015a) and globular clusters (GCs; e.g., Beasley et al. 2000; Peng et al. 2006; Hwang et al. 2008; Liu et al. 2011; Pastorello et al. 2015; Usher et al. 2015; Lim et al. 2017; Longobardi et al. 2018b; Ko et al. 2018, 2019, 2020; Powalka et al. 2018) with success to elucidate the stellar content of halo populations. In particular, colors of old GCs are indicative of metallicities, which have been broadly used to trace the mass assembly history of host galaxies. The GCs in massive galaxies commonly show bimodal color distributions, with the blue and red GCs corresponding to metal-poor and metal-rich subpopulations, respectively. The properties of red GCs show tighter correlations with the underlying stars in their host galaxies in terms of their spatial distributions (e.g., Park & Lee 2013; Wang et al. 2013) and kinematics (e.g., Peng et al. 2004; Pota et al. 2013). Several studies have suggested a scenario where the red (metal-rich) GCs were formed with the bulk of the stars during dissipative mergers, while the blue (metal-poor) GCs were formed in small galaxies and have been continuously accreted onto massive galaxies (e.g., Côté et al. 1998; Lee et al. 2010b).
We can better understand these possibilities by combining chemical and spatial information. Metallicity gradients, a common feature of galaxies, inform us about the role of the gravitational potential in the regulation of star formation and chemical enrichment, but they also highlight the role of major mergers in erasing these gradients (e.g., Spolaor et al. 2010). Gas-poor minor merging can also create metallicity gradients, as higher-mass galaxies with more metal-rich stars have a tendency to sink further into the host galaxy owing to dynamical friction (Amorisco 2019). For example, in the Milky Way's stellar halo, the existence or not of abundance gradients has been important in shaping our view of how the halo may have formed. The lack of a metallicity gradient for Milky Way GCs beyond 8 kpc led Searle & Zinn (1978) to propose that the outer halo was formed from the accretion of low-mass fragments.
The color gradients of GC systems have been explored in a relatively large sample of early-type galaxies in the Virgo and Fornax Clusters (Liu et al. 2011), finding that both the blue (metal-poor) and red (metal-rich) GC systems exhibit radial color gradients, showing equal strengths within the uncertainties. This result holds even when gradient measurements are extended to R ∼ 5R e –8R e , as in the literature compilation of Forbes & Remus (2018), but flattening of the gradients of metal-poor GCs was found at R ≳ 8R e .
Another chemical parameter for understanding the formation history of galaxies is the ratio between α-element abundances and iron, [α/Fe]. This provides information about the star formation timescale, i.e., a high [α/Fe] value indicates rapid star formation and early quenching. As an example, the correlation of [α/Fe] in Milky Way stars with their locations in the Galaxy has been used to understand the formation of thin and thick disks (e.g., Recio-Blanco et al. 2014). For GCs in massive early-type galaxies, it is known that they have a wide range of [α/Fe] values and often have supersolar [α/Fe] ratios (e.g., Puzia et al. 2005; Park et al. 2012), which indicates that they were formed over a short timescale. Studies of the spatial variation of [α/Fe] in extragalactic GCs, however, are rare, but they have the potential to reveal in detail how galaxy halos are assembled.
Recently, there have been attempts to understand the characteristics of the GC populations in the context of cosmological simulations. The E-MOSAICS project (Pfeffer et al. 2018; Kruijssen et al. 2019; Reina-Campos et al. 2019) combined the MOSAICS semianalytic models for star cluster formation and evolution with the EAGLE hydrodynamical simulation for galaxy formation, including 25 Milky Way–like disk galaxies. They studied the formation and coevolution of galaxies and their star clusters in terms of chemical abundances, progenitors, and accretion epochs of the GCs. In another approach, Ramos-Almendares et al. (2020) investigated the GC populations in nine simulated Virgo-like galaxy clusters. They tagged the GCs to the simulated galaxies at their infall time and traced the kinematics and spatial distributions of the GCs. Other groups have also started to tag GCs in cosmological simulations to study the spatially resolved properties of GC systems (Trujillo-Gomez et al. 2021; Chen & Gnedin 2022). These developments have made it possible to use the spatial gradients and distributions of GCs in meaningful comparisons to simulations.
The Virgo Cluster is an excellent target to investigate the mass assembly history of the central galaxy because of its proximity. It is a dynamically young system that shows many substructures of member galaxies both spatially and kinematically (Binggeli et al. 1985, 1987). In particular, the core region of the Virgo Cluster shows diffuse stellar features (Mihos et al. 2005, 2017) and extended structures of ionized gas filaments (Kenney et al. 2008; Boselli et al. 2018) and H i gas (Yoshida et al. 2002; Oosterloo & van Gorkom 2005) around galaxies, indicative of strong galaxy interactions. The most massive galaxy at the center of the Virgo Cluster, M87, has a weak cD envelope, and kinematic studies of GCs and PNe have been done to understand the connection between its faint stellar halo and the intracluster light (Longobardi et al. 2015a; Ko et al. 2017; Longobardi et al. 2018a, 2018b). In addition, the GCs and PNe in M87 also constitute kinematic substructures related to recent mergers (Strader et al. 2011; Romanowsky et al. 2012; Longobardi et al. 2015b; Oldham & Evans 2016; Lambert et al. 2020). Therefore, M87 and the Virgo Cluster are a good laboratory to study the ongoing assembly of the most massive galaxy in a cluster core.
In addition to the many aforementioned GC kinematic studies, the chemistry and ages of GCs have also been studied from photometry and spectroscopy in Virgo ellipticals (e.g., Cohen et al. 1998, 2003). Powalka et al. (2018) reported a substructure of relatively young GCs located ∼115 south of M87, possibly related to past mergers or accretions. More recently, Villaume et al. (2020) extended the study of chemical abundances in GCs out to ∼140 kpc from M87, showing the promise of stellar population studies in extragalactic GCs across a range of radius. They do not find any significant iron abundance gradients of both metal-poor and metal-rich populations within R < 40 kpc and conclude that the flat metallicity gradient of metal-poor GCs could be a sign that some metal-poor GCs were formed in situ at high redshift. Moreover, they found that the metal-poor GC population has a higher α-element abundance than the metal-rich GC population in both the inner (R < 40 kpc) and outer (R > 40 kpc) halos, which is not expected by current cosmological simulations.
In this study, we investigate mean trends of the stellar populations of the GCs in the core region of the Virgo Cluster out to ∼840 kpc (29), as a function of major-axis distance from M87, using a new homogeneous spectroscopic catalog. This survey mainly covers a 3° × 2° field including the massive elliptical galaxies M87, M86, and M84. This paper is organized as follows. The GC catalog used in this study and the sample selection are described in Section 2. In Section 3, we explain how we measure ages and metallicities of the GC subpopulations. We present the main results about the age and metallicity distributions of the GCs in Section 4. We discuss the assembly history of M87 and the Virgo Cluster core in Section 5. In Section 6, we summarize the results. We adopt a distance to the Virgo Cluster of 16.5 Mpc (Mei et al. 2007; Blakeslee et al. 2009), so 1 corresponds to 4.91 kpc.
2. Data and GC Sample Selection
2.1. Sample Definition
We produced a new homogeneous catalog of the GCs in the Virgo core region by combining two spectroscopic data sets. The first data set was part of a spectroscopic follow-up program to the Next Generation Virgo Cluster Survey (NGVS; Ferrarese et al. 2012). The NGVS is an extensive deep imaging survey using the MegaCam imager on the Canada–France–Hawaii Telescope (CFHT). It images a 104 deg2 field of the Virgo Cluster with five optical filters . The images were taken with excellent seeing conditions (<06 for i-band imaging), and many GCs in the Virgo Cluster are marginally resolved in the NGVS images. The -band "fuzziness" index (Δi4−8), defined as the difference between 4-pixel and 8-pixel -band aperture-corrected magnitude, has larger values for the GCs than for foreground stars on average (Durrell et al. 2014). In addition, a Ks -band imaging survey was carried out with WIRCam on the CFHT (NGVS-IR; Muñoz et al. 2014) for a 2° × 2° region around M87. Muñoz et al. (2014) showed that the color–color diagram allows us to clearly distinguish Virgo GCs from contaminants such as foreground stars and background galaxies. We identified GCs in the Virgo core based on a "fuzziness" index and various color combinations of point sources, adopting the "extreme deconvolution" (xD) classification algorithm (Bovy et al. 2011). GCs were selected using xD, but also using a broader set of color and morphology cuts, and targeted for spectroscopy. This program obtained spectra of 776 GCs in the Virgo Cluster core, some of which have been used in previous studies (Zhu et al. 2014; Li et al. 2020).
Second, we used spectra taken from Ko et al. (2017), who presented a spectroscopic catalog of the GCs in a central 2° × 2° field of the Virgo Cluster, covering the vicinity of M87, M84, and M86. They obtained 910 optical spectra of GC candidates. We reclassified their sample using the color–color diagram for the objects in the NGVS-IR field and confirmed 212 GCs in their sample. There are 88 sources in common between the NGVS follow-up program and Ko et al. (2017). Finally, we combined these two samples to construct a combined spectroscopic catalog of 903 GCs in total in the Virgo core.
Both programs were taken with the MMT/Hectospec (Fabricant et al. 2005) with a 270 mm−1 grating having a dispersion of 1.2 Å pixel−1. Data were obtained in the wavelength range 3650—9200 Å. Data reduction was done with the SAO pipeline, following the description in Caldwell et al. (2011), and flux calibration was done following the steps in Fabricant et al. (2008). We used the xcsao task in the IRAF RVSAO package (Kurtz & Mink 1998) to measure heliocentric radial velocities of targets, taking 10 different templates including stars, GCs, and galaxies. Combining radial velocity, color, and size information, we classify the objects into four groups: GCs, UCDs, foreground stars, and background galaxies. The details of the data reduction, flux calibration, heliocentric radial velocity measurements, and target classification will be described in a forthcoming paper about the spectroscopic catalog of the GCs (Y. Ko et al., 2022, in preparation). The median signal-to-noise ratio (S/N) of the spectra of GC candidates at 5000 Å is S/N ∼ 9, where the Poisson noise, sky noise, and readout noise are considered to calculate spectral noise as a function of wavelength.
We select 726 out of the 903 total GCs, excluding the GCs that belong to other Virgo galaxies in the core. The remaining GCs are either GCs that belong to M87 or intracluster GCs (IGCs), a population governed by the galaxy cluster gravitational potential (Williams et al. 2007; Firth et al. 2008; Lee et al. 2010a; Durrell et al. 2014). We did this because our interest for this paper is mainly to investigate the assembly history of M87 in the cluster environment. There are 154 Virgo galaxies that have radial velocity information in the survey region. We used both spatial and kinematic criteria to separate the GCs bound to Virgo galaxies, following Longobardi et al. (2018b): (1) R < 10R e , where R e is the effective radius of a given galaxy; and (2) ∣v r − vgal∣ < 3σgal where v r , vgal, and σgal are the radial velocities of GCs, galaxy systemic velocity, and galaxy central velocity dispersion, respectively. The radial velocities of galaxies are compiled from various literature sources (NASA/IPAC Extragalactic Database (NED) 2019; Ahn et al. 2014; Kim et al. 2014; Ferrarese et al. 2020). The velocity dispersions are from the GOLDMine database (Gavazzi et al. 2003) compiling data from McElroy (1995) and Scodeggio et al. (1998). For the galaxies where velocity dispersions are not known, mostly faint dwarf galaxies, we adopted σgal = 50 km s−1.
We expect that these cuts are very conservative and remove from our sample nearly all the GCs that belong to other galaxies. As an example, for larger galaxies, the effective radius of the GC system (R e ,GC) may be as big as ∼5 times that of the stars (e.g., Kartha et al. 2016; Caso et al. 2019), so our cut to eliminate GCs within 10R e ,stars corresponds to a cut of 2R e ,GC, which results in ∼ 85% of GCs being removed, leaving ∼ 15% to be potential contaminants. The velocity cut, however, removes all GCs within ± 3 × the velocity dispersion from the mean galaxy velocity. If the GC line-of-sight velocity distribution (LOSVD) is roughly Gaussian, then this selection removes all but 0.27% of the GCs. Combined with the spatial cut, only 0.04% of the GC system could potentially contaminate our sample. Our spectroscopic limit ( mag) and a Gaussian GC luminosity function shape (〈g〉 ≈ 24 mag, σ ∼ 1.3 mag; Jordán et al. 2007) mean that we only observe the brightest ∼ 12% of the GCs in any massive galaxy. Therefore, the fraction of GCs that could enter our sample is roughly 5 × 10−5. The massive nearby galaxies M86 and M84 have ≈ 3000 GCs each (Peng et al. 2008). The analysis above therefore leads us to expect that only 0.3 GCs from these galaxies could contaminate our sample. Of course, this analysis assumes that the GC spatial distribution follows a Sérsic profile at large radii and that the LOSVD is Gaussian in the wings. These simple assumptions may not hold in practice, but any GC that is significantly detached (spatially and kinematically) from its host is probably on its way to being accreted by the M87 system and arguably should be included in our analysis.
Figure 1(a) shows the radial velocities of 726 GCs as a function of major-axis distance from M87. The major-axis distance Rmaj is calculated as follows:
with a position angle (PA) of 155° and an ellipticity e of 0.4 for M87 (Ferrarese et al. 2006; Janowiecki et al. 2010). Most GCs have radial velocities consistent with a mean value of 1255 km s−1, similar to the systemic velocity of M87 (1284 km s−1; Cappellari et al. 2011). In addition, we identify another group of GCs at , that appears to be a broad feature in the velocity−position phase-space diagram connecting M87 with M86. The mean radial velocity of these GCs is v r = − 66 km s−1 , which is about 1370 km s−1 lower than the systemic velocity of M87.
Figure 1(b) shows the spatial distribution of the GCs in our study. The GCs are strongly concentrated around M87, although some of them are uniformly distributed in the survey region beyond . The GCs with low radial velocities are mostly located in the region west of M87 (between it and M86). Note that the low-velocity GCs that are strongly concentrated on M86 (systemic velocity of −224 km s−1) are considered to belong to it and are not plotted in this figure. For this analysis, we only used the 692 GCs that have higher radial velocities in order to focus on the main structure of the Virgo Cluster, assuming that the low-velocity GCs constitute an infalling group associated with M86.
2.2. GC Subsamples
Many studies have confirmed the existence of color bimodality for GCs in M87 (Whitmore et al. 1995; Elson & Santiago 1996; Kundu et al. 1999; Larsen et al. 2001; Peng et al. 2006, 2009; Tamura et al. 2006; Harris 2009), a feature typical of GCs in massive early-type galaxies. Since we expect that the blue and red GC subpopulations have different formation origins, we divided our GC sample into blue and red, using their colors, to search for any differences in their stellar populations.
We performed Gaussian mixture modeling (Muratov & Gnedin 2010) on the color distribution of the GCs in a 2° × 2° region around M87, identified in the NGVS photometric survey (see Section 2.1) within radial bins 5 in width, assuming heteroscedastic bimodal distributions. Figure 2(a) shows the color of the photometric sample as a function of major-axis distance from M87. We adopted the color at which GCs have a 50% probability of belonging to either component as the color that divides the two components. This color criterion changes as a function of major-axis distance. We adopted a broken linear fit for the dividing color for GCs with as follows:
Download figure:
Standard image High-resolution imageThe color criterion becomes almost flat at beyond . It is worth mentioning that the main radial gradient results are not dependent on the details of the color criterion and are little changed even if we adopt a single color at all radii, such as mag.
Figure 2(b) shows the color of our GC samples as a function of major-axis distance from M87. The majority of the GCs are located within where a strong sequence at and a weak red tail are shown. Most of the red GCs are located within , while the blue GCs extend out to . We have 509 blue and 183 red GCs in our sample with the color criteria described above.
3. Age, [Z/H], [Fe/H], and [α/Fe] Measurements
The colors of GCs are an efficient tool to measure their metallicities. However, there are inherent uncertainties because GC colors also depend on GC ages. A spectral analysis has the advantage of measuring ages and metallicities of GCs more independently.
Lick line indices have been commonly used to measure ages and metallicities of old simple stellar populations with their integrated spectra. They are defined as the strengths of 25 prominent absorption lines seen in optical spectra of old stellar systems, established and refined by Burstein et al. (1984), Worthey et al. (1994), Worthey & Ottaviani (1997), and Trager et al. (1998).
We used the lick_ew code in the EZ_Ages package (Graves & Schiavon 2008) to measure the Lick indices from our spectra, adopting the passbands defined by Worthey & Ottaviani (1997) and Trager et al. (1998). This code smooths spectra to the resolution of the Lick system and estimates index errors based on error spectra (Cardiel et al. 1998).
We measured the mean ages, metallicities, and α-element abundances of the GC subgroups using a χ2 minimization technique (Proctor et al. 2004). This technique is based on the residuals between the Lick index measurements and the prediction from stellar population models. We used the flux-calibrated stellar population models of Lick indices presented by Thomas et al. (2011). The models were constructed with the ages ranging from 0.1 to 15 Gyr, the metallicities [Z/H] from −2.25 to + 0.67, and the α-element abundances, [α/Fe], from −0.3 to + 0.5. We measured the iron abundances by adopting the relation from Thomas et al. (2003), [Z/H] = [Fe/H] + 0.94 [α/Fe].
We used all Lick indices except CN1, CN2, Ca4227, and NaD to calculate the χ2 values. The stellar population model we adopted is calibrated with Galactic GCs that show an anomaly in the nitrogen abundance, which might not be applicable to typical simple stellar populations. The CN1, CN2, and Ca4227 are sensitive to nitrogen abundances, so they were excluded from the fitting. The NaD index is also excluded from the fitting because it is strongly affected by interstellar absorption. We determined the χ2 values with the remaining Lick indices, using iterative 2σ clipping. Finally, 14–21 line indices of each co-added spectrum are used for the fitting.
Because our individual GC spectra typically do not have the required S/N for a robust stellar population analysis, we performed our analysis on co-added spectra. We stacked the spectra of various GC subgroups to obtain mean ages, metallicities, and α-element abundances by taking the error-weighted mean values for each wavelength bin after strong sky emission line regions were masked. All spectra were de-redshifted prior to the co-addition. Figures 3(a) and (d) show flux-calibrated spectra of typical blue and red GC spectra, with S/N ∼ 9 Å–1 at 5000 Å, which is the median S/N for our sample. Panels (b), (c), (e), and (f) show stacked spectra of blue and red GCs with S/N of 30 and 80 as examples. The stacking analysis enables us to determine stellar population parameters from high-S/N spectra for GC subgroups. We have made the stacked spectra used in our analysis publicly available on Zenodo (10.5281/zenodo.6390374).
Download figure:
Standard image High-resolution imageFor the following analysis, the measurement errors for the ages, [Z/H], [Fe/H], and [α/Fe] were calculated using the bootstrap. For a given GC subgroup, we chose the same number of GCs in the group allowing replacement, stacked their spectra, and measured the stellar population parameters. We repeated this process 1000 times and took the 16th and 84th percentiles of the parameter distributions for lower and upper errors. The following results show the stellar population analysis for the GCs stacked using various criteria.
4. Results
4.1. Radial Gradients
Figures 4(a)–(d) show the mean age, [Z/H], [α/Fe], and [Fe/H] of the GC population in each radial bin as a function of major-axis distance from M87. The measurement values are listed in Table 1. The co-added spectrum of the GCs in each radial bin has an S/N of 80. Note that the numbers of blue and red GCs in each radial bin are comparable within , but the blue GCs dominate in the outer region, where its number fraction rises up to ∼80%.
Download figure:
Standard image High-resolution imageTable 1. Mean Ages, [Z/H], [α/Fe], and [Fe/H] for the Entire GC Sample
Rmaj Range (kpc) | a (kpc) | Age (Gyr) | [Z/H] (dex) | [α/Fe] (dex) | [Fe/H] (dex) |
---|---|---|---|---|---|
8–18 | 14 | ||||
17–38 | 28 | ||||
38–58 | 48 | ||||
58–75 | 65 | ||||
75–90 | 83 | ||||
90–108 | 99 | ||||
109–139 | 126 | ||||
139–175 | 154 | ||||
175–232 | 197 | ||||
234–376 | 279 | ||||
381–837 | 598 |
Note.
a Luminosity-weighted mean major-axis distance of GCs in each radial bin.Download table as: ASCIITypeset image
The GC populations in all radial bins are older than 8 Gyr, but the ones in the outskirts have formal mean ages about 2 Gyr younger than those in the inner region. However, it is well known that changing horizontal branch morphology can affect age estimates of integrated stellar populations (e.g., Lee et al. 2000; Koleva et al. 2008; Cabrera-Ziri & Conroy 2022), especially for metal-poor GCs. Tests with our integrated spectra show that age estimates using only Hδ (and masking all other Balmer lines) are 1–2 Gyr younger than those using Hβ alone for the GC bins at large radii, which is a sign that blue horizontal branch stars may be having a moderate effect on age estimates (Schiavon et al. 2004). Therefore, although we report this apparent gradient in age, we urge some caution in their interpretation.
We find a steep radial gradient in the mean metallicities of the GCs within a certain radius and a flat gradient beyond. We derive broken linear fits for these data with the form as follows:
The coefficient a is the value of [Z/H] at and beyond the scale distance R0 (in kpc), and the coefficient b is the slope within R0. The best-fit values we estimated are listed in Table 4. The entire GC sample shows a steep metallicity gradient with a slope of within (109 kpc) and a constant metallicity of [Z/H] = −1.34 beyond. This value of b means that, within R0, [Z/H] changes by −0.84 dex for a factor of 10 in radius. The break point of Rmaj ∼ 109 kpc is consistent with the distance at which the radial velocity dispersion of the GC system rises, found in Longobardi et al. (2018b), where the IGC population becomes dominant.
The inner radial bins have estimated ages that saturate at the limit of the models (15 Gyr). This model age zero-point issue could potentially affect our metallicity measurements. To test this, we artificially limit the model ages to be below 10 Gyr, finding that the estimated [Z/H] increases by only ∼0.1 dex. This moderate increase is likely representative of the potential effects of age saturation and is much smaller than the magnitude of the detected gradient.
We also confirmed that the metallicities of individual GCs with S/N > 20 are derived down to [Z/H] ∼ − 2 with uncertainties of ∼0.1. In addition, integrated metallicities of the Milky Way GCs, measured using Schiavon et al. (2005) spectra with the same method, are distributed below [Z/H] = −1.6 and show a good agreement with the values from the literature. We therefore conclude that the flattening of the metallicity gradient in Figure 4(b) is real and not due to a lack of metallicity sensitivity.
The α-element abundances also show a decreasing trend out to and then become constant at larger distances. We fitted the data with the same functional form as used for [Z/H] (Equation (3)), with R0 fixed to the value derived for [Z/H]. The coefficients a and b are derived to be a = 0.15 ± 0.04 and b = − 0.34 ± 0.11. Described in words, this value of b reflects the fact that GCs in M87's inner region are roughly twice as abundant in alpha-elements relative to iron compared to GCs at and beyond R0. We perform a similar analysis for the radial trend of [Fe/H]. Here we also find a significant gradient within R0. Although we again fixed R0 to 109 kpc, allowing R0 to vary freely produces a very similar value. The coefficients a and b for [Fe/H] are estimated to be and . The magnitudes of the slopes for [α/Fe] and [Fe/H] are comparable within the uncertainties, showing that they each contribute roughly equally to the [Z/H] gradient.
We also investigate the radial trends of the stellar populations of GC subpopulations where the co-added spectrum in each radial bin has an S/N of 70. The mean ages, [Z/H], [α/Fe], and [Fe/H] for blue and red GC samples are listed in Tables 2 and 3, respectively. The blue and red GCs in most radial bins have mean ages older than 10 Gyr (Figure 4(e)). Figure 4(f) shows the metallicity gradients of the blue and red GCs. The red GCs are clearly more metal-rich than the blue GCs, which indicates that the color corresponds well to the average metallicity.
Table 2. Mean Ages, [Z/H], [α/Fe], and [Fe/H] for the Blue GC Sample
Rmaj Range (kpc) | a (kpc) | Age (Gyr) | [Z/H] (dex) | [α/Fe] (dex) | [Fe/H] (dex) |
---|---|---|---|---|---|
9–25 | 19 | ||||
25–54 | 40 | ||||
54–77 | 65 | ||||
77–94 | 85 | ||||
94–112 | 103 | ||||
112–142 | 128 | ||||
141–181 | 158 | ||||
181–241 | 206 | ||||
241–421 | 301 | ||||
423–829 | 615 |
Note.
a Luminosity-weighted mean major-axis distance of blue GCs in each radial bin.Download table as: ASCIITypeset image
Table 3. Mean Ages, [Z/H], [α/Fe], and [Fe/H] for the Red GC Sample
Rmaj Range (kpc) | a (kpc) | Age (Gyr) | [Z/H] (dex) | [α/Fe] (dex) | [Fe/H] (dex) |
---|---|---|---|---|---|
8–29 | 18 | ||||
30–62 | 48 | ||||
61–103 | 79 | ||||
104–253 | 165 | ||||
259–837 | 488 |
Note.
a Luminosity-weighted mean major-axis distance of red GCs in each radial bin.Download table as: ASCIITypeset image
We fit linear relations for the metallicity and α-element abundance gradients of the blue subpopulation in the same way as they were derived for the entire sample. For the red GCs, we do a single linear fit, as there is no convincing evidence of a break point. The best-fit coefficients are listed in Table 4. The blue GCs have a steep total metallicity gradient with a slope of within (165 kpc) and a constant metallicity of [Z/H] = −1.56 beyond, which is a similar trend to that of the entire sample. On the other hand, the red GCs have a much shallower gradient in [Z/H], with a slope of . A value of [Z/H] ≈ − 0.7 describes the red GCs across a large range in radius.
Table 4. Linear Fit Results for Radial Gradients of Mean [Z/H], [α/Fe], and [Fe/H] for the GC Subsample
Sample | Rmaj Range (kpc) | R0 (kpc) | a[Z/H] | b[Z/H] | a[α/Fe] | b[α/Fe] | a[Fe/H] | b[Fe/H] |
---|---|---|---|---|---|---|---|---|
All | 8 − R0 | |||||||
R0 − 837 | 0 | 0 | 0 | |||||
Blue | 9 − R0 | |||||||
R0 − 829 | 0 | 0 | 0 | |||||
Red | 8–837 | 100a |
Note. The functional form of linear fits is a + b × log(Rmaj/R0).
a For red GCs, R0 is not a break point of segmented linear fits, but simply a reference point.Download table as: ASCIITypeset image
Figure 4(g) shows the radial trends of α-element-to-iron abundance ratios of the blue and red GCs. The blue GCs show a noticeable gradient with a slope of within (165 kpc), while the red GCs have a constant [α/Fe] value of ∼0.3 in all radial ranges with a radial gradient slope of . The mean [α/Fe] value of the outermost blue GCs is + 0.09, similar to the solar abundance.
Recently, Villaume et al. (2020) presented a stellar population study of M87 GCs based on full spectrum fitting of individual GC spectra. They found no significant radial gradient in [Fe/H] for either blue or red GCs. A direct comparison to this work is complicated by the fact that the sample of Villaume et al. (2020) covers a smaller range in radius (R < 150 kpc) and a gap in their data at 50 kpc that motivated their choice of break point in their fitting. Our [Fe/H] gradient for the blue GCs (−0.32 dex per decade in radius) is consistent with their findings, within the uncertainties, although a more homogenous analysis between the two data sets is certainly desirable in the future. Villaume et al. (2020) also find a strong negative gradient in the blue GCs for certain α-elements, like [Mg/Fe], which is similar to what we see in this work.
4.2. Azimuthal Variations
The Virgo Cluster core around M87 is known to be a place of active galaxy and cluster assembly that is likely not well mixed (e.g., Mihos et al. 2017). Our wide spatial coverage of GCs in this region allows us to investigate the possibility that there may also be azimuthal variations of mean ages and total metallicities of the GC population in addition to radial trends. We take ellipses centered on M87 that have a position angle of 155° and an ellipticity of 0.4 (Ferrarese et al. 2006; Janowiecki et al. 2010). We divide these ellipses into four sectors azimuthally: northeast (NE), northwest (NW), southeast (SE), and southwest (SW) (Figure 5). We further subdivide these sectors into four regions along the major-axis distance (, , , and ). Figure 5 shows the S/N of the co-added spectra of the blue and red GCs in subregions, ranging from 6 to 68. We measured the mean ages and total metallicities for 16 and 14 co-added spectra for the blue and red GCs, respectively.
Download figure:
Standard image High-resolution imageThe blue and red GCs in 26 out of 30 subregions are older than 10 Gyr, not showing specific trends. In addition to the sectors, we investigated the substructure that Powalka et al. (2018) identified as the region that relatively young GCs populate, located at 115 south of M87 with a diameter of ∼6'. We found 41 GCs in that region and measured their mean age with the co-added spectrum. They have a mean age of 9 Gyr, which is indeed a few Gyr smaller than the azimuthally averaged age of other GCs at the same distance from M87. Further kinematic studies of the GCs in this region will be helpful to reveal the origin of any younger GC population.
The two GC populations show radial gradients in metallicity for each sector. Figure 6 shows the radial trends in metallicity of the blue and red GCs in four sectors. The radial trends of the blue GCs in all sectors are consistent with those derived with the entire blue GC sample (Figure 4(e)). The red GCs also show shallow radial gradients in general, but the red GCs in the NW region of M87—the direction of M86 and M84—have slightly elevated mean metallicities compared to those in other sectors (Figure 6(b)). There is a possibility that this metal-rich population originates from massive galaxies in the northwest region of M87. However, the difference is not significant enough for us to draw a definite conclusion.
Download figure:
Standard image High-resolution image5. Discussion: The Assembly History of M87 and the Virgo Cluster Core
We have estimated ages and metallicities from co-added spectra of GC subpopulations and have found several notable features, especially in GC metallicity variations as a function of projected distance from M87. The GCs in M87 and the intracluster region are mostly old, potentially allowing us to trace the mass assembly history of M87 to early epochs (>10 Gyr). In this section, we discuss the influence of two different parameters related to the assembly history of the Virgo core region and M87 itself: the masses of GC progenitor host galaxies and the dependence on environment.
5.1. Gradients and GC Progenitor Galaxy Masses
The chemical abundances of stars and star clusters are intimately tied to the mass of the galaxy within which they are formed. The stellar mass–metallicity relation of galaxies has connected these quantities since its discovery (Lequeux et al. 1979), and the relation was confirmed with gas-phase metallicities, especially for star-forming galaxies (e.g., Tremonti et al. 2004; Mannucci et al. 2010; Andrews & Martini 2013). A series of studies have also clarified and extended the relation with stellar metallicities for Sloan Digital Sky Survey star-forming galaxies (Gallazzi et al. 2005) down to Local Group dwarf spheroidal and irregular galaxies (Kirby et al. 2013). Although the metallicities of the GCs are about 0.8–1 dex lower than those of the stars in a galaxy (Jordán et al. 2004; Lamers et al. 2017), we expect that the metallicities of the GCs also represent masses of progenitors where the GCs were formed. Therefore, one interpretation of the metallicity radial gradients seen in the blue and red GCs is that both blue and red GCs in the inner region of M87 were formed in more massive galaxies than those in the outskirts.
If the GCs in massive galaxies like M87 largely originate from merged and accreted galaxies that were disrupted (e.g., Côté et al. 1998) and GC metallicity correlates with progenitor host mass, then a metallicity gradient suggests a radial gradient in the stellar masses of disrupted progenitors. This is a natural consequence of the efficiency of dynamical friction as a function of the merger ratio (Amorisco 2019). Merged galaxies will spiral into the center owing to dynamical friction until they are disrupted. Not only do more massive progenitors experience dynamical friction more efficiently, but they can also resist complete disruption until they are deeper in the potential well of the larger galaxy. In this way, one expects more metal-rich stars and star clusters to be deposited closer to the center of the accreting galaxy, creating a metallicity gradient.
To check whether this explanation can work in practice, we use the Illustris cosmological simulation (Vogelsberger et al. 2014a, 2014b), where we analyze the GC system of a simulated galaxy cluster with a virial mass consistent with that of the Virgo Cluster as an example. Ramos-Almendares et al. (2020) tagged GCs to individual galaxies in the simulated galaxy at their infall time, following a Hernquist density profile (Hernquist 1990). The infall time here is defined as the epoch when a given galaxy can no longer be identified as the central galaxy of its own friends-of-friends group. They adopted different values of the characteristic density and scale length in the Hernquist profile for blue and red GC systems. The relation between the total mass in GCs and host galaxy mass (Harris et al. 2015; Hudson et al. 2015) is taken into account to assign masses to GCs. The detailed GC tagging analysis is described in Ramos-Almendares et al. (2020). In our analysis, we only study "orphan" GCs—GCs for which the progenitor host galaxies are destroyed at the present day.
We investigate where orphan GCs are located as a function of their progenitors' masses. The mean magnitude of our spectroscopic sample is mag, with 99% of the sample having mag. The latter limit corresponds to a mass of M = 4 × 105 M⊙ with an assumption of a GC mass-to-light ratio in the g-band of 2. We therefore only select the simulated galaxies with GC system masses higher than 4 × 105 M⊙ for a fair comparison. These galaxies have stellar masses ranging from 3 × 108 M⊙ to 1011 M⊙. Assuming that the mass of a GC is 4 × 105 M⊙, we randomly chose GC particles until their total mass was equal to the total mass of GCs in each galaxy. We repeated this particle sampling 100 times to determine the mean trend of progenitors' masses as a function of projected distance from the central galaxy.
Figure 7 shows the mean stellar masses at infall time of the progenitor galaxies of orphan GCs as a function of their present-day projected clustercentric distance. We find a steep gradient in the mean progenitor masses for blue orphan GCs located at 50 kpc < Rproj < 300 kpc, ranging from 1010.5 to 109.3 M⊙. This implies that the blue orphan GCs in the inner region of the central galaxy originate from more massive satellites that are more likely to sink further into the center than those in the outer region. We cannot constrain the gradient in the innermost region (Rproj < 50 kpc) because in this analysis we do not include the "in situ" population that has formed at the center and would be the most metal-rich population in the central galaxy. On the other hand, the red orphan GCs do not show a steep gradient in the mean progenitor masses because the progenitor masses of the red GCs have a much narrower range than those of the blue GCs. The red GCs are rarely formed in galaxies with M*, infall < 1010 M⊙. This trend matches well with the observed radial gradient in metallicity for the red GCs being shallower than that of the blue GCs.
Download figure:
Standard image High-resolution imageIn addition, we also found that the orphan GCs, especially blue ones, from progenitors with the earliest infall time (z ∼ 3.6) are located at the innermost region in this simulation, while those from recent infallers with a mean infall time of z ∼ 2.5 are distributed in the outer region. This hints at the observed age difference of a few Gyr between the inner and outer GCs.
While the simulations have their limitations, this example illustrates how including GCs in cosmological galaxy simulations, in a way that allows for studies of GC properties as a function of their spatial distributions (e.g., Ramos-Almendares et al. 2020; Reina-Campos et al. 2021; Chen & Gnedin 2022), makes available new diagnostics for galaxy assembly history.
5.2. GC Occupation Fraction and Gradient Flattening
In addition to the steeper radial metallicity gradient for blue GCs, we have found that the gradient flattens beyond ∼165 kpc (Figure 4(f)). Since the mass–metallicity relation extends to quite low metallicities and masses, and the gravitational processes that govern dynamical friction and galaxy disruption are also continuous, it might seem surprising that we should find any break in these radial trends. We show below that this could be related to the GC occupation fraction—the fraction of galaxies that contain GCs—as a function of galaxy mass.
Sánchez-Janssen et al. (2019) presented the star cluster occupation fraction of Virgo early-type galaxies as a function of galaxy stellar mass (see their Figure 6). All galaxies with M* > 109 M⊙ host GCs with mag (roughly the mean of the GC luminosity function), while only half of the galaxies with M* ∼ 107 M⊙ do. Because of the flux limit of our survey, we only observed relatively bright GCs, with a mean sample magnitude of mag, and with 99% of the sample having mag. The relevant GC occupation fraction for determining whether an accreted progenitor contributes GCs to our observed sample is therefore whether it contains a GC with mag. Using the GC data from the ACS Virgo Cluster Survey (Côté et al. 2004; Jordán et al. 2009), we determine the GC occupation fraction for Virgo galaxies as a function of GC magnitude and galaxy stellar mass. For galaxies with M* ∼ 1010 M⊙, the occupation fraction for GCs with mag is ∼80%, but for galaxies with M* ∼ 109 M⊙, the occupation fraction drops to ∼25%. Galaxies with even lower stellar masses are then not expected to host many GCs with mag. Therefore, we expect that beyond ∼100–150 kpc, while many lower-mass galaxies are being accreted and disrupted, only a small fraction of them can contribute GCs. The flattening of the gradients at large distances is a natural by-product of the fact that, below a certain mass, the low-mass galaxies that would host even more metal-poor GCs have few to no GCs. We see a sign of this in the simulated GC system shown in Figure 7(a). Beyond a certain distance, the mean progenitor masses of blue orphan GCs become constant because in this simulated galaxy cluster the galaxies disrupted in the outskirts have lower masses and do not host any GCs more massive than 4 × 105 M⊙. The location of the break point likely depends on the mass and assembly history of individual galaxy clusters. While the flux limit of our survey may affect the exact [Z/H] or [α/Fe] value at which the gradient flattens, we expect that the general trend would be similar even if one could observe fainter GCs, and this could be tested with future observations.
We present this as an alternative to a more traditional explanation for stellar population gradients, which is that they occur during dissipational star formation, where more rapid star formation and gas recycling toward the galaxy center cause more rapid chemical evolution. Evidence for steep metallicity gradients in "relic" galaxies with little apparent accreted populations (Martín-Navarro et al. 2015; Beasley et al. 2018) suggests that in situ star formation likely plays a role in these gradients, especially for more metal-rich populations. However, given the complexity of modeling a dissipational collapse, we find it interesting to show that the trends we see in the M87 GC system can be explained using a picture dominated by dry merging. This is especially relevant because it is expected that the intracluster light is dominated by stars acquired through dissipationless accretion (Purcell et al. 2007). It is likely, however, that dissipational collapse plays some role in forming the gradients we observe, especially in the inner regions, and it will be interesting to explore that possibility with future simulations.
5.3. [α/Fe] and Dependence on Environment
The [α/Fe] ratios allow us to understand the star formation timescales in the host environments. There is a well-known relation between α-element abundances and central velocity dispersion values for galaxies. The larger velocity dispersion the galaxies have, the higher [α/Fe] values they have (Thomas et al. 2005; Annibali et al. 2007; McDermid et al. 2015). However, for low-mass galaxies, the local density could be a more critical parameter than the mass. Liu et al. (2016) measured the α-element abundances of 11 low-mass dwarf galaxies in the Virgo Cluster and found that the [α/Fe] values of dwarf galaxies closer to M87 are higher than those distant from M87.
Similarly, we find moderate negative radial gradients in mean [α/Fe] values, especially for the blue GC sample, a population that likely originates from low-mass galaxies. This suggests that the blue GCs in the central region of M87 were formed in progenitors where rapid star formation occurred. Figure 8 shows the comparison of [α/Fe] gradients of GCs from this study and dwarf galaxies from Liu et al. (2016). The gradients are similar to each other, although the GCs are more centrally concentrated (in projected radius) than the Virgo dwarfs.
Download figure:
Standard image High-resolution imageFollowing the simulation result above, the blue GCs in the inner region of M87 might originate from relatively massive galaxies that became satellites of the central galaxy at high redshift. Intense star formation and rapid quenching occurred in these massive progenitors in the high-z universe, resulting in high [α/Fe] abundances in their stars and star clusters. Similarly, dwarf galaxies in dense environments have enhanced α-element abundances (Chilingarian et al. 2008; Smith et al. 2009; Liu et al. 2016). There is a possibility that the most metal-rich of the blue GCs in the inner region originate from dwarf galaxies that were quenched quickly by the dense environment (Mistani et al. 2016).
5.4. Local Substructures
The Virgo Cluster is a dynamically young galaxy cluster with several subclusters (Binggeli et al. 1985, 1987, 1993). Mihos et al. (2005, 2017) presented the deep imaging survey showing the diffuse intracluster light in the Virgo core region, indicating the ongoing interaction between galaxies. Related to this dynamical status, it is expected that the GCs can be stripped not only from low-mass dwarf galaxies but also from massive galaxies. The observed metal-rich GCs imply the ongoing interaction between M87 and nearby massive galaxies.
We detected a hint that the metal-rich GCs located in the NW region of M87 (Figure 6) are slightly more metal-rich than in other areas of the Virgo core. This direction is toward massive early-type galaxies, M86 and M84. We suspect that these metal-rich GCs might be associated with M84 because we exclude the GCs that have lower radial velocities corresponding to the systemic velocity of M86 (Figure 1(a)). We defer to a future paper an investigation of the many GCs associated with known substructures, particularly the M86 group.
6. Summary
We present a stellar population analysis of the GCs in the central region of the Virgo Cluster, especially around M87, using the MMT/Hectospec spectra. We select 692 GCs that belong to the main body of the Virgo Cluster to investigate the stellar populations of GCs, excluding the GCs associated with substructures. The GCs are divided into several groups with various criteria for colors, major-axis distances, and position angles. We measured mean ages, total metallicities, [Fe/H], and α-element abundances of the GC subpopulations using their co-added spectra based on Lick indices. The main findings of this study are summarized as follows.
- 1.The GCs as a whole follow a two-component radial metallicity ([Z/H]) trend, with a steep gradient following a power-law slope of within (109 kpc), and a flat gradient with a mean metallicity of [Z/H] ∼ − 1.34 beyond. The trends for [Fe/H] and [α/Fe] are similar, with inner region slopes of and -0.34 ± 0.11 and outer flattenings at and [α/Fe] = 0.15 ± 0.04, respectively.
- 2.When analyzing gradients in GC color subpopulations (blue and red), we find that the two-component trend is clearly seen only in the blue GC system, while the red GCs are best fit by a single component. For [Z/H], the blue GCs exhibit an inner region power-law slope of within (165 kpc) and a mean metallicity of [Z/H] = − 1.56 ± 0.05 beyond. The red GCs show a much shallower radial gradient in [Z/H] than do the blue GCs, with a power-law slope of and no obvious flattening at large radius.
- 3.The [Fe/H] and [α/Fe] gradients also behave differently for the blue and red GC subpopulations. The blue GCs follow the same two-component trend, with inner region slopes of and and outer flattenings at and , respectively. These results show that the [Z/H] gradient seen in the blue GCs is driven in roughly equal portions by [Fe/H] and [α/Fe]. By contrast, the red GCs exhibit only a weak gradient in [Fe/H] (slope of ) and no gradient in [α/Fe].
- 4.We investigated the azimuthal variation of GC metallicities and detected a hint of local substructures toward the massive early-type galaxy M84, consisting of metal-rich GCs.
- 5.The GCs in the Virgo core are mostly old (>8 Gyr) regardless of their colors and locations. The outer GCs are formally about 2 Gyr younger than inner GCs, but there is a possibility that the presence of blue horizontal branch stars is causing this apparent age gradient.
We interpreted the radial gradient of GCs in metallicity and [α/Fe] in the context of the assembly history of M87 and the Virgo Cluster core. We used the Illustris cosmological simulations to investigate masses of disrupted GC progenitors in a simulated cluster.
- 1.According to the mass–metallicity relation of galaxies, the more metal-rich blue GCs at smaller galactocentric distances were probably formed in more massive galaxies. We expect that these more massive progenitors were better able to survive disruption as they approached the galaxy center due to dynamical friction. From the simulation, we found that the GCs from massive progenitors sink further into the central galaxy, while GCs from lower-mass progenitors are located preferentially at larger distances.
- 2.The flattening of metallicity gradients of blue GCs at large distances (R ≳ 100 kpc) is a natural consequence if the lower-mass dwarfs that could contribute GCs with even lower metallicity actually host few or no GCs.
- 3.We suggested that the dense environment around M87 could result in the high [α/Fe] values of the GCs in the central region, mirroring the environmental dependence of the [α/Fe] of Virgo dwarf galaxies (Liu et al. 2016).
Overall, we suggest a scenario where more massive progenitors accreted into the central galaxy (M87) contribute a significant fraction of the GCs in the inner region. Future simulation studies on various galaxy clusters will be helpful to understand the relation between the assembly history of a host galaxy cluster and chemical properties of the GCs residing in it.
Observations reported here were obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution. MMT telescope time was granted, in part, by NOAO, through the Telescope System Instrumentation Program (TSIP). TSIP is funded by NSF. E.W.P. thanks the staff of the MMT Observatory for their professional support over multiple observing runs. E.W.P. also thanks the Harvard-Smithsonian Center for Astrophysics and its Institute for Theory and Computation for hosting him during an unexpectedly long pandemic visit. We thank Daniel Fabricant for sharing his flux calibration routine for MMT/Hectospec data. We also thank the anonymous referee for useful comments.
This paper is based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/IRFU, at the Canada–France–Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.
C.L. acknowledges support from the National Natural Science Foundation of China (NSFC, grant Nos. 12173025, 11673017, 11833005, 11933003), 111 project (No. B20019), and Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education. A.L. is supported by Fondazione Cariplo, grant No. 2018-2329. A.L., T.H.P., and P.-A.D. acknowledge support from grant ANR-19-CE31-0022(POPSYCLE) of the French Agence Nationale de la Recherche. L.V.S acknowledges support from the NASA ATP 80NSSC20K0566 and NSF CAREER 1945310 grants. A.J. and S.E. acknowledge support from ANID—Millennium Science Initiative—ICN12_009. M.G.L. was supported by the National Research Foundation grant funded by the Korean Government (NRF-2019R1A2C2084019). H.S.H. acknowledges the support by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2021R1A2C1094577). S.L. acknowledges the support from the Sejong Science Fellowship Program through the National Research Foundation of Korea (NRF-2021R1C1C2006790).