Galaxy archaeology for wet mergers: Globular cluster age distributions in the Milky Way and nearby galaxies

Abstract

Context. Identifying past wet merger activity in galaxies has been a longstanding issue in extragalactic formation history studies. Gaia’s 6D kinematic measurements in our Milky Way (MW) have vastly extended the possibilities for Galactic archaeology, leading to the discovery of a multitude of early mergers in the MW’s past. As recent work has established a link between younger globular clusters (GCs; less than about 10–$11\text{\,}\mathrm{Gyr}$ old) and wet galaxy merger events, the MW provides an ideal laboratory for testing which GC properties can be used to trace extragalactic galaxy formation histories.

Aims. To test the hypothesis that GCs trace wet mergers, we relate the measured GC age distributions of the MW and three nearby galaxies, M 31, NGC 1407, and NGC 3115, to their merger histories and interpret the connection with wet mergers through an empirical model for GC formation.

Methods. The GC ages of observed galaxies are taken from a variety of studies to analyze their age distributions side-by-side with the model. For the MW, we additionally cross-match the GCs with their associated progenitor host galaxies to disentangle the connection to the GC age distribution. For the modeled GCs, we take galaxies with similar GC age distributions as observed to compare their accretion histories with those inferred through observations.

Results. We find that the MW GC age distribution is bimodal, mainly caused by younger GCs (10–$11\text{\,}\mathrm{Gyr}$ old associated with Gaia-Sausage/Enceladus (GSE) and in part by unassociated high-energy GCs. The GSE GC age distribution also appears to be bimodal. We propose that the older GSE GCs (12–$13\text{\,}\mathrm{Gyr}$ old) were accreted together with GSE, while the younger ones formed as a result of the merger. For the nearby galaxies, we find that clear peaks in the GC age distributions coincide with active early gas-rich merger phases. Even small signatures in the GC age distributions agree well with the expected wet formation histories of the galaxies inferred through other observed tracers. From the models, we predict that the involved cold gas mass can be estimated from the number of GCs found in the formation burst.

Conclusions. Multimodal GC age distributions can trace massive wet mergers as a result of GCs being formed through them. From the laboratory of our own MW and nearby galaxies we conclude that the ages of younger GC populations of galaxies can be used to infer the wet merger history of a galaxy.

1 Introduction

Over the course of galaxy formation and evolution, in situ formed structures mix with accreted matter, concealing their origin and thereby the formation history of the galaxy. Through precise measurements of star and gas properties, such as their distribution in space, their velocity, their chemical compositions, and the stellar ages, it is possible to disentangle many of the individual clues on the formation history. This is the aim of galaxy archaeology (e.g., Freeman & Bland-Hawthorn 2002; Binney 2013; Helmi 2020). In particular, stellar structures and overdensities, such as stellar streams and other tidal features, can reveal details on a galaxy’s history, where especially stellar streams are generally the remains of a tidally disrupted satellite galaxy. This has been done extensively for the Milky Way (MW; e.g., Helmi et al. 1999; Belokurov et al. 2006, 2007; Bell et al. 2008; Shipp et al. 2018), aided particularly by the Gaia mission (Gaia Collaboration et al. 2016, 2018a, 2023) in recent years with 6D phase-space data (e.g., Helmi et al. 2018; Helmi 2020; Prudil et al. 2022; Malhan et al. 2022; Ruiz-Lara et al. 2022). It is also possible to detect tidal features for other galaxies in a more limited way through photometric and integral field unit (IFU) observations. Such identified structures have also been connected to the merger history of their host galaxies in observations (e.g., Bílek et al. 2020, 2023; Chandra et al. 2023) and simulations (e.g., Bullock & Johnston 2005; Johnston et al. 2008; Amorisco 2015; Hendel & Johnston 2015; Pop et al. 2018; Karademir et al. 2019; Valenzuela & Remus 2024). However, all of these are tracers for stellar-dominated merger events and cannot trace the gas that has been accreted through such mergers.

For our own Galaxy, data are available in unprecedented detail, including measured proper motions. By detecting kinematic substructures of stars clustered in phase space, it has been possible to identify a number of past and ongoing mergers for the MW. The Sagittarius Dwarf Spheroidal Galaxy was the first merger discovered with the MW through positional and line-of-sight kinematic data by Ibata et al. (1994). Using Hipparcos data and line-of-sight distances and velocities, Helmi et al. (1999) discovered substructures in the inner halo now known as the Helmi streams. Through the Gaia mission and the 6D kinematic data that were obtained for stars in the MW, a large number of halo stars were found to have distinct kinematics in phase space, which were attributed to a major merger event that is expected to have taken place around $10\text{\,}\mathrm{Gyr}$ ago, Gaia-Sausage/Enceladus (GSE; Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Mackereth et al. 2019). It is the last major merger that the MW experienced, as well as the most massive one, forming a large part of the stellar halo. There are also other interpretations of the kinematic signatures, however, that the local halo cannot be dominated by GSE alone, but is composed of the remains of multiple mergers over a prolonged time (e.g., Donlon et al. 2020, 2022; Donlon & Newberg 2023). Finally, Myeong et al. (2019) find a second group of halo stars kinematically and chemically different from GSE stars with overall retrograde motions, which they attribute to a separate merger event, referred to as Sequoia. Further groups of stars in phase space have been found, which are likely the debris of disrupted galaxies that fell into the MW, though the connection to specific merger events is not yet clear (see Dodd et al. 2023 and Horta et al. 2023 for a current overview of structures found in phase space).

Through Gaia data, globular clusters (GCs) in the MW have also been used to further disentangle the Galactic formation history. Based on their 6D phase-space properties and their age-metallicity relations, some studies have linked the GCs with their likely progenitor host galaxies (Myeong et al. 2018; Massari et al. 2019; Forbes 2020; Horta et al. 2020; Callingham et al. 2022; Belokurov & Kravtsov 2024; Chen & Gnedin 2024b), such as to the MW itself as in situ GCs, or to some of the inferred accreted galaxies. These studies also reveal unassociated groups of GCs with low and high orbital energies, where it is proposed that a group of unassociated low-energy GCs are part of a further past merger event (Massari et al. 2019; Forbes 2020; Callingham et al. 2022). Similar structures in phase space with overlapping properties are identified through different methods (e.g., Kruijssen et al. 2019; Forbes 2020; Horta et al. 2021, 2023). However, it has also been shown that GCs can migrate in phase space over time, potentially making it difficult to disentangle the origin of GCs based on their phase-space properties alone (Pagnini et al. 2023). An overview of the positions on the sky and in phase space of the streams and GCs detected in the MW can be found in the figures of Riley & Strigari (2020), Malhan et al. (2022), and Mateu (2023).

Because of their intrinsic brightness, GCs have also been used as tracer populations in the outskirts of other galaxies to study their mass distribution, kinematics, and formation history (e.g., Peng et al. 2004; Foster et al. 2011; Coccato et al. 2013; Pota et al. 2013; Zhu et al. 2014; Dolfi et al. 2021; Veršič et al. 2024). Through their old age, GCs in particular have experienced a large part of a galaxy’s history, making them valuable tracers for past merger events. However, the formation of GCs themselves is still poorly understood. For this reason, models and simulations have been developed to help constrain the details of their formation process. Highly-resolved hydrodynamic simulations help study the resolved formation of individual GCs (e.g., Kravtsov & Gnedin 2005; Lahén et al. 2019; McKenzie & Bekki 2021), and sub-grid models for GCs applied to isolated or cosmological simulations allow one to follow GC properties and their spatial distribution through time, making comparisons with observations of nearby galaxies possible (e.g., Bekki et al. 2005; Kruijssen et al. 2011; Pfeffer et al. 2018; Chen & Gnedin 2022; De Lucia et al. 2024; Doppel et al. 2023; Reina-Campos et al. 2023; Chen & Gnedin 2024a). Finally, empirical and semi-analytic GC formation models applied to cosmological merger trees of galaxies enable one to test the parameter space of a limited number of free parameters for a large number of galaxies (e.g., Beasley et al. 2002a; Choksi et al. 2018; El-Badry et al. 2019; Valenzuela et al. 2021; Chen & Gnedin 2023). Such models can lead to a better understanding of the statistical properties of GC formation that are necessary to reproduce GC properties and relations as they are observed today (e.g., Spitler & Forbes 2009; Harris et al. 2015, 2017; Forbes et al. 2018; Burkert & Forbes 2020).

The stars of GCs can be individually observed in the MW, such that reasonably good measurements of their ages can be determined through the color-magnitude diagram (CMD), which has been done for various GC subsamples in the MW (e.g., Salaris & Weiss 1998; Dotter et al. 2008; Marín-Franch et al. 2009). For extragalactic GCs, the ages have much larger uncertainties associated with them because only the integrated GC properties can be measured. For this reason, stellar population models and evolutionary tracks have to be used for extragalactic systems to determine the ages. Ages have been estimated from photometric studies of the integrated light of extragalactic GCs (e.g., Chies-Santos et al. 2011; Georgiev et al. 2012; de Brito Silva et al. 2022), a technique which will also become increasingly important with the upcoming generation observatories such as Euclid, the Nancy Grace Roman Space Telescope, and the Vera C. Rubin Observatory (e.g., Lançon et al. 2021; Dage et al. 2023; Usher et al. 2023). Even for integrated spectroscopic measurements of extragalactic GCs, there are biases towards younger ages compared to the CMD method. This bias arises from hot horizontal branch (HB) stars, whose presence is degenerate with younger stars (e.g., Worthey 1994; de Freitas Pacheco & Barbuy 1995; Beasley et al. 2002b; Cabrera-Ziri & Conroy 2022) and even causes issues when using state-of-the-art models (e.g., Perina et al. 2011; Gonçalves et al. 2020). This has been done for GCs in galaxies in the Local Group (e.g., Beasley et al. 2005; Schiavon et al. 2013; Wang et al. 2021) and for selected nearby galaxies (e.g., Usher et al. 2019).

In this work, we used a recent empirical GC formation model with two formation pathways (Valenzuela et al. 2021; Valenzuela 2023, the first pathway forms GCs in small halos, the second forms GCs in gas-rich mergers) to study in what way its second pathway of forming GCs in gas-rich wet mergers (e.g., Ashman & Zepf 1992) can help shed light on the formation history of the MW and other nearby galaxies. The model has been shown to agree well with the observed numbers of GCs in galaxies from dwarf to galaxy cluster masses, where a linear relation has been found to exist between the number of GCs and the dark matter (DM) halo virial mass (e.g., Blakeslee et al. 1997; Harris et al. 2017; Forbes et al. 2018), as well as with GC age distributions of the MW and nearby galaxies. We introduce the GC model and observational data in Sect. 2. In Sect. 3, we present a bimodal feature found in the observed GC age distribution of the MW and how it could be related to the predictions of the GC model. We then test and discuss these predictions in detail for the MW in Sect. 4 and for other nearby galaxies in Sect. 5. Finally, we summarize and conclude the results in Sect. 6.

2 Data & method

In the following, the empirical GC model and the GC observational data used in this work are presented. The main property studied is the GC age distribution of galaxies.

2.1 Globular cluster model

In this work, we used the empirical GC formation model introduced by Valenzuela et al. (2021), which builds on previous models and investigations by Boylan-Kolchin (2017), Choksi et al. (2018), and Burkert & Forbes (2020). The model employs two formation pathways for GCs: The first, the small halo pathway, forms GCs in small haloes as soon as a halo’s virial mass surpasses a given threshold value, $M_{\mathrm{seedGC}}$. With equal probability, 0, 1, or 2 GCs are formed. The second, the wet merger pathway, is the formation pathway introduced by Choksi et al. (2018) and triggers GC formation when the relative halo mass accretion rate surpasses a given threshold value, $A_{\mathrm{min}}$. The formed number of GCs is then determined by converting the available cold gas mass $M_{\mathrm{gas}}$ to a total GC mass via a conversion factor, $\eta_{\mathrm{GC}}$ (Kravtsov & Gnedin 2005; Li & Gnedin 2014; Choksi et al. 2018; Valenzuela et al. 2021):

$$M_{\mathrm{GC}}=$1.8\text{\times}{10}^{-4}$\eta_{\mathrm{GC}}M_{\mathrm{gas}}.$$ (1)

By assuming a cluster initial mass function of

$$\frac{dN}{dM}\propto M^{-2},$$ (2)

the expected number of formed GCs is obtained as (combining eqs. 3, 6, and 7 of Valenzuela et al. 2021)

$$\langle N\rangle=\exp\bigg{(}W\Big{(}\frac{$1.8\text{\times}{10}^{-4}$\eta_{\mathrm{GC}}M_{\mathrm{gas}}}{M_{\mathrm{min}}}\Big{)}\bigg{)}-1,$$ (3)

where $W$ is the Lambert $W$ function, $\eta_{\mathrm{GC}}=0.5$ for the best-fitting model, and $M_{\mathrm{min}}=${10}^{5}\text{\,}\mathrm{M_{\odot}}$$ is the minimum mass that a GC needs at formation time to survive for a few $\mathrm{Gyr}$ (Li & Gnedin 2014; Choksi et al. 2018). For more details on the models, see Valenzuela et al. (2021). Note that recent work by Chen & Gnedin (2023) now uses a smaller value of $M_{\mathrm{min}}=${10}^{4}\text{\,}\mathrm{M_{\odot}}$$, although they note that GCs with low initial masses of for example ${10}^{4}\text{\,}\mathrm{M_{\odot}}$ have an estimated lifetime of around $1\text{\,}\mathrm{Gyr}$ at $3\text{\,}\mathrm{kpc}$ distance of the center of a MW-mass galaxy. The model only tracks the numbers of GCs per galaxy and their formation times, but does not include metallicities. This limits the comparison with observations to only the GC ages, though the available measured metallicities can be used as indicators for the formation sites of the observed GCs. In contrast, for the modeled GCs this information is already known.

The GC model was applied to the merger tree of a DM-only simulation of side length $30\text{\,}\mathrm{Mpc}$ with a DM particle mass of $m_{\mathrm{DM}}=$7.90\text{\times}{10}^{6}\text{\,}\mathrm{M_{\odot}}$$ that was run with the TreePM code Gadget-3 (Springel 2005). The empirical model emerge (Moster et al. 2018) provides the model with the baryonic matter content per galaxy. Because the model only tracks the number of GCs in each galaxy and at what times they were formed, the model parameters were fit to match observations of the numbers of GCs, since those are available for a sufficiently large sample: the GC numbers are taken from Burkert & Forbes (2020). The GC age distributions are consistent with those found by Usher et al. (2019) for the MW and three SLUGGS galaxies (SAGES Legacy Unifying Globulars and GalaxieS; Brodie et al. 2014), and the fractions of GCs formed through the wet merger pathway agree with the red GC fractions measured by Harris et al. (2015). For more information on the comparability to observations and how the model parameters affect the resulting GC properties, see Valenzuela et al. (2021).

2.2 Globular cluster observational data

A variety of observed GC age measurements from the literature are used in this work for selected galaxies in the Local Universe. These use different methods to obtain the GC ages and are presented in the following. For all of the measured samples, it is important to keep in mind that while the absolute ages have large uncertainties and are difficult to measure even in the MW itself (e.g., Ying et al. 2023, for M 92), within a given GC sample they are subjected to the same systematic uncertainties, resulting in much more precise relative ages. This is important for the study of GC age distributions, where the features in the distribution itself are given by the relative ages as opposed to the absolute ones.

2.2.1 Milky Way

The MW is the only galaxy besides its own satellites for which it is currently possible to obtain accurate CMDs of its GCs. This allows for a much more exact determination of their ages and has been done in multiple studies for different sized samples of GCs. In this work, we considered three of these as compiled by Kruijssen et al. (2019), and additionally investigated the mean ages of those three studies in Sect. A.1:

• Forbes & Bridges (2010) with a sample of 92 GCs,

• Dotter et al. (2010, 2011) with a sample of 68 GCs,

• VandenBerg et al. (2013) with a sample of 54 GCs,

• Kruijssen et al. (2019) with a sample of 96 GCs, which contains the mean GC ages from the previous three studies.

The sample from Forbes & Bridges (2010) is based on a number of previous age measurement studies (Salaris & Weiss 1998; Bellazzini et al. 2002; Catelan et al. 2002; De Angeli et al. 2005; Carraro et al. 2007; Dotter et al. 2008; Carraro 2009; Marín-Franch et al. 2009), of which 64 GCs were measured using the Advanced Camera for Surveys (ACS) from the Hubble Space Telescope (HST) through the ACS survey for Galactic GCs (Sarajedini et al. 2007; Marín-Franch et al. 2009) to obtain relative ages. They were normalized to absolute ages with the Dartmouth models of Dotter et al. (2007). While that sample is restricted to the inner $20\text{\,}\mathrm{kpc}$ of the MW, the age measurements of further GCs were supplemented from the other works. The sample from Dotter et al. (2010, 2011) is for the most part also based on the GCs measured through the ACS survey of Galactic GCs using the photometric catalog from Anderson et al. (2008), and the remaining GCs were observed with further HST/AST measurements. Lastly, the GC age measurements from VandenBerg et al. (2013) were computed from the same photometric catalog of the ACS survey of Galactic GCs as Marín-Franch et al. (2009) used, but employing the stellar evolutionary tracks from VandenBerg et al. (2012). It should be noted that while not all these objects are necessarily actual GCs, but in part also nuclear star clusters, metal complex clusters, or a combination thereof (e.g., McKenzie et al. 2022), we refer to all of them simply as GCs in this work.

In addition to the three CMD age samples of the MW GCs, we also included two GC age samples obtained through integrated measurements, as this is also what one is restricted to for other galaxies:

• Usher et al. (2019) with a sample of 61 GCs. Their measurements are obtained through combining photometry and spectroscopy, to which stellar population models are fitted using a Markov chain Monte Carlo (MCMC) method.

• Cabrera-Ziri & Conroy (2022) with a sample of 32 GCs, of which we removed the 3 spurious young GCs (see their section 6.1, where they detail how the spectral fit residuals indicate whether the best fit for a young GC is real or spurious, which is a result of the simple single population HB star model that they use). These have also been shown to be much older from CMD measurements. Compared to Usher et al. (2019), Cabrera-Ziri & Conroy (2022) used spectroscopy only, but additionally took the hot horizontal branch (HB) stars into account in their modeling of the integrated stellar population measurements.

For the MW, additional GC properties can be determined that are not possible to obtain for other galaxies at the moment. Six-dimensional phase space measurements have been made available for many MW GCs through Gaia (Gaia Collaboration et al. 2018b; Vasiliev 2019), which Massari et al. (2019) used to assign the likely origin of the individual GCs. Their list of progenitors for the GCs consist of the MW itself (i.e., in situ formed GCs in the disk or bulge), the GSE galaxy, the Sagittarius dwarf, the Helmi Streams, the Sequoia galaxy, and unassociated high- and low-energy GCs. Forbes (2020) used the age-metallicity relation (AMR) to improve these progenitor assignments and propose that the unassociated low-energy GCs belong to a single progenitor dwarf satellite, which they dub Koala and is likely related to or overlaps with Kraken (Kruijssen et al. 2019) and Heracles (Horta et al. 2021, 2023). Further work was done by Callingham et al. (2022), who used a chemo-dynamical model and hydrodynamical simulations of MW-like galaxies to associate the GCs with their progenitor hosts. Their assumed accretion events largely align with those used by Massari et al. (2019) and Forbes (2020). Chen & Gnedin (2024b) have recently applied clustering techniques to chemical, spatial, kinematic, and age properties of the GCs to associate the GCs with the MW, GSE, Sagittarius, and other ex-situ-formed GCs. It should be noted that due to the ongoing observations and work on this topic, this is a rapidly evolving field. The identifiers of the GCs were available to us for the CMD age samples and for the sample of Cabrera-Ziri & Conroy (2022), such that we could cross-match the ages for those four samples to the progenitor assignments. For this work, we use the assignments made by Forbes (2020), though using the assignments from Massari et al. (2019), Callingham et al. (2022), or Chen & Gnedin (2024b) do not change the statistical findings presented in this work when we did the analysis using those instead (also see Sect. A.2 for an analysis using the associations from Callingham et al. 2022 and Chen & Gnedin 2024b).

2.2.2 M 31

As the nearest more massive galaxy to the MW, the Andromeda galaxy (M 31) is an ideal galaxy for which GCs can be identified and analyzed, since observations have much better resolution and better signal-to-noise values than for more distant galaxies. For M 31, we used two different studies for the GC ages:

• Wang et al. (2021) with a sample of 343 clusters, of which we used the 293 old GCs ($t_{\mathrm{age}}>$1.5\text{\,}\mathrm{Gyr}$$) in this work,

• Usher et al. (2024) with a sample of 290 clusters,

• Cabrera-Ziri et al. (in prep.) with a sample of 286 GCs, of which we used the 136 GCs whose ages are sufficiently constrained.

The sample from Wang et al. (2021) was observed with the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Cui et al. 2012; Luo et al. 2015). The GC ages were then determined based on the obtained integrated spectra and multi-band photometry from Beijing-Arizona-Taiwan-Connecticut (BATC; Ma et al. 2015) and Sloan Digital Sky Survey (SDSS; Peacock et al. 2010). For the old clusters, they fit the spectra using an empirical stellar spectra library and the measured colors using an MCMC method. Their ages are all younger than $10\text{\,}\mathrm{Gyr}$, however, which is likely a consequence of their underlying models and almost certainly not actually the case for the majority of M 31’s GCs (e.g., Beasley et al. 2005; Caldwell et al. 2011; Schiavon et al. 2013). Their ages should therefore be taken with caution. The GC age determinations from Usher et al. (2024) were also made from photometric and spectroscopic measurements, where the photometry was obtained with SDSS (mostly from Peacock et al. 2010) or from the Pan-Andromeda Archaeological Survey (PAndAS) Canada France Hawaii Telescope MegaCam for the halo GCs not covered by SDSS (Huxor et al. 2014). Usher et al. (2024) used the same method as Usher et al. (2019), for which they fit the stellar spectra using an MCMC method.

The integrated spectra measurements of Cabrera-Ziri et al. (in prep.) use the same method as applied for the MW GCs presented by Cabrera-Ziri & Conroy (2022). Their GC sample is selected from the inner halo of M 31, of which 150 have only a lower bound on their age, leaving 136 GCs with sufficiently constrained ages to study their distribution. The reason for the large number of lower age bounds is that their method is able to determine that the metal-poor GCs with a horizontal branch are very old, but not what exact age they have (${>}$9\text{\,}\mathrm{Gyr}$$ in all cases, and ${>}$12.5\text{\,}\mathrm{Gyr}$$ is the median lower bound) due to degeneracies between the ages and horizontal branch properties. This also means that the resulting age distribution has a selection bias towards metal-rich GCs, which means that the metal-poor GCs typically accreted through smaller galaxies are removed. This should be kept in mind, though we believe that the consequences for this work are not severe since we focus on GCs formed during mergers, which produce younger GCs at generally higher metallicities. GCs formed through such mergers tend to be younger by on average 2–$3\text{\,}\mathrm{Gyr}$, having their formation peak at an age of roughly $10\text{\,}\mathrm{Gyr}$ (see fig. 2 of Valenzuela 2023).

2.2.3 NGC 3115 and NGC 1407

The age measurements for 116 GCs in NGC 3115 and 213 GCs in NGC 1407 were made by Usher et al. (2019) by combining photometry and spectroscopy to fit the stellar spectra using an MCMC method. The spectral data were obtained through the SLUGGS survey (Brodie et al. 2014) and the photometric data are from Arnold et al. (2011) using the Suprime-Cam of the Subaru telescope.

3 Age dating wet mergers with GCs in the Milky Way

Having the observational data, we investigated some reoccurring features in the GC age distributions with the aim of verifying if such features can be explained by the model. In the following, we present our initial findings from the observations and what predictions the model makes with respect to these.

3.1 GC age distribution observations in the Milky Way

As the galaxy for which the best data exists, we first considered the GC age distribution of the MW. In two of the available samples, the ages appear to have a bimodal or even multimodal distribution, which are shown in Fig. 1. For both Forbes & Bridges (2010) and Usher et al. (2019), the main peak of the distribution is at $13\text{\,}\mathrm{Gyr}$. There is also a slight second peak in the distribution of Forbes & Bridges (2010) at around $11\text{\,}\mathrm{Gyr}$, while for Usher et al. (2019) there are two small peaks between $8\text{\,}\mathrm{Gyr}10\text{\,}\mathrm{Gyr}$. Since the data from Forbes & Bridges (2010) are composed of multiple different studies, we show the age distribution only of the largest of the underlying GC age studies from Marín-Franch et al. (2009) in Sect. A.3 (it includes 64 of the 92 GCs), which also features the age bimodality. Thus, the bimodality seen for Forbes & Bridges (2010) is not the result of combining multiple studies with different systematic uncertainties. Studies of the age-metallicity relationship for MW GCs have shown that there are multiple branches of GCs corresponding to the in situ formed GCs and different progenitor galaxies that fell into the MW (e.g., Kruijssen et al. 2019; Forbes 2020), which could be what is seen as a multimodal age distribution. The peaks in the age distributions align well with the estimated infall times of Sagittarius (8–$9\text{\,}\mathrm{Gyr}$ ago) and Sequoia (${\sim}$10\text{\,}\mathrm{Gyr}$$ ago; Forbes 2020), and of the GSE merger event, the last major merger that the MW experienced (Helmi et al. 2018; Haywood et al. 2018): around $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ ago according to Belokurov et al. (2018) and around $10\text{\,}\mathrm{Gyr}$ ago according to Helmi et al. (2018). With an estimated merger mass ratio of 1:4 (Helmi et al. 2018; Gallart et al. 2019), GSE is expected to have been a gas-rich major merger that triggered a starburst in the disk of the MW (e.g., Helmi 2020; Ciucă et al. 2023). There is also evidence for this from stellar population measurements, where Gallart et al. (2019) found a clear peak of high star formation at around $9.5\text{\,}\mathrm{Gyr}$ ago.

Figure 1: Globular cluster age distributions in the MW from Forbes & Bridges (2010) and Usher et al. (2019), with sample sizes of 92 and 61 GCs, respectively. The boxy lines are the actual histograms, while the smooth curves show the distributions smoothed by the measurement uncertainties. These are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations. The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018).

Of course, the findings from the observed age distributions are accompanied by some caveats: First of all, GC age measurements are subject to large uncertainties, especially for the older ages (see the solid lines in Fig. 1, which are the distributions smoothed by Gaussian kernels given by the measurement uncertainties). Still, a large part of these uncertainties are systematic effects, such that the relative ages can still be trusted more than individual absolute ages. Second, the sample sizes are not large enough to make any kind of statistically significant statements, in particular for the few GCs that make up the minor peaks in the age distributions. Third, bringing together the estimated time of the GSE merger event and the time of the GC age distribution peaks is far from having established a causal relation between the two. However, since the features are present in the data and the time also aligns with that of the GSE merger event, GC formation models can be employed to investigate if there is a theoretical basis for a causal connection.

3.2 Model predictions

The model introduced by Valenzuela et al. (2021) consists of two GC formation pathways, of which one allows GCs to form through wet mergers containing a sufficient amount of cold gas. To study what the model predicts for galaxies such as the MW, we first extracted MW-like galaxies from the simulation to which the GC model was applied. For this, we selected MW-mass galaxies by their virial mass of $M_{\mathrm{vir}}=1$–$2\text{\times}{10}^{12}\text{\,}\mathrm{\mathnormal{M}_{\mathrm{vir}}}$, which applies to 21 simulated galaxies. Out of these, there are two with GC age distributions that best match the distribution measured by Usher et al. (2019) for the MW in terms of their cumulative distributions, which is shown in the top panel of Fig. 2. The analogs were selected using the measurements by Usher et al. (2019) instead of one of the CMD measurements since the model was originally calibrated by Valenzuela et al. 2021 to be in agreement with the ages measured by Usher et al. 2019 to have a comparison with multiple galaxies. A recalibration of the model to the CMD age distributions is unsuitable because it would surpass the scope of this work and would only provide one single galaxy with a sufficient number of CMD-measured GC ages. The conclusions drawn in the following are still applicable to the model in general, independent of the age calibration that was used. We show the GC age distributions and the accretion histories of the other 19 MW-mass galaxies in Appendix B.

Both of the simulated galaxies have roughly bimodal GC age distributions with a minor peak at around $9\text{\,}\mathrm{Gyr}10\text{\,}\mathrm{Gyr}$ (middle panel of Fig. 2). Assuming a measurement uncertainty of $0.75\text{\,}\mathrm{Gyr}$ for each GC age continues to show a clear bimodality in one case (orange distribution), but only leaves a hint for the other case (red distribution). This shows that even if there is a bimodal age distribution, the large measurement uncertainties for GC ages can make it difficult to actually confirm it in practice.

Figure 2: GC ages and virial mass evolutions of two selected modeled MW-mass galaxies compared to the observed GC ages in the MW. Top: Cumulative GC age distributions of the MW from Usher et al. (2019) (blue) and of two modeled GC populations in simulated MW-mass galaxies from Valenzuela et al. (2021) (red and orange). Middle: Age distributions of the same three GC populations as in the top panel. The smooth distributions are smoothed using an assumed uncertainty of $0.75\text{\,}\mathrm{Gyr}$ for the GC ages. These are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations. For the sample from Usher et al. (2019), the distribution was scaled by a factor of four to be comparable to the modeled populations. Bottom: Dark matter virial mass evolution of the two modeled galaxies, showing their accretion histories. The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018).

The underlying reason for these peaks in the age distributions becomes apparent when studying the accretion histories of the two galaxies: both of them experience a major merger in the same time period as the GC formation bursts (bottom panel of Fig. 2). The mergers occur at the same time as the GSE merger is estimated to have happened (8–$11\text{\,}\mathrm{Gyr}$ ago). The two accretion histories differ strongly in their later evolution, however: while the galaxy with the more pronounced bimodal GC age distribution only experiences mini to minor mergers afterwards (orange line in the bottom panel of Fig. 2), the other galaxy has a second major merger at a later time of around 2–$4\text{\,}\mathrm{Gyr}$ ago (red line). The lack of GCs having formed around that time is a clear indication that the merger was rather gas-poor (i.e., dry). The formation history of the first galaxy is therefore more similar to that inferred for the MW (?MW-analog?), while the other galaxy has had a much more violent recent history (?non-MW-analog?). The fact that the non-MW-analog GC age distribution seen in the middle panel of Fig. 2 is more similar to the observed one by Usher et al. (2019) only indicates that the major merger around $10\text{\,}\mathrm{Gyr}$ ago is more similar to the GSE merger in terms of their GCs formed. However, the second major merger around 2–$4\text{\,}\mathrm{Gyr}$ ago is in no way comparable to the MW, making it the non-MW-analog. As shown for the other 19 MW-mass galaxies in Appendix B, the peaks in their GC age distributions also correspond to gas-rich merger events, while those mergers that do not leave a strong imprint in the GC age distribution are the dry mergers that the galaxy experienced.

In conclusion, the model makes the following predictions: wet mergers with a large amount of cold gas are capable of producing a bimodal or even multimodal GC age distribution for a galaxy, providing an indication for a type of event that is generally difficult to trace through other means. However, dry mergers with little to no cold gas do not result in noticeable signatures in the GC age distributions. This is the case for the majority of late-time mergers, but at higher redshifts gas-rich mergers are increasingly more common and beyond $z\approx 2$ the most common kind of merger event (e.g., Bournaud et al. 2011). Thus, GC ages provide a means to probe the very early turbulent formation times of galaxies by tracing their massive wet merger events.

4 Discussion: GCs in the Milky Way

Having the model prediction that wet mergers leave an imprint on the GC age distribution of a galaxy, we tested it through the available GC measurements in the MW. For this we could make use of additional properties currently unavailable for GCs around other galaxies, such as kinematic phase-space information and more accurate age measurements.

4.1 Milky Way diagnostics

To address some of the caveats brought up in Sect. 3.1, we used further data available on the MW GCs to study to what extent the model prediction can be applied to the MW and its GC population. Using the GC progenitor assignments by Forbes (2020), the age distributions for four of the samples (Forbes & Bridges 2010; Dotter et al. 2010, 2011; VandenBerg et al. 2013; Cabrera-Ziri & Conroy 2022) can be split up by those assignments. Figure 3 shows the age distributions for the GCs formed in situ in the MW and those associated with GSE and the other GC host progenitors, for each of the four samples. This figure uses the age uncertainties to smooth the distributions. For the total GC age distributions see Fig. 7, which shows that for the smoothed distributions, only the total age distribution from Forbes & Bridges (2010) shows a hint at a bimodality. This bimodality is smoothed out when taking the mean GC ages from Kruijssen et al. (2019), which may lead to GC age biases (discussed in Sect. A.1). See Sect. A.4 for the unsmoothed age distributions plotted as histograms for the most relevant GC host progenitor groups that correspond to Fig. 3.

Figure 3: Globular cluster age distributions in the MW from Forbes & Bridges (2010), Dotter et al. (2010, 2011), VandenBerg et al. (2013), and Cabrera-Ziri & Conroy (2022), split up according to their likely progenitors from Forbes (2020). The shown classifications are the MW, GSE, unassociated high-energy GCs (H-E), unassociated low-energy GCs (L-E, which was given the name Koala by Forbes 2020), the Helmi Streams (H99), Sequoia (Seq), and Sagittarius (Sag). The total number of GCs in the respective sample is indicated in the top left. The distributions are computed from the GC ages and their uncertainties through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations. The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018). See Fig. 12 for the corresponding histograms of the most relevant GC progenitor groups without taking the uncertainties into account.

The two smaller samples from VandenBerg et al. (2013) and Cabrera-Ziri & Conroy (2022) do not reveal further information besides an indication for the low-energy GCs to have formed at slightly later times for the measurements by VandenBerg et al. (2013). In contrast, both of the other samples show that there is a physical reason for the bimodal GC age distribution shown in Fig. 1: in the largest sample by Forbes & Bridges (2010), it is clearly seen that the younger peak seen in Fig. 1 at $11\text{\,}\mathrm{Gyr}$ is dominated by GSE-associated GCs, with contributions from the Helmi Streams and the high-energy GCs. This is less pronounced in the sample by Dotter et al. (2010, 2011), where the high-energy GCs contribute most of the young GCs, albeit there is also a significant contribution from GSE. The fact that these GCs have a different origin than the other MW GCs is known through the age-metallicity relation (e.g., Leaman et al. 2013; Kruijssen et al. 2019; Forbes 2020), in which the GCs associated with GSE and the other progenitor galaxies follow their own tracks. While the Sagittarius dwarf is expected to have fallen into the MW around 8–$9\text{\,}\mathrm{Gyr}$ ago on a most likely rather circular orbit, its extended time range of GC formation could be a result of tidally induced star cluster formation during its orbit around the MW while there was still enough cold gas available (e.g., Williams et al. 2022, for a study on young clusters with ages of around $2\text{\,}\mathrm{Gyr}$ in the Small Magellanic Cloud, which are suggested to have formed through tidal interactions with the Large Magellanic Cloud).

It is curious, however, that there appears to be a hint at a non-unimodal age distribution within the GSE GCs, which can be seen more clearly for Forbes & Bridges (2010) than for Dotter et al. (2010, 2011) in Fig. 3. The bimodality is also present when only using the GC ages from the subsample of Marín-Franch et al. (2009) (Sect. A.3) as well as when using the mean GC ages from Kruijssen et al. (2019). The bimodality is even slightly more prominent when using those ages, despite the overall distribution being more smoothed out (Sect. A.1). Additionally, it is also present when using the GC progenitor associations from any of Massari et al. (2019), Callingham et al. (2022), or Chen & Gnedin (2024b), and it is actually much more clearly visible for their classifications (Sect. A.2). The GCs associated with GSE are therefore shown to have an age bimodality across multiple different studies of their ages and association models with the GSE. The reason the bimodality is not seen in the samples from Dotter et al. (2010, 2011) and VandenBerg et al. (2013) is that they each only contain 13 GCs associated with GSE, compared to 21 GCs from Forbes & Bridges (2010). Considering the raw age data without uncertainties, both distributions are actually bimodal (Fig. 12), with the caveat of there being a very small amount of GCs involved.

Still, if the signal is real, it can be brought together with the model prediction in a straightforward manner. As the model forms GCs in small halos at early times (first pathway as described in Sect. 2.1), a GSE satellite galaxy hosts its own GCs as it falls into the MW. We refer to these GCs as the accreted GCs in the following. Assuming there is enough cold gas available, the major merger event then triggers GC formation through the gas collision and tidal forces between the two galaxies (merger-induced GCs in the following, second pathway as described in Sect. 2.1; Ashman & Zepf 1992; Williams et al. 2022). This scenario leads to multiple properties arising for the GCs: (1) The age distribution of the combined early-formed accreted and merger-induced GCs is bimodal. (2) Assuming that GSE brings in its own cold gas, many of the merger-induced GCs are also formed from that gas, or a mixture of that and the MW’s gas and therefore have similar metallicities and phase-space properties as those of the accreted GCs. This would then also lead to such GCs being associated with GSE through an analysis of phase-space and the age-metallicity relation. (3) Assuming the merger leads to violent gas interactions, it is possible that some recently formed merger-induced GCs are ejected from the overall orbit of GSE, leading to unassociated high-energy GCs. In that case, the high-energy GCs could also be associated with GSE and thus the bimodality in the GSE GC age distribution would be even more enhanced. It is also possible for GCs to distribute themselves further apart in phase space through other dynamical processes as shown by Pagnini et al. (2023), potentially ending up in the high-energy regime. This may occur in a similar fashion to the Splash (e.g., Bonaca et al. 2017; Belokurov et al. 2020), a group of more metal-rich stars in the MW halo that appear to have been formed in situ and dynamically ejected around the time of the GSE merger (e.g., Belokurov et al. 2020; Ciucă et al. 2023). Of course, this would not result in a change in age and metallicity of the formed GCs, though Ciucă et al. (2023) argue that the GSE merger could have first driven down the metallicity, which was then again enriched by the induced starburst.

While it is not possible to prove this theory with the current state of GC age precision and the difficulties of kinematic associations, the GC formation model indicates that the GCs associated with GSE are not only those that were brought in by the accreted galaxy (e.g., Côté et al. 1998), but also those that were formed through the wet merger (e.g., Ashman & Zepf 1992). Additionally, it is possible that some of the unassociated high-energy GCs were also formed in the process of the GSE merger, though it should be noted that many of those GCs have lower measured metallicities than those associated with GSE (Forbes 2020). In turn, this could be a result of the metallicity lowering through the merger (Ciucă et al. 2023). However, analyzing and modeling the details of GC metallicities in such a gas-rich major merger scenario are beyond the scope of this work and will be addressed in a future study.

4.2 Model diagnostics

From Eq. 3 it is possible to obtain the number of GCs expected from the model to form through a merger event based on the cold gas mass, $M_{\mathrm{gas}}$. Vincenzo et al. (2019) estimated GSE to have brought in a cold gas mass of $6.62\text{\times}{10}^{9}\text{\,}\mathrm{M_{\odot}}$. For the MW, we estimated a range of possible cold gas masses by determining the cold gas masses of galaxies at $z=2$ in a hydrodynamical cosmological simulation since the simulation that our GC model was applied to only contains DM particles. For this, we used Magneticum Pathfinder Box4 (uhr) (with a side length of $68\text{\,}\mathrm{Mpc}$ and a gas particle mass of $m_{\mathrm{gas}}=$1.0\text{\times}{10}^{7}\text{\,}\mathrm{M_{\odot}}$$) because it has the necessary resolution to resolve MW-like galaxies and their formation well, and at the same time it is large enough to contain a large sample of MW-like galaxies, from which we can draw conclusions about the expected gas mass at $z=2$. The galaxies contained within it have been shown to agree well with observations across a broad range of properties (see Teklu et al. 2015 for details on the implementations and Valenzuela & Remus 2024 for an overview of comparisons to observations). For the galaxies at $z=2$ with virial masses $$5\text{\times}{10}^{11}\text{\,}\mathrm{M_{\odot}}$\leq M_{\mathrm{vir}}\leq$8\text{\times}{10}^{11}\text{\,}\mathrm{M_{\odot}}$$, the typical cold gas mass (which we define as star-forming gas particles with temperatures below ${10}^{5}\text{\,}\mathrm{K}$) is $M_{\mathrm{gas}}=$(2.1\pm 0.6)\text{\times}{10}^{10}\text{\,}\mathrm{M_{\odot}}$$. We computed the expected numbers of formed GCs for a range of different MW cold gas masses between $1.0\text{\times}{10}^{10}\text{\,}\mathrm{M_{\odot}}7.5\text{\times}{10}^{10}\text{\,}\mathrm{M_{\odot}}$, leading to 6 to 22 GCs being formed (Table 1).

Table 1: Predicted number of GCs, $N_{\mathrm{GC,formed}}$, formed through the GSE merger with the MW for different assumed MW total cold gas masses, $M_{\mathrm{gas,MW}}$. The GSE galaxy is assumed to have had a cold gas mass of $6.62\text{\times}{10}^{9}\text{\,}\mathrm{M_{\odot}}$ (Vincenzo et al. 2019).

$M_{\mathrm{gas,MW}}/$\mathrm{M_{\odot}}$$	$N_{\mathrm{GC,formed}}$
$1.0\text{\times}{10}^{10}$	$6$
$2.5\text{\times}{10}^{10}$	$11$
$5.0\text{\times}{10}^{10}$	$17$
$7.5\text{\times}{10}^{10}$	$22$

As seen in the GC age distribution from Forbes & Bridges (2010), there are 12 out of 21 GCs with available ages associated with GSE and 7 out of 9 unassociated high-energy GCs in the time range of GSE and the secondary peak of the overall MW GC age distribution. While the latter GCs are surely not all related to the GSE merger, the GC sample also does not include all GCs found in the MW (around 50–$60\text{\,}\mathrm{\char 37\relax}$). Overall, the number of GCs expected from the model to be formed through such a merger aligns well with the observed number of GCs found to be associated with GSE or to be unassociated with high energies, while also taking into account the statistical incompleteness of the GC sample. This supports the prediction that multimodal GC age distributions can be used to trace wet mergers.

5 Discussion: Extension to nearby galaxies

While the GC age measurements are by far the most accurate for the MW due to resolved measurements of the clusters, the general properties from integrated GC age measurements can still give indications about the wet merger history of the host galaxy, through the relative ages between the GCs. One of the galaxies studied by Usher et al. (2019) is NGC 3115, a fast-rotating S0 field galaxy with stellar mass $M_{*}=$9\text{\times}{10}^{10}\text{\,}\mathrm{M_{\odot}}$$ and virial mass $M_{\mathrm{vir}}=$1.2\text{\times}{10}^{12}\text{\,}\mathrm{M_{\odot}}$$ (Forbes et al. 2016, assuming an NFW profile such that the virial mass is a factor 10 larger than the DM mass within $8\,R_{e}$, where $R_{e}$ is the effective radius, the radius within which half the light of the galaxy is emitted). It has $550\pm 80$ GCs (Harris et al. 2013) and features multimodal age distributions (based on 116 measured GC ages from Usher et al. 2019; top panel of Fig. 4). The multimodal behavior of the ages is even retained when smoothing the histogram with a Gaussian kernel of $1.5\text{\,}\mathrm{Gyr}$ (smooth line), which is a value we selected to illustrate the effect of smoothing that can be caused by measurement uncertainties. There is in fact a galaxy in the simulation that has a very similar GC age distribution as NGC 3115, which can here be seen in the middle panel of Fig. 4 (the similarity of the distributions is seen especially well in the cumulative distributions shown in fig. 17 of Valenzuela et al. 2021). The simulated galaxy has a virial mass of $M_{\mathrm{vir}}=$1.8\text{\times}{10}^{12}\text{\,}\mathrm{M_{\odot}}$$, similar to that of NGC 3115, and a total of 525 GCs, consistent with NGC 3115, which places it $0.2\text{\,}\mathrm{dex}$ above the mean linear scaling relation from Burkert & Forbes (2020), but still within the observed scatter. Interestingly, both the observed and simulated ages feature three minor peaks in the distribution without smoothing over it, although this could very well be a coincidence given the large uncertainties for the observational measurements. For the simulation this means that there were overall three especially gas-rich mergers leading to an increased amount of GC formation. This occurred during a time of steady assembly between $6\text{\,}\mathrm{Gyr}12\text{\,}\mathrm{Gyr}$ ago (bottom panel of Fig. 4).

Figure 4: GC ages and virial mass evolution of the modeled NGC 3115-analog galaxy compared to the observed GC ages in NGC 3115. Top and middle: Age distributions of the GC population in NGC 3115 from Usher et al. (2019), with a sample of 116 GCs, and of the modeled GC population in an NGC 3115-analog galaxy from Valenzuela et al. (2021) with 525 GCs. The smooth distributions are smoothed using an assumed uncertainty of $1.5\text{\,}\mathrm{Gyr}$ for the GC ages. These are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations. Bottom: Dark matter virial mass evolution of the simulated NGC 3115-analog galaxy, showing its accretion history.

This agrees well with the inferred formation history of NGC 3115 obtained through kinematics from IFU data and tracer populations, the metallicity profiles, its mass distribution, and the study of its morphological components: it is believed to have experienced an early gas-rich accreting phase followed by a lack of significant mergers thereafter (Arnold et al. 2011; Guérou et al. 2016; Zanatta et al. 2018; Poci et al. 2019; Dolfi et al. 2020; Buzzo et al. 2021). In particular, this supports the prediction of GC ages being able to trace wet mergers, also for galaxies outside the Local Group.

Such behavior is not the norm, however. This could already be seen from the median age distributions presented by Valenzuela et al. (2021) for all the virial mass bins, which show that generally GC populations are dominated by old GCs like those of the MW. For most galaxies, GC formation bursts are too close in time to the oldest GC populations to be able to distinguish them without having further properties available like in the MW, or the bursts are not significant enough due to a lack of cold gas as gas-rich merger events become less and less frequent with time.

One such case is M 31, for which the GC age measurements of Wang et al. (2021), Usher et al. (2024), and Cabrera-Ziri et al. (in prep.) show no clear bimodal distribution (Fig. 5). There is a hint at some additional modes between $6\text{\,}\mathrm{Gyr}$ and $10\text{\,}\mathrm{Gyr}$ ago in the data from Usher et al. (2024) and Cabrera-Ziri et al. (in prep.) (middle and lower panel). The numbers of GCs in these modes are small for the latter study (6 around 6–$8\text{\,}\mathrm{Gyr}$ ago and 8 around 8–$10\text{\,}\mathrm{Gyr}$ ago) and the uncertainties of the GC age measurements are large compared to the individual age peaks (smooth curves in the middle and lower panel), so we cannot exclude that they originate from statistical uncertainties. If they do not, this would indicate only small gas-rich mergers given the small numbers of GCs produced. As discussed in Sect. 2.2.2, the absolute ages of the measured distributions should not be compared with those of Wang et al. (2021) due to differences in the measurement techniques, which led to the much younger determined ages for Wang et al. (2021). In terms of their relative ages, however, all of the samples show that significant features in the GC age distributions are not found for all galaxies and strong features are only visible for mergers that have a sufficient amount of gas. Additionally, which features are identifiable in GC age distributions will always depend on how large the measurement errors are, which smooth out the data, which first removes the signatures from smaller mergers or those with smaller gas fractions. Thus, the method is the most reliable in detecting wet mergers with large mass fractions. Finally, note that the sample from Cabrera-Ziri et al. (in prep.) is biased towards metal-rich GCs (Sect. 2.2.2). However, gas-rich major mergers are expected to involve the more metal-rich GCs as opposed to metal-poor ones, such that we believe the implications of our analysis to be unaffected.

Figure 5: Globular cluster age distribution in M 31 from Cabrera-Ziri et al. (in prep.), Usher et al. (2024), and Wang et al. (2021), with sample sizes of 136 GCs, 290 GCs, and 293 GCs, respectively. The ages in the sample of Wang et al. (2021) correspond to their determined old clusters ($t_{\mathrm{age}}>$1.5\text{\,}\mathrm{Gyr}$$). For further comparisons with previous GC age determinations in M 31, see fig. 8 of Wang et al. (2021). The smooth lines show the distributions smoothed by the measurement uncertainties. These are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations for Usher et al. (2024) and over log-normal distributions for the other two data sets. Note that for the logarithmic uncertainties, the visualization in linear space is skewed with respect to the peak.

Observationally, it has been proposed that M 31 experienced a major merger at around $2\text{\,}\mathrm{Gyr}$ ago, possibly with the progenitor of M 32 (D’Souza & Bell 2018a, b; Hammer et al. 2018). A large peak in star formation between $2\text{\,}\mathrm{Gyr}4\text{\,}\mathrm{Gyr}$ ago in the disk and outskirts of M 31 (Bernard et al. 2012, 2015; Williams et al. 2017) is likely related to this event. Using planetary nebulae, Bhattacharya et al. (2023) find further evidence for a wet major merger of M 31 $2.5\text{\,}\mathrm{Gyr}$ to $4\text{\,}\mathrm{Gyr}$ ago. In fact, the full sample by Wang et al. (2021) also includes young stellar clusters with ages below $1.5\text{\,}\mathrm{Gyr}$, which could potentially also have been formed as a result of the gas brought in by the merger and would not be referred to as GCs yet. Similarly, the distribution by Usher et al. (2024) shows GCs with these kind of young ages as well as GCs with ages around 3–$4\text{\,}\mathrm{Gyr}$, and Cabrera-Ziri et al. (in prep.) also find one GC with an age of $2.5\text{\,}\mathrm{Gyr}$, which could coincide with that merger.

It has also been determined that the star formation rate was very low before the recent peak, with most stars having formed prior to $8\text{\,}\mathrm{Gyr}$ ago (Williams et al. 2017). Due to the smooth and very massive stellar halo observed for M 31, it is estimated that the merger history was dominated by many smaller accretion events (e.g., Ibata et al. 2014; Mackey et al. 2019). This scenario is supported by the model prediction presented in this work, in which no single sufficiently massive merger exists that would lead to a significant GC formation burst.

Finally, the massive elliptical galaxy NGC 1407 also lacks a significant second peak in its GC age distribution, though there is a tail with a slight peak towards younger ages around 6–$9\text{\,}\mathrm{Gyr}$ ago and potentially another around $10\text{\,}\mathrm{Gyr}$ ago (Fig. 6). Due to the integrated age measurements, this could be the result of underestimated ages for those GCs as the peak again consists of only very few GCs. Overall, it appears that from the GC ages NGC 1407 has not experienced any massive wet mergers since the early buildup phase of the galaxy, and at most a later merger with a low cold gas fraction. In the latter case, it should be expected that there would also be a sign of late star formation activity in the stellar populations themselves. In fact, Spolaor et al. (2008) found from their stellar population measurements of NGC 1407 that the stars are uniformly old, having formed around $12\text{\,}\mathrm{Gyr}$ ago.

Figure 6: Globular cluster age distribution in NGC 1407 from Usher et al. (2019) with a sample size of 213 GCs. The smooth line shows the distribution smoothed using an assumed uncertainty of $1.5\text{\,}\mathrm{Gyr}$ for the GC ages. These are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations.

A kinematically decoupled core (KDC) hints at a major merger of NGC 1407 with gas fractions between $15\text{\,}\mathrm{\char 37\relax}$ and $40\text{\,}\mathrm{\char 37\relax}$ (Hoffman et al. 2010; Schulze et al. 2017). While Forbes & Remus (2018) find that a simulated galaxy from Magneticum with a similar size and mass as NGC 1407 had a late major merger around $8\text{\,}\mathrm{Gyr}$ ago, it was not selected based on further properties such as stellar ages or metallicity gradients. However, Ferré-Mateu et al. (2019) find in their observations of NGC 1407 that the KDC is slightly younger than the rest of the galaxy (we estimate a difference of around $1\text{\,}\mathrm{Gyr}$ based on their fig. 5), suggesting that a wet major merger could have occurred slightly later than the early buildup of the galaxy. This scenario is compatible with the GC age distribution found by Usher et al. (2019): the lack of late gas-rich mergers is consistent with no large amount of late GC formation, and the slight peak in the GC ages could be related to the wet major merger that formed the KDC. Such a merger is expected to have had a relatively high gas fraction (Hoffman et al. 2010), potentially forming many GCs as a result. However, due to the violent nature of such a merger, many of those systems would likely be disrupted in the same process, leading to a smaller peak in the age distribution. Since the GC sample analyzed by Usher et al. (2019) is biased towards central GCs due to the SLUGGS survey having been focused on the inner regions, it is expected that it would be more likely to pick up signatures from major mergers.

6 Conclusion

In this work, we have used the GC formation model from Valenzuela et al. (2021) with dual formation pathways to predict that massive wet mergers with enough cold gas can leave an imprint on the age distribution of the GC population in the host galaxy. This imprint results in a bimodal or even multimodal distribution, indicating when the wet mergers occurred. This prediction is in part also a consequence of the idea that red GCs tend to form through mergers that also induce star formation, thus resulting in properties that overall trace the underlying stellar component, such as spatial, kinematic, or chemical properties (e.g., Brodie & Strader 2006; Pota et al. 2013; Dolfi et al. 2021). In contrast, mergers with little to no gas are not traced by the GC ages since the lack of gas means that no significant number of GCs could be formed in the process.

The prediction is discussed for the MW in detail. We find that a hint at bimodality visible in the pure GC age distributions compiled by Forbes & Bridges (2010) can be further disentangled by combining the data with phase-space information on the GCs to map which galaxy progenitors the GCs are likely associated with (e.g., Massari et al. 2019; Forbes 2020; Callingham et al. 2022; Chen & Gnedin 2024b). We find the later peak in the GC age distribution to correspond to the GCs associated with GSE, the last major merger of the MW (Belokurov et al. 2018; Helmi et al. 2018), and in part also to the unassociated high-energy GCs and the GCs of the Helmi Streams. In fact, the age distribution of the GSE GCs appears to also have a bimodality. Since GSE is expected to have been a massive and gas-rich merger, we suggest that GSE not only brought in its own older GCs, but also formed a second group of GCs through the merger with the MW. The second group of GCs would be located near the first in phase space and lie on the same age-metallicity relation, as is also found for them in the observations. We further suggest that some of the unassociated high-energy GCs may also originate from the GSE merger, since it is possible that the violent merger dynamics would eject some GCs, or that GCs migrate away in phase space over time (Bonaca et al. 2017; Belokurov et al. 2020; Pagnini et al. 2023).

Two simulated MW-mass galaxies with similar GC age distributions as the MW are found to both have had a wet major merger around the same time as GSE (around $8\text{\,}\mathrm{Gyr}$ to $10\text{\,}\mathrm{Gyr}$ ago). One of the two then evolved only through smooth accretion and mini and minor mergers, as is believed to have been the case for the MW. In contrast, the other simulated galaxy encountered a dry merger at later times that is not traced by the GC ages, (for those kind of mergers, other kind of indicators have to be used).

We also tested the model with three other observed galaxies, NGC 3115, M 31, and NGC 1407, for which GC ages have been obtained. For NGC 3115, the observed GC age distribution is clearly multimodal with a considerable population of younger GCs that are 7–$11\text{\,}\mathrm{Gyr}$ old. We were able to find a simulated galaxy of comparable virial mass and a similar GC age distribution. It underwent an early phase of multiple wet mini and minor mergers that led to the formation of the younger GC population, after which it experienced no further significant mergers. This agrees well with the expected formation history of NGC 3115 inferred from other tracers (Arnold et al. 2011; Guérou et al. 2016; Poci et al. 2019; Buzzo et al. 2021) and supports the prediction that additional modes in the GC age distribution trace wet mergers.

The GC age distribution of M 31 features no bimodality, indicating that it experienced no significant late wet mergers. However, it shows two or even three small peaks in the GC age distributions at 2–$3\text{\,}\mathrm{Gyr}$ and at 7 and $9\text{\,}\mathrm{Gyr}$, indicating either small merger events with large gas fractions or large merger events with low gas fractions, with the former being more likely than the latter. Since these peaks are small and the uncertainties are large, this is not conclusive, however. Similarly, NGC 1407 features no clear bimodality, but a possible second late peak in the GC age distribution may be related to an early wet major merger that led to the kinematically distinct core found in the overall very old galaxy. These are examples for the most common case for galaxies: the simulation predicts that strong peaks in the GC age distribution with GCs younger than around $12\text{\,}\mathrm{Gyr}$ clearly indicate a massive wet merger event and that small wet mergers can still leave minor peaks in the GC age distribution. Dry mergers, however, leave no imprints in the GC age distributions.

To conclude, the age distribution of GCs can be used as a tracer for wet mergers of galaxies, which are generally more difficult to infer from observations: while the old age peak of the GC age distribution does not help in constraining the merger history as here the old populations from the main progenitor of a galaxy mix with the old GCs brought in through other merging galaxies of all kind, the young GCs with ages less than around $11\text{\,}\mathrm{Gyr}$ are formed in the otherwise untraceable wet merger events. However, due to the current large observational uncertainties in determining GC ages from integrated measurements, further development of the age determination techniques will be essential to better understand individual galaxies’ formation histories. Finally, increasing the number of galaxies with accurate GC age measurements will also help improve our understanding of GC formation and set more constraints on current GC formation models. We thus propose that this is a good method to infer the wet merger times of extragalactic galaxies from observations. With the advent of observatories such as Euclid, the Nancy Grace Roman Space Telescope, and the Vera C. Rubin Observatory, applying photometric methods to obtain GC ages from such observations has the great potential of providing deeper insight into the early gas-rich formation history of galaxies.

Appendix A Milky Way globular cluster age distributions

A.1 Age distributions of all globular clusters

The age distributions of the full GC samples of Forbes & Bridges (2010), Dotter et al. (2010, 2011), VandenBerg et al. (2013), Kruijssen et al. (2019), Usher et al. (2019), and Cabrera-Ziri & Conroy (2022) are shown in Fig. 7. The values from Kruijssen et al. (2019) are obtained through the mean ages from the first three samples, all of which used CMD measurements, which leads to the total distribution being further smoothed out. Since not all of the CMD samples contain all 96 GCs from that study, the mean ages of the individual GCs will be biased towards younger or older ages depending on which samples they are contained in due to differences in the systematic uncertainties between the samples.

Figure 7: Globular cluster age distributions in the MW from Forbes & Bridges (2010), Dotter et al. (2010, 2011), VandenBerg et al. (2013), Kruijssen et al. (2019), Usher et al. (2019), and Cabrera-Ziri & Conroy (2022), smoothed by the measurement uncertainties. The curves are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations. The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018).

Nevertheless, the bimodality in the GSE GC age distribution is still clearly visible (Fig. 8), if not even clearer than in the top left panel of Fig. 3 for the ages from Forbes & Bridges (2010). This strongly indicates that the found bimodality is an actual feature and not a statistical artifact.

Figure 8: Globular cluster age distributions in the MW from Kruijssen et al. (2019), split up according to their likely progenitors from Forbes (2020). The shown classifications are the MW, GSE, unassociated high-energy GCs (H-E), unassociated low-energy GCs (L-E, which was given the name Koala by Forbes 2020), the Helmi Streams (H99), Sequoia (Seq), and Sagittarius (Sag). The curves are smoothed by the measurement uncertainties, which are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations. The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018).

A.2 Alternative progenitor associations

While we used the GC progenitor associations from Forbes (2020) in this work, it makes no qualitative difference to use the associations determined by other studies. For instance, the data from Callingham et al. (2022) show a somewhat different distribution of GCs among GSE, the Helmi Streams, and the low-energy GCs, but the GC age distributions of GSE still features a bimodal distribution (Fig. 9). This bimodality is visible in a much more distinct way when compared with the top left panel of Fig. 3, which shows the distributions for the GC associations determined by Forbes (2020). Similarly, when using the associations from Chen & Gnedin (2024b), Fig. 10 shows that the age distributions change for GSE and Sagittarius (the only two ex-situ groups they specifically identify), but again the bimodality of the GSE-associated GCs persists. This time, the younger peak at $11\text{\,}\mathrm{Gyr}$ is significantly stronger than the ones seen from the associations of Forbes (2020) and Callingham et al. (2022), which would mean that potentially more GCs were formed through the GSE merger.

Figure 9: Globular cluster age distributions in the MW from Forbes & Bridges (2010), split up according to their likely progenitors from Callingham et al. (2022). The shown classifications are the MW, GSE, unassociated high-energy GCs (H-E), unassociated low-energy GCs (L-E, which was given the name Koala by Forbes 2020), the Helmi Streams (H99), Sequoia (Seq), and Sagittarius (Sag). The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018).

Figure 10: Globular cluster age distributions in the MW from Forbes & Bridges (2010), split up according to their likely progenitors from Chen & Gnedin (2024b). The shown classifications are the MW, GSE, Sagittarius (Sag), and other ex-situ formed GCs (other). The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018).

A.3 Age distributions of a single study

The GC sample of Forbes & Bridges (2010) is composed of multiple different studies (see Sect. 2.2.1), of which the largest subsample is of Marín-Franch et al. (2009), which includes 64 of the 92 GCs. As combining several GC studies could cause problems due to different systematic biases, we show the single GC study by Marín-Franch et al. (2009) in Fig. 11. Both the bimodality of the overall GC age distribution as well as in the GSE GC age distribution are clearly visible. This demonstrates that these bimodalities do not result from combining several GC studies, but in fact come from a single study.

Figure 11: Globular cluster age distributions in the MW from Marín-Franch et al. (2009), split up according to their likely progenitors from Forbes (2020). The shown classifications are the MW, GSE, unassociated high-energy GCs (H-E), unassociated low-energy GCs (L-E, which was given the name Koala by Forbes 2020), the Helmi Streams (H99), Sequoia (Seq), and Sagittarius (Sag). The curves are smoothed by the measurement uncertainties, which are computed through a summation over normal distributions with the respective ages as the means and their uncertainties as the standard deviations. The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018).

Figure 12: Globular cluster age distributions in the MW from Forbes & Bridges (2010), Dotter et al. (2010, 2011), VandenBerg et al. (2013), and Cabrera-Ziri & Conroy (2022), split up according to their likely progenitors from Forbes (2020). The shown classifications are the total population (all), GSE, the unassociated high-energy GCs (H-E), and the Helmi Streams (H99). The total number of GCs in the respective sample is indicated in the top left. The shaded region between $8\text{\,}\mathrm{Gyr}11\text{\,}\mathrm{Gyr}$ indicates the estimated time of the GSE merger (Belokurov et al. 2018; Helmi et al. 2018). See Fig. 3 for the corresponding age distributions smoothing the data points through their uncertainties.

A.4 Age histograms

The unsmoothed age distribution of the MW are shown in Fig. 12 for the four GC samples by Forbes & Bridges (2010), Dotter et al. (2010, 2011), VandenBerg et al. (2013), and Cabrera-Ziri & Conroy (2022), using the progenitor associations given by Forbes (2020). It becomes apparent that GSE, the unassociated high-energy GCs, and the GCs from the Helmi Streams are important contributors to the second peak in the total GC age distribution. The smoothed distributions according to the measured age uncertainties are found in Fig. 3.

Appendix B Other Milky Way-mass modeled globular cluster age distributions

As we found 21 MW-mass galaxies in the simulation with virial masses of $M_{\mathrm{vir}}=1$–$2\text{\times}{10}^{12}\text{\,}\mathrm{\mathnormal{M}_{\mathrm{vir}}}$, but only discussed the two of them with GC age distributions most similar to those observed in the MW (see Sect. 3.2), we here present the other 19 simulated galaxies. Figure 13 shows their GC age distributions (top of each of the four double panels) and their respective virial mass evolutions in the same color (bottom of the double panels). From the top left to the bottom right, the galaxies are sorted by their formation time, that is at what time the galaxy reached half of the virial mass that it has at $z=0$, where the galaxies formed the earliest are shown at the top left. It can immediately be seen that the galaxies formed the earliest also have the oldest GC populations, whereas the later-formed galaxies have increasingly young GC subpopulations. The time of the youngest GC formation peak is furthermore strongly correlated with the formation time as seen in Fig. 14. From the MW-analog, we also obtain a prediction for the formation of the MW itself of around 9–$10\text{\,}\mathrm{Gyr}$, which corresponds to $z=1.3$–1.7. Overall, the ages of the youngest GCs also directly mark the time of the last major wet accretion event, as expected from the GC model.

Figure 13: Age distributions of the globular clusters in 19 simulated MW-mass galaxies and their virial mass evolutions plotted below the respective age distributions in the same colors for each galaxy. The distributions are smoothed using an assumed uncertainty of $0.5\text{\,}\mathrm{Gyr}$ for the GC ages. The galaxies shown in the top left panels are the galaxies that reached half of their $z=0$ virial mass at the earliest time, followed by the galaxies in the top right panel. The galaxies in the bottom right panel are those that reach half of the current virial mass the latest. The galaxy with the dashed virial mass evolution line has an issue with the merger tree, but the GC population is correct.

Figure 14: Correlation between the age of the youngest GC peak in their age distribution ($t_{\text{last GC peak}}$) and the formation time of the MW-mass simulated galaxies ($t_{\mathrm{form}}$), colored by the virial mass at $z=0$. The youngest peak is required to have an absolute value above 10. The formation time of a galaxy is taken as the lookback time at which the galaxy has reached half of its virial mass at $z=0$. The MW-analog and non-MW-analog from Fig. 2 are additionally circled in orange and red, respectively.

Overall, the oldest GC population is almost always the dominant one (with the exception of the double-peaked age distributions in the upper left panel shown in pink and in the upper right panel shown in blue). In these two cases, there are nearly as many of the oldest GCs as the second youngest peak in the age distribution contains.