DATING THE TIDAL DISRUPTION OF GLOBULAR CLUSTERS WITH GAIA DATA ON THEIR STELLAR STREAMS

Dating the Tidal Disruption of Globular Clusters With Gaia Data on Their Stellar Streams

[    Idan Ginsburg    Abraham Loeb
Abstract

The Gaia mission promises to deliver precision astrometry at an unprecedented level, heralding a new era for discerning the kinematic and spatial coordinates of stars in our Galaxy. Here, we present a new technique for estimating the age of tidally disrupted globular cluster streams using the proper motions and parallaxes of tracer stars. We evolve the collisional dynamics of globular clusters within the evolving potential of a Milky Way-like halo extracted from a cosmological CDM simulation and analyze the resultant streams as they would be observed by Gaia. The simulations sample a variety of globular cluster orbits, and account for stellar evolution and the gravitational influence of the disk of the Milky Way. We show that a characteristic timescale, obtained from the dispersion of the proper motions and parallaxes of stars within the stream, is a good indicator for the time elapsed since the stream has been freely expanding away due to the tidal disruption of the globular cluster. This timescale, in turn, places a lower limit on the age of the cluster. The age can be deduced from astrometry using a modest number of stars, with the error on this estimate depending on the proximity of the stream and the number of tracer stars used.

(Galaxy:) globular clusters: general – proper motions – methods: numerical
Corresponding author: Sownak Bose

0000-0002-0974-5266]Sownak Bose \move@AU\move@AF\@affiliationHarvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

\move@AU\move@AF\@affiliation

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

\move@AU\move@AF\@affiliation

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

1 Introduction

Globular clusters (GCs) are dense, spherical, concentrations of stars with a characteristic mass of (see Brodie & Strader 2006 for a review) and half-light radii of 3–10 pc (e.g. Jordán et al., 2005; van den Bergh, 2008). It is universally agreed that GCs are among the oldest observable objects and are found in numerous galaxies, with giant elliptical galaxies harboring thousands of them (Peng et al., 2011). Surveys have discovered about 150 GCs around the Milky Way (e.g. Kharchenko et al., 2013). Understanding the nature and origin of these GCs has important implications for not only the formation history of the Milky Way, but also for models of structure and galaxy formation (e.g. Kravtsov & Gnedin, 2005; Boylan-Kolchin, 2017).

GCs are grouped into two sub-populations depending on metallicity. So-called “blue clusters” are metal-poor whereas, “red clusters” are metal-rich (Zinn & West, 1984; Usher et al., 2012; Roediger et al., 2014). Renaud et al. (2017) argued that blue clusters form in satellite galaxies and are accreted onto the Milky Way, whereas red clusters form in situ. Furthermore, Kundu & Whitmore (1998) noted that blue GCs appear to be 20% larger than their redder counterparts. However, it is unclear whether this is a real physical phenomenon or a projection effect (Larsen & Brodie, 2003).

GC ages are typically acquired from studies of color-magnitude diagrams in conjunction with stellar evolution models (Forbes & Bridges, 2010; Correnti et al., 2016; Powalka et al., 2017; Kerber et al., 2018). Such studies show that Galactic GCs (GGCs) typically have ages over 10 Gyr, often exceeding 12 Gyr. However, there are also young GGCs such as Whiting 1, which has an estimated age of Gyr (Carraro, 2005; Valcheva et al., 2015). As they traverse the Milky Way’s potential, GGCs evaporate over time and leave extended tidal tails (Gnedin & Ostriker, 1997; Fall & Zhang, 2001; Carlberg, 2017).

In this Letter we explore a novel technique to date the stream associated with a tidally disrupted GC. This, in turn, can be used to constrain the gravitational potential and the assembly history of the Milky Way (e.g. Johnston et al., 1999; Pearson et al., 2015; Bonaca & Hogg, 2018). Our method relies on measurements of positions and motions of stars in the plane of the sky, which will become available soon with the Gaia DR2 catalog. Gaia is expected to provide precise astrometric measurements for over a billion stars in the Milky Way, with proper motions and parallaxes determined to 1 percent or better up to 15 kpc (Pancino et al., 2017).

The remainder of this Letter is organized as follows: in Section 2, we discuss the numerical techniques used and the parameter space that we explore. In Section 3 we visualize the simulated streams as they might be detected by Gaia, and in Section 4, we discuss our method for age determination. Finally, we summarize our conclusions in Section 5.

2 Numerical methods

In the following subsections, we outline the details of our numerical setup.

2.1 The galactic potential of the Milky Way

We consider the combined effect of the time-evolving potential of the smooth dark matter halo component of a Milky Way-like galaxy, as well as the contribution of a central disk galaxy. For the halo component, we first create zoom-in initial conditions of a Milky Way-like halo () extracted from the dark matter-only Copernicus Complexio simulations (Bose et al., 2016; Hellwing et al., 2016). Here, is the mass contained within , the radius within which the mean density is 200 times the critical density of the Universe. The halo that we have chosen for re-simulation exhibits a fairly typical accretion history for halos of this mass, and was evolved from redshift to using the P-Gadget3 code (Springel et al., 2008).

At each simulation output, we compute the gravitational potential associated with the high-resolution particle distribution, and find that its evolution with redshift is captured well by an analytic expression of the form (Buist & Helmi, 2014; Renaud & Gieles, 2015):

 ϕ(r,z)=−GMs(z)rln(1+rrs(z)), (1)

where is Newton’s constant, is the scale radius of the halo, while is the mass contained within this radius. The redshift evolution of these parameters can be written in the form:

 Ms(z) = Ms(0)exp(−0.08z), (2) rs(z) = rs(0)exp(−0.05z),

where , kpc, and the coefficients of the exponentials are obtained by fitting the potential from snapshots of our re-simulation.

The contribution of a central disk galaxy is modeled using a superposition of three Miyamoto-Nagai disk potentials (Miyamoto & Nagai, 1975), following the procedure described by Smith et al. (2015). The disk is assigned a total mass of with a scale length of 3 kpc and a scale height of 300 pc. The disk potential does not evolve dynamically.

2.2 Globular cluster dynamics

Initial conditions for the GC simulations were generated using the McLuster code (Küpper et al., 2011). We sample the positions and velocities of star particles in the GCs according to a King profile (King, 1966). In each case, the total initial mass of the star cluster is set to , with a concentration parameter . In five out of the six GC simulations that we run, we set the half-mass radius pc, which is a relatively large value for typical GCs. In a final simulation, therefore, we additionally include an example with pc. Masses of individual stars are assigned by sampling a Salpeter (1955)-like initial mass function (IMF). For computational reasons, we truncate the low-mass end of the IMF to , resulting in a total of 21,730 star particles in the cluster.

The dynamical evolution of these star clusters is performed using the publicly available Nbody6++ code (Wang et al., 2015). Nbody6++ is a hybrid MPI-GPU accelerated version of the Nbody6 code (Aarseth, 2003; Nitadori & Aarseth, 2012), and contains routines accounting for the evolution of star clusters in background tidal fields developed as part of Nbody6-tt (Renaud & Gieles, 2015). The simulations (and, therefore, the halo potential in Eq. 1) are initialized at using the same random seed and evolved for 12.96 Gyr, corresponding to the time interval until . Mass loss through stellar evolution is included in all of our simulations. The initial configurations for all the simulations that we have performed are summarized in Table 2.2, where we have also listed the orientations of the position/velocity vectors with which the GCs are initialized. Note that the parameter space that we explore in our simulations is not designed to be representative of the population of GCs in the Milky Way.

Figure 2.2 shows projections of a subset of our simulations at the present time. As the initial GC structure is identical in each run, different stream morphologies can be attributed to the different choice of initial configurations as listed in Table 2.2. For example, while the cluster is initialized with the same orbital parameters in GC-2 (top row) and GC-6 (bottom row), the more compact initial structure of GC-6 ( pc) compared to GC-2 ( pc) means that the former is less susceptible to tidal disruption, resulting in a more prominent cluster core at . On the other hand, while GC-5 (middle row) shares identical structural properties to GC-2, it also experiences less tidal disruption for the simple reason that it undergoes fewer pericentric passages than GC-2 (see, e.g. Peñarrubia et al., 2009).

3 Simulated streams as observed by gaia

In what follows, we emulate mock Gaia observations of our simulated GC streams by projecting the stellar distribution onto all-sky maps. We convert between cartesian coordinates used in the simulation to heliocentric latitude, , and longitude, , using the astropy package (The Astropy Collaboration et al., 2018). For each star, we use its mass to calculate its luminosity, and use its distance from the observer to compute the equivalent apparent magnitude, . Gaia is expected to be complete down to a nominal magnitude limit of  mag; this threshold can be used to determine which sections of each simulated stream would be within Gaia’s detection capability.

Figure 2.2 displays all-sky projections of streams from three of our simulations (GC-2, GC-4 and GC-5) as they would be seen by an observer on the solar circle around the center of the Milky Way. Stars that would be observed by Gaia are shown in color; the range of colors corresponds to the radial velocity, , of the star relative to the observer. Unobserved portions of the stream (i.e., composed of stars fainter than  mag) are shown in gray. In each panel, we determine an ‘optimal’ observer position as the point on the solar circle that maximizes the number of observed stars.

As expected, the number of stars observed at present depends on the initial orbital configuration for each simulation. For example, both GC-2 and GC-4 are initialized 35 kpc from the center of the Milky Way, but owing to its lower initial orbital velocity (see Table 2.2), GC-4’s motion is more tightly bound by the disk potential, ending up at a galactocentric distance of  kpc at (rather than  kpc in the case of GC-2). The result is that almost three times as many stars in GC-4 fall within Gaia’s observable window compared to GC-2. Similarly, GC-5, which started off 60 kpc from the center of the Milky Way, displays the smallest number of observed stars. In the remainder of this Letter, we will be concerned only with the set of observed stars in each simulation.

4 A characteristic timescale using proper motions

Next, we demonstrate that a characteristic timescale that can be defined using the proper motions of tracer stars in a stream matches remarkably well the period that these stars have spent outside the tidal influence of the GC from which they originated.

Given a set of stars observed in a stream at present day, we can label their positions (velocities) in two orthogonal directions along the plane of the sky as and ( and ). The subscripts on each quantity refer to the fact that they are measured along the width and along the length of the stream – corresponding to measurements along the and axes, respectively, as shown in the upper panel of Figure 3. These quantities are directly related to the parallax and proper motion in and coordinates using the relations:

 Xw =r⊙−robscosbcosl, (3) Xl =robssinb, Vw =V⊙,x+Vrcosbcosl −robs[μbcoslsinb+μlcosbsinl], Vl =V⊙,z+Vrsinb+robsμbcosb,

where is the radial distance to the star from the observer, and , respectively, are the proper motions in the and directions, kpc is the distance to the Galactic center and kms is the solar motion relative to it (Schönrich et al., 2010; Bovy, 2015). For these tracers, we can then define the dispersion in position, , and velocity, , along the stream as:

 σx = √σ2Xw+σ2Xl, (4) σv = √σ2Vw+σ2Vl,

Finally, we define a characteristic timescale ascertained from these proper motions, , given by:

 tPM=σxσv. (5)

Once stars have escaped the tidal influence of the GC, their motions should be dominated by acceleration due to the external galactic potential these stars are embedded in, rather than the GC itself. For an external potential that evolves reasonably slowly in time, the timescale provides an estimate for the duration of the period that the escaped stars have been under the influence of the external potential. Consequently, if the age of a stream is defined as the time since the disruption of the GC resulting in the stream, the timescale defined by Eq. (5) should be correlated with the age of the stream itself.

To investigate whether this is indeed the case, we followed the evolution of the tracer stars backward in time in each simulation to find out when their dynamics are no longer dominated by the GC. This regime can be identified using the tidal radius, , which sets the radial distance from the center of the GC at which the potential of the cluster is balanced by the background potential. Theoretically, we estimate the tidal radius as a function of redshift as:

 r3t(z)≈GMc(z)V2circ(z)r2gc(z), (6)

(c.f., Dehnen et al., 2004) where is the circular velocity of the host halo (measured at ) at redshift, , while and , respectively, are the mass and galactocentric distance of the star cluster at this redshift. Each of these quantities are estimated at the output times of our simulations. Over the course of the simulations, increases while decreases, resulting in an overall decrease in as a function of time. We consider two possible criteria for a star to be within the tidal influence of the cluster: (i) and (ii) , where is the distance of a given star particle from the center of the cluster. More specifically, we determine the age of the stream using the two tidal radii criteria, , as the last epoch when at least 50% of the all stars in the stream (i.e., observed and unobserved) satisfy conditions (i) or (ii). Condition (i) is more commonly adopted in the literature (e.g. Wilkinson et al., 2003; Aarseth, 2012; Madrid et al., 2014); additionally, considering the more restrictive condition (ii) gives a handle of the uncertainty in measuring the time since tidal escape. By comparing with , we can test how accurately Eq. (5) can be used to determine the time since the disruption of the GC.

Figure 3 is an illustration of this comparison for the GC-1 simulation. The stars in the stream marked in orange (158 in total) are those that could be observed by Gaia, and act as the tracer population for estimating the age of the stream. Using Eqs. (4) & (5), we estimate the age to be Gyr, corresponding to a lookback time of Gyr. To estimate the “dynamical” age of the stream, we then trace the evolution of the entire set of stars (tracers and unobserved) through the GC-1 simulation, and find that at least half of this population satisfies condition (i) at Gyr and condition (ii) at Gyr.

The result of this analysis for our other simulations is shown in Figure 4. We observe streams in each simulation at two output times. In some cases, the same part of the stream is observed twice; these correspond to data points that have the same value of . Given the apparent magnitude of each tracer star in the simulated stream, we estimate errors on its proper motion and parallax as they would be measured by Gaia using the PyGaia package; these errors are combined to get a rough estimate of the error on . The error on any given measurement therefore depends on both the size of the tracer population (which we indicate using the color scale of each data point in Figure 4) and the brightness of the stars in this set (or, equivalently, the proximity of the stream).

In general, the age of the stream estimated using proper motions matches quite well the lookback at which these stars escaped the tidal radius of the GC. The agreement between and is typically within 20-40% in the majority of the experiments that we have carried out. Encouragingly, this level of agreement is also true for GC-6 ( pc), for which the tidal disruption rate (and, consequently, the distribution of energy and angular momenta of the disrupted stars) is different to the other examples that we have considered.

Figure 4 shows that underpredicts the “true” age of the stream, . The reason for this systematic difference can be ascribed to the population of stars chosen to trace the age using proper motions. Limiting the estimate of to only “observed” stars preferentially selects more massive stream members. Due to mass segregation within a GC, massive stars tend to sink toward the cluster center while low-mass stars move farther away – these low-mass stars are the first to exit the tidal radius of the cluster upon disruption. Using only the more massive, observed stars slightly underestimates the actual time since disruption, and should therefore be interpreted as a lower bound on the age of the stream.

It is worth highlighting that Eq. (5) does not necessarily reflect the age of the GC itself; instead, it is a measure of the time at which the tidal disruption of the globular cluster resulted in a given part of the observed stream. An observational determination of this time can be used to constrain the assembly history of the Milky Way.

5 Conclusions

By virtue of their ancient stellar populations, GCs represent the most interesting astrophysical entities found in galaxies. As these dense concentrations of stars orbit their host galaxy, they are continually stripped by the tidal field of their host, resulting in the formation of cold, extended stellar streams.

The Gaia mission will provide precise measurements of the parallax and proper motions for more than a billion stars in the Milky Way, paving the way for a deeper understanding of the kinematics of stars and the assembly history of our Galaxy. In this Letter we have shown that, using a sample of tracer stars in a stellar stream, it is possible to use the proper motions of these stars to infer the epoch at which a globular cluster was disrupted, resulting in the formation of the stream. Specifically, we find that a timescale defined using the dispersion in positions and velocities of stars in the plane of the stream (Eq. 5) provides a good estimator for how long these stars have spent outside the tidal radius of the cluster.

We verified our results by running a sequence of simulations evolving a cluster of mass in a time-evolving potential comprising the dark matter halo of the Milky Way (extracted from a cosmological zoom simulation) and a (static) central disk (Section 2). As shown in Figure 4, the age of the stream inferred using astrometric information of tracer stars at matches very well the epoch at which these stars were last contained within the tidal radius of the globular cluster. The procedure outlined in this Letter can therefore be used as an effective method for dating the tidal disruption of a globular cluster, which in turn serves as a lower bound on the age of the cluster itself.

In reality, correctly identifying stream members in the observed data is more challenging than simply assuming a magnitude cut as we have done in this work. In particular, misidentifying stream members can result in incorrectly measured dispersions, which would subsequently propagate as a large error in the inferred age of the stream. Radial velocities, where available, can be useful: GC streams are typically very cold, and “true” stream members are likely to be clustered in a phase diagram of radial distance and line-of-sight velocity. Jointly considering both the proper motions and color-magnitude diagram can also be informative in faithfully separating members of a stellar stream from foreground/background contaminants. The Gaia DR2 dataset will be particularly transformative in this regard.

We are thankful to the referee and the editor for making a number of suggestions that have significantly improved the overall quality of this manuscript. We are grateful to Long Wang for providing extensive support in the setup and use of Nbody6++, and to Adrian Jenkins for giving us access to his codes for generating cosmological initial conditions. We benefited from some useful conversations with Ana Bonaca. We are thankful to the developers of astropy and PyGaia for maintaining these codes and making them public. S.B. is supported by Harvard University through the ITC Fellowship. S.B. also acknowledges the hospitality provided by the Kavli Institute for Theoretical Physics in Santa Barbara during the “Small-Scale Structure of Cold(?) Dark Matter” programme, where part of this work was completed. This research was supported in part by the National Science Foundation under grant No. NSF PHY-1748958, and in part by the Black Hole Initiative, which is funded by a grant from the John Templeton Foundation.

References

• Aarseth (2003) Aarseth, S. J. 2003, Gravitational N-Body Simulations, by Sverre J. Aarseth, pp. 430. ISBN 0521432723. Cambridge, UK: Cambridge University Press, November 2003., 430
• Aarseth (2012) Aarseth, S. J. 2012, MNRAS, 422, 841
• Bonaca & Hogg (2018) Bonaca, A., & Hogg, D. W. 2018, arXiv:1804.06854
• Bose et al. (2016) Bose, S., Hellwing, W. A., Frenk, C. S., et al. 2016, MNRAS, 455, 318
• Bovy (2015) Bovy, J. 2015, ApJS, 216, 29
• Boylan-Kolchin (2017) Boylan-Kolchin, M. 2017, arXiv:1711.00009
• Brodie & Strader (2006) Brodie, J. P., & Strader, J. 2006, ARA&A, 44, 193
• Buist & Helmi (2014) Buist, H. J. T., & Helmi, A. 2014, A&A, 563, A110
• Carlberg (2017) Carlberg, R. G. 2017, ApJ, 838, 39
• Carraro (2005) Carraro, G. 2005, ApJL, 621, L61
• Correnti et al. (2016) Correnti, M., Gennaro, M., Kalirai, J. S., Brown, T. M., & Calamida, A. 2016, ApJ, 823, 18
• Dehnen et al. (2004) Dehnen, W., Odenkirchen, M., Grebel, E. K., & Rix, H.-W. 2004, AJ, 127, 2753
• Fall & Zhang (2001) Fall, S. M., & Zhang, Q. 2001, ApJ, 561, 751
• Forbes & Bridges (2010) Forbes, D. A., & Bridges, T. 2010, MNRAS, 404, 1203
• Gnedin & Ostriker (1997) Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223
• Hellwing et al. (2016) Hellwing, W. A., Frenk, C. S., Cautun, M., et al. 2016, MNRAS, 457, 3492
• Johnston et al. (1999) Johnston, K. V., Zhao, H., Spergel, D. N., & Hernquist, L. 1999, ApJL, 512, L109
• Jordán et al. (2005) Jordán, A., Côté, P., Blakeslee, J. P., et al. 2005, ApJ, 634, 1002
• Kerber et al. (2018) Kerber, L. O., Nardiello, D., Ortolani, S., et al. 2018, ApJ, 853, 15
• Kharchenko et al. (2013) Kharchenko, N. V., Piskunov, A. E., Schilbach, E., Röser, S., & Scholz, R.-D. 2013, A&A, 558, A53
• King (1966) King, I. R. 1966, AJ, 71, 64
• Kravtsov & Gnedin (2005) Kravtsov, A. V., & Gnedin, O. Y. 2005, ApJ, 623, 650
• Kundu & Whitmore (1998) Kundu, A., & Whitmore, B. C. 1998, AJ, 116, 2841
• Küpper et al. (2011) Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300
• Larsen & Brodie (2003) Larsen, S. S., & Brodie, J. P. 2003, ApJ, 593, 340
• Madrid et al. (2014) Madrid, J. P., Hurley, J. R., & Martig, M. 2014, ApJ, 784, 95
• Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
• Nitadori & Aarseth (2012) Nitadori, K., & Aarseth, S. J. 2012, MNRAS, 424, 545
• Pancino et al. (2017) Pancino, E., Bellazzini, M., Giuffrida, G., & Marinoni, S. 2017, MNRAS, 467, 412
• Peñarrubia et al. (2009) Peñarrubia, J., Walker, M. G., & Gilmore, G. 2009, MNRAS, 399, 1275
• Pearson et al. (2015) Pearson, S., Küpper, A. H. W., Johnston, K. V., & Price-Whelan, A. M. 2015, ApJ, 799, 28
• Peng et al. (2011) Peng, E. W., Ferguson, H. C., Goudfrooij, P., et al. 2011, ApJ, 730, 23
• Powalka et al. (2017) Powalka, M., Lançon, A., Puzia, T. H., et al. 2017, ApJ, 844, 104
• Renaud et al. (2017) Renaud, F., Agertz, O., & Gieles, M. 2017, MNRAS, 465, 3622
• Renaud & Gieles (2015) Renaud, F., & Gieles, M. 2015, MNRAS, 449, 2734
• Renaud & Gieles (2015) Renaud, F., & Gieles, M. 2015, MNRAS, 448, 3416
• Roediger et al. (2014) Roediger, J. C., Courteau, S., Graves, G., & Schiavon, R. P. 2014, ApJS, 210, 10
• Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
• Schönrich et al. (2010) Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829
• Smith et al. (2015) Smith, R., Flynn, C., Candlish, G. N., Fellhauer, M., & Gibson, B. K. 2015, MNRAS, 448, 2934
• Springel et al. (2008) Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685
• The Astropy Collaboration et al. (2018) The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, arXiv:1801.02634
• Usher et al. (2012) Usher, C., Forbes, D. A., Brodie, J. P., et al. 2012, MNRAS, 426, 1475
• Valcheva et al. (2015) Valcheva, A. T., Ovcharov, E. P., Lalova, A. D., et al. 2015, MNRAS, 446, 730
• van den Bergh (2008) van den Bergh, S. 2008, MNRAS, 385, L20
• Wang et al. (2015) Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070
• Wilkinson et al. (2003) Wilkinson, M. I., Hurley, J. R., Mackey, A. D., Gilmore, G. F., & Tout, C. A. 2003, MNRAS, 343, 1025
• Zinn & West (1984) Zinn, R., & West, M. J. 1984, ApJS, 55, 45
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters