Dating the Tidal Disruption of Globular Clusters with GAIA Data on Their Stellar Streams

Globular clusters (GCs) are dense, spherical, concentrations of stars with a characteristic mass of ∼10⁵ M_⊙ (see Brodie & Strader 2006 for a review) and half-light radii of ${r}_{h}$ ∼ 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 ∼6 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.

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 (M₂₀₀ = 1.5 × 10¹² M_⊙) extracted from the dark matter-only Copernicus Complexio simulations (Bose et al. 2016; Hellwing et al. 2016). Here, M₂₀₀ is the mass contained within r₂₀₀, 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 z = 127 to z = 0 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 2015a):

$$\begin{eqnarray}&&\phi (r,z)=-\displaystyle \frac{{{GM}}_{s}(z)}{r}\mathrm{ln}\left(1+\displaystyle \frac{r}{{r}_{s}(z)}\right),\end{eqnarray} \tag{ 1 }$$

where G is Newton's constant, r_s(z) is the scale radius of the halo, while M_s(z) is the mass contained within this radius. The redshift evolution of these parameters can be written in the form:

$$\begin{eqnarray}{M}_{s}(z) & = & {M}_{s}(0)\exp (-0.08z),\\ {r}_{s}(z) & = & {r}_{s}(0)\exp (-0.05z),\end{eqnarray} \tag{ 2 }$$

where M_s(0) = 2 × 10¹¹ M_⊙, r_s(0) = 20 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 5 × 10¹⁰ M_⊙ 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 10⁵ M_⊙, with a concentration parameter W₀ = 6. In five out of the six GC simulations that we run, we set the half-mass radius r_h = 8 pc, which is a relatively large value for typical GCs. In a final simulation, therefore, we additionally include an example with r_h = 3 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 ∼0.9 M_⊙, resulting in a total of 21,730 star particles in the 10⁵ M_⊙ 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 2015b). The simulations (and, therefore, the halo potential in Equation (1)) are initialized at z = 7 using the same random seed and evolved for 12.96 Gyr, corresponding to the time interval until z = 0. 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 1, 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.

Table 1.  A Summary of the Numerical Experiments Performed

Label	Initial Position	Initial Velocity	Half-light Radius
	(kpc)	(kms−1)	(pc)
GC-1	(20.0, 0.0, 0.0)	(84.7, 0.0, 84.7)	8
GC-2	(35.0, 0.0, 0.0)	(84.7, 0.0, 84.7)	8
GC-3	(60.0, 0.0, 0.0)	(84.7, 0.0, 84.7)	8
GC-4	(35.0, 0.0, 0.0)	(63.0, 0.0, 63.0)	8
GC-5	(60.0, 0.0, 0.0)	(35.0, 0.0, 35.0)	8
GC-6	(35.0, 0.0, 0.0)	(84.7, 0.0, 84.7)	3

Note. All distances and velocities are measured relative to the center of the galactic potential. The GC orbits, which for simplicity we have constrained to lie solely in the x–z plane (i.e., perpendicular to the plane of the galactic disk), exhibit a range of eccentricities sensitive to both the initial position and velocity of the GC. The potential of the dark matter halo and that of the central disk is included in all simulations.

Download table as:  ASCIITypeset image

Figure 1 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 1. 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 (r_h = 3 pc) compared to GC-2 (r_h = 8 pc) means that the former is less susceptible to tidal disruption, resulting in a more prominent cluster core at z = 0. 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).

Zoom In Zoom Out Reset image size

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, ${\mathscr{b}}$, and longitude, ${\mathscr{l}}$, 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, G. Gaia is expected to be complete down to a nominal magnitude limit of G ≈ 20 mag; this threshold can be used to determine which sections of each simulated stream would be within Gaia's detection capability.

Figure 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, V_r, of the star relative to the observer. Unobserved portions of the stream (i.e., composed of stars fainter than G = 20 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.

Zoom In Zoom Out Reset image size

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 1), GC-4's motion is more tightly bound by the disk potential, ending up at a galactocentric distance of ∼12 kpc at z = 0 (rather than ∼20 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.

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 X_w and X_l (V_w and V_l). 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 x and z axes, respectively, as shown in the upper panel of Figure 3. These quantities are directly related to the parallax and proper motion in ${\mathscr{l}}$ and ${\mathscr{b}}$ coordinates using the relations

$$\begin{eqnarray}{X}_{w} & = & {r}_{\odot }-{r}_{\mathrm{obs}}\cos {\mathscr{b}}\cos {\mathscr{l}},\\ {X}_{l} & = & {r}_{\mathrm{obs}}\sin {\mathscr{b}},\\ {V}_{w} & = & {V}_{\odot ,x}+{V}_{r}\cos {\mathscr{b}}\cos {\mathscr{l}}\\ & & -{r}_{\mathrm{obs}}[{\mu }_{{\mathscr{b}}}\cos {\mathscr{l}}\sin {\mathscr{b}}+{\mu }_{{\mathscr{l}}}\cos {\mathscr{b}}\sin {\mathscr{l}}],\\ {V}_{l} & = & {V}_{\odot ,z}+{V}_{r}\sin {\mathscr{b}}+{r}_{\mathrm{obs}}{\mu }_{{\mathscr{b}}}\cos {\mathscr{b}},\end{eqnarray} \tag{ 3 }$$

where r_obs is the radial distance to the star from the observer, ${\mu }_{{\mathscr{l}}}$ and ${\mu }_{{\mathscr{b}}}$, respectively, are the proper motions in the ${\mathscr{l}}$ and ${\mathscr{b}}$ directions, r_⊙ = 8.3 kpc is the distance to the Galactic center and ${V}_{\odot }=(11.1,232.24,7.25)$ km s⁻¹ 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, σ_x, and velocity, σ_v, along the stream as

$$\begin{eqnarray}{\sigma }_{x} & = & \sqrt{{\sigma }_{{X}_{w}}^{2}+{\sigma }_{{X}_{l}}^{2}},\\ {\sigma }_{v} & = & \sqrt{{\sigma }_{{V}_{w}}^{2}+{\sigma }_{{V}_{l}}^{2}}.\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&{t}_{\mathrm{PM}}=\displaystyle \frac{{\sigma }_{x}}{{\sigma }_{v}}.\end{eqnarray} \tag{ 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 Equation (5) should be correlated with the age of the stream itself.

Zoom In Zoom Out Reset image size

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, r_t, 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

$$\begin{eqnarray}&&{r}_{t}^{3}(z)\approx \displaystyle \frac{{\rm{G}}{M}_{c}(z)}{{V}_{\mathrm{circ}}^{2}(z)}{r}_{\mathrm{gc}}^{2}(z),\end{eqnarray} \tag{ 6 }$$

(c.f., Dehnen et al. 2004) where V_circ(z) is the circular velocity of the host halo (measured at r₂₀₀) at redshift, z, while M_c(z) and r_gc(z), 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, V_circ increases while M_c decreases, resulting in an overall decrease in r_t as a function of time. We consider two possible criteria for a star to be within the tidal influence of the cluster: (i) r < 2r_t and (ii) r < 5r_t, where r 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, t_tidal, as the last epoch when at least 50% of 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 t_PM with t_tidal, we can test how accurately Equation (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 Equations (4) and (5), we estimate the age to be 9.95 Gyr, corresponding to a lookback time of ${t}_{\mathrm{PM}}^{{lb}}\approx 3.01\,\mathrm{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 t = 1.53 Gyr and condition (ii) at t = 1.83 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 ${t}_{\mathrm{tidal}}^{{lb}}$. 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 ${t}_{\mathrm{PM}}^{{lb}}$. 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).

Zoom In Zoom Out Reset image size

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 ${t}_{\mathrm{PM}}^{{lb}}$ and ${t}_{\mathrm{tidal}}^{{lb}}$ 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 (r_h = 3 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 ${t}_{\mathrm{PM}}^{{lb}}$ underpredicts the "true" age of the stream, ${t}_{\mathrm{tidal}}^{{lb}}$. 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 ${t}_{\mathrm{PM}}^{{lb}}$ 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 Equation (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.

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 (Equation (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 10⁵ M_⊙ 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 z = 0 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 data set 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.