STAR STREAMS IN TRIAXIAL ISOCHRONE POTENTIALS WITH SUB-HALOS

R. G. Carlberg

doi:10.1088/0004-637X/808/1/15

1. INTRODUCTION

Tidal star streams pulled out from star clusters have been extensively studied in simulations where the galaxy is modeled as an analytic potential, often taken to be spherically or axially symmetric for simplicity (McGlynn 1990; Johnston et al. 1996; Dehnen et al. 2004; Küpper et al. 2008, 2012; Varghese et al. 2011; Fardal et al. 2014; Ngan & Carlberg 2014; Amorisco 2015; Carlberg 2015; Hozumi & Burkert 2015). These simulations have described how tidal heating ejects stars into a stream. There is a good analytic understanding of how sub-halos interact with thin streams in angular momentum conserving potentials to create velocity perturbations, which open up to gaps in a stream (Yoon et al. 2010; Carlberg 2012; Erkal & Belokurov 2015). Simulations that begin with a cosmological halo show considerably more orbit diffusion than smooth analytic potentials (Bonaca et al. 2014; Ngan et al. 2015; W. Ngan et al. 2015, in preparation; Pearson et al. 2015) which is not yet well understood. However, the increased orbital freedom is certainly related to the rich orbit structure in static triaxial potentials (de Zeeuw 1985; Binney & Tremaine 2008), which is bound to play a role in the aspherical and non-stationary cosmological galactic halos in which tidal streams develop.

In principle, kinematic measurements along a stream can tightly constrain the full phase space parameters of the highly redundant orbits of the stream particles and can serve as a precision probe of a galactic potential (Johnston et al. 1999; Binney 2008; Bovy 2014). The initial velocities of streams of stars that are tidally stripped from star clusters are calculated for a given orbit and potential with collisionless n-body codes. The stars emerge with highly correlated orbital velocity substructure (Carlberg 2015). The stream width scales with the tidal radius and the length of scales with the tidal radius and the age of the stream (Johnston 1998). In a spherical potential, angular momentum and radial action are always conserved, but tidal mass-loss features blur due to orbital overlap down the stream. However, as stars move away from the progenitor on slightly different orbits in which the gravitational field has variations due to aspherical structure in the halo, the streams can disperse relative to a simple smooth potential model (Bonaca et al. 2014; W. Ngan et al. 2015, in preparation). Relative to a spherical potential the streams are almost everywhere increased in width, with some portions greatly blurred. However, a few segments of high or even enhanced density can remain (W. Ngan et al. 2015, in preparation). The view of the stream in the dynamical action variables, which are adiabatically conserved quantities of a potential, can provide insight into what is causing the stream to change its appearance (Ibata et al. 2002).

The purpose of this paper is to quantitatively analyze the dynamical behavior of tidal streams in moderately aspherical potentials that also contain sub-halo populations. The analysis will be done within a simplified potential which is able to reproduce many of the features of a cosmological halo but allows the calculation of action-angle variables at the beginning and end of the simulation. A simplified n-body code that reproduces the tidal stream from a standard n-body code is used to speed up the simulations. The changes in angular momentum and the radial action along a stream as a function of the halo shape and sub-halo content appear as changes to the stream in position and velocity space. These results will be useful to advance the understanding of the mass structure and contents of galactic halos as more and better kinematic data become available and serve as a guide to more realistic simulations.

2. ORBITS IN TRIAXIAL ISOCHRONE POTENTIALS

We want to analyze the orbit changes of a stellar stream in a galactic halo with sub-halos. Action variables, as conserved quantities, provide direct insight into what has changed. Henon's isochrone potential has a simple algebraic expression for the actions (Henon 1959a, 1959b) and provides a reasonable match to a galactic rotation curve over the range of radii where we want to study stream dynamics. One fairly standard way of creating a triaxial potential is to apply a linear scaling to the Cartesian coordinates of a spherical potential to create a triaxial potential. The simplicity of the axis scaling comes at the cost of losing the analytic properties of the orbits in the spherical potential. A number of axisymmetric potentials with good analytic properties have been developed (Eyre & Binney 2011; Binney 2014) but these still do not have the triaxiality, or more generally, the asphericity, of cosmological halos.

With the goal of allowing the rich orbital structure of a triaxial potential, yet retaining some of the analytic properties of the isochrone, we introduce a time dependent potential, which should be a reasonable approximation for cosmological halos. The potential starts as a spherical isochrone, makes a smooth transition to a triaxial form, then remains as a triaxial for a relatively long interval, then reverses to become spherical again. The scale radius and the mass remain constant, the shape being the only thing that changes. We will leave the x axis unaltered, but compress the y–z axes, which alter the density profile, increasing the density in the central region and generally leading to negative densities at a large radius on the most compressed axis. We will also smoothly turn on and off a small fraction of the mass in an identically distributed sub-halo population.

The potential's spherical initial and final form is a device to allow the analytic actions to be calculated and is not intended to be realistic. The action variables of the spherical isochrone will be valid quantities to calculate at the beginning and end of the simulation, although they are not conserved quantities over the simulation. We expect that there will be smooth changes in the action variables associated with the turn-on and turn-off of the triaxiality and local changes as a result of sub-halo encounters, the two being sufficiently different that it might be possible to untangle the effects. Based on these considerations, the triaxial isochrone potential used here is

$\begin{eqnarray}&&{\rm{\Phi }}(x,y,z,t)=\displaystyle \frac{-{GM}(t)}{b+\sqrt{{b}^{2}+{x}^{2}+{y}^{2}{q}_{y}^{-2}(t)+{z}^{2}{q}_{z}^{-2}(t)}}.\end{eqnarray} \tag{ 1 }$

The q_y(t) and q_z(t) are dimensionless functions that begin and end at unity, at which time the potential is spherically symmetric. We largely focus on compressions where the q values are less than one, but cosmological halos also undergo small expansions during mergers with other halos. Our standard compression has q_y = 0.9 and q_z = 0.95, which leads to small negative densities on the y axis beyond y = 6.51, with a negative peak at y = 8.83 with about 10⁻⁴ of the central density. The negative densities occur outside of the orbits that will be studied. The gravitational force is always directed inward.

2.1. Time Variation of Potential Shape

The potential asphericity turn-on takes place over a time t_r, after which it remains on for a time t_f, using the arbritrary function,

$\begin{eqnarray}A(t) & = & 2{(t/{t}_{r})}^{2}\quad 0\leqslant t\leqslant {t}_{r}/2,\\ & = & 1-2{(1-t/{t}_{r})}^{2}\quad {t}_{r}/2\lt t\leqslant {t}_{r},\\ & = & 1\quad {t}_{r}\lt t\lt {t}_{f}+{t}_{r}.\end{eqnarray} \tag{ 2 }$

For $t\gt {t}_{f}+{t}_{r}$ , the rising function is reversed and A(t) returns smoothly back to zero. The rise time, t_r, is chosen to be reasonably slow, typically about three radial orbits, which is neither in the impulsive regime nor truly in the adiabatic regime. The time spent as a constant triaxial potential, t_f, is typically about 15 radial orbits, comparable to the time over which thin streams have been developing in the Milky Way halo.

All calculations here use units with GM = 1 and b = 1. We usually start the calculation with a cluster placed at an x position of four units. We characterize the starting conditions with the angular momentum, L, relative to the angular momentum of a circular orbit in the spherical potential, L_c. Orbits with $L/{L}_{{\rm{c}}}$ of 0.4 and 0.7 have radial orbital periods and are 29.2 and 35.3, respectively, in the spherical potential when started at r = 4. For the time variation function A(t) we select t_r = 80 and t_f = 440 so that the potential becomes spherical again at time 600 and the simulation terminates.

The coordinate distortions are evaluated as ${q}_{y}(t)=1+{a}_{y}A(t)$ and ${q}_{z}(t)=1+{a}_{z}A(t)$ . Positive ${a}_{y},{a}_{z}$ correspond to an expansion of the potential, that is, reduced gravitational attraction at locations off the x axis and negative values correspond to a compression. Both expansion and compression are possible in cosmological situations, though compression likely dominates and is relatively small, since accreted mass spends little time in the central region and largely deposits mass near the current virial radius of the halo. Cosmological halos can also have a rotation of the potential shape, which we will ignore at this time. Louis & Gerhard (1988) studied radial perturbations to the isochrone potential finding that higher order resonances quickly appeared along with complex orbits, leading us to expect behavior of that sort.

2.2. Orbiting Sub-halos

Sub-halos are present in large numbers in the dark matter halos created in cosmological simulations (Klypin et al. 1999; Moore et al. 1999; Diemand et al. 2007; Springel et al. 2008; Stadel et al. 2009). The differential mass distribution is $m\;{dN}/{dm}\propto {m}^{-0.9}$ hence it has a nearly constant mass per decade of mass, ${m}^{2}\;{dN}/{dm}\propto {m}^{0.1}$ . The mass fraction of sub-halos varies steeply with radius relative to the local dark matter density in cosmological halos (Springel et al. 2008). For the radial range around 30 kpc where the few known thin streams have been found, the models predict that about 1% of the mass is in sub-halos, with the sub-halo mass fraction dropping to about 0.1% at a nominal solar radius of 8 kpc.

For the interests of this paper, the practical issue is to have a well-defined sub-halo population over the radial range of the star cluster orbit. Conceptually, our approach is to convert a small fraction, f_sh, of the mass of the fixed isochrone potential into a set of sub-halos with the same collective massive distribution as the isochrone. The sub-halos are free to orbit within the isochrone potential as test particles, which requires that they be randomly drawn from a distribution function, which generates the isochrone potential. Conveniently, analytic forms are available for the isochrone's mass and energy distribution functions (Hénon 1960; Binney 2014). Therefore, on average, there is no change in the summed sub-halo plus diminished external potential and the global and local sub-halo fractions are equal everywhere for the initial spherical potential. Although the density distribution of the isochrone extends to infinity, 59% of the mass is within the r = 4 orbit where we start the star clusters. Hence, little computational effort is wasted on the relatively few sub-clusters at larger radii that do not have orbital pericenters within the range of radii that the star cluster explores. The sub-halos respond to the undiminished triaxial isochrone potential but not each other. The stars of the cluster and the tidal streams feel gravitational forces from all other stars, plus the isochrone reduced by $(1-{f}_{\mathrm{sh}}A(t))$ and the sum of the forces from all the sub-halos, which vary as ${f}_{\mathrm{sh}}A(t)$ , where the A(t) factor, Equation (2), which rises from zero and returns to zero at the end of the simulation.

All of the sub-halos are given identical masses, thereby allowing us to examine mass dependent effects. That is, once a number, N_sh, of sub-halos is chosen, the sub-halos all have masses of ${m}_{\mathrm{sh}}={f}_{\mathrm{sh}}/{N}_{\mathrm{sh}}$ . Each sub-halo is modeled as a Herquist potential, $-{Gm}/(r+a)$ (Hernquist 1990), with a scale radius of ${a}_{\mathrm{sh}}=\sqrt[3]{{m}_{\mathrm{sh}}}$ , as expected for cosmological sub-halos of constant scale density. The dependence of the results on the scale radius of the sub-halos are visible in the simplified analytic analysis (Carlberg 2015; Erkal & Belokurov 2015), where there is little change beyond the scale radius, but the induced velocities increase for stars that have encounters that take them within the scale radius of a sub-halo.

2.3. Test Particle Orbits without Sub-halos

Simple orbit integrations in the model potentials illustrate some of the potential's basic properties. We start a set of 10,000 particles with no self-gravity at x = 4, $y=z=0$ with a uniform distribution in v_y between zero and the circular value at the chosen starting radius. The value of x = 4 is chosen to help avoid the core region where there is little effective tidal force. The particle v_x and v_z velocities are given a small Gaussian random value, 0.002 units, about 0.5% of the circular velocity, to allow the particles to explore mildly off plane orbits. The orbits are integrated with a leapfrog integrator with step sizes of 0.004 units, which is nearly 10⁴ steps per orbit, conserving the actions to better than one part in 10⁶.

Orbits are evolved in four different versions of the triaxial isochrone, with ${a}_{y},{a}_{z}$ values of ±0.05 and ±0.1. The differences of the actions of the particles calculated at the start and end of the orbit integrations are plotted as a function of the initial orbital angular momentum in Figures 1–4. Figure 1 shows changes roughly, as one might expect, as a function of the orbital angular momentum (green points). That is, the range of initial orbital angular momentum means that particles have different orbital phasing relative to the changing axes. Some particles experience a net positive torque and others where the torques lead to a net loss. At high L, hence low eccentricity, the torques largely cancel out, leaving little net change.

**Figure 1.** Changes between the beginning and the end of the simulation in the particle's total angular momentum, L (blue), radial action, J_r (red), and angular momentum about the axis perpendicular to the initial orbital plane, J_z (green), as a function of the initial L in a potential with ${a}_{y}=-0.05$ and ${a}_{z}=-0.10$ . The particles are started at r = 4 with a uniform angular momentum distribution between zero and the circular value. All orbits are started in the plane with an initial vertical and radial velocity dispersion of 0.002 units. A point below the downward sloping line drawn from the origin would have had its orbital direction reversed.
Download figure:
Standard image High-resolution image

**Figure 1.** Changes between the beginning and the end of the simulation in the particle's total angular momentum, L (blue), radial action, J_r (red), and angular momentum about the axis perpendicular to the initial orbital plane, J_z (green), as a function of the initial L in a potential with ${a}_{y}=-0.05$ and ${a}_{z}=-0.10$ . The particles are started at r = 4 with a uniform angular momentum distribution between zero and the circular value. All orbits are started in the plane with an initial vertical and radial velocity dispersion of 0.002 units. A point below the downward sloping line drawn from the origin would have had its orbital direction reversed.
Download figure:
Standard image High-resolution image

**Figure 2.** Same as Figure 1, but for ${a}_{y}=-0.10$ and ${a}_{z}=-0.05$ .
Download figure:
Standard image High-resolution image

**Figure 2.** Same as Figure 1, but for ${a}_{y}=-0.10$ and ${a}_{z}=-0.05$ .
Download figure:
Standard image High-resolution image

In contrast to the simple behavior of Figure 1, Figures 2–4 show nearly chaotic behavior over a range of orbits (Siegal-Gaskins & Valluri 2008). The change in potential axis length remains the same as in Figure 1, but with the y–z axes switched, as in Figure 2, or made positive, the potential axis length was hence allowed to expand, as in Figures 3 and 4. The figures show that nearby orbits have dramatic differences in the changes to their action variables for higher eccentricity orbits. Low eccentricity orbits at high angular momentum are generally much less affected, although the vertical action can still be dramatically affected. This behavior is often seen for orbits that have higher order resonances between orbits and the rate of potential shape change.

**Figure 3.** Same as Figure 1, but for a_y = 0.05 and a_z = 0.10.
Download figure:
Standard image High-resolution image

**Figure 4.** Same as Figure 1, but for a_y = 0.10 and a_z = 0.05.
Download figure:
Standard image High-resolution image

If the potential changes are made very slowly, we expect that the action changes for regular orbits will become adiabatic and become exponentially small (Binney & Tremaine 2008). In Figure 5, the ${\rm{\Delta }}L$ and ${\rm{\Delta }}{J}_{r}$ actions diminish as expected but the changes in the ${\rm{\Delta }}{J}_{z}$ values remain large. Overall, the rich range of orbital behavior of the triaxial isochrone makes it a simplified, but interesting, stand-in for a cosmological galactic halo potential that can capture realistic star stream behavior beyond a spherical halo.

**Figure 5.** Same as Figure 2, with an increased vertical axis, but for t_r = 2400, which is 30 times slower than the standard rate of potential change. The green dots are the changes in J_z.
Download figure:
Standard image High-resolution image

2.4. Test Particle Orbits with Sub-halos

The effects of adding sub-halos to the simple orbit integrations of Figure 2 are shown in Figure 6, where 1% of the halo mass transfers to 1000 identical mass sub-halos. That is, each sub-halo has a mass of 10⁻⁵ of the host halo. The test particles are initiated with no additional random velocities to avoid the confusion of orbit diffusion. The outcome with no sub-halos is shown as the inset in Figure 6.

The sub-halo interactions occur on the scale of the radius of the sub-halos, which is approximately 0.02 units for the case studied. The sub-halos act to create fairly abrupt changes in the action variables, spread out over the orbits in accord with the density profile of the halo. Reducing the mass of each sub-halo by a factor of two, which reduces their scale radius by $\sqrt[3]{2}$ , we find that the sharp changes in J_z for $L\gtrsim 0.8$ are one half the height of the higher mass sub-halo results, indicating that single sub-halo encounters dominate for test particle orbits in this range over the duration of the simulation. For $L\lesssim 0.8$ , the particles encounter sub-halos multiple times leading to much larger J_z changes. The highest eccentricity orbits, those with $L\lesssim 0.3$ are largely unperturbed by sub-halos, likely because of their high velocities through the central regions. Overall, the amplitude of the changes in the actions due to sub-halos are comparable to the changes due to the triaxial potential. We can foresee that disentangling sub-halo effects from potential effects will not be straightforward as might be hoped.

3. SIMULATION OF A SMALL STAR CLUSTER

We create a small, self-gravitating satellite with a mass of $m\ll 1$ ( $m\sim {10}^{-8}$ , where the host halo has M_I = 1) from a King model distribution function (King 1966). The cluster is created with its own dimensionless units of velocity and position, then scaled down in size to fit its zero density radius within the estimated tidal radius. The velocities are scaled to the desired mass and new size. Our standard star cluster model has a central potential, $\psi (0)$ , relative to the central velocity dispersion, σ, of $-\psi (0)/{\sigma }^{2}=4$ . Increasing $-\psi (0)/{\sigma }^{2}$ decreases the core radius relative to the radius where the model density goes to zero. The outer density profile has a slope nearly independent of the core radius, so the mass loss into the tidal stream has a weak dependence on the concentration of the model star cluster. Consequently, the star cluster's core radius only begins to have a significant effect on the tidal stream when mass loss is so large that the core is eroded.

The outer radius of the star cluster model is scaled to be equal to the tidal radius, $r{[m/(3{M}_{I}(r))]}^{1/3}$ , where r is the starting radius of the orbit. The small adjustment to the tidal radius for a distributed mass potential (Binney & Tremaine 2008) is not included, nor is any allowance for a non-circular orbit, which leads to some modest internal changes of the star cluster model in response to the tidal field in its first orbit.

The standard model uses 10⁵ particles in the initial star cluster. The star particles in the model star clusters are subject to gravitational forces from other star particles, the host galaxy gravitational potential, and any orbiting sub-halos that may be present. The rate of mass ejection varies around the orbit and leads to complex, strongly correlated, velocity differences (or action variable offsets) between the stream and the star cluster (Carlberg 2015). The mass loss in the systems of interest is driven by tidal heating and is well-modeled with a collisionless n-body code. As a reference, we use the gravitational and orbital evolution components of the well-known and characterized Gadget2 code (Springel 2005).

3.1. An $n-{Body}$ Shell Code

The precise description of the evolution of a star cluster in the tidal field of a host galaxy remains an important research problem as a generalization of the dynamics of an isolated star cluster. Current work allows for a realistic range of star masses, their mass evolution and their remnants, the presence of binaries and higher multiplicity systems. The evolution of collisional star clusters in tidal fields generally either use some form of gravitational calculation splitting to handle the dynamics on different timescales (Fujii et al. 2007) or interactions via statistical approximations (Brockamp et al. 2014; Sollima & Mastrobuono Battisti 2014). Fukushige & Heggie (2000) demonstrate the complexity of star cluster orbits near the tidal radius that exists for cluster in a circular orbit.

Two well-known thin star streams are Pal 5, which emanates from an exceptionally low concentration cluster (Odenkirchen et al. 2001) and GD-1, which has no known progenitor (Grillmair & Dionatos 2006). Given the complexity and high computational costs of studying a collisional star cluster, and, the observational fact that the importance of star–star encounters in the progenitor clusters is either relatively low, or, unknown, it is currently conventional to study tidal streams as emanating from tidally heated collisionless star clusters (Dehnen et al. 2004). Given that our interest is in the subsequent dynamics of the star streams the assumption is appropriate for the time being, but it will eventually be important to incorporate a more complete dynamical description of the originating star cluster.

The evolution of a very low mass, collisionless, satellite orbiting in an external potential in which dynamical friction is not a factor is well-handled with multipole expansion codes. The relatively low cost and good accuracy of multipole n-body codes has a rich history as the first codes that could handle large numbers of self-gravitating collisionless particles (Aarseth 1967; van Albada & van Gorkom 1977; Fry & Peebles 1980; Villumsen 1982; White 1983; McGlynn 1984). Multipole methods remain popular for moderately distorted spheres with arbitrary density profiles for which there are now efficient parallel processing algorithms (Meiron et al. 2014) and have been generalized to the powerful self-consistent field method (Hernquist & Ostriker 1992). A multipole field provides substantial smoothing of the gravitational field and completely suppresses the complex astrophysics of the two and more body interactions in a dense globular cluster. If desired, star–star collisional physics can be added using a Monte Carlo method (Giersz 2006), though none is included here.

After testing potentials calculated with up to fourth order multipole expansions, we found that a simple shell code, which only calculates the spherical part of the gravitational field, gave results that were essentially indistinguishable, with the shell code being less subject to the problem of defining the center. This is not unexpected in the situation we are modeling since much of the mass of the satellite is in the central regions where the tidal field is weak. On the other hand, the degree of quantitative accuracy needs to be verified. The shell code calculates the radial gravitational acceleration of a particle at radius, r, with a mildly softened gravity, $-{GM}(\lt r)r/{({r}^{2}+{s}^{2})}^{3/2}$ where $M(\lt r)$ is simply the total mass interior to r relative to a central point. The softening s is set at 20% of the core radius of the King model cluster that we wish to simulate and is needed for the case where two particles are near the coordinate center. Making a good choice for the central point of the potential is important, since a poor choice can lead to undesirable oscillations and even an unstable cluster. For the situations of interest here, the simple approach of integrating the orbit of a fictitious particle placed at the position of the coordinate center of the cluster works extremely well and shows no signs of instability.

The calculation sorts particles in order of ascending radii and then calculates the gravitational field at each successive particle, simply as the mass interior to its radius. The sub-halo gravity is calculated with a direct sum. Our shell code has fixed time steps, ${\rm{\Delta }}t=0.004$ , that accurately capture the tidal heating in the star cluster and the interactions with sub-halos down the tidal stream. For many of our simulations, about half of the particles remain in the cluster so that increasing the time step size for particles outside of the cluster would not necessarily give rise to a large speed increase. The code was parallelized with a parallel quicksort and MPI routines to reduce the execution time. Other than the quicksort the algorithm is essentially trivially parallel.

3.2. Comparison of $n-{Body}$ and Shell Code Results

As a test of the accuracy of the shell code, Figure 7 shows the mass as a function of time in Gadget2 and the shell code for a cluster started at r = 2.2 in a spherical isochrone. Note that the difference between the full n-body code and the shell code is almost entirely one of slightly different orbital periods, not of mass lost per orbit. The shell code uses the high accuracy orbit of the center of the star cluster. The angular momentum and the radial action calculated at the end of the simulation when the potential is again spherical are measured relative to the star cluster as a function of their time of ejection from the cluster. The resulting distributions are shown for the n-body simulation in Figure 8 and for the shell code in Figure 9. The two are essentially identical, other than the slight difference in the orbital period, with the shell code having the correct period. We conclude that the shell code gives accurate results for the range of satellites, orbits, and potentials of interest here.

**Figure 7.** Mass in the cluster as a function of time on an $L/{L}_{{\rm{c}}}=0.4$ orbit evolved with an $m={10}^{-7}$ with a full n-body code (solid line) and with a shell code (dashed line).
Download figure:
Standard image High-resolution image

**Figure 8.** Sngular momentum (blue) and radial action (red) relative to the cluster center, ${\rm{\Delta }}L,{\rm{\Delta }}{J}_{r}$ , respectively, of particles as a function of their ejection time for the cluster evolved in the Gadget2 n-body code.
Download figure:
Standard image High-resolution image

**Figure 8.** Sngular momentum (blue) and radial action (red) relative to the cluster center, ${\rm{\Delta }}L,{\rm{\Delta }}{J}_{r}$ , respectively, of particles as a function of their ejection time for the cluster evolved in the Gadget2 n-body code.
Download figure:
Standard image High-resolution image

**Figure 9.** Dame model as in Figure 8, but evolved with the shell code.
Download figure:
Standard image High-resolution image

3.3. Streams in Triaxial Potentials without Sub-halos

Clusters were started at x = 4 with $L/{L}_{{\rm{c}}}=0.7$ and 0.4 in spherical and $[{a}_{y},{a}_{z}]=[-0.1,-0.05]$ triaxial isochrones with no sub-halos. A cluster in a low eccentricity orbit initiated with $L/{L}_{{\rm{c}}}=0.7$ loses 39% of their mass over the interval to time 600 nearly independent of triaxiality. The cluster undertakes approximately 17 radial orbits in the spherical potential. Over the same time, a cluster on a high eccentricity orbit loses 99.7% of its mass in the the triaxial potential and 94% in the spherical potential (21 radial orbits). Although the difference is easily measurable, when considered as a rate, the differences with potential shape are small compared to the dependence on orbital eccentricity. The mass loss does depend on the location of the surface used to measure contained mass. We use 0.03 units, approximately three times the initial tidal radius.

Figures 10 and 11 show the action variables relative to the star cluster of the particles in the stream as functions of their ejection age for two triaxial potentials. Two clear features are that the two tidal tails remain almost completely symmetric, largely a consequence of the small size of the cluster, which has a tidal radius of 0.007 units, about 0.3% of the orbital radius. Consequently, the tidal field is essentially symmetric across the cluster and the two tails are nearly the same, although moving away from the cluster in opposite directions. The second feature is that the radial action and angular momentum of the stars relative to the satellite vary slowly and smoothly over the duration for which particles are in the stream. The change of the angular momentum in the initial orbital plane, ${\rm{\Delta }}{J}_{z}$ , is largely a measure of the subsequent tilt of the orbits. ${\rm{\Delta }}{J}_{z}$ rises exponentially with time, the doubling time being approximately four radial orbits. The exponential growth of the orbit tipping is an important new feature that the triaxial potential brings into play for the tidal stream relative to an orbit started with the same velocities in a spherical or axisymmetric potential. Though the lower panels of Figures 10 and 11 show the same triaxility on different axes relative to the orbital plane they do not show the exponential growth. Hence, orbital structure will depend on the details of the potential and the initial orbits of particles in the potential.

**Figure 10.** Changes of the actions of a stream's particles from a cluster that was started with $L/{L}_{{\rm{c}}}=0.7$ in a smooth triaxial with $[{a}_{y},{a}_{z}]=[-0.1,-0.05]$ (top) and $[-0.05,-0.1]$ (bottom). The actions are measured at the end of the simulation in the spherical potential relative to the cluster center values. The differences are plotted as a function of the time since the particle was ejected from the cluster. The action variables are red for ${\rm{\Delta }}{J}_{r}$ , blue for angular momentum, ${\rm{\Delta }}L$ and green for ${\rm{\Delta }}{J}_{z}$ , whose values are scaled downward by a factor of 30 in this plot.
Download figure:
Standard image High-resolution image

**Figure 10.** Changes of the actions of a stream's particles from a cluster that was started with $L/{L}_{{\rm{c}}}=0.7$ in a smooth triaxial with $[{a}_{y},{a}_{z}]=[-0.1,-0.05]$ (top) and $[-0.05,-0.1]$ (bottom). The actions are measured at the end of the simulation in the spherical potential relative to the cluster center values. The differences are plotted as a function of the time since the particle was ejected from the cluster. The action variables are red for ${\rm{\Delta }}{J}_{r}$ , blue for angular momentum, ${\rm{\Delta }}L$ and green for ${\rm{\Delta }}{J}_{z}$ , whose values are scaled downward by a factor of 30 in this plot.
Download figure:
Standard image High-resolution image

**Figure 11.** Same as Figure 10, but for a cluster that was started with $L/{L}_{{\rm{c}}}=0.4$ . Note that the vertical scale is reduced relative to Figure 10.
Download figure:
Standard image High-resolution image

4. STREAMS IN TRIAXIAL POTENTIALS WITH SUB-HALOS

Two star clusters, started with the same velocites in the ${a}_{y}=-0.1$ , ${a}_{z}=-0.05$ halo of Figures 10 and 11, but now with 1000 sub-halos containing 0.5% of the mass of the halo, are shown in projection (both at two times) onto the orbital plane in Figure 12. The addition of the sub-halos leads to gaps in the stream and helps to push the progenitor and stream onto somewhat different orbits relative to the same potential without sub-halos.

**Figure 12.** Snapshots at times 580 and 600 of a low and high eccentricity stream. In the top panel, the low eccentricity stream is represented by the two blurry structures that come near the center, whereas the high eccentricity stream is represented by the thinner structures that cross them and nearly connect. The contrast scale goes from 0 to 10 particles per 0.01 unit square pixel in the top panel and from 20 to 100 particles per pixel in the bottom.
Download figure:
Standard image High-resolution image

4.1. Stream Visibilities

A star cluster on a high eccentricity orbit is subject to stronger tidal fields at pericenter than a low eccentricity cluster with the same apocenter. The high eccentricity orbit also takes the stream into regions where there are more sub-halos; therefore, we expect high eccentricity streams to have more orbital spreading than low eccentricity streams. To examine this idea, we compare the streams from clusters on $L/{L}_{{\rm{c}}}=0.4$ and 0.7 orbits in Figure 12. A cluster started with $L/{L}_{{\rm{c}}}=0.4$ loses mass at more than twice the rate of a cluster started with $L/{L}_{{\rm{c}}}=0.7$ . Over the somewhat arbitrary simulation duration here the high eccentricity mass cluster is reduced to a few percent of its initial mass and it is only a coincidence of timing that the final moment is caught. Figure 12 shows the projected view of the two streams at two snapshot times. The upper panel shows the location of all particles, showing that the high eccentricity stream has many stars scattered well away from the stream. In currently available observational data, once match-filtered to identify streams, it usually is difficult to identify stream locations where the local surface density is less than about a factor of five in density below the peak values (Carlberg & Grillmair 2013).

Figure 13 shows that the surface density distribution of high and low eccentricity streams has a different character, with low eccentricity streams having relatively more particles at higher surface densities. The distribution is not strongly orbital phase dependent. The peak values are similar, with about 100 particles per pixel (each pixel is 0.01 units square) for both low and high eccentricity streams, though we do note that the high eccentricity stream has somewhat higher peak densities. At a factor of five down from the peak, 20 particles per pixel, about 25% of the low eccentricity stream is visible but only 5% of the high eccentricity stream. The lower panel of Figure 12 shows the regions with surface densities above 20 particles per pixel. More quantitatively, from Figure 13, we see that 95% of the high eccentricity stream is below the 20% of peak surface brightness level, nearly independent of the orbital phase. Similar behavior is seen in the more general simulations of W. Ngan et al. (2015, in preparation), the results here clarifying the combined roles of sub-halos and triaxiality by allowing much of the stream to disperse to low densities.

**Figure 13.** Cumulative distribution, normalized to the total, of surface densities in the low eccentricity stream (green and turquoise lines) and the high eccentricity stream at times 580 and 600 (red and blue lines) in the $[-0.1,-0.05]$ halo shown in Figure 12.
Download figure:
Standard image High-resolution image

The visible segment of the Ophiuchus stream has a luminosity of $1.4\pm 0.6\times {10}^{3}\;{L}_{\odot }$ (Bernard et al. 2014) and is estimated to be on an e = 0.67 orbit (Sesar et al. 2015). If the visible section of Ophiuchus has a similar threshold for visibility as the high eccentricity stream here, then the total luminosity of the entire Ophiuchus stream could be 20 times higher, a total of $3\times {10}^{4}\;{L}_{\odot }$ . These exploratory ideas mainly point to a scenario to understand the very short and low luminosity Ophiuchus stream, which needs to be developed into a model to test its viability.

5. STREAM ACTION VARIABLES

Action variables are constants of the motion in a stationary potential. If the actions have the same distribution of values along the length of a stream, then they can be used along with the stream shape and any available velocities to put strong constraints on the mass and shape of the galactic dark matter potential. The actions of stars ejected into a stream have pulses that vary systematically around an orbit, but these details average together in position space a few orbital cycles further down the stream (Carlberg 2015). As seen in Figures 10 and 11, a time variable triaxial potential induces smooth systematic changes in the potential, which are symmetric relative to the satellite. Although such a smooth variation will complicate modeling, the constraining power of streams remains.

An important issue is how dark matter sub-halos affect stream actions. The effects of a modest sub-halo fraction is shown in Figures 14 and 15, where 0.5% of the halo mass of the same model potentials as in Figure 10 is converted to 1000 equal mass sub-halos. The addition of sub-halos largely leaves the local character of the action distribution in place, a consequence of the probability of a sub-halo hitting the stream just as it emerges being fairly low. However, with increasing distance down the stream, the cumulative action of the sub-halos is to create sharp offsets in the actions that can be up to a factor of two changes relative to the satellite values. However, the absolute size of the changes are typically only 1% of the underlying value, thus, within a variation of about 1%, the action variables along the stream can be considered approximately constant.

**Figure 14.** Actions of stream particles relative to the satellite center for the $L/{L}_{{\rm{c}}}=0.7$ satellite shown in Figure 12, as a function of ejection time. The triaxial potentials $[{a}_{y},{a}_{z}]=[-0.05,-0.1]$ (top) and $[-0.1,-0.05]$ (bottom) contain 1000 sub-halos with 0.5% of the mass. The ${\rm{\Delta }}{J}_{z}$ values are scaled down by a factor of 30 to fit on the plot.
Download figure:
Standard image High-resolution image

**Figure 15.** Same as Figure 14, but for satellites initiated on $L/{L}_{{\rm{c}}}=0.4$ orbits. Note that the vertical scale is a factor of two larger. The cluster essentially dissolves by the end of the simulation so the tidal radius is reduced for more recently lost star particles.
Download figure:
Standard image High-resolution image

6. LOW ECCENTRICITY STREAM KINEMATICS

Low eccentricity streams have less differential shear and less morphological change with orbital phase than high eccentricity streams (Carlberg 2015) making them better suited than high eccentricity streams as indicators of the sub-halo population of the galactic halo. Tidal streams started with $L/{L}_{{\rm{c}}}=0.7$ in an ${a}_{y}=-0.1$ , ${a}_{z}=-0.05$ halo with an increasing mass fraction in a varying number of sub-halos, ranging from 100 to 3000, along with zero for reference, are shown in Figure 16. The plots show the radial position of the particles relative to radius of the orbit of the center of the star cluster at that azimuthal angle. The reference orbit is calculated in the spherical isochrone potential at the end of the simulation. The figures show that the sub-halos are sufficiently numerous and massive to have significant encounters with the streams. At a 1% sub-halo mass fraction, the stream is becoming so badly contorted for any of the sub-halo numbers that 1% becomes an approximate upper limit to the sub-halo fraction that can be tolerated for thin streams, though this is consistent with the predicted local fraction. With the mass entirely in 100 or 300 sub-halos the stream encounters are rarer, but have a larger effect leading to large gaps. As the sub-halos become more numerous and less massive, the stream is able to largely retain its basic shape, with the sub-halos producing gaps in the stream, as is most clearly seen for the 3000 sub-halo case.

**Figure 16.** Stream locations relative to radius of the orbit of the satellite, ${\rm{\Delta }}r$ relative to the angular location of the satellite, ϕ. The particles from satellites on $L/{L}_{{\rm{c}}}=0.7$ orbits with 100, 300, 1000, and 3000 sub-halos. The mass in subhalos in each panel is 0 at the bottom then 0.2%, 0.5%, and 1% (top).
Download figure:
Standard image High-resolution image

The 1000 sub-halo case is repeated for a star cluster of mass $m={10}^{-7}$ in Figure 17 which is 10 times the cluster mass of Figure 16. The stream width and length scale with the tidal radius are thus a factor of 2.15 times larger. Two snapshots are shown to illustrate the variation with the orbital phase. The lower panel of Figure 17 can be compared to the second panel from the bottom of Figure 16, but note that the axes are not the same. Although the sub-halo effects remain constant, the figure also shows that sub-halo-induced gaps remain visible in this slightly more than twice as wide stream from a 10 times more massive progenitor cluster. Figure 17 also shows that the tidal "spurs" on the stream are much more evident, particularly near apocenter, as a consequence of the orbital phase dependence of the mass ejection from the satellite being initiated at pericenter with the particles drifting away at apocenter, as is the case for Pal 5 (Dehnen et al. 2004).

${\rm{\Delta }}r-\phi $ — **Figure 17.** ${\rm{\Delta }}r-\phi$ plot for a satellite that is 10 times as massive as in Figure 16, $m={10}^{-7}$ with no sub-halos at the bottom, then 1000 sub-halos with 0.2%, 0.5%, and 1% (top) of the halo mass, for time 580 (top) and 600 (bottom) to illustrate the dependence on orbital phase with a time near apocenter at the top and near pericenter at the bottom. The bottom panel is comparable to the second from the bottom panel of Figure 16, though the scales are different.
Download figure:
Standard image High-resolution image

The radial and vertical velocities of the streams emanating from the $m={10}^{-8}$ cluster of Figure 16 are shown in Figures 18 and 19. The effects of increasing mass fraction in sub-halos, f_sh, over a range of N_sh are clearly visible in the velocities of 18 and 19. In all cases the amplitude of the perturbations grows significantly with increasing f_sh, which gives a better defined "gap" along the stream. The radial velocities mix together the velocities of the tidal features and have considerable orbital phase dependence, whereas the velocities perpendicular to the orbital plane have little orbital phase dependence. The mean vertical velocity and its dispersion about the mean are shown in Figure 20. Along the stream, the dispersion rises as the orbits spread out in the aspherical potential (Helmi & White 1999; Siegal-Gaskins & Valluri 2008; W. Ngan et al. 2015, in preparation). Further down the stream, the stream has been out of the cluster for a longer time and has been subject to more sub-halo strikes. The mean velocity is knocked away from zero and there are local increases in the velocity dispersion. on a length scale comparable to the scale size of the dark matter halo. The Hernquist radius of the sub-halos is $0.017\sqrt[3]{({f}_{\mathrm{sh}}/0.005)(1000/{N}_{\mathrm{sh}})}$ . For an orbital apocenter of 30 kpc, the sub-halo radius scales to 0.13 kpc. Kinematic detection of sub-halo crossings requires enough velocity measurements per kiloparsec to be able to detect the expected velocity changes. The numbers can be worked out from simulations matched to a stream of interest, but it will be of the order of 100 measurements of stream velocity per kiloparsec to be able to get started. For 0.5% of the halo mass in sub-halos, excursions around the mean velocities reach 0.02 units, which is 5% of the circular velocity of 0.4 units. Otherwise, if scaled to a circular velocity of 220 km s⁻¹ , values of 11 km s⁻¹ are expected in the largest excursions from the mean. For 0.2% of mass in sub-halos, the velocity excursions are proportionally smaller. The velocity dispersion in the thin parts of the streams, which are most readily identified, is about 0.002 units, which scales to about 1km s⁻¹ ; thus, these excursions should be readily detectable given enough velocities along the stream.

**Figure 18.** Same points as Figure 16 except with the radial velocities relative to the orbital radial velocity at that angle.
Download figure:
Standard image High-resolution image

**Figure 19.** Same points as Figure 16 except with the vertical velocities relative to the orbital vertical velocity at that angle.
Download figure:
Standard image High-resolution image

**Figure 20.** Mean (blue) and dispersion (red) in vertical velocities along the stream of Figure 19 for 1000 sub-halos, with the sub-halo mass fraction rising from 0 at the bottom to 0.2%, 0.5%, and 1% at the top.
Download figure:
Standard image High-resolution image

7. DISCUSSION AND CONCLUSIONS

We have developed a framework that is able to reproduce a large range of properties of tidal streams as are seen in galactic halos extracted from cosmological simulations. Basing the potential on the isochrone puts these simulations into a well-studied regime and allows easy dynamical insight into the orbital mechanics at play. A central element of tidal stream analysis is that the distribution of star velocities in the stream can be accurately calculated, here with a spherical shell n-body code. Once free of the satellite potential, stars orbit freely, responding to the overall potential around the orbit and the presence of sub-halos orbiting within the host potential. A question of interest is how much do the orbital actions of stream stars change relative to their initial values as a result of the asphericity of the potential and any sub-halos. If the changes are large, it would diminish the value of streams as probes of the potential shape. We find that for realistic local sub-halo mass fractions of 1% or less that the orbital actions changes are about 1%–2% over approximately 15–20 radial orbits. The stars most distant from the progenitor cluster are the most affected, but the changes are small enough that streams can be approximated, to this precision, as having invariant orbital actions. The outcome is a desirable situation in which sub-halos should have easily measurable local effects in velocity space, but the mean stream motions are not greatly disturbed.

Tidal streams in spherical halo models are much thinner and simpler than those seen in more realistic cosmological halos, whether there are sub-halos or not (W. Ngan et al. 2015, in preparation). The underlying cause is that aspherical potentials have orbits that are not present in a spherical potential, opening up more ways for the stream to spread. Near the progenitor, the stream is sufficiently young that the number of sub-halo encounters is low and the stream dynamics mainly reflect the change in tidal field around the orbit. Further down the stream, where the stream has been orbiting freely longer, the cumulative effect of sub-halos and potential asphericity builds up. In the case of a high ellipticity, $e\gtrsim 0.6$ , star cluster progenitor, the stream orbit spreading can be so large that the majority of the stream would be hard to identify as a high density region on they sky. However, the dense segments of high eccentricity streams certainly are visible, even though they may represent only a few percent of the entire stream, leaving the rest as a ghost. It is therefore possible that a short stream segment like Ophiuchus could be the high density region of a stream with a total luminosity that is a factor of 10 larger.

Velocity measurements are the key to making progress with stream dynamics as indicators of both the shape of the halo and the substructure in the halo. It is now straightforward to reliably predict the velocity distribution of the particles newly injected into the stream, provided that the mass loss is predominantly driven by galactic tides. The presence of sub-halos leads to patches along a stream comparable in size to the characteristic radius of a sub-halo, typically 0.1–0.3 kpc, with an enhanced velocity dispersion and velocity offsets from the mean value for the stream, as shown in Figure 20. It will likely be easier to detect the mean velocity offsets along the stream, rather than local velocity dispersion increases. The large-scale shape of the galactic potential leads to smooth variations along the stream. Although there is already considerable incentive to acquire spectroscopy of stars in stellar streams for abundances and orbital analysis, these simulations show that once large numbers of precision velocities are available for stream stars they will provide definitive tests for the presence of sub-halos.

This paper has concentrated on how sub-halos and triaxial potentials effect action variables, orbits, and velocities and left aside an examination of gaps in streams (Yoon et al. 2010; Carlberg 2012; Carlberg & Grillmair 2013; Erkal & Belokurov 2015) as indicators of the sub-halo population. Gaps will be most useful as indicators in low eccentricity streams, being blurred out fairly quickly in high eccentricity streams (Carlberg 2015). An analysis of gaps and velocity offsets in these streams as sub-halo indicators will be presented elsewhere.

This research was supported by CIFAR and NSERC Canada.

STAR STREAMS IN TRIAXIAL ISOCHRONE POTENTIALS WITH SUB-HALOS

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION