# The shape of dark matter subhalos in the Aquarius simulations

###### Abstract

We analyze the Aquarius simulations to characterize the shape of dark matter halos with peak circular velocity in the range , and perform a convergence study using the various Aquarius resolution levels. For the converged objects, we determine the principal axis () of the normalized inertia tensor as a function of radius. We find that the triaxiality of field halos is an increasing function of halo mass, so that the smallest halos in our sample are rounder than Milky Way-like objects at the radius where the circular velocity peaks, . We find that the distribution of subhalo axis ratios is consistent with that of field halos of comparable . Inner and outer contours within each object are well aligned, with the major axis preferentially pointing in the radial direction for subhalos closest to the center of their host halo. We also analyze the dynamical structure of subhalos likely to host luminous satellites comparable to the classical dwarf spheroidals in the Local Group. These halos have axis ratios that increase with radius, and which are mildly triaxial with and at kpc. Their velocity ellipsoid become strongly tangentially biased in the outskirts as a consequence of tidal stripping.

###### keywords:

methods: numerical - galaxies: dwarf - cosmology: dark matter.^{†}

^{†}pagerange: The shape of dark matter subhalos in the Aquarius simulations–Acknowledgements

^{†}

^{†}pubyear: 2019

## 1 Introduction

In the cold dark matter cosmological paradigm structures build hierarchically, via the mergers of smaller objects (Press & Schechter, 1974; Gott & Rees, 1975; White & Rees, 1978; Blumenthal et al., 1984). As mergers proceed, the innermost regions of some of the progenitors survive, resulting in non-linear structures where a wealth of substructure orbits the center of an otherwise monolithic dark halo. Early -body simulations showed that halos could host dozens of substructures, down to masses near the numerical resolution limit (Tormen, 1997; Tormen, Diaferio & Syer, 1998; Ghigna et al., 1998; Klypin et al., 1999a, b; Moore et al., 1999). For systems like the Milky Way, current numerical simulations have extended the dynamical range of resolved substructures by 4-5 orders of magnitude (Diemand, Kuhlen & Madau, 2007; Diemand et al., 2008; Springel et al., 2008; Stadel et al., 2009).

The properties of these substructures are of great interest since luminous satellites, such as the population of dwarf spheroidal (dSph) galaxies in the Local Group, are expected to be embedded in them (Stoehr et al., 2002; Strigari et al., 2007; Boylan-Kolchin, Bullock & Kaplinghat, 2012; Vera-Ciro et al., 2013). Furthermore, the large mass-to-light ratios of dSph, which range from 10 s to 1000 s (Mateo, 1998; Gilmore et al., 2007; Walker, 2013), indicate that their internal dynamics is dominated by the dark matter. This suggests that the predictions of pure dark matter simulations may be directly confronted with observations of these systems. For instance, it has been suggested that they provide an optimal place to look for signals of dark matter self-annihilation processes (Kamionkowski, Koushiappas & Kuhlen, 2010) due to the natural enhancement in density and the lack of significant contamination from the baryonic component.

The availability of large samples of line-of-sight velocities for individual stars in dSph galaxies offers new tests of the predictions of cold dark matter (CDM) regarding the structure of dark matter subhalos. For instance, studies of -body numerical simulations have shown that the inner slope of the dark matter density profile is expected to be cuspy in CDM models (Navarro, Frenk & White, 1996, 1997). This seems to contrast with the somewhat shallower slopes and even constant density cores proposed to explain the motions of stars in local dwarf spheroidals (Amorisco & Evans, 2011; Walker & Peñarrubia, 2011; Amorisco & Evans, 2012; Jardel & Gebhardt, 2012). This, however, is a subject of active debate, since various authors have shown that the stellar kinematics of Milky Way dwarfs are also consistent with the NFW cuspy profiles (Battaglia et al., 2008; Walker et al., 2009; Strigari, Frenk & White, 2010; Breddels et al., 2013; Breddels & Helmi, 2013).

Most of the dynamical modeling performed in the studies of Local Group dwarf spheroidals relies on simple assumptions about the structure of their dark matter component. In particular, spherical symmetry and specific anisotropy profiles have been extensively assumed. The orbital anisotropy has been taken to be constant (Richstone & Tremaine, 1986; Łokas, 2002; Łokas, Mamon & Prada, 2005; Łokas, 2009; Walker et al., 2009), or radially dependent (Kleyna et al., 2001; Wilkinson et al., 2002; Battaglia et al., 2008; Strigari et al., 2008; Wolf et al., 2010; Amorisco & Evans, 2011), while in Schwarzschild modeling or in made-to-measure -body methods it does not need to be assumed (Long & Mao, 2010; Jardel & Gebhardt, 2012; Breddels et al., 2013). The selection of geometric shape for the dark matter potential can also be relaxed. For example, Hayashi & Chiba (2012) considered axisymmetric dark matter halos to model some of the Milky Way dSph galaxies.

For isolated galaxies, numerical experiments of CDM have clear predictions for these quantities. The shapes of (isolated) dark matter halos in the mass range M, are generally found to be triaxial, with axis ratios depending on the mass of the object (Frenk et al., 1988; Dubinski & Carlberg, 1991; Warren et al., 1992; Cole & Lacey, 1996; Thomas et al., 1998; Jing & Suto, 2002; Bailin & Steinmetz, 2005; Kasun & Evrard, 2005; Hopkins, Bahcall & Bode, 2005; Allgood et al., 2006; Bett et al., 2007; Hayashi, Navarro & Springel, 2007; Kuhlen, Diemand & Madau, 2007; Stadel et al., 2009; Diemand & Moore, 2011; Muñoz-Cuartas et al., 2011; Schneider, Frenk & Cole, 2012). In terms of their internal kinematics, the velocity ellipsoid is close to isotropic near the center of halos and becomes mildly radial towards the outskirts (Wojtak et al., 2005; Hansen & Moore, 2006; Ludlow et al., 2011).

For subhalos, however, less is known because of the demanding numerical resolution needed to model properly low mass halos orbiting within hosts of Milky Way mass. This situation has recently improved with simulations that are now able to successfully sample the mass function on these scales, such as the Via Lactea, CLUES, GHALO and Aquarius simulations (Diemand, Kuhlen & Madau, 2007; Diemand et al., 2008; Springel et al., 2008; Stadel et al., 2009; Libeskind et al., 2010). For instance, using the Via Lactea simulations Kuhlen, Diemand & Madau (2007) found that subhalos are also not spherical, although the effect of tides tends to make subhalos rounder than comparable objects in the field. These results prompt questions about the validity of some of the assumptions involved in the mass modeling of stellar kinematics in dwarfs. And although the orbital anisotropy of the stars in a dSph is likely unrelated to that of dark matter (and associated to the formation history), it is nonetheless valuable to explore the dynamical structure of subhalos because they provide the main potential.

A detailed study of the shape of the Milky Way mass Aquarius halos was presented in Vera-Ciro et al. (2011). Here we extend this analysis to lower mass objects, both subhalos of the main central halo and field halos, up to from the center of the main Aquarius halo. The paper is organized as follows. In Section 2 we describe the numerical simulations, introduce the methods we use to determine halo shapes and explore the convergence of the results. In Section 3 we compare the properties of subhalos and isolated objects of similar mass. We analyze subhalo shapes in the context of the kinematic modeling of dwarfs around the Milky Way in Section 4 and summarize our main results in Section 5.

## 2 Shape measurements and convergence

We use the suite of cosmological -body simulations from the Aquarius project (Springel et al., 2008). These consist of six CDM halos (Aq-A to Aq-F), re-simulated with five different levels of resolution within the cosmology , , , . The simulations use the zoom-in technique, with a high-resolution region that extends at up to Mpc from the center of each main halo. This exceeds the typical virial radius of the Aquarius halos by 5–10 times and allows us to identify isolated halos that have been unaffected by tidal forces (see Springel et al., 2008, for further details).

In most of the analysis that follows we focus on the level-2 resolution runs, with a mass per particle . However, we use the other Aquarius levels to test the convergence of our results. Halos and subhalos are identified using the subfind algorithm (Springel et al., 2001). We keep all structures identified with at least 20 particles. We will call the central subhalo of a group a field/isolated halo. In this work we consider field halos up to a distance of Mpc from the center of the main Aquarius halo, to avoid contamination of low-resolution particles.

To measure the shape of halos in the simulations we follow the same approach as Vera-Ciro et al. (2011) and iteratively compute the inertia tensor in ellipsoidal regions. At a given radius, the algorithm begins with a spherical contour which is reshaped and reoriented according to the principal axis of the normalized inertia tensor for the encompassed material, until convergence is reached (Allgood et al., 2006). More specifically, we define the normalized inertia tensor as

(1) |

where is a distance measure to the -th particle and is the set of particles of interest. Representing dark matter halos as ellipsoids of axis lengths , the axis ratios and are the ratios of the square-roots of the eigenvalues of , and the directions of the principal axes are given by the corresponding eigenvectors. Initially, the set is given by all particles located inside a sphere which is re-shaped iteratively using the eigenvalues of . The distance measure used is , where and are updated in each iteration. In practice we find that the algorithm converges (i.e. the variation in the shape between successive iterations is %) when there are as few as 200 particles in set . Notice that this is a more stringent criteria than required by plain identification of bound objects in subfind, which is here taken to be only 20 particles.

In Vera-Ciro et al. (2011) we showed that shapes can be robustly measured from the convergence radius, outwards (Power et al., 2003; Navarro et al., 2010). is defined so that the ratio between the local relaxation time and the dynamical time at the virial radius equals where:

(2) |

where is the number of particles inside the radius , is the spherically averaged density and the critical density. We adopt the value because this guarantees that the circular velocity profiles of the main halos are accurate to better than (Navarro et al., 2010). Note that this equation has to be numerically solved for with fixed for each object in the simulation. As a rule-of-thumb, we find that requiring a minimum of particles enclosed within the radius of interest (i.e. or ) ensures that more than of the subhalos satisfy this constraint.

Fig. 1 shows that the same criteria applied to our sample of subhalos also ensure convergence of halo shape estimates. We compare the results for the Aq-A run in all resolution levels 1–5 (red to black, respectively). The left panels show, as a function of halo maximum circular velocity , the mean axis ratios computed at , here defined as the ellipsoidal contour enclosing the most bound particles identified by subfind. The thin lines correspond to the entire sample of subhalos, whereas the thick curve shows only “converged” objects (those where ). At level 2, the one used for most of our analysis, the mean axis ratios agree with the highest resolution run Aq-A-1 to better than across the full spectrum of “converged” subhalos.

A similar conclusion is reached for the inner regions of subhalos, as shown by the right panels of Fig. 1. Here, and are computed at the radius of the peak circular velocity , which is typically approximately nine times smaller than for field halos and approximately six times smaller than for subhalos. The number of objects for which is roughly ten times lower than those with . This explains the relatively more noisy behavior of the curves on the right column compared to those on the left (especially for the lowest two resolution runs, where typically less than ten objects satisfy the convergence condition). In general, a subhalo whose shape has converged at the radius has not necessarily converged at the radius.

Besides the mean trends shown in Fig. 1, we have also explored the convergence of halo shapes on an object-by-object basis. In order to do this, we identify the same subhalos in several resolution levels of the Aq-A halo by matching the Lagrangian positions of all the particles assigned to a substructure by subfind back in the initial conditions (see Section 4.2 Springel et al., 2008, for further details). In addition to this criterion, we impose a maximum deviation on the orbital path of matched objects in different level runs. This is to ensure that the evolution of each subhalo has been comparable in the different resolution runs also in the non-linear regime. More specifically, we define

(3) |

with is the position^{1}^{1}1The positions of all
halos and subhalos are defined by the particle with the minimum potential
energy. of the subhalo with respect to the main halo center at the -th
resolution level. This is computed for every snapshot from the time the object
is first identified in the simulation box until present day. We consider only
structures for which .

A total of 260 substructures are successfully matched in all levels 1,2 and 3 of the Aq-A halo by this procedure. For each object, we define , where measured at and the lower indices indicate the resolution level (1 and 3 in the example above). By construction, for well-converged objects. We show the distribution of in Fig. 2. The light gray histogram corresponds to all matched objects, and is significantly broader than the distribution for the converged sample (i.e. the subset of the matched sample formed only with the subhalos that have ), which is shown in dark gray.

We illustrate this more clearly on the right panels in Fig. 2, which show the behavior of as a function of radius for three subhalos in the sample. Small colored dots indicate their value on the histogram on the left. The various curves correspond to the results for different resolution levels as indicated by the labels. For each subhalo the convergence radius is marked with a vertical arrow and the position of the contour by a vertical thick gray line. The top panel shows a typical example of an unconverged object: the peak of the circular velocity occurs at a smaller radius than for levels 2 and 3. On the other hand, subhalos in the middle and bottom panels have and have therefore converged (according to our criterion) at all these levels.

Notice that a large number of particles do not guarantee convergence. For instance, the unconverged object on the top right panel of Fig. 2 has and particles in levels 3 and 2, respectively. These are significantly larger than the values previously used in the literature (e.g., Kuhlen, Diemand & Madau, 2007; Knebe et al., 2008a, b), and highlights the need to impose a second criterion to measure individual shapes reliably. With our criterion, for only % of the halos with particles have the shapes at the contour converged (i.e. ). The situation improves significantly for the contour, where 99.6% of such objects have converged.

The distribution of axis ratios for converged objects shown in the left panel of Fig. 2 has a standard deviation , meaning that 68.3% of the objects shapes determined at the Aquarius level-3 deviate less than 8% from their value in the highest resolution run. Since we focus on the level-2 runs for the analysis that follows, we expect resolution effects in our sample to be negligible.

The above discussion shows that our criterion for convergence is relatively strict. There are 21403 subhalos with at least 200 particles within the radius in all the Aquarius simulations, and we find that the inertia tensor algorithm converges for 11483 subhalos at the contour, and for 13970 at the contour. If we now impose our convergence criteria, there remain 412 and 6072 subhalos with well-determined shapes at the and contours, respectively. For halos in the field our convergence criteria leads to a reduction of 96% and 35% for the and contours respectively. As expected, there is a larger proportion of field objects whose shapes can be measured at the contour. However, despite this significant reduction in sample size, we have gained in the reliability of the shape determination for individual halos.

Therefore, in the next section we focus on those halos which satisfy our convergence criteria.

## 3 Halo shapes as a function of mass and environment

We proceed to characterize the variations in the axis ratios and of dark matter halos according to their mass or, equivalently, their maximum circular velocity. The left column of Fig. 3 shows and for isolated objects measured at the radius (labeled contour) and at ( contour) in the top and bottom panels, respectively. A thick line indicates the median trend of our sample and the open symbols at the high mass end show the results for the main Aquarius halos from Vera-Ciro et al. (2011). In agreement with previous work, we find that axis ratios tend to decrease gently with (Allgood et al., 2006; Macciò et al., 2007; Hahn et al., 2007; Bett et al., 2007; Muñoz-Cuartas et al., 2011), although we now explore a different mass regime.

Inspection of the top and bottom panels shows that the dependence of the axis ratios with circular velocity is somewhat steeper when measured at than at the contours. Typically, our lowest mass objects have inner axes that are rounder by than those of Milky Way-like halos. Nevertheless, the scatter from object to object at fixed circular velocity also is larger at , as indicated by the shaded regions.

A comparison between the left and the right column of Fig. 3 reveals that there are only small differences between subhalos and isolated objects. To make this comparison easier we overplot in the panels on the right the linear fits obtained for field halos. This shows that, on average, subhalos are slightly more spherical than field halos at a given , but differences are well within the scatter in the samples. The number of converged objects in the case of subhalos is 385 and 1522 for and contours, respectively.

Could the differences between field halos and subhalos be caused by measuring shapes at different physical radii? It has been shown in the literature that tidal evolution can significantly decrease in satellites while affecting significantly less (Hayashi et al., 2003; Kravtsov, Gnedin & Klypin, 2004; Peñarrubia, Navarro & McConnachie, 2008). In that case, the measurement of the halo shape at would be at a smaller radius for a subhalo than for a halo in the field with the same , and the same holds for the contour. We address this in Fig. 4, where we show the minor-to-major axis profiles for individual field halos (left) and subhalos (right) of similar mass ( km s). The solid dots show the location of for individual objects; they indicate that the radii of the peak circular velocity are comparable in both samples and therefore show that this can not be reason for the different trends reported in Fig. 3. We thus confirm that, on average, subhalos of a given are slightly more spherical than comparable field halos at all radii, particularly in the outskirts. Kolmogorov–Smirnov tests indicate that the difference between both samples is statistically significant only at the contours (the Kolmogorov–Smirnov probability is 0.09 in that case versus 0.42 at ). However, the differences are well within the object-to-object scatter (see the bottom panel Fig. 4).

The similarity between the subhalo and field populations apparent in Figs 3 and 4 explains the lack of appreciable trends as a function of distance to the center of the main Aquarius halos, shown in Fig. 5. The typical axis ratios measured at the as well as contours do not depend on the distance to the host center up to distances . We have explicitly checked that this is not due to the averaging over random orbital phases of individual subhalos; using the apocenter distances instead of the instantaneous position leads to very similar results. On the other hand, the tidal field of the host seems to imprint a significant radial alignment on the subhalos (see the bottom panel of Fig. 5), which tend to orient their major axis radially to the center of the host, albeit with a large scatter. The signal is stronger close to the center of the main host halos and decreases steadily until it disappears at , where the distribution is consistent with random (Pereira, Bryan & Gill, 2008; Pereira & Bryan, 2009). Interestingly, inner (at ) and contours are well aligned within each object, as shown in the bottom right panel.

Notice that, although our findings suggest only small differences between the shapes of subhalos and field halos, the evolution of single objects under the effects of tidal disruption can differ significantly from the behavior of the population as a whole (Barber et al., in preparation.) In general, the analysis of a population of subhalos such as that shown in Figs 3–5 will be dominated, in number, by objects with recent infall times (and therefore not largely exposed to tidal effects), minimizing the differences between subhalos and field halos in good agreement with our results.

## 4 Application to the Modeling of Local Group Satellites

As discussed in the Introduction, the Local Group satellite galaxies are expected to inhabit dark matter subhalos comparable to those studied in the previous section. Since the contribution of the baryons to the gravitational potential of these systems is thought to be sub-dominant, the shape, dynamics and orbital structure of their host dark halos may be compared in a reasonably direct way to those of a suitable subset of the subhalos in the Aquarius simulations.

To select subhalos likely to host luminous satellites comparable to local dwarfs we use the semi-analytical model of Starkenburg et al. (2013). The semi-analytic model includes physical prescriptions for the treatment of relevant processes such as radiative cooling, chemical enrichment, star formation, supernova feedback, etc. The parameters in the model are tuned to simultaneously reproduce the luminosity function and spatial clustering of bright galaxies as well as the properties of satellites in the Local Group (De Lucia & Blaizot, 2007; De Lucia & Helmi, 2008; Li, De Lucia & Helmi, 2010).

Parameter | Median | ||
---|---|---|---|

0.27 | 0.07 | 0.10 | |

1.21 | 0.42 | 0.58 | |

9.06 | 1.87 | 1.72 | |

0.70 | 0.07 | 0.10 | |

0.83 | 0.10 | 0.09 | |

1.72 | 1.26 | 2.86 | |

1.60 | 0.13 | 0.09 | |

-2.42 | 0.15 | 0.20 |

### 4.1 The shapes of subhalos hosting luminous satellites

Fig. 6 shows the axis ratios as a function of distance along the major-axis for our sample of subhalos. This consists of subhalos within the virial radius of their hosts at and that resemble the classical satellites of the Milky Way in their luminosity, i.e. their -band absolute magnitudes are in the range . Each curve is plotted from the convergence radius out to the contour, and the color scale indicates the luminosity assigned by the semi-analytic model to the satellites.

Fig. 6 shows that the dwarf galaxies in the model are surrounded by subhalos that are triaxial, with axis ratios and typically increasing from the inner regions to the radius. The scatter from object to object is large, but the overall trend with radius is similar for all subhalos. The median profile, and -equivalent percentiles of the sample are given, respectively, by the black solid line and the gray shaded area. These dark matter subhalos have on average and at a radius of kpc, and turn more spherical close to the radius, where and . Individual inner shapes of halos/subhalos can be clearly seen in Fig. 7, where we show a scatter plot of and ratios measured at kpc. Only converged objects have been included. Different symbols are used for different samples: blue circles for field halos, red diamonds for subhalos and black circles for luminous subhalos. The fact that subhalos and field halos are well mixed in this plane confirms the lack of any significant trend between shape and distance to the main halo, in agreement with Fig. 5.

We may use Poisson’s equation to derive a relation between the and of the density (which our method measures), and those of the potential and . Following Vogelsberger et al. (2008), we introduce a generalized radius:

(4) |

where is the Euclidean distance, is the ellipsoidal radius and a characteristic scale. With this definition, for and for . Assuming that the potential at any point is

(5) |

where is the spherically symmetric potential associated with the Einasto profile (Einasto, 1965), we find for our sample a median and . The median and error of the parameter fits for the density and axis ratios obtained in this way are given in Table 1.

### 4.2 The internal kinematics of subhalos hosting luminous satellites

#### 4.2.1 Behavior along the major axis

Fig. 8 shows the radial velocity dispersion (top) and the orbital anisotropy (bottom), both as a function of distance along the major axis. These quantities are computed in ellipsoidal coordinates that follow the axis ratios of the mass density at each radius. In practice, we calculate the component of the velocity in the direction tangential to a given ellipsoid, , and the radial component is derived by subtraction in quadrature of from the total velocity dispersion. Therefore this corresponds very closely to the velocity dispersion along the spherical radial direction as it is computed along the ellipsoid’s major axis.

As in the previous figure, individual objects are shown with thin lines color-coded according to their -band absolute magnitude. The median trend of the sample is indicated by a black solid curve, -equivalent percentiles in gray shading. The spread in seen in the top panel of Fig. 8 is due to the difference in mass of the subhalos, which span a range .

A declining velocity dispersion profile as seen in the top panel of Fig. 8 is generally expected for halos in CDM because of the relation between the density and pseudo phase-space density, namely (Taylor & Navarro, 2001):

(6) |

with a normalization constant, , and , where , and are the characteristic
density, velocity and radius for the Einasto density profile,
respectively^{2}^{2}2See Appendix A, where we show that a
power-law fit is a reasonable description of the pseudo-phase-space of
subhalos just like it is for field halos.. When we fit the -profiles
individually for each subhalo, we obtain median values of and
, the latter indicating a slightly shallower fall off than for
field halos.

In the bottom panel of Fig. 8 we plot the ellipsoidal velocity anisotropy profiles. Here , where we calculate the component of the velocity tangential and radial to a given ellipsoid as explained above. The velocity anisotropy profiles of dark matter subhalos tend to decline with radius. In the inner regions the motions are slightly radially-biased ( at kpc), while the ellipsoid becomes increasingly tangential () at larger radii. This behavior is markedly different from the radially-biased ellipsoids found in isolated CDM halos, particularly in the outskirts (Cole & Lacey, 1996; Taylor & Navarro, 2001; Wojtak et al., 2005; Ludlow et al., 2010). This difference is a result of tidal forces, which preferentially remove particles with large apocenters on radial orbits. Fig. 8 also shows that subhalos rarely have a constant profile.

We may derive an expression for using the spherical Jeans Equation which relates the density, anisotropy and radial velocity dispersion of a system (Binney & Tremaine, 2008). In the case of an Einasto profile with logarithmic slope we find:

(7) |

where as before and is given by Eq. (11) (see Appendix A for a more detailed derivation of this expression). The velocity anisotropy is therefore dependent on the logarithmic slopes of the mass density and of the pseudo phase-space density, and , respectively.

This expression provides a reasonable fit out to . Whereas the limiting behavior of the anisotropy in the inner regions is similar for all subhalos, beyond a radius of kpc large variations are seen from object-to-object. These variations are still accounted for by Eq. (7) when each profile is fitted individually. We find that the exact shape of the anisotropy profile depends most strongly on , while determines where the velocity ellipsoid becomes tangential at large radii. On the other hand, variations in have a very minor effect.

#### 4.2.2 Axisymmetric description of the internal kinematics

Fig. 6 shows that subhalos are not significantly triaxial having in average , therefore, an axisymmetric approach may be sufficient to describe their internal kinematics. Moreover, an oblate approximation seems favoured by the overall distribution of inner shapes shown in Fig. 7, albeit with a large scatter.

To this end, in Fig. 9 we explore the velocity structure in the cylindrical radial and vertical directions. Each object is rotated such that the minor axis coincides with the -direction. The velocity dispersion components along the vertical () and radial () axes provide information on how dynamically hot a system is in both directions. We compute this along the two preferential axis, minor (left) and major (right), using at each radius a sphere that contains 400 particles. The velocity dispersion and are then computed within these volumes and displayed as a function of distance along the axes. As before, individual halos are shown in thin blue lines, the median values with solid thick lines.

The velocity dispersion profiles show a steady decline with radius, and as expected with a trend similar to the radial velocity dispersion shown in Fig. 6. As in that case, there is scatter from object-to-object mostly due to the range in subhalo mass, with the most luminous objects (light green in Fig. 9) having larger velocity dispersion than the fainter ones (red). The anisotropy is defined as and is shown in the bottom panels of Fig. 9. Interestingly, along the minor axis we find , in agreement with the assumptions by Hayashi & Chiba (2012) when modeling stars in the dSph. However, this is not true for the velocity dispersions along the major axis as shown in the right panel, where we find positive values of close to the center and at the radius. This trend, although systematic, is quite weak, with large scatter amongst the individual objects. It is important to realize that our measurements are subject to significant noise as we are now restricted to consider a small region around a given radius (i.e. we do not average over entire ellipsoidal shells as before and so each sphere contains fewer particles). The typical uncertainty for the individual curves is shown by the black vertical error bar, computed as the dispersion obtained from drawing 100 samples with replacement in each bin.

The fact that the velocity dispersions and are not equal implies that the subhalos’ distribution functions are a function of a third integral. Although we have shown this to be the case for dark matter satellites, it could also be true for the stars embedded in these systems. Therefore dynamical models of dSph may need to take this into account (Battaglia, Helmi & Breddels, 2013), since neglecting this fact can lead to unrealistic estimates of the shapes of the host dark matter halos (see Hayashi & Chiba, 2012).

The shapes of the subhalos are consistent with dynamical support by the velocity ellipsoid, as shown by Fig. 10. The vertical axis shows the local anisotropy , where is again the direction of the major axis and points along the minor axis. These quantities are calculated in a sphere with 400 particles located at kpc and the individual points in this figure correspond to the different subhalos. In the axisymmetric case, the virial theorem in tensor form (Binney, 2005) gives:

(8) |

where are the components of the potential-energy tensor (Binney & Tremaine, 2008). For ellipsoidal systems is a function of the ellipticity and is independent of the radial density profile (Roberts, 1962). is the streaming velocity along the -axis and is the velocity dispersion in the direction. The solid line in Fig. 9 indicates the prediction from Eq. (8) in the case of a dispersion supported system with . This prediction provides a reasonable representation of the simulated objects that agrees well with the very little rotation that we find: halos and subhalos show an average and more than of the sample has in their inner regions. The scatter, however, is large and cannot be explained solely on the basis of rotation of subhalos at small radii. Further factors such as departure from axisymmetry or the lack of dynamical equilibrium generated by tides may also contribute to the scatter seen in Fig. 10.

## 5 Conclusions

We have used the Aquarius simulations to study the shapes of field and satellite dark matter halos with emphasis on the mass range expected for the hosts of the dwarf galaxies in the Local Group. We have used an iterative method based on the normalized inertia tensor to characterize the principal axis lengths of halos and subhalos as a function of radius. In particular, we have explored in detail halo shapes measured in the inner regions (radius of maximum circular velocity ) and in the outskirts or contour. Although stars are more centrally concentrated than the dark matter, our resolution allows us to characterize halo shapes at radii as small as , starting to probe the regime traced by the outer stars in dwarf galaxies.

Through a comparison of objects in common between the different resolution levels of the Aquarius simulations, we have noticed that simple number of particles cuts do not guarantee convergence in the measured halo shapes, especially in the inner regions. We find that instead the convergence radius (defined as the threshold in the ratio between the local relaxation time and the dynamical time at the virial radius (Power et al., 2003)) provides a good estimate of the radius where the axis ratios are robustly determined (with an error %).

We have found that the typical axis ratios of isolated halos in the Aquarius simulations decrease with increasing mass, or equivalently maximum halo circular velocity , i.e. low mass objects tend to be more spherical than Milky Way-like objects. These trends are well approximated by a relation between the axis ratio measured at the radius and , i.e. , while this relation is slightly steeper if the axis ratio is measured at contour in which case . The differences in the shapes of field versus satellite halos are small and within the intrinsic scatter of the samples. Nonetheless, at a fixed , subhalos tend to have larger axis ratios than isolated objects in the field.

The similarity between subhalos and field objects is also apparent in the lack of significant trends in the axis ratios with distance to the main host halo, . We find, however, that the alignment of the ellipsoids varies with : dark matter halos at close distances from the host center tend to be oriented preferentially with their major axis pointing radially. The signal disappears only for , where the orientations are consistent with random.

We have also focused on the properties of subhalos likely to host analogs of the classical satellites of the Milky Way (), according to the semianalytic model of galaxy formation run on the Aquarius suite by Starkenburg et al. (2013). Our analysis indicates that these galaxies are hosted by mildly triaxial dark matter objects with minor-to-major axis ratios and intermediate-to-major in the first kiloparsec with a clear trend towards becoming axisymmetric in the outskirts. Their internal orbital structure shows evidence of being affected by tidal forces from their hosts (i.e. the main Aquarius halos), since the velocity anisotropy becomes tangential with radius, in clear contrast to what is found for isolated systems. We have also found that this orbital structure may be modeled in the axisymmetric context, where the velocity anisotropy along the minor axis, and declines with distance along the major axis. These results may be used to motivate more realistic models of the subhalos hosting satellite galaxies like those observed around the Milky Way.

## Acknowledgements

We thank Else Starkenburg for access to the semi-analytical catalog of the Aquarius simulations. AH gratefully acknowledges financial support from the European Research Council under ERC-Starting Grant GALACTICA-240271. The authors thank the hospitality of the Kavli Institute for Theoretical Physics, Santa Barbara during the program “First Galaxies and Faint Dwarfs: Clues to the Small Scale Structure of Cold Dark Matter”, where part of this work was completed under the support of the National Science Foundation Grant No. NSF PHY11-25915. This work was supported in part by the National Science Foundation under Grant No. PHYS-1066293 and the hospitality of the Aspen Center for Physics.

## References

- Abramowitz & Stegun (1972) Abramowitz M., Stegun I. A., 1972, Handbook of Mathematical Functions. Dover Publications, New York
- Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
- Amorisco & Evans (2011) Amorisco N. C., Evans N. W., 2011, MNRAS, 411, 2118
- Amorisco & Evans (2012) Amorisco N. C., Evans N. W., 2012, MNRAS, 419, 184
- Bailin & Steinmetz (2005) Bailin J., Steinmetz M., 2005, ApJ, 627, 647
- Battaglia, Helmi & Breddels (2013) Battaglia G., Helmi A., Breddels M., 2013, New Astron. Rev., 57, 52
- Battaglia et al. (2008) Battaglia G., Helmi A., Tolstoy E., Irwin M., Hill V., Jablonka P., 2008, ApJ, 681, L13
- Bett et al. (2007) Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215
- Binney (2005) Binney J., 2005, MNRAS, 363, 937
- Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press, Princeton
- Blumenthal et al. (1984) Blumenthal G. R., Faber S. M., Primack J. R., Rees M. J., 1984, Nat., 311, 517
- Boylan-Kolchin, Bullock & Kaplinghat (2012) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2012, MNRAS, 422, 1203
- Breddels & Helmi (2013) Breddels M. A., Helmi A., 2013, A&A, 558, A35
- Breddels et al. (2013) Breddels M. A., Helmi A., van den Bosch R. C. E., van de Ven G., Battaglia G., 2013, MNRAS, 433, 3173
- Cole & Lacey (1996) Cole S., Lacey C., 1996, MNRAS, 281, 716
- De Lucia & Blaizot (2007) De Lucia G., Blaizot J., 2007, MNRAS, 375, 2
- De Lucia & Helmi (2008) De Lucia G., Helmi A., 2008, MNRAS, 391, 14
- Dehnen & McLaughlin (2005) Dehnen W., McLaughlin D. E., 2005, MNRAS, 363, 1057
- Diemand, Kuhlen & Madau (2007) Diemand J., Kuhlen M., Madau P., 2007, ApJ, 657, 262
- Diemand et al. (2008) Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J., 2008, Nat., 454, 735
- Diemand & Moore (2011) Diemand J., Moore B., 2011, Adv. Sci. Lett., 4, 297
- Dubinski & Carlberg (1991) Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496
- Einasto (1965) Einasto J., 1965, Trudy Astrofiz. Inst. Alma-Ata, 5, 87
- Frenk et al. (1988) Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507
- Ghigna et al. (1998) Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 1998, MNRAS, 300, 146
- Gilmore et al. (2007) Gilmore G., Wilkinson M. I., Wyse R. F. G., Kleyna J. T., Koch A., Evans N. W., Grebel E. K., 2007, ApJ, 663, 948
- Gott & Rees (1975) Gott, III J. R., Rees M. J., 1975, A&A, 45, 365
- Hahn et al. (2007) Hahn O., Porciani C., Carollo C. M., Dekel A., 2007, MNRAS, 375, 489
- Hansen & Moore (2006) Hansen S. H., Moore B., 2006, New Astron., 11, 333
- Hayashi, Navarro & Springel (2007) Hayashi E., Navarro J. F., Springel V., 2007, MNRAS, 377, 50
- Hayashi et al. (2003) Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn T., 2003, ApJ, 584, 541
- Hayashi & Chiba (2012) Hayashi K., Chiba M., 2012, ApJ, 755, 145
- Hopkins, Bahcall & Bode (2005) Hopkins P. F., Bahcall N. A., Bode P., 2005, ApJ, 618, 1
- Jardel & Gebhardt (2012) Jardel J. R., Gebhardt K., 2012, ApJ, 746, 89
- Jing & Suto (2002) Jing Y. P., Suto Y., 2002, ApJ, 574, 538
- Kamionkowski, Koushiappas & Kuhlen (2010) Kamionkowski M., Koushiappas S. M., Kuhlen M., 2010, Phys. Rev. D, 81, 043532
- Kasun & Evrard (2005) Kasun S. F., Evrard A. E., 2005, ApJ, 629, 781
- Kleyna et al. (2001) Kleyna J. T., Wilkinson M. I., Evans N. W., Gilmore G., 2001, ApJL, 563, L115
- Klypin et al. (1999a) Klypin A., Gottlöber S., Kravtsov A. V., Khokhlov A. M., 1999a, ApJ, 516, 530
- Klypin et al. (1999b) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999b, ApJ, 522, 82
- Knebe et al. (2008a) Knebe A., Draganova N., Power C., Yepes G., Hoffman Y., Gottlöber S., Gibson B. K., 2008a, MNRAS, 386, L52
- Knebe et al. (2008b) Knebe A., Yahagi H., Kase H., Lewis G., Gibson B. K., 2008b, MNRAS, 388, L34
- Kravtsov, Gnedin & Klypin (2004) Kravtsov A. V., Gnedin O. Y., Klypin A. A., 2004, ApJ, 609, 482
- Kuhlen, Diemand & Madau (2007) Kuhlen M., Diemand J., Madau P., 2007, ApJ, 671, 1135
- Li, De Lucia & Helmi (2010) Li Y.-S., De Lucia G., Helmi A., 2010, MNRAS, 401, 2036
- Libeskind et al. (2010) Libeskind N. I., Yepes G., Knebe A., Gottlöber S., Hoffman Y., Knollmann S. R., 2010, MNRAS, 401, 1889
- Łokas (2002) Łokas E. L., 2002, MNRAS, 333, 697
- Łokas (2009) Łokas E. L., 2009, MNRAS, 394, L102
- Łokas, Mamon & Prada (2005) Łokas E. L., Mamon G. A., Prada F., 2005, MNRAS, 363, 918
- Long & Mao (2010) Long R. J., Mao S., 2010, MNRAS, 405, 301
- Ludlow et al. (2010) Ludlow A. D., Navarro J. F., Springel V., Vogelsberger M., Wang J., White S. D. M., Jenkins A., Frenk C. S., 2010, MNRAS, 406, 137
- Ludlow et al. (2011) Ludlow A. D., Navarro J. F., White S. D. M., Boylan-Kolchin M., Springel V., Jenkins A., Frenk C. S., 2011, MNRAS, 415, 3895
- Macciò et al. (2007) Macciò A. V., Dutton A. A., van den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55
- Mateo (1998) Mateo M. L., 1998, ARA&A, 36, 435
- Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJL, 524, L19
- Muñoz-Cuartas et al. (2011) Muñoz-Cuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584
- Navarro, Frenk & White (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
- Navarro, Frenk & White (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
- Navarro et al. (2010) Navarro J. F. et al., 2010, MNRAS, 402, 21
- Peñarrubia, Navarro & McConnachie (2008) Peñarrubia J., Navarro J. F., McConnachie A. W., 2008, ApJ, 673, 226
- Pereira & Bryan (2009) Pereira M., Bryan G. L., 2009, in Bulletin of the American Astronomical Society, Vol. 41, American Astronomical Society Meeting Abstracts #213, p. 315.03
- Pereira, Bryan & Gill (2008) Pereira M. J., Bryan G. L., Gill S. P. D., 2008, ApJ, 672, 825
- Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
- Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
- Richstone & Tremaine (1986) Richstone D. O., Tremaine S., 1986, AJ, 92, 72
- Roberts (1962) Roberts P. H., 1962, ApJ, 136, 1108
- Schneider, Frenk & Cole (2012) Schneider M. D., Frenk C. S., Cole S., 2012, J. Cosmol. Astropart. Phys., 5, 30
- Springel et al. (2008) Springel V. et al., 2008, MNRAS, 391, 1685
- Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
- Stadel et al. (2009) Stadel J., Potter D., Moore B., Diemand J., Madau P., Zemp M., Kuhlen M., Quilis V., 2009, MNRAS, 398, L21
- Starkenburg et al. (2013) Starkenburg E. et al., 2013, MNRAS, 429, 725
- Stoehr et al. (2002) Stoehr F., White S. D. M., Tormen G., Springel V., 2002, MNRAS, 335, L84
- Strigari et al. (2007) Strigari L. E., Bullock J. S., Kaplinghat M., Diemand J., Kuhlen M., Madau P., 2007, ApJ, 669, 676
- Strigari et al. (2008) Strigari L. E., Bullock J. S., Kaplinghat M., Simon J. D., Geha M., Willman B., Walker M. G., 2008, Nat., 454, 1096
- Strigari, Frenk & White (2010) Strigari L. E., Frenk C. S., White S. D. M., 2010, MNRAS, 408, 2364
- Taylor & Navarro (2001) Taylor J. E., Navarro J. F., 2001, ApJ, 563, 483
- Thomas et al. (1998) Thomas P. A. et al., 1998, MNRAS, 296, 1061
- Tormen (1997) Tormen G., 1997, MNRAS, 290, 411
- Tormen, Diaferio & Syer (1998) Tormen G., Diaferio A., Syer D., 1998, MNRAS, 299, 728
- Vera-Ciro et al. (2013) Vera-Ciro C. A., Helmi A., Starkenburg E., Breddels M. A., 2013, MNRAS, 428, 1696
- Vera-Ciro et al. (2011) Vera-Ciro C. A., Sales L. V., Helmi A., Frenk C. S., Navarro J. F., Springel V., Vogelsberger M., White S. D. M., 2011, MNRAS, 416, 1377
- Vogelsberger et al. (2008) Vogelsberger M., White S. D. M., Helmi A., Springel V., 2008, MNRAS, 385, 236
- Walker (2013) Walker M., 2013, Dark Matter in the Galactic Dwarf Spheroidal Satellites, Springer Netherlands, p. 1039
- Walker et al. (2009) Walker M. G., Mateo M., Olszewski E. W., Peñarrubia J., Wyn Evans N., Gilmore G., 2009, ApJ, 704, 1274
- Walker & Peñarrubia (2011) Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20
- Warren et al. (1992) Warren M. S., Quinn P. J., Salmon J. K., Zurek W. H., 1992, ApJ, 399, 405
- White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
- Wilkinson et al. (2002) Wilkinson M. I., Kleyna J., Evans N. W., Gilmore G., 2002, MNRAS, 330, 778
- Wojtak et al. (2005) Wojtak R., Łokas E. L., Gottlöber S., Mamon G. A., 2005, MNRAS, 361, L1
- Wolf et al. (2010) Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., Muñoz R. R., Simon J. D., Avedo F. F., 2010, MNRAS, 406, 1220

## Appendix A Einasto Profiles and The spherical Jeans equation

Following the method of Vera-Ciro et al. (2013) we fit a Einasto profile to the circular velocity profile of each individual subhalo of our sample. For each object we take 20 bins equally spaced in logarithmic space between . We compute the cumulative circular velocity profile and fit an Einasto model by minimizing the merit function against the free parameters , and . Here, is the circular velocity corresponding to an Einasto profile with a scale radius (the radius at which the density profile has a slope ), a characteristic density at equal to and a shape parameter that controls the overall slope of the profile, . We deliberately chose the circular velocity profile over the more widely used density profile which is more sensitive to shot noise in each bin (see Vera-Ciro et al., 2013, for more details).

The density profile can be written as,

(9) |

where

(10) |

is therefore a dimensionless function of the dimensionless variable . In this spirit, it is possible to define a set of scaling factors in which we can express the dynamics of the system, namely, a characteristic mass and characteristic velocity [note that is not ]. The enclosed mass within a radius is therefore , where

(11) | |||||

with the lower incomplete gamma function. In a similar fashion we can define a dimensionless version of the radial velocity dispersion .

It has been previously reported in the literature that the pseudo phase-space density profile of isolated dark matter halos can be well modeled by a single power law , (Taylor & Navarro, 2001; Dehnen & McLaughlin, 2005; Navarro et al., 2010; Ludlow et al., 2011). Fig. 11 shows measured for our sample of subhalos hosting luminous satellites. We have found that the slope is slightly shallower than for objects in the field ( versus for isolated halos). The best-fitting values of the parameters, and their variance, are given in Table 1.

Notice that the pseudo-phase-space profiles start to deviate from a power-law in the outer regions, likely induced by ongoing tidal stripping. This typically occurs for , which is roughly the same scale at which our fit for the radial velocity dispersion deviates from the mean subhalo trends shown in Fig. 8.

Using the dimensionless quantities introduced before, we can now write the spherical Jeans equation as:

(12) | |||||

We can use Eq. (10) to further reduce this expression:

(13) |

The limiting value of this expression at small radii can be obtained from a Taylor expansion around zero. For we may use that (Abramowitz & Stegun, 1972). Finally

which leads to

(14) |