The distribution of globular clusters in kinematic spaces does not trace the accretion history of the host galaxy

Context. Reconstructing how all the stellar components of the Galaxy formed and assembled over time by studying the properties of the stars that form it is the aim of Galactic archaeology. Thanks to the launch of the ESA Gaia astrometric mission and the development of many spectroscopic surveys in recent years, we are for the ﬁrst time in the position to delve into the layers of the past of the Galaxy. Globular clusters play a fundamental role in this research ﬁeld since they are among the oldest stellar systems in the MW and thus bear witness to its entire past. Aims. As a natural result of galaxy formation, globular clusters did not necessarily all form in the Galaxy itself. Indeed, a fraction of them could have been formed in satellite galaxies accreted by the Milky Way over time. In recent years, there have been several attempts to constrain the nature of clusters (accreted or formed in the Milky Way itself) through the analysis of kinematic spaces, such as the E − L z , L perp − L z , eccentricity − L z , and the action space, as well as attempts to reconstruct the properties of the accretion events experienced by the Milky Way through time from this kind of analysis. This work aims to test a widely used assumption about the clustering of the accreted populations of globular clusters in the integrals of motions space. Methods. In this paper we analyse a set of dissipationless N -body simulations that reproduce the accretion of one or two satellites with their globular cluster population on a Milky Way-type galaxy. Results. Our results demonstrate that a signiﬁcant overlap between accreted and ‘kinematically heated’ in situ globular clusters is expected in kinematic spaces for mergers with mass ratios of 1:10. In contrast with the standard assumptions made in the literature so far, we ﬁnd that accreted globular clusters do not show dynamical coherence, that is, they do not cluster in kinematic spaces. In addition, we show that globular clusters can also be found in regions dominated by stars that have a di ﬀ erent origin (i.e. a di ﬀ erent progenitor). This casts doubt on the association between clusters and ﬁeld stars that is generally made in the literature and is used to assign them to a common origin. By means of Gaussian mixture models, we demonstrate that the overlap of clusters is not only a projection e ﬀ ect on speciﬁc planes but is also found when the whole set of kinematic properties (i.e. E , L z , L perp , eccentricity, radial, and vertical actions) is taken into account. Overall, our ﬁndings severely question the recovered accretion history of the Milky Way based on the phase-space clustering of the globular cluster population.


Introduction
Our Galaxy, the Milky Way, is a collection of hundreds of billions stars mostly distributed in a disc which extends over distances of at least 20 kiloparsecs from the Galactic centre.The disc is surrounded by a diffuse stellar, oblate-shaped halo, which contains only few percent of the total stellar mass of the Galaxy (Hartwick 1987;Deason et al. 2011;Belokurov et al. 2018;Deason et al. 2019).A stellar bulge-bar sits at the centre of the Galaxy and dominates the mass, and light distribution, in the inner kiloparsecs (e.g.Wegg & Gerhard 2013;Ness & Lang 2016).All these components are made of stars which have, on average, different properties, in term of ages, chemical abundances, and kinematics, thus suggesting different scenarios for their formation.Reconstructing how all these stellar components formed and assembled over time, by studying the properties of the stars which make it, is the aim of Galactic archaeology.
In these last years, this field of research has seen a succession of fascinating new discoveries, thanks to the launch of the ESA Gaia astrometric mission1 and the delivery of its catalogues comprising full 6D positions and motions for millions of stars (33 million stars have full 6D phase-space information in the Gaia Data Release 3, see Gaia Collaboration 2022).In addition, many spectroscopic surveys such as the APOGEE survey with SDSS (Prieto et al. 2008), and the soon operational WEAVE@WHT (Dalton et al. 2012), MOONS@VLT (Cirasuolo et al. 2014;Cirasuolo et al. 2020) and 4MOST (de Jong et al. 2012) surveys aim at complementing Gaia by providing detailed chemical abundances for millions of stars, each.We are, for the first time, in the position to delve into the layers of the far and recent past of our Galaxy, bringing to light the timeline of events which helped make the Milky Way the galaxy that we observe today.Having survived billions of years to Milky Way's changes, globular clusters (hereafter GCs) observed today in the Galaxy bear witness of this entire past.GCs are dense, gravitationally bound stellar systems distributed from the Galactic centre to the most remote regions (up to 150 kpc from the galactic centre, see Harris & Racine 1979;Harris 1996;Meylan & Heggie 1997).They are very luminous objects, with typical masses of few 10 5 M , and very compact, having sizes of few pc only (Harris et al. 2013;Baumgardt & Hilker 2018;Gratton et al. 2019).Being very old, typically older than 10 Gyr (Marín-Franch et al. 2009;Carretta et al. 2010;VandenBerg et al. 2013), GCs are tracers of the early epochs of the Galaxy and, as such, they keep precious information about its formation and evolution.At present, more than 160 GCs are known in the Milky Way (see, for example, Baumgardt & Vasiliev 2021, for a recent derivation of their distances).
According to the ΛCDM cosmological paradigm, galaxy formation proceeds in a bottom-up scenario, as small structures merge together to build up the larger galaxies we observe today (White & Rees 1978).The Milky Way is a prime example of this formation mechanism, as demonstrated, for example, by the discovery of the Sagittarius dwarf spheroidal galaxy in the process of being accreted by our Galaxy and disrupted by its gravitational tidal field (Ibata et al. 1994;Newberg et al. 2002;Majewski et al. 2003).As a natural result of galaxy formation, not only field stars but also globular clusters may have been accreted (Penarrubia et al. 2009): a fraction of them can indeed have been formed in satellite galaxies accreted by the Milky Way over time, while the remaining fraction would have formed in-situ, that is in the Galaxy itself.This double nature of the Galactic globular cluster system is now fully accepted, and evidence of this duality come from the analysis of observational data (Marín-Franch et al. 2009;Forbes & Bridges 2010;Dotter et al. 2010;Leaman et al. 2013;VandenBerg et al. 2013) and of models, as well (Renaud et al. 2017;Kruijssen et al. 2019Kruijssen et al. , 2020;;Pfeffer et al. 2020;Chen & Gnedin 2022).In the recent years, there have been several attempts to constrain the origin of clusters (accreted or formed in the Milky Way itself), and to reconstruct the properties (mass, numbers) of the accretions events experienced by the Milky Way through time.Concerning globular clusters originally belonging to satellite galaxies that are still being accreted (e.g. the Sagittarius dwarf galaxy), it is possible to reconstruct their origin from the similarity between their positions, velocities and orbits, and those of stars belonging to the tidal stream tracing the on-going accretion (Bellazzini et al. 2020).For globular clusters accreted in the past instead, it is not feasible to proceed in this way since systems accreted long ago are expected to be spatially mixed with the in-situ population and, as a consequence, they should have lost all spatial coherence.Since the publications of the Gaia DR2 catalogue (Gaia Collaboration et al. 2018a), and the access to phase-space information for almost the entire clusters population (Gaia Collaboration et al. 2018b;Vasiliev 2019b), different works have tried to establish the accreted or in-situ nature of Galactic globular clusters using kinematic spaces, such as the energy-angular momentum space (hereafter E − L z , E being the total orbital energy and L z being the z component of the angular momentum space in a reference frame with the Galactic disc in the x − y plane).Under the hypothesis that accreted clusters conserve a dynamical coherence even several billion years after their accretion onto the Galaxy -i.e. that they cluster in kinematic spaces -and that different groups are related to different accretion events (Helmi & de Zeeuw 2000), Galactic globular clusters located in different regions of the E − L z space have been associated to different progenitors (see for example Massari et al. 2019).Since the models in Helmi & de Zeeuw (2000) assume an analytic and static Galactic potential and no dynamical friction is taken into account, the energies and angular momenta of the accreted satellites are overall conserved by definition.Thus, the fact that accreted globular clusters stand out as clumps in the E − L z space corresponding to different satellites, is, generally, a direct consequence of the assumptions made in the models, rather than an intrinsic feature of the accretion event.In addition to the E − L z space, the L perp − L z space (L perp being the projection of the total angular momentum onto the Galactic plane) has been used several times as a natural space to search for the presence of stellar currents, ever since the initial work of Helmi & de Zeeuw (2000).This has been done by assuming that L perp , although not strictly conserved, varies only marginally during the merging.In particular, the Helmi stream was identified in the L perp − L z space as a localised overdensity in the prograde region (L z > 0) with a high orbital inclination (high L perp ) (see, for example Helmi et al. 1999;Kepley et al. 2007;Koppelman et al. 2019).As a consequence, this space has been applied in literature also to retrieve the origin of globular cluster together with the E −L z space (see for instance Massari et al. 2019).Eccentricity is another important parameter in describing an orbit which has been used for instance to separate the Gaia-Sausage-Enceladus debris from the rest of the stellar halo (e.g.Belokurov et al. 2018;Mackereth et al. 2019;Naidu et al. 2020).In particular the eccentricity − L z plane has been suggested by Lane et al. (2022) as another tool to identify accreted systems (see also Cordoni et al. 2021).Another kinematic space that has been considered in literature is the action space where the horizontal axis is the (normalized) azimuthal action (J φ /J tot ≡ L z /J tot , where J tot = J 2 R + J 2 z + J 2 φ ), while the vertical axis is the (normalized) difference between the vertical and radial actions ((J z − J R )/J tot ).The action space has been suggested as the ideal plane to identify and study possible substructures and debris from accretion events (e.g.Myeong et al. 2018;Vasiliev 2019b;Myeong et al. 2019;Lane et al. 2022) since actions are nearly conserved under the hypothesis that the potential is slowly evolving (Binney & Spergel 1982).Dynamical coherence is the common thread that has guided the interpretation of the kinematic spaces just listed, leading to several tentative reconstructions of the merger history of the Milky Way.For instance, by making use of the globular clusters classification proposed by Massari et al. (2019), Kruijssen et al. (2020) inferred the stellar masses and accretion redshifts of five satellite galaxies accreted by our Galaxy over time, namely Kraken (Kruijssen et al. 2019), the Helmi stream (Helmi et al. 1999), Gaia-Sausage Enceladus (Nissen & Schuster 2010;Belokurov et al. 2018;Haywood et al. 2018;Helmi et al. 2018), Sequoia (Myeong et al. 2019), and the still accreting Sagittarius galaxy (Ibata et al. 1994).Malhan et al. (2022) recovered other past mergers as Cetus (Newberg et al. 2009), Arjuna/Sequoia/I'itoi (Naidu et al. 2020), LMS-1/Wukong (Yuan et al. 2020), and discovered a possible new merger called Pontus (see also Malhan 2022).The classification by Massari et al. (2019), slightly revised, has been also used by Forbes (2020) who combined it with ages, metallicities and [α/Fe] abundances of globular clusters to reconstruct the properties of their progenitor galaxies, while the classifications by Kruijssen et al. (2019) and Malhan et al. (2022) have been used by Hammer et al. (2023) to estimate their infall time.
However, the assumptions used in the previous cited works -namely the hypothesis that accreted clusters should show a dynamical coherence -has not been proven even if several works have started doubting it.For example, various simulations showed that even a single merger debris span over the large range in the E − L z space (Jean-Baptiste et al. 2017;Grand et al. Pagnini et al.: The distribution of globular clusters in kinematic spaces does not trace the accretion history of the host galaxy 2019; Simpson et al. 2019;Koppelman et al. 2020), where, e.g.chaotic mixing may further affect the coherence of the tidal debris (Vogelsberger et al. 2008;Price-Whelan et al. 2016).Moreover, a single merger could result in a number of features or overdensities in the E − L z (Jean-Baptiste et al. 2017;Grand et al. 2019;Koppelman et al. 2020); finally, both E and L z are not conserved quantities in case of evolving and non-axisymmetric galaxy, because of the substantial mass growth of the galaxy over time (Panithanpaisal et al. 2021) and pericentric passages of massive satellites (Erkal et al. 2021;Conroy et al. 2021) perturbing both accreted and in-situ components of the halo.
In this paper, we thus wish to test the hypothesis of dynamical coherence in the E − L z and the other above cited spaces, by analysing N-body simulations of the accretion of satellite galaxies, and their globular cluster systems, onto a Milky Waytype galaxy, containing its own system of clusters.We focus on mergers with mass ratios of about 1:10, because Milky Waytype galaxies are expected to have experienced mergers of similar mass (relative to the Milky Way at the time of the accretion, see for example Deason et al. 2016;Renaud et al. 2021;Khoperskov et al. 2022a) and also because this is roughly the mass ratio estimated for the Gaia Sausage Enceladus accretion event (see Helmi et al. 2018, but also Kruijssen et al. (2020) for a lower estimate).It is worth mentioning that some studies considered even higher Gaia Sausage Enceladus masses (e.g. a 1:5 mass ratio in Naidu et al. (2020)).As detailed in the next sections, our analysis demonstrates that, for such mass ratios, accreted globular clusters do not show the claimed grouping in energy-angular momentum space nor in the other kinematic spaces, in agreement with the results of numerical simulations modelling the population of field stars (see Jean-Baptiste et al. 2017;Amarante et al. 2022;Khoperskov et al. 2022b,a).This finding has strong implications for the reconstruction of their kinematic-based galaxy membership, and questions the accretion history of our Galaxy, as traced so far in the literature (Massari et al. 2019;Myeong et al. 2019;Forbes 2020;Kruijssen et al. 2020;Malhan et al. 2022).The paper is organised as follows.In Sect.2, we present the numerical method, and the characteristics of the simulations.In Sect.3, we present the results of this study.We focus our analysis, firstly, on the spatial distribution of accreted and insitu globular clusters (Sect.3.1), secondly on their distribution in the E − L z plane (Sect.3.2), and we generalise our results showing also how GCs redistribute in other kinematic spaces, such as the L perp − L z and the eccentricity − L z planes, as well as the action space (see Sec. 3.3).We finally apply a Gaussian Mixture Model (GMM) to the outcomes to check whether such a procedure is able to retrieve the real accretion history of the simulated galaxies (Sec.3.4).Motivated by our findings, we reconsider the classification of Galactic globular clusters made on the basis of their distribution in the E − L z plane (Massari et al. 2019), showing that this classification cannot be supported by the age-metallicity relation described by these clusters (see Sect. 4).Finally, we present our conclusions (Sect.5).

Main simulations
In this paper we analyse 14 dissipationless, high-resolution, Nbody simulations of a Milky Way-type galaxy accreting one or two satellites.These simulations are similar to those presented in Jean-Baptiste et al. (2017).As comprehensively explained in that paper, in these simulations both the satellite(s) and the Milky Way galaxy can react to the interaction, experiencing kinematical heating, tidal effects and dynamical friction.Each satellite has a mass which is 1:10 of the mass of the Milky Way-like galaxy and this is motivated by the fact that mergers with these mass ratios are inevitable for Milky Way-type galaxies in ΛCDM simulations (Read et al. (2008); Stewart et al. (2008);De Lucia & Helmi (2008); Deason et al. (2016)).Both the main galaxy and the satellite(s) are embedded in a dark matter halo and contain a thin, an intermediate and a thick stellar disc -mimicking the Galactic thin disc, the young thick disc and the old thick disc, respectively (see Haywood et al. (2013); Di Matteo (2016)).A population of globular clusters initially redistributed in a disc is also taken into account and it is represented as point masses 2 .The total number of particles used in these simulations is N tot = 27 500 110, for the case of a single accretion, and N tot = 30 000 120, when two accretions are modeled.The total number of stellar (thin, intermediate, thick disc) particles in the main galaxy is 20 000 000, the number of globular cluster particles is 100 and the number of dark matter particles is 5 000 000.The satellite galaxies are rescaled versions of the main galaxy, with masses and total number of particles 10 times smaller and sizes reduced by a factor √ 10.All these values are listed in Table 1.The discs are modeled with Miyamoto-Nagai density distributions (Miyamoto & Nagai (1975)) of the form: where masses (M * ), characteristic lengths (a * ) and heights (h * ) vary for the thin, intermediate and thick disc populations (see Table 1); the system of disc globular clusters has scale length and scale height equal to those of the thick disc and the dark matter halo is modeled as a Plummer sphere: (1) 2 Resolving both the simulated galaxies and the globular clusters as N-body systems is currently not possible due to large range of spatial scales that should be modelled -from sub-parsec resolution to hundreds of kpc.In Renaud & Gieles (2015) a numerical method to describe the dynamical evolution of clusters in tidal fields has been carried out.We report to future works that employ the more generic method described in Renaud et al. (2011), which is particularly adapted to study the dynamics of some of our clusters -modelled as N-body systems -in the evolving tidal field generated by some of these accretion events.
A&A Table 2: Initial positions and velocities of the barycentres of the simulated satellite galaxies.The initial inclination of the satellite orbital plane is also given.Distances are given in kpc, masses in units of 2.3 × 10 9 M , velocities in units of 100 km/s and G = 1.Energies are thus given in units of 10 4 km 2 /s 2 and time is in unit of 10 7 years.
where the characteristic mass (M dm ) and radius (a dm ) are listed in Table 1.The choice of using a core dark matter halo for both the Milky Way-type galaxy and the satellite(s) comes from a number of observational evidence that seem to be more consistent with a dark halo profile with a nearly flat density core (Flores & Primack (1994) 2008)).The final N-body models of the main disc galaxy and of the satellite are then built by using the iterative method described in Rodionov et al. (2009).Note that we do not include any significant spheroidal bulge in our model, since the contribution of a classical bulge to the total stellar mass of the Milky Way has been shown to be limited to few percent (see, for example Shen et al. 2010;Kunder et al. 2012;Di Matteo et al. 2015;Gómez et al. 2018), nor any population of halo clusters before the interaction(s) (as previously stated, GCs are all initially confined in the disc).The inner regions of the Galaxy are indeed dominated by a stellar bar, whose impact on the kinematics and distribution of globular clusters has recently started to be explored (see for instance Pérez-Villegas et al. 2020).
Once the N-body systems which represent the main galaxy and the satellite are generated, we translate the barycentre of the satellite galaxy to the (x sat , y sat , z sat ) positions and assign it velocities (v x,sat , v y,sat , v z,sat ) in order to have a parabolic orbit in the case of a 1x(1:10) interaction, with the satellite initially at a distance of 100 kpc from the Milky Way-type galaxy.For each 1x(1:10) interaction, we also vary the initial orientation of the satellite orbital plane (Φ orb ), by rotating it about the y-axis, in such a way to span a range of possible orientations, from prograde orbits (Φ orb = 0 • , 30 • , 60 • , with Φ orb = 0 • being the planar, prograde orbit) to a polar orbit (Φ orb = 90 • ) to retrograde orbits (Φ orb = 120 • , 150 • , 180 • , with Φ orb = 180 • being the planar, retrograde orbit).For the 2x(1:10) interactions, we use as initial orbital conditions for the two satellites a combination of those adopted for single 1x(1:10) mergers.All these initial values, for the 14 simulations (seven 1x(1:10) and seven 2x(1:10)) are reported in Table 2.In this table, each simulation has been given an identifier of the form "MWsat_nN_Φα 1 (−α 2 )", with N being the number of satellites in the simulation and α 1 (α 2 ) the angles of the satellite(s) orbital plane(s) relative to the x − y plane.As for the initial inclination of the satellite disc, each satellite has an internal angular momentum whose initial orientation has been assigned parallel to the z-axis of the host.

Additional simulations
A couple of simulations with mass ratio 1:10 (namely the simulations with ID MWsat_n1_Φ60 and MWsat_n2_Φ30-150) have been rerun in a static Galactic mass distribution, artificially imposing that this latter is not influenced by the accretion of the satellite(s).This experiment allows us to study how the results of the present work would change if dynamical friction exerted by the Milky Way-like galaxy on the satellite and on the clusters was neglected (see Appendix A).Finally, the 1:10 mass ratio simulations are complemented by a set of 7 simulations with mass ratio 1:100, whose initial conditions and analysis are presented in Appendix B, and which are aimed to show the distribution in kinematic spaces of clusters from low-mass galaxies.
All simulations have been run by making use of the TreeSPH code described in Semelin & Combes (2002).Gravitational forces are calculated using a tolerance parameter3 θ = 0.7 and include terms up to the quadrupole order in the multiple expansion.A Plummer potential is used to soften gravitational forces, with constant softening lengths for different species of particles.In all the simulations described here, we adopt = 50 pc.The equations of motion are integrated using a leapfrog algorithm with a fixed time step of ∆t = 2.5 × 10 5 yr.
In this work, we make use of the following set of units: distances are given in kpc, masses in units of 2.3×10 9 M , velocities in units of 100 km/s and G = 1.Energies are thus given in units of 10 4 km 2 /s 2 and time is in unit of 10 7 years.With this choice of units, the stellar mass of the Milky Way-type galaxy, at the beginning of the simulation, is 8.4 × 10 10 M .

Results
In this section, we present and discuss the effect of one or two 1:10 mass ratio accretions on the Milky Way galaxy, focusing our attention on the resulting globular cluster system.We will concentrate our analysis to understand: how accreted and in-situ clusters distribute themselves spatially in the final galaxy; what information about the in-situ or accreted nature of clusters can be extracted from the so-called integrals-of-motion spaces; whether we can use the kinematic information contained in the Galactic globular cluster system today to trace the accretion history of the Galaxy.We show and discuss, firstly, the spatial distribution of accreted and in-situ globular clusters (Sect.3.1), secondly their distribution in the E − L z plane (Sect.3.2), and finally we generalise our results showing also how GCs redistribute in other kinematic spaces and applying a Gaussian Mixture Model to the outcomes to check whether such an approach is able to provide robust information about the accretion history experienced by the remnant galaxy (Sec.3.3).In the following of this analysis all quantities are evaluated in a reference frame whose origin is at the centre of the Milky Way-type galaxy with the spin of the Milky Way-type galaxy being oriented as the z-axis and positive.The centre is evaluated, at each snapshot of the simulations, as the density centre, following the method described in Casertano & Hut (1985).

Spatial distribution of accreted and in-situ clusters
Fig. 1 shows the globular clusters projections in the (x, y) and (x, z) planes, from one of the 7 single-accretion simulations (Φ orb = 60 • , i.e. simulation ID = MWsat_n1_Φ60), at different times: the initial (T = 0 Gyr), two intermediate times (T = 1 Gyr, 2.5 Gyr), and at the final time (T= 5 Gyr).In all the plots, the in-situ globular cluster system, i.e. originally in the Milky Way-like galaxy, is represented by grey circles, and the globular cluster system initially linked to the satellite by magenta circles.The distribution of field (in-situ and accreted) stars is also shown in the background.
At the beginning of the simulation, the Galactic and satellite globular cluster systems have both an axisymmetric disc distribution, as the corresponding field stellar particles.This initial axisymmetric configuration of the clusters and field stars is rapidly disturbed by the first satellite passage (occurring at T 0.3 Gyr, see Fig. 2): at time T = 1.0 Gyr, the satellite has disturbed the outer parts of the disc of the Milky Way-type galaxy, which has developed an extended spiral-like structure.In the meantime, the outer parts of the satellite are tidally disrupted by the main galaxy and also one of the clusters originally associated with the satellite has escaped to join the outer parts of the massive galaxy.We can also appreciate that the vertical structure of the disc of the massive galaxy and the corresponding cluster system start to be kinematically "heated" at this time: seen from the edge, the disc has thickened, it shows tidal plumes in the outer parts, partly dragged by tidal pulls of the satellite, partly made of in-situ stars, disturbed by the interaction.This kinematic heating is the result of the redistribution of a fraction of the orbital energy of the satellite into internal energy of the stars in the merger remnant (see, for example Quinn et al. 1993;Walker et al. 1996;Villalobos & Helmi 2008;Zolotov et al. 2009;Purcell et al. 2010;Di Matteo et al. 2011;Qu et al. 2011;McCarthy et al. 2012;Cooper et al. 2015;Jean-Baptiste et al. 2017).At T = 2.5 Gyr the merging process is toward the end and this is also visible in the plots: the stellar distribution of the satellite is no longer clearly sepa-Fig.2: Top panel: Time evolution of the distances of the satellite (black line) and its globular clusters (coloured lines) from the Milky Way-type galaxy, for the simulation MWsat_n1_Φ60.Each globular cluster is colour-coded according to its escape time from the progenitor satellite.Bottom panel: Time evolution of the apocentres of the orbits of all in-situ clusters in the Milky Way-type galaxy, for the simulation MWsat_n1_Φ60.Thick and coloured lines indicate clusters whose orbital apocentres, at the final time of the simulation, are at least 1.5 larger than their corresponding initial values; these are colour-coded according to the ratio of the final apocentre over the initial one.rated from that of the main galaxy, and also most of its globular clusters are now part of the bulk of the Milky Way population.However, we can still appreciate a series of tidal plumes and two globular clusters originating from the satellite that have not yet been trapped by the Milky Way-type galaxy.At the final time of the simulation (T = 5 Gyr), we can see that the majority of the clusters (in-situ and accreted) are located in the inner densest regions of the final galaxy and if it was not for the different colors, they would not be distinguishable from each other.The decay of part of the accreted clusters in the innermost regions of the remnant can be understood because these are the clusters that remain gravitationally bound to the satellite until the final phases of the merger.It is primarily the dynamical friction experienced by the satellite to which these clusters belong that make them decay towards the centre of the main galaxy.It is only when these clusters become gravitationally unbound to the satellite, that their motion from their progenitor galaxy decouples.This can be appreciated in Fig. 2, where we show -for the case of the simulation MWsat_n1_Φ60 -the temporal evolution of the distance of the satellite to the main galaxy, together with the corresponding evolution of all the clusters initially bound to the satellite.We can see that soon after the first pericentre passage, the satellite loses one of its clusters, which is trapped on an orbit with apocentre at about 40 kpc, while the satellite itself moves to its first apocentre, at about a distance of 95 kpc from the main galaxy centre.A second cluster is lost by the satellite at its next pericentre passage (T = 2 Gyr).The cluster lost at this time is ejected on a high energy orbit, which has an apocentre at more than 140 kpc from the centre of the Milky Way-type galaxy.Over the entire duration of the merging process, this loss of the satellite clusters from their progenitor leads to their redistribution over a large ranges of distances (pericentres and apocentres) from the merger remnant at the final time.While clusters associated to the satellite galaxy are lost at different times of the interaction, and hence contribute to populate both the outer and inner regions of the remnant galaxy, clusters initially in the Milky Way-type galaxy -which, we remind the reader, are initially distributed in a disc-like configuration -have their orbits perturbed by the interaction.In Fig. 2 (bottom panel), we show indeed the evolution with time of the orbital apocentres -computed as the maxima of the 3D distance -of all the 100 clusters in the Milky Way-type galaxy.GCs whose orbital apocentres, at the end of the simulation, are at least 1.5 larger than their corresponding initial values are shown as thick and colour coded lines while the others as grey lines.We can see from this plot that 14% of the clusters have their orbital apocentre increased of a factor 1.5 at least, and for some of them this heating coincides with the last phases of accretion of the satellite (around T = 2 Gyr).The number of kinematically heated clusters increases to 27% when we consider GCs that have their final orbital apocentre increased of at least a factor of 1.2 with respect to the initial one.

Energy -angular momentum space
All the results presented in this Section concern the distribution of globular clusters in the energy -angular momentum space (E− L z ) in the case of a Milky Way-type galaxy accreting one or two satellites, since this space has been proposed as the natural one where to look for the signatures of past accretion events (Helmi & de Zeeuw 2000).

Accreted clusters
Let us start with the case of a single accretion.In Fig. 1 (right column), we show the distribution in E − L z space of the accreted clusters in the case of the Φ orb = 60 • simulation (ID=MWsat_n1_Φ60), at different times, T= 0, 1, 2.5, 5 Gyr (these times correspond to those for which the x − y and x − z maps are shown in the left and middle columns of the same figure).In these E − L z plots, we also report, for each cluster, the escape time from its progenitor satellite.This escape time has been estimated by means of a simple spatial criterion, i.e. it is defined as the time when the distance between the GC and the centre of mass of the parent satellite is larger than 15 kpc.For comparison, the distribution in E − L z space of field stars of the satellite is also plotted in the background.At the beginning of the interaction, the distribution of satellite stars and globular clusters is clumped in the E − L z space and characterized by high energy.This is understandable because the satellite is, at the beginning of the simulation, a gravitationally bound system.At T = 1 Gyr the satellite is close to an apocentre passage (see Fig. 2), resulting in a broad extent in L z groups: the first consisting of GCs with t esc < 1 Gyr, the second with 1.8 Gyr < t esc < 3.6 Gyr and the last with t esc = 5 Gyr.The first two groups are associated with globular clusters lost at the first and subsequent pericentric passages, respectively.Globular clusters with t esc = 5 Gyr are those which did not escape from the satellite before the end of the merging process.From the color-coding of globular clusters in Fig. 3, it is possible also to notice that the clusters that escape the earliest from the progenitor satellite are those which are initially less bound to their progenitor.They are indeed the clusters with the highest values of E int , where E int is the sum of the kinetic energy of clusters in the satellite and of their potential energies, both estimated relative to the satellite centre.As expected, globular clusters that are more tightly bound to the satellite tend to escape later and have a lower final energy relative to the Milky Way-type galaxy reference system.However, from Fig. 1 (bottom, right panel) we can see that this trend shows a number of exceptions: for instance, the cluster lost at 2.59 Gyr, at the end of the simulation has lower energy (E ≈ −9.5) with respect to the clusters lost later at 3.05 Gyr (E ≈ −8.0).The globular cluster lost at 0.58 Gyr also ends up at higher energy than the cluster lost at 0.42 Gyr.This happens because after leaving their parent satellite, the energy and angular momentum of globular clusters may change due to changes in the gravitational potential induced by the final phases of accretion of the same satellite.
When two satellites are accreted, the interpretation of the energy-angular momentum space becomes even more tricky (see also Trelles et al. 2022).Figure 4 (left, middle panel) shows the distribution in the E − L z space of globular clusters originally belonging to the two satellites (magenta and white colors, respectively) for one of the two accretion simulations (MWsat_n2_Φ30 − 150).The background density maps show the fractional contribution of stars from satellite 1 (left panel) and satellite 2 (middle panel) relative to the totality of stars.As we can see from the figures, the globular cluster populations of the two satellites mix and overlap with each other.Furthermore, the cluster distribution does not necessarily coincide with the distribution of stars of the same satellite.For example, nearly half of the globular cluster population belonging to satellite 1 is redistributed in a region of the E − L z space dominated by stars of satellite 2 (see middle panel).The contribution of stars belonging to satellite 1 is still significant in this region, but not exclusive.Therefore, the stellar halo substructures and accreted GCs located in the same E − L z regions can not be directly associated with the same dwarf galaxy progenitor accreted in the past: regions of the E − L z dominated by stellar populations of a given progenitor can indeed contain a significant fraction of clusters (and stars) originating from a different progenitor.
Simulation ID   Table 3: Mean and standard deviation of the initial and final L z and E for the populations of accreted GCs within the different single accretion simulations.

In-situ clusters
The right panel of Fig. 4 shows the same distribution in the E − L z space of globular clusters originally belonging to the two satellites with the addition of GCs belonging to the Milky Way-type galaxy (grey circles).The density map illustrates the fractional contribution of stars from the two satellites relative to the totality of stars.This panel allows us to show an additional result: part of the in-situ cluster population overlaps with the accreted clusters in E − L z space.This in-situ population is made of disc clusters which have been perturbed enough by the interaction to be "pushed" into the halo (see also Fig. 2, bottom panel), following the same dynamical mechanism extensively discussed for in-situ field stars (see, for example, Zolotov et al. 2009;Purcell et al. 2010;Qu et al. 2011;Jean-Baptiste et al. 2017;Khoperskov et al. 2022b).We emphasize that, by construction, our modelled Milky Way-type galaxy does not contain, before the interaction(s), any population of halo clusters (they are all initially confined in the disc).This implies that the overlap of in-situ GCs with satellites GCs could be even more significant, if part of the in-situ population had halo-like kinematics already before the accretion(s).We further address this issue in Sect.3.5.Interestingly, if we compare the distribution of globular clusters with the distribution of stars in the right panel of Fig. 4, we can notice that a part of the in-situ GC population ends up in regions of E − L z space dominated by stars accreted from the two satellites.All this has as a consequence that, in the interpretation of the E − L z space, we cannot simply associate stars and clusters to the same origin (i.e. the same progenitor galaxy) simply because they are found at similar values of E and L z : clusters from a satellite can be found in regions where the stellar density distribution is dominated by another satellite (middle panel, Fig. 4), and clusters originally in the main galaxy can be found in regions of the E − L z space dominated by accreted stars (right panel Fig. 4).
About the possibility to separate in-situ from accreted clusters, it is interesting to cite the numerical work by Callingham et al. (2022), who constructed mock catalogues of accreted and in-situ clusters from the AURIGA simulations, reporting that "the total in-situ population is, in general, well recovered with very high purity.Rarely does the methodology misidentify an accreted GC as an in-situ one in our mock tests, with a median purity of 98%." This conclusion is clearly in contradiction with our results, but can be explained by the fact that in Callingham et al. ( 2022) work, while their accreted GC populations are drawn from the AURIGA simulations, the in-situ GCs is added a posteriori (that is, it is not extracted self-consistently from the AURIGA simulations, as for the accreted population).In their work, the in-situ GCs are constructed to be either in the bulge or in the disc, that is they do not allow the possibility that part of the in-situ GC population can have a halo-like kinematics (as it is the case, as we have shown, when part of a pre-existing disc GCs population is heated by a merger).The reason why they can separate so well in-situ disc GCs from accreted GCs in kinematic spaces (see for example the E − L z diagram shown in their Fig 3) is due to the way they construct the in-situ GC population, and it is not the result of a self-consistent evolution of the in-situ population itself during the merger(s).It would have been interesting, and more realistic, to draw also the in-situ GC population from their AURIGA simulations, by randomly extracting stars formed in the AURIGA MW-like progenitors.If this was the case, we argue that the separation between in-situ and accreted components would have been more challenging also in their models.

Overlap in the E − L z space
Figure 5 shows the final globular cluster distribution in the E −L z space for the whole set of 1x(1:10) accretion simulations.The distribution of the in-situ globular clusters from all the simulations at T = 5 Gyr has been stacked and is shown by grey circles.Accreted clusters are colour-coded circles according to the different Φ orb and linked by a dashed line.This figure summarise the results presented so far: regardless of the initial inclination of the satellite orbital plane, the overlap between accreted and in-situ globular clusters is clear and becomes critical in the most gravitationally bound regions.This may be partly caused by the narrowing of the L z range at high binding energies, but also may reflect the tendency of high-mass-ratio mergers to radialize their orbits (e.g.Amorisco 2017; Naidu et al. 2020;Vasiliev et al. 2022).The figure also shows that, in the case of a single accretion, the angular momenta of clusters at high energies (E −10 in our units) correlate with the inclination of the orbital plane of the infalling satellite at early times: the higher the orbital inclination Φ orb of the parent satellite, the less prograde (lower L z ) the orbits of the clusters are.We will see in the following that this trend is less prominent in the case of two accretions.
Fig. 6 shows the distribution in the energy-angular momentum plane of the whole set of globular clusters in the 2x(1:10) accretions simulations.The distribution of the in-situ population of globular clusters has been stacked and is shown by grey circles, while accreted GCs (belonging to satellite 1: left panel, to satellite 2: right panel) are shown as colour-coded circles according to the different Φ orb and linked by a dashed line.As we can see in Fig. 6, in-situ and accreted globular clusters overlap almost everywhere in the E − L z space.Only for E −6 and/or L z −3 in principle we can distinguish the two groups.Even in that case, however, it is not possible to identify from which of the two satellites the clusters originate because, as mentioned above, these clusters appear mixed, that is clusters originating from different satellites can end up having similar energies and angular momenta.Moreover, the trend between Φ orb and the final L z found for high-energy clusters in the case of a single accretion is clearly washed out in the case of two mergers.Looking, for example, at the case of satellites accreted with an initial orbital inclination of Φ orb = 180 • (as for the simulations MWsat_n2_Φ0-180 and MWsat_n2_Φ180-90), the final distribution of their clusters in the E − L z plane is different for the two 2x(1:10) simulations (compare the top-left and top-right panels of Fig. 6).This distribution appears also different from the one generated by clusters whose progenitor is on a Φ orb = 180 • inclination orbit, in the case of a single accretion (Fig. 5).This indicates that -unless there are strong reasons to think that the Galaxy experienced only a main massive accretion (and not a few) -we cannot infer from the current distribution of its accreted globular cluster population the inclination of the progenitor satellite.
From this analysis, we conclude that the only merger events that in principle could be distinguished in the E − L z plane are: (1) those that are still occurring and for which the population of GCs that the satellite is bringing with itself populates a rather delimited region at very high energies (as it is the case for the Sagittarius dwarf galaxy and its globular cluster system, see Bellazzini et al. 2020); (2) those that occurred in the past but which were characterised by a smaller mass ratio, typically 1:100 and lower (as discussed in Jean-Baptiste et al. 2017).In this case, stellar debris tend to be distributed in the high-energy regions of the E − L z plane (see Sect. 4.4 in Jean-Baptiste et al. 2017, and also Pfeffer et al. (2020)), and a similar behaviour is found also for the associated globular cluster population.We refer the reader to Appendix B for the distribution in the E − L z space obtained when considering the Milky Way accreting four 1:100 mass ratio satellites.Let us now examine how the populations of globular clusters in our simulations redistribute in other kinematic spaces.In addition to the E − L z space, we have analysed the L perp − L z space where L perp is the projection of the total angular momentum onto the Galactic plane and is defined as: L perp = L 2 x + L 2 y (L x , L y being respectively the x and y component of the angular momentum space in a reference frame with the Galactic disc in the x − y plane).Eccentricity (e), defined as e = R apo −R peri R apo +R peri (R apo and R peri being respectively the apocentre and the pericentre of the orbit) is another important parameter in describing an orbit.Here we use the eccentricity by combining it with L z as suggested in Lane et al. (2022) (see also Cordoni et al. 2021).The last kinematic space we analysed is the action space where the horizontal axis is the (normalized) azimuthal action (J φ /J tot ≡ L z /J tot , where , while the vertical axis is the (normalized) difference between the vertical and radial actions ((J z − J R )/J tot ).The right and left points of the space (|L z | = J tot , see bottom right panel in Fig. 7) are in-plane prograde and in-plane retrograde orbits respectively.The top point (J z = J tot ) is a circular polar orbit, and the bottom point (J R = J tot ) is a radial orbit.The bottom-right and bottom-left edges are prograde and retrograde in-plane orbits.The top-right and top-left edges are prograde and retrograde circular orbits.
Figures 7, 8 show the final globular clusters distribution in the four kinematic spaces described above i.e. the E − L z , L perp − L z , eccentricity − L z and action space for the two previously examined simulations with one or two accretions respectively (MWsat_n1_Φ60 and MWsat_n2_Φ30-150).The population of globular clusters originally belonging to the MW-type galaxy is represented by grey circles while the accreted one by magenta (and orange) circles for satellite 1 (and satellite 2).Actions are computed using the Stäckel fudge approach (Binney 2012), as implemented in the AGAMA code (Vasiliev 2019a).The addition of L perp (see top right panel in Fig. 7) to the previous analysis of the E − L z space only confirms the conclusions drawn in the previous paragraph: at the end of the simulation, overall the superposition between the in-situ and accreted population is considerable everywhere.Only GCs lost by the satellite at the first pericentric passage (t esc 0.5) which are therefore at high energies (E −5) are in principle distinguishable from the rest.Furthermore, the accretion event generates a population of in-situ halo clusters (consisting of disk GCs heated by the interaction), whose distribution extends towards high L perp .From the eccentricity − L z and action spaces (see bottom panels in Fig. 7), it is not possible to make any difference between accreted or insitu globular clusters.In the action space we can notice a higher density in the right-hand corner related to the more prograde orbits where in-situ and accreted GCs definitely overlap.Overall clusters from the satellite are preferentially located towards the bottom-right edge, i.e. towards prograde in-plane orbits, but are not found in distinct groups from in-situ clusters.The strength of the action space should lie in the fact that distinct types of orbits occupy different loci in the space (see Lane et al. 2022), but this is not sufficient to reconstruct the origin of the globular clusters (i.e.accreted or in-situ), since during the merger process they are scattered and mixed over most of the space.
As we have already mentioned, when we consider the merging with two satellites (see Fig. 8), the interpretation of kinematic spaces obviously does not improve.We can confirm that overall GCs that end at higher energy in the E − L z plane, are accreted.In this specific simulation (MWsat_n2_Φ30-150) they also stand out as the most retrograde (see top and bottom left panels of Fig. 8).Despite this, it is still hardly feasible to distinguish clusters coming from the first or second satellite and it is evident that also in-situ clusters originally in the disc can be found at these energy levels.If the points were not coloured, it would not be possible to assign a membership to the different GCs, on the basis of their location in the E−L z space.Other kinematic spaces added in the analysis seem not to add any relevant information for this purpose.In the L perp −L z and eccentricity−L z spaces, in fact, the only clusters that are detached from the rest of the mixed population of accreted and in-situ clusters are those with a more retrograde orbit comprising 2 MW GCs, two clusters originating from satellite 2 and one GC from satellite 1.In the action space, again, the clusters are scattered and mixed over most of the plane; the accreted clusters are predominantly in the region characterised by radial in-plane orbits, but we cannot identify groups with the same galactic membership: the small clumps detectable by eye are all composed of GCs of different origins.In this respect, it is worth recalling the work by Lane et al. (2022) which extensively investigates the best kinematic spaces to separate radially anisotropic from isotropic halo populations, concluding that, to this aim, the "action diamond" space is superior to other spaces used in the literature.The problem is that clusters of different nature (accreted or in-situ, and/or be- longing to different satellite progenitors) can end up having similar orbital anisotropy (see, for example Figs 7 and 8), thus a kinematic space can be very good at separating stars (or clusters) with specific orbital properties, but without being able, however, to assess their origin.This is the limitation of the analysis presented in the literature so far: to associate a specific region of a kinematic space, and hence specific orbital characteristics, with a specific accretion origin.
In Fig. 9 and 10 we have stacked together the GCs distributions in L perp − L z , eccentricity − L z and action − space for the whole set of 1x(1:10) and 2x(1:10) simulations respectively as in Fig. 5 and 6.All in-situ globular clusters are shown with gray circles while accreted GCs are colour-coded according to the initial inclination of the satellite of membership with respect to the Milky Way-like reference frame, Φ orb .These figures allow us to generalise the conclusions just drawn for the two examples of single and double accretion, the most important of which concerns the fact that the kinematic spaces considered in addition to the E − L z space do not add any useful insight for discriminating between globular clusters originating from satellites accreted over time by the MW or formed in the MW itself.Regarding the L perp − L z plane, we note that in the case of a single accretion, if we select the region with L perp ≥ 20 and L z 30, we will certainly find satellite clusters, an argument that no longer holds if we analyse simulations with double accretion: here indeed we find several in-situ clusters heated so much kinematically that they end up having L perp > 20.Looking at the eccentricity − L z plane, the only noteworthy point is the fact that with both one and two accreted satellites, at fixed |L z | > 0, we find accreted GCs with more eccentric orbits than those in-situ, having e ≥ 0.8 while for |L z | 0 accreted and in-situ clusters share the same more extended range of eccentricities.As far as action space is concerned, we can say that in general accreted clusters tend to have in-plane prograde or radial orbits.The distribution of insitu globular clusters is denser in the right corner corresponding to prograde in-plane orbits and this makes sense since, by hypothesis, the clusters formed in the MW-type galaxy initially have disc orbits.Despite this, we also find in-situ GCs at the opposite extreme -i.e. with strongly retrograde orbits -and indeed, the most retrograde clusters are actually in-situ in all cases.If we had also considered in-situ GCs with halo-like orbits, this feature would have been even more evident (see Sect. 3.5).

Does a clustering analysis allow to recover the history of accretion?
To make our analysis more quantitative, we exploited Gaussian Mixture Models both as a clustering algorithm but mainly for its fundamental purpose, which is to model the overall distribution of the input data.The purpose of this analysis was to investigate whether the algorithm can grasp any properties of the final distributions of GCs in kinematic spaces not visible by eye and in a more quantitative manner.Lately, several works have started to use this type of tool both to distinguish between accreted and in-situ stars/globular clusters and to retrieve the number of accretions undergone by the Milky Way (see for instance Donlon ).We then proceeded as follows.For each simulation we considered a six-dimensional space consisting of the globular clusters' kinematic quantities previously examined, namely: the z component of the angular momentum L z , the total orbital energy E, the projection of the angular momentum onto the galactic plane L perp , the eccentricity e, the normalized azimuthal action L z /J tot , and the difference between the vertical and radial actions (J z − J R )/J tot .We tried to fit this with a two/three-component GMM (depending on the number of accreted satellites considered in addition to the insitu component) viewed as a clustering model, but the results were not particularly useful as the real memberships were not retrieved.Therefore, we fitted several models with increasing number of components and determined the optimal number of components for our dataset by minimising the Bayesian information criterion (BIC, Schwarz (1978)).true labels given by the simulations (MW, satellite or sat1 and sat2) while each column represents the predicted group identified by the model.Values are normalised to the total number of GCs in each true class.Interestingly, for both types of simulation, there is a predicted class that includes both the majority of accreted and in-situ clusters.This class (group 3 in Fig. 7 and group 7 in Fig. 8) covers the region of the E − L z space where we have the largest overlap between GCs from different galactic progenitors, namely the region with −17 E −12.Overall, the groups that do not contain accreted clusters are those that trace highly prograde disc orbits in E − L z space, while the rest of the groups into which the MW is divided also contain a high fraction of satellite(s) clusters.The GMM, as it is designed, prac-tically cuts the different kinematic spaces (with the exception of the action space) into separate slices that overlap only slightly.This behaviour therefore tends to support the interpretation suggested so far based on grouping in kinematic spaces but fails in grasping the underlying physical process.We acknowledge the good prospects in using this type of tool to describe, in a quantitative manner, the distribution of stars and clusters in kinematic spaces, but also suggest caution in interpreting the results, as the components retrieved by the algorithm are not directly related to the accretion events experienced by our Galaxy: the number of groups does not trace the number of satellites accreted over time, and groups are in general made of a mixture of in-situ and accreted populations.

Adding an in-situ halo population
As we have stressed several times in previous sections, by construction, our modelled Milky Way-type galaxy does not contain any population of halo clusters before the interaction(s), since insitu GCs are all initially confined in a Miyamoto-Nagai disc (see Sec. 2).In literature, the presence of an in-situ halo is still being debated.For instance the only in-situ population that Haywood et al. (2018);Di Matteo et al. (2019) find in great proportions in Gaia DR2 and APOGEE data, is the thick disc i.e. the early disc of the Galaxy heated to hot kinematics.On the other hand, Belokurov & Kravtsov (2022) support a period of chaotic predisc evolution when stars are born in dense clumps scattered on all kinds of orbits and thus populating a hot halo.If the latter scenario were the case (both for stars and globular clusters), this would imply that the overlap between in-situ and satellites GCs in kinematic spaces could be even more significant.To test this claim with our simulations, we randomly selected a group of 20 particles from the dark matter halo of the MW-type galaxy, initially modelled with a Plummer distribution, and assumed these to be globular clusters formed in-situ, originally with halo kinematics.This was feasible since, as shown in Table 1, the mass of dark matter particles in our simulations is 7.4 × 10 4 M and so consistent with the mass range of Milky Way GCs (see, for example, Baumgardt et al. 2019).
Figure 13 shows the final globular clusters distribution in the four kinematic spaces analysed i.e. the E − L z , L perp − L z , eccentricity − L z and action space for the simulation MWsat_n1_Φ60 -as Figure 7 -with the addition of the mock halo in-situ GCs as empty circles.As expected, the population of in-situ GCs that initially has a halo kinematics, maintains the same type of kinematics: in E − L z space they end up in a rather high energy region since the dynamical friction on the individual clusters is not strong enough to allow them to lose energy and we also find many GCs with negative L z .As a result, the different kinematic spaces -in particular the E − L z space -become even more indecipherable and at this point we are no longer sure whether even high-energy or highly retrograde clusters are accreted.In this scenario therefore, with only kinematic information, not even GCs coming from less massive satellites -as for instance with a 1:100 mass ratio (see App. B) -would be distinguishable from in-situ GCs.

Discussion: a new look at Galactic globular clusters
The results presented in the previous section show that we do not expect globular clusters accreted with their own parent satellite to clump in specific regions of the E − L z plane nor in the other analysed kinematic spaces, unless there are reasons to assume that the Galaxy has not experienced any significant (i.e.mass ratio ∼ 1 : 10) merger -assumption which would contradict the current estimates of the mass ratio of the Gaia Sausage-Enceladus merger (see, for example, Helmi et al. 2018).Moreover, if the Galaxy has experienced more than one significant merger during its evolution (see, for example, Kruijssen et al. 2020), clusters associated to these mergers can overlap and mix in E − L z space, particularly in the most gravitationally bound regions of the diagram (i.e.low energies).
A suite of mergers would also muddle the traces left by previous accretion events in this space.Finally, a region of the E − L z space dominated by stars associated to a progenitor satellite is not necessarily dominated by clusters associated to the same accretion event.However, the current reconstruction of the merger tree of the Galaxy is exactly based on these assumptions: (1) dynamical coherence of globular clusters in the kinematic spaces, (2) negligible overlap of clusters originating in different satellite galaxies in all the kinematic spaces, (3) no kinematic heating of the in-situ population, (4) straightforward correspondance between field stars and clusters in the E − L z space, meaning that the regions of this space where field stars associated to Sequoia, Gaia Sausage-Enceladus, Helmi stream -and other possible accretions -have been identified are also the regions used as boundaries to assign clusters to each of these progenitors.
Our results show that there is no physical reason to proceed in this way.If we want to look for the remnants of these accre-tion events, we cannot require that the associated clusters satisfy any dynamical coherence in kinematic spaces.In this respect, let us recall the results by Kruijssen et al. (2019), who identified traces of an ancient accretion event in the Galaxy, called Kraken, studying the system of Galactic globular clusters.More specifically, Kruijssen et al. (2019) based their conclusions on the study of the age-metallicity relation of these clusters which was compared to that of a suite of cosmological zoom-in simulations of Milky Way-mass galaxies from the E-MOSAICS project.This association of the possible clusters associated to Kraken was then questioned by Massari et al. (2019), who noted that most of the clusters reported in Kruijssen et al. (2019) as possible mem-bers of Kraken were not dynamically coherent5 , since associated to different regions of the E − L z diagram.They thus reconsidered the initial suggestion of Kraken-like globular clusters made by Kruijssen et al. (2019), proposing a new one where Krakenlike clusters are a group of dynamically-coherent clusters found in the low-energy part of the E − L z plane, and concluded that ".. taking into account the dynamical properties is fundamental to establishing the origin of the different GCs of our Galaxy."Note that it is on the basis of this new classification that Kruijssen et al. (2020) subsequently based the reconstruction of the merger tree of the Milky Way (see also Forbes 2020).Our work urges to reconsider this overall approach, since it shows that the remnants of massive accretion events are not expected to show any global clustering in this space, as already shown to be the case also for field stars (Jean-Baptiste et al. 2017;Amarante et al. 2022;Khoperskov et al. 2022a): it is not because of a lack of dynamical coherence in the E − L z space that a subset of Galactic globular clusters cannot be associated to the same progenitor.If we return to the classification of the Galactic globular clusters made by Massari et al. (2019) on the basis of their kinematics (energies, and angular momenta), we can take a new look at the age-metallicity relation(s) (AMRs) of clusters which in their study are associated to different progenitors.It is important to rediscuss these age-metallicity relations briefly here, because they have been used in the literature to further justify kinematicallybased classifications.For example, Massari et al. (2019) concluded that quite remarkably, "the dynamical identification of associations of GC results in AMRs that are all well-defined and depict different shapes or amplitudes".However, in their work, no actual fit to the data was done (as acknowledged by the authors themselves).It is worth thus reconsidering whether the kinematic classification by Massari et al. (2019) effectively leads to groups whose AMRs have different shapes and amplitudes, especially considering that all age estimates come with associated errors, that each sample has a finite and limited number of clusters for which ages are available (typically less than 10) and that different ages and metallicity estimates exist in literature.Do kinematically-based classifications effectively lead to select clusters with different enrichment histories?And how dependent are the retrieved enrichment histories on the chosen set of ages and metallicities?These are fundamental questions, since the apparent difference between AMRs of groups of clusters identified on the basis of their kinematics has been used as an additional probe of the robustness of the classification itself.To rediscuss this issue, we make use of ages and metallicities from literature data, as reported by Marín-Franch et al. (2009) and VandenBerg et al. (2013).We recall the reader some main differences and similarities among these studies.Marín-Franch et al. (2009) measured relative ages of a sample of 64 Galactic globular clusters, observed in the framework of the HST/ACS Survey of Galactic globular clusters.The corresponding ages, and errors, are reported in Table 4 of their paper, for a set of different theoretical isochrones, and metallicities for two abundance scales.We adopt in the following the ages corresponding to the theoretical isochrones of Dotter et al. (2007) using the Zinn & West (1984)  Gaia-Sausage Enceladus, Sagittarius, Helmi Streams, Sequoia, together with two additional groups, the Low-Energy clusters, and High-Energy clusters (see Table 4).Despite the differences in the shape and extension of the age-metallicity relations resulting from these two datasets, we notice that: 1. some of the clusters classified as accreted clusters, either associated to the Low-Energy group or to the Helmi stream -namely E3, NGC 6441, and NGC 6121 -have ages and metallicities compatible with being in-situ disc clusters heated to halo-kinematics (see also Forbes 2020;Horta et al. 2020, for similar conclusions).They are indeed on the old, metal-rich branch of the age-metallicity relation (see for example the left panel of Fig. 14), where most of the disc-bulge clusters lies, but have hotter kinematics than that expected from disc-bulge systems (and this is the reason why they have not been classified by Massari et al. (2019) 2019) identified.Note that we restricted this analysis only to the groups for which the association was considered robust (i.e.Sag, H99, G-E, L-E, but not groups as H99/G-E, H99?, G-E?, etc) excluding Sequoia (Seq) and the high-energy group (H-E) since for them the number of clusters for which ages are provided by Marín-Franch et al. (2009) 2009), we find that the chemical evolutions of the Sag and H99 groups have similar yields (respectively p = 0.72 ± 0.14 and p = 0.79 ± 0.15) and formation times (respectively t f = 1.06 ± 0.07 and t f = 1.10 ± 0.06).The L-E and G-E groups have slightly higher p and t f (respectively p = 0.80 ± 0.01 and p = 0.83 ± 0.06, t f = 1.20 ± 0.01 and t f = 1.18 ± 0.03), however the evolution of G-E is still compatible -within 1σ -with those found for Sag and H99 groups.If we use the AMRs by VandenBerg et al. (2013), the retrieved chemical enrichment histories are all indistinguishable within 1σ (H99: p = 0.71 ± 0.12, t f = 12.94 ± 0.40, Sag: p = 0.68 ± 0.09, t f = 13.20 ± 0.34, L-E: p = 0.64 ± 0.06, t f = 13.40 ± 0.27, G-E: p = 0.75 ± 0.04, t f = 13.48 ± 0.19).Thus, the conclusion that kinematically-based classifications lead to groups of clusters with different chemical enrichment histories is not supported by a robust analysis of the data, and this for two independent datasets of ages and metallicities.
To summarise, we have demonstrated that the assumption of "dynamical coherence" for the interpretation of globular clusters in kinematic spaces is not supported by physical arguments (unless a very specific merger history for our Galaxy is assumed), and indeed we do not find in the observational data any confirmation that merger histories based on this assumption identify clusters with specific age-metallicity relations, and hence star formation histories of their progenitor systems.There is room for a new interpretation of AMR and/or chemical abundance spaces.Dynamical coherence should not be the "a priori" assumption to analyse globular cluster data, as done essentially by all studies published in the literature so far (see, for example, Forbes 2020; Callingham et al. 2022) with the exception of Kruijssen et al. (2019); Horta et al. (2020).Even if uncertainties on ages are still significant, and even if possible overlaps among different evolutionary sequences can be challenging 7 , it is by finding features in 7 See for example the overlap between the in-situ and accreted branches in the AMRs or abundance planes at low metallicities.
the age-metallicity relations or in abundance-spaces -which are conserved quantities through time -that we can hope to solve the question of the accretion history of the Galaxy.

Conclusions
In this paper, we have analysed dissipationless N-body simulations of a Milky Way-type galaxy accreting one or two satellites with mass ratio 1:10, as well as some 1:100.Each galactic system has a population of disc globular clusters represented by point masses.We have analysed this set of simulations to investigate the possibility to make use of kinematics information to find accreted globular clusters in the Galaxy, remnants of past accretion events and which have lost their spatial coherence.In particular, we have examined the energy -angular momentum (E − L z ) space, and shown that: (1) clusters originating from the same progenitor generally do not group together in E − L z space and thus GCs populating a large range of energy and/or angular momentum could have a common origin; (2) if several satellites are accreted, their globular cluster populations can overlap, in particular in the region E − L z plane where the most gravitationally bound clusters are found (i.e.low energies); (3) the clusters distribution does not necessarily match that of the field stars of the same progenitor galaxy: in several cases we find regions of E−L z space populated by clusters originating from one satellite, but dominated by field stars from another satellite.This implies that the correspondence between field stars and associated clusters in the E − L z diagram is not necessarily trivial: if a certain region of this space is dominated by the stellar debris of a satellite, this does not imply that clusters found in the same region all originated from the same satellite; (4) the in-situ population of globular clusters (that in our models is initially confined into a disc) can be heated up by the accretion, and a fraction of it acquires halo-like kinematics, thus becoming indistinguishable -on the basis of the kinematics alone -from the accreted populations.
We find that accreted clusters are confined in a tight range of energies only when they originate in low mass progenitors (mass ratio of 1:100).For such mass ratios, the distribution of angular momenta is still very extended, and these clusters tend to be all in the high energy part of the E − L z diagram, since dynamical friction is not able to bring them to the inner regions of the Milky Way-type galaxy before their parent satellite is severely affected by the gravitational tides exerted by the main galaxy.Because clusters associated to low-mass progenitors are all found at highenergies in our simulations, a significant overlap is found also for these systems.Note that, besides the effects studied in this paper i.e. dynamical friction, perturbations induced by the other ongoing accretions (see also Garrow et al. 2020), other processes can contribute to the non-conservation of energy and angular momenta of globular clusters, as for example the mass growth of the Galaxy with time, which is expected to be significant especially in the first 4 − 5 Gyr of its evolution (Snaith et al. 2014).Interestingly, simulations run in a cosmological context (Khoperskov et al. 2022a,b)   For these reasons, we are confident that these results are robust and would find confirmation also in a cosmological framework.
We have also exploited other kinematic spaces suggested in the literature to reconstruct the accretion history of our Galaxy, namely L perp − L z , eccentricity − L z and action space, and this additional analysis only confirms the problems encountered in disentangling GCs in the E − L z space: clusters with different origins appear scattered and mixed together also in those spaces.By means of Gaussian Mixture Models, we have demonstrated that the overlap of clusters is not only a projection effect on specific planes but it is found also when the whole set of kinematic properties (i.e.E, L z , L perp , eccentricity, radial and vertical actions) is taken into account.Consequently, applying algorithms such as Gaussian Mixture Models -a method lately used to identify groups of Galactic GCs in kinematic spaces -is conceptually wrong if the results of such methods are used to infer the merger tree of the Galaxy, and the interpretation of the results of these models tends to support the classification based on dynamical coherence.These findings, together, question the history of accretions experienced by the Galaxy, as it has been reconstructed so far by analysing its globular cluster system (see Massari et al. 2019;Kruijssen et al. 2020;Malhan et al. 2022).Indeed, this reconstruction usually assumes a dynamical coherence for clusters originating by the same accreted galaxy -assumption theoretically motivated by the Helmi & de Zeeuw (2000) work, which however was based on a number of oversimplifications, as extensively discussed in Jean-Baptiste et al. (2017).When more realistic simulations are analysed, indeed, the dynamical coherence of the accreted stars and globular clusters is not guaranteed anymore.When we test one of this merger histories (specifically the first one, proposed by Massari et al. 2019) we find no confirmation that the association that has been made so far in the literature between clusters and their progenitor galaxies is confirmed by the analysis of the age-metallicity relation of Galactic globular clusters.
To understand the origin of the globular clusters population in the Milky Way, we need to exploit the information from other dimensions like detailed chemical abundances and ages coming from the spectroscopic surveys such as the APOGEE survey and the soon operational WEAVE@WHT and MOONS@VLT sur- Table 4: Metallicities, ages and errors on ages (when available) for all Galactic globular clusters studied in Massari et al. (2019).Ages and metallicities are taken from (1) Marín-Franch et al. (2009), (2) VandenBerg et al. (2013).Note that in the case of (1) ages are relative ones, while for (2) they are expressed in units of Gyr.For each cluster, its progenitor galaxy, as reported by Massari et al. (2019), is also given.

Appendix A Static Milky Way potential
To illustrate how the distribution in the E − L z plane changes when the dynamical friction experienced by the satellite during the interaction with the Milky Way is not taken into account, we show here the results obtained by considering a static potential for the MW.To this aim, we have run the same simulations with one or two accreted satellites this time keeping the positions of the MW particles fixed at the initial values.
Figure 16 is the analogue of Fig. 1 thus showing the globular clusters projections in the (x, y) and (x, z) planes (left and central columns), for the single-accretion simulation with Φ orb = 60 • (i.e.simulation ID = MWsat_n1_Φ60), at different times: the initial (T = 0 Gyr), two intermediate times (T = 1 Gyr, 2.5 Gyr), and at the final time (T= 5 Gyr).In all the plots, the in-situ globular cluster system is represented by grey circles, and the globular cluster system initially linked to the satellite by magenta circles.The distribution of field (in-situ and accreted) stars is also shown in the background.Fig. 17 shows the temporal evolution of the distance of the satellite to the main galaxy, together with the corresponding evolution of all the clusters initially bound to the satellite.Due to its quasi-parabolic orbit, in this scenario we can see the satellite approaching the fixed MW-type galaxy at T 0.4 Gyr and then moving away again.At this approach, the satellite loses a globular cluster that is trapped by the Milky way on an orbit with apocentre at about 35 kpc and, as a result of its consequent departure, it retains the rest of its globular clusters.This can also be seen in Fig. 16 (left, middle columns) where indeed we notice that at T = 1 Gyr the satellite is receding (at T = 2.5, 5 Gyr it is no longer visible within the 100 kpc) after having deposited its globular cluster and part of its stars in the Milky Way, which at the end of the simulation (T = 5 Gyr) constitute the halo of the MW.
The corresponding evolution of the distributions in the E − L z space (right column of Fig. 16) of globular clusters (magenta circles) and field stars (background density map) of the accreted satellite are very different from the full N-body case (see Fig. 1, right column).At the beginning of the simulation the satellite, being still a bound system, appears to be clumped.With the passing of time, however, since in this case the positions of the MW particles are fixed and energy redistribution is not possible, i.e. there is no dynamical friction, the orbital energy of the satellite does not decrease as the satellite does not penetrate deeper in the potential well of the main galaxy.As a consequence, at the end of the simulation, we do not see a distribution of satellite stars and GCs in the E − L z space with a funnel-like shape elongated towards low energies as in the full N-body scenario.The energy fluctuates instead of decreasing, as we can also see in Fig. 17 (bottom panel) showing the temporal evolution of the orbital energy of the 10 satellite globular clusters.In fact the mean energy of the GCs remains equal to -3.3 as the initial one.Angular momentum is not conserved since the mass distribution of the system is not axisymmetric and changes over time.Moreover, as we have seen, the satellite does not merge with the Milky Way and this implies a drift in time of the L z for the stars and clusters that remain attached to it.This trend is clearly visible in Fig. 17 (middle panel) where, except for the one GC trapped by the MW for which L z fluctuates, for all the satellite's clusters the L z increases over time as it moves away.The mean L z in fact goes from being equal to 34.5 at the beginning of the simulation to being equal to 203.1 at the end, and the spread in L z also rises (as the standard deviation increases from 27.0 to 92.7).
As a simulation with two accretions, we re-run the one with ID = MWsat_n2_Φ30 -150 with a fixed MW potential.Figure 18 shows the final distribution in the E − L z space of the MW (grey circles) and satellite(s) GCs (magenta and orange circles) for the single and double accretion simulations with fixed galactic potential.These plots summarise the arguments just presented as the satellites' GCs are generally positioned at high energies and drifted towards large angular momenta.We do not find welldefined groups of clusters at different E−L z levels resembling the regions populated by different galactic progenitors in Fig. 4 in Helmi & de Zeeuw (2000).Furthermore, in this scenario, MW's clusters end up in a typical disc-like distribution since they are not heated up to halo-type kinematics by the redistribution of the energy, and this prevents us to account for a more realistic overlap between accreted and in-situ GCs.

Appendix B 4x(1:100) mergers
In this section we show how globular clusters distribute in the E − L z space in the full N-body case when the MW-type galaxy accretes four satellites with relative mass ratios of 1:100 with respect to the Milky Way.In these simulations, the Milky Waytype galaxy has the same properties (in terms of masses and sizes of its components, and of number of particles adopted) as those of the main simulations described in Sect. 2. As for the satellites, each of them is made of 25 000 particles and contains a population of 5 globular clusters.We have explored 7 different inclinations for the initial orbital plane of the satellite relative to the Milky Way disc.Here we present the results of two such simulations following the accretion over a period of 5 Gyrs that differ, as before, in the initial inclination of the various satellites Φ orb and so are referred as  show the initial (left panel) and final (right panel) distributions of GCs and stars belonging to each accreted satellite for the two simulations.As we can see, all the satellites in both the simulations start with a clumpy distribution of stars and clusters.On the other hand, the final stellar distributions are different as we can identify three types of shapes: high-energy and angular-momentum spread distributions (see satellites 1,2 in Fig. 19 and satellites 1,3 in Fig. 20), which correspond to those satellites that do not end up merging but move away from the MW after an initial approach, more clumped distributions in energy-angular momentum (see satellite 2,4 in Fig. 20) corresponding to satellites that are at the beginning of the merging, and more elongated distributions in energy (see satellite 3,4 in Fig. 19) corresponding to satellites that have orbited the MW for a while and are towards the end of the merging.For some satellites (see for example satellite 2,3,4 in Fig. 20), we can also see structures that start from the main distribution and elongate towards higher and lower energies that correspond to the trailing and leading tails, respectively.If we look at how the clusters are distributed, we can see that, in all the accreted systems, most of them are positioned at the lower extremes in energy of the corresponding stellar distributions.
Fig. 21 shows the final E − L z distributions of globular clusters belonging to all the four accreted satellites (colour coded circles) and to the MW-type galaxy (grey circles) respectively for the MWsat_n4_Φ150-60-0-30 and MWsat_n4_Φ180-90-30-120 simulations.As a general trend, we can see that accreted GCs remain at high energies compared to the cases involving mergers with satellites of 1:10 mass ratio.This happens because the progenitors to which they belong, being less massive, cannot penetrate deep enough into the MW-type galaxy potential.The only satellites that are able to lose more energy and carry a part of their clusters at energies lower than -5 (see satellites 3,4 in  Fig. 19), deposit some GCs also at high energies thus resulting in a final distribution that is not clustered as expected.In fact, even in this case we cannot identify clumps of clusters belonging to the same progenitor, being overall mixed and spread over an extended range of angular momentum.

Fig. 1 :
Fig. 1: Left, middle columns: projections of the simulated globular clusters positions on the xy and xz planes, for different times (increasing from top to bottom) of the single-accretion simulation with Φ orb = 60 • (MWsat_n1_Φ60, see Tab. 2 for the initial parameters).The in-situ clusters are represented by grey circles and the accreted clusters by magenta circles.In the background, the surface density of the totality of the stars of the simulation is also shown.Right column: Distributions of accreted GCs and field stars of the same satellite in the E − L z space, for different times (increasing from top to bottom) of the single-accretion simulation MWsat_n1_Φ60.Each globular cluster is colour-coded according to its escape time from the progenitor satellite which is also specified by the number on the top right and increases from 0.42 Gyr for GC1 to 5 Gyr for GC10 (see Sec 3.2 for the definition of t esc ).

Fig. 3 :
Fig. 3: Final energy of all the accreted clusters in the single accretion simulations as a function of their escape time from their satellites (t esc ); the color-coding indicates the clusters' initial energy in the reference frame of their progenitor satellite.

Fig. 4 :
Fig. 4: Accretion of two satellites on the Milky Way-type galaxy (MWsat_n2_Φ30-150).Left, middle panel: Distribution in the E − L z space of globular clusters originally belonging to the two satellites (magenta and white colors, respectively).The density maps show the fractional contribution of stars from satellite 1 (left panel) and satellite 2 (middle panel) relative to the totality of stars.Righ panel: Distribution in the E − L z space of globular clusters originally belonging to the two satellites (magenta and white colors, respectively) and to the Milky Way-type galaxy (grey circles).The density map shows the fractional contribution of stars from the two satellites relative to the totality of stars.

Fig. 5 :
Fig. 5: E − L z distribution for the whole set of 1x(1:10) merger simulations at T = 5 Gyr.The distribution of the in-situ population of globular clusters has been stacked and is shown by the grey circles.The color-coding of accreted GCs is different according to the initial inclination of the satellite orbital plane with respect to the Milky Way-like reference frame, Φ orb .

Fig. 6 :
Fig. 6: E − L z distribution for the whole set of 2x(1:10) merger simulations at T = 5 Gyr.The distribution of the in-situ population of globular clusters has been stacked and is shown by the grey circles.The color-coding of satellites GCs (satellite 1: left panel, satellite 2: right panel) is different according to the initial inclination of the satellite orbital plane with respect to the Milky Way-like reference frame, Φ orb .

Fig. 7 :
Fig. 7: Final globular clusters distribution in kinematic spaces namely E − L z , L perp − L z , eccentricity − L z and action space for the simulation MWsat_n1_Φ60.The in-situ population of globular clusters is represented by grey circles while the accreted one by magenta circles.
Fig. 8: Final globular clusters distribution in kinematic spaces namely, E − L z , L perp − L z , eccentricity − L z and action space for the simulation MWsat_n2_Φ30-150.The in-situ population of globular clusters is represented by grey circles while the accreted one by magenta and orange circles respectively for satellite 1 and satellite 2.

Fig. 9 :
Fig. 9: Distribution in L perp − L z , eccentricity − L z and action space for the whole set of 1x(1:10) merger simulations at T = 5 Gyr.The distribution of the in-situ population of globular clusters has been stacked and is shown by the grey circles.The colorcoding of satellites GCs is different according to the initial inclination of the satellite orbital plane with respect to the Milky Way-like reference frame, Φ orb .

Fig. 10 :
Fig. 10: Distribution in L perp − L z , eccentricity − L z and action space for the whole set of 2x(1:10) merger simulations at T = 5 Gyr.The distribution of the in-situ population of globular clusters has been stacked and is shown by the grey circles.The color-coding of satellites GCs (satellite 1: left column, satellite 2: right column) is different according to the initial inclination of the satellite orbital plane with respect to the Milky Way-like reference frame, Φ orb .

Fig. 11 :
Fig. 11: Top, middle row: Final globular clusters distribution in kinematic spaces (E − L z , L perp − L z , eccentricity − L z and action space) for the simulation MWsat_n1_Φ60.The colour-coding of GCs is related to the different components retrieved applying the minimum BIC criterion in the Gaussian Mixture Model.Truly accreted globular clusters (i.e. with the true label given by our simulation) are bordered by magenta circles.Bottom panel: Confusion matrix obtained by the GMM with values normalised to the total number of GCs in each true class.

Fig. 12 :
Fig. 12: Top, middle row: Final globular clusters distribution in kinematic spaces (E − L z , L perp − L z , eccentricity − L z and action − space) for the simulation MWsat_n2_Φ30-150.The colour-coding of GCs is related to the different components retrieved applying the minimum BIC criterion in the Gaussian Mixture Model.Truly accreted globular clusters (i.e. with the true label given by our simulation) are bordered by magenta and orange circles respectively for satellite 1 and satellite 2. Bottom panel: Confusion matrix obtained by the GMM with values normalised to the total number of GCs in each true class.

Fig. 13 :
Fig. 13: Same as Fig. 7 with the addition of in-situ GCs with halo kinematics at the beginning of the simulation, represented as empty circles.
abundance scale.VandenBerg et al. (2013) analyse 55 clusters -many of which are also in the Marín-Franch et al. (2009) -whose ages and metallicities are reported in Table 2 of their paper.The resulting age-metallicity relations are shown in Fig. 14, where colours indicate the different galaxy progenitors to which Massari et al. (2019) associate the clusters: Fig. 14, and possibly NGC 5139, red empty circle in the same figure) we find a cluster identified as a disc cluster (NGC 6752, blue point), two clusters which have been associated to the Helmi stream (NGC 5272, NGC 6981, orange points, and possibly NGC 5904, orange empty circle), two High-Energy clusters (NGC 6584, and NGC 6934, cyan points), one cluster associated to Sagittarius (NGC 6715, green point) and one possible to Sequoia (NGC 3201, brown empty circle).At higher [M/H]> −0.8, the Pal 1 cluster, which is classified as a disc cluster, lies in between two clusters associated to the Sagittarius galaxy, Palomar 12 and Terzan 7. The only three clusters in this region that seem to slightly separate from the bulk of the distribution are NGC 4147, Arp 2 and NGC 6535 (relative ages lower than 1.0 and −1.6 < [M/H] ≤ −1.5), but their are still compatible within 2σ with the other clusters; 3. at [M/H]< −1.4,no distinction can be made, on the basis of the age-metallicity relation, among the clusters that have been classified as Sequoia, Gaia Sausage Enceladus, Bulge/Disc or Sagittarius clusters.To further probe that the accreted groups in the Massari et al. (2019) classification -which, we remind the reader, constitutes the backbone of many other classifications that have been published afterwards in the literature (see also Forbes 2020; Pfeffer et al. 2020; Callingham et al. 2022) -do not depict specific AMRs in terms of shapes or amplitude, in the bottom panels of Fig. 14 we report a fit to the AMR of each of the groups Massari et al. ( have recently confirmed the results of tailored N-body simulations (Jean-Baptiste et al. 2017; Amarante et al. 2022) about the large spread of accreted stars in E − L z plane.

Fig. 14 :
Fig. 14: Top panels: Age-metallicity relations of Galactic globular clusters from the literature: (Left panel): Marín-Franch et al. (2009); (Right panel): VandenBerg et al. (2013).Different colours and symbols in each plot indicate the different galaxy progenitors of these clusters, following the classification given by Massari et al. (2019): Gaia Sausage Enceladus (G-E), Sequoia (Seq), Helmi stream (H99), Sagittarius (Sag), Low-Energy (L-E), High-Energy clusters (H-E) are shown by solid circles, while tentative associations are shown by empty circles.Clusters classified as in-situ by Massari et al. (2019) (M-D and M-B groups) are shown with blue dots.Errors on ages are also reported.Bottom panels: Fits to the age-metallicity relations of accreted globular clusters in the G-E, H99, Sag and L-E groups.For each group, the mean fit to the data, as a function of metallicity, is shown, together with the corresponding standard deviation.Metallicities and ages as provided by Marín-Franch et al. (2009) are used for the fits shown in the left panel, while in the right panel we make use of ages and metallicities as reported by VandenBerg et al. (2013).The in-situ clusters are also shown for comparison (blue points).The two clusters identified by a purple cross in the bottom-left panel are NGC 6441 and NGC 6121, which have been excluded from the fit, as motivated in the text.Error bars are not reported in these two panels, to avoid having overly complex figures, but they have been taken into account in the fitting procedure (see text).

Fig. 15 :
Fig. 15: Time, t f , when the chemical enrichment started, as a function of the effective yield p for the four different groups reported in Fig. 14.Top panel: p and t f resulting from fitting the AMR of Marín-Franch et al. (2009); Bottom panel: p and t f resulting from fitting the AMR of VandenBerg et al. (2013).

Fig. 16 :
Fig. 16: Left, middle columns: projections of the simulated globular clusters positions on the xy and xz planes, for different times (increasing from top to bottom, same as Fig. 1) of the single-accretion simulation with Φ orb = 60 • (MWsat_n1_Φ60, see Tab. 2 for the initial parameters).The in-situ clusters are represented by grey circles and the accreted clusters by magenta circles.In the background, the surface density of the totality of the stars of the simulation is also shown.Right column: Distributions of accreted globular clusters and field stars of the same satellite in the E − L z space, for different times (increasing from top to bottom) of the single-accretion simulation with Φ orb = 60 • .

Fig. 17 :
Fig. 17: Top panel: time evolution of the distances of the satellite (black line) and its globular clusters (coloured lines) from the Milky Way-type galaxy, for the simulation with fixed MW potential MWsat_n1_Φ60.Middle panel: time evolution of the angular momenta of the satellite globular clusters.Bottom panel: time evolution of the orbital energies of the satellite globular clusters.

Fig. 18 :
Fig. 18: Distribution of globular clusters in the E − L z space at the end of the single (MWsat_n1_Φ60, top panel) and double (MWsat_n2_Φ30-150, bottom panel) accretions simulations.Milky Way's globular clusters are identified as grey circles while satellite(s) GCs are represented as magenta and orange circles.

Fig. 19 :
Fig. 19: Initial (left panel) and final (right panel) distributions of accreted globular clusters and stars of the same satellite in the E −L z space for the four accreted satellites having 1:100 mass ratio with respect to the MW-type galaxy (ID = MWsat_n4_Φ150-60-0-30).Note that the ranges on the x and y axes are the same for all panels.

Table 1 :
Masses, characteristic scale lengths and heights, number and mean masses of particles for the different components of the Milky Way-type galaxy and the satellite(s).All masses are in units of 2.3 × 10 9 M , distances in kpc.
(see right panel) is made of clusters which have been associated to different progenitors, but which have age and metallicities that are indistinguishable one from another.If we look at the left panel of Fig. 14, indeed, in between the group of clusters associated to Gaia Sausage-Enceladus (such as NGC 5286, NGC 6205, NGC 7089, NGC 2808, NGC 288, NGC 362, NGC 1261, NGC 1851, red colours in or VandenBerg et al. (2013) is less than 3, making any fit meaningless.In the case of the AMRs from Marín-Franch et al. (2009) data, we have also excluded two clusters which were assigned by Massari et al. p being the effective yield, t the look-back time and t f the time in the past where the chemical enrichment started.The result of this analysis are shown in the bottom panels of Fig. 14 and in Fig. 15.If we use the AMRs by Marín-Franch et al. (