HALO ORBITS IN COSMOLOGICAL DISK GALAXIES: TRACERS OF FORMATION HISTORY

The MUGS sample of galaxies were produced using an N-body+SPH simulation of the formation of galaxies in a fully cosmological context. Dark matter halos of mass ∼10¹² M_☉ at z = 0 were identified in a dark-matter-only (N-body) ΛCDM (WMAP3 parameters; Spergel et al. 2007) simulation within a box of comoving length 50/h Mpc, without regard to halo spin or mass accretion history. These halos were then re-simulated using GASOLINE (Wadsley et al. 2004) at high resolution using the volume renormalization technique to allow for high resolution in the region of interest while simultaneously including the full tidal field. The simulations include metal cooling, heating from UV background radiation, star formation, stellar energetic and metal feedback, and metal diffusion. Star formation and stellar feedback are included via recipes based on the "blastwave model" (Stinson et al. 2006, 2010). In the high resolution region of these simulations, a dark matter particle has a mass of 1.1 × 10⁶ M_☉, a gas particle has an initial mass of 2.2 × 10⁵ M_☉, and star particles form with masses of <5.5 × 10⁴ M_☉. All particles have a gravitational softening length of = 312.5 pc. Full details of the MUGS simulations and properties of the sample galaxies can be found in Stinson et al. (2010). The disk galaxy analyzed in this paper (g15784) is the largest of the disk galaxies from the MUGS sample.

The galaxy g15784 contains a bulge which is fitted by a de Vaucouleurs R^1/4 profile with an effective radius of r_e = 0.67 kpc, and a disk with an exponential scale length of 3.38 kpc. The maximum halo circular velocity is 360 km s⁻¹ and its virial radius (r_vir) is 290.8 kpc. We use the definition of r_vir (Bryan & Norman 1998) which corresponds approximately to the radius within which the average density of the halo is ∼100 times the critical density.⁷ This galaxy experienced its last major merger at z = 2.5 about 10 Gyr prior to the time at which we examine the simulation.

Cosmological simulations predict that the dark matter halos of all galaxies contain hundreds to thousands of dark subhalos (Moore et al. 1998), although only about 10% of the total mass of the halo is in such subhalos. Subhalos in g15784 were identified by the AMIGA halo finder AHF⁸ (Gill et al. 2004; Knollmann & Knebe 2009). To determine the frozen potential in which to integrate orbits of stars and dark matter particles in g15784, we included mass within four virial radii to avoid an unphysical discontinuous density change beyond the virial radius. This region included several small satellite galaxies and numerous dark subhalos and a small mass contribution at the edge of the simulation volume from a neighboring 10¹² M_☉ galaxy (Nickerson et al. 2011). Both the central stellar spheroid and luminous subhalos are unrealistically massive due to overcooling. Both dark and luminous satellites can perturb orbits of halo particles. Subhalos found by AHF were removed; however, a few subhalos of mass ∼10⁷ M_☉ remained in the inner regions of the galaxy, where automated halo finders have the most difficulty distinguishing substructure (Knebe et al. 2011), making it difficult to fully assess the effects of the subhalos on the orbits of stars and dark matter particles. We tested the characteristics of orbits integrated both with and without the gravitational effect of subhalos included in the frozen potential and found a higher fraction of chaotic orbits when subhalos were included. We also carried out controlled simulations where subhalos were allowed to evolve (Debattista et al. 2008—hereafter D08). We find that both frozen and evolving subhalos increase the fraction of chaotic orbits. However, the population that is most strongly affected is the resonant orbit family. These orbits which are "thin" (i.e., they occupy only a small volume of physical space), were shown to be important in halos where disks grow adiabatically (V12). When subhalos are "frozen" in place, a resonant orbit may repeatedly encounter a subhalo that lies close to that volume, causing it to become chaotic. In this paper, we only consider orbits integrated in the frozen potential from which subhalos more massive than 10⁷ M_☉ were removed. This does not reduce the potential for scattering from smaller subhalos or irregularities such as tidal streams however, and we will see that this scattering probably contributes to the large fraction of chaotic orbits.

We compare the orbits of particles from the cosmological simulation to orbits drawn from two controlled simulations of disk galaxies grown adiabatically inside dark matter halos (from V10 and V12).

In the controlled simulations, the disks form quiescently, and hence alter the halos adiabatically. They also do not have subhalos and other irregularities arising from hierarchical structure formation. The first controlled simulation SA1 (see D08, V10, and V12 for details) is a triaxial dark matter halo in which a disk of particles representing a baryonic mass fraction of 2.5% was grown (starting from nearly zero initial mass) adiabatically and linearly on a timescale t_g = 5 Gyr with the disk's symmetry axis parallel to the short-axis of the halo. The triaxial N-body halo used in simulation SA1 was produced by multi-stage merger of spherical NFW halos. In the second controlled simulation, SBgs, a baryonic stellar disk forms self-consistently from 2.8 × 10¹¹ M_☉ of initially hot gas (constituting 10% of the total mass) in a prolate NFW halo (V12, Debattista et al. (2013)). In SBgs, the hot halo gas initially has the same spatial distribution as the dark matter and slowly cools to form a disk with symmetry axis along the short axis of the halo. The halo and gas particles are given an initial specific angular momentum j, determined by an overall cosmological spin parameter λ = (j/G)(|E|/M³)^1/2 = 0.039, which is motivated by cosmological N-body experiments (Bullock et al. 2001). The simulation closely follows that described in Roškar et al. (2008) and is evolved with the parallel N-body+SPH code GASOLINE (Wadsley et al. 2004) for 10 Gyr. Table 1 summarizes some of the key properties of the different simulations used in this paper.

Table 1. Details of Simulations

Run	rvir	Mvira	Mbb	DMc	hd	Run Description	Referencee
Name	(kpc)	(1012 M☉)	(1011 M☉)	(pc)	(kpc)
g15784	290.8	1.4	2.1	312.5	3.38	MUGS disk galaxy	Stinson et al. (2010)
SA1	491.9	6.9	1.75	100	3.0	Triaxial halo+short-axis stellar disk	D08, V12
SBgs	317.9	1.9	0.9f	100	1.9	Prolate halo+hot gas → disk w/SF	Debattista et al. (2013)

Notes. ^aM_vir: mass within the virial radius (r_vir). ^bM_b: mass in baryons. ^c_DM: gravitational softening length. ^dh: exponential radial scale length of disk. ^eReference where the simulation is described. ^fMass of baryons in stars only.

Download table as:  ASCIITypeset image

In each simulation, we selected 1–2 × 10⁴ halo particles, randomly distributed within some volume. The "global sample" refers to particles selected at random within a spherical radius r_g centered on the model's galactic center (potential minimum).

In the MUGS galaxy g15784, we identify halo stars by excluding stellar orbits associated with the massive central spheroid (orbits with apocenter radius r_apo < 5 kpc) and stellar disk (orbits whose maximum distance from the disk plane Z_max < 3 kpc). The global sample comprised ∼10⁴ halo stars with r_g < 50 kpc. Since the controlled simulations do not form stellar halos, we follow the orbits of dark matter particles selected with r_g < 200 kpc in these two simulations.

In all three, potentials orbits were integrated forward from their z = 0 initial conditions for a duration of ∼50 Gyr in a frozen potential resulting from the full mass distribution of the simulation including dark matter and baryons using an integration scheme based on the PKDGRAV tree. The assumption made here (and by most studies that involve the kinematics of halo stars in Milky Way) is that the potential of the galaxy at z = 0 is nearly in steady state. While this assumption is not strictly true, there is evidence to support the view that halos of mass similar to the Milky Way evolve less at the present epoch than they have in the past. Busha et al. (2007) show that halos with M₂₀₀ ≲ 10¹³ M_☉ in ΛCDM simulations evolved for 64 h⁻¹ Gyr after the big bang (see their Figure 8) grow rapidly till t/h = 10 (corresponding to 13.8 Gyr after the big bang), after which there is little change in their virial mass. This does not mean that secular evolution in the disk and halo stops at z = 0, merely that the rate of mass accretion slows down. We emphasize that the long orbital integration times only serve to improve the accuracy with which orbits in the potential at z = 0 can be characterized. We do not imply that the potential will have been static for the duration of integration or that the long integration time is physically meaningful.

In this section, we provide a brief description of the Laskar's frequency analysis method. A more detailed description of the application to orbits in N-body simulations is given in V10 and V12.

Orbits in galaxies are approximately quasi-periodic (Binney & Tremaine 2008; hereafter BT08), and hence their space and velocity coordinates can be represented by a time series of the form: $x(t) = \sum _{k=1}^{k_{{\rm max}}} A_{k} e^{i\omega _kt}$, and similarly for y(t), v_x(t), etc. This means that a Fourier transform of such a time series will yield the spectrum of frequencies, ω_k, and amplitudes, A_k, that define the orbit (Binney & Spergel 1982, 1984). In a three-dimensional potential, only three frequencies in the spectrum are linearly independent. The remaining ω_k are given by ω_k = l_kΩ₁ + m_kΩ₂ + n_kΩ₃. For a regular orbit, the frequencies Ω_i, i = 1...3 are constant "fundamental frequencies" that are related to three angle variables θ_i(t) = Ω_it and three integrals of motion (more precisely actions J_i; Binney & Tremaine 2008). Since the full phase-space trajectory in the action-angle coordinates (J_i, θ_i) resembles a three-dimensional torus, the mapping from (x(t), v(t))⇔(J, θ) is referred to as torus construction. Laskar (1990, 1993) developed a very accurate numerical technique, "Numerical Analysis of Fundamental Frequencies," to recover frequencies by taking Fourier Transforms of complex time series of the form x(t) + iv(t). We use the implementation of Valluri & Merritt (1998), which uses integer programming to recover orbital fundamental frequencies with extremely high accuracy in ∼20–30 orbital periods (significantly faster than 100s of orbital periods needed by other codes; e.g., Carpintero & Aguilar 1998). Orbits were integrated and orbital frequencies computed in Cartesian coordinates and in cylindrical coordinates, as described in V12.

Orbital fundamental frequencies can be used to (1) classify orbits into major orbit families (Carpintero & Aguilar 1998); (2) quantify the elongation of an orbit relative to the shape of the potential (V10); (3) identify chaotic orbits (Laskar 1993); (4) identify important resonant orbit families (Robutel & Laskar 2001) which have trapped a large number of resonant orbits yielding insights into the importance of secular or adiabatic processes; and (5) determine the self-consistent equilibrium potential from which a population of orbits was drawn (V12).

To determine the fraction of chaotic orbits in a potential, the orbital time series is divided into two equal segments and the orbital fundamental frequencies computed in each time segment. Since regular orbits have fixed frequencies that do not change in time, the change in the frequency measured in the two time segments can be used to measure the drift in frequency space. V10 showed that even for orbits in N-body potentials (which are inherently noisy), it is possible to distinguish between N-body jitter and true chaos by a quantitative measurement of frequency drift by defining a diffusion rate parameter log (Δf) which measures the logarithm of the change in the frequency of the leading term in the orbit's frequency spectrum, measured in two consecutive time segments. V10 showed, using orbits in N-body simulations of spherical NFW halos with gravitational softening = 100 pc, that orbits with log (Δf) < −1.2 were regular. In g15784, the fraction of orbits with log (Δf) > −1.2 (chaotic by the criterion established in the spherical NFW halo) is over 50%. However, visually, many of the orbits have almost regular appearances. g15784 has a larger gravitational softening ( = 312.5 pc) which appears to result in a slightly smoother potential, and hence orbits with log (Δf) < −0.5 appear more "regular" (i.e., they fill a well defined volume in configuration space) than they do in the controlled simulations. We therefore use this more relaxed criterion for identifying regular orbits in the cosmological simulation, but it must be kept in mind that in N-body simulations, all orbits are inherently "noisy" resulting in a continuous distribution of diffusion parameters log (Δf) and any cutoff between "regular" and "chaotic" is arbitrary. Furthermore, we have found that decreasing the cutoff to log (Δf) < −1.2 does not significantly affect the fractions of short- and long-axis tubes, but decreases the fraction of box orbits while increasing the fraction of chaotic orbits. We therefore caution readers that in this cosmological simulation, the boundary between chaotic and regular box orbits is rather fuzzy and although we distinguish between them for consistency with our previous work, they are better thought of as a single family of "centrophilic" orbits with no net angular momentum.

3.1.1. Halo Shape

We computed the shapes of the spatial matter distributions of dark matter and star particles in g15784 within ellipsoidal shells centered on the minimum of the potential. The shape is determined at each radius r by finding the ellipsoidal shell of width ∼5 kpc whose principal axes a, b, and c self-consistently (1) have a geometric mean radius $\sqrt{abc}=r$, and (2) are the eigenvectors of the second moment tensor of the mass distribution within the ellipsoidal shell (J. Bailin et al. 2013, in preparation). In practice, this is done by starting with a spherical shell, calculating the second moment tensor, deforming the shape of the shell to match the eigenvectors of the tensor, and iterating until the solution has converged; this is very similar to the method advocated by Zemp et al. (2011).

To measure the shape of the stellar halo (at z = 0), we select halo stars by excluding all stars within 5 kpc of the center of the galaxy (defined as the "bulge"), as well as all stars which lie within 3 kpc of the disk plane (defined as a thin+thick disk) for particles within 20 kpc (the radius within which most disk particles reside). This cut to identify "halo stars" is not perfect since it incorrectly excludes the small fraction of halo stars that happen to lie close to the disk plane or within 5 kpc of the galactic center at z = 0. This cut may potentially bias c/a values to be somewhat higher than they should be within the 20 kpc inner region of the halo, but is unlikely to significantly alter b/a. We use a more sophisticated method (relying on kinematic properties) to exclude disk stars when we examine the orbital properties of stars.

Following standard nomenclature a, b, and c are defined as the semi-axis lengths of long (x), intermediate (y), and short (z) axes, respectively. The short axis is perpendicular to the plane of the disk in the inner region (the disk warps beyond about 30 kpc). Figure 1 shows axis ratios as a function of radius b/a (top panel) and c/a (middle panel), and the triaxiality parameter T = (1 − b²/a²)/(1 − c²/a²) (bottom panel). The triaxiality parameter T is often used to characterize the shape of ellipsoidal figures that are termed "oblate" if T < 1/3, "triaxial" if 1/3 < T < 2/3, and "prolate" if T > 2/3. Thin dashed horizontal lines marking these shape regimes are marked in the bottom panel of Figure 1.

Zoom In Zoom Out Reset image size

For the dark matter halo, b/a and c/a vary only slightly with radius (1 < b/a < 0.9, c/a ∼ 0.8). The bottom plot shows that the halo is oblate within ∼15 kpc, becoming increasingly triaxial to prolate-triaxial beyond this radius. The shape of the stellar halo is highly oblate (T ∼ 0.1) within 40 kpc, becoming mildly triaxial at larger radii. The stellar halo is also much flatter (c/a is smaller) than the dark matter halo at all radii. Note that at radius ∼17 kpc, the sharp dip in b/a and the corresponding spike in T are the result of a subhalo of mass ≲ 10⁷ M_☉ that was not identified by the Amiga halo-finder. Zemp et al. (2011) show that when subhalos are not subtracted (or improperly subtracted), they produce spurious spikes in the shape parameters. The overall changes in the shape parameters (e.g., T) of the dark matter halo and stellar halo as a function of radius are consistent with the results of previous studies (e.g., Zemp et al. 2012).

3.1.2. Orbital Type Distributions

To assess how orbital structure depends on shape, we examined orbits of ∼10⁴ star and dark matter particles in the global sample of the galaxy g15784, and for comparison, 10⁴ dark matter orbits from each of the two controlled simulations SA1 and SBgs. Orbits were classified into four types (box orbits, SATs, LATs, and chaotic orbits). Their distributions as a function of orbital pericenter radius (r_peri) and orbital eccentricity were examined in order to understand how the shapes of the two components (dark matter and halo stars) are related to the types of orbits that self-consistently reproduce their shapes.

Figure 2 shows orbit populations as a function of r_peri in the two controlled simulations SA1 and SBgs (top panels), and in g15784 (bottom panels). We choose to plot the orbit fractions as a function of r_peri rather than some other measure of the "radius" of an orbit (such as the time-averaged radius) because we are specifically interested in understanding if the SATs contribute to the oblateness in the inner halo. Since r_peri defines the inner boundary of a tube orbit, it defines the region interior to which the orbit contributes no mass. Hence, this quantity is best suited to addressing this question. In theory, box orbits and chaotic orbits can travel arbitrarily close to the center of the potential and so, formally, they have r_peri = 0 but, in practice, r_peri > 0 during the 50 Gyr (and at least 20 orbital period) integration time. Each panel shows kernel density distributions of each of the four orbit families as indicated by the line legends. The area under each curve is proportional to the total number of orbits of that family. It is important to note that in the top two panels, the gravitational softening length of dark matter particles is ∼0.1 kpc, while in the bottom panels the gravitational softening length is ∼0.3 kpc; distributions at pericenter radii smaller than the softening lengths are not reliable.

Zoom In Zoom Out Reset image size

The progenitor halo of SA1 was triaxial with 86% of particles on box orbits and 11% on LATs (V10). Figure 2 shows that this population was transformed (due to the growth a stellar disk) to a halo that is still dominated by boxes (49%) but the fraction of SATs (32%), LATs (12%), and chaotic orbits (7%) has increased significantly. In contrast, the progenitor halo of SBgs was prolate, with 78% LATs and 15% box orbits (V10). Following the condensation of hot gas and the growth of a stellar disk, the dark matter halo orbit population evolves to an overall potential in which SATs dominate (51%) with LATs (35%) also being important. In the cosmological simulation g15784, halo star particles (Figure 2, bottom right) are primarily on SATs (47%) and LATs (29%) and dark matter particles also have similar orbit fractions (40% are on SATs and 36% are on LATs).

Table 2 shows the orbit fractions in the inner and outer regions of all three halos following the growth of disks. We computed the relative fractions of different types of orbits with r_peri less than or greater than 20 kpc. Since the inner ∼20 kpc region is where the disk dominates, this is the region where the shape of the halo is close to oblate axisymmetric in all three potentials. In none of the halos do SATs constitute more than 50% of the orbit fraction in the oblate inner halo (r_peri < 20 kpc). In all four models, the "triaxial orbit families" (box orbits and chaotic orbits and LATs) together constitute 50% or more of the population. In SBgs and g15784, the outer parts of the halo (r_peri > 20 kpc) are mildly triaxial. Here, we find that LATs are the dominant population. In SA1—which is quite strongly triaxial beyond 20 kpc—SATs constitute 65% of the population, with LATs constituting 33%. Thus, in no case do we find evidence that an oblate density distribution implies a distribution function dominated by axisymmetric SATs, and conversely, we find that the presence of a large fraction of SATs does not imply axisymmetry.

Table 2. Orbit Fractions in Inner and Outer Halo Regions

Run	rperi ⩽ 20 kpc				rperi > 20 kpc
	Box	Chaotic	LAT	SAT	Box	Chaotic	LAT	SAT
SA1	0.51	0.09	0.11	0.29	0.01	0.05	0.33	0.65
SBgs	0.10	0.10	0.30	0.50	0.03	0.02	0.52	0.43
Stars	0.12	0.13	0.32	0.43	0.03	0.02	0.55	0.41
DM	0.12	0.12	0.34	0.42	0.09	0.00	0.59	0.32

Download table as:  ASCIITypeset image

Although the distributions of orbit types as a function of r_peri differ in detail for dark matter particles and halo stars in g15784, their overall distributions are significantly more similar to each other than they are, for instance, to the orbit distribution of dark matter particles in model SA1 (Figure 2, top left). Furthermore, g15784 shows a striking resemblance to the orbit populations in SBgs. The similarity is probably coincidental in view of the fact that g15784 formed in a cosmological context and has experienced significant perturbations due to mergers and substructure. Minor mergers are entirely absent from the history of SBgs, which formed from a single major merger between two spherical halos. We expect the halos of galaxies of the mass of g15784 to have experienced multiple major mergers and to therefore be more like SA1, rather than SBgs which has only had one major merger. It appears that since g15784 had its last major merger nearly 10 Gyr ago, this galaxy's stellar halo was primarily formed via minor majors (rather than major mergers). Since the orbits at large radii are primarily tubes, we hypothesize that the satellites that formed the stellar halo were mainly accreted on tube orbits. We will see in Section 3.3 that there is support for this hypothesis.

Recently, Bryan et al. (2012) analyzed orbits in halos of mean mass ∼6 × 10¹³ M_☉ (at z = 0) and 7 × 10¹¹ M_☉ (at z = 2) from the OverWhelmingly Large Simulations (OWLS) to investigate the effects of various feedback prescriptions on the orbital properties of stars and dark matter particles in a cosmological context. Bryan et al. (2012) selected 50 halos from 5 different simulations (with different feedback prescriptions) and selected 500 particles per halo. They integrated orbits in gravitational potentials computed with the Self-Consistent Field (SCF) method and classified them using the method of Carpintero & Aguilar (1998).

Bryan et al. (2012) find that as the central baryonic fraction increases, the fraction of box orbits in the inner part of the halo decreases. Since we only examine one simulation, we are unable to verify their result that the fraction of box orbits in the inner region decreases with increasing central baryonic mass fraction; however, the baryonic mass fraction in the inner 20 kpc region of g15784 (where the disk dominates) is obviously much higher than in the outer region, and we find no evidence that box orbits have been depleted in the inner regions as a consequence of baryonic condensation. The difference between our result and those of Bryan et al. (2012) is most probably a consequence of differences in the way in which the total galactic potential (in which orbits are evolved) is computed. The SCF method (used by Bryan et al. 2012) uses a low-order multipole expansion code to compute the potential, which consequently is much smoother and has more symmetry than our potential, which is computed with the PKDGRAV tree-code and therefore includes all the irregularities of the full cosmological simulation (except subhalos more massive than 10⁷ M_☉ which were removed).

3.1.3. Orbital Shapes: Elongation versus Eccentricity

In any self-consistent potential, the shapes of the majority of the orbits must match the overall shape of the three-dimensional matter density distribution. In a triaxial potential, the elongation along the major axis is provided either by box orbits or inner LATs. The ratios of the fundamental frequencies of orbits can be used to characterize their overall shape. For a triaxial potential, the fact that semi-axes lengths a > b > c implies that the oscillation frequencies obey the condition |Ω_x| < |Ω_y| < |Ω_z| for any (non resonant) orbit with the same over-all shape as the density distribution (we consider only the absolute values of the frequencies since their signs only signify the sense of oscillation). V10 use this property to define an average "elongation" parameter χ_s for any orbit. When an orbit is elongated in the same way as the potential,

$$\begin{eqnarray*} & & |\Omega _z| > |\Omega _y| > |\Omega _x| \\ &\Rightarrow & {\frac{|\Omega _y|}{|\Omega _z|}} > {\frac{|\Omega _x|}{|\Omega _z|}}. \end{eqnarray*}$$

Therefore, we can define an elongation parameter χ_s by

$$\begin{eqnarray} \chi _s &\equiv & {\frac{|\Omega _y|-|\Omega _x|}{|\Omega _z|}}, \end{eqnarray} \tag{ 1 }$$

such that χ_s > 0 for orbits elongated along the major axis, with larger values of χ_s implying a greater the degree of elongation (V10). V12 showed that some systems can be elongated along the y-axis at small radii but elongated along the x-axis at larger radii; also, some outer LAT orbits can have a greater extent along the y-axis than along the x-axis. In both these cases, χ_s is slightly negative. Orbits for which all frequencies are almost equal enclose a volume that is almost spherical have χ_s ∼ 0 because Ω_x ∼ Ω_y ∼ Ω_z. Also, orbits that are close to axisymmetric about the short (z) axis (i.e., SATs) have Ω_x ∼ Ω_y, and hence χ_s ∼ 0, regardless of the value of Ω_z.

Figure 3 shows the kernel density distribution of orbits of the four different types as a function of orbital elongation χ_s. SA1 shows a significantly larger fraction of elongated orbits (χ_s > 0.1) than the other models. The SATs (dashed lines) have χ_s = 0 because they are close to axisymmetric. Model SBgs (which is prolate at large radii) has a significant fraction of elongated LAT orbits. In contrast, for g15784, all four families of orbits for both dark matter (bottom left) and stars (bottom right) have very small elongation (i.e., they are quite "round" or axisymmetric).

Zoom In Zoom Out Reset image size

V10 showed for several controlled simulations (including SA1) that the orbital elongation parameter χ_s in the inner regions decreased from χ_s ∼ 0.4–0.5 in the triaxial or prolate models to χ_s ∼ 0.2 after baryonic components were grown. They showed that although some orbits do change from boxes or LATs in the prolate systems to SATs, all orbits adapt to the more axisymmetric baryonic component by becoming less elongated, especially at small r_peri. However, only a small fraction of all orbits (between 4% and 25%) actually transform themselves to become SATs. Figure 4 shows the median value of χ_s in 15 bins in r_peri for orbits in the three simulations. The top two panels show that when simulations are controlled (adiabatic), the baryons make box and chaotic orbits with smaller pericenter radii that are less elongated, but at larger pericenter radii these orbits remain elongated (LATs always tend to be elongated but as shown by V10, these too are more elongated when baryons are absent). The lower two panels show the elongations of orbits in g15784. In this simulation, all types of orbits are essentially "round" over the entire range in r_peri. Although this confirms the findings of V10 that the orbits of all three families can become "round" to support the shape of the potential, it is also clear that overall, the shapes of all the orbits are similar and close to axisymmetric at all radii.

Zoom In Zoom Out Reset image size

Another commonly used metric for characterizing individual orbits is the eccentricity parameter:

$$e = {\frac{r_{\rm apo}- r_{\rm peri}}{r_{\rm apo}+r_{\rm peri}}}.$$

Strictly speaking, orbital eccentricity e is defined for a two-dimensional orbit in terms of its apocenter radius r_apo and pericenter radius r_peri. It is important to emphasize that eccentricity does not measure the degree of axisymmetry; e.g., a perfectly planar orbit in a spherical potential is always an axisymmetric rosette, but can have an arbitrary eccentricity that depends on its angular momentum. At a given energy, orbits with very high angular momentum have an eccentricity approaches zero, while orbits with very low angular momentum have eccentricity approaches unity. Figure 5 shows the kernel density distributions of orbital eccentricity for particles in the three simulations. We see that all four orbit families have larger fractions of orbits with high eccentricity than with low eccentricity. Even SATs (dashed curves) and LATs (dot-dashed curves) are very eccentric, implying that they have fairly low angular momentum. The absence of a low eccentricity box and chaotic orbits is not unexpected—these orbits are highly radial (formally having zero time averaged angular momentum) and have r_peri approaching zero. The average orbital eccentricities of halo stars in each orbit family are = 0.88 (boxes), = 0.95 (chaotic), = 0.76 (SATs), and = 0.72 (LATs). For dark matter particles, the average orbital eccentricities of each orbit family are = 0.82 (boxes), = 0.95 (chaotic), = 0.71 (SATs), and = 0.63 (LATs). These average orbital eccentricities for the various orbit families are somewhat larger than the average orbital eccentricities ( = 0.6) of dark matter subhalos and particles in dark-matter-only simulations (van den Bosch et al. 1999). It is particularly interesting to reflect on the fact that for g15784 halo stars, even those on SATs are more likely to have high eccentricity than low eccentricity. If a significant fraction of the stellar halo was formed in situ in an early disk that was then heated to form the halo, then one would expect this population to have a higher angular momentum and lower eccentricity than an accreted population. The fact that the majority of the stellar halo orbits in this global sample have a relatively high eccentricity reflects the fact that they were primarily accreted.

Zoom In Zoom Out Reset image size

The main conclusions from Figures 2–5 are as follows.

• 1.  
  The types of orbits that characterize a halo (either in a controlled simulation or in a cosmological simulation) reflect both the merger/formation history and the effect of baryons.
• 2.  
  All four types of orbits (boxes, SATs, LATs, and chaotic) can be present in significant proportions even when the shape of the potential is nearly oblate.
• 3.  
  All four types of orbits are able to become less elongated (i.e., more axisymmetric) to adapt to the nearly oblate shape
• 4.  
  All four types of orbits have a high fraction of eccentric orbits reflecting the fact that the halos are products of mergers or radial infall.

Since the mass of the stellar halo is only a tiny fraction of the mass of the dark matter halo, the dynamics of halo stars are determined by the shape of the dark matter halo, the orbital initial conditions at the time of accretion, and the evolution of the halo potential since the epoch of accretion (Knebe et al. 2005). All of these factors also influence the orbital types, orbital shapes and eccentricity distributions of orbits. The similarities in the distributions of orbits of dark matter and halo stars with r_peri, χ_s and eccentricity (Figures 2–5) show qualitative similarities that are quite significant, especially compared with the orbit populations of the controlled simulation SA1, which formed from multiple major mergers. Bryan et al. (2012) also found in their analysis of orbits from the OWLS cosmological simulations that star particles and dark matter particles have similar orbital distributions as a function of radius.

V12 showed that a useful way to map the entire distribution function of a self-consistent system was via frequency maps. A frequency map is obtained by plotting the ratios of fundamental frequencies: in Cartesian coordinates, Ω_x/Ω_z versus Ω_y/Ω_z, or in cylindrical coordinates, Ω_ϕ/Ω_R versus Ω_ϕ/Ω_z. They showed that frequency maps give a pictorial representation of the different orbit families, their relative importance, and how they are distributed with energy. In the frequency maps in Figure 6, each point represents a single orbit with color representing the orbital binding energies in 3 bins (blue: the 1/3rd most tightly bound; red: 1/3rd least tightly bound; green: intermediate energy).

Zoom In Zoom Out Reset image size

Figure 6 shows the frequency maps for orbits in galaxy g15784. The top left panel shows a map (in Cartesian coordinates) for dark matter particles, while the bottom left shows the same for halo stars. SATs lie along the diagonal Ω_x/Ω_z ∼ Ω_y/Ω_z. The LAT family is clustered along the horizontal line Ω_y/Ω_z ∼ 1, box and chaotic orbits are mostly the blue points scattered around the map. The dashed lines mark the major resonances (satisfying conditions like lΩ_x + mΩ_y + nΩ_z = 0), with numbers representing the integers l, m, n. In this Cartesian representation, no orbits are seen to be associated with box-orbit resonances (diagonal lines on the left-side of the graph). The right hand panels of Figure 6 show the frequency maps in cylindrical coordinates for the same sets of orbits. A prominent resonance is seen at Ω_z/Ω_R = Ω_ϕ/Ω_R ∼ 0.5 in orbits of both the halo stars and dark matter particles. V12 showed that numerous resonances characterize the phase space of a self-consistent distribution function evolved adiabatically in controlled simulations. This resonance is associated with a thin shell resonance of the SAT family and it is impressive that such resonant orbits persist despite the hierarchical accretion history of this system. The frequency map representations of the halo stars and dark matter particles are very similar to each other both in Cartesian maps (left) and in the cylindrical maps (right), providing additional evidence that their phase space distribution functions are similar. This implies that the phase space distribution function of this sample of stellar halo orbits probably originated in a manner similar to that of the dark matter orbits.

Since the prescient work of Eggen et al. (1962), it has been recognized that correlations between orbital properties of halo stars and their metallicities and ages can be used to infer the formation history of the Galaxy. If the orbital integrals of motion of halo stars (e.g., total energy E, the total angular momentum L, and angular momentum about the z-axis L_z) change only slowly as the potential evolves and the debris disperses, then correlations in integral-of-motion space can lead to the recovery of the progenitor satellite and also enable recovery of the underlying potential from tidal debris even in time-dependent potentials (Johnston et al. 1996; Helmi & de Zeeuw 2000; Bullock & Johnston 2005). Searching for correlations in angle-action or frequency space may improve the accuracy of such recovery (McMillan & Binney 2008; Gómez & Helmi 2010). Using semi-analytic prescriptions to construct stellar populations for accreted halos from cosmologically motivated pure N-body simulations, Font et al. (2006) showed that combining phase space information with chemical abundances and stellar ages can aid the recovery of individual satellites, even when the stars in a satellite have a small spread in ages and chemical characteristics. To date, there have been no studies that have tried to identify halo stars associated with individual satellite progenitors in a cosmological hydrodynamical simulation.

Figure 7 shows orbital energy versus one component of angular momentum for halo stars from the global sample in each of the four orbit families studied in Section 3.1.2. For box orbits, SATs, and chaotic orbits, the abscissa is the angular momentum about the z axis (L_z), while for LATs the abscissa is L_x since x is the axis about which this tube family rotates. Points corresponding to orbits in each family are fairly smoothly distributed and almost no substructure typical of tidal debris associated with distinct individual satellites (e.g., Helmi & de Zeeuw 2000) is seen. The ∼10⁴ particles plotted in Figure 7 were selected at random within a spherical volume of r_g = 50 kpc. However, even a local sample selected within a small volume around the location of a fiducial "sun" showed no substructure in E, L, L_z, L_x space. It is important to emphasize that we only plot E versus L for comparison with previous work, but that in fact none of the components of angular momentum are expected to be integrals of motion for any of the orbit families in this triaxial potential. Figure 7 clearly shows that in this simulation, the hierarchical growth of the halo, accompanied by dissipative evolution, has erased all correlations between these pseudo integrals of motion. It should be noted that our choice of orbital initial conditions (randomly chosen within 50 kpc of the galaxy's center) is not ideal for searching for coherent tidal debris. Typically for such studies, a local volume representing the solar neighborhood is examined. Our examination of such a volume (i.e., R_s = 4 kpc) yielded too small of a fraction of bona fide halo stars (because the disk in this simulation is quite thick) to perform this exercise. We also searched for substructure in plots of energy versus metallicity, energy versus stellar age, and angular momentum versus metallicity and stellar age. None of these plots revealed any detectable substructure that would enable us to identify satellite progenitors.

Zoom In Zoom Out Reset image size

J. Bailin et al. (2013, in preparation) analyze the age-metallicity relation for all stars (halo and disk) in the MUGS cosmological disk galaxy simulation g15784. They find that while the most metal-poor stars in the galaxy were formed at the earliest times, the ISM is rapidly enriched, and hence the formation of highly metal enriched stars follows very quickly thereafter. Bailin et al. show that when [Fe/H] is plotted versus formation time, numerous tight "stream-like" features that correspond to continuous star formation within individual galactic subunits, are seen. A study of the relative spatial, metallicity, and age distributions of stars that were accreted onto g15784 and stars that formed in the disk and were subsequently kicked into the halo will be presented in Bailin et al. Here, we focus only on 10⁴ halo stars whose orbits were analyzed in previous sections, and examine how their orbital characteristics depend on their metallicities and ages.

In Figure 8, we plot total metallicity (in solar units) for stars in the halo orbit sample studied in previous sections, versus their formation time in the simulation (t_form). The formation time is scaled such that stars that form at the present time (z = 0) have t_form = 13.7 Gyr. The panels shows stars on different types of orbits: boxes, SATs, LATs, and chaotic orbits as indicated by the labels. Each point represents a star particle from the simulation.

Zoom In Zoom Out Reset image size

Most striking are the thin stream-like features (identified by Bailin et al.) which arise from progressive enrichment and star-formation inside individual sub-galactic units. Note that even though we have randomly selected ∼10⁴ halo stars within 50 kpc, with no other constraints other that they have r_apo > 5 kpc and Z_max > 3 kpc the thin streams associated with individual progenitors are clearly visible. Star formation appears to stop at some point in all the streams. While more detailed analysis with all the MUGS snapshots is needed to determine the precise reason for this, a plausible interpretation is that star formation in a subhalo can continue for some time after it is accreted by the galaxy but stops once it is tidally disrupted. Consequently, the end of the "stream" serves as a proxy for when it is significantly disrupted by tidal forces from the main galaxy.

The thin stream-like features which were disrupted early (t_form < 8–10 Gyr) appear in all four panels, indicating that at z = 0, stars associated with a single progenitor satellite traverse the halo on all four types of orbits. For example, the features labeled "A" consists of stars with a wide range of metallicity and t_form and are clearly seen in four panels. This implies that these two galactic sub-units were completely disrupted by tidal forces and then scattered onto different types of orbits by mixing processes.

We believe that this is probably "chaotic mixing" rather than "phase mixing" since the latter conserves orbital integrals of motion in a static potential and changes them adiabatically in a slowly varying one. Phase mixing is not capable of changing the orbital type, but chaotic mixing is (Merritt & Valluri 1996). Therefore, it is unlikely that the phase space coordinates of stars (at z = 0) in the two units labeled "A" can be used to fully reconstruct the progenitor or its orbital properties at the time it was accreted. For dynamically older streams like "A," the integrals of motion, especially energy, are expected to change significantly due to the change in the halo potential due to mass accretion, both dark (Knebe et al. 2005) and baryonic (V10). Additionally, chaotic mixing in the triaxial halo causes tidal streams to disperse more rapidly (Vogelsberger et al. 2008).

In contrast, the t_form-metallicity feature labeled "B" is most clearly seen in the top right panel (SATs). Its stars also span a range of metallicities (−1.5 < [Z/Z_☉] < −0.5) and formation epochs (6–13.5 Gyr); however, most of the stars of this satellite are on SATs (although a faint hint of it is seen in the LATs panel), while there is no distinct counterpart in the box or chaotic orbit panels.⁹ This indicates that the orbits of this satellite are not mixed, and hence the progenitor's phase space coordinates and the halo potential can probably be extracted by backward evolution of stellar orbits. The stream labeled "C" in the bottom left panel consists of stars on LATs and is also only seen in this panel. The near-exclusivity in orbit-type for streams "B" and "C" implies that the orbits of these stars have not experienced much mixing (probably because they were accreted/tidally disrupted more recently). Figure 8 shows that the thin steam-like features have stars with a wide range of metallicities and star formation occurs over an extended period of time. We see that stars with a range of metallicities and ages that were once associated with the same progenitor are easy to identify because they form thin features in metallicity-formation time plots despite sometimes having very different orbital characteristics.

Our preliminary analysis of the age-metallicity signatures of our selected halo stars in Figure 8 shows that satellites leave the stars in long and short-axis tube orbits as they accrete. Satellites that accreted >5 Gyr ago have chaotically mixed these stars into chaotic and box orbits, but more recently accreted satellites have not yet become well mixed. Such chaotic mixing of early accreted satellites is consistent with previous studies of satellite accretion (Bullock & Johnston 2005; Font et al. 2006; Gómez & Helmi 2010).

It is important to note that star formation appears to continue in this halo population until quite late times—well beyond what is expected based on observations of the stellar halos of the Milky Way and M31. This is in part due an artifact of the detailed sub-grid physics in the simulations: the cooling of gas in to dark matter subhalos results in their becoming more massive than in the real Universe. This also allows subhalos in these simulations to continue to form stars after they are accreted by the main galaxy. This (presumed) artifact is not unique to our simulations and has been previously observed (Tissera et al. 2012, 2013). Tissera et al. (2013) refer to this component of the stellar halo as "endo-debris," and distinguish it from "debris" stars which were formed in a satellite before it was accreted. There is tentative evidence that such a population has also been discovered in the Milky Way halo by Sheffield et al. (2012), who have identified a population of halo stars with low [Fe/H] but high [α/Fe], which they argue could have been formed in situ in the halo (i.e., inside satellites after they were accreted). Tissera et al. (2013) also suggest that one population of Carbon Enhanced Metal Poor stars found by Carollo et al. (2012) could be such "endo-debris" stars. If star formation in satellites terminates earlier than in these simulations, then narrow features would not extend to late times, but would still remain detectable at early times as in Figure 8.

Together, Figures 7 and 8 imply that the tidal debris in cosmological hydrodynamical simulations experiences significantly more chaotic evolution than in collisionless simulations, making it much harder to identify individual progenitors using phase space coordinates alone. But the identification of distinct progenitors is quite easy in age-metallicity plots, and such plots when combined with orbital information can be used to probe the accretion history of the halo. Current cosmological hydrodynamic simulations, including the one studied here, do not have adequate resolution to identify the types of stellar streams observed in the Milky Way and M31. The study of the correlations between kinematics, metallicities, and spatial distributions of halo stars and their formation history will become increasingly important in the future as the resolution of simulations and the baryonic physics prescriptions improve.