Star Streams and the Assembly History of the Galaxy

The current globular cluster population is the remnant of a larger population that lost its lower-mass members to evaporation driven by internal two-body relaxation and tides (Fall & Rees 1977; Gnedin & Ostriker 1997; Fall & Zhang 2001). However, at least some fraction of those remnants are visible as the thin stellar streams found in the halo (Grillmair & Carlin 2016). Combining globular clusters and their remnant stellar streams into one dynamical picture offers the possibility of providing additional, less "survivor"-biased, insights into the early-epoch globular cluster population and the assembly history of the Milky Way halo.

Two-body relaxation in globular clusters causes stars outside the core to gain energy, which drives a gradual evaporation of the cluster in the tidal field of the galaxy. The evaporating stars form tidal star streams that are useful for measuring the shape of the galactic potential (Binney 2008; Law & Majewski 2010; Eyre & Binney 2011; Sanders & Binney 2013; Bovy et al. 2016) and detecting the many dark subhalos that are expected to be present (Klypin et al. 1999; Moore et al. 1999). The dynamical analysis of stream formation usually starts with star clusters instantaneously inserted into a fixed background potential, although that potential may contain a small fraction of its mass in the form of orbiting subhalos. Some simulations use an approximation that prescribes the release of stars near the tidal surface (Küpper et al. 2012). In n-body simulations of stream formation particle softening makes the clusters collisionless, although tidal heating in a reasonably elliptical orbit is sufficient to cause mass loss. Although the mass-loss rate varies around the orbit, it is periodic so that the the mean mass-loss rate is constant and the dynamical properties of the stars along the stream repeat nearly exactly along the stream (Dehnen et al. 2004; Carlberg 2015a). Much valuable insight results from the analysis of streams and subhalos orbiting in a nonevolving potential (Küpper et al. 2012; Carlberg 2015a; Bovy et al. 2017), but it does overlook the fact that the overall halo potential is assembled over comparable timescales.

Observed stellar streams, or stream segments, are typically about 10 kpc, such as the well-studied Pal 5 (Odenkirchen et al. 2001; Grillmair & Dionatos 2006b) and GD-1 (Grillmair & Dionatos 2006a) streams. There length is somewhat puzzling because globular clusters are 10 Gyr old systems that release stars into the tidal stream at typical velocities of about 5 $\mathrm{km}\,{{\rm{s}}}^{-1}$, which should lead to long streams, ∼50 kpc. Examining how the process of assembly of the globular clusters into our galactic halo affects the properties of the stellar streams and, conversely, how stellar streams might provide insight into the history of the halo and its globular cluster population is the goal of this paper.

The internal dynamical evolution of globular star clusters in a tidal field is fairly well understood (Spitzer 1987; Davies 2013). The mass evolution of globular clusters in the galaxy halo has been discussed at length (Fall & Rees 1977, 1985; Gnedin & Ostriker 1997; Fall & Zhang 2001; Gieles et al. 2011; Lamers et al. 2013). Renaud et al. (2011) have developed an approach to the problem of a star cluster evaporating in a merging system, concentrating on the star cluster evolution. The rate of evaporation depends on the internal two-body relation time and the local tidal field. The tidal field changes from the site of cluster formation to its accretion into the galactic halo. As a cluster evaporates, the rate and velocity at which stars are ejected into the tidal streams change, which in turn affects the density of the tidal stream created.

This paper explores how the observed star streams and remnant halo globular clusters are related to their history of assembly into the galaxy. A toy model of accretion infall of a satellite galaxy containing model globular clusters is developed. The satellite is represented as an evolving analytic potential that falls into the host with dynamical friction. The particles in the model clusters are subjected to a simple Monte Carlo heating procedure that mimics the relaxation processes of globular clusters with the expected mass and size scaling. The evolution of the clusters with time and the resulting properties of the streams on the sky are presented. As an initial practical application, the number of star streams in the recent Pan-STARRS1 volume is compared to the number of globular clusters.

An n-body simulation of the dynamical evolution of a star cluster starting in a satellite galaxy that merges with its host ideally requires a large numbers of particles and a high-precision n-body code that follows stellar two-body interactions but has no two-body interactions with dark matter particles. Here the problem is split into two parts, one being to model the dark matter merger of a satellite and host galaxy and the other to provide internal evolution of the star cluster. The aspect of galaxy merging needed here is a satellite galaxy potential that loses orbital energy and blends into the host galaxy, releasing its globular clusters to orbit in the host galaxy potential. Analytic potentials for the host and satellite, combined with dynamical friction, can meet these requirements. There is no need in this exploratory study to have a full dynamical model of the merger process or detailed cosmological starting conditions.

The the initial host galaxy model is Milky Way 2014 of Bovy (2015). The halo mass component of the model potential absorbs mass from the satellite as it merges, but there is no violent relaxation. A population of dark matter subhalos is included in the host galaxy as fixed mass Hernquist spheres following the prescriptions of Carlberg & Grillmair (2016). The subhalos extend out to 360 kpc. Subhalos below a mass of $0.9\times {10}^{7}\,{M}_{\odot }$ are not included because they make essentially no difference to the streams studied here.

2.1. The Satellite Galaxy Starting Conditions

2.2. The Merger Model

An accurate violent relaxation model is not essential for this investigation, which focuses on the star clusters and their tidal streams. Two analytic time-varying potentials moving together with dynamical friction give an approximate description of the merger. Chandrasekhar's formula (Binney & Tremaine 2008; Boylan-Kolchin et al. 2008) is used to evaluate the dynamical friction, with $\mathrm{ln}\,{\rm{\Lambda }}=10$ adopted. The NFW dark matter halo in the MW2014 potential is used to calculate the local mass density for friction. The NFW is given a small core radius equal to the MW2014 bulge core radius to avoid the singularity when the satellite is dragged to the center of the host galaxy. As the satellite merges into the host, it loses mass and expands to spread over the host halo potential. The key dynamical quantity is ${{\rm{\Omega }}}_{s}$, the local circular frequency of the satellite orbit, which is estimated from the local acceleration in the host galaxy, ${{\rm{\Omega }}}_{s}=\sqrt{{\boldsymbol{a}}\cdot {\boldsymbol{r}}}/r$. The spread in circular frequency of the orbit at radius r across the scale radius of the satellite, a_s, is approximated as $\delta {{\rm{\Omega }}}_{s}\simeq {{\rm{\Omega }}}_{s}{a}_{s}/r$. Therefore, the rate of satellite mass decrease due to orbital divergence is modeled as

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{s}}{{dt}}=-\displaystyle \frac{{M}_{s}{a}_{s}}{{r}^{2}}\sqrt{{\boldsymbol{a}}\cdot {\boldsymbol{r}}},\end{eqnarray} \tag{ 1 }$$

with the lost mass being added to the host halo. The scale radius of the satellite, a_s, expands as

$$\begin{eqnarray}&&\displaystyle \frac{{{da}}_{s}}{{dt}}=\displaystyle \frac{{a}_{s}}{r}\sqrt{{\boldsymbol{a}}\cdot {\boldsymbol{r}}}.\end{eqnarray} \tag{ 2 }$$

These ideas are implemented in the simulation, with the maximum scale radius of the satellite limited to 5 times the scale radius of the NFW of the host galaxy. The mass taken from the satellite is added to the host galaxy halo without changing the scale radius. Under this process the satellite orbits into the center of the host galaxy and expands, and its mass goes to zero. An example merger orbit is shown in Figure 1.

Zoom In Zoom Out Reset image size

2.3. Initial Conditions and Gravity Code

2.4. A Simple Model of Cluster Internal Dynamics

The King model clusters are effectively collisionless over a model Hubble time. Hence, with only weak tidal heating, evaporation of star particles will be very slow, with even the most eccentric orbits losing less than 10% per Hubble time, as shown in Figure 2.

Zoom In Zoom Out Reset image size

The internal dynamics of globular clusters is a rich subject that continues to develop with ever more realistic n-body simulations and sophisticated dynamical analysis. The dynamical evolution of star clusters is reviewed in Spitzer (1987), Binney & Tremaine (2008), and Davies (2013). In brief, two-body relaxation leads to the core shrinking and envelope expansion until heating from central binary- and multiple-star systems intervenes to limit core collapse and sometimes cause core oscillations. The gravitational collisions, largely in the central region, cause stars in the envelope to gradually move to larger orbits until they reach a tidal surface and become unbound from the cluster. Gieles et al. (2011) find that clusters with mass greater than $(8\,\mathrm{kpc}/R)4\times {10}^{4}\,{M}_{\odot }$ are in the expansion-dominated phase of evolution, where R is the galactocentric radius of the cluster. The vast majority of thin Milky Way streams are located outside the solar circle and well beyond. A typical halo cluster at 15 kpc that is more massive than $2\times {10}^{4}\,{M}_{\odot }$ will be in the expansion phase of evolution, which simplifies the modeling requirements.

The controlling timescale for internal evolution of a cluster of half-mass ${M}_{1/2}$ is the relaxation time at the half-mass radius, ${r}_{1/2}$,

$$\begin{eqnarray}&&{t}_{{rh}}\propto {M}_{1/2}^{1/2}{r}_{1/2}^{3/2}\propto \sqrt{\langle {v}^{2}\rangle }{r}_{1/2}^{2}\end{eqnarray} \tag{ 3 }$$

(Spitzer 1969, 1987), where $\sqrt{\langle {v}^{2}\rangle }$ is a representative half-mass velocity dispersion. The relaxation process causes the envelops to expand in a self-similar manner (Hénon 1961; Freitag et al. 2006a, 2006b).

The goal here is to understand how a somewhat more complete cluster and galaxy assembly history affects the stream properties. The heating process is simply modeled with the addition of small velocity kicks to stars in the clusters. The required size of the velocity increments is estimated from the diffusion equation for stellar velocities (Binney & Tremaine 2008) in a cluster with collisions,

$$\begin{eqnarray}&&\displaystyle \frac{d\langle {v}^{2}\rangle }{{dt}}=\zeta \displaystyle \frac{\langle {v}^{2}\rangle }{{t}_{{rh}}({M}_{1/2},{r}_{1/2})},\end{eqnarray} \tag{ 4 }$$

where Gieles et al. (2011) find that $\zeta \simeq 0.1$. In a dynamically correct Monte Carlo simulation an energy increase to outer stars would be extracted from inner stars and binaries (Giersz 1998). This initial investigation does not require an accurate internal model of the cluster, so energy is simply added to the envelope at the predicted rate.

The velocity increments, $\delta v$, will be applied at discrete intervals, $\delta {t}_{r}$. Over some simulation time T, the velocity gains are ${\rm{\Delta }}\langle {v}^{2}\rangle =T/\delta t{(\delta v)}^{2}$. Comparing Equations (3) and (4) gives the required stepwise velocity changes,

$$\begin{eqnarray}&&{(\delta v)}^{2}\propto \zeta \displaystyle \frac{\delta t}{T}\displaystyle \frac{{M}_{1/2}^{1/2}}{{r}_{1/2}^{5/2}}.\end{eqnarray} \tag{ 5 }$$

A practical implementation of Equation (5) is, every $\delta t$ in time, to add random velocities in all three directions drawn from a Gaussian velocity distribution of width,

$$\begin{eqnarray}&&\delta v=C\sigma {\left(\displaystyle \frac{\delta t}{T}\right)}^{1/2}{\left(\displaystyle \frac{{M}_{1/2}}{{10}^{-7}}\right)}^{1/4}{\left(\displaystyle \frac{{r}_{1/2}}{{10}^{-3}}\right)}^{-5/4},\end{eqnarray} \tag{ 6 }$$

in the dimensionless units of the calculation, where C is a constant around unity that absorbs ζ and σ is a characteristic velocity dispersion in the cluster. The $\delta v$ value is used to generate random velocity offsets in the xyz directions from a Gaussian distribution. In this paper $\delta v$ is fixed for any model cluster; however, it could be normalized to the evolving rms velocity dispersion of the cluster.

A little experimentation is required to fine-tune the parameters controlling the heating process. An outer radius for the evaluation of the rms velocity and the half-mass radii of the clusters is selected to r_m = 0.02 units (160 pc), which is somewhat larger than the initial tidal radius of the clusters. No heating is done for particles beyond r_m or inside the radius containing 1/3 of the mass inside r_m, the latter to roughly replicate the results of detailed collisional n-body models. The heating calculation is done every $\delta t$ of 100 steps ($0.71\times {10}^{7}$ yr). Equation (6) is evaluated with $C\sigma \sqrt{\delta t/T}$ set equal to $3\times {10}^{-4}$ for T = 400, the duration of the simulation. The simulations have 200,000 time steps, so this random walk process happens 2000 times.

The simulations have 24 star clusters arranged in two rings of 12 clusters in the satellite galaxy. Simulations with all clusters having an initial mass of $0.9\times {10}^{4}\,{M}_{\odot }$ and $0.9\times {10}^{5}\,{M}_{\odot }$ are done. The simulations reported use a satellite galaxy with M_s = 2, which can be considered a 22% mass satellite, relative to the MW2014 model to its nominal virial radius. The simulations run for 400 time units, or 14.2 Gyr, where one time unit is 0.0355 Gyr. All simulations use the same heating parameters and merger model. Simulations have been run for many other satellite masses, satellite starting conditions, subhalo contents of the host galaxy, and the cluster internal relaxation parameters; however, the set presented provides a good guide to understanding the general outcomes.

Figure 3 shows the fractional cluster mass as a function of time for initial cluster masses of $\sim {10}^{4}\,{M}_{\odot }$ (top panel) and $\sim {10}^{5}\,{M}_{\odot }$ (bottom panel). As expected, the $\sim {10}^{4}\,{M}_{\odot }$ clusters largely, but not completely, evaporate in a Hubble time, whereas the $\sim {10}^{5}\,{M}_{\odot }$ clusters, on average, lose less than half their mass.

Zoom In Zoom Out Reset image size

The mass-loss rates are shown in Figure 4, with the clusters separated into the inner ring (blue) and the outer ring (red). The peaks in the mass-loss rates with time occur at orbital apocenter, where stars that were unbound earlier move away from the cluster. The inner clusters initially have higher mass-loss rates, most clearly shown in the top panel of Figure 4, which is a consequence of the higher tidal fields nearer to the center of the satellite galaxy. The outer clusters tend to be launched into more eccentric orbits, where they experience stronger tidal fields at orbital pericenter and have higher mass-loss rates. Two of the 24 clusters in the bottom panel have sufficiently elliptical orbits that tidal fields help drive complete evaporation.

Zoom In Zoom Out Reset image size

Although the higher-mass clusters have a lower fractional rate of evaporation, their absolute mass-loss rate is slightly higher than that of the lower-mass clusters, as visible in Figure 3, with the bottom panel scale increased by a factor of 10 to allow for the mass difference. At around a time of 200 time units (7 Gyr) the mass-loss rate is around $8\times {10}^{3}\,{M}_{\odot }\,{\mathrm{Gyr}}^{-1}$ in the top panel and $1.2\times {10}^{4}\,{M}_{\odot }\,{\mathrm{Gyr}}^{-1}$ in the bottom panel, about 1.4 times higher. Detailed calculations and the resulting scaling relations (Baumgardt 2001; Baumgardt & Makino 2003) expect an approximately ${M}^{0.25}$ dependence of the mass-loss rate, which would give an increase of 1.8.

The $\sim {10}^{5}\,{M}_{\odot }$ clusters here are comparable to the direct n-body simulations for a similar-mass cluster on a similar range of orbital eccentricities (Zonoozi et al. 2017). All the clusters have orbital apocenters near 250 kpc. The two clusters that dissolve have pericenters of 4 and 6 kpc, whereas the two that drop to about 1/3 and 1/2 of their initial particle count (the masses drop a little more due stellar) have pericenters of 8 and 10 kpc, respectively.

The characteristic dynamical quantities of a single cluster and its stream are illustrated in Figure 5 as a function of time for individual clusters at masses $\sim {10}^{4}$ and $\sim {10}^{5}\,{M}_{\odot }$, in the top and bottom panels, respectively. The two clusters are on essentially the same orbit, with an eccentricity of 0.86. The half-mass radius of the clusters, around 20 pc, is well inside the tidal radius of the clusters, which varies around the orbit and is typically about 100 pc. The particle density distributions have no crisp cutoff. The half-mass radius expands modestly, approximately 50%, over the course of the cluster evolution as long as the cluster retains significant mass. The rms velocity dispersion of the cluster stars is largely dependent on the cluster mass, dropping to zero if the cluster evaporates completely. The spray of (green) points in Figure 5 shows the current angular distance from the progenitor cluster (or its remnant center) of stream particles as a function of the time at which the star particle went through a fixed radius near the tidal surface. The decline of velocity dispersion as the cluster evaporates will be directly reflected in the properties along the stream. The lower-mass cluster loses relatively more stars when the cluster was still in the satellite galaxy. These stars spread to relatively large angular distances away from the progenitor cluster. The particles that evaporate late do so at such low velocity dispersions that they remain relatively close to the cluster compared to the higher drift velocities of particles released earlier. The more massive clusters lose relatively few stars while orbiting in the satellite galaxy, after which the cluster maintains a fairly uniform low rate of mass loss, ejecting stars at velocities that only allow particles to reach angular separations of 20° or so from the cluster. Stream properties are cluster mass dependent.

Zoom In Zoom Out Reset image size

The mean particle density summed over all 24 streams from a set of equal-mass clusters is plotted as a function of azimuthal angle with the orbit projected onto the $x-y$ plane in Figure 6. The average largely removes the orbital phase variations of individual streams. The particle density in the lower-mass streams is higher because the clusters essentially completely evaporate. The cluster models all have 50,000 star particles, so a conversion to stellar density would lead to a relative increase of a factor of 10 for the higher initial mass cluster streams of the bottom panel. That is, the higher-mass clusters here will, on average, be easier to find on the sky, for the same orbital parameters. Yet higher-mass clusters evaporate so slowly that their streams are very low density. An implication is that star clusters with initial masses around $\sim {10}^{5}\,{M}_{\odot }$ will produce the highest sky density star streams.

Zoom In Zoom Out Reset image size

The density along individual streams from the inner dozen clusters at $\sim {10}^{4}$ and $\sim {10}^{5}\,{M}_{\odot }$ clusters done in a simulation with subhalos present is shown in Figure 7. There is a substantial variation in stream density profiles depending on the details of the orbits and the subhalo interactions. The considerable variation of the stream density profiles from one orbit to the next is clearly visible, much of it due to orbital phase where the stream near apocenter increases in linear density on the sky and declines near pericenter. Subhalo-induced gaps are visible in most of the stream densities. The orbits of both clusters and subhalos are the same for the two sets of cluster masses, and the gaps are nearly identical, as are the subhalo interactions, so the differences between the two mass sets can largely be attributed to cluster mass dependence. The most prominent gaps range in size from a degree up to about 10°, with the most visible gaps being around 5° in length.

Zoom In Zoom Out Reset image size

The $x-y$ projections of streams from $\sim {10}^{4}\,{M}_{\odot }$ and $\sim {10}^{5}\,{M}_{\odot }$ clusters are shown in Figure 8. To a good approximation the cluster orbits are independent of the cluster masses, so for every cluster in the top panel there is a matching one at 10 times the initial mass in the bottom panel. All of these clusters are on orbits with apocenters near 240 kpc, but with a wide range of pericenters. The lower-mass clusters on the top make narrower streams, and the sky density is less concentrated to the location of the progenitor cluster, as also shown in Figure 6. The satellite galaxy model potential can be adjusted to produce smaller apocenters.

Zoom In Zoom Out Reset image size

The simulations here lead to the conclusion that streams from star clusters at masses around ${10}^{5}\,{M}_{\odot }$ should be the most visible on the sky, when converted to number of stars per unit angular length of the stream. The density and width vary along the stream, with the ends being relatively low linear density and relatively wide and hence reduced surface brightness relative to near the progenitor cluster.

Streams are very low sky density features that in photometric data alone become detectable through the use of color–magnitude diagram filtering (Rockosi et al. 2002), which optimally weights some desired common age, metallicity, and distance population of stars in preference to field stars. Even with such procedures currently available, data generally lead to streams with a signal-to-noise ratio per degree length of stream of a few. If the stream is at a sufficiently large distance that a large number of stars near the turnoff from the main sequence are not within the sample, about 15–20 kpc for current images, then the signal-to-noise ratio drops further. The outcome is that finding star streams is difficult at best, and it is essentially impossible with current data to detect lower-density regions of streams.

In Figure 9 the sky density is shown for all 24 streams in the bottom panel of Figure 8. In the top panel of Figure 9 any pixel on the sky with a star is black, irrespective of the numbers. A pixel is 0.002 units in the Aitoff–Hammer coordinate, which has a maximum of $2\sqrt{2}$. Scaling to 360° indicates that the pixels are approximately 013 on a side. To more accurately represent sky density, in the middle panel the pixel density is de-weighted, with the inverse distance squared for stars beyond 16 kpc and the resulting pixel density plotted with a gray scale from 0 to 4 units. The bottom panel shows the expected sky densities for the $\sim {10}^{4}\,{M}_{\odot }$ streams on a gray scale of 0–40 to mimic the conversion to the brightness of stars.

Zoom In Zoom Out Reset image size

The weighted map of the sky density of streams is a much more realistic view of what should be expected in current sky maps (Ngan et al. 2015). The streams are thinner and shorter than they would be if all stars were detected, which will eventually be possible. However, the streams do not disappear entirely. That is, even though the length and width of streams are likely being underestimated, the counts of streams are likely to be accurate for streams within double the color–magnitude main-sequence turnoff detection distance. The details of the stream selection function will be studied more quantitatively in a future paper.

For every completely dissolved globular cluster, there will be a thin star stream with its stars distributed close to the orbit of the cluster. The streams also provide an indicator of when their progenitor cluster fell in through their length. A cluster that recently merged into the galaxy will have a short tail, and old additions should have very long tails. Creating a sample of both streams and globular clusters with well-understood sampling is not yet possible, but an interesting start can be made.

The PS1-3π survey provides nearly uniform sky coverage over a large sky area north of decl. −30. The Bernard et al. (2016) stream search is defined within a maximum distance of 35 kpc, which here is taken as a galactocentric radius. The minimum galactocentric radius is 10 kpc, inside of which streams are not found, possibly because they encounter the stellar disk and its massive gas clouds and are quickly eroded (Amorisco et al. 2016). Within the PS1-3π sky Bernard et al. (2016) recover nine previously known thin streams (wider, dwarf galaxy streams are excluded here) and find five new ones, for a total of 14. A few more streams are known in this sky area from other work (compiled in Grillmair & Carlin 2016), but PS1-3π provides the most uniform and large area of sky available at the moment. A few of the streams have mass or luminosity estimates, such as GD-1, estimated to have a mass of $2\times {10}^{4}\,{M}_{\odot }$ (Koposov et al. 2010). There is no estimate of how much mass associated with any stream may be missed at low surface brightness or dust obscuration. The simulations presented here indicate that the streams are likely in the mass range around ${10}^{5}\,{M}_{\odot }$. Streams above this mass limit should be readily visible.

The Harris (1996) catalog is used to construct a sample of globular clusters in the same sky volume. Remnant halo clusters are those visible within the same sky volume as the star streams. To a first approximation, the obscuration of the galactic disk should be about the same for streams and the halo population of globular clusters. An upper limit on metallicity of [Fe/H] of −1 is imposed to focus on halo clusters, although this only eliminates three very low luminosity clusters out of 27 in the PS1-3π volume. The selected clusters have a mean galactocentric distance of 19 kpc. Masses are assigned to the globular clusters using a uniform M/L of 2 solar units and plotted in Figure 10. The pincushion symbol shows the 14 PS1 streams, assigned a mass of ${10}^{5}\,{M}_{\odot }$, although that is likely an upper limit, as indicated with the arrows.

Zoom In Zoom Out Reset image size

The mass distribution of current-epoch dense star clusters is a steep power law, approximately ${M}^{-2}$ (Portegies Zwart et al. 2010). Since the high-mass clusters in the halo are expected to have largely retained their initial mass, the ${M}^{-2}$ can be fit to the higher-mass clusters to predict the numbers at lower masses. In Figure 10 the dashed line is $7{(M/{M}_{* })}^{-1}\exp (-M/{M}_{* }){\rm{\Delta }}M/M$, with ${M}_{* }={10}^{6}\,{M}_{\odot }$ and ${\rm{\Delta }}{\mathrm{log}}_{10}\,M=0.2$, which approximately fits the mass distribution of the half dozen clusters above $3\times {10}^{5}\,{M}_{\odot }$. The dotted histogram is the cumulative difference between the clusters observed to be present and their assumed progenitor numbers, that is, the number of "lost" clusters. The difference between the 28.6 clusters expected from the dense star cluster progenitor distribution and the 14 clusters more massive than ${10}^{5}\,{M}_{\odot }$ is 14.6. Given the uncertainty on the stream masses, this primarily means that the progenitor cluster distribution could not have been steeper than ${M}^{-2}$. It also suggests that essentially all of the massive streams in the PS1 volume have been found, but there could be many more lower-mass streams.

The processes of galaxy assembly and star cluster evaporation need to be included to provide an accurate dynamical history of the development of stellar streams. This paper develops a toy model with a simple merger model and a cluster heating function that approximates the internal dynamics of star clusters. Clusters expand and cool as they lose mass, causing the stream linear density to vary with distance from the progenitor clusters. The density decline of streams in an evolving potential contrasts with the effectively constant linear density of streams produced in static background potentials (Küpper et al. 2012; Bovy 2014; Carlberg 2015a), setting aside orbital phase variation. The portions of the tidal streams created prior to merging of the satellite galaxy into the host are often very widely dispersed.

The heated model clusters have mass loss in rough agreement with detailed cluster dynamical studies. The clusters having a mass around ${10}^{5}\,{M}_{\odot }$ lose about one-third of their mass in a Hubble time, and ${10}^{4}\,{M}_{\odot }$ clusters largely evaporate in about 2/3 of a Hubble time. The linear density of the resulting stellar streams declines by about a factor of 10 between 10° and 50° from the progenitor cluster, although there is considerable case-to-case variation. When projected onto the sky, the decline of the mean density away from the progenitor helps explain why the readily visible parts of most streams are relatively short, even though the simulated streams often wrap a good fraction of their orbit around the galaxy. The observational signature of this explanation will be quite different from short stream segments, due to halo triaxiality (Pearson et al. 2015), chaotic orbits (Price-Whelan et al. 2016), or preferential escape of low-mass stars (Balbinot & Gieles 2017).

Dark matter subhalos always create gaps along the streams, but because gaps are largely changes parallel to the orbital motion (Yoon et al. 2011), the streams remain robust structures (Helmi & White 1999) and not readily blurred out. The result of these models is an improved conceptual understanding of the development of streams and their sky density within an accreting halo, but almost every aspect of the simulation can be improved in more detailed future investigations.

The globular clusters around ${10}^{5}\,{M}_{\odot }$ that dissolve should all produce streams sufficiently bright to be detectable with optimized color–magnitude filtering methods if they are in an accessible part of the sky. Therefore, when combined with the remnant globular clusters in the same volume, streams can provide an estimate of the total cluster population present at early epochs. The current data do not have well-determined stream masses, and the stream selection function is not well understood. However, an interesting star can be made with the PS1-3π survey, which finds 14 streams (Bernard et al. 2016) in its volume. In the same sky volume there are 14 globular clusters to ${10}^{5}\,{M}_{\odot }$. Making the assumption that all the streams are all at or above ${10}^{5}\,{M}_{\odot }$ in total mass then leads to the conclusion that there were 28 clusters in the progenitor cluster population to ${10}^{5}\,{M}_{\odot }$, which would be in accord with the expectations of an ${M}^{-2}$ initial cluster mass distribution. The minimum mass of the current stream sample is almost certainly somewhat lower than ${10}^{5}\,{M}_{\odot }$, which then makes ${M}^{-2}$ the upper limit to the shape of the progenitor mass distribution. A possibly interesting implication is that if binary black holes (BBHs) at least partially form in clusters with low-metallicity masses around ${10}^{5}\,{M}_{\odot }$ (Chatterjee et al. 2017) within the first few gigayears prior to assembly into the Milky Way, then a prediction of BBH numbers based on the current globular cluster population is an underestimate, by a factor of two or less, based on the stream counts.

This research was supported by CIFAR and NSERC Canada. An anonymous referee provided constructive criticism and pointed to an error in the mass scaling in an earlier version.