The evolution of dwarf galaxy satellites with different dark matter density profiles in the ErisMod simulations. I. The early infalls
We present the first simulations of tidal stirring of dwarf galaxies in the Local Group carried out in a fully cosmological context. We use the ErisDARK cosmological simulation of a MW-sized galaxy to identify some of the most massive subhalos () that fall into the main host before . Subhalos are replaced before infall with high-resolution models of dwarf galaxies comprising a faint stellar disk embedded in a dark matter halo. The set of models contains cuspy halos as well as halos with ”cored” profiles (with asymptotic inner slope ) consistent with recent results of hydrodynamical simulations of dwarf galaxy formation. The simulations are then run to with as many as 54 million particles and resolution as small as pc using the new parallel N-Body code ChaNGa. The stellar components of all satellites are significantly affected by tidal stirring, losing stellar mass and undergoing a morphological transformation towards a pressure supported spheroidal system. However, while some remnants with cuspy halos maintain significant rotational flattening and disk-like features, all the shallow halo models achieve and round shapes typical of dSph satellites of the MW and M31. Mass loss is also enhanced in the latter, and remnants can reach luminosities and velocity dispersions as low as those of Ultra Faint Dwarfs (UFDs) In two cases a merger transforms the objects quickly into a spheroidal system already before infall. We argue that cuspy progenitors must be the exception rather than the rule among satellites of the MW since all the MW and M31 satellites in the luminosity range of our remnants are dSphs, a result matched only in the simulation with ”cored” models. This suggests an early and efficient transformation of cusps into cores for Local Group satellites.
Dwarf galaxy satellites of the Milky Way (MW) and M31 are the main focus of ”near-field cosmology”, providing a plethora of useful constraints and tests for galaxy formation theories as well as for the nature of dark matter in the CDM paradigm. Yet their origin is still not completely understood. Nearly all the known satellites, located within the expected virial radius of the bright spirals, i.e at kpc, are classified as dwarf spheroidal galaxies (dSphs). They have fairly round stellar components with low or negligible rotation but relatively high velocity dispersions (Mateo (1998); McConnachie (2012)). At the same time, nearly all of the dwarf galaxies located outside such radius and far from M31 are gas-rich dwarf irregular galaxies (dIrrs), which have ongoing or recent star formation, flattened disky shapes, and exhibit substantial rotational support (). The different spatial distribution of dwarf galaxy types is known as the morphology-density relation, and reflects an analogous effect seen in galaxy clusters and rich groups on larger mass scales (Dressler (1980)).
Environmental mechanisms have been often invoked to transform dIrrs into dSphs (Mateo (1998); McConnachie (2012)). One of the mechanisms proposed to achieve such transformation is tidal stirring, namely the cumulative effect of repeated tidal shocks occurring near the pericenter of the eccentric orbits of satellites (Mayer (2010)).
The process has been investigated thoroughly with collisionless N-Body simulations adopting high resolution models of dwarf galaxies interacting with rigid or symmetric models of hosts on orbits with high eccentricities as those of cosmological subhalos (Mayer et al. (2001a); Mayer et al. (2001b) ;Klimentowski et al. (2007); Kazantzidis et al. (2011a)). It has also been investigated in combination with ram pressure mass removal using hydrodynamical simulations of dwarf-main host interaction in which a diffuse gaseous halo consistent with observational constraints is added to the dark halo of the main host (Mayer et al. (2006); Mayer (2010)).
These studies established tidal stirring as a major mode of morphological evolution of dwarf galaxies, turning them from faint rotationally supported disks into even fainter pressure supported spheroidals, thus providing a natural explanation for the morphology-density relation ubiquitously observed in galaxy groups, including our own Local Group (Grebel et al. (1999)).
Direct tidal heating as well as indirect heating stemming from bar and buckling instabilities triggered by tidal shocks have been both shown to play a role in tidal stirring, although the latter requires initial disks with mass-to-light ratios at the upper end of those observed in dwarf galaxies (Mayer et al. (2001b); Kazantzidis et al. (2011a)). These works also determined that ram pressure mass removal aided by tidal shocks can rapidly remove the gas content of satellites explaining why dSphs are devoid of gas and why they can have very high mass-to-light ratios.
Recently the effect of the dark matter density profile of dwarf galaxies on tidal stirring has been studied, motivated by the new findings of cosmological hydrodynamical simulations of field dwarf galaxy formation which have shown that feedback from supernovae and massive stars can turn the cuspy halos predicted by CDM into halos with shallower inner density slopes ( in [0.3, 0.7] range), depending on the star formation rate achieved in the galaxy (Governato et al. (2010); Pontzen & Governato (2012); Shen et al. (2014); Governato et al. (2012); Di Cintio et al. (2013); Trujillo-Gomez et al. (2015)). There have also been attempts to show the existence of this effect specifically for dwarf galaxy satellites, albeit in this case resolution is a potential issue ((Brooks et al., 2013)). Łokas et al. (2012) and Kazantzidis et al. (2013) have thus updated the previous tidal stirring simulations by carrying out N-Body simulations in which disky dwarf galaxies have shallower dark matter profiles as those arising in the cosmological hydro simulations, and have checked results against varying resolution. They have found that the transformation by tidal stirring is enhanced significantly when a shallower density profile is adopted ( or ). They have also found that, since mass loss is also enhanced, remnants of dwarfs with shallow dark matter profiles can evolve into very faint objects such as Ultra Faint Dwarfs (Łokas et al. (2012)). The enhanced mass loss and transformation follow from the lower binding energy and more impulsive dynamical response of the system in the case of objects with a shallower halo profile.
However, these recent simulations, exactly like those of the last decade, are fairly idealized since they completely lack the cosmological context. Hence orbits are also idealized, being disconnected with the epoch of fall inside the host and with the shape of its gravitational potential, which is just assumed to be spherically symmetric. Instead, cosmological simulations show that the host is triaxial and grows substantially in mass and radius over the last ten billion years (Diemand et al. (2007)), which implies satellites can be accreted on increasingly wider orbits over time. Furthermore, in the idealized simulations subhalos can only interact with the main host, and not between them. Also, for simplicity these simulations have always explored a fixed or limited mass range for the satellites, rather than considering the wide range of masses from those of actual subhalos. It is expected that, the initial mass, and also the stellar mass relative to that of the halo, ought to have an important effect on the response of the system to tidal shocks. Furthermore, other studies have shown that mergers or strong interactions, while rare, can occur among dwarf galaxies just before infall, turning them into dwarf spheroidals very quickly (Kazantzidis et al. (2011b)), but leaving it completely open to identify the distinct signature of a remnant of one of such mergers as opposed to one from tidal stirring. Ideally, all these difficulties should be overcome in a fully cosmological hydro simulation of a MW sized galaxy, but with current resolutions around 100 pc and in the baryons numerical results cannot be put on firm grounds for objects that should be only a few hundreds parsecs in size.
Therefore, in order to make progress beyond the idealized simulation we devise a numerical strategy that in the past has only been attempted for simulations of galaxies in clusters ( Dubinski (1999); Mastropietro et al. (2005)). In this approach, known as the ”replacement technique”, cosmological subhalos are replaced with high resolution models of dwarf galaxies within the original cosmological volume. Models are constructed based on observational and theoretical constraints, and matched as closely as possible to the total virial mass of the parent subhalo. After the replacement the simulation is rerun to the present epoch. In this way galaxy satellites are automatically placed on realistic orbits and are subject to the full tidal interaction history of the parent subhalos. Due to the computational facilities and improved speed and parallelism of the codes, the simulations can span the equivalent of more than 10 Gyr of evolution with this technique, while in previous cluster simulations only a small time span, from onwards, could be covered (Mastropietro et al. (2005)). In this paper we present the first replacement simulation carried out for a MW sized halo, using the ErisDARK cosmological simulation carried out in the concordance CDM model (Kuhlen et al. (2013); Pillepich et al. (2014)). The simulations are carried out with state-of-the-art massively parallel N-Body codes, PKDGRAV2 (Stadel (2001)) and ChaNGa (Jetley et al. (2008); Jetley et al. (2010); Menon et al. (2015); http://www-hpcc.astro.washington.edu/tools/changa.html) on the CRAYXE6 ”Monte Rosa” of the Swiss National Supercomputing Centre.
The paper continues with Section 2 where the replacement technique and numerical simulations are presented in details. Next the analysis methods and results are enumerated in Section 3. Section 4 contains an extensive discussion regarding the implications of the results for the origin of dSphs, including their relation to observational data. Finally Section 5 includes the conclusion and outlook of the presented work.
2 Numerical simulations and replacement technique
2.1 Initial conditions and selection procedure
The initial conditions of the simulations of this paper are based on the ErisDark run (Pillepich et al. (2014)), which is a dark matter only cosmological gravity run of 40 million particles using the CDM cosmological model ( = 73 km , ,, and ). The cosmological volume contains a 1 Mpc scale zoom-in region with three levels of refinement that lead to a minimum mass of the dark matter particle of 1.2 and a minimum softening length of 124 pc (Pillepich et al. (2014)). In the zoom-in region an isolated MW-sized halo assembles, with mass consistent with the lower limit for the estimated halo mass of the Milky Way, which provides a good fit with current kinematical data on the Milky Way stellar halo (Rashkov et al. (2013)). The companion hydrodynamical simulations, Eris and ErisLE, have shown to produce realistic spiral galaxies by , which reproduce nearly all the main properties of the disk and bulge of the Milky Way (Guedes et al. (2011); Mayer (2012); Guedes et al. (2013); Bird et al. (2013)).
Satellite galaxies falling into the main host after the last major merger, which takes place at , are considered for replacement with hi-resolution dwarf galaxy models including stars. We analize the entire available set of outputs of the original simulation with the use of the Amiga Halo Finder (AHF) (Gill et al. (2004);Knollmann & Knebe (2009)) with a virial overdensity criterion that varies with redshift according to the cosmology (Bryan & Norman (1998)). In this paper, in particular, we analyze the earliest infalling population of satellites among which there are survivors at . In the next paper we will present larger populations of dwarf galaxies which are accreted at later epochs.
We thus proceed in the following way. First of all, among subhaloes
with bound masses greater than at z = 0, the satellite with the earliest infall is identified.
The adopted threshold mass ensures that we are considering well resolved objects, with at least particles
(presumably their progenitors are even better resolved at higher redshift, before tidal mass loss begins), to
minimize numerical effects on their internal dynamical evolution (Moore et al. (1996)).
The object enters the host’s virial radius just before , when it is still unrelaxed, being the result of
a merger between two subhalos. We trace the satellite back to
when the two progenitor merging subhalos are still well separated, namely their virial volumes do not overlap.
Therefore the epoch corresponding to z = is set as the epoch of replacement of the original ErisDark’s infalling subhalos with their
high resolution versions.
Next we identify other subhalos that at this epoch are located between 1 and 4 of the host. We restrict the sample to the most massive ones, having pre-infall masses around and above. Among these we reject 60% that cannot be properly identified and replaced because they are either too elongated (axis ratio less than 0.65, possibly tidally disturbed by companions) or have a substantial fraction of their mass in subhalos. Four objects remain after this rejection procedure, with masses in the range . Note that at this translates into satellites having masses in the range 80-300 times smaller than the host, which has a virial mass of and a virial radius of 55.6 kpc at the corresponding epoch (both relatively high since ErisDARK has a fast assembly history with no more major mergers after , Pillepich et al. (2014)). We remark that the lower end of the satellite-to-host mass ratio is still in a regime in which dynamical friction is effective in eroding the orbital energy over many Gyrs (Colpi et al. (1999)), while this is a negligible effect in published idealized simulations that assume satellites of similar mass but within a host having a fixed mass (Kazantzidis et al. (2013)). Finally, we observe that none of the four objects has a mass above at z = 0 in the original ErisDark run. Therefore they are not part of the largest subhaloes at the final epoch, which can well exceed this mass at low z (Paper II, Pillepich et al. (2014)). In total we thus replace 6 subhalos (figure 1), of which 2 are merging during infall and the other four are accreted individually. The halo properties are listed in table 1.
2.2 Model galaxy generation
The replacement models are generated using the GalactICS code (Kuijken & Dubinski (1995), Widrow & Dubinski (2005), Widrow et al. (2008)), which generates multi-component equilibrium models of galaxies. Each model galaxy is comprised of two components: the dark matter halo and the stellar disk. It is assumed by construction that the initial model has a stellar disk since this stems from the assumption that disky dwarfs are the progenitors of dSphs, the core idea of tidal stirring. The code creates halos with negligible angular momentum and a spherically symmetric density profile as described by the relation:
In the previous equation the scale radius and velocity scale are denoted with and , and stands for the cusp coefficient.
In order to restrain the model to finite dimensions a complementary error function with cut-off radius and width is used.
The velocity distribution is computed by solving the Eddington equation for the prescribed density profile.
For each removed halo two replacements are created: one with and another with . The cut-off radius and width are set to be linear functions of the virial radius of the original object:
The scale radius is determined indirectly by assuming that the concentration of an average halo of mass at redshift are given by the empirical relations found by Klypin et al. (2011):
With the concentration and virial radius set, the scale length can be determined.
The two variants corresponding to shallow and steep dark matter density profile have by construct the same scale radius. In order to have the same mass enclosed in the same virial volumes with different cusp coefficients the scale velocity is adjust accordingly. The values of the characteristic quantities associated to dark matter halos with the virial quantities presented in table 1 are listed in table 2.
The previously mentioned code can generate a disk component of a galaxy according to the cylindrical spatial distribution:
and cylindrical distribution of the radial dispersion:
The reader who is interested in more details on how the code constructs the galactic model is referred to the publication Widrow et al. (2008).
There are free parameters that have to be set using educated guesses based on available constraints on galactic structure at low mass scales. Most importantly, the mass of the stellar disk for each object is of its virial mass Such ratio is chosen to reflect the stellar-to-halo mass ratio, , inferred from the abundance matching analysis extrapolated to lower masses (Behroozi et al. (2013)). Past simulations of disk galaxy formation presented in the Kaufmann et al. (2007) paper associate to objects with maximum circular velocities in the 24-53 km range, scale height values in the 0.08-0.18 kpc range. Considering the previously mentioned results the values of the scale heights of all the objects is set to the minimum of 0.08 kpc. The choice of constant scale height is motivated by the expectation that at the masses of interest the disk thickness acquires minimum threshold value. The value is dictated by the balance between the thermodynamics of disk formation, such as the effect of the cosmic ionizing flux providing a minimum temperature floor of K, corresponding to a bulk velocity dispersion km/s, and the gravitational potential of the halo. For all the stellar disks introduced in the simulation the scale height was set to 0.08 kpc.
With the scale height and stellar mass to total mass ratio assigned, the scale length remains to be determined. The square of the scale length is naturally constrained to be proportional to the virial mass of the respective object. An observational study provided by Swaters et al. (2009) finds several dwarf galaxies with disk scale lengths as low as 0.33 kpc. Considering the inferred galaxy masses for the latter sample and the values assumed in previous work on tidal stirring based on various constraints (Kazantzidis et al. (2011a)), the reference value of was set for an object with a virial mass of . Using such reference value as a normalization the following relation is used for the calculation of the scale length as a function of galaxy mass:
For each pair of objects the parameters of the stellar component are listed in table 2.
From the original cosmological simulation, quantities such as the virial mass and virial radius are measured (table 1). The virial quantities are used for determining all the parameters required for the generation of the pairs of high resolution galaxies that will be used as replacements. Namely the following quantities are determined: () the dark matter halo quantities and () the stellar disk quantities. The remaining free parameters are the number of dark matter particles and star particles. For all the models created both numbers were set to one million.
2.3 Stability of the galaxies
Each of the twelve high resolution generated objects evolve in isolation for 1 Gyr prior to their introduction in the cosmological environment. The state of the objects at the end of the isolation run can be seen in figure 2 and 3. The purpose of these simulations is to verify the stability of the models and allow the system to remove instabilities prior to replacement. Initial discreteness noise is unavoidable, and it is known to seed waves in cold disk-like systems which could lead to mass redistribution (eg Mayer (2011)). These often manifest themselves as transient spiral waves, although the method adopted by our initial condition generator combined with the low self-gravity of the light disks embedded in a massive halo should minimize amplification of the waves. Furthermore, two-body relaxation between more massive dark matter particles and lighter stellar particles may also lead to some spurious evolution (Mayer et al. 2001). All isolation runs were performed using the pkdgrav2 code (Stadel (2001)). The softening lengths of the particles were set according to the equation:
And are listed in table 3.
The baryonic distributions of the constructed galaxy models with a cusp coefficient equal to unity are presented in figures 2 and 3. The previously mentioned images correspond to the objects at the end of the isolation run and at the moment of replacement. Figure 4 shows the cylindrically averaged surface density at the beginning of the isolation run and at the end. Small fluctuations are visible as small amplitude waves, which cause departure from perfect axisymmetry. Despite the existence of the local density fluctuations for both versions there are no significant changes happening on scales comparable to the scale length during the isolation run. Moreover, the analysis of the global mass distribution including the dominant dark matter component as shown by the circular velocity curves in figure 5 along with the radial distribution of the cumulative mass present in figure 6, do not show any appreciable changes. We conclude that the stability of the models is satisfactory.
2.4 Replacement method and numerical simulations
At the end of their evolution in isolation the galactic models replace the original objects in four steps. All bound particles of the original objects are removed after measuring their centre of mass (COM) velocity and position from the cosmological box at redshift 2.87. The particles of the newly created pairs of galaxies have their residual COM velocities and positions removed. Afterwards their individual positions and velocities are rotated with random matrices generated with the Arvo (1991). Finally the COM phase space coordinates are translated to match the COM phase space coordinates of the original objects that they replace. The resulting modified versions of the cosmological volume containing ErisDark are used as initial conditions for the gravity calculation runs down to redshift 0. ChaNGa was used to carry out these new cosmological runs. The code ChaNGa uses a fast hierachical oct-tree gravity solver, and is based on the Charm++ parallel programming infrastructure (Kale & Krishnan (1996)). It leverages the object-based virtualization and data-driven style of computation inherent in Charm++ to adaptively overlap communication and computation and achieve high levels of resource utilization on large systems (Jetley et al. (2008), Jetley et al. (2010)). The Charm++ load balancing infrastructure is used to distribute pieces of this tree across processors. The leaf nodes are buckets that contain several particles (usually 8 to 32) whose force calculations are collectively optimized. At each level of the tree, multipoles are calculated to speed distant force evaluations. Time adaptivity is achieved by assigning individual timesteps to particles. A special load balancer was developed by the UIUC Charm++ team to efficiently handle the resulting processor load imbalance (Menon et al. (2015)).
3 Results and data analysis
The resulting raw simulation data is analized using the AHF in the same manner as performed with the original ErisDark simulation data.
From the set of halos generated by the AHF code, the six pairs of satellites are identified through the indices of bound star particles.
As output of the AHF code, we use in our satellite analysis only the list of bound dark matter and star particles.
There are cases when the halo finder does not detect one of the objects at a particular epoch either due to the pericenter shock or
temporary interference with other large subhalos. Such events will appear as unavoidable interruption of the flow of data points on the analysis plots.
As soon as the object is identified again by the AHF code, the flow of data points in the analysis plots resumes.
The center of the satellite is defined as being the position of the stellar density peak of the respective object. From this point a 3D half light/mass radius can be defined as the radius of the sphere centered on the previously mentioned density peak which includes and excludes half of the mass amount of bound stellar particles. All the values of the half light radius associated with each object at the available evaluation epochs are presented in figure 7. The next step in the preparation of the satellite data for analysis is the alignment of the object. This is performed by defining a set of particles that include both dark matter and stars inside a sphere centered on the stellar density peak and with a radius equal to two half light radii. All the bound particles of the object including the ones exterior to the two half light radii sphere are rotated such that the moment of inertia tensor, corresponding to the set of stellar particles interior to the previously mentioned sphere, is diagonalized. This process is repeated for each individual object at each epoch of measurement.
Once the satellite is identified and analysed there are a number of diagnostics that are routinely measured, as described in the following subsections. In addition to dynamical masses, dark matter masses and stellar masses (see tables 4 and 5 for the results at the starting and final time), the key diagnostics used to quantify the morphology of the model galaxies are the shape as measured by the axis ratios, the ratio between 1D velocity dispersion and rotational velocity of the stars (note that we refer here to the actual stellar rotational velocity rather than to the circular velocity), and the surface density profiles of the stars. The characteristic velocities are measured inside () and around () the half-mass radius, which is defined for the whole 3D mass distribution.
Our classification scheme to decide whether or not the final remnant is a dSph is more stringent than the one used in (Kazantzidis et al. (2011a);Kazantzidis et al. (2013)). There two criteria were used simultaneously, namely that the object had to have and . With the latter criterion outliers from the typical structural properties of dSphs, such as the Tucana dSph, which has , would be included. Here instead we use only the criterion , which is indeed the case for nearly all observed dSphs and that, as we will see, automatically guarantees also that in our simulated remnants at . Finally, for brevity we will refer to ”cored models” when we talk about the models with inner slope of the dark halo , and to ”cuspy models” when we discuss models with an NFW profile.
3.1 Dimensions along the principal axes
With the object aligned and principal axes identified three maps of the projected stellar mass corresponding to
the three pairs of principal axes can be generated. Perpendicular to each of the vectors defining the axes directions
a square grid with size and a million square cells is defined.
For each cell a projected mass is defined that equals the mass of all bound star particles whose projection on
the grid is included in the respective cell.
Using the information in the projected mass maps, the dimensions of the object along the principal axes can be determined. In the process of dimensions calculations ellipses containing cells with different levels of mass density are being fitted. The distribution of these ellipses around the center of the image allows for the calculation of the ratio of the dimensions. For each of the three planes the ratio between the shortest and longest dimensions and is measured: .
After the previously presented analysis is performed the first available scalar indicator of the transformation of the object can be determined:
The time evolution of the previously defined quantity is presented in figure 8.
3.2 Velocity dispersion and tangential velocity
With the directions of the principal axes determined, the velocity dispersion is computed in the direction of each axis. For the previously mentioned quantities the stars considered for the calculations reside in a spherical volume of a radius of interest, centered on the satellite. The velocity dispersion of the system associated with radius is:
In order to compute the rotational velocity associated to the system, three cylindrical shell volumes are defined. Each of the volumes is defined with a central axis constrained to pass through the maximum density region of the object and to be parallel to a principal axis dimension. The inner and the outer radii of the volumes are 0.95 and 1.05 .
For all the stars bound to the object and included in the previously defined volume,
the component of their velocity tangential to any concentric circle centered on the axis of the cylinder is considered for the average.
The resulting average is the velocity , where .
3.3 Projected mass surface densities and cylindrically averaged surface mass densities
In addition to the quantities determined based on the bound amount of particles for the respective dwarf galaxies, the projected surface mass densities (PSMD) and cylindrically averaged surface mass densities (CASMD) have been determined. Unlike with the previous quantities where the nearby unbound particles are excluded from the calculations, in the PSMD and CASMD case they are included. The plots present in figures 2, 3, 10, 11, 12 and 13, are generated by projecting along the principal axes of the bound satellite all stellar particles inside a spherical volume of 10 kpc radius centered on the stellar density maximum of the galaxy. In the case of the CASMD calculations, the depth of the projection was reduced in order to include the stars with distances from the plane of projection smaller than 2.5 kpc. The resulting plots are available in figures 4, 14.
3.4 Evolution of satellites
3.4.1 Galaxy A
Originally the most massive of the studied objects, satellite A enters the host halo on an orbit which brings it very close to the ErisDark’s center.
Soon after the first pericenter passage, owing to the effect of dynamical friction
the dwarf galaxy settles on an eccentric orbit with pericenter distances smaller than 5 kpc and apocenter distances less than 50 kpc.
Before the end of the simulation the cuspy version of the object completes 15 orbits, while the cored version completes 14 orbits, as it can be seen from the top left panel of
figure 15. In both cases the satellite is found at z = 0 close to its maximum distance from the host.
The three dimensional half mass radius decreases from the original 1300 pc values to 220 pc and 110 pc for the NFW version and nonNFW version respectively.
The time evolution of the half-mass radius can be inspected in figure 7.
The tides reach the stellar disc and unbinds 10 of the stellar mass before the second pericenter is reached.
90 of the stellar mass of the cored dwarf is lost prior to completion of the fourth orbit (figure 16).
In the case of the counterpart with steep density profile, seven orbits have to be completed before loosing the same amount of stars.
Based on the projected mass distribution of the satellite, it can be observed that the difference in the halo profile has no clear impact on the evolution of the shape of the satellite (figures 10 and 11). The ratios between the shortest and longest dimension corresponding to both versions reaches the value 0.5 during the third pericenter passage and remains above this value throughout the simulation (figure 8).
From a kinematic perspective, the ratios between the rotational velocity and the dispersion associated with the respective half-light radius decrease to the 0.5 level differently for the two versions (figure 9). Before the first orbit is completed the cored satellite is described by a of 0.5. At the same time the ratio of the NFW variant decreases to values in the 1.0 - 1.7 range. Later during the third orbit, the value 0.5 is reached as well. It should be noted that the final size of the object with = 0.6 is comparable to the softening length of its dark matter particles. Therefore the very final stages in the evolution of the object should be regarded with caution as the response of the object to tides is affected by the limited resolution ((Moore et al., 1996)).
3.4.2 Galaxy B
Unlike the previous dwarf galaxy, the second object has a significantly wider orbit.
The apocenter distance varies around the 100 kpc value with a pericenter distance less than 10 kpc.
During the 11 Gyr timespan of the simulation the satellite completes six orbits, and at the epoch corresponding to z = 0 it is aproaching its 7th apocenter (figure 16).
The value of half mass radius shows large fluctuations in the case of the variant with shallow dark matter density profile, with large variations when the
object is in the proximity of the periapsis. The half light radius of the alternative presents smaller fluctuations.
Each of the versions of dwarf galaxies start with a half light radius of 1160 pc and end with a value of 1010 pc in the case of the cored dwarf and 900 pc
in the case of the cuspy dwarf (figure 7).
The stellar discs of the two versions start to lose significant amounts of stars during the second pericenter passage. By the end of the simulations the version with the higher central density has lost half of the initial stellar mass. In the case of the other variant this level of loss is reached before the 4th apocenter epoch. Moreover at z=0 90 of the original stellar mass is lost (figure 16).
The velocity ratio corresponding to the NFW profile reaches the level of 1.5 at redshift 0. At the same epoch the cored dwarf reaches a value of 0.5 (figure 8).
Spatially, the two versions of the object, maintain similar ratios of the minor and major dimensions until epoch 8 Gyr. Afterwards in the case of the variant with a cusp coefficient of 0.6, the ratio grows faster. The divergence in the velocity ratio appears at the same time; namely after the epoch corresponding to the 4th pericenter. Finally at z =0 , c/a fluctuates around 0.8 for the non-NFW variant and around 0.5 for the NFW variant (figure 9). When studying the projected density maps it is evident that the cupsy object maintains a two-armed spiral pattern while the cored dwarf approaches spherical symmetry (figures 10 and 12).
3.4.3 Galaxy C
The third most massive object is originally on a trajectory with a pericenter distance of 40 kpc and an apocenter distance of 110 kpc.
Soon after one orbit is completed it decays to a trajectory closer to the host. The following orbits are characterised by pericenter
distances in the range 10-20 kpc and apocenter distances in the 90-100 kpc range. As the previous satellite before the last epoch it completes six orbits (figure 15).
Kinematically, the characteristics of the two versions diverge after epoch 7.5 Gyr which roughly corresponds to the third apocenter passage (figure 8).
The indicator of the spheroidal transformation c/a clearly reflects the large difference in the transformation level (figure 9).
While the cuspy object maintains a ratio around 1.5, the cored object reaches values below 0.2.
The shape indicators of the two variants diverge earlier during the third pericenter passage and eventually converge to 0.7 values at z=0,
roughly the time of the last apocenter passage.
A visual inspection of the density map discloses a more complex case (figures 10 and 12). The shape of the stellar distribution in the inner regions resemble case by case, the distribution of the satellite B pair. At the same time the outer, low density and spatially more extended regions are very different. Particularly they are more isotropically distributed in the case of object C than in the case of object B for both versions of the respective object.
The stellar masses of the two versions of satellite C remain approximately constant until the 7 Gyr epoch. From this moment on the stellar mass decreases from 1.7 to 1.0 in the case of the cuspy dwarf, and to 4 in its cored counterpart. The original value of the half light radius for both versions was 810 pc, while the final values are 680 pc and 650 pc in the cuspy and cored dwarf, respectively (figure 8). Fluctuations similar to ones observed with object B are present, although with smaller amplitudes. Moreover they start to be visible at later epochs than in the case of object B.
3.4.4 Galaxy D
The object D is bound to ErisDark on an orbit with a pericenter distance of about 20 kpc and an apocenter distance in the 70-80 kpc range.
During the simulation time frame the dwarf galaxy passes through 7 pericenters, and at z = 0 it is approaching the 8th (figure 15).
The tides of the host galaxy start to significantly affect the stellar component just before 7 Gyr ago.
From this epoch onwards the cuspy object loses mass. During the same period its cored
counterpart loses 8 (figure 16).
Similarly to the previous cases the stellar component in the cored model begins to lose substantial mass
when the onset of large fluctuations in the values of the ratio is also observed,
reflecting the more impulsive response to tidal shocks relative to the cuspy case (figure 9).
The fluctuations are particularly enhanced for the cored version.
Close to the present epoch reaches for both variants of the model.
Similarly, the velocity ratio starts to diverge at the 7 Gyr epoch and converge to similar values of 0.5 at z = 0 (figure 8). Shapewise the final state of the two versions appear to be similar (figures 11 and 13). Moreover the evolution of the half light radius is similar for both cases, showing small fluctuations and a small decrease of the value from the original 640 pc values to 570 pc and 580 pc for the cuspy and cored object, respectively (figure 7).
3.4.5 Galaxy EF
The Dwarf galaxy EF, which is the result of a merger between two objects comparable in size, evolves on a wide orbit with the pericenter distance varying from 30 to 50 kpc and apocenter distances with values as high as 160 kpc (figure 15). Since infall the satellite completes three orbits. After the merger, the stellar mass of the object remains virtually constant for both variants (figure 16). The ratio / stabilizes around for the object resulting from the merging of two cuspy dwarfs. For the other variant, / continues to decrease until it reaches values close to 0.1 at 13 Gyr epoch (figure 9). In the case of the shape ratio c/a the object resulting from the merger of the cuspy subhalos has on average a higher value; after 12 Gyr the average c/a is around 0.8 for the cuspy pair and 0.7 for the cored pair (figure 8).
The original values of the half mass radius prior to the merging of the two objects were 420 pc for object E and 390 pc for object F. Immediately after the merger the radius of the resulting objects is 550 pc and during the remaining simulation time it increases to 610 pc, with no significant differences between the cuspy and core variants (figure 7). This highlights a clear difference with tidally stirred dwarfs, in which the half-light radius always decreases.
In the current work the evolution of five pairs of satellites entering a MW-sized halo has been presented.
Each pair consisted of two objects differentiated by a single parameter, namely the cusp coefficient of the dark matter halo.
In each pair, each of the objects was composed of a stellar disk component and a dark matter halo.
The virial mass of the objects spanned one order of magnitude from 3.3 to 3.83 .
The simulation commences at z =2.87 when all dwarf galaxies are in the proximity of their future host halo, yet outside its virial volume.
The original distances to the host center ranged from 75 to 165 kpc. Since the last major merger experienced by the host halo
takes place before z = 3 the described population of objects is representative of the earliest infall objects and thus the most tidally affected
population. Galaxies with later infall epochs will be affected by the tides on average to a lesser extent.
All of the objects with cored density profiles match our definition of dSphs before z=0. At the same time only two out of four objects with cuspy density profiles satisfy the criterion. The more efficient transformation of cored models is expected as their response to tides is more impulsive owing to their lower gravitational binding profile for comparable mass (Kazantzidis et al. (2013)). In particular the binding energy decreases much faster towards the center allowing the tides to disturb the stars found closer to the center of the satellite (Łokas et al. (2012)).
Figures 17, 18, 19 and 20 show the time evolution of various structural properties of all models, highlighting how the cored models in general fit better the range of values that classical dSphs and UFDs span for such properties (McConnachie (2012)). We caution, though, that there is a lack of objects with small half light radii, which may reflect the relatively small sample of initial conditions adopted here (see next section).
The half mass radius, velocity and shape ratios have significant fluctuations on orbital timescales. Their fluctuations appear to be in phase with orbital position and grow in amplitude as more stellar mass is unbound from the respective satellite. In general the objects with an original show significantly greater amplitudes in the fluctuations than the ones with , again reflecting the different binding energy profiles which place cored models in the more impulsive regime for the response to tides (which will be maximal near pericenter). The minimum of the fluctuations in half mass radius is reached when the object is near its periapsis, while the maxima is reached a few hundred million years afterwards.
In general the satellites with are also transformed faster into dwarf spheroidals than their counterparts. The final dark matter to stellar mass ratio are significantly higher in the cuspy remnants (see table 5), yet the final values in cored remnants are still high enough to match those of classical dSphs and UFDs. Assuming a stellar mass-to-light ratio , appropriate for an old stellar population, (table 5) would yield values in the range for the cored remnants. Instead velocity dispersions and dynamical masses of cuspy remnants indeed appear to exceed observed values in model B and D, because of the very high central dark matter density. Note this result is not coincidental rather it reflects the fact that, while the core remnants have radii of less than 1 kpc, the initial dark matter-to-stellar mass ratio was very high even inside the central kpc in order to match the constraints from abundance matching with the chosen density profiles of dark matter.
All the tidally stirred remnants are found at kpc from the host at . Within such a distance from the MW and M31 all the galaxies in the mass range of the remnants (stellar masses ) are dSphs (or UFDs, which structurally appear similar to classical dSphs, see McConnachie (2012)). One is thus tempted to conclude that this is evidence of the fact that the most massive satellites falling early inside the host, which is the type of objects we are studying here, have predominantly shallow dark matter profiles. In other words, dwarf galaxies with the properties of the cuspy remnants that did not transform fully into dSphs are not observed in the inner halo of the MW and M31, which means cuspy progenitors must be rare (the only irregular galaxies in the inner halo, the LMC and SMC, are a few magnitudes more luminous than any of our remnants). Note that the latter inference is well grounded even with the small sample in this paper since specifically this sample contains the objects that have experienced the longest tidal evolution possible inside the main host (early infall population). We thus expect a larger number of such incompletely transformed objects with disk-like features in the much larger sample with varying infall time that will be presented in Paper II.
Note that idealized simulations already hinted at such a more efficient transformation of the cored models, still they never found such a strong effect because the mass ratio adopted between the satellite and the main host was at least an order of magnitude higher than in this work. The reason is that here the host halo is relatively low in mass (closer to than to when the most massive satellites fall in at , an effect that only a fully cosmological simulation can account for. The case for cored progenitors is suggested by the notion that dwarf galaxy satellites which had the highest star formation rates, as required to reach a stellar mass in the range of already at , are those in which feedback was most effective to turn cusps into cores (Governato et al. (2010); Governato et al. (2012); Shen et al. (2014); Oñorbe et al. (2015)). These are the objects in which the star formation (SF) efficiency, defined as ratio between stellar mass assembled and dark halo mass, was particularly high (Madau et al. (2014)). A high star formation efficiency in classical dSphs versus similarly low mass field dIrrs is indeed seen in SF histories based on hi-res color magnitude diagrams (eg Skillman et al. (2014); Aparicio et al. in preparation).
Interestingly, gas-poor dwarfs with disk features have been detected in clusters (Lisker et al. (2006)) and have been explained as harassed galaxies that have accreted onto the cluster core only recently, thus have not been able to transform fully into dwarf spheroidals. It is tempting to speculate that such a population would be present in clusters even if the progenitors have dark matter cores before infall since a much larger fraction of galaxies accrete at in clusters as opposed to for satellites into galaxy-sized halos. This is strongly suggested by the evolution of some of our models, which also in the cored variant transform fully into dSphs only after nearly 10 billion years (eg see models B and D in figure 9).
Moreover in order for the velocity ratio of a tidally affected object to reach the threshold for a spheroidal it must lose more than 25 of its original stellar mass (figure 17). The total fractional mass loss was 1/3 at a minimum (in the cuspy version of model satellite D). This occurs for both cuspy and core variants, and implies a prominent tidal stream has to form if the object is transformed into a dSph. Detecting tidal features or tails might be difficult though because, as our simulations show, tails and streams show up at a surface density that is between and times lower than the central surface density of the object (14, 21), which corresponds to 7-10 magnitudes lower relative to the central surface brightness. The lowest values apply to cuspy models in which stellar removal is less efficient (figure 14 and 21). This confirms earlier results of Mayer et al. (2002). The transition to a marginally bound tidal tail and fully unbound stream following the orbital path is marked by the flattening of the profile at large radii, typically occurring at 3-6 kpc, or 3-6 times the half-light radius (figures 14 and 21), confirming past work with a variety of simulation techniques and models (eg Johnston et al. (1999); Mayer et al. (2002)). This means a very wide field of view would be needed to unequivocally detect the tidal debris surrounding the dwarf in addition to very deep photometry. The actual shape of the stellar profile varies also with projection and over time (figure 14 ), but the basic facts just highlighted hold. We also notice that in some cases the change of slope of the surface density profile between the inner bound and outer unbound stellar component is more evident in the cuspy variant, probably because in the cored models tidal mass loss occurs all the way to the very center, leading to more continuous removal of stars as a function of radius.
The evolution of the observable structural parameters for the simulated objects shown in figures 18, 19 and figure 20 also highlights that cored models can evolve enough in stellar mass and radius to end up in the region populated by Ultra Faint Dwarfs (UFDs), confirming a previous claim of Łokas et al. (2012) obtained with idealized simulations. In general the objects evolve from the proximity of the locus populated by dwarf irregulars towards the locus of ultra-faint dwarf spheroidals by passing through the classical dwarf spheroidal locus. Except for the satellite A pair all galaxies can be found in the locus of classical dwarf spheroidals and the transition region between the faintest irregular dwarf galaxies and spheroidals. Rather than being necessarily reionization fossils, as it has been claimed in the literature (Ricotti (2010)), some UFDs might be thus just more extreme outcomes of tidal stirring. This scenario has the attractive aspect of allowing to explain naturally why there is a substantial continuum between the properties of classical dSphs and those of UFDs. Such a population of UFDs would have cored profiles, while on the contrary a population of pristine UFDs formed by gas collapse in mini-halos would likely maintain its original cuspy profile as feedback from star formation would have a negligible effect (Peñarrubia et al. (2012)).
Despite the low final stellar mass such a population would yield no tension with the proposed core formation mechanism via feedback (Peñarrubia et al. (2012)) because it would be the remains of a population of much more massive dwarf galaxies that were severely affected by the tides of the host. Indeed among our remnants the UFD-like object originates from one of the most massive infalling dwarfs, whose initial luminosity would be comparable to Fornax, exactly the mass scale at which the effect of feedback should be important (Governato et al. (2012); Shen et al. (2014)). Finding UFDs with cored profiles is thus a potential test of our scenario. UFDs forming from an early infall population of satellites would also be consistent with the notion that these galaxies typically lack an intermediate or young population (Brown et al. (2012)) as star formation would be truncated soon after infall as gas is rapidly removed by tides and ram pressure when the cosmic ionizing background is still at its peak amplitude (Mayer et al. (2007)).
The smallest two objects replaced merge while orbiting the host galaxy at significantly larger pericenter and apocenter distances in comparison to the other satellites. The effect of the tides is minimal with insignificant stellar mass loss and small effects in the shape and velocity ratios for the merged objects with = 1.0. In the case of the = 0.6 mergers the velocity ratio further decreases due to the gravitational interaction with the host. If one considers the condition / as the definition of a dwarf spheroidal, then the merging pair with = 0.6 clearly reaches this state at the end. At the same time the merger remnant of galaxies with = 1.0 fails the velocity criterion, although the later appears shapewise to be closer to a spheroid than the former. The formation of a significantly round object albeit with rotation still dynamically important cannot be caused by tidal stirring only since in that case the transformation of the initial disky shape is correlated with the redistribution of angular momentum due to tidally induced instabilities and tidal heating.
Hence finding a dwarf with photometric properties analogous to those of classical dSphs but with a significant rotational velocity might be the signature of dwarf-dwarf mergers. Indeed low luminosity ellipticals, whose kinematics is reproduced well by mergers of disks with intermediate masses, do retain significant rotation (Cappellari et al. (2011)). One of such rare objects might indeed be the Tucana dSph, which is on a very wide trajectory around the primary and exhibits appreciable rotation (Fraternali et al. (2009)).
Figure 14 also shows that the change of slope in the case of the merger remnant occurs much more gradually, with a clear flattening not appearing before 6-7 kpc from the center. This is not surprising as in in this case the morphological evolution is not accompanied by significant tidal mass removal. The different shape of the profile relative to tidally stirred remnants is very interesting because it hints at a further possible way to distinguish between a dSph formed by tidal stirring and one formed via dwarf galaxies mergers. The presence of a swift change from one exponential to another should be the signature of a tidally affected stellar system.
5 Conclusion and Caveats
To the authors’ knowledge the current work addresses for the first time the importance of the dark matter density profile for the dynamical evolution and morphological transformation of disky satellites into spheroidals in a cosmological context. Therefore the work bridges the gap between detailed simulations of dwarf-host interaction and fully cosmological hydrodynamical simulations of MW-sized galaxies which do not have at the moment the resolution to reliably address this subject. This is an important step forward since the assembly history of the host, the infall orbits of the satellites and the structure of the potential well of the host, which all concur to determine the effect of tidal shocks, are all self-consistently taken into account. As a by product, dwarf-dwarf interactions are also take into account, and we witnessed a particular example of an interaction which gave rise to a merger.
An important overall conclusion of this work is that tidal stirring is confirmed to be an effective mechanism for transforming disky field dwarfs into objects with properties resembling dwarf spheroidals and even Ultra Faint Dwarfs (UFDs). Most importantly, it now appears clear that the morphological transformation induced by tidal stirring is a natural and inevitable consequence of hierarchical structure formation. The transformation is enhanced when the satellites have shallow dark matter density profiles, as in our cored models, confirming previous results based on idealized simulations of individual dwarf-host interaction (Kazantzidis et al. (2013)). In the latter case, both tidal mass loss of the stellar component as well as tidal heating and internal redistribution are augmented as the response of the galaxy is impulsive down to its central regions. As a consequence remnants can become as faint and puny as UFDs, which does not happen among the set of cuspy models.
We also found that major mergers of satellite galaxies with shallow density profiles can effectively evolve into spheroidals on wide orbits. Since the merger happens before the satellites are accreted by the host, when they have lower velocities,it naturally leads to a dwarf spheroidal orbiting far from the center of the host. Especially in the cuspy case the remnant retains significant rotation despite achieving a spheroidal shape, which would not happen in tidally stirred remnants This is a plausible scenario for the origin of distant dSphs such as Tucana and Cetus, and confirms a previous claim based on idealized simulations (Kazantzidis et al. (2011a)).
Another important conclusion of this work is the identification of two signatures that could separate dwarf spheroidals formed due to a merger on large orbits from the ones formed due to the strong tides generated by the host galaxy. The first signature is the presence of a tidal stellar stream with a total mass comparable to the associated satellite. The second signature can be found in the stellar density/surface brightness profile of the object. It should appear as a fast transition on lengths of a few hundreds parsecs from the steep central exponential to a significantly shallower one. The objects that are not the source of such a stream or do not have a fast departure from the central brightness exponential decay most probably did not become spheroidal due to the gravity of their current host. They were probably formed through the merger of comparable galaxies.
There are of course important caveats in the simulations, which primarily stem from the fact that they are carried out in the dissipationless approximation, namely without treating hydrodynamics, star formation and feedback. One caveat is the inexistence of the baryonic disk of the primary. This would increase significantly the central density of the host, within 10-20 kpc, enhancing tidal mass loss and tidal stirring for satellites with orbits having pericenters smaller than kpc D’Onghia et al. (2010). Given the orbits of the satellites analysed in this paper, this effect would play a role at least for some of them (eg model A and model B, which have the smallest pericenters). The neglected effect should enhance the tidal transformation, but may lead to complete tidal disruption of objects that would have otherwise ended up as UFDs, or produce UFDs out of cuspy progenitors (Peñarrubia et al. (2010)). Another caveat is the lack of a a gaseous component in the disky dwarf progenitors, but previous work has shown that this would be readily removed by ram pressure in the galactic diffuse corona over 1-3 Gyr for the mass range considered here ( km/s), especially at , when the cosmic photoionizing background maintains the ISM warm and loosely bound, thus suppressing further star formation even for the gas that has not been removed yet (Mayer (2010)).
Furthermore, our modelling approach relies on the assumption that the progenitors of dwarf spheroidal satellites started out their life as disky dwarf galaxies. Cosmological hydrodynamical simulations show that this is a well justified assumption for dwarf galaxies forming in halos in the range (Governato et al. (2010); Shen et al. (2014)), while at lower masses field dwarfs can have a rather thick and turbulent disk, with a ratio (Shen et al., in preparation), which appears consistent with the kinematics of the faintest dIrrs in the Local Group (Mateo (1998); McConnachie (2012)). A thicker, kinematically hotter disk could have a different response to tidal shocks. On the one hand the triggering of non-axisymmetric instabilities, such as bars and spiral arms, would be weaker due to the reduced self-gravity, on the other hand the disk binding energy will be lower, so that the direct effect of tidal heating could be stronger. Kazantzidis et al. (2011a) have considered thicker disks in the wide range of initial conditions they explored in their tidal stirring simulations, finding small differences in the transformation as the initial disk thickness was varied by a factor 3, with only a slight tendency of thicker disks to transition faster to a nearly isotropic configuration. However, future cosmological hydrodynamical simulations with enough resolution to study the morphological evolution of satellites with the same level of detail as done here will be required to assess the impact of disk thickness and kinematics, of both gas and stars, on the efficiency of the transformation.
Finally, here we have considered a relatively small sample of satellites as we focused on the earliest infalling objects with a sizable mass, large enough to host luminous dwarfs and be sufficiently well resolved. The advantage of the small sample is that we could study in detail the different dynamical evolution between cuspy and cored variants of the same object. In addition, since these early infallers are exposed for the maximum amount of time possible to the tides of the primary we can safely assume that the results provide an upper limit on how deeply the structure of the dwarfs can be modified by tidal shocks in cored versus cuspy variants. As a result, we inferred that the population of satellites of the Milky Way should be dominated by remnants of cored objects otherwise we would see a population of dwarfs with intermediate properties between dwarf irregulars and dwarf spheroidals at kpc galactocentric distances, which is not seen neither in the MW or the M31 halo. In a forthcoming paper we will study a much larger sample of satellites which includes objects that are accreted at low redshifts. The size of the sample that we will present will be of the same order of the number of luminous satellites known so far to lie within the putative virial radius of the Milky Way, thus allowing to make statistically meaningful statements on the evolution of the entire population of satellites.
- Arvo (1991) Arvo, J. 1991, Graphics Gems II, 355
- Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
- Bird et al. (2013) Bird, J. C., Kazantzidis, S., Weinberg, D. H., et al. 2013, ApJ, 773, 43
- Brooks et al. (2013) Brooks, A. M., Kuhlen, M., Zolotov, A., & Hooper, D. 2013, ApJ, 765, 22
- Brown et al. (2012) Brown, T. M., Tumlinson, J., Geha, M., et al. 2012, ApJL, 753, L21
- Bryan & Norman (1998) Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80
- Cappellari et al. (2011) Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 416, 1680
- Colpi et al. (1999) Colpi, M., Mayer, L., & Governato, F. 1999, ApJ, 525, 720
- Di Cintio et al. (2013) Di Cintio, A., Knebe, A., Libeskind, N. I., et al. 2013, MNRAS, 431, 1220
- Diemand et al. (2007) Diemand, J., Kuhlen, M., & Madau, P. 2007, ApJ, 667, 859
- D’Onghia et al. (2010) D’Onghia, E., Springel, V., Hernquist, L., & Keres, D. 2010, ApJ, 709, 1138
- Dressler (1980) Dressler, A. 1980, ApJ, 236, 351
- Dubinski (1999) Dubinski, J. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 182, Galaxy Dynamics - A Rutgers Symposium, ed. D. R. Merritt, M. Valluri, & J. A. Sellwood, 491
- Fraternali et al. (2009) Fraternali, F., Tolstoy, E., Irwin, M. J., & Cole, A. A. 2009, A&A, 499, 121
- Gill et al. (2004) Gill, S. P. D., Knebe, A., Gibson, B. K., & Dopita, M. A. 2004, MNRAS, 351, 410
- Governato et al. (2010) Governato, F., Brook, C., Mayer, L., et al. 2010, nature, 463, 203
- Governato et al. (2012) Governato, F., Zolotov, A., Pontzen, A., et al. 2012, MNRAS, 422, 1231
- Grebel et al. (1999) Grebel, E. K., Seitzer, P., Dolphin, A. E., et al. 1999, in Bulletin of the American Astronomical Society, Vol. 31, American Astronomical Society Meeting Abstracts, 1380
- Guedes et al. (2011) Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742, 76
- Guedes et al. (2013) Guedes, J., Mayer, L., Carollo, M., & Madau, P. 2013, ApJ, 772, 36
- Jetley et al. (2008) Jetley, P., F., G., Mendes, C., Kale, L., & Quinn, T. 2008, IEEE International Parallel and Distributed Processing Symposium, 1
- Jetley et al. (2010) —. 2010, IEEE International Conference for High Performance Computing, Networking, Storage and Analysis (SC), 1
- Johnston et al. (1999) Johnston, K. V., Sigurdsson, S., & Hernquist, L. 1999, MNRAS, 302, 771
- Kale & Krishnan (1996) Kale, L. V., & Krishnan, S. 1996, in Parallel Programming using C++, ed. G. V. Wilson & P. Lu (MIT Press), 175–213
- Kaufmann et al. (2007) Kaufmann, T., Wheeler, C., & Bullock, J. S. 2007, MNRAS, 382, 1187
- Kazantzidis et al. (2011a) Kazantzidis, S., Łokas, E. L., Callegari, S., Mayer, L., & Moustakas, L. A. 2011a, ApJ, 726, 98
- Kazantzidis et al. (2013) Kazantzidis, S., Łokas, E. L., & Mayer, L. 2013, ApJL, 764, L29
- Kazantzidis et al. (2011b) Kazantzidis, S., Łokas, E. L., Mayer, L., Knebe, A., & Klimentowski, J. 2011b, ApJL, 740, L24
- Klimentowski et al. (2007) Klimentowski, J., Łokas, E. L., Kazantzidis, S., et al. 2007, MNRAS, 378, 353
- Klypin et al. (2011) Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102
- Knollmann & Knebe (2009) Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608
- Kuhlen et al. (2013) Kuhlen, M., Guedes, J., Pillepich, A., Madau, P., & Mayer, L. 2013, ApJ, 765, 10
- Kuijken & Dubinski (1995) Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341
- Lisker et al. (2006) Lisker, T., Grebel, E. K., & Binggeli, B. 2006, AJ, 132, 497
- Łokas et al. (2012) Łokas, E. L., Kazantzidis, S., & Mayer, L. 2012, ApJL, 751, L15
- Madau et al. (2014) Madau, P., Weisz, D. R., & Conroy, C. 2014, ApJL, 790, L17
- Mastropietro et al. (2005) Mastropietro, C., Moore, B., Mayer, L., et al. 2005, MNRAS, 364, 607
- Mateo (1998) Mateo, M. L. 1998, ARAA, 36, 435
- Mayer (2010) Mayer, L. 2010, Highlights of Astronomy, 15, 193
- Mayer (2011) Mayer, L. 2011, in EAS Publications Series, Vol. 48, EAS Publications Series, ed. M. Koleva, P. Prugniel, & I. Vauglin, 369–381
- Mayer (2012) Mayer, L. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 453, Advances in Computational Astrophysics: Methods, Tools, and Outcome, ed. R. Capuzzo-Dolcetta, M. Limongi, & A. Tornambè, 289
- Mayer et al. (2001a) Mayer, L., Governato, F., Colpi, M., et al. 2001a, ApJ, 559, 754
- Mayer et al. (2001b) —. 2001b, ApJL, 547, L123
- Mayer et al. (2007) Mayer, L., Kazantzidis, S., Mastropietro, C., & Wadsley, J. 2007, nature, 445, 738
- Mayer et al. (2006) Mayer, L., Mastropietro, C., Wadsley, J., Stadel, J., & Moore, B. 2006, MNRAS, 369, 1021
- Mayer et al. (2002) Mayer, L., Moore, B., Quinn, T., Governato, F., & Stadel, J. 2002, MNRAS, 336, 119
- McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4
- Menon et al. (2015) Menon, H., Wesolowski, L., Zheng, G., et al. 2015, Computational Astrophysics and Cosmology, 2, 1
- Moore et al. (1996) Moore, B., Katz, N., & Lake, G. 1996, ApJ, 457, 455
- Oñorbe et al. (2015) Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, ArXiv e-prints, arXiv:1502.02036
- Peñarrubia et al. (2010) Peñarrubia, J., Benson, A. J., Walker, M. G., et al. 2010, MNRAS, 406, 1290
- Peñarrubia et al. (2012) Peñarrubia, J., Pontzen, A., Walker, M. G., & Koposov, S. E. 2012, ApJL, 759, L42
- Pillepich et al. (2014) Pillepich, A., Kuhlen, M., Guedes, J., & Madau, P. 2014, ApJ, 784, 161
- Pontzen & Governato (2012) Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464
- Rashkov et al. (2013) Rashkov, V., Pillepich, A., Deason, A. J., et al. 2013, ApJL, 773, L32
- Ricotti (2010) Ricotti, M. 2010, Advances in Astronomy, 2010, 33
- Shen et al. (2014) Shen, S., Madau, P., Conroy, C., Governato, F., & Mayer, L. 2014, ApJ, 792, 99
- Skillman et al. (2014) Skillman, E. D., Hidalgo, S. L., Weisz, D. R., et al. 2014, ApJ, 786, 44
- Stadel (2001) Stadel, J. G. 2001, PhD thesis, UNIVERSITY OF WASHINGTON
- Swaters et al. (2009) Swaters, R. A., Sancisi, R., van Albada, T. S., & van der Hulst, J. M. 2009, A&A, 493, 871
- Trujillo-Gomez et al. (2015) Trujillo-Gomez, S., Klypin, A., Colín, P., et al. 2015, MNRAS, 446, 1140
- Widrow & Dubinski (2005) Widrow, L. M., & Dubinski, J. 2005, ApJ, 631, 838
- Widrow et al. (2008) Widrow, L. M., Pym, B., & Dubinski, J. 2008, ApJ, 679, 1239