This chapter presents the material published in Ferrone et al. (2023). The aim of this work is to provide the first prediction of the expected tidal debris from all Milky Way globular clusters and to characterize their spatial distribution. For this study, we used GCsTT, a predecessor of tstrippy, which follows the same philosophy of representing stars as test particles that experience the gravitational field but do not contribute to it.

Introduction

Globular clusters are the oldest gravitationally bound stellar systems in the Galaxy (Meylan and Heggie 1997). About 170 are currently known in the Milky Way (Vasiliev and Baumgardt 2021) and the census is still incomplete, particularly in the inner regions of the bulge and disk of our Galaxy, where dust extinction and high stellar number density limit detections. It is in these regions in particular that new globular cluster candidates have been recently discovered, especially thanks to the analysis of near-infrared surveys

(D. Minniti et al. 2011, 2021; Moni Bidin et al. 2011; Dante Minniti, Palma, et al. 2017; Dante Minniti, Geisler, et al. 2017; Dante Minniti et al. 2018, 2021; Gran et al. 2019, 2022; Garro et al. 2020, 2021; Garro, Minniti, Alessi, et al. 2022; Garro, Minniti, Gómez, et al. 2022). The current population of globular clusters is likely to merely represent the leftovers of an initially more numerous and more massive one that had been depopulated as a result of many disruptive processes (Gnedin and Ostriker 1997; Murali and Weinberg 1997b, 1997a; Vesperini and Heggie 1997; Fall and Zhang 2001). One of the main processes affecting the globular cluster population and its evolution in number, mass, and size is tidal stripping.

As all stellar systems are characterized by a finite size and defined orbit of the Galaxy, globular clusters are subject to tidal effects, which arise because the opposite sides of these systems experience a different gravitational acceleration. The long-term effect of this process strips the system of its most loosely bound stars, which redistribute themselves onto orbits similar to those of their progenitor, forming so-called “tidal tails” or streams around it (see Carl J. Grillmair et al. 1995; Leon, Meylan, and Combes 2000, for some of the earliest studies). Some spectacular tails have been discovered and studied over the past twenty years around Milky Way globular clusters, ranging from the long tails (of roughly \(30^\circ\) degrees) departing from the Palomar 5 cluster (Odenkirchen et al. 2001, 2003, 2009; C. J. Grillmair and Dionatos 2006a; Guillaume F. Thomas et al. 2016; Starkman, Bovy, and Webb 2020; R. Ibata et al. 2021) to those of NGC 5466 (Belokurov et al. 2006), Palomar 14 (Sollima et al. 2011), and the GD-1 stream, whose parent cluster has still to be discovered (or has already been completely destroyed, leaving behind the stream as the only vestige of its past existence; see C. J. Grillmair and Dionatos 2006b; Webb and Bovy 2019; Ana Bonaca, Conroy, et al. 2020). These studies have been boosted in the last few years thanks to the publication of the ESA Gaia mission catalogues (Gaia Collaboration et al. 2016, 2018, 2021a, 2021b) which has been delivering parallaxes, proper motions, and magnitudes for about 1.4 billion stars, as well as the radial velocities for several million, and thus allowing for searches of stars with coherent distances and motions in the Galaxy, revealing the existence of a number of new and spectacular streams, as well as rediscovering and confirming already known ones (Navarrete, Belokurov, and Koposov 2017; Malhan, Ibata, and Martin 2018; Malhan et al. 2018; Rodrigo A. Ibata et al. 2018, 2019; Piatti 2021, 2021, 2022; N. Shipp et al. 2018; Rodrigo A. Ibata, Malhan, and Martin 2019; Kaderali et al. 2019; Bianchini, Ibata, and Famaey 2019; Malhan, Ibata, Carlberg, Bellazzini, et al. 2019; Malhan, Ibata, Carlberg, Valluri, et al. 2019; Palau and Miralda-Escudé 2019, 2021; Piatti and Carballo-Bello 2019, 2020, 2020; Caldwell et al. 2020; R. Ibata et al. 2020, 2021; Piatti et al. 2020, 2021; Nora Shipp et al. 2020; Guillaume F. Thomas et al. 2016; Boldrini and Vitral 2021; Malhan, Valluri, and Freese 2021; Jensen et al. 2021; Yuan et al. 2022; Zhang, Mackey, and Da Costa 2022; Nie et al. 2022). For a general overview, Mateu (2023) provides a recent compilation of known stellar streams.

All these studies are unraveling a very complex and rich set of stellar structures in the Milky Way that are mainly distributed in the halo, where their identification is the easiest because of the low density of the background stellar field.

From a numerical and theoretical point of view, many studies over the years have been focused on the formation and evolution of tidal streams around globular clusters (Keenan and Innanen 1975; Oh, Lin, and Aarseth 1992; Oh and Lin 1992; C. J. Grillmair 1998; Combes, Leon, and Meylan 1999; R. A. Ibata et al. 2002; Johnston, Spergel, and Haydn 2002; Yim and Lee 2002; Capuzzo Dolcetta, Di Matteo, and Miocchi 2005; Di Matteo, Capuzzo Dolcetta, and Miocchi 2005; Montuori et al. 2007; Siegal-Gaskins and Valluri 2008; Küpper et al. 2010; Lane et al. 2010; Küpper, Lane, and Heggie 2012; Mastrobuono-Battisti et al. 2012; Sanders and Binney 2013; Bovy 2014; Amorisco et al. 2016; Erkal et al. 2016; Sanders, Bovy, and Erkal 2016; Pearson, Price-Whelan, and Johnston 2017; Carlberg 2018, 2020; G. F. Thomas et al. 2018; Vitral and Boldrini 2022). These studies have contributed to understanding how these structures form and evolve, to what extent they trace the globular cluster orbit, and how their shape, extension, and morphology depend on the orbital phase and characteristics of the Galactic potential, as well as on the potential tidal shocks experienced by the cluster itself when it crosses the Galactic disk.

Some works have presented models and simulations for specific streams (Dehnen et al. 2004; Mastrobuono-Battisti et al. 2012; Banik and Bovy 2019; Ana Bonaca et al. 2019; Banik et al. 2021; A. Bonaca 2021), contributing to an understanding of their morphology, density variations, and their extent. From these works, it is clear that the tidal loss of stars from globular clusters and the formation of related structures are important for several reasons: (1) in quantifying to what extent globular clusters have contributed to the field stellar populations, from the halo to the disk to the bulge, and to what extent they still do; (2) reconstructing the properties (in terms of numbers and masses) of the early Galactic globular clusters, through their current mass loss; and (3) using globular cluster streams as a probe of the Galactic potential and, more generally, of the physical laws governing gravity (see, e.g., G. F. Thomas et al. 2018; Bianchini, Ibata, and Famaey 2019; Naik et al. 2020; Banik et al. 2021).

In this chapter, we wish to contribute to the current discourse on this matter by presenting the first complete catalog of simulated extra-tidal features around globular clusters. We emphasize that we are speaking generically on their features, rather than specifically on tails or streams, because the latter are but one of the morphologies that extra-tidal material can reveal, as we go on to show in this work. This project is motivated, on the one hand, by the aforementioned discoveries of many numerous new streams and tails in the Galaxy and, on the other hand, by the availability of the full 6D phase information and internal parameters (masses and sizes) for more than 150 Galactic globular clusters (Baumgardt and Hilker 2018; Baumgardt and Vasiliev 2021; Vasiliev, Belokurov, and Erkal 2021). The aims of this project are manifold: (1) to obtain a complete view of the expected distribution of globular clusters tidal structures in the sky; (2) to inform the interpretation of recent and future discoveries; (3) to support the search for new extra-tidal features in the data; (4) to offer the community a repository of all these models to be compared to other theoretical and numerical predictions, which adopt different Galactic potentials and/or gravity laws.

Numerical method

To model the formation and evolution of extra-tidal features around Galactic globular clusters, we use a set of codes, called Globular Clusters’ Tidal Tails (GCsTT) developed by our group. It comprises two python codes, for the backward and forward integration of a stellar system, made of N test-particles (see Sect. 2.1). These codes are separated for data organization and management, while the (computationally) most expensive part, namely, the calculation of the accelerations acting on the N particles and the orbits integration, is realized by means of a Fortran module written by our group. This module is interfaced to python by means of f2py directives from NumPy. The use of test-particle methods for modeling the tidal stripping process is widespread in the literature, where these methods are usually applied to one or few clusters at a time (see, e.g., Lane, Küpper, and Heggie 2012; Mastrobuono-Battisti et al. 2012; Palau and Miralda-Escudé 2019; Piatti 2021; Carl J. Grillmair 2022). In this work, we apply a test-particle methodology to the whole set (159) of Galactic globular clusters for which this is currently possible, also taking into account, for each cluster, the errors on astrometry, line-of-sight velocities¹ and distances. In the following, we describe the two main steps of the procedure used by GCsTT to simulate the tidal stripping process (Sect. 2.1), the initial conditions adopted for the clusters’ parameters and their mass distribution (Sect. 2.2), as well as the Galactic potentials (Sect. 2.3).

Simulations of the tidal stripping process: Two-step procedure

To model the formation and evolution of extra-tidal features around Galactic globular clusters, and predict their current properties, we proceed as follows:

Step i: Backward integration. Reconstructing the globular cluster orbit over the last 5 Gyr: First, for each Galactic globular cluster for which the distances from the Sun, proper motions, line-of-sight velocities, and structural parameters are available (see Sect. 2.2), we determine their current positions and velocities in a Galactocentric reference frame, in which the Sun is at \((x_\odot, y_\odot, z_\odot) = (-8.34, 0., 0.027)\) kpc (Chen et al. 2001; Reid et al. 2014) and at a given velocity for the local standard of rest, \(v_{LSR}= 240\) km/s (Reid et al. 2014), and a peculiar velocity of the Sun with respect to the LSR, \((U_{\odot}, V_{\odot}, W_{\odot}) = (11.1, 12.24, 7.25)\) km/s (Schönrich, Binney, and Dehnen 2010). We then integrate the orbit of a single point mass, representing the cluster barycenter, backwards in time for 5 Gyr, and in this way, we retrieve its position and velocity at that time in the chosen Galactic potential (see Sect. 2.3). We notice that other choices for the Sun’s position or velocity with respect to the Galactocentric frame would have been possible. For example, Piatti (2021) adopted the same values as ours for the \(v_{LSR}\) and for the peculiar velocity of the Sun but with a different distance to the Galactic center (8.1 kpc in their work, see GRAVITY Collaboration et al. 2018). The difference in the adopted position of the Sun is, however, generally smaller than the uncertainties affecting our knowledge of the distance of Galactic globular clusters to the Sun. For this reason, we do not to explore the dependency of the results presented in this paper with regard to these choices.

Step ii: Forward integration. Test-particle streams from the past to the present day: Once the positions and velocities of the barycenter of each cluster, 5 Gyr ago, have been determined, we build the corresponding \(N\)-body system, with N = 100 000 particles. The phase-space coordinates of these particles are generated following a Plummer distribution, with the total mass and half-mass radius as described in Sect. 2.2. The barycenter of this \(N\)-body cluster is then assigned initial positions and velocities in the Galactic model, as those retrieved at step \((i),\) and the cluster is then integrated forward in time until the present day. Particles in this \(N\)-body system are modeled as test-particles, that is, they experience the gravitational field exerted by the globular cluster itself (see Sect. 2.2) and by the Galaxy (see Sect. 2.3), but do not generate any gravitational field themselves. This allows us to maintain a computational time which scales as \(O(N)\) and not as \(O(N^2)\), as would be the case for a direct \(N\)-body self-consistent computation.

In the following, we refer to these simulations, made by using the most probable values on distances, proper motions, and line-of-sight velocities, as the “reference simulations.” In addition, for each globular cluster, we also take into account the errors on its distance, proper motions, and line-of-sight velocity, assuming Gaussian distributions of the errors, treated independently, and by generating 50 random realizations of these parameters. For each of these realizations, we repeat the steps described above, that is: (step i) we determine the associated current positions and velocities in the chosen Galactocentric reference frame, we integrate the orbit of the single-point mass (representing the cluster barycenter) backwards in time, retrieving the corresponding values 5 Gyr ago, (step ii) we build an \(N\)-body cluster containing \(N\)= 100 000 particles, with total mass and half-mass radius as those used for the reference simulation, and then we integrate the \(N\)-body cluster forwards in time until the present-day position.

To summarize, for a given Galactic potential, we run \(159\times (50+1)=8109\) simulations, where 159 is the total number of clusters for which we currently have both 6D phase-space information and structural parameters. As we discuss in the following section, the whole set of globular clusters has been evolved in three different Galactic potentials, which implies that a total of 24 327 simulations have been run.

For the orbit integration, a leap-frog algorithm is used, with a fixed time-step, \(\Delta t\), and a total number of steps, \(N_{steps}\), such that the total simulated time is \(\Delta t \times N_{steps}=5\) Gyr. The choice of the value of \(\Delta t\) adopted to simulate each cluster in the Galactic potential has been based on the energy conservation of the corresponding cluster evolved in isolation (i.e., without the effect of the Galactic gravitational field for 5 Gyr). For the majority of the clusters (109/159), this value was set to \(\Delta t = 10^5\) yr (for a corresponding value of \(N_{steps}=50\,000\)), while for the remaining clusters (50/159) a \(\Delta t = 10^4\) yr (for a corresponding value of \(N_{steps}=500\,000\)) was used. We refer to Supplemental Material 5.1 (and in particular to Table [tcross-energy]) for additional details on the choice of \(\Delta t\) for the whole set of clusters. As for the total simulated time, while globular clusters are much older than 5 Gyr, we chose this time limit because the longer back in time we could go, the less certain we would be of the Galactic environment. In addition, the last significant mergers in the Galaxy happened between 9 and 11 Gyr ago (see Belokurov et al. 2018; Helmi et al. 2018; Di Matteo et al. 2019; Gallart et al. 2019; Kruijssen et al. 2020) – well before the time interval simulated in this study. Other more recent interactions, such as the accretion of Sagittarius and of the Magellanic Clouds, may perturb the Galactic potential as well (see, e.g., Vasiliev, Belokurov, and Erkal 2021) and we plan to investigate their impact on the properties of globular cluster streams in the future.

For each realization, we generate an output file in an hdf5 format² containing the values for the right ascension (\(\alpha\)), declination (\(\delta\)), distance from the Sun (\(D\)), along with the components for proper motion in the equatorial coordinate system (\(\rm \mu_{\alpha}\cos(\delta)\) and \(\rm \mu_\delta\)), the line-of-sight velocity (\(\rm v_{\ell os}\)), longitude (\(\ell\)), latitude (\(b\)), as well as the components for proper motion in the Galactic coordinate system (\(\rm \mu_{\ell}\cos(\mathit{b})\) and \(\rm \mu_b\)) and the Galactocentric positions (\(x, y, z\)), velocities (\(v_x, v_y, v_z\)) and energy, \(E\), of each particle in the simulated system. We used Astropy (Astropy Collaboration et al. 2013, 2018) to convert the Galactocentric positions and velocities in the equatorial and Galactic quantities \(\alpha, \delta, D, \rm \mu_{\alpha}\cos(\delta), \rm \mu_\delta, \rm v_{\ell os}, \ell, b, \rm \mu_{\ell}\cos(\mathit{b})\), and \(\rm \mu_b\).

For each particle, we also save its escape time \(t_{\rm esc}\), defined as the time at which the particle escapes from the cluster, that is, the time, \(t,\) at which the particle satisfies the relation³: \[E_{GC}= 0.5 \times \left( (v_x-v_{x, GC})^2+(v_y-v_{y,GC})^2+(v_z-v_{z,GC})^2\right)+\Phi_{GC} > 0,\] with \(E_{GC}\) being the total specific energy of the particle relative to the cluster, that is, the sum of the potential energy, \(\Phi_{GC}\), due to the gravitational field of the cluster (see Eq. [gcpot]), and of the kinetic energy, relative to the cluster barycenter, \(T_{GC}=0.5 \times \left( \left(v_x-v_{x, GC}\right)^2+\left(vy-v_{y,GC}\right)^2+\left(vz-v_{z,GC}\right)^2\right)\), where \(v_x, v_y\), and \(v_z\) are its velocity components at time, \(t\), and \(v_{x,GC}, v_{y,GC}\), and \(v_{z,GC}\) of the cluster barycenter at the same time. A positive value of \(E_{GC}\) implies that the particle is no longer gravitationally bound to the cluster and, hence, it is lost in the field. Overall, the total volume of the whole set of 24 327 simulations, saved in hdf5 format, amounts to about 370 Gb.

Simulations of the tidal stripping process: Globular clusters’ current and initial conditions and their gravitational potential

Steps \((i)\) and \((ii)\) described in the previous section require some input conditions to be adequately executed. The current distances from the Sun, proper motions, and line-of-sight velocities, as well as the related uncertainties, of all 159 globular clusters considered in this study are taken respectively from Baumgardt and Vasiliev (2021) and Vasiliev and Baumgardt (2021). These values are then converted into Galactocentric positions and velocities by making use of Astropy and used as initial conditions to execute step \((i)\).

Step \((ii)\) requires generating an \(N\)-body system, representing the globular cluster, whose initial total mass and half-mass radius are assigned on the basis of their current values, as given by Baumgardt and Hilker (2018)⁴ and reported in Table 1. As anticipated at step \((ii)\) in Sect. 2.1, the phase-space coordinates of each \(N\)-body cluster are generated by assuming a Plummer distribution of total mass, \(M_{GC}\), and half-mass radius, \(r_{h}\), for which the corresponding potential is: \[\label{gcpot} \Phi_{GC}(r) = -\frac{GM_{GC}}{\sqrt{r^2+{r_c}^2}},\] where \(r_c\) is the cluster scale radius and it is related to the half-mass radius, \(r_{h}\), through \(r_{h} \simeq 1.305 r_c\) (Heggie and Hut 2003). The variable \(r\) here indicates the distance of the test particle from the center of the cluster. For each cluster, the same Plummer distribution used to generate the \(N\)-body system is also used to calculate the accelerations exerted on each particle as the system moves through time. The Plummer sphere, representing the cluster potential, indeed moves through the Galaxy along the orbit retrieved at step (i), traveling this time in the opposite direction, from 5 Gyr ago to the present day.

It might be noted that this implies that the globular cluster density profile and its internal parameters (total mass and characteristic radius) are constant over time in these models. This is, of course, a crude approximation, because in reality both the internal parameters and the density profile itself can change over time. We consider these assumptions to be acceptable within the scope of our work given that we are primarily interested in the distribution of extra-tidal stars, which had once escaped from the cluster have dynamics primarily dictated by the Galactic potential rather than the globular cluster itself. Of course, the density of stars along the extra-tidal structures, as well as the total mass lost, depend on these assumptions. That is to say that if the mass of the cluster was not assumed constant over time, but could possibly decrease, the gravitational attraction exerted by the cluster itself on its stars would be weaker and this would lead to an increasing mass loss and density along the tails. We could have proceeded with diminishing the mass over time, based on some assumptions on the temporal behavior of this relation, however, we did not find this approach satisfying. In this way. we would have taken into account a temporal evolution of the mass, but not of the size of the cluster, adding a supplementary hypothesis to the problem. For these reasons, we decided to maintain the simplest approach. We emphasize that other groups have followed the same methodology, maintaining masses and sizes that remain constant over time (see, e.g., Palau and Miralda-Escudé 2019).

The summary tables giving both the current internal parameters of the clusters (total mass and half-mass radius), their astrometric quantities of relevance for this study and the line-of-sight velocities are publicly available⁵. We have made use of these tables for our work and we report them in a unique table in our paper for the sake of the completeness and self-consistency of the data used (see Table 1).

Simulations of the tidal stripping process: Galactic potentials

As for the Galactic mass distribution, we make use of the two axisymmetric Galactic mass models presented in Pouliasis, Di Matteo, and Haywood (2017) and of an asymmetric mass model, containing a central stellar bar, and we present it here for the first time. We recall the main properties of the two models of Pouliasis, Di Matteo, and Haywood (2017) below and we describe the asymmetric Galactic mass model, presented here for the first time, in more detail.

Model I by Pouliasis, Di Matteo, and Haywood (2017): An axisymmetric mass model for the Galaxy including a spherical bulge

Model I by Pouliasis, Di Matteo, and Haywood (2017) (abbreviated name: PI) consists of four components: two disks (thin and thick), both described by Miyamoto & Nagai potentials, a dark matter halo, and a central bulge. Its total potential is: \[\Phi_{tot}(R, z) = \Phi_{thin}(R, z) + \Phi_{thick}(R, z) + \Phi_{halo}(r)+ \Phi_{bulge}(r),\] with \(r=\sqrt{R^2 + z^2}\), \[\begin{aligned} \Phi_{thin}(R,z)&=&\frac{-GM_{thin}}{\left(R^2+\left[a_{thin}+\sqrt{z^2+b_{thin}^2}\right]^2\right)^{1/2}},\\ \Phi_{thick}(R,z)&=&\frac{-GM_{thick}}{\left(R^2+\left[a_{thick}+\sqrt{z^2+b_{thick}^2}\right]^2\right)^{1/2}}, \end{aligned}\] \[\begin{split} \Phi_{halo}(r)=&\frac{-GM_{halo}}{r}-\frac{M_{halo}}{1.02a_{halo}}\times\\ & \Bigg[\frac{-1.02}{1+\left(\frac{r}{a_{halo}}\right)^{1.02}}+\ln{(1+\left(\frac{r}{a_{halo}}\right)^{1.02})}\Bigg]_R^{100}, \end{split}\] and \[\Phi_{bulge}(r) = -\frac{GM_{bulge}}{\sqrt{r^2+{b_{bulge}^2}}},\] where \(M_{thin}, M_{thick}\), \(M_{halo}\), and \(M_{bulge}\) are the masses of the disks, halo, and bulge. Also, \(a_{thin}, b_{thin}, a_{thick}, b_{thick}, a_{halo} ,b_{bulge}\) are the characteristic scale lengths of the thin and thick disks, the halo, and the central bulge, respectively (see Table [PII]).

This model is a modification of the classical Allen and Santillan (1991) model, made to include also the presence of a thick disk. As it has been discussed in detail by Pouliasis, Di Matteo, and Haywood (2017), the choice to include a massive spheroid in this model, as well as in the original Allen and Santillan (1991) model, is dictated by the need to reproduce CO/HI-based velocity curves, as those provided by Sofue (2012), which show a rise and then a sudden decrease of the velocity curve in the inner Galactic regions (\(R \le 2-3\) kpc). In an axisymmetric model, such a rise can be reproduced only if a central spheroidal component, with a typical mass greater than 10% of that of the disk(s), is added. However, as shown by Chemin, Renaud, and Soubiran (2015), the central rise observed in the rotation of the molecular gas in the inner Galaxy may be an effect of non-circular motions generated by large-scale asymmetries such as the bar. Moreover, this feature is not reported in all the observational studies (see, e.g., Reid et al. 2014, on which model PII is based). In other words, if we do not assume that the mass distribution of the inner Galaxy is axisymmetric, the need for a massive spheroidal component to reproduce velocity curves, such as those from Sofue (2012), no longer persists. In addition to that, in the last decade, a number of works have shown that if a spheroidal bulge exists in the central regions of our Galaxy, it has to be small (few percents of the mass of the disk at the most, see among others Shen et al. 2010; Kunder et al. 2012; Di Matteo et al. 2015; Gómez et al. 2018). All these arguments suggest to employ this model (as well as all models including a massive central spheroid; see, e.g., Irrgang et al. 2013) with care when dealing with the central parts of the Galaxy. Since models with a massive central spheroid, however, are still used in the literature, we have included model PI here, as a term of comparison.

Model II by Pouliasis, Di Matteo, and Haywood (2017): An axisymmetric, bulge-less mass model for the Galaxy

Model II by Pouliasis, Di Matteo, and Haywood (2017) (abbreviated name: PII) consists of a spherical dark matter halo, with the same functional form adopted in the Allen and Santillan (1991) model, and of two disk components (a thin and a thick disk), with same functional form as PI. This model does not include any central spheroid (i.e., it is a bulge-less model) and thus its total potential is the sum of three components only: \[\Phi_{tot}(R, z) = \Phi_{thin}(R, z) + \Phi_{thick}(R, z) + \Phi_{halo}(r),\] with the thin, thick disks, and dark matter halo having the same functional forms adopted in PI.

As it has been shown in Pouliasis, Di Matteo, and Haywood (2017), this model satisfies a number of observational constraints, such as the stellar density at the solar vicinity, thin- and thick-disk scale lengths and heights, the rotation curve as provided by Reid et al. (2014) and the absolute value of the perpendicular force, \(K_z\), as a function of distance to the Galactic centre (see Sect. 2.5 in Pouliasis, Di Matteo, and Haywood 2017). As it is, however, an axisymmetric model, it fails to accurately describe the inner few kpc of the Galaxy, where the stellar mass distribution has been shown to be asymmetric.

Model II with a massive, slowly rotating stellar bar

The third mass model (abbreviated name: PII-0.3-SLOW) that we use in this paper is a version of PII by Pouliasis, Di Matteo, and Haywood (2017) modified to include a rotating stellar bar, whose mass has been assigned to be 30% of the (thin+thick) disk mass of PII. We assume that the bar rotates with a constant pattern speed of \(\Omega_{bar}=38\rm km~s^{-1}kpc^{-1}\) and that it is currently inclined of \(25^\circ\) with respect to the Sun-Galactic center direction (see Bland-Hawthorn and Gerhard 2016). We model it as a triaxial distribution, whose gravitational potential is given by Long and Murali (1992): \[\Phi_{bar}(x,y,z)=\frac{GM_{bar}}{2a_{bar}}\ln\left( \frac{x-a_{bar}+T_{-}}{x+a_{bar}+T_{+}} \right),\] with \(T_{\pm}=\left[ (a_{bar}\pm x)^2+ y^2 + (b_{bar}+ \sqrt{c_{bar}^2+z^2})^2 \right]^{1/2}\) and \(a_{bar}, b_{bar}, c_{bar}\) the characteristic bar parameters. The total gravitational potential generated by this model thus takes the form: \[\Phi_{tot}(x, y, z) = \Phi_{thin}(R, z) + \Phi_{thick}(R, z) + \Phi_{halo}(r)+ \Phi_{bar}(x, y, z).\] with all characteristic values given in Table [PII]. Practically, to include the bar, we reduced the mass of the disks in such a way to maintain the total stellar mass of this model as that of PII. Long and Murali (1992) provide the formulas of the accelerations generated by this triaxial distribution in the reference frame of the bar. To calculate and add them to the accelerations generated by the disks and dark matter halo, at each time step, we converted the positions of all particles in the rotating, non-inertial reference frame of the bar, computed the corresponding accelerations on each particle, and then transformed these accelerations back in the inertial reference frame described in Sect. 2.1. In this way, the accelerations due to the bar are added to those generated by the other terms of the Galactic mass distribution.

We emphasize that we do not consider this model as the best possible representation of the Galactic mass distribution, especially in the central region. It can, however, provide a first indication on how the inclusion of a rotating asymmetric component in the inner Galaxy can affect the globular cluster streams, near and far from the Galactic center.

Moreover, since the exact characteristics of the Milky Way bar are still subject to debate (see, e.g., Bland-Hawthorn and Gerhard 2016), it is important to explore how varying the parameters adopted in this paper, such as the pattern speed, the mass, or the length of the bar, can affect the characteristics of the whole set of streams. More complex shapes for the bar can also be explored, for example, by substituting the inner parts of the triaxial bar with a boxy-peanut-shaped morphology, which has been shown to characterize the inner Milky Way (see, e.g., Wegg and Gerhard 2013; Wegg, Gerhard, and Portail 2015). These topics are, however, beyond the scope of the present paper. In sum, given the uncertainties on the bar’s physical extent and how it can change over the time span investigated here, its affect on the streams presented here are purely indicative.

Results

A few basic premises:

Given the large number of simulations carried out and the wealth of information contained in them, it is not possible to exhaust all possible applications of this simulation database in this paper. We have therefore chosen to proceed as follows.

In Section 3.1, we present an overview of the distribution of all streams in Galactic coordinates. This coordinate space will be the one used in the remainder of the entire article. This first section allows us to show qualitatively how the global distribution of streams varies, depending on the Galactic potential used.

We then move on (Section 3.2) to present the global system of streams as a function of their distance from the Sun. In this section, we also show the kinematic properties of the streams, such as proper motions and line-of-sight velocities, that can be directly compared to Gaia data or other astrometric and spectroscopic surveys. This section also allows us to show the variety of morphologies that the stars which escaped from globular clusters can take. In Section 3.3, we explore this issue in more detail, showing how these morphologies depend primarily on the orbital characteristics of the clusters and their distance from the Galactic center. For the most interesting cases, we compare the tidal structures predicted by our simulations with streams found in observational data. For this purpose, we make use of the galstreams library of stellar streams in the Milky Way (Mateu 2023), which constitutes a unique and public database summarizing angular positions, distances, proper motions, and line-of-sight velocity tracks for nearly a hundred Galactic stellar streams. Any stream that is not included in this library is not be compared to our simulations in the context of this paper.

We note that the tidal features associated with each of the 159 simulated clusters are presented in Supplemental Material 5.2. To avoid making this Supplemental Material too long, the tidal features above are presented only in the case of the potential PII. However, all the data are available on a dedicated site⁶ where it is possible to the way the characteristics of these streams change, for any cluster, in the three chosen potentials.

A sky full of streams

Figure [galleryLB] shows the number and mass density distributions of the whole set of simulated globular clusters and their extra-tidal features in Galactic longitude and latitude for the three Galactic potentials.

For all Galactic mass models adopted, a striking characteristic among the plots in Figs. [galleryLB] is the variety of features that our models predict, which are reminiscent of the tidal tails, stellar streams, and shells that are produced by interacting and merging galaxies in the process of mass assembly (see, e.g., Mancillas et al. 2019). Some clusters have very thin and elongated streams, which describe arcs that can extend up to \(180^\circ\) in longitude, or tens of kpc in physical space. In some other cases, extra-tidal features appear shorter (a few up to ten degrees) and sometimes also thicker (about \(10^\circ\) in the sky) than others. Finally, in some cases, clusters are surrounded by extended structures, such as halos, rather than coherent and thin streams.

This variety of properties depends on several factors: the distance of the stream to the Sun (due to projection; i.e., for a given physical thickness, the closer the stream is to the Sun, the more extended it appears in the \((\ell, b)\) plane), its orbital phase (towards the peri-center or the apo-center of the orbit), and the orbital properties of the parent globular cluster. We also see from these figures that stellar particles stripped from their parent clusters do not only redistribute in coherent structures, but in some cases, they can also contribute to a more diffuse density distribution.

Because of the large number of simulated clusters, we have chosen not to present the corresponding extra-tidal features one-by-one in the main part of this paper, but we have rather decided to describe these extra-tidal features with a global approach, by first adopting a criterion based on the distance of these features to the Sun (see Sect. 3.2) and then discussing the types of distributions tidal debris can have and how these depend on the cluster orbital parameters (see Sect. 3.3). All the extra-tidal structures generated by the 159 globular clusters simulated in this paper, and their corresponding uncertainties, are reported in Supplemental Material 5.2. Among them, we include clusters with thin and elongated tails, such as IC 4499, NGC 3201, NGC 4590, NGC 5024, NGC 5053, Pal 5, to cite a few, as well as clusters such as AM 1, Pal 14, Pal 4, and Pal 15, whose extra-tidal material shows a halo-like configuration, and clusters such as NGC 1261, NGC 4147, NGC 6356, and UKS 1, whose stripped stars show a remarkable diffuse distribution in the field.

Finally, even when our models are not tailored to accurately reproduce the mass loss from globular clusters, since (as described in Sect. 2.2) we have adopted a test-particle approach with a time-constant globular cluster potential, it is however tempting to estimate, to a first order, the total mass associated with the tidally stripped population and compare it to the current mass. By calculating the mass lost in the field in the past 5 Gyr as the sum of the mass of all particles⁷ that have escaped the cluster (\(t_{\rm esc}> 0\)), we find that the PI model sheds \(2.1\times 10^{7} M_{\odot}\), which is 55% of the GC population’s current mass. The PII model shed is \(2.7\times 10^{6} M_{\odot}\), which is 7% of the current mass. Similarly, the PII-0.3-SLOW model lost \(3.7\times 10^{6} M_{\odot}\), which is 10% of the current mass and gives a half mass radius of 6.3 kpc.

This mass roughly constitutes one-hundredth to one-tenth of the total stellar halo mass (Bland-Hawthorn and Gerhard 2016) and it is probably only a lower limit on the mass of escaped stars in the field, since a number of clusters initially in the Galaxy must have been destroyed over time (see introduction) and are thus no longer identifiabl as globular clusters today. It is also interesting to note that escaped stars are mostly redistributed in the inner Galaxy, with the half-mass radius of the PI model being 4.0 kpc and that of the PII and PII-0.3-SLOW models being 6.3 kpc.

The total mass lost from the clusters, as well as its spatial distribution, thus depends on the Galactic potential adopted: the variations between the PI model and both the PII and PII-0.3-SLOW models are of course caused by the PI’s inclusion of the bulge, which leads to larger tidal forces in the center of the galaxy and subsequently drives larger mass loss.

Despite the differences in the modeling approach, it is interesting to compare our results to those of Baumgardt (2017). Briefly, our experiments differ in the following ways: they employ N-body simulations while we use our test-particle approach; they have an integration time of 12 Gyrs compared to our 5 Gyrs; lastly, their clusters have circular orbits in a Galactic potential modeled as an isothermal sphere as compared to the more realistic orbits and Galactic potentials considered in this work. Interestingly, the authors find that over 12 Gyrs their population of globular clusters lose two-thirds of their initial mass. This is roughly consistent with our application of the PI model, whose globular clusters shed 35% of their initial mass in 5 Gyrs, which is roughly half the mass found by Baumgardt (2017) in a period that is also about half as long.

From the nearest to the furthest extra-tidal structures

The analysis presented in this section along with the corresponding information are given in Figures [D0-10], [D10-20], [D20-30], and [D30-300] are restricted to the PII model.

In this section, we analyze structures located at different distances from the Sun. We identify structures as main or secondary on the basis of the fraction of stars they represent. If, for example, a cluster contributes more than 10% of its stripped stars in a given distance range, we consider the associated structure to be significant and we define it as a main structure in this distance range. If, on the other hand, the fraction of stripped stars from a cluster is between 1% and 10%, we define the associated structure as being “secondary”. Structures that constitute less than 1% of the cluster mass are considered insignificant in that range. The main and secondary structures for each distance bin are reported in Table [summarySTREAMS], where they are named by their progenitor cluster. Below, we describe the structures encountered at different distances from the Sun, from the closest to the most distant. For the following discussions, we refer to Figs. [D0-10] to [D30-300], where we report the different streams in a given distance range (top row in each figure), the corresponding mass density generated by all the structures found in that bin (second row), the longitudinal and latitudinal proper-motions map (third and fourth rows), and the line-of-sight velocities (bottom row).

The [0–2] kpc distance range:

Among the 159 simulated clusters, only NGC 6121 is currently at a distance of less than 2 kpc from the Sun. In addition to stars stripped from this cluster within this distance-to-the-Sun range, we have identified very diffuse and low density extra-tidal stars associated with five other clusters, which are currently several kpc away from the Sun, as reported in Table [summarySTREAMS]. Of these, only UKS 1 and NGC 6121 have a significant fraction of their stripped mass in this distance range. All the others contribute with only a few percent. None of these clusters seem to show well-defined, stream-like features in the (\(\ell,b\)) plane. This extra-tidal material appears indeed quite uniformly redistributed, in a latitude range mostly inside \(-30^\circ\) to \(30^\circ\), and is mostly found at negative longitudes. Because in this interval range there is no remarkable extra-tidal structure found on the sky, we decided to plot this distance bin together with the [2–5] kpc bin (see Fig. [D0-10], left column).

The [2–5] kpc distance range:

As the distance from the Sun increases, many more tidal structures are intercepted. Eleven globular clusters are found at a distance between 2 and 5 kpc from the Sun; and together with structures emanating from these clusters, we also find extra-tidal material associated with other 48 clusters, 19 of which are significant (as listed in Table [summarySTREAMS]). Some extra-tidal structures are clearly identifiable, even when over-plotted with all the others: this is the case, for example, for the tidal material associated with the E 3 cluster, which appears as an extended thick stream (as shown in Fig. [D0-10]) at approximately \(-100^{\circ}\) longitude. This is also the case for the complex tidal structure associated with BH 140, which resembles a ribbon in the sky, with a bifurcation at negative longitudes whose edges extend to about \(-30^\circ\) and \(30^\circ\) latitude. In addition to these, there are a series of circular halos concentric about the Galactic center, with abrupt drop-offs in density tracing the furthest extent of diffuse debris emanating from a variety of clusters. Overall, few streams are immediately recognizable in this range, although plenty of debris is present. This is expected given that this range samples a bite of the Sun-side of the Galactic disk and that this range has a relatively small volume. In this, as in the following distance ranges, the \(v_{\ell os}\) maps show that the tidal material at negative Galactic longitudes has, on average, positive \(v_{\ell os}\) while material at positive Galactic longitudes has, on average, negative \(v_{\ell os}\), which is due to the solar reflex velocity. This is the same trend observed for the whole set of radial (i.e., line-of-sight) velocities in Gaia DR3; for instance, we refer to the bottom panel of Fig. 5 from Katz et al. (2023), although less extreme velocities are reported in their plot since their dataset is dominated by disk stars, whereas our maps have a high proportional contribution of halo stars (additionally, they use median values in their bins, while we use an average). In order to quantify the net rotation of the system of streams, we calculated the mean angular momentum about the Galactic pole as \(\langle L_z \rangle =\sum_i L_{z,i} m_{p,i} / \sum_i m_{p,i}\), where \(m_{p,i}\) is the mass of each star particle indexed by \(i\) (as discussed in footnote 7) and \(L_{z,i}\) is the corresponding particle’s angular momentum. The mean angular momentum is found to be \(\langle L_z \rangle =-300\) kpc km s\(^{-1}\), which shows a slight co-rotation of the system of streams with the disk though there is much dispersion about this value – as shown in bottom right panel of Fig. [orbparam] in the Supplemental Material.

The [5–10] kpc distance range:

The [5–10] kpc distance range, which includes the Galactic center, contains much more material. A total of 79 clusters are found in this range, together with tidal structures associated with 131 different progenitors redistributed among main and secondary structures. Some tiny streams are visible in the density maps as well as in proper motions and line-of-sight velocity spaces (see Fig. [D0-10], right column): the trailing portion of the tail of the globular cluster Pal 1, at \((\ell, b)\sim(140^\circ, 25^\circ)\); the most extreme portion of the trailing tail of NGC 3201 at \((\ell, b)\sim(150^\circ, -37^\circ)\); the waterfall-like shape of NGC 288, particularly evident at \(b \lesssim -60^{\circ}\); the thin inverted U-shape of NGC 4590 at positive latitudes spanning a large longitude extent from \(\ell \simeq -60^\circ\) to \(100^\circ\); the portion of the E 3 tails the closest to the cluster at \((\ell, b)\sim(-75^{\circ},-15^{\circ})\), which continues from the more easily recognizable portion in the [0–5] kpc bin.

The [10–15] kpc distance range:

In the [10–15] kpc range, we find tidal structures associated with 134 progenitors, as listed in Table [summarySTREAMS] and reported in Fig. [D10-20] (left column), 27 of which are related to globular clusters that are also found in this distance bin. We note that at these distances from the Sun, the distribution of tidal features in directions towards the Galactic center, from \(-30^{\circ}\) to \(30^{\circ}\) longitude, appears less fuzzy than the one characterizing the [0–5] and [5–10] kpc distance bins. Tidal features here are beyond the Galactic center, and are mostly associated with disk or halo clusters.

Among the thinnest structures, we find the stream associated with Pal 1, at \((\ell, b)\sim(120^\circ, 15^\circ),\) which is also visible in the distance bin [5–10] kpc (see previous discussion), but which is even more elongated here. NGC 6101 shows the nearest portion of its long thin diagonal tidal tail that spans negative longitudes and ranges from \(-15^{\circ}\) to \(45^\circ\) latitude. Additionally, this stream is also unique against its counter parts in proper motion space. NGC 5053’s nearest portion appears as a vertical tidal tail at \(-80^{\circ}\) longitude. Similarly, NGC 5466 is shown vertically at \(25^{\circ}\) longitude.

Among the thickest structures, we can recognize general diffuse and bowtie-like shapes. There are also spoke-like structures departing radially from the Galactic center. For instance, we can associated the extra-tidal material with NGC 7078 at \((\ell, b)\sim(60^\circ,-28^{\circ}),\) as well as NGC 7089, which is nearly parallel to the previous structure, but at lower latitudes at \((\ell, b)\sim(50^\circ,-40^{\circ})\).

The [15–20] kpc distance range:

At larger distances (\([15--20]\) kpc range, see Fig. [D10-20]), some of the most striking features are associated with the clusters NGC 5024 and NGC 5053, whose long thin tails essentially overlap in this distance, with the latter covering the former, and appearing at high latitudes spanning a longitudinal range from \(\ell \sim -90^{\circ}\) to \(\ell \sim 45^{\circ}\). Again, the long thin stream of NGC 5466 appears in this range (as will be the case for the next) and is at high latitudes at roughly \(85^{\circ}\) and positive longitudes. There is also the thicker extended structure of NGC 4147 whose diffuse structure emanates from about \((\ell, b)\sim(-100^{\circ},80^{\circ})\).

In this distance bin, we find long tidal tails emanating from globular clusters associated with the Sagittarius dwarf galaxy, which are particularly visible at negative longitudes: a long thin stream is associated with Pal 12 at positive longitudes and latitudes \(b \le-15^{\circ}\), as well as two overlapping structures at \(0^\circ \lesssim \ell \lesssim 30^\circ\) longitude, namely, Ter 7 (and Ter 8 in the next distance bin). A word of caution is needed here: the mass loss from these clusters may be incorrect, since we do not include the presence of the Sagittarius dwarf galaxy itself. The potential well associated with this latter could change the tidal effects experienced by clusters associated with Sagittarius, especially in the case of NGC 6715, which sits at the center of this dwarf galaxy. The inclusion of the Sagittarius dwarf will be the subject of future investigations. Overall, in this distance bins, we find 11 clusters and 61 streams, all listed in Table [summarySTREAMS].

The [20–25] and [25–30] kpc distance ranges:

In the following distance bins (at [20–25] kpc and [25–30]; see Fig. [D20-30]), globular clusters and extra-tidal structures become less numerous, although some are still visible, such as Pal 5 at \((\ell, b)\sim(0^{\circ},45^{\circ})\). In more detail, in the [20–25] kpc bin we find tidal features associated with 37 different progenitors, 8 of which are associated with globular clusters whose current positions are in the same distance bin; in the [25–30] kpc bin, 7 clusters are found, together with tidal features associated with 30 other progenitor clusters which do not lie in this same distance range. In both bins, the streams emanating from globular clusters associated with the Sagittarius dwarf galaxy are still visible, as well as the most extreme portion of the tail associated with NGC 5466.

The [30–35] and [35–300] kpc distance ranges:

Finally, in the last distance bins (see Fig. [D30-300]), thin streams become rare. Some small streams are visible: Pyxis at \((\ell, b)\sim(-100^{\circ},0^{\circ})\); NGC 2419 at \((\ell, b)\sim(-180^{\circ},30^{\circ})\); Pal 4 at \((\ell, b)\sim(-160^{\circ},75^{\circ})\); Pal 3 at \((\ell, b)\sim(-120^{\circ},45^{\circ})\). Many more have a diffuse and halo-like structure. For instance, the blob associated with Pal 15, centered at \((\ell, b)\sim (15^\circ, 20^\circ\)); AM 1 at \((\ell, b)\sim(-100^{\circ},-55^{\circ})\); Eridanus at \((\ell, b)\sim(-140^{\circ},-45^{\circ})\); Pal 14 at \((\ell, b)\sim(30^{\circ},45^{\circ})\); Laevens 3 at \((\ell, b)\sim(65^\circ,-20^\circ)\), which is completely enveloped by NGC 7006. In total, in these two distance bins, we find 4 and 12 clusters, respectively, along with their associated streams, together with extra-tidal material associated with 14 and 19 progenitors in total.

Disks of inner and outer halo clusters: A variety of morphologies and shapes for extra-tidal structures

The analysis presented in the previous section allows us to appreciate the variety of morphologies found for extra tidal structures, from padlocks to “Easter eggs,” disks, ribbons, and canonical streams. Moreover, some structures are limited in latitude and longitude, while some others fill nearly the entire sky.

To more easily capture the similarity and differences in the morphology of the extra-tidal features surrounding Galactic globular clusters, we can group the latter on the basis of their orbital parameters⁸(see Supplemental Material 5.4 for more details), as follows: Disk clusters: A cluster is classified as a disk cluster if arctan(\(z_{max}/R_{max}\)) \(\le 10^\circ\), where \(z_{max}\) and \(R_{max}\) are, respectively, the maximum height above or below the Galactic plane reached by its orbit in the last 5 Gyr and its maximum in-plane distance from the Galactic center. Inner clusters: All clusters with \(r_{max} \le R_{\odot}\) that are not classified as disk clusters enter this group. Contrary to \(R_{max}\), which is the maximum in-plane distance that a cluster reaches from the Galactic center, \(r_{max}\) is the maximum 3D distance, that is, \(r_{max} = \textrm{max}(\sqrt{R^2+z^2}),\) with the maximum calculated over the whole cluster orbit. Outer clusters: All clusters with \(r_{max} > R_{\odot}\) that are not classified as disk clusters are included in this group.

By using the orbital radius of the Sun as the criterion for inner and outer clusters, debris from outer clusters can span the whole sky while inner clusters must be restricted in longitude and latitude. With these definitions, 21 clusters are disk clusters, 71 are inner clusters, and 67 are outer clusters (see Table [classification] in Supplemental Material 5.4). We emphasize that this classification does not aim to suggest any specific origin for these systems (e.g., whether they are in-situ or accreted, see Massari, Koppelman, and Helmi 2019), but it is uniquely based on their current orbital characteristics and helps in capturing some of the properties in the extension (projected in to the sky) and shape of their extra-tidal material, as we discuss in the following.

Extra-tidal features originating from disk clusters: ribbons in the Galactic plane

Disk clusters are defined on the basis of the flatness of their orbits (i.e., on their \(z_{max}/R_{max}\) ratio). As a result, they typically are restricted to low latitudes, though the exact distribution depends on the relationship of their orbit to the solar radius. To specify, clusters whose \(R_{max}\) are interior to the Solar radius generate tidal debris in a limited range in longitude and latitude. For instance, the material associated with clusters as Ter 1, Ter 5, Ter 6, and Ter 9 has a disk-like shape and is completely confined to \(|\ell| <30^{\circ}\) and \(|b|<10^{\circ}\). If \(R_{max}\) is greater than the solar radius, material can cover the full longitude space and most of the material will still appear at low latitudes. This is the case, for example, of BH 140, whose escaped stars diffusely occupy all longitudes and most of them are found at \(|b| \le 30^\circ\), while for Pal 2 and Pal 10, their extra-tidal stars have a very limited latitudinal extension and appear as ribbons in the sky. In the following, we discuss some of the structures associated with NGC 6121, Pal 2, and Pal 10. We refer to Supplemental Material 5.2 for the tidal structures generated by the whole set of disk clusters.

NGC 6121:

With a current position at \(x=-6.58\), \(y=-0.28\) and \(z=0.53\) kpc, NGC 6121 is the closest globular cluster to the Sun in our list. This cluster has a remarkably planar orbit, with a maximal vertical excursion from the Galactic plane of only 0.5 kpc (see Fig. [ngc6121_stream]), and an eccentricity \(e=0.80 \pm 0.01\), which makes it oscillate between an apo-center at \(R_{\rm max}=6.81 \pm 0.02\) kpc and a peri-center at \(R_{min}=0.76 \pm 0.04\) kpc. Because this cluster lies inside the solar circle, its orbit is limited to a longitude interval from \(-60^{\circ}\) to \(60^{\circ}\); because the cluster currently lies very close to the Sun, and is at its highest height above the Galactic plane, the orbit forms a hook-like pattern in longitude-latitude space. This hook-like portion of the orbit, nearest to the cluster, is traced by the recently stripped tidal material (see Fig. [ngc6121_stream]), with a leading tail oriented mostly in a vertical direction in the \((\ell,b)\) plane, from the current cluster location, up to approximately \(0^\circ\) latitude. This portion of the stripped material lies at less than 2 kpc from the Sun, and it constitutes the nearest stream found in our simulations.

To our knowledge, no extra-tidal structure has been discovered around NGC 6121 thus far. Recently, Kundu, Minniti, and Singh (2019) used RR-Lyrae stars to trace the extra-tidal material around NGC 6121, without finding any clear evidence of structures. The current position of the cluster in the sky, at a latitude of roughly \(20^\circ\) and at a longitude close to \(0^\circ\), makes this search difficult due to the strong contamination of field disk stars, despite the fact that this portion of the stream is expected to be very close to the Sun.

Pal 10 and Pal 2:

Pal 10 and Pal 2 are two disk clusters whose orbit crosses the solar radius. While for Pal 10, the maximal in-plane distance, \(R_{\rm max}\), is approximately \(12\) kpc, in the case of Pal 2, the orbit can reach about 40 kpc from the Galactic center. The fact that both these clusters have a radial excursion of the orbit which is beyond the solar radius implies that their stripped stars can redistribute over the whole longitudinal range and, thus, is also in the anti-center direction. The fact that both clusters have orbits confined close to the disk plane implies that the escaped material redistributes in very thin structures (i.e., confined in a limited latitude interval), which resemble typical “ribbons” in the sky. Our models predict that both clusters are surrounded by a long stream of tidal material, which is however probably very difficult to identify because in both cases, these extra-tidal stars are confined close to the Galactic plane. No tidal streams emanating from these two clusters, to our knowledge, have been identified in the observational data so far. Because the tidal structures associated with these two clusters have similar properties, we report only the case of Pal 10 in Fig. [pal10_stream].

Extra-tidal features originating from inner clusters: Bow-ties and more complex shapes

We have defined inner globular clusters as systems that are not disk clusters (their orbit is not confined close to the Galactic plane) but that are confined inside the solar radius. Seventy-one clusters are found in this category (see Table [classification]). We discuss some of them in the following.

NGC 5946 & NGC 5986:

These are inner-non disk clusters whose escaped stars redistribute in a characteristic “bow-tie” shape. These stars are all confined in a relatively narrow longitudinal range (typically within \(-30^\circ\) to \(30^\circ\)). Towards the edges of the longitude interval, the distribution of extra-tidal stars tends to flare, whereas it instead shrinks at zero longitude. These trends can be explained as an effect of the projection of the orbits of these clusters in the \((\ell, b)\) plane. Moreover, because these clusters always stay in the inner region of the Galaxy, where the dynamical timescales are short, their orbit - and consequently their stripped stars - can experience many disk crossings over the whole duration of the simulation, filling the whole \((\ell, b)\) space allowed by their orbital parameters. An example of such a distribution is given in Fig. [ngc5986_stream] for the extra-tidal material associated with the cluster NGC 5986.

There are also clusters such as NGC 104 (47 Tuc) for which the morphology of the extra-tidal material can take on more complex shapes. This cluster has an orbit confined inside the solar radius, but which, at its apocenter, can reach a distance of less than 1 kpc from the Sun (see Fig. [ngc104_stream]). The projection of the past and future orbit in the (\(\ell, b\)) plane gives rise to a kind of figure-eight shape, with the stars stripped more recently from the cluster tracing the portion of this shape that extends to negative latitudes. Our model of NGC 104 is in agreement with the conclusions reached by Lane, Küpper, and Heggie (2012), who suggested that two clear tidal tails should emanate from this cluster, however, the orientation of these tails in the (\(\ell, b\)) plane differs slightly between the two works. The models of Lane, Küpper, and Heggie (2012) do indeed suggest that the leading tail of NGC 104 should extend a bit beyond \(\ell < -70^\circ\) (see Fig. 4 in their paper, as an example), while our model predicts a minimum longitude of about \(-60^\circ\). The exact comparison between these two works is however difficult, since the model of the Galactic potential, the distance to the Sun, proper motions, and line-of-sight velocities of NGC 104 used in their study are different from ours. To date, clear tails around NGC 104 have not been found yet, with the most recent observational works pointing to the possibility of the presence of a diffuse extended halo-like structure around this cluster (see Piatti 2017). We note that we only find a more diffuse distribution of extra-tidal material in the case of the PII-0.3-SLOW model. Additional work for comparing the current observational data with simulations will be needed to resolve this apparent discrepancy between theoretical predictions and observational findings.

Other shapes found in this category include “Easter eggs”, which are generated by clusters whose orbits are confined to the innermost kpc of the Galaxy and show significant variations in the z-coordinates that are at least as large as those found in the radial direction. Among clusters whose extra-tidal material shows these peculiar shapes we have found: HP1, NGC 6093, NGC 6273, NGC 6293, NGC 6723, and NGC 6809.

Extra-tidal features originating from outer clusters: “canonical” tidal tails

Caveat: In this group, there are some elongated streams emanating from such clusters as NGC 6715, Pal 12, Ter 7, and Ter 8, which are associated with the Sagittarius dwarf galaxy. As a caution, we add that for these clusters the elongation and shape of the streams may be severely modified if the gravitational potential generated by the Sagittarius dwarf galaxy itself was included in the model. For completeness, we decided to include these streams in the paper, and to report them in Supplemental Material 5.2. In future works, we plan to investigate how the inclusion of the Sagittarius dwarf may alter these streams and possibly affect also those of other clusters that are not necessarily associated with this dwarf galaxy.
Among the extra-tidal structures emanating from outer globular clusters, we find some of the most beautiful and elongated tidal tails, of which those associated with Pal 5 were the first to be discovered (Odenkirchen et al. 2001). In this category, we note the stream associated with the E 3 cluster that extends about \(120^\circ\) in longitude based on our models prediction; the thin stream emanating from IC 4499, which we predict has an extension of about \(150^\circ\) in longitude. To recap, the finest and thinnest stellar streams have been predicted from: AM 4, Arp 2, IC 4499, NGC 1261, NGC 3201, NGC 4590, NGC 5024, NGC 5053, NGC 5272, NGC 5466, NGC 5694, NGC 5824, NGC 5904, NGC 6101, NGC 6426, NGC 6584, NGC 6934, Pal 1, Pal 5, Pyxis, Rup 106, Sagittarius II, Ter 7, Ter 8, and Whiting 1.

Many of the above-cited streams have been discovered and also found in observational data, but in many cases, the extent of the tails, as predicted by our models, is larger. Of these observed streams, many tracks are available in the galstreams (Mateu 2023) library. In the following, we compare our model predictions to some of these tracks in the three different Galactic potentials adopted in this paper. More specifically, we compare the projected density distribution of the simulated stream to observations in the \((\ell, b)\), \((\ell, \rm \mu_{\ell}cos(\mathit{b}))\), \((\ell, \rm \mu_b),\) and \((\ell, \rm D)\) planes. As for the projected density distributions derived from the models, we calculated them by taking into account all the 51 simulations realized for each cluster, that is, both the reference simulation and those realized by a Monte-Carlo sampling of the uncertainties. Overall, the agreement between models and observational data is excellent for all Galactic models used. By looking at more extended regions in the sky than those covered by current observations, it should be possible to better constrain the streams, as well as favor or otherwise disfavor some of these models.

NGC 3201:

NGC 3201 is an outer cluster at a distance of about 4.7 kpc from the Sun. It has received a great deal of attention in the last couple of years, followng the suggestion by Riley and Strigari (2020) that part of its tidal tails could be associated with the Gjöll stream, discovered by Rodrigo A. Ibata, Malhan, and Martin (2019). The association of the Gjöll and NGC 3201 streams has been further confirmed by Hansen et al. (2020) on the basis of the similarity of the chemical abundances of stars in the Gjöll stream and in NGC 3201. More recently, Palau and Miralda-Escudé (2021) conducted an extensive study of the tidal tails emanating from this cluster, discovering a long stream, with an overall length of \(140^\circ\) in the sky. Our models suggest that the tails emanating from NGC 3201 may, in fact, be even more extended than those found by Palau and Miralda-Escudé (2021). To further illustrate this point, in Fig. [NGC3201_comp] we compared our model predictions to the tracks available for this stream in the galstreams library as well as those taken from R. Ibata et al. (2021), from Palau and Miralda-Escudé (2021), and from the Gjöll stream, as reported by R. Ibata et al. (2021). All models represent the portion of the stream discovered so far very well, specifically, in terms of the distribution of the stream in the sky, distance, and proper motion spaces. Our models indeed predict that for this cluster, the differences in the stream properties change very little with the mass model adopted. The most striking difference is found for the elongation of the stream at \(\ell > 0\): in the case of the barred potential, the stream associated with NGC 3201 does extend to smaller values of \(\ell\) (up to \(\ell \sim 70^\circ\)), which are not reached in the case of the axisymmetric models. This portion of the stream is expected to be found at distances greater than 20 kpc from the Sun (see bottom panels in Fig. [NGC3201_comp]).

NGC 4590:

NGC 4590 (M 68) is an outer cluster at a current distance of about 10 kpc from the Sun. This cluster is surrounded by a very extended stream (Palau and Miralda-Escudé 2019; R. Ibata et al. 2021), a long portion of which is represented by the Fjörm stream, discovered by R. Ibata et al. (2021). The comparison between our model predictions and the tracks available in galstreams is shown in Fig. [NGC4590_comp]. Interestingly (and differently from the case of NGC 3201), not all the Galactic potentials adopted in this paper seem to represent the stream distribution in the sky equally well. While the axisymmetric models PI and PII predict generally a good match – with model PI describing the stream at positive longitudes even more accurately than model PII – the model PII-0.3-SLOW fails to reproduce the stream in its observed extension: the modeled stream in this case appears quite thick in the (\(\ell, b\)) plane and, moreover, it is much shorter than the stream found in the observational data. Interestingly, the axisymmetric models capture very well also the proper motions and distance-to-the-Sun trends as a function of longitude, as found by R. Ibata et al. (2021) (for Fjörm) and by Palau and Miralda-Escudé (2019), while the NGC 4590 as reported by R. Ibata et al. (2021) (and named “M68-I21” in Fig. [NGC4590_comp]) tends to be off in all models – and also off when compared to the other observational tracks. As suggested by Mateu (2023), it is possible that the stellar stream that was associated with NGC 4590 by R. Ibata et al. (2021)does indeed have a different progenitor. Finally, the failure of the barred model to reproduce the extension of the NGC 4590 stream (and in particular Fjörm) might be due to the choice of the pattern speed adopted or that of the bar length. We will explore these topics in future work. Here, we simply note that given the sensitivity of NGC 4590 stream to the choice of the barred potential, this stream is potentially very interesting for determining the parameters of the latter.

NGC 5466:

NGC 5466 is another cluster known to be surrounded by a thin and very extended stream, as shown by C. J. Grillmair and Johnson (2006) and Belokurov et al. (2006). More recently, Jensen et al. (2021) used Gaia DR2 data to study the stream, finding an extension of about 30 degrees on the sky and a somewhat different orientation than that suggested by C. J. Grillmair and Johnson (2006). Also, R. Ibata et al. (2021) confirmed the existence of an elongated stream around this cluster, even if it appears somewhat less extended than what was found by C. J. Grillmair and Johnson (2006). In Fig. [NGC5466_comp], we show the comparison of our model predictions to the tracks available in galstreams, taken from the works by C. J. Grillmair and Johnson (2006) and R. Ibata et al. (2021). The comparison is excellent for all the Galactic potentials used in this paper. We note that there is a slight offset between the modeled streams and the track by R. Ibata et al. (2021) in the (\(\ell-D\)) plane (of about 1 kpc at a fixed longitude) – which is maybe less evident for the case of the PII-0.3-SLOW potential than in the axisymmetric cases. We emphasize that all our models predict that the stream should be more extended than that discovered so far in the observations. In particular, there is a portion of the leading tail at \(0 < \ell < 42^\circ\) that has not yet been discovered, lying at distances from the Sun closer or similar to that of the known stream.

Pal 5:

About 20 years after the discovery of its tidal tails (Odenkirchen et al. 2001, 2003), Pal 5 still represents the prototype cluster surrounded by thin and extended streams of stars. The extension, morphology, kinematics, and chemical composition of its tails have been extensively studied, both observationally and numerically (Odenkirchen et al. 2003, 2009; Rockosi et al. 2002; Dehnen et al. 2004; Koch et al. 2004; C. J. Grillmair and Dionatos 2006a; Mastrobuono-Battisti et al. 2012; Küpper et al. 2015; Kuzma et al. 2015, 2015; Fritz and Kallivayalil 2015; Rodrigo A. Ibata, Lewis, and Martin 2016; Ishigaki et al. 2016; Guillaume F. Thomas et al. 2016; Koch and Côté 2017; Rodrigo A. Ibata et al. 2017; Pearson, Price-Whelan, and Johnston 2017; Price-Whelan et al. 2019; Starkman, Bovy, and Webb 2020; Ana Bonaca, Pearson, et al. 2020; R. Ibata et al. 2021; Phillips et al. 2022). In Fig. [Pal5_comp], we report the comparison of our model predictions to the tracks available for this cluster in galstreams and those coming from the work by Price-Whelan et al. (2019), Starkman, Bovy, and Webb (2020), and R. Ibata et al. (2021). Projected in the \((\ell, b)\) plane, the observed streams are in excellent agreement with the model predictions for all the Galactic potentials adopted. Interestingly, all potentials suggest more extended streams than those discovered so far. In the case of the barred potential, we note that the prediction of the stream position in the sky at large angular distances from the cluster center is still highly uncertain. This is due to the fact that the uncertainties are still affecting the current distance of Pal 5 to the Sun, combined with the impact of the rotating bar, which can be more or less efficient in perturbing the stream depending on the torques experienced by the latter at pericenter; these, in turn, depend on the orbital parameters of the cluster itself. Pearson, Price-Whelan, and Johnston (2017) already explored the effect of a rotating bar on the extension and morphology of Pal 5 tails. While we can confirm the Pearson, Price-Whelan, and Johnston (2017) findings, namely, that both characteristics are affected by a rotating bar, we also emphasize that both the extension of the leading tail (i.e., the portion of the stream at negative longitudes) and its morphology depend on the choice of the bar pattern speed. For example, we do not find a density drop in the leading tail as reported by Pearson, Price-Whelan, and Johnston (2017). These authors adopted a higher bar pattern speed than the one adopted in this work (\(\Omega_b=60\) kms\(^{-1}\)kpc\(^{-1}\) for the example discussed in Fig. 2 of their work, right panel, while \(\Omega_b=38\) kms\(^{-1}\)kpc\(^{-1}\) in our PII-0.3-SLOW model). Additionally, taking into account the uncertainties in the cluster distance, proper motions, and line-of-sight velocity is also important for obtaining robust predictions on the stream characteristics. Some of our solutions, for example, predict a very extended leading tail that is significantly more extended than those found for the axisymmetric potentials.

Conclusion

In this work, we have presented the first simulated catalogue of all Galactic globular clusters for which 6D phase-space information, along with masses and sizes are available. A total of 159 globular clusters has been simulated in three Milky Way-like potentials, modeling the process of tidal stripping that these clusters have experienced over the past 5 Gyr. As a result, for all clusters, we can predict the distribution of the extra-tidal material in the sky, their proper motions, distances to the Sun, and line-of-sight velocities. Errors on 6D phase-space information have been taken into account by generating 50 complementary simulations for each cluster, with a Monte-Carlo sampling of the uncertainties. This catalogue currently contains 24 327 simulations, for a total volume of about 370 Gb. It will be made publicly available⁹, with the intention to provide the community with an instrument that allows for: a more complete view of the expected distribution of globular clusters tidal structures in the sky, helping in formulating the interpretation of recent and future discoveries, supporting the search for new extra-tidal features in the data, and offering the community a repository of all these models to be compared to other theoretical and numerical predictions, which employ different Galactic potentials and/or gravity laws.

In this first paper, we have presented the distribution in the \((\ell, b)\) plane of all the simulated extra-tidal features. A striking result is the variety of extra-tidal shapes that globular clusters can give rise to. The canonical tidal tails à la Palomar 5 are only some of the multiple morphologies these structures can have. Ribbons, bow-ties, padlocks, and halo-shapes are also common. This variety of shapes that these stripped stars can show in the sky depends on the characteristics of the cluster orbit, on its current position in the orbit itself, and on the distance of the extra-tidal material from the Sun. Any search for the left-overs of globular clusters in the field should take into account this richness of distributions and morphologies.

These simulations have also allowed us to derive an estimate of the expected mass of stars that have escaped from the clusters in the last 5 Gyr and now in the field. Although these are approximations, given the limitation of our models, our estimate of the mass lost from the Galactic globular cluster system over the last 5 Gyr is between \(2-21\times 10^6 M_{\odot}\), which is comparable to up to half of the total stellar mass found nowadays in the Galactic globular cluster system itself.

This work is intended to be the first of a series which investigates the properties of globular clusters streams in a variety of realistic Galactic potentials, including the perturbations induced by close dwarf satellites (Sagittarius and the Large and Small Magellanic Clouds), as well as more complex and time-varying distributions for the dark matter component.

Supplementary material

Choice of the time-step for orbit integration

To choose the optimal time-step \(\Delta t\) for the simulations, we quantified the energy conservation of the orbit integration, in the case where the 159 clusters evolve in isolation, that is, the case where each test-particle in a cluster feels the gravitational attraction of the cluster itself – but not that of the Galaxy. Since in the case of an isolated cluster, the gravitational potential is time-independent, the total energy, \(E_i\), of each particle in the system (sum of the particle kinetic and potential energy) must be conserved. Any departure from energy conservation is thus a check of the quality of the integration and in particular of the choice of the adopted time-step.

For the isolated simulations, we thus proceeded as follows. We modeled each cluster as a set of \(N\)=100 000 test-particles, subject to the cluster potential, only. To model this latter, for each cluster we adopted the same Plummer sphere distribution, with same characteristic radius \(r_c\) and total masse \(M_{GC}\), as those adopted for the simulations described in Sect. 2 (see in particular Sect. 2.2). A first set of 159 isolated simulations was run adopting a \(\Delta t=10^5\) yr and a total number of steps \(N_{steps}=50 000\) for all clusters, for a total simulated time interval of 5 Gyr. We then quantified the energy conservation by calculating the median error per time-step, \(m_{err}\), of \((\Delta E/E)_i=|(E_{fin,i}-E_{ini,i})/E_{ini,i}|\), where \(E_{ini,i}\) and \(E_{fin,i}\) are, respectively, the initial (time \(t=0\)) and final (time t=5 Gyr) energy of each particle in the system.

While for most of the simulated clusters (109 over 159), this choice of the time-step was sufficient to guarantee an excellent energy conservation (with \(m_{err}\) typically of the order of \(10^{-11}-10^{-12}\)) for 50 clusters the corresponding \(m_{err}\) values were found to be above \(10^{-10}\). This was the case for all clusters with crossing times \(t_{cross}=\sqrt{{r_c}^3/GM_{GC}}\) lower than \(2\times 10^5\) yr, that is about twice the time-step. For these clusters, we hence reduced the \(\Delta t\) of a factor 10, rerunning the simulations with \(\Delta t=10^4\) yr and \(N_{steps}=500 000\). With such a choice, the corresponding energy conservation turned out to be excellent (below \(10^{-10}\) per step). In Table [tcross-energy], we summarize the result of this study, reporting the cluster name, crossing time, and median error, \(m_{err}\), in energy conservation obtained for all isolated cluster simulations. Clusters for which a \(\Delta t=10^4\) yr, and a corresponding number \(N_{steps}\) have been used, are indicated in the Table with an asterisk. The values of \(\Delta t\) adopted for the isolated simulations, and the associated number of time steps, \(N_{steps}\), are also those used to run the simulations of the same clusters orbiting in the gravitational field of the Milky Way.

Extra-tidal features generated by all the simulated clusters

In this section, we report all the extra-tidal features as predicted by our models. For each cluster, we show (from Figs. [stream1] to [stream20]) the probability density of finding associated extra-tidal features in the sky, by calculating the 2D histogram of the escaped particles. Each particle per cluster is present, meaning that all 100,000 particles originating each of the 50 Monte-Carlo realizations are included, plus the case with the best values (see Sect. 2.1). Thus, each map is a 500x500 histogram that bins 5.1\(\times10^6\) particles. The retrieved cumulative 2D histogram is then normalized to its maximum value and plotted in the following figures in logarithmic scale.

Globular cluster initial conditions

Globular cluster classification

In Fig. [orbparam] (top panel), we show the distribution of the arctangent of the \(z_{max}/R_{max}\) ratio, with \(z_{max}\) and \(R_{max}\) being (respectively) the maximum height above or below the Galactic plane and their maximum in-plane distance from the Galactic center that was reached in the past 5 Gyr of orbital evolution in the Galactic potential adopted in this paper (see Sect. 2.3). As was already noticed for field stars (see Haywood et al. 2018), the GC distribution also shows a dip at about \(10^\circ\), which separates clusters with flattened orbits (arctan(\(z_{max}/R_{max}\)) \(\le 10^\circ\)) from thicker ones. We thus define a first set of clusters (the disk GCs) as that containing all globular clusters with arctan(\(z_{max}/R_{max}\)) \(\le 10^\circ\). This first set contains 21 clusters. Of the remaining 138, we distinguish between an “inner” GCs sample and an “outer” GCs sample, on the basis of the maximum 3D distance (\(r_{max}\)) that the cluster reaches from the Galactic center. Inner GCs are those with \(r_{max} \le 10\) kpc and outer GCs are those with \(r_{max} > 8.34\) kpc, which is the value of the distance of the Sun to the Galactic Center used in this experiment. Such a value allows us to discriminate between two classes of tidal debris, for inner clusters are necessarily restricted in latitude and longitude, whereas outer cluster can fill the sky.

Finally, the third and bottom panels of Fig. [orbparam] show the distribution of these three defined groups in the \(R_{max}-z_{max}\) and \(E-L_z\) planes. We note that since disk clusters are uniquely defined on the basis of the ratio between the maximum vertical and in-plane orbital excursion – and not on the circularity of their orbits (as seen, e.g., in Massari, Koppelman, and Helmi 2019), some of our disk GCs have elongated (i.e., radial) orbits (\(L_z \sim 0\)) or even retrograde ones (\(L_z > 0\)). Our definition of disk clusters is purely related to a morphological criterium: disk GCs are those whose orbits are confined close to the Galactic plane, independently of their eccentricity.