© ESO, 2017

1. Introduction

In the last two decades much progress has been made on the discovery and characterisation of tidal streams around our Milky Way and in other nearby galaxies (see e.g. Koposov et al. 2012; Martin et al. 2014). Tidal streams consist of stars stripped from satellites (dwarf galaxies and globular clusters) that move on nearby almost parallel orbits. As such, they constitute extremely sensitive probes of the mass distribution in the host system (Johnston et al. 1999; Ibata et al. 2001; Johnston & Bullock 2004; Law et al. 2005; Law & Majewski 2010; Koposov et al. 2010; Vera-Ciro & Helmi 2013; Sanders & Binney 2013b; Sanderson et al. 2014). This is one of the drivers of the observational and theoretical studies of streams, as one of the ultimate goals is to establish the mass distribution in the dark halo of galaxies like the Milky Way, which in turn will lead to a better understanding of the nature of dark matter (see e.g. Strigari 2013).

The recently launched Gaia satellite (Perryman et al. 2001) will provide the phase-space coordinates of a vast sample of stars in the Milky Way in the next decade. Together with spectroscopic follow-up surveys of the fainter Gaia stars such as 4MOST (de Jong et al. 2012) and WEAVE (Dalton et al. 2012), these datasets will provide an unprecedented detailed view of our Galaxy. Our understanding of the dynamics of the halo and its streams needs to be sharpened to maximally exploit the wealth of data that will soon become available (see e.g. Johnston et al. 1999; Law & Majewski 2010; Eyre & Binney 2011; Bonaca et al. 2014).

Action-angle coordinates provide an excellent tool to describe the evolution of streams (Tremaine 1999; Helmi & White 1999). For example, the evolution of stars in angle space is linear with time for a static potential, while the actions are adiabatic invariants. The difficulty lies in finding the necessary coordinate transformations to map observables into action-angle space. Only for specific types of potentials, including spherical and those of Staeckel form, can we directly compute the angles and actions because the Hamilton-Jacobi equation is separable (Goldstein 1950; de Zeeuw 1985; Binney & Tremaine 2008). However, in recent years several approximate schemes have been developed to overcome this problem. For example, an appropriate toy potential can be used to compute the true action-angles, based on the work of McGill & Binney (1990) and Kaasalainen & Binney (1994; see e.g. McMillan & Binney 2008; Fox 2012; Sanders & Binney 2014; Bovy 2014). Another method is to approximate the potential locally by a Staeckel potential (axisymmetric or triaxial), a procedure known as the “Staeckel fudge” (Binney 2012; Sanders & Binney 2015).

The availability of full phase-space information for a large number of stars in the Gaia dataset will assist greatly in exploiting the power of action-angle coordinates for streams (McMillan & Binney 2008; Gómez et al. 2010). Actions (and integrals of motion) may be used as well to derive the accretion history of the halo of the Milky Way, because even when a stream is fully phase mixed it is still clumped in this space (Helmi & de Zeeuw 2000; also see Gómez & Helmi 2010). We do require the potential to evolve adiabatically slowly for the actions to remain clumped. This clumpiness in action space can also be employed to determine the gravitational field in which the stars have evolved because the largest degree of clustering occurs when the actions are computed in the true potential (Sanderson et al. 2014; Magorrian 2014; Peñarrubia et al. 2012). The angles and frequencies can also be used to this end, because streams should lie along straight lines that have the same slope in angle and in frequency space for the correct gravitational potential under the condition that it is static (Sanders & Binney 2013a,b). In Buist & Helmi (2015) we argued that an adiabatically growing potential will cause a small difference in these slopes, or an “angle-frequency misalignment”. We also pointed out that in angle and frequency space there are several other indicators that the potential used in the computation may be incorrect. This will allow us to determine the characteristic parameters of the true potential to be separated from its time-evolution.

We continue along that line of study in this Paper, where we explore the behaviour of streams evolved in fully cosmological N-body simulations. In particular, we characterise the streams in the angle and frequency spaces obtained using an approximate spherical potential. We also ran several test-particle simulations to understand and guide our interpretation of the different features we see in the cosmological simulations. Our goal is both to establish which perturbations of the action-angle behaviour of streams exist because of the approximations made on the potential, but also to derive to what degree we can use these coordinates to study cosmologically evolved streams.

The structure of this Paper is as follows. In Sect. 2 we describe the cosmological N-body simulations and the stream catalogue used. In Sect. 3 we discuss the behaviour of these streams in action-angle coordinates, computed using an appropriately chosen spherical potential. To guide our understanding of the behaviour of these streams we present in Sect. 4 a set of test-particle streams evolved in an axisymmetric potential, loosely based on the cosmological simulations. In Sect. 5 we gain further insight into the generic behaviour of streams by determining the impact of computing the action-angles of the test-particle simulations of Sect. 4 in several incorrect potentials, and discuss the effect of self-gravity. We end in Sect. 6 with a discussion and conclusions.

2. Streams in cosmological simulations

2.1. Description of the Aquarius project and its stellar haloes

The Aquarius project (Springel et al. 2008; Navarro et al. 2010) consists of a set of six re-simulations of Milky-Way mass (~ 10¹² M_⊙) dark matter haloes extracted from a larger cosmological parent simulation (Boylan-Kolchin et al. 2009). These haloes were selected to have no close massive neighbours at z = 0, and form late-type galaxies when evolved using semi-analytic galaxy formation models. Haloes Aq-A to Aq-E are believed to be representative of the Milky Way, while Aq-F experienced a recent major merger and hence is less suitable (Wang et al. 2011). We focus our study of streams to two of the haloes, Aq-A and Aq-D, to avoid flooding this article with examples while providing a flavour of the variations found in different systems.

We extracted streams from the accreted component of stellar haloes modelled using the Durham semi-analytic model GALFORM. Cooper et al. (2010) have associated stellar populations with dark matter particles in the simulations via a “tagging” scheme. Lowing et al. (2015) have taken this a step further and generated individual stars from these populations by re-sampling the dark matter particles and using stellar population synthesis modelling. Another difference between these works is that Cooper et al. used the Bower et al. (2006) version of GALFORM, while Lowing et al. used the Font et al. (2011) version which has improved physics on dwarf galaxy scales that makes model satellites more similar to those observed around the Milky Way. Here we use the public catalogue that Lowing et al. offer online¹, which has all stars with magnitude M_g< 7.

The Lowing dataset re-samples the dark-matter-particle positions and velocities in a way that aims to preserve their distribution in phase-space. To observe streams in the halo of the Milky Way, it is common to select the Red Giant Branch (RGB) stars because they are intrinsically bright. We noticed that some of the thinner simulated streams look quite clumpy when only using the RGB stars, most likely because of the re-sampling of the dark matter particles. Since our interest is only in the dynamical features of the streams, we instead decided to use the source (“tagged”) dark matter particles. To this end, we have matched the dark-matter-particle ID’s to those in the outputs of the Aquarius simulations and found the positions and velocities of the source particles from the Lowing dataset. Whether we would have used the RGB stars or the dark matter particles does not matter much for the number statistics as both sets have similar sizes.

We aligned the coordinate system to that of the parent dark halo by using the principal axes determined at 100 kpc from the centre by Vera-Ciro et al. (2011). These authors used the reduced inertia tensor method (Allgood et al. 2006) which closely follows isodensity contours. The z-direction is set along the major axis of these haloes. Further details on the shapes are discussed in Sect. 4.1.

2.2. Mass distribution of the Aquarius haloes

Although intrinsically triaxial, with triaxiality parameters in the range 2/3 to 1 (see Sect. 4.1), the Aquarius haloes can be fitted reasonably well with a spherical mass profile. Here we use Navarro-Frenk-White functional form (Navarro et al. 1996, 1997, hereafter NFW), which provides a relatively good description in the regions where we study streams (r ~ 50–100 kpc, see Springel et al. 2008; Navarro et al. 2010), and we prefer it because of its simplicity and computational efficiency compared to the slightly better fitting Einasto profile (Einasto 1965).

In Fig. 1 we show the results obtained when fitting the spherically-averaged circular-velocity profile with two free parameters, the scale mass M_s and the scale radius r_s, for halo Aq-A and halo Aq-D (see also Navarro et al. 2010). The scale mass is defined as M_s = M(r_s) and is related to the virial mass as (1)where the concentration c ≡ r_vir/r_s. We note that halo Aq-D is much better fitted by an NFW profile than halo Aq-A, which appears to have a bump in its circular-velocity profile. Aq-A is also more triaxial than Aq-D (Vera-Ciro et al. 2011).

2.3. Selection of streams

Not all objects in the Aquarius stellar haloes are apparent as stream-like structures. For example this may be the case if the object is too small or if it has not been significantly disrupted. Figure 2 shows that large streamy structures exist out to large radii (see also Cooper et al. 2010). Each of the two Aquarius stellar haloes have in the order of 100–200 individual progenitors, of which about 20% have produced discernible stream-like features by z = 0, for example, a visible piece of a loop. We selected these by eye in physical space and verified this selection by checking that they also appeared streamy in their approximate angle coordinates, as discussed later. For comparison we also included several massive and seemingly phase-mixed objects that most likely represent debris on quite radial orbits. We also imposed a lower limit of at least 500 “tagged” dark matter particles for each progenitor to exclude the really small streams. This corresponds to a “tagged” dark matter mass of ~ 7 × 10⁶ M_⊙, and yields a total of 35 streams for halo Aq-A and Aq-D together. In the main part of this Paper, we focus on ten streams from each halo to illustrate their behaviour, a selection which includes a few of the highly phase-mixed shell-like structures that appear near the long axis (here the z-axis) and which move on more radial orbits. The remaining 25 streams are shown in the Appendix.

2.4. The morphology of selected streams

Streams consist of groups of stars that have similar orbits, and the relatively small variance in their orbits creates structures whose trajectory follows closely the orbit of the progenitor (Jin & Lynden-Bell 2007; Binney 2008), although not exactly (Choi et al. 2007; Eyre & Binney 2009; Sanders & Binney 2013b). We can therefore analyse the streams in terms of the orbits permitted by the potential.

Our selection of individual streams for halo Aq-A and Aq-D are shown in Fig. 3. To some extent, the consideration of structures in these two different haloes helps us gauge the variety and similarities in the characteristics of their streams. The IDs of the streams given in this figure correspond one-to-one with tree-IDs in the Lowing catalogue. We provide the progenitor’s masses of these streams in Appendix A. The colours in the figure represent the binding energy computed in the best fitting spherical NFW potential (yellow is the most bound, blue is the least bound). Most of the streams depict at least one clear stream-like feature or loop (almost by construction, and for as far as a 2D projection allows us to show this), except for the debris for A98, A104, D56, and D98, which we selected to demonstrate massive and seemingly very phase-mixed objects on radial orbits.

Typically the streams consist of one or more petals, but they are not the clean rosette-like figures seen in the case of spherical potentials (see e.g. Buist & Helmi 2015), because in a triaxial potential many more orbital families exist. In a triaxial potential, orbits are typically box or tube orbits around the major, intermediate, and minor axes. Box orbits get arbitrarily close to the centre of the potential, a property they share with purely radial orbits in a spherical potential, and they do not have a sense of rotation. Tube orbits circulate about one of the axes of the potential and never get to the centre of the potential, which they have in common with loop orbits in a spherical potential. For example, stream D56 and stream A98 are distributed on a structure that appears similar to that defined by a box orbit, while stream A108 resembles a “fish”-orbit (3:2 resonance, see e.g. Miralda-Escude & Schwarzschild 1989; Merritt & Valluri 1999; Binney & Tremaine 2008). We also notice that the streams in halo D also seem to be a bit better defined (i.e. smoother and with longer loops) than those in halo A.

As a first step in our characterisation of streams, and for simplicity, we explore how orbits integrated in the best fitting spherical potential follow the trajectories delineated by the streams. For the relatively massive progenitors, we naturally expect worse agreement, as discussed below. The initial conditions of the orbit were taken from a particle located in the highest density portion of the stream (typically a bound particle in the progenitor). We integrate this particle forward and backwards for 4 Gyr in our best fitting NFW potential to approximately match the streams’ length.

In Fig. 4 we show the resulting orbits on the streams, but in an orientation where the z-axis is in the direction of the mean angular momentum of the stream. An overview of the apocentre and pericentre distances of the spherical orbits is given in Appendix A. For thinner streams, with only one or two wraps, the orbits seem to follow the stream, such as for A158 and A151; but also the much thicker A164 has at least one loop reasonably matched. Many of the heavier streams that have more wraps show a big difference in radial extent when compared to the orbits, with stream A104 giving one of the least satisfactory results. This is not too surprising given that the spherical potential only supports a very limited range of orbit families. Also, a single orbit cannot fit the stream’s radial extent because this depends on the range of energies of the particles, and this is particularly large for some of our objects.

In Fig. 5 we show the streams in the r-v_r plane with the corresponding orbits overlaid. For all objects, individual stream wraps are discernible, even for those objects that are very massive and well phase-mixed. For these it is clear that the radial extent is not well matched by the orbits, as discussed earlier. Part of the complex structure seen in this phase-space projection may also be attributed to the triaxiality of the true potential. Nonetheless, in all cases the spherically integrated orbit seems to match reasonably at least a single wrap. Often the still bound part of the progenitor appears as a vertical diamond-shaped object in the stream. Its location coincides quite well with the location of the highest density of the stream, which was chosen to define the initial conditions for the orbital integration and is indicated by the red cross. Near the progenitors of D82 and D72, we see very long arms that at the ends seem to dissolve in the stream. Most likely these particles are not yet completely unbound from the progenitor (see e.g. Gibbons et al. 2014)

We conclude that our orbits do reproduce some of the wraps of the streams, but typically they lack the complexity that is seen in the Aquarius simulations’ streams. This is partially because a single orbit cannot represent these streams well, but also importantly because the spherical potential does not support the exact same types of orbits.

3. Action-angle behaviour of Aquarius streams in spherical potentials

3.1. General behaviour in a static potential

We now investigate the properties of our streams in action-angle space, where their behaviour is expected to be particularly simple. For an individual particle in a static potential, angles evolve linearly with time as (2)where Ω_i are the orbital frequencies, which only depend on the adiabatically invariant actions and hence are constant. For an ensemble of particles this implies the spread in angles evolves as (3)For streams that have evolved for long enough (), we can ignore the initial angle spread Δθ_i(0), and therefore the appearance of the stream in angle and frequency is rather similar, with the stream in angle space being stretched out in time.

We can understand the shape in frequency space more quantitatively using a Taylor expansion of the orbital frequencies Ω_i with respect to the actions J_j (Helmi & White 1999): (4)with the spreads measured with respect to the centre of mass of the progenitor, and the Hessian of the Hamiltonian H_ij also evaluated at this point. If we diagonalise the Hessian (and assuming the Einstein notation convention), then (5)with V the matrix of eigenvectors e_i of H_ij (as its columns), D_kl the diagonal matrix with the eigenvalues λ_i, and , the action spreads in the eigenspace. In vector notation, this expression takes the simpler form of (6)Generally one of the eigenvalues of the Hessian is much larger than the others² , and this is responsible for the stream being a thin 1D structure in frequency space (and therefore also in angle space), that is elongated mostly in the direction associated to the corresponding eigenvector (Tremaine 1999). We note here that also the relative magnitude of the action spreads is relevant for the direction in which the frequency distribution is elongated.

The action-angle coordinates of course depend on the underlying potential. For example, in a spherical potential there are two independent frequencies Ω_φ and Ω_r (and Ω_φ = Ω_ϑ, apart from a possible sign difference), and therefore for sufficiently long times Δθ_r(t) = ΔΩ_rt and Δθ_φ(t) = ΔΩ_φt. This implies that in this regime in the space of angles, streams follow a straight line with slope S(Δθ) = Δθ_φ/ Δθ_r which is the same slope as in frequency space S(ΔΩ) = ΔΩ_φ/ ΔΩ_r.

3.2. Computing the angles and frequencies

In this section, and as a first approximation, we assume the mass distribution of the Aquarius haloes is represented by our best fitting spherical NFW potential. For this potential we compute the actions, angles, and frequencies (see Goldstein 1950; Binney & Tremaine 2008, for the general procedure). As noted above, in a spherical potential, Ω_φ and Ω_ϑ are equal by definition, however, the corresponding angles θ_ϑ and θ_φ will be different for the particles in our streams. This is because they depend on the current positions and velocities of the particles which are the result of evolution in a different potential than that used to compute the angles³.

In Figs. 6 and 7 we show the distribution of particles in the θ_r−θ_φ and θ_r−θ_ϑ spaces respectively. The most striking feature in these figures is that each of the streams is distributed along more or less straight lines (Sanders & Binney 2013a; Buist & Helmi 2015), even though the host potential is really not spherical. For all streams, the behaviour in θ_r−θ_ϑ space is generally cleaner and the individual streams are seen more clearly than in θ_r−θ_φ space. That such a difference exists is a reflection of the haloes being non-spherical.

We fitted straight lines to the streams, in analogy to what we did in Buist & Helmi (2015), and the results are given also in the figures and in Table 1. Our method to fit lines to the distributions of streams in angle space is discussed in Appendix D. We fitted the distributions in θ_r−θ_φ and θ_r−θ_ϑ separately because as we just saw, the behaviour in these spaces is different. Particles that are still bound to the progenitor are removed from the fitting procedure because they are not expected to follow the general action-angle behaviour of the stream. The result of this removal is clearly visible, for example in streams D72 and A164 where we see a gap in the middle of the distribution in angle space.

The fitting procedure makes evident several noticeable distortions in angle space because the different wraps of a stream are not always on parallel lines, such as for stream A108 and A116. Much more subtle is the deviation seen in the bottom-left of stream A158. These distortions clearly pose a challenge to the determination of the slope. For this reason we have to interpret the slopes determined with our method generally with some care (especially in the case of A116) as the quoted errors do not account for such systematic uncertainties. A quick visual inspection usually helps to evaluate the outcome and the reliability of the fit. The slopes in θ_r−θ_φ space are much more varied than those in θ_r−θ_ϑ space, which, especially for halo Aq-D, all cluster around 0.5. Generally, we also notice that the streams in halo Aq-D seem less distorted than those in halo Aq-A.

In our selection of streams, we also included several objects on very radial orbits and whose debris is distributed in an hour-glass shape in configuration space, such as A98, D56, and D98. In θ_r−θ_φ space, they seem very mixed and do not show very distinct structures, yet their behaviour in θ_r−θ_φ is remarkable, with streams along straight lines being clearly apparent (see e.g. D56 and D98 in Fig. 7). These objects are typically more massive and have deposited debris close to the centre of the halo and therefore are more phase mixed.

3.3. Comparing the slopes in frequency and angle space

In Fig. 8 we show the frequency distributions of the streams. These follow closely a straight line, although sometimes they are quite thick. For A108 and A116 the width in frequency space is not everywhere the same. This is one of the signatures expected when using the wrong potential to compute the frequencies. The determination of the fitted slopes is more robust in frequency space, and these show a considerably smaller range than the slopes in angle space.

The slopes in angle space and frequency space differ, especially when comparing the θ_r−θ_φ space to frequency space. They are expected to be equal in the true (static) potential (Sanders & Binney 2013a,b). In Buist & Helmi (2015), we found that in a time-dependent potential, the slope in angle space is steeper than in frequency space (i.e. S(Δθ) >S(ΔΩ)). Therefore, since the potentials in the Aquarius N-body simulations have grown in time, we might expect the angle space slope to be larger. For some of the streams in our sample this is indeed the case, but there are quite a few streams for which this does not hold, such as D116 and A112, for which the magnitude of the angle-frequency difference is far greater than what can be expected from (adiabatic) evolution of the halo. Also in θ_r−θ_ϑ space, many streams are not in line with our expectations. Overall we find that in θ_r−θ_ϑ, 28 of the 35 streams have a larger slope in frequency than in angle space, that is, they do not follow the expected behaviour, while for θ_r−θ_φ, this is the case for 17 out of the 35 streams.

The simulations of Buist & Helmi (2015) that follow the evolution of streams in a spherical potential, show that differences between the slopes can be obtained when the wrong potential is assumed. For example, S(Δθ)−S(ΔΩ) > 0 when the enclosed mass is too high, and S(Δθ)−S(ΔΩ) < 0 when the enclosed mass is too low. We therefore may attribute the differences in slope to having assumed an incorrect potential (mostly in shape), and not to time-dependence.

3.4. Influence of the potential

We also expect an energy gradient to be present along a stream (Buist & Helmi 2015) and, although this is visible in frequency space, it is less clear in angle space. This is not unexpected as the frequencies depend on energy (via the actions) and more energetic particles move faster, while, on the other hand, the angles depend also on other phase-space coordinates. Even small (of the order of 10%) differences in the characteristic parameters of the potential can lead to the energy gradient being lost in angle space, even when the stream itself has a normal appearance (Buist & Helmi 2015). An exception is A164 in θ_r−θ_φ-space, which does seem to have a continuous energy gradient when following the stream along its various wraps.

To see if the behaviour in angle space and the energy gradient could be improved, we experimented by changing one of the characteristic parameters, namely the enclosed mass for the particular case of A158. However, we were unable to remove the bend seen at the bottom-left in θ_r−θ_φ space, a behaviour that in the spherical case is known to be indicative of wrong values of the characteristic parameters of the potential. It seems clear from the analysis presented in this section that it is the shape of the potential that is wrong rather than the value of the enclosed mass.

We conclude that streams can still look rather regular in angle and frequency space when making an incorrect assumption about the shape of the potential, even if they evolved in a potential that grew via accretion and merging. In later sections of this Paper, we will investigate the conditions for streams to be distributed along straight lines such as those seen in Figs. 6–8.

4. Test-particle simulations of streams in axisymmetric potentials

The Aquarius haloes are not spherically symmetric, and it is important to understand which of the deviations seen in the streams’ behaviour in angle space arise from having assumed a spherical potential. Rather than attempting to study their dynamics in a full triaxial potential, we add only one degree of complexity and now analyse a set of example streams evolved in the axisymmetric Kuzmin-Kutuzov Staeckel potential (Dejonghe & de Zeeuw 1988) with similar density axis-ratios as halo Aq-D of the Aquarius simulations⁴. For the streams we have used as initial orbital conditions those extracted from the present-day positions and velocities of particles in halo Aq-D (as described in Sect. 2.4). An axisymmetric Staeckel potential allows for a relatively straightforward computation of the actions and angles. This implies that we can directly compare their behaviour using true action-angles to those computed assuming a spherical approximation to this potential.

4.1. Potential set-up

To motivate our choice of the characteristic parameters of our axisymmetric potential, we use the Aquarius haloes. These haloes are triaxial, and may be described with two axis ratios: q = c/a (long to short axis ratio) and s = b/a (intermediate to long axis ratio, i.e. a ≥ b ≥ c). A useful quantity is the triaxiality parameter (Franx et al. 1991), (7)which is zero for an oblate halo and equals unity for a prolate halo. Haloes Aq-A and Aq-D have triaxiality parameters between T = 2 / 3 and 1 at a radius of 50 kpc, making them somewhat more prolate, with halo Aq-A being more triaxial than halo Aq-D. Therefore, these haloes are not in the range of being extremely triaxial (T = 1 / 3 to 2 / 3, see also Warren et al. 1992), although at larger radii their triaxiality increases (Vera-Ciro et al. 2011).

The functional form of the Kuzmin-Kutuzov potential is given by Dejonghe & de Zeeuw (1988): (8)where M_K, a_K, and c_K are characteristic mass and scale parameters, respectively. A coordinate transformation to prolate ellipsoidal coordinates (λ, ν, φ) leads to a simple expression, (9)where the relations between λ, ν and R, z can be expressed as (10)We take λ>ν. Depending on the ratio c_K/a_K, this coordinate system can be used to represent a prolate (q_K> 1) or an oblate (q_K< 1) mass distribution (see also the discussion in Dejonghe & de Zeeuw 1988). The potential in Eq. (9) reduces in the spherical limit (a_K = c_K) to Henon’s isochrone potential (Dejonghe & de Zeeuw 1988; Henon 1959). For convenience we link the parameters M_K, a_K, and c_K to the isochrone scale radius r_iso and scale mass M_iso, and introduce the flattening parameter q_K ≡ c_K/a_K. We define the isochrone scale mass as M_iso = M_isochrone(r_iso), such that (11)These relations satisfy a_K + c_K = 2r_iso, that is, the isochrone scale radius is the average of the Kuzmin-Kutuzov axis lengths. For the determination of the numerical values of the parameters, we refer to Appendix C, where we set q_K = 1.87, M_iso = 1.58 × 10¹¹ M_⊙, and r_iso = 22.84 kpc.

4.2. Streams’ set-up

As mentioned earlier, we use a subset of the streams’ positions and velocities extracted from the Aquarius simulations to evolve our test-particle simulations. The stream progenitor is composed of 10 000 particles that follow an isotropic Gaussian distribution in position and velocity, characterised by dispersions σ_pos = 0.3 kpc and σ_vel = 10 kpc/Gyr respectively. This can be roughly translated into a progenitor mass M using that and where we take R = σ_pos, which results in . The orbits are integrated for 10 Gyr in the Kuzmin-Kutuzov potential with the parameters derived from Aquarius halo D.

Figure 9 shows our streams, where the inset labels are the IDs of the corresponding Aq-D streams (e.g. S140 is based on the position and velocity of one particle in D140). We note it is not our goal to reproduce the original streams from the Aquarius haloes in the axisymmetric limit, but to understand the behaviour of streams in the non-spherical regime.

4.3. Action-angles in the true potential

For the Kuzmin-Kutuzov potential, the Hamilton-Jacobi equation separates in prolate spheroidal coordinates (λ, ν, φ). This means that we can compute the true actions (J_λ, J_ν, J_φ), the corresponding frequencies, and the angles in a straightforward manner⁵. At large distances from the centre λ behaves like the radial coordinate r from the spherical case, while ν follows almost the longitudinal angle ϑ. For this reason, we will generally compare Ω_λ with Ω_r, and Ω_ν with Ω_ϑ (and the same for their corresponding angles).

The results for each of the streams in Fig. 9 are shown in Figs. 10 and 11 for both projections of frequency space (top) and angle space (bottom). In angle space, the streams are distributed along straight lines (just as in frequency space), spreading out after 10 Gyr of evolution. The insets in the figures show the slopes obtained from fitting straight lines to the distributions, where the error is estimated from bootstrapping the fits 200 times. These results are also summarised in Table 2.

In the top panels, we also give the angle-frequency misalignment ΔS = S(Δθ)−S(ΔΩ), with the corresponding error. We note that typically | ΔS | ~ 0.01–0.02, which is larger than found by Buist & Helmi (2015) for the streams evolved in a static spherical potential. This is mostly due to the larger progenitor used here which causes the streams to be wider, and this has an impact both on the initial spreads and on the fitting procedure. This is explicitly demonstrated in Fig. 12, where we show a stream on the same orbit as S140 but now for the Carina-like progenitor (σ_pos = 0.1 kpc and σ_vel = 5 kpc/Gyr) used in Buist & Helmi (2015) and evolved for 10 Gyr. We see that the angle-frequency misalignment is now consistent with being zero.

5. Action-angle behaviour of streams in approximate potentials of varying shape

Thus far we have discussed, mostly from a theoretical point of view, what to expect in angle and frequency space in the true potential, and why streams are on straight lines. In this section, we explore what happens when the wrong potential is used to compute the angles and frequencies. To this end we employ the test-particle simulations of Sect. 4 run in the axisymmetric Staeckel potential, but we will assume a spherical potential (Sect. 5.1), a Staeckel potential but with a different flattening parameter (, Sect. 5.2), and a spherical potential with a more dissimilar radial mass distribution (Sect. 5.3). We also explore the effect of self-gravity in Sect. 5.4.

5.1. Spherical approximation

We now compute the angles and frequencies for the test-particle simulations assuming the isochrone potential. This potential constitutes the spherical limit of the Kuzmin-Kutuzov potential, and its characteristic parameters were defined in Sect. 4.1. The resulting angle and frequency distributions are shown in Fig. 13, where we do not include the Ω_r−Ω_ϑ space as it is redundant in the spherical case. Instead, the behaviour in the angles θ_φ and θ_ϑ is relevant since these additionally depend on the physical location of the particles (which encodes information about the potential in which they were evolved as discussed earlier in Sect. 3). The characteristic parameters of the fitted straight lines are also listed in Table 3.

In the top panel of Fig. 13 we show the resulting frequency distributions. We find that quite a few streams have become thicker and distorted in the spherical approximation to the potential in comparison to their behaviour in the true axisymmetric potential. Some of the streams seem quite irregular such as S140, a characteristic that was also seen for some of the Aquarius streams (see e.g. stream A107 in Fig. 8). As a consequence of the distortions, it is generally more difficult to fit straight lines to the frequency distribution, and therefore also the determination of their slopes is less reliable (see e.g. S122 and S138).

To understand this behaviour, it is interesting to explore how orbits map onto frequency space if the frequencies are computed at each time step in the approximate spherical potential. This is shown in Fig. 14 for two different orbits. The trajectories oscillate in frequency, while if the potential had been truly spherical they would collapse into a single point. These oscillations are therefore directly related to how well the isochrone potential represents the mass distribution. We take from this that the thickening of streams is caused by the degree to which the orbits of those particles are described by a spherical potential.

The central panel of Fig. 13 shows the θ_r−θ_φ space. Streams such as S103 and S121 appear quite similar to their counterparts in Fig. 10, plotted in their natural (true) angle space (θ_λ−θ_φ). When the straight line fits look reasonable, the slopes in this spherical angle space and in the corresponding Staeckel angles are similar, differing by ~ 0.05. As for the Aquarius streams, the slopes in θ_r−θ_φ span a range between ~ 0.4–0.7. In all cases, the streams have become thicker. The energy gradient along the streams seems especially discontinuous at some locations, while for some streams, such as S140 and S121, it is somewhat better retained.

The bottom panels of Fig. 13 show the angles θ_r−θ_ϑ. The comparison of θ_λ−θ_ν to θ_r−θ_ϑ is less straightforward because Ω_ϑ and Ω_ν differ by a factor of two in the spherical limit.

We conclude that using a spherical approximation to the potential leads to an increase of the spread in frequency space, by making the streams longer, wider, and/or distorted. However, the energy gradient remains quite intact in this space. The streams in angle space often look much more well-behaved, although they do become thicker and the energy gradient is not preserved. This good behaviour probably reflects that the spherical potential has the right average enclosed mass (as argued by Buist & Helmi 2015). We find that the fitted slopes in angle and frequency space deviate significantly from each other, indicating that the true potential in which the streams have evolved has not been used for the angle and frequency computations. We have explored possible correlations between slopes’ differences and orbital characteristics in our stream sample that could perhaps be used to infer the asphericity of the true potential. However, we found no specific trends in the differences in slopes between the angle and frequency spaces, nor between the angles in the spaces θ_θ−θ_r) and θ_φ−θ_r).

5.2. Axisymmetric Staeckel potential with different flattening

In this section we compute the angles and frequencies for the Kuzmin-Kutuzov potential, but now assuming a flattening , which reverses the axis lengths a_K and c_K compared to the true potential. For the simulations of Sect. 4 we used a prolate potential, and inverting q_K results in an oblate shape.

The resulting distributions in frequency and angle space are shown in Figs. 15 and 16, and are listed in Table 4. They are qualitatively very similar to what we saw in Sect. 5.1 and in previous figures. The frequency distributions are even more broadened than when assuming a spherical potential, and for some of the streams such as S56, also more extended (i.e. beyond the boundaries of the box, which we retained for easier comparison, and which is centred on the progenitor’s centre of mass). This broadening appears to be asymmetric with respect to the progenitor. As in the previous section, the energy gradient in frequency space remains almost intact while in angle space it almost cannot be discerned.

The slopes fitted to the angle and frequency distributions differ considerably, which we for a large part can attribute to the difficulty of the fitting in frequency space. In both angle spaces, we see that some of the streams depict small scale wiggles, such as streams S86, S118, and S140, which are more apparent in θ_λ−θ_ν space. This is probably because the shape of the potential is significantly different from the true form, and from the spherical shape assumed in the previous section.

5.3. Spherical approximation with incorrect radial form

We now focus on the impact of computing the actions and angles in a spherical potential whose radial dependence is quite different from the Kuzmin-Kutuzov potential in the spherical limit. We explore an NFW potential that has the same enclosed mass and mass slope at r_fix = 50 kpc (see also Fig. C.1).

In Fig. 17 and Table 5 we show the resulting distribution of these angles and frequencies for the streams evolved in the axisymmetric Kuzmin-Kutuzov potential of Sect. 4. At first sight, many of the streams look very similar to those in Fig. 13, which corresponds to the computations assuming a spherical isochrone mass distribution.

In frequency space, we see that some of the streams are longer and thinner (see e.g. S140 and S118). The slopes derived from fitting straight lines to these distributions are different from those computed in the isochrone potential, such as for S86 and S103. In angle space, the differences are more subtle, and the typical difference between the fitted slopes of the isochrone and NFW cases is ~ 0.05, which is comparable to the estimated error in the fit, but also the incorrect shape of the potential may contribute here. The energy gradient in angle space for the isochrone and NFW potentials is not always the same, with the isochrone potential showing more fluctuations in the θ_r direction. Overall we find that the frequencies are the most sensitive to the potential, but the differences are small because the slope and the enclosed mass at r_fix = 50 kpc are equal for the NFW and isochrone profiles. The behaviour in frequency and angle space is only slightly worse than for the spherical isochrone potential.

5.4. The effects of self-gravity

In reality, the progenitors of streams will initially be bound by their own self-gravity. Particles that become unbound with time will typically be released at specific points along the orbit (close to the pericentre), rather than continuously as modelled thus far. This results in the leading and trailing arms being offset from each other in configuration space (Johnston 1998) and in energy-angular momentum space (Gibbons et al. 2014). The process of disruption also causes particles to define a bow-tie structure in action (and energy-angular momentum) space. Since one of our goals is to understand the properties of streams in the cosmological Aquarius N-body simulations, we attempt here to establish what the effect of self-gravity is on the distribution of particles in frequency and angle spaces.

To this end, we simulated the evolution of a 3.7 × 10⁸ M_⊙ progenitor for 10 Gyr on an orbit with r_apo = 23.6 kpc and r_peri = 10.2 kpc in a spherical NFW potential with M_s = 1.5 × 10¹¹ M_⊙ and r_s = 12 kpc using an N-body code that uses a quadrupole expansion to model the internal gravitational potential of the system (Helmi & White 2001). Figure 18 shows snapshots for t = 0 and t = 3.5 Gyr. The left panels plot the projection of the stream on the progenitor’s orbital plane, where the red particles correspond to those that are (still) bound. The stream has already spread out significantly by t = 3.5 Gyr because the progenitor is large and the orbit is rather confined to the inner regions of the host, which results in fast evolution. By this time, the progenitor has almost completely dissolved (i.e. only a few particles are marked in red).

In the second panel we show the angle distribution initially and at 3.5 Gyr, with the particles bound to the progenitor marked in red. By 3.5 Gyr many wraps fill up the angle space, to which we fitted parallel straight lines and whose slope is shown in the inset. The frequency distribution at both times is shown in the third column panels, which is overlaid onto the distribution at t = 10 Gyr in grey. We learn from this that at t = 3.5 Gyr, the distribution has reached its final (bow-tie like) shape, and the evolution of the particles is simply governed by θ ~ Ωt by this point in time.

The expected gap near the progenitor in action space forms slowly and is only subtle at the final time (see Gibbons et al. 2014, for a more detailed discussion of the process). Furthermore, we find no significant offset in configuration space between the leading and trailing arms in our simulation, nor any epicyclic oscillations in the stream. This is most likely related to the fast disruption of the progenitor as a consequence of its low density contrast with respect to the host. This is quite different to what is seen for N-body simulations of globular clusters whose disruption process is slower because these are strongly bound gravitational systems (Küpper et al. 2010, 2012).

At 3.5 Gyr, the distribution of angles and frequencies follow each other very closely, as quantified by the slope of the fitted straight lines in the bottom panels, and this is in agreement with the findings of Sanders & Binney (2013b). The difference between test-particle and N-body simulations is especially seen in the bow-tie shape in action space and frequency space, but it is almost absent in angle space. The dynamics of streams is otherwise mostly the same. Furthermore, the slopes of the straight lines along which the stream is distributed in angle and in frequency space are the same, in agreement with the results by Sanders & Binney (2013a). Only for very massive progenitors (~ 10¹⁰ M_⊙) might this picture change because of interactions between the stream and the progenitor (Choi et al. 2009), or because the overall potential of the halo changes while the progenitor system is in the process of being disrupted (Vera-Ciro & Helmi 2013; Gómez et al. 2015).

This analysis leads us to believe that self-gravity has a negligible impact on the results presented in earlier sections of the paper.

6. Discussion and conclusions

We have studied the behaviour of streams in fully cosmological N-body simulations of the formation of stellar haloes from the Aquarius project (Springel et al. 2008). These stellar haloes were produced by tagging dark matter particles in these otherwise dark-matter only simulations according to the GALFORM semi-analytic galaxy formation model (Cooper et al. 2010; Lowing et al. 2015). From these haloes we selected a set of stream-like objects and used the “tagged” particles to build up a catalogue of stellar streams. Our interest was to understand the behaviour of these streams, evolved in a fully cosmological time-dependent framework, particularly when studied in action-angle space.

Since the Aquarius haloes’ potential is not analytic, we explored as a first approximation their behaviour when a spherical NFW potential was assumed. For the best fitting NFW mass distribution, we computed the streams’ angles and frequencies. We found that many of the streams in the Aquarius haloes show several wraps in angle space that appear to be on relatively straight lines, as reported in other works on streams evolved in static or evolving smooth potentials (Sanders & Binney 2013b; Buist & Helmi 2015). However, in many cases these lines are not parallel. We also found patchy features and wiggly behaviour in angle space. In frequency space, often the structures are very broad but relatively linear and depict some amount of irregularity. The width of these streams and features are typically larger than what we have seen before in simple simulations of streams evolved in a smooth spherical potential.

To understand the nature of the various features, we proceeded to explore how various deviations from spherical symmetry could be affecting the behaviour of streams. We have been able to demonstrate that, independently of the form of the host potential, if the angles and frequencies are computed self-consistently then the streams are expected to be along straight lines in frequency and angle space. This is because the Hessian of the Hamiltonian generally has one eigenvalue that is much larger than the others and this dictates to a large extent the direction in which the streams will expand (Tremaine 1999). The exact direction depends also on the action distribution. These results are valid provided the progenitor of the stream is relatively compact in phase-space.

We next focused on why streams evolved in a particular potential but whose angles and frequencies have been computed in a different approximate potential are still on straight lines, as we found for the Aquarius simulations. To this end we ran a set of test-particle simulations in an axisymmetric prolate Kuzmin-Kutuzov potential, which is of Staeckel form. We computed the angles and frequencies in this potential, and found again the characteristic linear appearance of streams in these spaces. Next, we assumed different forms of the potential and computed the angles and frequencies for those cases as well. We found that even if this procedure is not self-consistent, streams are still distributed along relatively straight lines. However, in frequency space, streams became typically thicker and somewhat distorted, and in angle space they depict wiggly behaviour. For example, we found that using a potential with the wrong flattening (spherical or oblate, instead of prolate) has a strong effect on the size of streams in frequency space. On the other hand, differences in the angles and frequencies’ distributions for spherical potentials of different radial form remained subtle provided the enclosed mass was approximately correct within the radial extent probed by the streams’ orbits.

In all cases, the energy gradient along the stream seems almost intact in frequency space (as seen for the Aquarius streams) but clearly distorted or broken in angle space. The straight lines that we fitted to the angle and frequency distributions differ in slope, even when the potential is assumed to be spherical, contrary to expectations. This is the clearest indication that the shape assumed for the angle computation is incorrect and does not correspond to that of the potential in which the streams were evolved. Finally, we also investigated what happens to the actions and angles in a simulation with self-gravity. The largest difference is that during the disruption process of the progenitor, the action distribution of the particles that eventually form the stream is altered, but otherwise the dynamics of streams are the same as in test-particle simulations.

In conclusion, we have been able to reproduce and understand most of the features seen in the approximate angles and frequencies for the Aquarius streams, with the exception of the “noisy” and “patchy” appearance of the streams in angle and configuration space. We believe these can be attributed to interactions of a stream with dark matter substructures, which are known to give rise to disturbed morphologies (Bonaca et al. 2014). Such interactions may also introduce non-adiabatic time-dependent effects on streams that lead to the formation of gaps (Yoon et al. 2011; Carlberg 2013; Ngan & Carlberg 2014).

Finally, since the angle-frequency misalignments found for the Aquarius streams can mostly be attributed to using the wrong potential, this implies that they cannot be used to determine the mass growth history of the Aquarius dark matter haloes as we had proposed in Buist & Helmi (2015). This may be resolved with approximate schemes to compute the actions in a triaxial potential (e.g. Sanders & Binney 2014; Bovy 2014) and by using the distortions of sufficiently thin streams in angle and frequency space to further constrain the present-day potential.

A similar conclusion may be drawn for the determination of the time-evolution of the Milky Way’s gravitational potential, although the challenge in this case is greater. As discussed in Buist & Helmi (2015), the ability to measure time dependence also requires the presence of nearby thin and long streams, as only for such streams will Gaia be able to determine their full phase-space coordinates precisely. Clearly the first step is to have an accurate model for the present-day mass distribution in the Milky Way. Once this has been constructed and we are fortunate enough to be able to exploit the presence of suitable streams, measuring the angle-frequency misalignment to determine the evolution of our Galaxy’s gravitational potential may become feasible.

---

Acknowledgments

We are grateful to Volker Springel, Simon White, and Carlos Frenk for generously allowing us to use the Aquarius simulations. Carlos Vera-Ciro is acknowledged for his help in the analysis of these simulations. H.J.T.B. and A.H. gratefully acknowledge financial support from ERC-Starting Grant GALACTICA-240271 and a Vici grant. We wish to thank Tim de Zeeuw for his helpful suggestions on an earlier draft, and the anonymous referee for their constructive comments on this manuscript.

References

Appendix A: Stream properties

In Table A.1 we list the apocentre and pericentre distances of the orbits integrated in a spherical NFW potential as described in Sect. 2.4. We also provide the total dark matter masses of the progenitors with more than 500 dark matter particles for all the streams considered in that section.

Appendix B: More streams

In Sect. 2 we showed only ten of the streams present in halo Aq-A and Aq-D with more than 500 dark matter particles in their progenitors. In total there are 35 streams that match our selection criteria, and the remaining 15 (5 in Aq-A and 10 in Aq-D) are shown in Figs. B.1 and B.2. As in Figs. 4 and 5, we have overplotted the orbits integrated into the corresponding best fitting NFW spherical potential.

Appendix C: Numerical parameters for Kuzmin-Kutuzov potential

To set the numerical values of the parameters of the axisymmetric Staeckel potential used in the experiments of Sect. 4 we proceed as follows: we first set up the isochrone potential by choosing a radius r_fix at which the enclosed mass and its slope equals that of the best-fitting NFW potential for the Aquarius haloes. This ensures that around r_fix the (spherically averaged) mass distributions are similar. We chose r_fix = 50 kpc because most of the streams we study are located around 50–100 kpc from the halo centre. Having chosen r_fix, we first match the slope of the mass profile and then proceed to set the enclosed mass. The logarithmic slope κ is (C.1)with x = r/r_s and r_s the scale radius of the NFW potential such that M_NFW(r_s) = M_s. The condition of equal slopes at r_fix is (C.2)which we can invert to find r_iso. The scale mass M_iso ensures the enclosed mass is equal at r_fix(C.3)where A_i are the normalised radial mass profiles for the used potentials. In Fig. C.1 we show the result of our fitting procedure for halo Aq-D⁶. The overall mass and velocity profiles differ significantly outside r ~ r_fix.

The next step in setting up the axisymmetric Kuzmin-Kutuzov potential is to use Eqs. (11) to determine its characteristic parameters. This requires some measurement of the shape of the Aquarius haloes, and we use the axis ratios determined using the reduced moment of inertia tensor at r_fix = 50 kpc by Vera-Ciro et al. (2011), and which closely follow the isodensity contours. Since Aq-A and Aq-D haloes are more prolate, we define an axisymmetric equivalent of the axis ratio of the density q_ρ as (C.4)where a, b, and c are the major, intermediate, and minor axis lengths. For halo Aq-D q_ρ(50 kpc) = 1.53 (q_ρ = 1.39 when the axis ratios b/a and c/a are determined at 100 kpc).

The density profile of the Kuzmin-Kutuzov potential is given by (Dejonghe & de Zeeuw 1988): (C.5)We can use this expression to solve numerically for the value of q_K such that the isodensity-contour at R = 50 kpc has the desired axis ratio q_ρ. This results in q_K = 1.87 for halo Aq-D. We note that for this value of q_K, we find q_ρ,K = 1.44 at 100 kpc for the Kuzmin-Kutuzov potential, which is also not too far off the value of q_ρ(100 kpc) = 1.39 measured for Aq-D.

Appendix D: Fitting algorithm

When the individual wraps of a stream are sufficiently distinct in angle space, we can fit straight lines to them. Their slopes can then be compared to those obtained when fitting in frequency space. In Sect. 3, we derived the angle distributions for streams assuming a spherical potential. These distributions are rather noisy and make the process of fitting straight lines non-trivial and complex. Therefore, we proceed to describe now in detail the steps that we take to measure the slopes in angle space.

In Fig. D.1 we show two examples: stream A164 and the test-particle simulations stream S140 as they are being fitted. In the case of the Aquarius haloes, we first use the binding energies of the particles in the stream to find the most bound particle and then centre angle and frequency space on this particle. If there is no bound structure, the centre is put at the location of the highest density in angle space. We generously remove the progenitor and particles in its surroundings by computing the total binding energy with a much higher mass per particle (~ 4). These bound particles are marked in red in the top-left panel of Fig. D.1.

In the next step, we fit a straight line to the remaining particles in each independent projection of frequency space using a simple least-squares method. We show this line for the r−φ projection of angle space in the left panels of Fig. D.1. We then bin the data of this angle-space projection in N bins horizontally and N/S(ΔΩ_φ) vertically, where N = 40 for the test-particle streams and N = 30 for the Aquarius streams, unless there are fewer than 100 particles remaining. We then “clean” this image by emptying the bins that have fewer counts than the median of non-empty pixels. This is done twice to generate enough contrast with respect to the “background” of particles between the streams for the Sculptor progenitors, and only once for the thinner Carina progenitors (e.g. Fig. 12).

We group the particles by connecting pixels that fall within the pattern of Fig. D.2, which is elongated in the same direction as the stream, which is on an angle of 45 degrees if it follows the frequency distribution. Some parallel streams may be “connected” at a single pixel, and the algorithm keeps them separate by not combining groups that individually had more than 25 connected pixels. This procedure results in the groups shown in the third column of Fig. D.1 with different colours. It is clear from the top row that sometimes structures are grouped together that should not be connected, but it is difficult to completely prevent this from happening without removing too many particles from the streams.

We then use a least-squares fit for parallel lines with different offsets. The result is optimised by running the full fitting procedure twice, because we can then use the slope fitted to the angles found in the first iteration when binning angle space⁷. To estimate the errors in the fitting we bootstrap the data 200 times, which was found to be enough to get a reasonable estimate of the errors, although this does not fully reflect the error when the wraps in angle space overlap or if they are not on parallel lines.

All Tables

All Figures