Catalogue of model star clusters in the Milky Way and M31 galaxies

Detailed understanding of the formation and evolution of globular clusters (GCs) has been recently advanced through a combination of numerical simulations and analytical models. We employ a state-of-the-art model to create a comprehensive catalogue of simulated clusters in three Milky Way (MW) and three Andromeda (M31) analogue galaxies. Our catalogue aims to connect the chemical and kinematic properties of GCs to the assembly histories of their host galaxies. We apply the model to a selected sample of simulated galaxies that closely match the virial mass, circular velocity profile, and defining assembly events of the MW and M31. The resulting catalogue has been calibrated to successfully reproduce key characteristics of the observed GC systems, including total cluster mass, mass function, metallicity distribution, radial profile, and velocity dispersion. We find that clusters in M31 span a wider range of age and metallicity, relative to the MW, possibly due to M31's recent major merger. Such a merger also heated up the in-situ GC population to higher orbital energy and introduced a large number of ex-situ clusters at large radii. Understanding the impacts of galaxy mergers and accretion on the GC populations is crucial for uncovering the galaxy assembly histories.


INTRODUCTION
Over the past decade, significant observational efforts have been made to uncover the origins of globular clusters (GCs) in the local universe.Spectroscopic surveys like APOGEE (Majewski et al. 2017) have offered valuable insights into the chemical compositions of these ancient stellar clusters.Particularly, the advent of the Gaia mission (Gaia Collaboration et al. 2016Collaboration et al. , 2018Collaboration et al. , 2023) ) has revolutionized our understanding of the GC spatial and kinematic properties, and their stellar populations.Studies such as Massari et al. (2019) and Malhan et al. (2022) took advantage of the Gaia data to paint a comprehensive picture of the origins of Galactic GCs, shedding light on key questions related to the formation and evolutionary history of GCs, including where GCs formed and how they were brought to their current locations in the galaxy.
At the same time, numerous observational and theoretical studies have improved our understanding of the formation history of the Milky Way (MW) and Andromeda (M31) galaxies.For example, Deason et al. (2015) investigated the ratio of blue straggler stars to blue horizontal branch stars in the MW halo and suggested the accretion of massive satellite galaxies as progenitors of the stellar halo.Inspired by this, subsequent chemical and kinematic studies focused on disk and halo stars (Belokurov et al. 2018;Helmi et al. 2018;Deason et al. 2018) and confirmed the likely major merger at lookback time ≳ 10 Gyr.The progenitor galaxy of this merger is referred to as the Gaia-Sausage/Enceladus (GS/E).Analysis of the kinematics of MW stars (Belokurov & Kravtsov 2022) showed that after this early bursty star forming epoch, the Galaxy transitioned to a steady stage of disk formation within 1 − 2 Gyr.Cosmological simulations ★ E-mail: ybchen@umich.edu of galaxy formation also supported the early bursty stage and the subsequent transition to a steady state (Yu et al. 2023;Semenov et al. 2023).On the other hand, recent studies found that M31 galaxy had a distinct assembly history from the MW, characterized by a massive merger with an M32-like progenitor around 2 Gyr ago (D'Souza & Bell 2018).
Our understanding of the formation and evolution of star clusters throughout cosmic time has also been greatly advanced over the past decade.Some of the young massive clusters formed in high-redshift galaxies may survive tidal disruption, becoming GCs at the present.In this framework, theoretical modelling of GCs has become possible by embedding necessary recipes of cluster formation and evolution as sub-grid prescriptions into cosmological simulations.For example, Li et al. (2017Li et al. ( , 2018) ) and Li & Gnedin (2019) treated forming star clusters effectively as sink particles centered on giant molecular clouds in a suite of zoom-in hydrodynamical simulations.The E-MOSAICS project (Pfeffer et al. 2018;Kruĳssen et al. 2019) applied the MOSAICS model for cluster formation and evolution (Kruĳssen & Lamers 2008;Kruĳssen 2009;Kruĳssen et al. 2011) to a re-run of the EAGLE simulations (Schaye et al. 2015;Crain et al. 2015).The EMP-Pathfinder project (Reina-Campos et al. 2022) further updated the E-MOSAICS physics recipes with new prescriptions for cluster formation and the multi-phase interstellar medium (ISM).These works were able to match various observational properties of both young and old clusters.They also offered insights to the links between GC properties and the assembly history of their host galaxy (Kruĳssen et al. 2019;Pfeffer et al. 2020;Reina-Campos et al. 2020;Trujillo-Gomez et al. 2021).
In addition to these direct simulations, post-processing methods have also gained traction in the field.Such methods employ analytical prescriptions for GC formation and evolution, applying them to the merger trees and particle outputs of existing cosmological simulations.Unlike full galaxy formation simulations, these models are not sensitive to the specific baryonic prescriptions used in simulations, and do not require costly re-runs.Such robustness and flexibility make them particularly advantageous when applied to large samples of galaxies.Examples of these methods are Renaud et al. (2017); Creasey et al. (2019); Phipps et al. (2020); Halbesma et al. (2020); Valenzuela et al. (2021Valenzuela et al. ( , 2023)), including previous versions of our model (Muratov & Gnedin 2010;Li & Gnedin 2014;Choksi et al. 2018;Choksi & Gnedin 2019a,b).
Our model takes the halo merger tree as the input and triggers GC formation when the specific mass accretion rate exceeds a predefined threshold.Utilizing a sequence of scaling relations, the model analytically calculates the mass and metallicity of the newly-formed GCs.It also accounts for mass loss due to stellar evolution and tidal disruption.Starting in Chen & Gnedin (2022, hereafter CG22), we have improved the model by including the spatial and kinematic information for GCs using tracer star particles from the simulations.This enabled a comprehensive comparison between model predictions and the 9-dimensional (mass + age + metallicity + 3D positions + 3D velocities) characteristics of MW GCs.A subsequent extension by Chen & Gnedin (2023, hereafter CG23) incorporated dark matter (DM) particles to represent GCs in collisionless simulations.Our refined model successfully replicated key statistics of observed GC systems, including the distributions of GC mass, metallicity, distance from the galaxy center, velocity dispersion and anisotropy.By applying this model across a wide range of galaxy mass, we reproduced global scaling relations such as the effective radius-galaxy mass and the nearly linear GC mass-galaxy mass correlations.
In this work we use the model results to produce a catalogue of model GCs with the properties that match the observations of the GC systems in MW and M31.The primary use of this catalogue is to help the analysis of observations to reveal the relationship between GC features and the assembly histories of the MW and M31.To achieve this, we start by selecting a set of simulated galaxies that meet all major observational constraints, including the virial mass, circular velocity profile, and the defining assembly events, such as the GS/Elike merger for the MW and the M32-progenitor-like merger for M31.Next, we calibrate our model parameters by comparing observed GC attributes with those in our model systems.These attributes include total GC number, mass function, metallicity distribution, radial profile, and velocity dispersion.The resulting catalogue can enhance our understanding of how galaxy assembly events influence the current distribution of GCs in the property space, such as the age-metallicity plane and the integral of motion (IOM) space.The catalogue enables us to assess the accuracy of various classification algorithms in these spaces to identify the original progenitors for ex-situ GCs and, ultimately, reconstruct the host galaxy's assembly history.
The rest of the paper is organized as follows.We outline the background simulations and detail all updates made to the CG22 and CG23 models in Sec. 2. Next, we define the criteria used to select the best-matching MW and M31 analogues within our simulations in Sec. 3. In Sec. 4, we present the GC systems in the selected analogues and compare them with the observational data.We also explore the influence of the last major merger in M31 on the spatial distribution of GCs.The discussion in Sec. 5 focuses on the potential application of the model to uncover galaxy assembly histories and addresses known caveats.We conclude and summarize the study in Sec. 6.In addition, we present a new functional form for the time-dependent galactic stellar mass-metallicity relation (MZR), used to determine the metallicity of model GCs, in Appendix A.

MODEL SETUP
To create model catalogues of GC systems in the MW and M31, we apply our GC formation model to 6 carefully selected galaxies from the cosmological Illustris TNG50-1 simulation (Nelson et al. 2019;Pillepich et al. 2019;Nelson et al. 2021, hereafter TNG50) and from a suite of collisionless simulations of a Local Group (LG) environment (CG23).These simulated galaxies have assembly histories similar to the MW or M31 and reproduce most of the observable properties of their GC systems.In this section, we provide an overview of the simulations and the model setup.We describe the criteria used to select the best MW/M31 analogues in the next section.
TNG50 applies the Friends-of-Friends algorithm and the subfind code (Springel et al. 2001) to identify halos and subhalos.We use the term 'galaxy' to refer to subhalos throughout the paper.Additionally, TNG50 provides galaxy merger trees using the sublink code (Rodriguez-Gomez et al. 2015).
The second suite contains two collisionless zoom-in simulations, which we ran in CG23.We performed the simulations with the adaptive refinement tree (ART, Kravtsov et al. 1997) code on the modified initial conditions (ICs) from the ELVIS suite (Garrison-Kimmel et al. 2014).The ICs are Thelma & Louise and Romeo & Juliet, each producing a galaxy environment similar to the LG.We refer to this simulation suite as the 'LG simulations' hereafter.
We run the rockstar (Behroozi et al. 2013a) halo finder and the consistent tree code to construct halo catalogues and merger trees for the LG simulations, respectively.The two LG simulations each produces two main galaxies.We find the Louise, Romeo, and Juliet galaxies have a quiescent 'MW-like' mass assembly history after  ∼ 5 with no major merger with a mass ratio less than 4:1.This feature is similar to the formation history of the MW (Hammer et al. 2007).On the other hand, Thelma has more major mergers at later times.

Modeling cluster formation and evolution
Our GC formation and evolution model is based on CG22, where we modified the previous versions of the model (Muratov & Gnedin 2010;Li & Gnedin 2014;Choksi et al. 2018) to include positional and kinematic information by tagging simulation particles as 'GC tracer particles'.In this section, we briefly describe the setup of the CG22 model and further modifications we make in this work.
Our model consists of four steps: 1) cluster formation, 2) cluster sampling, 3) particle assignment, and 4) cluster evolution.The formation of GCs is triggered by rapid mass growth of the host galaxy (e.g., major mergers or intense mass accretion), quantified by the specific mass accretion rate  h / h exceeding a threshold value  3 , which is an adjustable model parameter.We then apply the empirical stellar mass-halo mass (SMHM) relation (Behroozi et al. 2013b) to compute the stellar mass  ★ from the halo merger history; the gas mass-stellar mass relation (Lilly et al. 2013;Genzel et al. 2015;Tacconi et al. 2018;Wang et al. 2022) to calculate gas mass  gas from stellar mass; and finally the linear gas mass-cluster mass relation (Kravtsov & Gnedin 2005),  tot = 1.8 × 10 −4  2  gas , to compute the total GC mass  tot , where  2 is another model parameter quantifying the cluster formation rate.We also use the time-dependent MZR to assign the host galaxy metallicity to its population of GCs forming at a given epoch.We make some modifications to the scaling relations employed in CG22 to match the updates in the theoretical and observational results and to improve the modeling of scatter evolution.We provide details of these modifications in Sec.2.2.1.
Next, we sample the masses of individual clusters from the Schechter (1976) initial cluster mass function (ICMF) derived from observations of young star clusters.We extend the range of cluster masses down to 10 4 M ⊙ instead of 10 5 M ⊙ in CG22.This allows us to capture some surviving low-mass GCs in the outskirts of the galaxies.We do not model clusters with initial mass below 10 4 M ⊙ for two primary reasons.First, the lowest-mass galaxies we consider have halo masses around 10 8 M ⊙ (resolved by ∼ 200 particles in TNG50), and typically produce clusters with masses close to 10 4 M ⊙ , see Fig. 5 in CG23.Second, our tidal disruption model is effective in disrupting clusters of this mass within a few Gyr.According to calculations in CG23, a 10 4 M ⊙ cluster located at a galactocentric radius = 3 kpc would likely survive less than 1 Gyr.Therefore, extending the mass function below this limit is meaningless.
After obtaining the list of newly formed clusters, we assign them to certain types of simulation particles depending on the simulation.For the hydrodynamic simulations like TNG50, we first select young (age < 10 Myr) stellar particles within an initial radius of 3 kpc from the galactic center.We have tested different initial radii and found that a larger radius leads to a final spatial distribution more extended than observations.However, reducing the initial radius does not significantly impact the final spatial distribution of the clusters.When there are not enough newly formed stellar particles, we also use older stellar particles formed between the adjacent snapshots.In the rare cases (∼10%) when there are still not enough stellar particles, we use DM particles near the galactic center as GC tracer particles.We only use the positions and velocities of these particles to passively track individual GCs but calculate all other properties analytically.This minimizes the model dependence on the baryon physics employed in the hydrodynamic simulation.
On the other hand, for the collisionless LG simulations we follow CG23 to select collisionless particles in local density peaks near the galaxy center, corresponding to surviving dense cores of satellite galaxies or other galactic structure with deep potential wells where massive clusters are more likely to form than elsewhere.We identify peaks within the scale radius of the best-fit Navarro-Frenk-White (NFW) halo profile.We also require the peak to be denser than any of the 16 closest grid cells and 30 times the mean density enclosed within the scale radius.
Finally, we compute the GC mass loss due to the stellar evolution and dynamical disruption based on the local tidal field along their orbits.Different from CG22, we update the prescription for tidal disruption following Gieles & Gnedin (2023), where the disruption rate is a multivariate power-law function of the cluster initial mass and current mass.In Sec.2.2.2, we describe the modifications made to the cluster evolution step in greater detail.
To choose the best model parameters, we calibrate the model to match the metallicity distribution, mass function, total GC mass-halo mass relation, and the spatial distribution.We describe the calibration in Sec.3.2.

Modifications to the scaling relations
Some of the scaling relations in the CG22 model are modified.First, we modify the gas mass-stellar mass relation by updating the upper bound of the gas mass constrained by the extragalactic ultraviolet background after reionization.Following the more recent CG23 model, we require the sum of the gas fraction  gas =  gas / h and the stellar fraction  * =  ★ / h to be less than the total accreted baryon fraction  in parametrized by Kravtsov & Manwadkar (2022): where (2 /3 −1)  ] −3/ is a soft step-function, and  ch is the characteristic mass scale at which  in = 0.5  b : where We adopt the reionization epoch at  rei = 6 and  = 15 as in Kravtsov & Manwadkar (2022).If  gas +  * >  in , we enforce  gas =  in −  * by setting Next, although we still use the Behroozi et al. (2013b) SMHM relation, we change the way we model and evolve the scatter.In our previous models since Choksi et al. (2018), the scatter was modeled as the cumulative scatter of short-term star formation periods.Such a technique guarantees the stellar mass  ★ to be a monotonic function of time but may underestimate the resultant scatter at present-day.Also, this method computes the final  ★ as the summation of a series of log-normal distributions at each epoch, which may not result in a final log-normal distribution to match the present-day observations.Instead, in this work we model the evolution of scatter by a Gaussian process, log  ★ ( h , ) ∼ GP {log SMHM( h , ),  ( 1 ,  2 )} . (4) For simplicity, we drop the '10' subscript in the base-10 logarithm for this expression and hereafter.A Gaussian process is a probabilistic model that represents a function as a probability distribution over all possible functions that are consistent with a given relation.It is fully specified by its mean function (in our case, log SMHM( h , )) and covariance function, or kernel,  ( 1 ,  2 ).The mean function provides the expected value of the function at any epoch, and the covariance function determines how the values at different epochs are related to each other.We choose a squared exponential kernel with the scatter of the Behroozi et al. (2013b) relation:  () = 0.218+ 0.023 /(1 + ).The parameter  characterizes the autocorrelation timescale.The limit of  → 0 represents a pure Gaussian noise.
We set  = 2 Gyr to reflect a typical gas depletion timescale in galaxies.However, we have verified that any value in the range  = 1 − 4 Gyr does not noticeably affect any property we analyze.The Gaussian process technique solves all the caveats of the previous methods.In addition, the non-zero autocorrelation preserves memory of  ★ at previous epochs, leading to a smoother evolution of  ★ .
This provides a monotonic function  ★ () during the peak of GC formation at  = 1 − 6.At  ≲ 1 when the galactic star formation rate (SFR) drops and the growth of  h slows, the calculated value of  ★ () may also decrease.In this case, we set  ★ to equal the previous value as we do not expect the stellar mass to decrease in reality (aside from the mass loss due to stellar evolution, which is already incorporated in the SMHM relation).
Lastly, we update the galactic MZR to the following expression: This relation is updated from the previous version (Choksi et al. 2018) by re-calibrating with the new observational data spanning a broad range of mass ( ★ = 10 7 − 10 11 M ⊙ ) and redshift ( = 0.7 − 12).The low- ★ and high- data are mainly from the recent programs obtained with the James Webb Space Telescope (JWST), see Appendix A for a detailed description.
Similarly to SMHM, we model the scatter of [Fe/H] and its evolution via the Gaussian process, where  g is also a squared exponential function, The scatter of MZR is characterized by the parameter  g , which we set  g = 0.3 dex to be consistent with the observed scatter (Appendix A).Note that [Fe/H] g refers to the mean metallicity of the galaxy.Clusters formed within this galaxy do not necessarily inherit exactly the same metallicity because of spatial variation and gradients within the ISM.We add an additional Gaussian noise to Eq. ( 7) to account for the internal metallicity dispersion: where [Fe/H] c stands for the metallicity of individual clusters formed within a galaxy of stellar mass  ★ at redshift .The internal scatter is quantified by the parameter  c .While  g can be measured directly from galaxy surveys, we must calibrate  c to match the observations of nearby GC systems.We achieve this by running the model with different  c to search for the model realizations that reproduce the total metallicity dispersion observed in the Virgo Cluster Survey (Peng et al. 2006).We find that  c = 0.2 − 0.35 dex can match observations within one standard deviation.We adopt the smaller  c = 0.2 dex in this work.

Modifications to cluster evolution
In the cluster evolution step, we follow the trajectories of GC particles taking into account two main processes of mass evolution: stellar mass loss and tidal disruption.Since most mass loss due to stellar evolution happens in the first tens of Myr, which is much shorter than the lifetime of a typical GC (∼10 Gyr), we treat stellar evolution as an instantaneous mass loss at formation: where  0 is the cluster mass at formation and  i is the cluster mass after stellar evolution.The remaining mass fraction  sev is a function of the stellar initial mass function (IMF) and metallicity.
Here we assume that the IMF is constant for all clusters.Metallicity affects the exact duration of the stellar evolution but not the final remaining fraction (only by ∼ 1% for −4 < [Fe/H] < 1).Therefore,  sev is close to a constant.We follow Gieles & Gnedin (2023) to set  sev = 0.55.After accounting for stellar evolution, we can consider  i as the 'initial' mass before tidal disruption.Following CG23 we express the tidal disruption rate of a cluster with mass  as with the parameters  = 2/3 and  = 4/3 that match the low-density -body models of Gieles & Gnedin (2023).This is different from CG22 where we used  =  = 2/3 (which is appropriate only for very concentrated clusters).
The disruption rate also depends directly on the local tidal field strength (Gieles & Baumgardt 2008).We parameterize the tidal field strength by the angular frequency Ω tid via the effective eigenvalue  1,e that takes into account the centrifugal, Euler, and Coriolis forces (Renaud et al. 2011): where  1 ,  2 , and  3 are the eigenvalues of the tidal tensor in descending order.This expression describes the mass loss rate more accurately than the expression in CG22.
We derive the tidal tensor numerically following the same method as in CG22: we compute the second order derivatives of the gravitational potential on a 3 × 3 × 3 cubic grid centered on the GC tracer particle, with side length = 300 pc.Although the side length is still too large compared to a typical tidal radius of GCs, we cannot adopt a lower value because of limited spatial resolution of simulations.To distinguish between the true eigenvalues  and those we derive numerically, we use the notation λ for the latter.To correct for the systematic deviation of the derived value from the true Ω tid , we use the third adjustable model parameter  as a correction:

GALAXY SELECTION
To find GC systems in TNG50 and the LG simulations that can match the observational properties of the GC systems in the MW and M31, we first select two samples of galaxies that have properties similar to the MW and M31, respectively.Next, we apply our model to these galaxies to obtain the model GC systems at present-day.Finally, by comparing the model GC systems to the MW and M31 systems specifically, we rank the model GC systems and output the three best analogues to the MW and M31.In the following subsections, we describe the criteria to select the galaxy samples and rank model GC systems.

Galaxy samples
We apply the following criteria to select MW analogues in TNG50: • Galaxies with total mass1 between 10 11.9 and 10 12.3 M ⊙ .
• Galaxies with maximum circular velocity  c,max between 210 and 270 km s −1 .

M31 M32 progenitor
Figure 1.Mass growth histories of the main progenitor branch for the samples of 10 MW analogues (left) and 14 M31 analogues (right).The dashed vertical lines label the epochs for key assembly events required by our selection criteria, including the GS/E-like merger (10 − 12 Gyr) and the subsequent quiescent stage (< 10 Gyr) for MW, and the M32-progenitor-like merger (< 6 Gyr) for M31.We show each galaxy as a thin gray curve and highlight the three best-matching galaxies in each set with thick colored curves.We keep the same color scheme in all plots hereafter when referring to these galaxies.
• Galaxies with at least one major merger between 10 − 12 Gyr ago to match the accretion of the GS/E satellite.
• Galaxies with no major merger in the last 10 Gyr.
• Galaxies formed 25 − 35% of their present-day stellar mass (calculated from halo mass via the SMHM relation of Behroozi et al. 2013b) at  lookback = 10 Gyr.
Our virial mass range includes the result (9×10 11 M ⊙ ) of the most recent modeling of the Sagittarius stream (Vasiliev et al. 2021), but some deviation at the virial radius is expected because the simulated halos are not likely to have exactly the same density profile as the MW.We keep the mass variation within a factor of ∼ 2. We also add the circular velocity criterion to match the inner mass distribution, which may be more relevant for modeling star and cluster formation.The range of  c,max is chosen to match the observed rotation curve (Eilers et al. 2019) with a variation ±30 km s −1 .The circular velocity value at Sun's location at 8.5 kpc from the center is constrained to be within 200 − 240 km s −1 .The stellar mass constraint at 10 Gyr is from Leitner (2012).
We define a major merger as the mass ratio less than 4:1.That is, the total mass of the incoming galaxy is greater than 1/4 of the main galaxy when the incoming galaxy reaches its maximum mass.Our two merger criteria select galaxies that assembled early.To illustrate this, we plot the mass growth histories of the main progenitor branch of these galaxies in the left panel of Fig. 1.In TNG50, there are 6 galaxies that match the above criteria.Additionally, we include one more galaxy (ID 519311 2 ) from the sample of 'early-spin-up' MWlike disk galaxies by Semenov et al. (2023).Our best MW analogue (523889) is also included in their sample.
We include the Louise, Romeo, and Juliet galaxies from the LG simulations since they have MW-like mass assembly histories.The small group environment of the LG simulations provides a more realistic background for the satellite accretion.In total, the MW sample has 10 (7 from TNG50 + 3 from LG simulations) galaxy candidates with similar properties to the MW. 2 The galaxy ID refers to the SubhaloID in TNG50 subfind catalogue.
For the M31 sample, we follow the criteria below: • Galaxies with total mass between 10 12 and 10 12.5 M ⊙ .
• Galaxies with at least one major merger in the recent 6 Gyr to resemble the recent merger with the M32 progenitor as suggested by D' Souza & Bell (2018).
Since M31 contains a much richer system of GCs and its total mass is less certain than in the MW, we allow a wider range of mass.We find 14 analogues of M31 satisfying the criteria.In the right panel of Fig. 1 we show the mass growth histories of the main progenitor branch of these galaxies.The M31 analogues are typically assembled later and have more variable mass growth at late times compared to the MW analogues.This is directly linked to the last selection criterion.
When selecting samples of the MW and M31 analogues, we only focus on the mass assembly history and do not take into account any baryonic properties such as the luminosity, surface brightness, or metallicity (even the stellar mass criterion for selecting MW analogues is derived from  h using SMHM relation).This minimizes the dependence on the specific prescriptions used in the TNG50 model.However, as we show later, even the differences only in the mass assembly history can lead to a wide range of properties of the GC systems.Some of these realizations correctly match the observed properties of GC systems in the MW and M31, suggesting that GC formation is strongly related to the hierarchical assembly of galaxies.

Model calibration
Before going into details about how we calibrate the model parameters, we introduce the observational data with which we compare the model systems.
The mass and spatial distributions for the MW GCs are from the third version of the Hilker et al. ( 2019) catalogue3 .The metallicities are from the 2010 version of the Harris (1996) catalogue.
The magnitude and position data for the M31 GCs are from the Revised Bologna Catalogue 4 (RBC) of Galleti et al. (2004, version V.5, August 2012).We compute the cluster mass from the V magnitudes using a mass-to-light ratio / V = 1.83 M ⊙ / V,⊙ (Baumgardt et al. 2020).The metallicities are from the LAMOST spectroscopy survey of star clusters in M31 (Chen et al. 2016) using spectral fitting with the models of Vazdekis et al. (2010).
We limit the MW sample to clusters with  > 10 4 M ⊙ to avoid the potentially incomplete low-mass clusters.For the M31 sample we require the clusters to have  > 10 4.5 M ⊙ and locate outside the central 1 kpc (projected radius) annulus to avoid contamination near the center of M31.We also apply the same criteria to the model GC systems to consistently compare them with the observations.Since the inclination angle of the M31 galaxy is 77 • (Simien et al. 1978), we project the model coordinates for M31 analogues using the same inclination angle with respect to the disk plane.We define the disk plane using the direction of the angular momentum of all stellar particles in the galaxy.We always refer to this inclination angle for projected radius throughout the rest of the work, unless specified otherwise.
Our model has three adjustable parameters:  2 ,  3 , and .To find the best values for these parameters, we need to calibrate the model to match important observable features.For this purpose we use a merit function that quantifies the discrepancy between the model and observations: where G  is a gauge function of the -th simulated galaxy.We calculate the merit function separately for the MW and M31 samples.We take into account 5 features of the GC system: total GC mass, the metallicity distribution, mass function, radial distribution, and velocity dispersion.They correspond to the 5 independent terms: The first term is the  2 function to evaluate the total mass of surviving clusters, where  GC represents the total GC mass of the -th galaxy, and  GC,obs = 3.7 × 10 7 M ⊙ when compared with the MW.This number comes from the sum of individual MW cluster masses greater than 10 4 M ⊙ .However, since the M31 data are incomplete below 10 4.5 M ⊙ , we cannot simply sum up the mass of GCs greater than 10 4 M ⊙ .Instead, we fit the M31 mass function greater than 10 4.5 M ⊙ with a log-normal function, and integrate this function down to 10 4 M ⊙ to get the total GC mass.This gives a correction ≲ 2% and yields  GC,obs = 9.5 × 10 7 M ⊙ for M31.The denominator is the uncertainty of the number of clusters.We set   = 0.2 dex to approximate the uncertainty of total GC mass.
The second term is similar to the first term but evaluates the velocity dispersion of GCs, where  GC represents the 3D velocity dispersion of all surviving Inside the summation,   ∈ {  ,   ,   } stands for the Kolmogorov-Smirnov (KS) -value for the mass function, metallicity distribution, and radial distribution.The radial distribution refers to the 3D galactocentric radius when calibrating for the MW.It refers to the projected radius for the M31.The Heaviside  function returns 1 if   is greater than a threshold value  th , otherwise 0. We employ  th ∈ {0.1, 0.03, 0.01}.In other words,   > 0.1 contributes 3 points to the final G, 0.03 <   < 0.1 contributes 2 points, and 0.01 <   < 0.03 contributes 1 point.For example, the greatest possible G is 9 if all   ,   ,   > 0.1,  GC =  GC,obs and  GC =  GC,obs .We employ three threshold values instead of one (as in our previous model versions) to provide a finer grading hierarchy with more possible grades contributed by the three KS tests.This finer grading system avoids the discreteness that too many galaxies receive the same grade, and thus allows more nuanced calibration and selection.
For any set of model parameters, we can compute G of all MW analogues by comparing their properties to the MW.The average of G yields the merit function for the MW sample.Similarly, we can obtain the merit function for the M31 sample by comparing to M31.We then span the parameter space and select the parameters that maximize the merit functions.
It is worth noting that due to the different numerical resolution and output frequencies in TNG50 and the LG simulations, the best model parameters for the two sets of simulations are not necessarily the same.Therefore, we calibrate the model parameters separately for the two simulation sets.For TNG50, we first compute the merit functions on a grid of model parameters.The best parameters for the MW and M31 differ slightly but systematically: the merit function of the M31 maximizes at greater  2 and smaller  compared to that of the MW.However, the merit function flattens near the maximum, allowing us to vary the parameters near the peak values without significantly affecting the model performance.We then select a single set of best parameters between the peak values of both merit functions for the MW and M31 samples: (  2 ,  3 , ) = (18, 0.5 Gyr −1 , 1.5).
Since the LG simulations only contain MW analogues, we simply compute the merit function for the three MW-like galaxies and select the parameters that maximize the merit function.This gives (  2 ,  3 , ) = (14, 0.5 Gyr −1 , 1.5).These values differ from our previous values (CG23, which has 3 = 0.7 Gyr −1 ) because we have updated the model setups, especially the new prescription for cluster evolution, which effectively leads to stronger tidal disruption.To balance the decrease of clusters due to disruption, we need a smaller  3 to enable more GC formation at later times.

SELECTED MW AND M31 ANALOGUES
We run our model 50 times on each model galaxy with different random seeds to generate an ensemble of 50 realizations.Having multiple random realizations offers two major advantages.First, the model randomness, which includes the scatter in scaling relations and the stochasticity involved in sampling clusters from the ICMF and We highlight the three best-matching galaxies in each panel with the colored curves.For comparison, we overplot the observational data for the MW and M31 systems as thick black curves.
assigning clusters to simulation particles, can result in GC systems with diverse properties when the same model is rerun for the same galaxy.Analyzing the average properties of these 50 realizations reveals their systematic dependence on the model's input and the model itself.Furthermore, having multiple random realizations enables us to search for the best MW/M31 analogues from significantly larger samples.For each MW/M31 analogue, we calculate G for the 50 realizations and identify the one with the highest G as the best representative.In this section, we present the three best-matching representatives from all sample galaxies, considering these galaxies as the best MW/M31 analogues.

Properties of globular clusters in model galaxies
The three best MW analogues are Romeo from the LG simulations and 523889 and 519311 from TNG50.The IDs of the three best M31 analogues are 532301, 474008, and 441709.In Fig. 1, we plot the mass growth history of these analogues as colored curves.The G function of the three M31 analogues is generally lower than that of their MW counterparts, indicating that our model is more effective in matching the MW system compared to M31.We note that the two TNG50 MW analogues and three M31 analogues are also included in the TNG50 MW/M31 sample by Pillepich et al. (2023), who selected 198 galaxies based on stellar mass, stellar morphology, and environment.
Next, we study the distributions of important GC properties.In the upper panel of Fig. 2 we plot the cumulative mass function for clusters in each galaxy.The GC systems of MW analogues in TNG50 have slightly more abundant low-mass clusters ( ≲ 10 5 M ⊙ ) compared to the observations.This is likely because the tidal disruption is underestimated for the low-mass clusters due to the limited spatial resolution of TNG50 ∼ 100 pc.We do not find such a deviation for Romeo because the LG simulations have higher spatial resolution ∼ 60 pc, and because the surviving GCs in Romeo are formed mainly in-situ, which subjects them to stronger tidal disruption.On the other hand, the mass functions of all M31 analogues are consistent with the mass function of M31 for  > 10 4.5 M ⊙ .
In the middle panel of Fig. 2, we plot the cumulative cluster metallicity distribution.Our model can match the observed metallicity distributions of both MW and M31 very well.The observed MW metallicity distribution completely falls into the 1- region of all MW sample galaxies.Our model also predicts that the M31 clusters have systematically broader metallicity distribution with higher abundance of metal-rich clusters compared to the MW.This is likely due to the younger GC population in M31 caused by the last major merger.
The abundance of [Fe/H] ≲ −1.5 clusters in the full M31 sample is slightly higher but still within the 1- region of the observations.Since [Fe/H] ≲ −1.5 corresponds to GCs formed in smaller galaxies and at early epochs ( ≳ 3, see Eq. ( 6)) such a mismatch may indicate the need for additional constraints on the early assembly history of M31 (our M31 selection criteria have no constraint at  lookback > 6 Gyr).
Finally, we plot the cumulative distribution of distance from the galaxy center in the lower panel of Fig. 2.Although the radial profile of Romeo can closely match the observed profile of the MW GCs with KS -value > 0.1, our model applied to TNG50 tends to produce more spatially extended GC distribution than the MW system.The half-number radii for the two best-matching MW analogues in TNG are 8 − 9 kpc, which is close to twice the observational value ∼5 kpc.Similarly, we find that only 532301 can match the radial profile of M31 GCs.The other two best-matching M31 analogues are also spatially more extended than the M31 GC system: the half-number radius of M31 GCs is ∼6 kpc, whereas the two TNG50 galaxies yield 8 − 9 kpc.The only model setting directly related to the spatial distribution is the initial boundary radius to distribute newly formed clusters.However, as we have verified in Sec.2.2, changing this parameter does not improve the final result significantly.As we discuss later in Sec.5.2, this mismatch is likely because TNG50 galaxies have more abundant mergers than the galaxies in the LG environment.

Cluster catalogues
Based on the above comparison, the three MW analogues and three M31 analogues are consistent with all the observable properties of their GC systems.We release the catalogues of model GCs on our model site www.github.com/ognedin/gc_model_mw.The catalogues include the following properties: galaxy ID, cluster formation/disruption/accretion time, cluster mass at formation/current, position, velocity, orbital actions, pericenter/apocenter radii, gravitational potential, metallicity, progenitor galaxy IDs, galaxy total mass at formation/current, and galaxy stellar mass at formation/current.We summarize these properties in Table 1 with a brief description for each entry.
Except for the orbital parameters (orbital actions, pericenter/apocenter radii, and gravitational potential) all other properties are direct outputs of the model.We describe the calculation of orbital parameters in the following section.

Orbital properties
We calculate the orbital properties of GCs using the same methods as in CG22.We first employ the agama code (Vasiliev 2019) to model the galactic potential with the multipole expansion and spline approximation methods.Since the potential of MW can be described by spheroids and disks (e.g., McMillan 2017), we model the presentday potentials of TNG50 galaxies with these two components.A largely spheroidal DM potential is modeled by the axisymmetric spherical harmonic expansion, where    are the real-valued spherical harmonics.We only take into account the axisymmetric terms with  = 0.The -dependent coefficient Φ  () is calculated on 20 grid points evenly spaced in log .The agama code employs a quintic spline to connect Φ  () values on grid points.
The disky baryonic (stars + gas) potential is modeled as an axisymmetric form in a cylindrical coordinate system Φ(, ).We use a 2D quintic spline to approximate the potential on a 20 × 20 grid evenly spaced in the log -log  space.
We find  max = 2 is sufficient to describe the simulation-provided potential to < 2% accuracy.Such an error is so small that we can ignore its influence on the subsequent calculation of the orbital parameters.
Next, we input the present-day positions and velocities of model GCs to agama to calculate the orbital actions and pericenter/apocenter radii.The orbital actions of a closed orbit are defined as where  ∈ {, , } is the radial, azimuthal, and vertical coordinates.  and   characterize the oscillations in radial and azimuthal directions, whereas   equals the  component of angular momen-tum   .In separable potentials5 , the actions are functions of the IOMs,   =   ( 1 ,  2 ,  3 ), where  1 =  is the total energy.The agama code computes actions using the Stäckel fudge method (see Sec. 3.2 in Vasiliev 2019).This approach assumes the potential is a Stäckel form (which is not necessarily true) so that we can analytically find the second and third integrals: the second integral in this case is   , while the third integral is non-classical.agama speeds up the calculation by using a pre-computed interpolation table of   =   ( 1 ,  2 ,  3 ).Vasiliev (2019) showed that this approximation is accurate to 90 − 99% even for extremely eccentric orbits.We calculate the pericenter/apocenter radii by computing the closest/farthest possible distances the cluster can reach with the current energy and angular momentum.They are the two roots of the equation  = Φ(,  = 0) +  2 /2 2 , where Φ(,  = 0) is the in-plane axisymmetric potential.
We do not calculate orbital properties for Romeo because it is from a collisionless simulation with no baryon physics.Since the baryonic matter should become dominant near the galaxy center where most in-situ clusters reside, the orbital parameters would be severely miscalculated.Therefore, all entries for the orbital properties in the Romeo catalogue are set to zero.

Progenitor branch ID
In the simulation merger trees, it is common for the same galaxy at different snapshots to be assigned different IDs.It is thus more appropriate to refer to this galaxy as a 'branch' rather than a 'galaxy' within the terminology of merger trees.To establish a clear and unique identifier for each branch, an effective approach is to label this branch with the ID of the 'main leaf' galaxy.The main leaf galaxy corresponds to the first instance of the galaxy in this branch.This labeling technique guarantees a one-to-one mapping between each branch and its main leaf ID without any ambiguity.
We provide the host_id_form and host_id_accrete entries to identify the progenitor branches for the ex-situ clusters.The host_id_form specifies the main leaf ID of the host galaxy at the time of cluster formation.This host galaxy may either merge directly into the central galaxy or first merge with another satellite, which is subsequently accreted onto the central galaxy.In the first case, the host galaxy at cluster formation is the same as the galaxy that ultimately delivered the cluster to the central galaxy.However, in the second scenario, the two are different.To account for this, we introduce host_id_accrete as the main leaf ID of the last satellite responsible for bringing the cluster into the central galaxy.We also provide the t_accrete entry specifying its accretion time.
In a major merger event, the accreted satellite is likely to have its own merger history prior to joining the central galaxy.
Since we expect all GCs from such a satellite to exhibit similar kinematics at present, regardless of their specific progenitors, the host_id_accrete entry provides a more straightforward way to query all GCs brought to the main galaxy by that satellite.Conversely, during the early stages of galaxy formation when the main galaxy was not dominant in its local environment, distinguishing all satellites becomes important.In such a scenario, the host_id_form entry is more relevant.This entry is also useful for studying the multiple populations of GCs in the same progenitor galaxy, under the assumption that GCs originating from different satellites lead to distinct populations.

Impact of the last major merger of M31 on the cluster spatial distribution
In our selection, the M31 sample is characterized by a major merger in the last 6 Gyr, whereas the MW analogues did not experience any major merger in the last 10 Gyr.The former correspond to the 'late-assembled' galaxies, while the latter correspond to the 'earlyassembled' galaxies.In this section, we investigate the impact of the last major merger on the observable properties of the resulting GC population.Without loss of generality, we study the recent merger event of the best-matching M31 analogue 532301 as an example.All other M31 analogues lead to similar conclusions.In Fig. 3, we plot the face-on (not 77 • inclination angle) projection of surviving GCs (defined as  > 10 4.5 M ⊙ ) before the merger, at the satellite apocenter after the first passage, and at the end of the merger.The merger happens at  lookback = 3 − 4 Gyr, with a merger ratio close to unity.It can be easily noticed that the clusters brought by the recent merger have broader spatial distribution than the in-situ clusters at the end of the merger.
To quantify such a difference in the spatial distribution, we split the clusters into three categories: clusters formed in the satellite prior to the merger (or 'satellite' clusters), clusters formed or accreted into the central galaxy prior to this merger ( form > 6 Gyr, or 'central' clusters), and new clusters formed during this merger ( form < 6 Gyr, or 'new' clusters).We plot the evolution of the face-on radial profiles for the three categories in Fig. 4. The radial profiles are obtained by kernel density estimation using an Epanechnikov (1969) kernel:  () = 0.75[1 − (/ℎ) 2 ] with ℎ = 0.25 dex.The peak at  ∼ 1 kpc in the middle panel is contributed by just one cluster captured by the central galaxy during the first passage.Since most satellite clusters are still bound to the satellite at this stage, we do not analyse this outlier further.
As the satellite approaches the main galaxy around  lookback = 3.6 Gyr, it does not largely alter the spatial distribution of central clusters ≲ 10 kpc (approximately the pericenter distance of the satellite's first passage).In contrast, it pulls the clusters in the outer region ≳ 10 kpc outwards, as the satellite is massive enough that it perturbs the potential of the central galaxy significantly in this region.The formation of about 110 new clusters increases the number density within ≲ 3 kpc.However, strong tidal disruption quickly brings the numbers down by the end of the merger ( lookback = 3.2 Gyr).In contrast to such dramatic early disruption, the number of new clusters only drops gradually in the next 3.2 Gyr until the present.Similarly, the early tidal disruption removes around 40 central clusters during the merger, whereas the subsequent evolution only disrupts less than 30 central clusters in the next 3.2 Gyr, leaving ∼90 central clusters surviving until the present.
Despite the change in the normalization, the radial profile of the central clusters does not change significantly during the merger.The merger only raises the effective radius of the central clusters from 2.5 kpc to 3.7 kpc.However, the satellite clusters have a radial profile significantly more out-spread than the centrals: the effective radius of the satellite clusters is 6.1 kpc at the end of the merger, raising the overall effective radius to 4.6 kpc.
In the subsequent 3.2 Gyr, the radial profile does not change much, except for the inner part being lowered by tidal disruption, slightly increasing the final effective radius to 5.2 kpc.
By adding a large number of ex-situ clusters with distinct radial distribution from the central cluster population, the recent major merger can significantly enlarge the GC effective radius by a factor of 2 in less than 1 Gyr.However, the merger does not significantly alter the spatial distribution of central clusters ≲ 10 kpc.Also, the merger triggers the formation of a young in-situ GC population within ≲ 3 kpc.Combining these two factors may explain why many M31 clusters are still located around the disk without being stripped away by the merger.But as we show below in Sec.5.1, the last major merger can still heat up the in-situ clusters to higher energy.

Investigating galaxy assembly with the catalogue
The main goal of constructing this catalogue is to investigate connections between the galaxy assembly history and the properties of GCs.For example, we show the age-metallicity relations (AMR) of surviving clusters for a typical MW analogue (523889) and a typical M31 analogue (532301) in Fig. 5. GCs in the MW analogue are in general older than the M31 counterparts because of the early assembly history of MW.This leads to almost no GC formation after  lookback = 10 Gyr.In contrast, GCs in M31 span a wider range of ages.For 532301 specifically, there are four bursts of GC formation: 1) ∼ 40 ex-situ GCs formed at  lookback ≈ 12 Gyr, corresponding to active cluster formation in the early universe when galaxy mergers are frequent; 2) subsequent 'bottom-up' assembly of satellites contributed a significant amount of mass to the main galaxy, leading to the formation of majority of in-situ GCs at ≈ 10 Gyr; 3) similar Each dot represents a surviving cluster either formed in-situ (red) or ex-situ (blue).We calculate the contours using Gaussian KDE with bandwidth = 0.1 dex for horizontal axis, and 0.3 Gyr for vertical axis.The contours from dark to light enclose 10, 50, and 90% of total number, respectively.
growth of the largest companion galaxy at ≈ 8 Gyr, forming ∼ 130 GCs; 4) the major merger between the main galaxy and its companion boosted in-situ GC formation at ≈ 4 Gyr (as illustrated in Sec.4.3).
Other M31 analogues also have more extended and discrete GC formation histories, but the timing and order for each GC population may differ among galaxies.
Another noticeable difference is the larger metallicity separation between in-situ and ex-situ clusters in the MW analogue.Because the assembly of MW is less hierarchical than that of M31, the main progenitor galaxy of MW is more dominant in the GC populations.This enlarges the gap between the metallicity of clusters formed in the main galaxy and the satellites.In contrast, most of ex-situ clusters in the M31 analogue formed in the largest companion galaxy which had mass and metallicity comparable to the main galaxy.The metallicity of these ex-situ clusters is thus similar to the in-situ clusters formed before the merger.
In addition to the AMR, galaxy assembly histories shape the orbits of GCs.For example, Fig. 6 shows the normalized IOM space for the MW and M31 analogues (523889 and 532301).The circularity parameter  is defined as   / circ (), where  circ () is the angular momentum of an in-plane circular orbit with the same energy  as the cluster.The normalized energy  is defined as /| 0 |, where  0 is the gravitational potential at the galaxy center.We have  ∈ [−1, 1] and  ∈ [−1, 0).Perfectly circular and in-plane orbits have  = 1 (prograde) or  = −1 (retrograde), while purely radial or polar orbits have  = 0.
We find that the ex-situ populations in both MW and M31 analogues have similar distribution around (, ) = (0, −0.3).However, the in-situ populations have quite different distributions: in the MW in-situ GCs are located near (, ) = (0.5, −0.6), indicating significantly lower-energy and more disky orbits.The in-situ GCs of the M31, however, have a more similar distribution to the ex-situ population and center around (, ) = (0.2, −0.Each dot represents a surviving cluster either formed in-situ (red) or ex-situ (blue).We calculate the contours using Gaussian KDE with bandwidth = 0.2 for horizontal axis, and 0.1 for vertical axis.The contours from dark to light enclose 10, 50, and 90% of total number, respectively.
ulation, the distributions of in-situ and ex-situ clusters in the M31 analogue largely overlap in the normalized IOM space.This is likely due to the recent major merger which heated up the in-situ clusters to higher energy and moved some of them away from the disk.Since GCs formed in different progenitor galaxies are located in distinct regions within the chemical and kinematic property space, numerous attempts (e.g., Massari et al. 2019;Malhan et al. 2022;Belokurov & Kravtsov 2023) have been made to identify GCs of various origins by classifying isolated clusters/groups in the AMR, IOM, and chemistry spaces.However, these studies often yield conflicting results.Our catalog can serve as a tool to assess the accuracy of the classification algorithms employed in these studies since our model provides 'true' labels from the simulations (as indicated by the host_id_from and host_id_accrete entries in Table 1).Furthermore, our catalog offers a unique opportunity to investigate the spread of properties among GCs originating from the same progenitor and to explore how much GCs retain their original kinematic signatures during galaxy assembly.Finding the answers to these questions is crucial for characterizing the accuracy of current classification schemes.

Why are TNG50 GC systems more extended?
As shown in Sec.4.1, the GC systems in TNG50 galaxies tend to have more radially extended distribution than the observations: the halfnumber radii of the two best-matching MW analogues in TNG50 are around 1.7 times the observational value ∼5 kpc.In fact, all 7 TNG50 galaxies in the MW sample have greater half-number radii, with a median value ∼12 kpc.We find a similar trend for the M31 sample.In contrast, all three MW-like galaxies in the LG simulations can match the radial profile of MW GCs with half-number radii ∼5 kpc.
These extended distributions of the TNG50 systems may be a general outcome of the large-scale environment probed by the TNG50 simulation.For example, Semenov et al. (2023) showed that approximately 90% of TNG50 galaxies developed their disk later than the MW.This difference arises from TNG50 lacking a quiescent, merger-free environment, which is necessary for the survival of an early-formed disk.While we limit our MW sample galaxies to those without major mergers with a mass ratio less than 4:1 in the last 10 Gyr, we do not impose constraints on higher merger ratios (∼10:1) which can still play a significant role in structure formation.
Additional evidence indicating that MW may be less mergerdominated than typical TNG50 galaxies lies in the high in-situ fraction of MW GCs.Observational classification studies suggest that the in-situ fraction may range from about 40% (Massari et al. 2019) to 56% (Malhan et al. 2022) and up to 66% (Belokurov & Kravtsov 2023).In contrast, the TNG50 galaxies in the MW sample typically have in-situ fractions between 20% and 50%, once again suggesting that even the selected MW sample galaxies may be too merger-dominated compared to the real MW.
The relatively low fraction of in-situ clusters in TNG50 leads to two discrepancies between the model and observed data.First, the low in-situ fraction accounts for the too extended distribution of TNG GCs, as shown in the lower left panel of Fig. 2.This is because the median radius of ex-situ clusters is approximately four times larger than that of the in-situ clusters (see CG22).Second, the low in-situ fraction in TNG50 also explains why the AMR of TNG GCs lacks a sequence of old, metal-poor in-situ clusters (as shown in the left panel of Fig. 5), which is a robust feature of the MW assembly history.To align the simulations more closely with observed data, a more realistic simulation of the MW assembly is needed, particularly to match the observed AMR of the in-situ clusters.

SUMMARY
We introduce a catalogue of model star clusters in the MW and M31 analogues, drawn from the TNG50 and LG simulations.By applying criteria based on galaxy virial mass, circular velocity profile, and defining assembly events, we select a sample comprising 10 MW analogues and 14 M31 analogues.We then apply an analytical model of GC formation and evolution to the merger trees and particle outputs of these simulated galaxies.This enables us to obtain key observables such as the mass, age, metallicity, positions, and velocities of the surviving cluster population.We calibrate the model parameters to optimize agreement with the observed total cluster mass, mass function, metallicity distribution, radial profile, and velocity dispersion.From these model results we select the three best-matching GC systems each for MW and M31.One of the best-matching MW analogues is from the LG simulations; the other five are from TNG50 (see Fig. 1 for the mass growth histories of the six galaxies).
The GC systems in the three best-matching MW analogues successfully reproduce the observed mass function and metallicity distribution (Fig. 2) after model calibration.The galaxy from the LG simulations also matches the observed radial distribution, yielding a KS -value > 0.1.However, the radial profiles of the other two GC systems from TNG50 are more extended, with effective radii approximately 1.7 times the observed value.Likewise, the three bestmatching M31 analogues have mass and metallicity distributions that agree with observational data.However, only one M31 analogue reproduces the observed radial distribution, while the remaining two exhibit more extended profiles, with effective radii ∼1.5 times the observed value.This discrepancy may arise from the higher rate of mergers in TNG50 compared to the LG environment.Such mergers typically introduce a significant population of ex-situ clusters, which tend to be located ∼4 times father away from the galaxy center than their in-situ counterparts.
An in-depth study of the last major merger in M31 analogues also supports the above conclusion.We find that even a 1:1 merger has a limited impact on the spatial distribution of in-situ clusters within ≲ 10 kpc (Fig. 4).However, such a merger can notably enlarge the overall GC effective radius by a factor of 2 in a short period ≲ 1 Gyr, primarily by bringing in a large number of ex-situ clusters.Moreover, the merger triggers the formation of a new population of young insitu clusters near the galactic center.This population, in addition to older clusters brought in by the merger, forms distinct groups in the AMR (Fig. 5).These groups have lower age and higher metallicity compared to the central clusters formed prior to the merger, resulting in an AMR spanning a wide range in both age and metallicity.
Furthermore, we show that the galaxy assembly history significantly influences GC kinematics.The MW analogues have a quiescent formation history over the past 10 Gyr.Such a long period preserves the original location of clusters in the IOM space (Fig. 6).In this space, in-situ clusters tend to have lower energy and angular momentum   closer to  circ (), whereas ex-situ clusters generally have higher energy and   ≈ 0. This distinction is particularly helpful when using the GC properties to decode the merger history of their host galaxy.Conversely, the last major merger of M31 elevates the in-situ populations to higher energy, making the in-situ and ex-situ populations less distinguishable in the IOM space.
We make our catalogue available in a machine-readable format.Along with all the variables directly output by the model, the catalogue also includes several orbital parameters: the gravitational potential, actions, and pericenter/apocenter radii.A comprehensive list of keys and descriptions for these variables is provided in Table 1.The catalogue can be accessed at www.github.com/ognedin/gc_model_mw.Points with black errorbars show data for individual galaxies (Jung et al. 2023;Curti et al. 2022;Boyett et al. 2023;Williams et al. 2023;Hsiao et al. 2023).The observational compilations that fit the data with functional forms are shown as colored curves or bands (Maiolino et al. 2008;Zahid et al. 2014;Lewis et al. 2023;Strom et al. 2022;Li et al. 2022;Sanders et al. 2021;Curti et al. 2023;Nakajima et al. 2023;Faisst et al. 2016;Matthee et al. 2023;Heintz et al. 2023).We convert the observed oxygen abundance 12 + log(O/H) to iron abundance [Fe/H] using Eq.(A5).The mean MZR adopted in this paper is represented by two black lines in each panel, evaluated at the lower/upper redshift bounds of the panel.The gray shaded ranges show the scatter  g = 0.3 dex around the bounding mean relations.

4
http://www.bo.astro.it/M31/GCs in -th galaxy, and  GC,obs = 200 km s −1 (calculated from the Hilker et al. 2019 catalogue of Galactic GCs) when compared with to MW.In M31, due to background contamination near the galaxy center, most velocity measurements are limited to the outer region.For example,Mackey et al. (2019) studied the kinematics of clusters ≳ 25 kpc and found  GC,obs = 134 km s −1 at  = 25 kpc.To properly compare the model M31 analogues with the observations, we only calculate  GC of the 32 GCs closest to 25 kpc.The intrinsic scatter of the dispersion can be approximated as   = 0.3 dex.

Figure 2 .
Figure 2.Cumulative distributions of the cluster mass (upper row), metallicity (middle row), and distance from the center (projected for M31; lower row) for the 50 realizations of 10 MW analogues (left column) and 14 M31 analogues (right column), represented by the gray shaded regions as the 16 − 84th percentiles.We highlight the three best-matching galaxies in each panel with the colored curves.For comparison, we overplot the observational data for the MW and M31 systems as thick black curves.

Figure 3 .
Figure3.Face-on projection plots of GCs during the recent major merger of the best-matching M31 analogue 532301.From left to right, the three panels correspond to the epochs before the merger, at the satellite apocenter after the first passage, and at the end of the merger, respectively.Each panel in the top row corresponds to a (200 comoving kpc) 3 cube centered on the main galaxy.We plot the DM column density as background using the quadtree-based projection code prj_plotter(Chen 2023).The in-situ and ex-situ clusters are shown as red and blue circles, respectively.The size of the circles increases with cluster mass.We also plot the trajectory (smoothed with cubic splines) of the incoming satellite as the blue curve.In the bottom row, we zoom-in to the central (20 ckpc) 3 region of each galaxy to show the GC distribution near the galactic center.

Figure 4 .Figure 5 .
Figure 4. Evolution of the GC number density profile during the recent major merger of 532301.From left to right, the three panels correspond to the epochs before the merger, at the satellite apocenter after the first passage, and at the end of the merger, respectively.We plot the profiles of all surviving clusters at each epoch as the solid black curve, with individual contributions of clusters brought by the satellite as blue, clusters formed in the central galaxy or accreted prior to this merger as gray, and new clusters formed during this merger as red.For comparison, we plot the present-day profile as the dashed curve.

Figure 6 .
Figure 6.Normalized integral of motion space for clusters in a typical MW analogue (523889, left panel) and a typical M31 analogue (532301, right panel).Each dot represents a surviving cluster either formed in-situ (red) or ex-situ (blue).We calculate the contours using Gaussian KDE with bandwidth = 0.2 for horizontal axis, and 0.1 for vertical axis.The contours from dark to light enclose 10, 50, and 90% of total number, respectively.

Table 1 .
Cluster properties in the catalogue.