THE CATERPILLAR PROJECT: A LARGE SUITE OF MILKY WAY SIZED HALOS

The Caterpillar suite was run using P-Gadget3 and Gadget4, tree-based N-body codes based on Gadget2 (Springel 2005). For the underlying cosmological model we adopt the ΛCDM parameter set characterized by a Planck cosmology given by, ${{\rm{\Omega }}}_{m}=0.32$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.68$, ${{\rm{\Omega }}}_{{\rm{b}}}=0.05$, n_s=0.96, ${\sigma }_{8}=0.83$, and Hubble constant, H₀ = 100 h km s⁻¹ Mpc⁻¹ = 67.11 km s⁻¹ Mpc⁻¹ (Planck Collaboration et al. 2014). All initial conditions were constructed using music (Hahn & Abel 2011). We identify dark matter halos via rockstar (Behroozi et al. 2013) and construct merger trees using Consistent-Trees (Behroozi et al. 2012). rockstar assigns virial masses to halos, ${M}_{\mathrm{vir}}$, using the evolution of the virial relation from Bryan & Norman (1998) for our particular cosmology. At z = 0, this definition corresponds to an overdensity of 104× the critical density of the universe. We have modified rockstar to output all particles belonging to each halo so we can reconstruct any halo property in post-processing if required. We have also improved the code to include iterative unbinding (see Section 2.5). In this work, we restrict our definition of virial mass to include only those particles which are bound to the halo.

Initially a parent simulation box (see Figure 1) of width 100 ${h}^{-1}$ Mpc was run at ${N}_{{\rm{p}}}={1024}^{3}$ (${m}_{{\rm{p}}}=8.72\times {10}^{7}\ {h}^{-1}$ ${M}_{\odot }$) effective resolution (see music/P-Gadget3 parameter files on project website) to select viable candidate halos for re-simulation (i.e., ∼10,000 particles per host). The candidates for re-simulation were selected via the following mass and isolation criteria.

• 1.  
  Halos were selected between 0.7 × 10¹² ${M}_{\odot }$ $\;\leqslant \;$ ${M}_{\mathrm{vir}}$ $\;\leqslant$ 3 × 10¹² ${M}_{\odot }$ (Smith et al. 2007; Xue et al. 2008; Tollerud et al. 2012; Boylan-Kolchin et al. 2013; Sohn et al. 2013; Piffl et al. 2014; Guo et al. 2015; Penarrubia et al. 2015; see Wang et al. 2015 for a review).
• 2.  
  No halos with ${M}_{\mathrm{vir}}$ ≥ 7 $\times \ {10}^{13}\;$ ${M}_{\odot }$ within 7 Mpc (Li & White 2008; van der Marel et al. 2012). We avoid halos near large clusters which would greatly enhance our Lagrangian volumes, making our ability to run simulations at our desired resolution impossible.
• 3.  
  No halos with ${M}_{\mathrm{vir}}$ ≥ 0.5 × M_host within 2.8 Mpc (Karachentsev et al. 2004; Tikhonov & Klypin 2009). We currently avoid pairs in our sample owing to the difficulty of running them at our desired resolution at the present time. We have nevertheless selected equivalent pairs of our current isolated sample but will be examining those in future work.

Zoom In Zoom Out Reset image size

This results in 2122 candidates being found (from an original sample of 6564 within the specified mass range). We use an extremely weak selection over merger history such that we require no halo to have had a major merger (1:3 mass ratio) since z = 0.05 (<5%). Our overall aim is to construct a representative sample of 10¹²${M}_{\odot }$ halos and not specifically require Milky Way analogs a priori as has been done in previous studies (Diemand et al. 2007; Springel et al. 2008; Stadel et al. 2009; Garrison-Kimmel et al. 2014b; Sawala et al. 2014). This also allows us to apply statistical tools to constrain semi-analytic models in future work (e.g., Gómez et al. 2014). We place our halos into three mass bins, with the largest number of halos centered on the most likely mass for the Milky Way (M₂₀₀ = ${1.6}_{-0.4}^{+0.5}\times {10}^{12}$ ${M}_{\odot }$, Piffl et al. 2014);

$$\begin{eqnarray*}{M}_{i}=\{\begin{array}{ll}0.7{\rm{\mbox{--}}}1.0\times {10}^{12} & {M}_{\odot }:7\ \mathrm{halos}\\ 1.0{\rm{\mbox{--}}}2.0\times {10}^{12}\ & {M}_{\odot }:46\ \mathrm{halos}\\ 2.0{\rm{\mbox{--}}}3.0\times {10}^{12}\ & {M}_{\odot }:7\ \mathrm{halos}.\end{array}\end{eqnarray*}$$

For this paper we are only considering a subset of the total sample in preparation, specifically 21 halos within the 1–2 × 10¹² ${M}_{\odot }$ mass range and 3 halos within the $0.7-1.0\times {10}^{12}\;$ ${M}_{\odot }$ mass range.

As has been highlighted by Onorbe et al. (2014), a great deal of care has to be taken when carrying out re-simulations of this kind so as to avoid contamination of the main halo of interest by low-resolution particles at z = 0. If mass from low resolution particles contributes more than ∼2% of the total host mass there can be offsets to estimates of the halo profile, shape, spin, and especially gas properties in hydrodynamic runs. To avoid contamination in our sample we custom built a Python GUI (using TraitsUI), Caterpillar Made Easy (cme), for running and analyzing cosmological simulations (for both single and multi-mass simulations). This tool allowed us to carry out an extensive contamination study (i.e., using ∼264 low resolution test halos with a particle mass of ∼10⁷ ${M}_{\odot }$) specifically for the halos to be re-simulated. We have automated the monotony of constructing hundreds of qualitatively similar but quantitatively distinct cosmological simulations with the added benefit of being able to interactively select over initial condition parameters, cosmologies, halo finders and merger trees. This procedure was carried out self-consistently across all runs, allowing for a systematic study of which simulation parameters produce the most computationally inexpensive to run, uncontaminated halos.

Using cme we tested 11 Lagrangian geometries (e.g., convex hull, ellipsoid, expanded ellipsoids, cuboids, and expanded cuboids) so as to ensure a sphere of radius ∼1 ${h}^{-1}$ Mpc exists of purely uncontaminated (high-resolution) particles centered on the host halo at z = 0. Our need to run 11 different Lagrangian geometries for each halo is motivated by the fact that the geometries vary substantially from halo to halo (due partially to their varied merger histories) and we wished to minimize the computation cost while achieving our contamination goals. It must also be highlighted that unlike many other studies, we did not select one Lagrangian geometry for all halos but used a specific geometry for a given halo depending on the needs of its simulation.

In Table 1 we show the various geometries we used for constructing our initial conditions. We modified music to be able to produce expanded Lagrangian volumes rather than the bounded volumes with which it was originally published. In Figure 2 we show four examples of Lagrangian geometries for halos selected for re-simulation. In some cases the geometries are reasonably compact, allowing for the traditional minimum cuboid enclosing to be used. Some larger regions, however, are extremely non-spherical (e.g., bottom right of Figure 2) and so a minimum ellipsoid was used. For each geometry we take the enclosed volume at z = 0 denoted by either $4\times {R}_{{\rm{vir}}}(z\quad =\quad 0)$ and $5\times {R}_{{\rm{vir}}}(z=0)$. We run each of these halos to z = 0, run our modified rockstar and determine at what distance the closest low resolution or contamination particle (type = 2) resides in each case. With the knowledge that the high-resolution volume distance decreases at higher levels of refinement (Onorbe et al. 2014), we ensure a minimum contamination distance of ∼1 ${h}^{-1}$ Mpc at z = 0 at our lowest resolution re-simulation, with the desire to have uncontaminated spheres of radius, ∼1 Mpc at our highest resolution re-simulation. In cases where four times the virial radius enclosure created contaminated halos but five times the virial radius created too large a simulation to run, we opted for an expanded ellipsoid of the minimum enclosing ellipsoid. In some cases there were a handful of offending particles far away from the primary Lagrangian volume (e.g., Figure 2 top right and bottom left) making no standard geometry enclosure feasible, expanded or otherwise. Here we trimmed the Lagrangian volume by hand and simulated the new geometry to z = 0 to ensure it had no contamination. Traditionally these types of halos are avoided but since we do not want to bias our sample, we dealt with complicated geometries in this specialized manner and have included them in our sample. Using this tailored approach, our highest resolution runs obtain very large, high-resolution regions with spheres of radius ∼1.4 ± 0.4 Mpc of solely high-resolution particles.

Zoom In Zoom Out Reset image size

Table 1.  The Contamination Suite Used on the First Refinement Level (i.e., levelmax = 11) for Every Halo in the Caterpillar Suite

Name	${{nR}}_{\mathrm{vir}}(z=0)$ a	Geometry	Factorb
CA4	4	Convex Hull	⋯
CA5	5	Convex Hull	⋯
EA4	4	Bounded Ellipsoid	⋯
EA5	5	Bounded Ellipsoid	⋯
EX4	4	Expanded Ellipsoid	1.05
EX5	5	Expanded Ellipsoid	1.05
EB4	4	Expanded Ellipsoid	1.1
EC4	4	Expanded Ellipsoid	1.2
BA4	4	Minimum Cuboid	⋯
BA5	5	Minimum Cuboid	⋯
BB4	4	Expanded Cuboid	1.1

Notes.

^aThe multiple of the z = 0 virial radius that we used to construct the Lagrangian volume. ^bThe factor by which we increased the original volume (e.g., 1.2 means the ellipsoid was expanded by 20% in size. A dash represents the minimum ellipsoid/cuboid/hull exactly). These values were arbitrarily chosen with the only requirement being that the initial condition files were not overly large in size (i.e., a few hundred megabytes at lx11).

Download table as:  ASCIITypeset image

In Figure 3 we show box plots of the median contamination distance and respective quartiles for all 264 of our test halos using each of our selected geometries. Typically, the best performing geometry (i.e., the largest uncontaminated volume with the cheapest computational expense) was the expanded ellipsoid which enclosed all particles within 4 or 5 times the virial radius of the host in the parent simulation at z = 0.

Zoom In Zoom Out Reset image size

Starting from our parent simulation resolution, we re-simulated each halo at iteratively higher resolutions (a factor of 8x increase in particle number for each level) to ensure we did not obtain contaminated particles within the host halo (the uncontaminated volume shrinks with an increase in the ratio of the zoom-in resolution to that of the parent simulation resolution). In Figure 4, we show the dark matter distributions for iteratively higher resolutions of the same halo. One can clearly identify the same subhalos across all resolutions, indicating the qualitative success of our numerical techniques. Regarding computational resources, each halo at our highest resolution took between 150 and 300 K hours on TACC/Stampede and occupies ∼5–10 TB of storage for both the raw hdf5 snapshots and halo catalog. Table 2 shows the mass and spatial resolution for each of our refinement levels. Our softening length is $\epsilon \;\sim$ 76 ${h}^{-1}$ pc for our fiducial resolution.

Zoom In Zoom Out Reset image size

Table 2.  The Resolution Levels of the Caterpillar Suite

lx	Np	mp
		(${M}_{\odot }$)	(${h}^{-1}$ pc)
15	327683	$3.7317\times {10}^{3}$	38
14	163843	$2.9854\times {10}^{4}$	76
13	80963	$2.3883\times {10}^{5}$	152
12	40963	$1.9106\times {10}^{6}$	228
11	20483	$1.5285\times {10}^{7}$	452
10	10243	$1.2228\times {10}^{8}$	904

Note.

lx represents represents the effective resolution (N_p =(2^X)³) of the high resolution region given by parameter, levelmax in music. m_p is the particle mass and is the Plummer equivalent gravitational softening length. The parent simulation parameters are also shown in the last row. Only a select sample of halos are run at resolution level lx15. These runs will be presented in future works.

Download table as:  ASCIITypeset image

We space our snapshots (320 per simulation) in the logarithm of the expansion factor until z = 6 (∼5 Myr/snapshot) and then linear in expansion factor down to z = 0 (∼50 Myr/snapshot). The motivation for this piecewise stitching of the two temporal schemes is twofold. At z > 6, we enter the era of mini-halo formation and the reionization epoch. If we wish to model the transport of Lyman–Werner (LW) radiation semi-analytically from the first mini-halos (e.g., Agarwal et al. 2012), we require a temporal resolution on par with the mean free path of LW photons in the intergalactic medium and the lifetime of a massive Population III star ($\leqslant 10$ Myr). Second, we also wish to resolve the disruption of low mass dwarf galaxies at low redshift, which requires a temporal resolution of order ∼50 Myr (e.g., Segue I has a disruption timescale of ∼50 Myr). These timescales are also required if one is attempting to determine subhalo orbital pericenters that can be input into semi-analytic models of tidal disruption (e.g., Baumgardt & Makino 2003). While we intend to use the capabilities offered by finely sampled snapshots in future work, this initial paper primarily focuses on the z = 0 halo properties.

rockstar is able to find any overdensity in 6D phase space, including both halos and streams. To distinguish gravitationally bound halos from other phase space structures, rockstar performs a single-pass energy calculation to determine which particles are gravitationally bound to the halo. Overdensities where at least 50% of the mass is gravitationally bound are considered halos, with the exact fraction a tuneable parameter (unbound_threshold) of the algorithm (Behroozi et al. 2013).

This definition is generally very effective at identifying halos and subhalos—but it fails in two important situations. First, if a subhalo is experiencing significant tidal stripping, the 50% cutoff can remove a subhalo from the catalog that should actually exist. We have found that changing the cutoff can recover the missing subhalos, but the best value of the cutoff is not easily determined. Second, rockstar is occasionally too effective at finding substructure in our high resolution simulations. In particular, it often finds velocity substructures in the cores of our halos that are clearly spurious based on their mass accretion histories and density profiles. Importantly, these two issues do not just affect low mass subhalos, but they can also add or remove halos with ${V}_{\mathrm{max}}$ > 25 km s⁻¹.

Both of these problems can be alleviated by applying an iterative unbinding procedure. We have implemented such an iterative unbinding procedure within rockstar. At each iteration, we remove particles whose kinetic energy exceeds the potential energy from other particles in that iteration. The potential is computed with the rockstar Barnes-Hut method (see Appendix B of Behroozi et al. 2013). We iterate the unbinding until we obtain a self-bound set of particles. Halos are only considered resolved if they contain at least 20 self-bound particles. All halo properties are then computed as usual, but with the self-bound particles instead of the one-pass bound particles. The iterative unbinding recovers the missing subhalos and removes most but not all of the spurious subhalos. Across 13 of our Caterpillar halos, we recover 52 halos with subhalo masses above 10⁸ ${M}_{\odot }$ which would have otherwise been lost using the conventional rockstar. Figure 5 demonstrates how these large haloes can be recovered when our iterative unbinding procedure is used.

Zoom In Zoom Out Reset image size

To remove the remaining spurious subhalos, we also remove halos if ${R}_{\mathrm{max}}$ (i.e., the radius at which the velocity profile reaches its maximum) of the subhalo is larger than the distance between the subhalo and host halo centers. The downside to adding iterative unbinding is that it increases the run time for rockstar by ∼50%. In the rest of this paper, we only consider subhalos with at least 20 self-bound particles passing the ${R}_{\mathrm{max}}$ cut. We define the subhalo mass, ${M}_{\mathrm{sub}}$, as the total gravitational bound mass of a subhalo which is obtained after the complete iterative unbinding procedure has been carried out.

In Table 3, we provide the basic properties of our first 24 re-simulated halos. This includes the simulation name, the halo virial mass, the halo virial radius, concentration, maximum circular velocity, the radius at which the maximum circular velocity occurs, the formation time (defined as when the halo reaches half its present day mass), the redshift of the last major merger (1:3 mass ratio), the fraction of the host mass contained within subhalos, the axis ratios defining the halo shape, and the distance to the closest contamination particle from the host. We adopt a simple naming convention based on when the halos were post processed (1 – ${N}_{\mathrm{halos}}$). Where required, we use a shorthand reference to the resolution of the simulation. These refer to the parameter levelmax inside the IC generation code music (e.g., levelmax = 14 is simply lx14). This means lx14 represents an effective resolution of ${N}_{{\rm{p}}}={({2}^{14})}^{3}$ ($\sim {Aquarius}$ level-2), lx13 $\;\to \;$ ${N}_{{\rm{p}}}={({2}^{13})}^{3}$ ($\sim {Aquarius}$ level-3 or ∼ELVIS resolution), lx12 $\;\to \;$ ${N}_{{\rm{p}}}={({2}^{12})}^{3}$ and lx11 $\;\to \;$ ${N}_{{\rm{p}}}={({2}^{11})}^{3}$. Unless otherwise stated, all halos in the analysis of this paper are the lx14 halos (i.e., our flagship resolution). All halos have similar z = 0 properties except Cat-7 whose properties can be tied to the fact that it has recently undergone a massive merger (1:3 mass ratio at z = 0.03). We obtain extremely large uncontaminated volumes (∼1.4 Mpc) in all but one of our halos (Cat-18 is ∼3% contaminated by mass). The fraction of mass held in subhalos across our sample is ${f}_{m,\mathrm{subs}}=0.121\pm 0.041$ (1σ), though this excludes Cat-7.

Table 3.  The Halo Properties of the First 24 Caterpillar Halos

Halo	${M}_{\mathrm{vir}}$	${R}_{\mathrm{vir}}$	ca	${V}_{\mathrm{max}}$	${R}_{\mathrm{max}}$ b	zformc	zlmmd	${f}_{m,\mathrm{subs}}$ e	c/a	b/a	${R}_{\mathrm{hires}}$ f
Name	($\times {10}^{12}$ ${M}_{\odot }$)	(kpc)		(km s−1)	(kpc)						(Mpc)
Cat-1	1.559	306.378	7.492	169.756	34.083	0.894	2.157	0.207	0.841	0.869	0.998
Cat-2	1.791	320.907	8.374	178.851	55.268	0.742	0.731	0.148	0.636	0.719	1.463
Cat-3	1.354	292.300	10.170	172.440	31.701	0.802	0.802	0.136	0.865	0.927	1.894
Cat-4	1.424	297.295	8.573	164.344	53.466	0.936	0.922	0.175	0.671	0.739	1.531
Cat-5	1.309	289.079	12.108	176.399	32.103	0.564	0.510	0.069	0.552	0.815	1.608
Cat-6	1.363	292.946	10.196	171.647	33.632	1.161	1.295	0.153	0.508	0.528	1.295
Cat-7	1.092	272.099	1.757	134.148	157.438	0.070	0.032	0.735	0.151	0.207	1.477
Cat-8	1.702	315.466	13.507	198.564	40.819	1.516	2.235	0.078	0.605	0.787	1.540
Cat-9	1.322	289.987	12.401	177.414	30.336	1.255	1.236	0.094	0.513	0.762	2.080
Cat-10	1.323	290.119	11.714	174.989	39.721	1.644	2.010	0.103	0.559	0.703	1.775
Cat-11	1.179	279.187	12.522	172.723	53.187	1.059	4.368	0.215	0.597	0.867	1.135
Cat-12	1.763	319.209	11.402	191.259	52.717	1.336	9.616	0.073	0.584	0.645	1.162
Cat-13	1.164	277.938	12.850	171.222	33.757	1.161	11.092	0.090	0.578	0.645	1.566
Cat-14	0.750	240.119	9.135	137.437	26.660	1.144	4.258	0.113	0.705	0.859	2.178
Cat-15	1.505	302.787	8.983	174.124	37.043	1.144	3.165	0.126	0.849	0.877	1.119
Cat-16	0.982	262.608	11.737	155.362	28.768	1.315	3.165	0.106	0.618	0.792	0.671
Cat-17	1.319	289.800	12.765	179.056	38.329	1.846	1.976	0.093	0.664	0.881	1.299
Cat-18	1.407	296.099	7.887	163.920	57.217	0.493	0.435	0.159	0.676	0.816	0.397
Cat-19	1.174	278.770	10.468	164.726	29.112	1.541	2.118	0.169	0.664	0.937	1.712
Cat-20	0.763	241.484	13.324	149.672	30.417	1.492	5.427	0.099	0.601	0.733	1.311
Cat-21	1.881	326.206	10.618	190.683	50.954	1.126	1.198	0.118	0.482	0.611	1.453
Cat-22	1.495	302.116	10.666	180.647	35.860	0.841	29.488	0.080	0.512	0.694	1.744
Cat-23	1.607	309.524	12.489	190.705	32.421	1.161	9.616	0.094	0.607	0.784	1.207
Cat-24	1.334	290.866	11.378	176.911	36.800	1.144	3.608	0.090	0.689	0.734	1.102
Meang	1.368	291.791	10.903	173.167	38.886	1.144	4.410	0.121	0.634	0.771	1.402
±1σ	0.285	21.610	1.761	13.441	9.530	0.329	6.112	0.041	0.103	0.102	0.409

Notes.

^aConcentration defined by ratio of the virial radius and the scale radius; ${R}_{\mathrm{vir}}/{R}_{{\rm{s}}}$. ^bThe radius at which the ${V}_{\mathrm{max}}$ occurs. ^cRedshift of host formation defined as when the host main branch progenitor mass equals 0.5 ${M}_{\mathrm{vir}}(z=0)$. ^dRedshift of last major merger defined as a halo with 1/3 mass merging into the main branch of the host. ^eFraction of the host mass in subhalos. ^fDistance to the closest contamination particle from the host. ^gMeans and deviations were calculated over all halos except Cat-7 as it has undergone a very recent major merger.

Download table as:  ASCIITypeset image

In Figure 6 we plot the concentration–mass (c–M) relation of the parent simulation for similarly sized halos ($11.5\lt {\mathrm{log}}_{10}\ {M}_{\odot }\lt 12.5$, gray band indicating the 1σ dispersion) and overlay the concentration (${R}_{\mathrm{vir}}$/${R}_{{\rm{s}}}$) and host mass of the high resolution halos. This shows that for nearly all of the halos, we are sampling within 68% of the average c–M relation at a fixed halo mass. Again, Cat-7 is an outlier with an extremely low concentration because it recently underwent a major merger and has an extremely large substructure mass fraction, so its concentration is not meaningful. For this reason we do not include it in the quantitative analysis in terms of determining average halo profile shapes or the mass function slopes. Its properties are still shown and plotted in the various tables and figures, however. Recently, Buck et al. (2015) found that the thickness of planes of satellites depends on the concentration of the host halo. Specifically, they found the thinnest planes are only found in the most concentrated, and hence earliest formed halos. The fact that we sample relatively average concentrations for halos of this mass range means that it is less likely that these hosts will contain planes of satellites, or if they do, their thicknesses will be quite large (B. Griffin et al. 2016a, in preparation). As the Caterpillar sample grows, we will eventually sample overly concentrated halos, enabling us to see in better detail how this concentration-plane relation holds.

Zoom In Zoom Out Reset image size

In Figures 7 and 8 we show images of the dark matter distribution in each of our 24 high-resolution halos at redshift z = 0. The brightness of each pixel is proportional to the logarithm of the dark matter density squared (i.e., log(${\rho }^{2}$)projected along the line of sight. To enhance the density contrast, each panel has a different maximum density. We note that a similarly colored pixel in one panel does not necessarily mean the density is the same for another panel. The panel width is 1 Mpc and the local dark matter density of the particles in each pixel is estimated with an SPH kernel interpolation scheme based on the 64 nearest neighbor particles. Upon first inspection it is clear that each halo is littered with an abundance of dark matter substructures of varied shapes and sizes. In some cases, there are reasonably large neighbors (e.g., Cat-4, 7, 11, 24). By virtue of our selection criteria these neighbors are no larger than 0.5 × ${M}_{\mathrm{vir}}$ of the central host. In any case, in under a Gyr, these SMC/LMC sized systems (${M}_{\mathrm{peak}}\gt {10}^{11}$ ${M}_{\odot }$) will likely undergo a major merger with the host galaxy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 9 we show the mass evolution of each of the halos. As highlighted by the inset which shows the normalized mass evolution, there is a wide variety of formation histories. In our initial catalog of 24 halos, 6 halos (Cat-2, 3, 4, 5, 7, 18) have had major mergers since z = 1. The halos going above a normalized mass ratio of 1.0 have had a halo pass through them relatively recently, which momentarily gives them extra mass such that it is larger than their z = 0 mass (e.g., Cat-2). This indicates that many of the halos are yet to reach an equilibrium state.

Zoom In Zoom Out Reset image size

Adopting the same criteria as Neto et al. (2007) we assess whether the host halos are relaxed. If their substructure mass fraction is below 0.1, their normalized offset between the center of mass of the halo (i.e., computed using all particles within ${R}_{\mathrm{vir}}$) and the potential center (${x}_{\mathrm{off}}\;\equiv | {r}_{c}-{r}_{\mathrm{cm}}| /{R}_{\mathrm{vir}}$, ${r}_{c}\;\equiv$ center of the potential well, ${r}_{\mathrm{cm}}\ \equiv$ center of mass) is below 0.07 and their virial ratio ($2\;T/| U|$) is below 1.35, then the host is considered relaxed. In Table 4 we provide the relaxed state of the halo. Many of the halos are in fact unrelaxed under this definition which is by design—we are sampling a wide range of assembly histories and so halos with recent merger events that prevent the halos from being fully virialized naturally make up part of our sample.

Table 4.  The Relaxed Nature and Einasto Profile Parameters for the First 24 Caterpillar Halos

Halo	Relaxeda	${\rho }_{-2}$ b	${r}_{-2}$ c	αd	Qmin	Qmin
Name			(×105)		Ein.	NFW
Cat-1	⨯	9.929	27.182	0.128	0.039	0.103
Cat-2	✓	3.531	46.846	0.185	0.040	0.033
Cat-3	✓	11.604	25.309	0.151	0.028	0.067
Cat-4	⨯	5.504	34.962	0.154	0.039	0.080
Cat-5	✓	13.305	24.152	0.167	0.018	0.050
Cat-6	⨯	9.254	28.299	0.164	0.027	0.058
Cat-7	⨯	0.482	90.528	0.075	0.082	0.162
Cat-8	✓	8.027	34.555	0.236	0.033	0.036
Cat-9	✓	12.159	25.543	0.186	0.011	0.034
Cat-10	✓	14.482	23.492	0.168	0.029	0.049
Cat-11	⨯	10.617	25.450	0.173	0.049	0.071
Cat-12	✓	8.897	31.157	0.160	0.030	0.062
Cat-13	✓	11.842	25.010	0.187	0.009	0.036
Cat-14	✓	10.455	21.347	0.139	0.018	0.078
Cat-15	✓	8.905	28.744	0.139	0.025	0.078
Cat-16	✓	10.445	24.162	0.166	0.026	0.056
Cat-17	✓	11.283	26.799	0.195	0.019	0.026
Cat-18	⨯	6.998	31.012	0.150	0.023	0.082
Cat-19	⨯	9.079	26.881	0.164	0.032	0.054
Cat-20	✓	11.415	22.158	0.199	0.018	0.020
Cat-21	✓	6.682	36.486	0.175	0.025	0.043
Cat-22	✓	10.857	26.957	0.156	0.017	0.063
Cat-23	✓	13.486	26.024	0.172	0.018	0.039
Cat-24	✓	7.040	31.734	0.181	0.044	0.050
Meane	⋯	9.817	28.446	0.169	0.027	0.055
±1σ	⋯	2.6106	5.5677	0.023	0.010	0.020

Notes. For comparison,  Q_min  fits for NFW halo profiles are also listed.

^aRelaxed criteria are based on those of Neto et al. (2007). If the substructure mass fraction is below 0.1, their center of mass displacement (${x}_{\mathrm{off}}=| {r}_{c}-{r}_{\mathrm{cm}}| /{R}_{\mathrm{vir}}$) is below 0.07 and their virial ratio ($2\;T/| U|$) is below 1.35, then the host is considered relaxed. ^bThe density at the characteristic radius in units of: ${10}^{10}{h}^{2}$ ${M}_{\odot }$ kpc⁻³. ^cThe characteristic radius or "peak" radius of the ${r}^{2}\rho$ profile in units of ${h}^{-1}$ kpc. ^dEinasto slope parameter of the form $\rho (r)\propto \mathrm{exp}({-\mathrm{Ar}}^{\alpha })$. ^eMeans and deviations were calculated over all halos except Cat-7, as it has undergone a very recent major merger.

Download table as:  ASCIITypeset image

In Figure 10 we plot the spherically averaged halo profiles for each of our 24 simulated halos. We draw the measured density profile as a thick line of a given color and continue the fit beyond the smallest radius possible set by Power et al. (2003) as a vertical black dashed line. We truncate each fit at this radius. There is a clear diversity in the profile shapes owing in part to the assembly histories of each halo. Halos which have undergone a recent major merger whose substructure mass fractions are higher than average are primarily dominated by a single subhalo (e.g., Cat-7 has a large subhalo at 200 kpc). The fitting formula we have used to describe the mass profile of our simulated halos follows the method of Navarro et al. (2010) and is given by the following Einasto form (over all particles within the virial radius):

$$\begin{eqnarray}&&\mathrm{log}[\rho (r)/{\rho }_{-2}]=(-2/\alpha )[{(r/{r}_{-2})}^{\alpha }-1].\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

The ${r}_{-2}$ is the scale length of the halo which can be obtained without resorting to a particular fitting formula. We compute the logarithm of the slope profile and identify where a low-order polynomial fit to it intersects the isothermal value ($\gamma =2$). Unlike the Navarro–Frenk–White profile (NFW) the peak parameter in the Einasto profile, α, is allowed to vary and thus provides a third parameter for the fitting formula. The best fitting parameters are found by minimizing the deviation between model and simulation at each bin. Specifically we minimize the function Q², defined as

$$\begin{eqnarray}&&{Q}^{2}=\displaystyle \frac{1}{{N}_{\mathrm{bins}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{bins}}}{(\mathrm{ln}({\rho }_{i})-\mathrm{ln}({\rho }_{i}^{\mathrm{model}}))}^{2}.\end{eqnarray} \tag{ 2 }$$

In this manner we find a function which clearly illustrates the deviation of the simulated and model profiles. In Table 4 we show our minimum Q² parameter (Q_min), characteristic scale radius ${r}_{-2}$, and their corresponding densities for each halo. For our relaxed halos, Q_min for our Einasto fits are 0.027 ± 0.010 indicating reasonable agreement between the simulated and model Einasto profiles. This is better than our NFW profile fits for which we obtain ${Q}_{\mathrm{min}}=0.055\pm 0.020$. The peak parameters for our Einasto fits are 0.169 ± 0.023, which is comparable to those of the Aquarius halos (α = 0.145–0.173) studied in Navarro et al. (2010). For halos which significantly deviate from the mean, it is important to note that those halos are not relaxed and so by definition will not provide meaningful Einasto/NFW fits. A more detailed study of the halo density profiles is reserved for future work.

In Figure 11(a) we show the cumulative abundance of subhalos as a function of their maximum circular velocity for each Caterpillar halo. Since we achieve excellent convergence (see the Appendix), we reliably resolve halos with circular velocities of ∼4 km s⁻¹, which is crucial for identifying the sites of first star formation. At the high ${V}_{\mathrm{max}}$ end we find a variety of different sized subhalos for each host. Some hosts have only one 10 km s⁻¹ subhalo, whereas another host halo has a large 70 km s⁻¹ subhalo within the virial radius. Between 5–20 km s⁻¹ all halos are very similar in their ${V}_{\mathrm{max}}$ function slopes within a slight offset owing to normalization stemming from the differences in host halo mass. At low ${V}_{\mathrm{max}}$ values (∼3 km s⁻¹) we begin to lose completeness of our host subhalo sample due to lack of resolution. We additionally include in this figure the ${V}_{\mathrm{max}}$ function for subhalos at infall (i.e., when a subhalo first crosses the virial radius of the host). Since dynamical friction affects the highest mass subhalos the fastest, the biggest difference in the functions occurs at the high-mass end whereby several LMC sized systems (${M}_{\mathrm{peak}}$ > ${10}^{11}\ {M}_{\odot }$) have been destroyed (over a timescale of 1–2 Gyr) between infall and z = 0. These large LMC sized-systems at infall can host anywhere from 4%–30% of the Milky Way sized halo's subhalos at z = 0 depending on their orbit and infall time (B. Griffen et al. 2016b, in preparation). In solid black we also plot the Aquarius Aq-A2 halo from Springel et al. (2008) (using the same version of rockstar that we used for the Caterpillar halos). We find the differences in the cosmology (${\sigma }_{8}=0.9$) and the slightly higher resolution of Aquarius leads to systematic differences in subhalo abundance.

Zoom In Zoom Out Reset image size

In Figure 11(b), we show the subhalo mass functions for each of the halos. Our results are best fit by the power law dN/d$M\propto {M}^{-1.88\pm 0.10}$, which is less steep than that found in the Aquarius halos of Springel et al. (2008). This slope is the best fit over the ranges ${10}^{5}\mbox{--}{10}^{8}\;$ ${M}_{\odot }$. We do observe a scatter in the subhalo abundances. This can be explained by the subtle concentration-subhalo-abundance relation whereby for fixed halo mass, there are more (less) subhalos belonging to hosts which are less (more) concentrated (e.g., Zentner et al. 2005; Watson et al. 2011; Mao et al. 2015). Indeed, we find halos that are more concentrated (see Figure 6) have lower normalizations than those less concentrated at fixed ${V}_{\mathrm{max}}$. This is simply because halos that are more concentrated have formed earlier and so subhalos have spent substantially longer undergoing dynamical disruption within the host compared to similar sized subhalos orbiting less concentrated hosts.

In Figure 11(c) we show the subhalo radial mass fraction which indicates high variability in the contribution to the total halo mass from substructure as a function of galactocentric distance. For example, at 0.1${R}_{\mathrm{vir}}$, the total mass contributing to the host halo mass from substructure varies by a factor of 10 or more when normalized by mass. At ${R}_{\mathrm{vir}}$ our substructure mass fraction varies by $\sim 10\%$ (see Table 3 for exact fractions). Cat-7 has a large component of the halo mass in substructure at low radii because it has recently undergone a major merger (z = 0.03). Those halos with a large substructure mass fraction generally have had a recent major merger and are in the process of disrupting the recently accreted systems. On average, for a fixed fraction of the virial radius, the Caterpillar halos have less mass in substructure than that found in the Aq-A Aquarius halo (see solid black line, calculated using the exact same code). In Figure 11(d) we plot the normalized number of subhalos as a function of radius scaled by the virial radius of the host. We find the scatter in the number of subhalos as a function of galactocentric distance is a factor of 3 across all halos except within the inner 10% of the host halo where we are subject to noise in the halo finding produced by rockstar. Again, in solid black we also plot the Aquarius Aq-A2 halo from Springel et al. (2008). We find the differences in the cosmology (${\sigma }_{8}=0.9$) and in particular the slightly higher resolution of Aquarius lead to this systematic difference in the subhalo number density.

We also examine halos which are massive enough to form stars but have no luminous counterpart in the nearby universe (i.e., the too big to fail problem, hereafter TBTF. Boylan-Kolchin et al. 2011). To do this we select halos with ${V}_{\mathrm{peak}}$ > 30 km s⁻¹, which are subhalos large enough to retain substantial gas in the presence of an ionizing background and therefore theoretically should form stars. We follow the same definition as in Garrison-Kimmel et al. (2014a) to count two classes of halos. Strong massive failures (SMFs) are too dense to host any of the currently known bright MW classical dwarf spherioidals (dSph) galaxies. MFs include all SMFs as well as all massive subhalos which have densities consistent with the high-density dSphs (i.e., Draco and Ursa Minor) but cannot be associated with them without allowing a single dwarf galaxy to be hosted by multiple halos (i.e., assuming every observable dSph galaxy is hosted by exactly one halo). Most subhalos in the range of ${V}_{\mathrm{max}}$ = 25–30 km s⁻¹ could host a low-density dwarf and as such are not defined as MFs.

In Figure 12, we plot a sample of the rotation curves for three different Caterpillar halos (Cat-19, Cat-13 and Cat-18). We adopt the Boylan-Kolchin et al. (2012) Einasto correction for $R\lt 291$ pc, which differ from the ELVIS profile fits in that they extrapolate their entire profile from ${R}_{\mathrm{max}}$ with various density profile shapes. Black squares depict circular velocities of the classical dwarf galaxies (with luminosities above 2 × 10⁵${L}_{\odot ,V}$), as measured by Wolf et al. (2010). Dashed–dotted cyan lines are LMC analogs (i.e., ${V}_{\mathrm{max}}$ > 60 km s⁻¹, which are excluded from our failure analysis), and blue dashed lines are MFs and red solid lines are SMFs. Thin black solid lines are subhalos which pass the test of having at least one observed dwarf with a comparable circular velocity, i.e., a circular profile goes through one of the observed dwarf galaxy data points. The cumulative numbers of profiles above and below the observed classical dwarfs for Cat-19 are 10 MFs, 5 SMFs and 1 LMC analogue. Similarly we find Cat-13 has 11 MFs and 8 SMFs. Cat-18 has the most failures of any halo with 21 MFs and 14 SMFs with 1 LMC analogue.

Zoom In Zoom Out Reset image size

In Figure 13 we plot the number of strong and massive failures across all Caterpillar halos. Specifically we plot the fraction of hosts with fewer than N MFs and N SMFs within 300 kpc of each host as a function of N (black lines). Averaging over the entire Caterpillar sample (excluding Cat-7 as it has recently had a massive major merger), we predict 8 ± 3 (1σ) SMFs and 16 ± 5 (1σ) MFs within 300 kpc. If the Milky Way was well described by such an average we would expect to have these failures.

Zoom In Zoom Out Reset image size

For comparison, we plot the ELVIS MF and SMF counts (red lines). The lower resolution in these simulations requires an extrapolation of the velocity profile from ${V}_{\mathrm{max}}$ and ${R}_{\mathrm{max}}$ using an analytic Einasto profile. Qualitatively, both simulation suites agree that there are a significant number of both MFs and SMFs. Quantitatively, there are several differences, which we now describe. The Caterpillar suite has many more MFs. This is due to our ability to better resolve high ${V}_{\mathrm{peak}}$ subhalos that have been tidally stripped. In particular, the iterative unbinding procedure described in Section 2.5 removes the need for the rockstar unbound_threshold parameter (Behroozi et al. 2013). We can simulate the effect of the standard unbound_threshold = 0.5 cut by removing halos whose bound masses are less than 50% of their mass prior to unbinding. The MF counts with this cut are shown in Figure 13 by the blue lines, which are very similar to the ELVIS MF counts.

Zoom In Zoom Out Reset image size

The Caterpillar suite also has significantly fewer SMFs compared to ELVIS. This discrepancy is likely due to the fact that we have measured rather than extrapolated the subhalo density profiles. Variations in the Einasto shape parameter (α) greatly affect the MF count (see Figure 4 of Garrison-Kimmel et al. 2014a). The Einasto fits to our density profiles have α typically closer to 0.2, which is less discrepant.

While TBTF is a prevalent problem in pure N-body simulations, many authors have indicated that the tension between the circular velocities of observed classical dwarfs and simulated subhalos can be alleviated with the addition of supernovae feedback and ram pressure stripping (e.g., Pontzen & Governato 2012; Zolotov et al. 2012; Brooks et al. 2013; Gritschneder & Lin 2013; Arraki et al. 2014; Del Popolo et al. 2014; Elbert et al. 2015; Maxwell et al. 2015) or by making dark matter self-interacting (e.g., Vogelsberger et al. 2012; Zavala et al. 2013). Our results are within 1σ of the number of failures found by Garrison-Kimmel et al. (2014a) (i.e., 12 MFs within 300 kpc), even when using a better density profile estimation.

Table 5.  Halo Properties for the First Set of 12 Caterpillar Halos at Each Resolution

Name	Geometry	n${R}_{\mathrm{vir}}$	lx	${M}_{\mathrm{vir}}$	${R}_{\mathrm{vir}}$	c	${V}_{\mathrm{max}}$	${R}_{\mathrm{max}}$	${z}_{0.5}$	zlmm	${f}_{m,\mathrm{subs}}$	b/c	c/a	${R}_{\mathrm{res}}$
				(1012${M}_{\odot }$)	(kpc)		(km s−1)	(kpc)						(Mpc)
Cat-1	EA	5	11	1.579	307.690	7.762	172.293	35.049	0.881	2.118	0.161	0.810	0.828	1.391
			12	1.560	306.491	7.494	171.060	34.292	0.894	2.118	0.171	0.824	0.869	1.248
			13	1.560	306.458	7.647	170.707	36.451	0.894	2.157	0.197	0.842	0.883	1.138
			14	1.559	306.378	7.492	169.756	34.083	0.894	2.157	0.207	0.841	0.869	0.998
Cat-2	EB	4	11	1.807	321.876	7.621	176.924	64.931	0.742	0.719	0.092	0.596	0.724	1.577
			12	1.782	320.357	8.575	179.069	53.856	0.742	0.719	0.112	0.607	0.716	1.522
			13	1.792	320.970	8.382	178.753	54.360	0.742	0.731	0.137	0.643	0.731	1.480
			14	1.791	320.907	8.374	178.851	55.268	0.742	0.731	0.148	0.636	0.719	1.463
Cat-3	EB	4	11	1.343	291.538	10.763	175.066	26.554	0.802	0.790	0.079	0.961	0.971	1.966
			12	1.355	292.387	10.489	175.142	29.083	0.802	0.790	0.100	0.868	0.915	1.926
			13	1.355	292.400	10.523	172.946	31.565	0.802	0.802	0.117	0.850	0.905	1.906
			14	1.354	292.300	10.170	172.440	31.701	0.802	0.802	0.136	0.865	0.927	1.894
Cat-4	EB	4	11	1.503	302.676	8.308	169.309	91.132	0.894	0.908	0.128	0.749	0.825	1.791
			12	1.415	296.632	9.208	168.170	85.989	0.922	0.922	0.120	0.681	0.762	1.594
			13	1.434	298.009	8.434	164.999	59.225	0.936	0.922	0.156	0.673	0.743	1.561
			14	1.424	297.295	8.573	164.344	53.466	0.936	0.922	0.175	0.671	0.739	1.531
Cat-5	EB	4	11	1.306	288.846	11.897	173.913	33.844	0.584	0.519	0.025	0.551	0.835	1.676
			12	1.318	289.714	11.896	174.223	36.356	0.574	0.519	0.041	0.547	0.765	1.657
			13	1.314	289.450	12.324	176.818	29.346	0.574	0.519	0.055	0.556	0.825	1.617
			14	1.309	289.079	12.108	176.399	32.103	0.564	0.510	0.069	0.552	0.815	1.608
Cat-6	EB	4	11	1.371	293.516	10.373	172.100	32.944	1.144	1.275	0.094	0.495	0.534	1.708
			12	1.347	291.848	10.522	172.873	30.622	1.161	1.275	0.116	0.496	0.525	1.495
			13	1.366	293.186	10.086	170.858	32.794	1.161	1.295	0.138	0.510	0.529	1.300
			14	1.363	292.946	10.196	171.647	33.632	1.161	1.295	0.153	0.508	0.528	1.295
Cat-7	EB	4	11	1.142	276.168	2.513	139.055	140.859	0.065	0.057	0.693	0.191	0.301	1.756
			12	1.111	273.686	2.487	136.803	145.085	0.074	0.057	0.615	0.170	0.288	1.520
			13	1.091	272.009	1.674	133.574	162.291	0.065	0.036	0.693	0.168	0.235	1.510
			14	1.092	272.099	1.757	134.148	157.438	0.070	0.032	0.735	0.151	0.207	1.477
Cat-8	EB	4	11	1.729	317.150	13.081	198.577	46.800	1.541	2.195	0.032	0.602	0.768	1.690
			12	1.716	316.337	13.154	198.229	39.671	1.315	2.195	0.053	0.594	0.775	1.597
			13	1.701	315.450	13.340	197.637	39.810	1.516	2.235	0.066	0.599	0.791	1.550
			14	1.702	315.466	13.507	198.564	40.819	1.516	2.235	0.078	0.605	0.787	1.540
Cat-9	EB	4	11	1.330	290.616	12.568	177.522	32.309	1.236	1.217	0.050	0.493	0.762	2.383
			12	1.331	290.654	11.616	175.047	27.903	1.236	1.236	0.070	0.486	0.754	2.101
			13	1.329	290.538	12.132	176.808	30.297	1.255	1.236	0.085	0.500	0.754	1.833
			14	1.322	289.987	12.401	177.414	30.336	1.255	1.236	0.094	0.513	0.762	2.080
Cat-10	EB	4	11	1.319	289.809	11.902	175.553	41.894	1.699	2.010	0.052	0.561	0.709	1.983
			12	1.332	290.764	11.439	174.479	29.806	1.516	2.010	0.069	0.551	0.679	1.870
			13	1.328	290.477	11.714	175.124	25.839	1.644	2.010	0.088	0.559	0.703	1.740
			14	1.323	290.119	11.714	174.989	39.721	1.644	2.010	0.103	0.559	0.703	1.775
Cat-11	EB	4	11	1.194	280.361	10.551	165.980	62.881	1.059	4.368	0.175	0.527	0.719	1.490
			12	1.196	280.471	10.044	163.290	70.202	1.043	4.368	0.200	0.525	0.703	1.408
			13	1.190	280.043	12.272	173.893	45.727	1.059	1.644	0.199	0.590	0.868	1.192
			14	1.179	279.187	12.522	172.723	53.187	1.059	4.368	0.215	0.597	0.867	1.135
Cat-12	EA	5	11	1.786	320.627	11.723	191.564	59.256	1.336	2.542	0.034	0.592	0.724	1.664
			12	1.749	318.388	11.824	192.085	56.859	1.336	2.542	0.042	0.572	0.686	1.342
			13	1.767	319.441	11.663	191.320	49.435	1.336	2.542	0.062	0.571	0.703	1.239
			14	1.763	319.209	11.402	191.259	52.717	1.336	9.616	0.073	0.584	0.645	1.162

Note. The resolution details for each refinement level (i.e., 11, 12, 13, 14) can be found in Table 2 and the geometry definitions can be found in Table 1.

Download table as:  ASCIITypeset image