Creating mock catalogues of stellar haloes from cosmological simulations
Ben Lowing, Wenting Wang, Andrew Cooper, Rachel Kennedy, John Helly, Carlos Frenk, Shaun Cole
MMon. Not. R. Astron. Soc. , 1–19 (2013) Printed 15 October 2018 (MN L A TEX style file v2.2)
Creating mock catalogues of stellar haloes fromcosmological simulations
Ben Lowing (cid:63) , Wenting Wang † , Andrew Cooper , , Rachel Kennedy , John Helly ,Shaun Cole and Carlos Frenk Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham, DH1 3LE, UK National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang, Beijing 10012, P.R. China
Accepted ... Received ... in original form ...
ABSTRACT
We present a new technique for creating mock catalogues of the individual stars thatmake up the accreted component of stellar haloes in cosmological simulations andshow how the catalogues can be used to test and interpret observational data. Thecatalogues are constructed from a combination of methods. A semi-analytic galaxyformation model is used to calculate the star formation history in haloes in an N -body simulation and dark matter particles are tagged with this stellar mass. Thetags are converted into individual stars using a stellar population synthesis model toobtain the number density and evolutionary stage of the stars, together with a phase-space sampling method that distributes the stars while ensuring that the phase-spacestructure of the original N -body simulation is maintained. A set of catalogues basedon the ΛCDM Aquarius simulations of Milky Way mass haloes have been created andmade publicly available on a website. Two example applications are discussed thatdemonstrate the power and flexibility of the mock catalogues. We show how the richstellar substructure that survives in the stellar halo precludes a simple measurementof its density profile and demonstrate explicitly how pencil-beam surveys can returnalmost any value for the slope of the profile. We also show that localized variations inthe abundance of particular types of stars, a signature of differences in the compositionof stellar populations, allow streams to be easily identified. Key words: methods: numerical, galaxies: haloes
Extensive diffuse stellar haloes are known to surround boththe Milky Way and M31 (e.g. Deason et al. 2011; Ibata et al.2014). The ΛCDM model of hierarchical structure formationpredicts that these haloes should be a ubiquitous outcomeof galaxy formation resulting from the accretion of satellitegalaxies that fall into the potential of the larger host galaxyand are subsequently disrupted (Bullock & Johnston 2005;Cooper et al. 2010). The dynamical time of the stellar halois very long compared to the age of the host galaxy and so itpreserves a memory of its initial conditions. From the identi-fication of distinct populations within haloes, and the studyof properties such as kinematics or chemical composition, avast amount of information about the assembly history ofthe galaxy can be inferred.The outer regions of stellar haloes are dominated bymultiple, extensive tidal streams, while the inner regions, (cid:63)
E-mail: [email protected] † E-mail:[email protected] where dynamical timescales are shorter, experience a morecomplete destruction of substructure and therefore have asmoother stellar distribution. In addition to the accretedcomponent, stellar haloes are thought to contain a secondcomponent made up of stars that formed in situ . Recentobservational claims of a two-component halo in the MilkyWay have focused attention on this possibility (e.g. Carolloet al. 2010; Beers et al. 2012; Kafle et al. 2013), while theevidence for two components is disputed by Sch¨onrich et al.2014. From a theoretical perspective there is a wide vari-ety of possible formation mechanisms for in situ halo stars,many of which overlap with possible origins of the Galacticthick disc. These include the scattering of stars from the thindisc, the collapse of the disc at high redshift, star formationin gaseous tidal streams and thermal instabilities in a hydro-static gas halo (e.g. Abadi et al. 2006; Zolotov et al. 2010;Font et al. 2011a; Tissera et al. 2013, Parry et al. in prep.).Most of these mechanisms result in in situ components thatdominate in the inner few kiloparsecs of the galaxy and re-main important out to ∼ c (cid:13) a r X i v : . [ a s t r o - ph . GA ] O c t Ben Lowing et al
Over the last few years there has been much effort ded-icated to modelling stellar haloes (e.g. Bullock & Johnston2005; Cooper et al. 2010; Zolotov et al. 2010; Font et al.2011a; Tissera et al. 2012) motivated by improvement inobservations. Recent surveys such as sdss segue (Yanny etal. 2009) and p an- starrs (Kaiser et al. 2010) have startedto probe deeper into the Milky Way halo and are starting tobuild up a picture of its complex structure. The gaia mis-sion will lead to further major improvements in the MilkyWay data (Jordi et al. 2010; de Bruijne 2012) and new pho-tometric surveys are extending these studies to other nearbygalaxies (McConnachie et al. 2009; Mart´ınez-Delgado et al.2010). We are now at the point where careful comparisonsbetween the haloes predicted by theory and those observedin the universe are becoming possible.In order to make a direct comparison of theoretical pre-dictions with observations of stellar haloes, a model of stellarhalo formation is required, along with a way to convert theoutput into observable quantities. Furthermore, to obtainrealistic accretion histories, a complete model of galaxy for-mation in a cosmological context is essential. Various studieshave focused on supplementing N -body simulations with an-alytical modelling of star formation and chemical enrichment(e.g. Bullock et al. 2001; Robertson et al. 2005; Font et al.2006; Johnston et al. 2008; De Lucia Helmi 2008; Tumlin-son 2010). Alternatively, cosmological hydrodynamical sim-ulations have also been used to study the overall structureand assembly of stellar haloes (e.g. Crain et al. 2009; Zolo-tov et al. 2010; Libeskind et al. 2011; Tissera et al. 2012;Pillepich et al. 2014), although these do not yet have the ca-pability to resolve the detailed kinematic and spatial struc-ture of individual halo stellar populations. Even “zoom” sim-ulations of individual galaxies represent the stellar halo witha relatively small number of star particles, since this com-ponent amounts to only a few per cent of the total stellarmass.In this paper we employ the particle tagging methoddeveloped by Cooper et al. (2010, hereafter C10) to con-struct stellar halo models. The method is based on high-resolution dark matter simulations in which the star forma-tion occurring in haloes is calculated using a semi-analyticgalaxy formation model. At every output time in the simu-lation, an appropriate set of dark matter particles is taggedwith the stellar mass that formed since the previous outputtime. The advantages of this technique are that it can cre-ate models of much higher resolution and at a much lowercomputational cost than full hydrodynamical simulations.C10 demonstrated that the method produces complex stellarhaloes containing tidal streams, shells and other substruc-ture, whose ages and metallicities are in broad agreementwith recent observations.A limitation of the C10 tagging method is that “starparticles” are not individual stars but rather represent en-tire stellar populations which can contain thousands of solarmasses of stars. The challenge is to compare these sparselysampled haloes to actual observational surveys which oftenfocus on a particular type of star, such as blue horizontalbranch (BHB) stars. These may not trace the overall stellarmass distribution in an unbiased way. To overcome this limi-tation, we use a method similar to that employed in galaxia (Sharma et al. 2011) to construct customizable mock cata-logues of halo stars. Each star particle is treated as a sep- arate simple stellar population using theoretical isochronesto determine the abundance of each type of star expectedin a population of the corresponding age and metallicity. Tocreate positions and velocities of individual stars, the starparticle distribution is oversampled using a phase-space ker-nel method based on the EnBiD code (Sharma & Steinmetz2006) which maintains the underlying phase-space structure.Our mock catalogues allow us to compare the structureof simulated stellar haloes to observational data in a more re-alistic way. For example, obscuration by the Milky Way disc,survey limits imposed by the high cost of deep photometryand spectroscopy and the low surface density of halo tracersmean that the structure of the Galactic stellar halo is onlywell sampled to large distances in small patches of the sky.Nevertheless, attempts are often made to infer the profile ofthe Milky Way stellar halo from these limited-area obser-vations (e.g. Watkins et al. 2009; Sesar et al. 2011; Akhteret al. 2012). In this paper we show how the complicated un-even structure of the outer accreted stellar halo results inlarge fluctuations among measurements of its global prop-erties based on observations of small sky patches. We alsoshow how different satellites that build up the halo createobservable variations in the stellar populations across thesky, for example in the ratio of stellar types. This variationcan potentially help to detect new halo substructures in theMilky Way.In Section 2 we outline the methods used to create stel-lar halo models. The final output of the process, a set ofmock catalogues of individual stars and their properties, isdescribed in further detail in Section 3. Sections 4 and 5present two simple example applications of the mock cata-logues to study the overall structure of stellar haloes in oursimulations and the distribution of different stellar typesover the sky. Finally, in Section 6 we discuss some potentialadvanced applications of the mock catalogues and how theycan be used to interpret and test observational strategies. Creating a model of stellar haloes based on a consistent the-ory of galaxy formation in a cosmological context requiresa combination of various techniques. The foundation of ourmodels is a set of N -body dark matter only “zoom” cos-mological simulations of Milky Way mass haloes onto whichwe have grafted the Durham semi-analytic model of galaxyformation, galform (Cole et al. 2000). galform predictsthe amount of stars that form in every halo and subhaloat each output time, along with an estimate of the metal-licity of each stellar population. This stellar mass is thenassigned to dark particles within the N -body simulationsusing the C10 particle tagging technique, which allows usto track the accretion and disruption of satellites and followthe fate of the stripped stars. Finally, to convert the massivesimulation particles into mock catalogues (mocks for short)of individual stars, we use a method based on stellar pop-ulation modelling and phase-space sampling. Each of thesesteps is briefly outlined in the following subsections. c (cid:13) , 1–19 aking Mocks of the Stellar Halo N -body Simulations For this work we use the haloes from the Aquarius project,which simulates six dark matter haloes of mass ∼ M (cid:12) at multiple resolution levels (Springel et al. 2008a,b; Navarroet al. 2010). The simulations assume the standard ΛCDMcosmology with parameters chosen to be consistent with theresults from the WMAP M = 0 .
25; cosmologi-cal constant, Ω Λ = 0 .
75; power spectrum normalization, σ = 0 .
9; spectral index, n s = 1; and Hubble parameter h = 0 .
73. The six haloes were selected randomly from aset of isolated ∼ M (cid:12) haloes identified in a lower resolu-tion (900 -particle) parent simulation of a 100 h − Mpc cubicvolume (Gao et al. 2008). The isolation criterion requires ahalo to have no neighbours exceeding half its mass within1 h − Mpc. This weak selection criterion ensures that thehaloes are not members of any massive groups or clusters.The Aquarius simulations are labeled Aq-A to Aq-E;in each case we use the second highest resolution (“level-2”in the Aquarius notation) simulation, with a particle massof at most 10 h − M (cid:12) . We have used only five of the sixAquarius haloes, omitting halo Aq-F which undergoes twomajor mergers at z ∼ . z = 0. For each simulation we have128 output times, with a constant spacing of 155 Myr after z ≈ . Dark matter N -body simulations are ideal for followingthe non-linear growth of cosmic structure, including thedetailed accretion history of a halo. However, including aseparate, realistic treatment of baryons in a cosmologicalsimulation at sufficient resolution for stellar halo studies isstill computationally prohibitive. Semi-analytic galaxy for-mation modelling permits baryonic physics to be includedin much higher resolution N -body simulations at low com-putational cost, and is thus ideal for this purpose. Based onsimple theoretical treatments and empirical prescriptions,semi-analytic models can provide a good match to a remark-able number of observations of the local and distant galaxypopulations (Baugh 2006). These models, however, are notdesigned to follow the internal structure of galaxies in detail,hence the need for the particle tagging extension describedbelow.The Durham semi-analytic model, galform , is used inthis work to postprocess the Aquarius N -body simulations. galform computes the mass of stars that forms in eachhalo between every two successive simulation outputs, alongwith the properties of these stars, such as their metallicity.Rather than using the Bower et al. (2006) model on whichC10 was based, we use the more recent Font et al. (2011b)model, which includes a modified treatment of physics rel-evant to satellites. The only significant changes comparedto the C10 model are the use of a higher yield, a modifiedfeedback model for supernovae in which the mass ejectionefficiency saturates in haloes with circular velocity V circ (cid:54) − , and an earlier epoch of hydrogen reionization. TheFont et al. (2011b) model produces a satellite luminosityfunction that matches that of Milky Way and also generates stellar populations in satellites which match the observedluminosity-metallicity relation. This eliminates a shortcom-ing of C10 where the overall metallicity of the stellar halowas found to be significantly lower than that of the MilkyWay’s stellar halo. The C10 particle tagging method is a means to associate thestellar mass calculated by galform with six-dimensionalphase space volumes, defined by carefully chosen sets of rep-resentative particles in high-resolution dark matter only sim-ulations. This allows the 3D spatial distribution and kine-matics of stellar populations to be followed, albeit approxi-mately, without including gas physics in the original simu-lation.After applying galform to the outputs of the Aquariussimulations, C10 associated each newly formed stellar pop-ulation (discretized in age by simulation output times) withthe 1% most tightly bound dark matter particles in their cor-responding dark matter haloes. These sets of ‘tagged’ darkmatter particles are chosen at the time when the populationforms. Their diffusion in configuration and velocity spacecan then be tracked to the present day. This method is idealto study the formation of the outer accreted component ofstellar haloes, which are thought to be created by the tidaldisruption of satellites. The fraction of most-bound particlestagged, the free parameter introduced by the method, has aneffect on the scale radii of the resulting galaxies. C10 fixedthis value at 1 per cent in order to reproduce the distributionof observed half-mass radii for surviving satellites. Havingfixed this scale, satellite profiles and velocity dispersions alsomatch well to observations (see C10).Most surviving satellites are highly dark matter domi-nated (e.g. Walker et al. 2007). Since the N -body potentialis not altered by the gravity of the ‘stars’ represented bythe tags, the technique is likely to be much less accurate inbaryon-dominated systems, such as the centre of the MilkyWay analogue (see C10, Cooper et al. 2013 and Cooper etal. 2014 for further details and discussion). Satellite growthrates, survival times and orbits may also be affected by amore self-consistent treatment baryons (e.g. Sawala et al.2013). The relative importance of these effects is uncertain,however. For example, supernova feedback can create darkmatter cores in haloes (e.g. Navarro et al. 1996; Pontzenet al. 2013). This reduction in central density may reducetheir survival time when they are subjected to tidal strip-ping (e.g. Pe˜narrubia et al. 2010). On the other hand, thecentral concentration of stars may also bind these galaxiesmore tightly in a self-consistent simulation – in this contextwe note that numerical softening of the gravitational forcebetween particles creates artificial cores in our simulatedhaloes on scales of ∼
100 pc.To explore the consequences of such limitations, Bailinet al. (2014) compared tagged dark matter particles inhydrodynamical simulations to star particles formed self-consistently in the same simulations. They found moderatedisagreement in the distributions of halo stars representedby these two sets of particles, and concluded that tightlybound dark matter is not a reliable proxy for the dynamicsof stars in satellite galaxies. However, there is a fundamentaldifference between the tagging method tested by Bailin et c (cid:13)000
100 pc.To explore the consequences of such limitations, Bailinet al. (2014) compared tagged dark matter particles inhydrodynamical simulations to star particles formed self-consistently in the same simulations. They found moderatedisagreement in the distributions of halo stars representedby these two sets of particles, and concluded that tightlybound dark matter is not a reliable proxy for the dynamicsof stars in satellite galaxies. However, there is a fundamentaldifference between the tagging method tested by Bailin et c (cid:13)000 , 1–19 Ben Lowing et al al. (2014) and that used by C10. Bailin et al. (2014) carriedout a single tagging procedure for each satellite at the timeof its infall into the main halo, whereas C10 tag each stel-lar population at the time of its formation. Star particlesundergo considerable diffusion in energy over time, particu-larly if their host halo is disturbed, and hence would not beexpected to occupy the same region of phase space as tightlybound dark matter chosen long after they form. Dark matterparticles with similar binding energy chosen at the time ofstar formation, on the other hand, will undergo similar dif-fusion to star particles. This diffusion is therefore capturedin the C10 approach, but not in the method used by Bailinet al. (2014). Two forthcoming papers will test the taggingtechnique of C10 in a similar way (Le Bret et al. and Cooperet al. in prep.). Within the limitations and approximationsof the method, mentioned above and in C10, these new stud-ies find good correspondence between stars and tagged darkmatter particles in hydrodynamical simulations.
The final step is to turn the output of the tagged simulationsinto full mock catalogues of individual stars. The output ofthe model is a collection of N -body particles, each taggedwith some amount of stellar mass, and having both an agefor when the stars were born and a metallicity at that time.When a single dark matter particle is tagged multiple times,multiple ‘tags’ are created with the same phase space trajec-tory. Turning these tags into a full mock catalogue requirestwo steps. The first is to determine the total number of starsof different types that each tagged particle represents byconsidering each as a single-age stellar population, and thesecond step is to distribute the stars over the phase-spacevolume represented by the corresponding dark matter par-ticle. The following subsections discuss each step in detail. Unfortunately, even with the high-resolution Aquarius sim-ulations we are a long way from being able to represent eachindividual star by a single N -body particle. Instead, in theC10 simulations, each stellar particle can represent up to ∼ M (cid:12) . In order to convert these massive stellar parti-cles into individual stars we consider each tagged particle asa stellar population (SSP), in which all the stars formed atthe same time from gas with the same metallicity; the the-oretical justification of this is discussed in detail in Pasettoet al. (2012). The abundances of the different star types thatmake up the population can then be obtained using stellarpopulation synthesis modelling.The initial mass function (IMF) determines the numberof stars forming within each stellar population as a functionof initial stellar mass. Here we adopt a Kennicutt IMF withno correction for brown dwarfs (Kennicutt 1983), just aswas used in the semi-analytical model, to ensure consistency.Theoretical stellar isochrones then specify whether each typeof star still survives at the final time, along with the currentproperties of these stars.Stellar isochrones allow us to assign properties suchas temperature, magnitude and colour to stars of a giveninitial mass, age and metallicity. We employ the parsec isochrones (Bressan et al. 2012), as these are recent, up-to-date isochrones that extend over a wide range of metallicities(0 . (cid:54) Z (cid:54) . . < log t/ yr < . t ) = 0 . . (cid:54) Z (cid:54) . The next step is to split the massive stellar particles repre-senting SSPs into stars with individual positions and veloci-ties. Although the massive stellar particles provide a coarsesampling of the phase-space distribution of the stellar halo,it is desirable to smooth the distribution and increase thesampling. This is done by distributing stars over the entirephase-space volume that each tagged particle represents. Us-ing the phase-space volume rather than just the real-spacevolume preserves coherent structures such as tidal streams.It must be noted that increasing the sampling does not in-crease the resolution and cannot reveal any more detailedstructure than was in the original simulation.To estimate the phase-space volume each stellar massparticle occupies we have used the entropy-based binary de-composition (EnBiD) code of Sharma & Steinmetz (2006).EnBiD numerically estimates the densities of discretely sam-pled data based on a binary space partitioning tree. En-BiD uses an entropy-based node splitting criterion to decidewhich subspace (real or velocity space) to split next. Thismakes the method metric free, as there is no need to specifyin advance how velocity and configuration spaces relate toeach other, and ensures the subspace with greater variationis subdivided more.Once the binary tree has been built, the leaf nodes canbe used to estimate the phase-space volume of each particle.The volume of the hypercube of the leaf nodes can vary sig-nificantly even between close neighbours due to Poisson sam-pling noise of the original particle distribution. Simply usingthe volumes of raw nodes would result in a large dispersion,but this can be reduced by using a smoothing scheme. Wehave modified EnBiD’s anisotropic kernel density estima-tion scheme to return a volume as well as a density. Thescheme works by first calculating a local metric about eachparticle based on the shape of the hypercube of each leafnode and using the length of the six sides to scale the phase c (cid:13) , 1–19 aking Mocks of the Stellar Halo − −
100 0 100 200 x [ kpc ] − − y [ k p c ] − −
100 0 100 200 x [ kpc ] − −
100 0 100 200 x [ kpc ] − − −
200 0 200 400 600 v r [ kms − ] v θ [ k m s − ] − − −
200 0 200 400 600 v r [ kms − ] − − −
200 0 200 400 600 v r [ kms − ] Figure 1.
Real-space projection of an example satellite (upper row) and velocity-space projection (bottom row). The left column showsthe original set of satellite particles, while middle column is the down-sampled version with just 10 percent of the number of particles.The right column shows the result of resampling 10 per cent of the particles back to the original resolution using phase space kernelsbased on EnBiD density estimates. space in each direction. A covariance matrix is then calcu-lated based on the nearest 64 neighbours in the scaled space.Diagonalizing the covariance matrix gives a set of eigenval-ues and eigenvectors that can be used to perform a furthertransformation to a coordinate system in which the eigen-vectors define the principal axes. In this new rotated andscaled space the local particle distribution will have unit co-variance in each direction. Finally, the phase-space volumeof the particle is defined as a hypersphere with 1/40 of thevolume of a hypersphere based on the distance to the 40thnearest neighbour, R . It therefore has a radius R of R = (cid:18) (cid:19) / R . (1)The modified version of EnBiD returns a 6x6 matrix thatdescribes the phase-space volume of each particle.The actual sampling is performed by choosing pointsfrom an 6d isotropic Gaussian distribution in which eachcomponent has zero mean and variance σ = γR . Eachgenerated point is then transformed by applying the inverseoperations to convert from the sampled space back to po-sitions and velocities in the original simulation. Some carehas to be taken with the choice of the scale parameter γ which relates the scale of the Gaussian kernel to the scaleof the hypersphere returned by EnBiD. Choosing too largea value of γ will result in the phase space structure of eachcomponent of the halo being made hotter and too diffuse.Choosing too small a value of γ will result in an artificiallyclumpy phase space distribution. As a compromise we havechosen γ = 1 /
48. In a separate study (Wang et al in prepa-ration) we have developed a maximum likelihood methodfor modelling the distribution function of the stellar halo inorder to estimate the mass of the host dark matter halo.Using this approach we have found that the recovered halomass is biased high if γ is too large. With γ = 1 /
48 the levelof the bias is controlled to be less than about 10%, which ismuch smaller than typical measurement errors. Smaller val-ues of γ would remove this bias entirely, but would negatethe benefit of distributing stars over a phase space kernel,because in the limit γ → γ , we attempt to reconstruct a debris streamfrom one of our simulations after artificially degrading its c (cid:13)000
48 the levelof the bias is controlled to be less than about 10%, which ismuch smaller than typical measurement errors. Smaller val-ues of γ would remove this bias entirely, but would negatethe benefit of distributing stars over a phase space kernel,because in the limit γ → γ , we attempt to reconstruct a debris streamfrom one of our simulations after artificially degrading its c (cid:13)000 , 1–19 Ben Lowing et al resolution by sampling only a fraction of its particles. Weextracted from the simulation a satellite in the process ofundergoing significant tidal stripping, such that its particleshave an extended distribution in both configuration and ve-locity space. This object comprises 15,250 particles. A lowerresolution version was created by randomly selecting 10 percent of the particles. Finally, the low resolution version wasthen resampled using our EnBiD method, with each parti-cle being converted back into 10 particles, thus creating anobject with the same resolution as the original satellite. Acomparison between the two then demonstrates the successof our phase space sampling method.Fig. 1 shows projections of phase space for the differentversions of the satellite. The top row is a real-space projec-tion and the bottom row a velocity projection. The originalsatellite (left column) has a large amount of structure inboth real and velocity space, with several complete wrapsof the tidal tails visible. In the down-sampled version (thesecond column from the left) the primary features are stillpresent but the subtle ones are missing, while the remnant ofthe satellite core is more prominent. The EnBiD up-sampledversion (the third column from the left) is remarkably simi-lar to the original satellite, though the subtle features tendto be smoothed out. In particular the arc in velocity spacecorresponding to an orbital apocentre at v θ ≈
400 km s − ,while still present, is less well resolved. This is expected, asinformation is lost during down-sampling and no resamplingscheme will be able to recover it perfectly.A further test of how well the original and EnBiD ver-sions match is a comparison between the volume distributionfunction of phase-space density, v ( f ), as shown in the upperpanel of Fig. 2. We have used the technique of Arad et al.(2004), based on Delaunay tessellation, as an independentmeasurement of the phase-space density. The phase spacevolume distribution function is found to be similar in allthree cases, suggesting that it is robustly defined by evena small subset of the particles and that the EnBiD resam-pling process does not significantly alter it. The differencescan be seen more clearly in the bottom panel of Fig. 2. Therelative difference between the down-sampled volume distri-bution function and the original is shown by the red curve,while the relative difference between the EnBiD up-sampledversion and the original is show by the green curve. We havealso tested how much noise is introduced to v ( f ) simply byrandomly sampling the phase space distribution defined bythe original particles (as estimated by EnBiD) at a ratio of1 : 1 (grey line). The relative difference in all these cases ismostly within 10 to 50 per cent, as indicated by the two hor-izontal black dashed lines. There is no systematic trend inthe relative difference across 9 orders of magnitude in phasespace density. A few spikes correspond to larger differencesin particular regions, which appear to be almost entirelydue to the shot noise introduced by the discreteness of theparticle distribution.We also notice that at the low density end of Fig. 2the red dashed curve based on the down-sampled particlesstarts to plunge at f ∼ − . , whereas the blue curvebased on the original sample continues to increase down to f ∼ − . . The down-sampled particles are more sparselydistributed in phase space and thus sample low density re-gions poorly, as well as increasing the Poisson noise. Thegreen dot-dash curve, which is based on the particles resam- − − − − − − − − − − f − − l og | ∆ v ( f ) | / v ( f ) | Sampled − Original | / Original | EnBiD − Original | / Original | − Original | / Original l og v ( f ) OriginalSampledEnBiD1to1
Figure 2.
Top: Comparison of the volume distribution functionof phase-space density, v ( f ), for different phase space samplingschemes applied to a single satellite. Both axes have arbitraryunits. Blue (labeled ‘original’): the original high-resolution parti-cle distribution. Red (‘sampled’): a sample of 1 /
10 of the originaldistribution. Green (‘EnBiD’): the same number of particles asthe original distribution, drawn from the phase space kernels de-fined by the down-sampled (red) distribution. Grey(‘1-to-1’): eachparticle resampled once from the phase space kernels defined bythe original distribution. Bottom: Relative difference between thevolume distribution functions of the ‘sampled’, ‘EnBiD’ and ‘1-to-1’ results to the ‘original’ volume distribution function. Blackdashed lines mark the 10% and 50% level of difference. pled with EnBiD, recovers the same small scale behaviourof the original particle distribution. The grey curve, whereeach particle is resampled once from the phase space ker-nel of the original sample, have almost the same small scalebehaviour of the original sample as well.When applying this resampling technique to the taggedparticles an important consideration is how the particles aredivided into subsets on which to run EnBiD. The particularsubset that is chosen defines which particles are consideredas phase-space neighbours. We choose to apply the EnBiDvolume estimation to all the tagged particles that form alonga given branch of the halo merger tree. This separates thetagged particles into sets based upon the halo they formedin. It is also possible to split the tagged particles by in-dividual assignments, i.e. only particles from an individ-ual single-age stellar population are considered as potentialphase-space neighbours. This has the advantage that thephase-space volumes of particular populations are not di- c (cid:13) , 1–19 aking Mocks of the Stellar Halo luted by stars forming in the same halo at later times. Theproblem with this method is that sets of particles from in-dividual assignments can be very small in many cases, suchthat their phase-space volumes are too poorly sampled toaccurately constrain the structure in six dimensions. Whileour implementation allows either method to be used, wehave chosen to use the first method in the rest of this paper.EnBiD requires a minimum of 64 particles for the co-variance matrix calculation. We therefore exclude assign-ments with smaller numbers of particles than this. We havechecked that this eliminates a negligible amount of stellarmass. In this section we describe the range of information providedby our mock catalogues, and present a set of example mocksgenerated from the Aquarius simulations. The combinationof the various techniques used results in a rich dataset.
Five complete mock catalogues for the Aquarius haloesA-E are publicly available at the following website http://galaxy-catalogue.dur.ac.uk:8080/StellarHalo . Eachmock is stored in a separate table within the StellarHalodatabase and can be accessed via the web interface usingSQL queries. The database will contain all stars with ab-solute magnitude M g <
7. This cut is intended to selectbright stars most useful to studies of the distant halo, in-cluding main sequence turn-off stars, red giant branch starsand blue horizontal branch stars (see 3.4). It excludes a verylarge number of faint main sequence stars which are not ac-cessible to surveys of the outer halo of the Milky Way.Table 1 lists the information that is provided. Positions,velocities and apparent magnitudes are given relative toboth the centre of the halo and to an observer placed in aposition similar to that of the Sun. In order to do this wefirst define a plane that a galactic disc would likely lie inand then choose a position for the Sun within that plane.We choose to orient the disc plane perpendicular to the mi-nor axis of the moment of inertia of the dark matter inthe inner 10 kpc, as the minor axis is a stable orbital axis.While the inclusion of baryons might be expected to changethe shape of the centre of the halo due to adiabatic con-traction (Kazantzidis et al. 2004), their presence is unlikelyto alter the halo’s orientation. We have taken the Sun tobe 8 kpc from the centre of the galaxy and moving with avelocity ( U, V, W ) (cid:12) = (11 , ,
7) km s − (Sch¨onrich et al.2010). The observer has been placed on the x -axis, thoughany orientation on the solar circle would be an equally validchoice.In addition to the magnitudes, we provide a few intrin-sic properties of stars, such as their age, metallicity, stellarmass, surface gravity and effective temperature. The semi-analytical galaxy formation model galform uses physically-motivated prescriptions to follow the evolution of the galax-ies within each halo. While galform outputs many proper- The magnitudes are based on sdss u, g, r, i, z filters. ties for each galaxy, it is only the star formation and chem-ical history that are used to build the mocks. These twomodel predictions for the baryonic components of satellitesare essential quantities required to populate the halo withstars. The final input are the theoretical isochrones that pro-vide the detailed properties of the stars. For a given age andmetallicity, together with an IMF, the isochrones providea breakdown of these stellar properties. The type columnis the original ‘stage’ flag in the PARSEC isochrones thatdenotes different phase of stellar evolution (see Table 1).Beyond the spatial and kinematical information for in-dividual stellar populations, the underlying dark matter N -body simulations provide the history of halo mass assembly.Each star is associated with a dark matter particle and thesubhalo to which that particle is currently bound, throughthe dark matter ID and Subhalo ID. This information al-lows the dynamical history of the particle to be traced backthrough the simulation. The satellite which brought the starinto the main halo can be identified, along with the time ofaccretion of that satellite and the time at which the particlewas stripped from the satellite. It is possible to find stellarstructures and streams easily by grouping stars according totheir satellite of origin, using the TreeID field. The ‘infall’redshift of the parent satellite has been provided, defined asbe the redshift when the satellite reached its maximum mass(for most objects this correlates closely with the actual timeof entry into the main FOF group and hence with crossingthe virial radius of the central halo). For generality and simplicity, the observables presented inour catalogue have not been convolved with sources of obser-vational error. In order to simulate a particular survey, theuser should therefore convolve physical quantities such asmagnitude, velocity and metallicity with appropriate errorfunctions. The user may also wish to impose a foregrounddust extinction model, and superpose additional source dis-tributions to account for contaminations of halo star selec-tions by faint foreground stars and the extragalactic back-ground.
While the mock catalogues we have generated are based onsophisticated cosmological models, they are not without lim-itations. The biggest limitation is that our model only gen-erates the accreted component of the stellar halo. It doesnot include processes that create in situ halo components inhydrodynamical simulations, such as the scattering of discstars or star formation in gaseous streams. This means thatthe model is most suitable for studying the outer regionsof a galaxy ( r >
20 kpc). There are various possible waysto add in additional components to create a more completemodel of the entire galaxy, and we will explore these in fu-ture papers.Stars belonging to the galactic disc are not included inthe mocks and the effect of having a disc present withinthe halo during its evolutionary history is also missing. Sub-haloes that pass near or through a disc are likely to be dis-rupted more rapidly (D’Onghia et al. 2010). The importance c (cid:13)000
20 kpc). There are various possible waysto add in additional components to create a more completemodel of the entire galaxy, and we will explore these in fu-ture papers.Stars belonging to the galactic disc are not included inthe mocks and the effect of having a disc present withinthe halo during its evolutionary history is also missing. Sub-haloes that pass near or through a disc are likely to be dis-rupted more rapidly (D’Onghia et al. 2010). The importance c (cid:13)000 , 1–19 Ben Lowing et al
Field Data type Unit DescriptionID long ID of the star, unique within the mock cataloguePosition ( x, y, z ) float kpc Star position, relative to halo centreVelocity ( v x , v y , v z ) float km s − Stars velocity, relative to halo rest frameGalactic coordinates ( l, b ) float o Galactic longitude and latitude relative to solar observerRadial distance ( r ) float kpc Distance from solar observerRelative velocity ( v r , v θ ) float km s − Radial and tangential velocity relative to solar observerlog (Age) float log (yr) Age of starMetallicity float Fraction of star mass in metalsMass float M (cid:12) Current mass of starlog ( g ) float log(cgs) Surface gravity of the starlog ( T eff ) float log (K) Effective temperature of the starType int The original PARSEC evolutionary stage flag that labels the type of star(0=Pre-MS, 1=Main sequence, 2=Subgiant branch, 3=Red giant branch,4,5,6=stages of core He burning, 7,8=stages of Asymptotic Giant Branch) M ( u,g,r,i,z ) float Absolute rest frame ( u , g , r , i , z ) band ( sdss ) magnitude of the star m ( u,g,r,i,z ) float Apparent rest frame ( u , g , r , i , z ) band ( sdss ) magnitude of the staras seen by the solar observerDark Matter ID long The tagged dark matter particle this star was spawned fromSubhalo ID long The current subhalo to which the associated dark matter particlebelongs to. -1 if not currently a member of any subhalo, 0 if partof the main haloTree ID long The subhalo tree that this star came from. All stars that arrive in themain halo via the same satellite have the same IDInfall Redshift float Accretion redshift of the subhalo this star once belonged to. Defined as theredshift at which the satellite reaches its maximum mass (see text). Table 1.
Information provided in our online mock catalogues that may be accessed at http://galaxy-catalogue.dur.ac.uk:8080/StellarHalo . SQL can be used to query any of these fields. of such effects is unclear because the disc has been grow-ing as the halo grows; even today, it presents a small crosssection to infalling satellites. This effect is therefore mostrelevant to the small number of satellites that pass close tothe centre of the halo after at late times. The consequencesof neglecting enhanced disruption by the galactic disc arelikely to vary among our five simulations, according to theirparticular mass assembly histories and subhalo orbital dis-tributions. While the properties of individual streams maychange significantly, we expect that, in most cases, the over-all properties of the halo should be robust to the inclusionof a realistic disc.As discussed in section 2.3, the C10 technique uses darkmatter particles to trace the distribution and dynamics ofstars. In reality, stars may form on phase space trajectoriesthat are not well sampled in collisionless simulations (forexample, centrifugally supported discs). For this reason, thespatial configuration of stars within surviving satellites maynot be reliable in detail. The importance of the initial con-figuration is greatly diminished once a satellite has been dis-rupted and the stars phase-mixed into the halo, but studiesof chemical and kinematic gradients along coherent streamsshould keep this limitation in mind.
We now describe a set of example mock catalogues gener-ated from the Aquarius simulations and use them to illus-trate how particular samples of stars can be easily selected.Each sample is defined by a set of colour and magnitudecuts that selects a particular stellar type commonly used forsurveying the stellar halo. The overall colour-magnitude di-agram for all the stars in the stellar halo of Aq-A (with the − . . . . . g − r − M g
12 3 . . . . . . . . l og N Figure 3.
Colour-magnitude diagram for the overall stellar pop-ulation of the accreted stellar halo of Aq-A. Stars still boundto satellites are not included. The colours encode the number ofstars within a bin and the boxes illustrate the approximate regionscovered in our three mock catalogue selection functions (see textfor details of the actual selection criteria, which make use of ad-ditional photometric bands and metallicity information). Box 1marks our selection for MSTO stars, box 2 marks the Bell et al.(2010) selection for BHB stars and box 3 the Xue et al. (2014)selection for K giants. Photometric magnitudes are given in theSDSS ugriz filter system. contribution from stars still bound in satellites removed) isshown in Fig. 3. The selections for our three samples of halotracers are marked by boxes.The vast majority of stars are low mass stars that havenot yet evolved off the main sequence and account for a c (cid:13) , 1–19 aking Mocks of the Stellar Halo significant fraction of the total stellar mass. However, asmentioned above, these are very faint. They can only beobserved in the Solar Neighborhood and are of little interestto studies of the distant halo or its global properties. A clearred giant branch and horizontal branch are visible; thesecontain the older, brighter stars which can be observed tolarge galactocentric distances and which are thus are mostuseful when studying the structure of the halo. Main sequence turn off (MSTO) stars can be unambiguouslyselected using the colour cut 0 . < g − r < . M g >
4, in order toonly select dwarf stars and remove giants. MSTO are rela-tively numerous and therefore serve as a good proxy for thegeneral stellar content of the stellar halo. Approximately40% to 50% of the stellar mass brighter than M g = 7 isfound in MSTO stars and there are 81,506,853 such stars inthe Aq-A halo. The volume currently accessible using MSTOstars is a relatively small fraction of the total halo, but theirhigh surface density in this volume allows a great deal of in-formation to be gathered by wide-field photometric surveyssuch as SDSS (e.g. Juri´c et al. 2008). Blue horizontal branch (BHB) stars are very often used astracers of the distant stellar halo. They are luminous andhave a very nearly constant magnitude, M g ∼ .
7, so theycan be located accurately to very large distances along theline of sight even with photometry alone. Unfortunately,theconditions under which BHB stars form are still not wellunderstood. This makes their abundance in our mock cat-alogues highly sensitive to the assumptions made by thetheoretical isochrones we use.To select BHB stars in our mock catalogues we use thecolour criteria proposed by Bell et al. (2010): 0 . < u − g < . − . < g − r < − .
06 and excluding the region([ u − g − . / . + ([ g − r + 0 . / . <
1. With thisselection, BHB stars make up a small fraction of the totalmass in stars brighter than M g = 7, ranging from 1-3% inthe five mocks. There are 239,950 BHBs in the Aq-A halo.The main source of stellar contamination in photomet-ric BHB selection is fainter blue straggler (BS) stars. Theformation of blue stragglers is also not well understood, andthey are not included in the parsec isochrones. Hence, evenhalo BS (and white dwarfs) are not a source of contam-ination in the mocks. This should be kept in mind whencomparing to observed photometric BHB counts, as shouldthe absence of quasar contamination at faint magnitudes. K giants are another stellar type frequently used to mapthe Milky Way’s stellar halo (Morrison et al. 1990, 2000;Starkenburg et al. 2009). They have the advantage that theyare found in predictable numbers in old populations of allmetallicities. However, unlike other standard candles, theirluminosities range by two orders of magnitude, which makesobtaining accurate distances more challenging. K giants are
10 100 r [kpc]10 − − − ρ [ M (cid:12) k p c − ] Aq-AAq-BAq-CAq-DAq-E
Figure 4.
The spherically averaged density profiles of the fiveAquarius stellar haloes. The arrows show the location of the breakin a best-fit broken power law. Three out of the five haloes haveclear breaks, a fourth has a weak break while Aq-E is well fittedby an unbroken power law. one of the types that have been targeted for spectroscopy bythe Sloan Extension for Galactic Understanding and Explo-ration ( segue , Yanny et al. 2009) and recently, Xue et al.(2014) have estimated distances to a sample of 6036 K giantsusing a combination of observed g − r colours and spectro-scopic metallicity.In order to isolate K giants in our stellar halo modelswe have used the colour cuts 0 . < g − r < . . .
087 [Fe / H] + 0 .
39 [Fe / H] + 0 . M g < −
7% of the total stellar halo mass in stars brighterthan M g = 7 is found as K giants and there are 956,246 Kgiants in the Aq-A halo. One of the most basic, but still unanswered, questions ingalactic structure is the density profile and shape of theMilky Way stellar halo. While there is now general consensusthat the stellar halo is oblate and follows a broken power-law density profile, the exact form is still ill constrained.Evidence for a break has been found by studies using RRLyrae variables as standard candle tracers of the distantmetal-poor halo, with measured break radii of ∼
25 kpc(Watkins et al. 2009), ∼
30 kpc, (Sesar et al. 2010) and ∼
45 kpc (Keller et al. 2008; Akhter et al. 2012). Measure-ments using other types of tracers have also suggested the c (cid:13)000
45 kpc (Keller et al. 2008; Akhter et al. 2012). Measure-ments using other types of tracers have also suggested the c (cid:13)000 , 1–19 Ben Lowing et al
10 100 r [kpc]10 − − − − − − ρ [ M (cid:12) k p c − ] OverallMSTOBHBK giants
Figure 5.
The spherically averaged density profile of the Aq-A stellar halo obtained using different types of tracer stars asindicated in the legend. The black line shows the density profileof all stars in the halo. existence of a break. Deason et al. (2011), using a large sam-ple of BHB and blue straggler stars from sdss
Data Release8 (DR8), again found the profile to follow a broken powerlaw, with a break radius at ∼
27 kpc. Similarly, from a sam-ple of near main-sequence turn off (MSTO) stars from theCanada-France-Hawaii Telescope Legacy Survey ( cfhtls ),Sesar et al. (2011) measured the break to be at ∼
28 kpc. Incontrast, no break or slope steepening was seen by De Pro-pris et al. (2010), who analyzed a sample of 666 BHB starsin two different directions separated by about 150 ◦ on thesky. Instead they found a smooth stellar distribution obey-ing a power-law density profile of index ∼ -2.5. Akhter et al.(2012) suggested that an explanation for these differences inthe measured location of the break could be that, at largerdistances, the outer halo is not be smooth but dominated byclumps and debris. Other studies of the power law index (orindices) of the Milky Way halo include Chiba Beers (2000);Miceli et al. (2008); Carollo et al. (2010).In Fig. 4 we show the density profiles of the accretedstellar haloes in our five simulations, and we have excludedstars in surviving bound subhaloes. Three of them showclear breaks, although generally at larger distances than es-timated for the Milky Way. To check whether the type oftracer makes a difference we also show in Fig 5 the stel-lar density profile of the Aquarius A halo measured fromMSTO, BHB and K giant stars. The overall density profileis that of all stellar mass in the tags from the C10 method,while the tracer profiles are from samples taken from themock catalogued described in Section 3.4. All of the curveshave the same shape, with the break occurring at the sameradius, but with differing normalizations, depending on theabundance of the stellar type.It has been suggested that the presence of a break in theMilky Way’s density profile occurs at the transition betweenan inner region where the halo is dominated by in-situ starsand an outer region which is primarily composed of accretedhalo stars. However, this cannot be the case for the breaksseen in the profiles of the mock catalogues, since they donot contain any in-situ component. Recently, Deason et al.
10 100 r [kpc]10 − − − − − − ρ [ M (cid:12) k p c − ] Overall1st-12.70Gyr2nd-10.67Gyr3rd-12.17Gyr4th-9.73Gyr5th-12.71Gyr
Figure 6.
The density profile of the Aq-A stellar halo brokendown into the five satellites which each contribute the largestamount of mass to the stellar halo. Numbers next to the legendshow the look-back infall time of the five satellites. (2013) have proposed that the breaks in stellar halo profilesare related to the satellite accretion history of their galax-ies and that a massive accretion event occurring 6 to 9 Gyrago is likely to have caused the break in Milky Way densityprofile. Satellites that are just falling into the halo and arestill in the process of being tidally disrupted, may have al-ready lost a significant fraction of their stars through tidalstripping. Although these stripped stars are unbound andnow form part of the halo, their tidal debris is still in theform of large coherent structures roughly tracing the orbitsof their progenitor satellites. The merit of constructing agalactocentric, spherically averaged density profile for thesedebris streams and clouds is questionable, as they are farfrom being spherical or dynamically relaxed. Such profilesgenerally appear flat in the central regions and break to avery rapid falloff at their outer edge.In Fig. 6 we show the density profiles of the stars thatwere brought in by the five largest contributors to the Aq-A stellar halo. Numbers beside the legends correspond tothe look-back time of infall for these satellites. Except forthe very largest contributor, the profiles of the others allhave very strong breaks. A visual inspection reveals that thelargest contributor was accreted earliest at about 12.7 Gyrago. By the end of the simulation, it has been completelydisrupted and mixed, with a correspondingly smooth profile.In contrast, the 4th largest contributor has fallen in morerecently 9.73 Gyr ago and is still in the form of a coher-ent, massive stream. It has a sharp break around ∼
80 kpc.The other three contributors are in intermediate stages ofdisruption; all show clear caustic/shell structures (visible aslocalized density peaks). Notice the fifth most massive con-tributor shows these structures despite having been accretedat roughly the same time as the largest contributor, which iswell mixed. The largest contributor has about 4.5 times themass of the fifth most massive, and hence the decay of itsorbit through dynamical friction may have been more rapid.In Aq-A, breaks in the individual contributor profilescorrespond to these caustic structures (this is true also forthe other four haloes, which we do not plot here). This find- c (cid:13) , 1–19 aking Mocks of the Stellar Halo ing is consistent with the results of Deason et al. (2013), thatbreaks are related to the apocentres of accreted satellites.These few most massive contributors determine the over-all profile – over certain ranges of radius, almost the entirestellar mass comes from a single progenitor. For example,around 50 kpc the overall density profile is almost entirelydetermined by the 2nd largest contributor. It is clear thatthe break in the overall profile in this region is driven bystars from this object. It can be clearly seen that the num-ber of breaks in the overall profile, and their sharpness, willdepend strongly on the relative distribution of stars from thefew satellites that contribute significant amounts of debristo the halo.It can be argued that when measuring the sphericallyaveraged stellar halo density profile, coherent structuressuch as those just described should be removed. However,while a massive stream such as the Sagittarius stream in theMilky Way might be easy to exclude, objects in later stagesof disruption are more difficult to identify. In the simula-tions we are able to use the subfind halo finder to identifybound structures and follow the accretion histories of starsin order to find candidate streams but this cannot be donewith observational data. In addition, even if structures havebeen identified, the question is then how to decide whetherthey should be excluded from the halo. For these reasons,as the overall profile is so sensitive to individual accretionevents, and thus is likely to change rapidly with time, it isnot clear how useful it is to measure the density profile of astellar halo, other than for the purpose of gaining informa-tion about recent large accretion events. As shown in Fig. 5, the spherically averaged density profileis not sensitive to the type of tracer used. What is muchmore important is the area of the sky over which the profileis measured. In order to measure the profile to large radii,narrow pencil-beam surveys of bright tracers are often used,probing deep into the stellar halo over a limited area of sky.In order to quantify the intrinsic differences in the slope ofthe profile in different directions we employ a similar ob-servational strategy to Sesar et al. (2011) and measure thedensity the simulated halo stars over in patches of small solidangle. Sesar et al. (2011) used four ‘pencil beams’ from theCanada-France-Hawaii Telescope Legacy Survey ( cfhtls )to select MSTO stars in order to measure the slope of thestellar halo density profile between 5 and 35 kpc in heliocen-tric radius. A direct comparison with Sesar et al. (2011) islimited by the lack of in situ halo stars in our model; it maybe that these dominate the inner part of the halo. Moreover,Sesar et al. (2011) attempted to excise the most obvious ac-creted substructures from their analysis (the Sagittarius andMonoceros streams). This is a subjective procedure and im-plicitly assumes an underlying ‘smooth’ halo. We do not tryto excise visible streams, or even bound satellites . Our aimin this section is to demonstrate how substructure in the haloaffects the determination of the density profile, rather thanto compare directly with the results of Sesar et al. (2011). Note that we remove bound satellites from the analysis in allother sections. . . . . . . . f ( n ) Aq-A . . . . . . . f ( n ) Aq-B . . . . . . . f ( n ) Aq-C . . . . . . f ( n ) Aq-D − − − − − n . . . . . . . f ( n ) Aq-E
Figure 8.
The distribution of volume density profile power-lawslope measured for halo stars in our mock catalogue selected usingthe criteria of Sesar et al. (2011, see text) in patches of area8 sq. degree tiling the whole sky. Each panel shows a differentmock catalogue. Bound subhaloes and other overdensities have not been removed; they contribute a small fraction of the pixels inthese area-weighted distributions. The dashed vertical lines markthe overall slope from a power law fit to the (mass-weighted)spherically averaged density profile in the same volume.c (cid:13)000
The distribution of volume density profile power-lawslope measured for halo stars in our mock catalogue selected usingthe criteria of Sesar et al. (2011, see text) in patches of area8 sq. degree tiling the whole sky. Each panel shows a differentmock catalogue. Bound subhaloes and other overdensities have not been removed; they contribute a small fraction of the pixels inthese area-weighted distributions. The dashed vertical lines markthe overall slope from a power law fit to the (mass-weighted)spherically averaged density profile in the same volume.c (cid:13)000 , 1–19 Ben Lowing et al n slope map(5kpc Top: A heliocentric all-sky map in which the colour of each 8 sq. degree pixel corresponds to a value of n , the best-fit(dimensionless) power law slope of the galactocentric density profile of stars from our mock catalogues, selected according to the MSTOcolour-magnitude criteria of Sesar et al. (2011, see text) in the heliocentric distance range 5 < r h < 35 kpc. The Galactic Centre isat the centre of this map and our fiducial disc plane is oriented along the equator. Unlike previous figures, bound satellites have not been excised in this analysis, although there are only three in this distance range (black circles, diameter equal to 10 times the apparentangular size of the corresponding subhalo). Bottom: The logarithm of the projected surface number density of stars selected and observedin the same way as the top panel. Labels A and B mark the location of satellites discussed in the text. We select stars in the same colour-apparent magnituderange as Sesar et al. (2011) (0 . < g − r < . g > 17, 17 Heliocentric sky map similar to Fig. 7, but showing the power law density slope (top) and surface number density (bottom)of all stars in our mock catalogue with galactocentric distance 20 < r g < 80 kpc. Bound satellites and other overdensities have not beenexcised. Black solid circles mark subhaloes that host satellites, with the angular size of circle equal to 10 times the half mass radius ofthe subhalo as viewed from our fiducial Solar position. expressed in terms of galactocentric radius, r g , hence pencilbeams in different directions sample different ranges of r g .We tile the entire sky with equal-area pencil beams andcompute the parameters of the best-fit power law densityprofile, ρ ∝ r n , in each beam. The top panel of Fig. 7 showsa sky map of Aq-A with pixels colour-coded by the value of n along a given beam. There are large variations in n acrossthe sky, with profiles ranging from very steep ( n ∼ − 5) insome areas to flat ( n ∼ 0) in others. The large-scale gradient in the value of n with galacticlatitude is likely to result from the non-spherical overall dis-tribution of halo stars. This is seen clearly in the projectednumber density of MSTO stars (including those bound tosatellites) shown in the lower panel of Fig. 7. The surfacedensity of accreted stars is high towards the galactic centreand falls off towards the galactic poles. This panel also showsthat variations of n on smaller scales (hundreds of square de-grees) correspond to localized density substructures in the c (cid:13)000 0) in others. The large-scale gradient in the value of n with galacticlatitude is likely to result from the non-spherical overall dis-tribution of halo stars. This is seen clearly in the projectednumber density of MSTO stars (including those bound tosatellites) shown in the lower panel of Fig. 7. The surfacedensity of accreted stars is high towards the galactic centreand falls off towards the galactic poles. This panel also showsthat variations of n on smaller scales (hundreds of square de-grees) correspond to localized density substructures in the c (cid:13)000 , 1–19 Ben Lowing et al halo. Some of the most dramatic fluctuations occur aroundsurviving galaxies that are being tidally stripped; for exam-ple the density ‘hot spots’ in the lower panel, labeled ‘A’ and‘B’, correspond to regions in the upper panel with n ∼ < r h < 35 kpc; we have overplotted black cir-cles on Fig. 7 to highlight all those in Aq-A, with relativediameters scaled to 10 times the angular size of their sub-halos (as seen by our fiducial Solar observer). Most smallscale fluctuations in the measured density slope, are due toinhomogeneities in the stellar halo itself rather than boundsubstructures.In order to quantify how much the slope varies over thesky, the distributions of n are shown in Fig. 8 for all fivehaloes. These are ‘area weighted’ distributions, in that eachbin shows the fraction of the total number of 8 deg × mode of the distri-bution, but for Aq-D and E the overall slope is significantlydifferent, lying at the steep end of the distribution. This canbe explained by the fact that, as we saw earlier, a substruc-ture or stellar stream can completely dominate the overallprofile, making it radically different from that measured ina typical patch of the halo. It should be noted that, so far,we have been considering the slope of the profile to no morethan 35 kpc from the Sun. In the Milky Way, it has beensuggested that an in situ component makes a significantcontribution to the halo out to galactocentric distances of r g ∼ 20 kpc. This component may also be much smootherthan the accreted halo originating from disrupted satellites.A smooth in situ component may reduce the variation inthe density profile.To investigate the structure of the outer halo, we takeall stars brighter than M g = 7 in the mock catalogue withgalactocentric radii in the range 20 < r g < 80 kpc. Thisupper limit is chosen because the overall density profiles ofFig. 4 do not show breaks out to 80 kpc. We again constructpencil beams to a fixed distance from the Sun and measurethe slopes of the density profiles from the galactic centre.The all sky map of slopes is presented in the top panel ofFig. 9 and the surface density of stars in the lower panel.Compared with the top panel of Fig. 7, the fluctua-tions across the sky in Fig. 9 are more numerous and pro-nounced. The large scale variations in n no longer trace largescale variations in the surface density of stars (bottom plot).However, we can see many of the small scale variations inthese two maps correspond. Since stars within 20 kpc to thegalactic centre are excluded, the central overdensity in Fig. 7disappears.The more complex fluctuations in Fig. 9 arise becausethere are more substructures at these larger galactocentricradii, and the local variations in the stellar halo are morepronounced. Again, we mark all subhaloes that are moremassive than 10 M (cid:12) with black circles. Most of these sub-haloes are distant and are thus small in their apparent size.The three most prominent local patches are associated withthree massive and extended substructures. Interestingly, thetwo patches that are close to the north galactic pole and onthe left hand side of the map, which can also be identified < / / / R c / a Aq-E < / / / R Aq-D c / a Aq-C KgiantMSTO Aq-B dark matter c / a Aq-A AllBHB Figure 10. The overall axis ratio (minor axis versus major axis)of different stellar tracers and the underlying dark matter, re-ported for four shells centred on the galaxy with inner radius20 kpc in all cases and outer radii of 1 / / 3, 1 / × R respectively. in Fig. 7, now give very steep slopes in contrast to their flatslopes in Fig. 7. This is because the associated subhaloes areclose to the Solar observer. The radial range of 5 to 35 kpcfrom the Sun has a very flat density distribution, beyondwhich the density profile drops very quickly. In Section 4.1 we showed that the shape of the sphericallyaveraged density profiles of different stellar tracers in ourmock catalogues are similar when averaged over the wholesky. In this section we examine the sky distribution of dif-ferent stellar tracers in more detail. We first look at theellipsoidal shape of different populations, and then considerthe degree of variation amongst different tracers in smallpatches of the sky. In this section we aim to quantify the overall ellipsoidalshape of the stellar halo. The degree of ‘flattening’ (or tri-axiality) of the stellar halo is of great interest, although tobe meaningful, measurements based on ellipsoidal fits to the3d distribution of stars require that a significant fraction ofthe halo be ‘smooth’ and symmetric.Fig. 10 gives the ratio between the ellipsoidal minor andmajor axes, c/a , obtained by diagonalizing the standard in-ertia tensor of different stellar populations (coloured lines,see legend) and dark matter particles (black lines) in eachsimulation. We consider stars and dark matter in four cu-mulative radial bins, the outer edges of which correspond to1 / 4, 1 / 3, 1 / × R . We exclude stars within 20 kpc c (cid:13) , 1–19 aking Mocks of the Stellar Halo < / / / R Aq-E < / / / R Aq-D a n g l e b e t w ee n m a j o r a x i s o f s t a r s a n d D M Aq-C KgiantMSTO Aq-B Aq-A AllBHB Figure 11. The angle (in units of degrees) between the majoraxis of stellar tracers and the major axis of dark matter particlesreported for four shells centred on the galaxy with inner radius20 kpc in all cases and outer radii of 1 / / 3, 1 / × R respectively. of the galactic centre and we do not include stars in anybound satellites.The overall shape of dark matter particles and stars aresomewhat different from each other. In the outer halo, withthe exception of Aq-E, the distribution of stars is signifi-cantly flatter than that of the dark matter. This reflects thetendency for stars even in the outer stellar halo to be con-centrated towards the main symmetry plane of the halo, asnoted by Cooper et al. (2011). As particles at larger radiiare included in the measurement, c/a remains approximatelyconstant for stars and dark matter. The large change in theoutermost point for Aq-E is likely to be because the stellardistribution at large radius is dominated by only one coher-ent tidal stream. Interestingly, we find that BHB stars havea slightly more isotropic distribution than the other tracers.Fig. 11 illustrates how the major axes of stellar tracerscorrelate with the major axes of dark matter particles . Forall the five haloes, we see the misalignment angle is typicallysmall (less than ∼ 20 degrees) particularly within R / We have checked the direction of the major axis of dark matterparticles is almost constant with radius dark matter out to R , despite the significant differencesin the flattening of these components seen in Fig. 10. It has been proposed that localized variations in the mixof stellar populations may provide a method for discoveringnew stellar structures in the Milky Way halo (e.g. Ibata et al.2007; McConnachie et al. 2009). Significant fluctuations inthe ratio of BHB to MSTO stars have been found in the stel-lar halo (Bell et al. 2010). Some of these, such as the “low-latitude stream” are almost devoid of BHB stars, while otherstructures are rich in them. The Sagittarius tidal streameven shows a variation in the BHB/MSTO ratio along itsextent (Bell et al. 2010), possibly related to a variation inmean metallicity that has been interpreted as the result ofa population gradient in the progenitor dwarf galaxy (e.g.Chou et al. 2007). In this section we perform a similar testwith the BHB, MSTO and K giant samples from our mockcatalogues and attempt to understand the origin of the dif-ferences.The simplest test is to search for variations by compar-ing sky maps of the projected map for our different tracersas seen by a hypothetical observer located on the solar circle.Fig. 12 (left column) shows a number of such maps for theAq-C halo. (Aq-D and Aq-B show similar features; to avoidrepetition we focus on Aq-C only.) It contains many clearlydefined streams that are visible in all the tracers. What ismost interesting, however, is that there are two streams thatdo not appear in the BHB map but appear in the two othermaps. In particular, the large continuous horizontal streamlying close to the equatorial plane that are clearly seen inthe MSTO and K giant map is partially missing in the BHBmap, while the shorter stream located near the top of theMSTO map, towards the galactic north pole, is totally ab-sent in the BHB map. T he absence of these streams in theBHB map helps to explain why BHB stars show distinct fea-tures in the behaviour of axis ratio and major axis at largeradii.Comparing the ratio of BHBs/MSTO stars (Fig. 12right column), the two missing streams are easily identifi-able. While the large, equatorial stream is at least partiallypresent in the BHB map, it has a lower than average abun-dance resulting in a higher MSTO/BHB ratio. In contrast,the stream near the top of the MSTO map is almost com-pletely absent in the BHB map and has a huge MSTO/BHBratio. Other streams with more standard stellar populationshave BHB to MSTO ratios similar to that of the overall haloand thus do not show up in the ratio maps. There are alsodifferences in the MSTO/K giant ratio, but in the oppositesense. From the MSTO/K giant map it can be seen thatthe two streams have a higher than average abundance of Kgiants. One other feature is that the BHB stars are slightlymore abundant in the galactic centre.These differences in the relative abundance of the threetypes of tracer can be explained by differences in the ageand metallicity of the stellar populations that make up thefeatures. In Fig. 13 we show a map of the average age (left)and the average metallicity (right) of the Aq-C stellar halo.These have been generated from the mean age/metallicity ofall stars within each equal-area patch of the sky. It is imme-diately obvious that the streams missing from the BHB map c (cid:13)000 20 degrees) particularly within R / We have checked the direction of the major axis of dark matterparticles is almost constant with radius dark matter out to R , despite the significant differencesin the flattening of these components seen in Fig. 10. It has been proposed that localized variations in the mixof stellar populations may provide a method for discoveringnew stellar structures in the Milky Way halo (e.g. Ibata et al.2007; McConnachie et al. 2009). Significant fluctuations inthe ratio of BHB to MSTO stars have been found in the stel-lar halo (Bell et al. 2010). Some of these, such as the “low-latitude stream” are almost devoid of BHB stars, while otherstructures are rich in them. The Sagittarius tidal streameven shows a variation in the BHB/MSTO ratio along itsextent (Bell et al. 2010), possibly related to a variation inmean metallicity that has been interpreted as the result ofa population gradient in the progenitor dwarf galaxy (e.g.Chou et al. 2007). In this section we perform a similar testwith the BHB, MSTO and K giant samples from our mockcatalogues and attempt to understand the origin of the dif-ferences.The simplest test is to search for variations by compar-ing sky maps of the projected map for our different tracersas seen by a hypothetical observer located on the solar circle.Fig. 12 (left column) shows a number of such maps for theAq-C halo. (Aq-D and Aq-B show similar features; to avoidrepetition we focus on Aq-C only.) It contains many clearlydefined streams that are visible in all the tracers. What ismost interesting, however, is that there are two streams thatdo not appear in the BHB map but appear in the two othermaps. In particular, the large continuous horizontal streamlying close to the equatorial plane that are clearly seen inthe MSTO and K giant map is partially missing in the BHBmap, while the shorter stream located near the top of theMSTO map, towards the galactic north pole, is totally ab-sent in the BHB map. T he absence of these streams in theBHB map helps to explain why BHB stars show distinct fea-tures in the behaviour of axis ratio and major axis at largeradii.Comparing the ratio of BHBs/MSTO stars (Fig. 12right column), the two missing streams are easily identifi-able. While the large, equatorial stream is at least partiallypresent in the BHB map, it has a lower than average abun-dance resulting in a higher MSTO/BHB ratio. In contrast,the stream near the top of the MSTO map is almost com-pletely absent in the BHB map and has a huge MSTO/BHBratio. Other streams with more standard stellar populationshave BHB to MSTO ratios similar to that of the overall haloand thus do not show up in the ratio maps. There are alsodifferences in the MSTO/K giant ratio, but in the oppositesense. From the MSTO/K giant map it can be seen thatthe two streams have a higher than average abundance of Kgiants. One other feature is that the BHB stars are slightlymore abundant in the galactic centre.These differences in the relative abundance of the threetypes of tracer can be explained by differences in the ageand metallicity of the stellar populations that make up thefeatures. In Fig. 13 we show a map of the average age (left)and the average metallicity (right) of the Aq-C stellar halo.These have been generated from the mean age/metallicity ofall stars within each equal-area patch of the sky. It is imme-diately obvious that the streams missing from the BHB map c (cid:13)000 , 1–19 Ben Lowing et al BHBs log ( N BHB ) K giants / BHBs N Kgiants /N BHB K giants log ( N Kgiants ) MSTOs / BHBs 100 1000 N MSTOs /N BHB MSTOs log ( N MSTOs ) MSTOs / K giants 70 300 N MSTOs /N Kgiants Figure 12. The distribution of different stellar tracers over the sky for the Aq-C stellar halo as seen by our fiducial heliocentric observer.The Galactic Centre is at the centre of the map and the disc plane is oriented along the equator. Stars within 20 kpc of the galactic centrehave been excluded. Left: the projected logarithmic number of BHBs, K giants and MSTO stars. Right: the ratio of K giants/BHBs,MSTO/BHBs and MSTO/K giant stars. are composed of the very youngest stars in the halo. As wellas being a few gigayears younger than the mean, they arealso slightly more metal rich. In these populations BHB starshave not yet had time to form, which explains their deficit.The actual age of these structures may be even younger, butthe average is increased by the other non-associated halostars along the same line-of-sight. Interestingly, there is apatch of intermediate age stars that stands out in the lowerright hand corner of the age map, but is not visible in anyof the tracer maps. This suggests that while differences in the make up of stellar populations will help identify somestructures, not all will be found with this technique. In this paper we have described a new way to construct mockcatalogues of stellar haloes based on cosmological N-bodysimulations. We have adopted the particle tagging method ofCooper et al. (2010) and extended it to turn the output into c (cid:13) , 1–19 aking Mocks of the Stellar Halo Age 10 13 Gyr [Fe / H] -2 -1 Figure 13. All-sky maps of the mean properties of the Aq-C stellar halo when averaged along the line of sight. Left map: average age. Right map: average metallicity map. Stars within 20 kpc of the galactic centre have been excluded. a catalogue of individual stars that can be directly comparedto observations.Using the Aquarius N -body simulations of galactichaloes as the basic framework, we employed the semi-analytical galaxy formation model, galform , to calculatethe evolution of the baryonic component of the Universe.This predicts the amount and properties of the stars form-ing in each halo throughout the simulation. These stellarpopulations were used to tag dark matter particles in the N -body simulation. Finally, these tags are turned into a mockcatalogue using theoretical isochrones to convert the fun-damental properties of the stellar populations, such as ageand metallicities, into a set of observables such as colour andmagnitudes. A phase-space sampling technique was appliedto “explode” the massive particles into numerous individ-ual stars while maintaining the phase-space structure of theoriginal simulation.The output of our method represents, within the limita-tions of the model, a perfect dataset containing precise infor-mation about the entire galaxy. Actual observations usuallyprovide an incomplete, limited view of a galaxy. They maybe subject to selection criteria, be limited in coverage, or,in the case of surveys of the stellar halo, be obstructed byother structures such as the galactic disc. To make a faircomparison between simulations and observations, it is nec-essary to convert the simulation data into a form that alsoreproduces the limitations and biases of the observationaldata. The selection criteria of an actual survey, encoded ina mask or selection function, can be readily applied to ourmocks. Observational errors, foreground dust extinction andcontaminating sources can also be included easily. Thus, inprinciple, the mock catalogues can be viewed and analyzedin a manner that mimics the observational data. In this pa-per, we have used the Aquarius simulations to make mockcatalogues of galactic stellar haloes that can be used to helpinterpret data from upcoming surveys such as those to becarried out by the GAIA mission. However, there is no rea-son why this technique could not be applied to other simula-tions of both smaller haloes, such as dwarfs, or larger haloes,such as clusters.Using our mock catalogues we carried out simple anal- yses to explore whether current observations of the MilkyWay halo provide an accurate picture. We have found that: • The accreted stellar halo is mainly built up from a fewmassive objects. Those that are in early stages of disrup-tion still maintain a coherent structure and are locally con-centrated. They dominate the spherically averaged densityprofile at their corresponding radii. Since the overall densityprofile is so sensitive to these few substructures it provideslimited information about the structure of the underlyingdark matter halo. Instead, it tells us more about the recentaccretion history (Deason et al. 2013). Three of the five haloswe have investigated show clear breaks in their sphericallyaveraged density profiles while two are well described by apower-law. In our model, these breaks have nothing to dowith the distinction between accreted and in-situ compo-nents, as has often been claimed for the Milky Way halo.Rather, they are produced by the accretion of satellites, al-though not necessarily by a single one. • The structure of the accreted stellar halo is clumpy andthus varies considerably over the sky. Attempts to measureproperties of the halo, such as its density profile, based onlimited regions of the sky, provide biased estimates. For ex-ample, if the profile is modeled as a power-law, the slopesmeasured in different pencil-beam surveys through our mockstellar haloes vary from very steep to flat depending on thedirection of the survey. • Beyond 20 kpc, the flattening of dark matter and halostars (as measured by c/a , the major to minor axis ratioof an ellipsoidal fit) is approximately constant out to R .The stellar distribution is significantly flatter than the darkmatter. The minor axis of the stellar distribution is approx-imately aligned with that of the inner dark matter halo,which we have also assumed to be the normal vector to themain symmetry plane. • Differences in the ages and metallicities of the satellitesthat build up the stellar halo are visible as local variationsin the composition of stellar halo populations. These arereflected in differences in the abundance of different stellartypes. For example young, more metal rich objects tend tolack BHB stars. Searching for such fluctuations in halo starsurveys offers a way to identify streams and structure in thehalo. c (cid:13)000 Right map: average metallicity map. Stars within 20 kpc of the galactic centre have been excluded. a catalogue of individual stars that can be directly comparedto observations.Using the Aquarius N -body simulations of galactichaloes as the basic framework, we employed the semi-analytical galaxy formation model, galform , to calculatethe evolution of the baryonic component of the Universe.This predicts the amount and properties of the stars form-ing in each halo throughout the simulation. These stellarpopulations were used to tag dark matter particles in the N -body simulation. Finally, these tags are turned into a mockcatalogue using theoretical isochrones to convert the fun-damental properties of the stellar populations, such as ageand metallicities, into a set of observables such as colour andmagnitudes. A phase-space sampling technique was appliedto “explode” the massive particles into numerous individ-ual stars while maintaining the phase-space structure of theoriginal simulation.The output of our method represents, within the limita-tions of the model, a perfect dataset containing precise infor-mation about the entire galaxy. Actual observations usuallyprovide an incomplete, limited view of a galaxy. They maybe subject to selection criteria, be limited in coverage, or,in the case of surveys of the stellar halo, be obstructed byother structures such as the galactic disc. To make a faircomparison between simulations and observations, it is nec-essary to convert the simulation data into a form that alsoreproduces the limitations and biases of the observationaldata. The selection criteria of an actual survey, encoded ina mask or selection function, can be readily applied to ourmocks. Observational errors, foreground dust extinction andcontaminating sources can also be included easily. Thus, inprinciple, the mock catalogues can be viewed and analyzedin a manner that mimics the observational data. In this pa-per, we have used the Aquarius simulations to make mockcatalogues of galactic stellar haloes that can be used to helpinterpret data from upcoming surveys such as those to becarried out by the GAIA mission. However, there is no rea-son why this technique could not be applied to other simula-tions of both smaller haloes, such as dwarfs, or larger haloes,such as clusters.Using our mock catalogues we carried out simple anal- yses to explore whether current observations of the MilkyWay halo provide an accurate picture. We have found that: • The accreted stellar halo is mainly built up from a fewmassive objects. Those that are in early stages of disrup-tion still maintain a coherent structure and are locally con-centrated. They dominate the spherically averaged densityprofile at their corresponding radii. Since the overall densityprofile is so sensitive to these few substructures it provideslimited information about the structure of the underlyingdark matter halo. Instead, it tells us more about the recentaccretion history (Deason et al. 2013). Three of the five haloswe have investigated show clear breaks in their sphericallyaveraged density profiles while two are well described by apower-law. In our model, these breaks have nothing to dowith the distinction between accreted and in-situ compo-nents, as has often been claimed for the Milky Way halo.Rather, they are produced by the accretion of satellites, al-though not necessarily by a single one. • The structure of the accreted stellar halo is clumpy andthus varies considerably over the sky. Attempts to measureproperties of the halo, such as its density profile, based onlimited regions of the sky, provide biased estimates. For ex-ample, if the profile is modeled as a power-law, the slopesmeasured in different pencil-beam surveys through our mockstellar haloes vary from very steep to flat depending on thedirection of the survey. • Beyond 20 kpc, the flattening of dark matter and halostars (as measured by c/a , the major to minor axis ratioof an ellipsoidal fit) is approximately constant out to R .The stellar distribution is significantly flatter than the darkmatter. The minor axis of the stellar distribution is approx-imately aligned with that of the inner dark matter halo,which we have also assumed to be the normal vector to themain symmetry plane. • Differences in the ages and metallicities of the satellitesthat build up the stellar halo are visible as local variationsin the composition of stellar halo populations. These arereflected in differences in the abundance of different stellartypes. For example young, more metal rich objects tend tolack BHB stars. Searching for such fluctuations in halo starsurveys offers a way to identify streams and structure in thehalo. c (cid:13)000 , 1–19 Ben Lowing et al In this paper we have developed a method of convert-ing massive stellar particles into stars based on the outputof the Cooper et al. (2010) N -body particle tagging method.However, with only minor modification it should be possi-ble to use the same conversion technique on hydrodynamicalgalaxy formation simulations. By combining these two ap-proaches, components that are missing from our treatment,such as the galactic disc or bulge, could be included, whileother components, such as the stellar halo, could be followedwith much higher resolution than is attainable with currenthydrodynamic simulations.The five complete mock catalogues for theAquarius haloes A-E, are publicly available at http://galaxy-catalogue.dur.ac.uk:8080/StellarHalo ,in a database that can by queried via SQL. The propertiesthat are available in each catalogue are described in theTable 1. ACKNOWLEDGMENTS We thank the anonymous referee for their careful readingof our paper and helpful suggestions. WW is grateful foruseful discussions with Jiaxin Han and Yanchuan Cai. Thiswork was supported by the European Research Council [GA267291] COSMIWAY and Science and Technology FacilitiesCouncil [ST/F001166/1,ST/L00075X/1]. WW is supportedby the Durham International Junior Research Fellowship[RF040353]. APC acknowledges a Chinese Academy of Sci-ences International Research Fellowship and NSFC grantno. 11350110323. The simulations for the Aquarius Projectwere carried out at the Leibniz Computing Centre, Garch-ing, Germany, at the Computing Centre of the Max-Planck-Society in Garching, at the Institute for Computational Cos-mology in Durham, and on the STELLA supercomputerof the LOFAR experiment at the University of Groningen.Sky maps were produced with the healpy implementationof the healpix algorithms ( http://healpix.jpl.nasa.gov ;G´orski et al. 2005). REFERENCES Abadi, M. G., Navarro, J. F., & Steinmetz, M. 2006, MN-RAS, 365, 747Akhter S., Da Costa G. S., Keller S. C., Schmidt B. P.,2012, ApJ, 756, 23Arad I., Dekel A., Klypin A., 2004, MNRAS, 353, 15Bailin, J., Bell, E. F., Valluri, M., et al. 2014, ApJ, 783, 95Baugh C. M., 2006, Reports on Progress in Physics, 69,3101Beers, T. C., Carollo, D., Ivezi´c, ˇZ., et al. 2012, ApJ, 746,34Bell E. F., Xue X. X., Rix H.-W., Ruhland C., Hogg D. W.,2010, AJ, 140, 1850Bower R. G., Benson A. J., Malbon R., Helly J. C., FrenkC. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS,370, 645Bressan A., Marigo P., Girardi L., Salasnich B., Dal CeroC., Rubele S., Nanni A., 2012, MNRAS, 427, 127Bullock J. S., Johnston K. V., 2005, ApJ, 635, 931 Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2001,ApJ, 548, 33Carollo D., et al., 2007, Nature, 450, 1020Carollo D., et al., 2010, ApJ, 712, 692Chou, M.-Y., Majewski, S. R., Cunha, K., et al. 2007, ApJ,670, 346Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000,MNRAS, 319, 168Colless M., et al., 2001, MNRAS, 328, 1039Cooper A. P., et al., 2010, MNRAS, 406, 744Cooper A. P., Cole S., Frenk C. S., Helmi A., 2011, MN-RAS, 417, 2206Cooper, A. P., D’Souza, R., Kauffmann, G., et al. 2013,MNRAS, 434, 3348Cooper A. P., Gao L., Guo Q., Frenk C. S., Jenkins A.,Springel V., White S. D. M., 2014, arXiv, arXiv:1407.5627Crain R. A., et al., 2009, MNRAS, 399, 1773Chiba M. & Beers T. C., 2000, AJ, 119..2843CCarollo D., Beers T. C., Chiba M., Norris J. E., FreemanK. C., Lee Y. S., Ivezi´c ˇZ., Rockosi C. M., Yanny B., 2010,ApJ, 712..692CDe Lucia G. & Helmi A., 2008, MNRAS, 391..14Dde Bruijne J. H. J., 2012, Ap&SS, 341, 31De Propris R., Harrison C. D., Mares P. J., 2010, ApJ, 719,1582Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS,416, 2903Deason A. J., Belokurov V., Evans N. W., Johnston K. V.,2013, ApJ, 763, 113D’Onghia E., Springel V., Hernquist L., Keres D., 2010,ApJ, 709, 1138Font A. S., Johnston K. V., Bullock J. S., Robertson B. E.,2006, ApJ, 638..585FFont A. S., McCarthy I. G., Crain R. A., Theuns T., SchayeJ., Wiersma R. P. C., Dalla Vecchia C., 2011a, MNRAS,416, 2802Font A. S., et al., 2011b, MNRAS, 417, 1260Gao L., Navarro J. F., Cole S., Frenk C. S., White S. D. M.,Springel V., Jenkins A., Neto A. F., 2008, MNRAS, 387,536G´orski K. M., Hivon E., Banday A. J., Wandelt B. D.,Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ,622, 759Helmi A., 2008, A&ARv, 15, 145Ibata R., Martin N. F., Irwin M., Chapman S., FergusonA. M. N., Lewis G. F., McConnachie A. W., 2007, ApJ,671..1591Ibata R. A., Lewis G. F., McConnachie A. W., MartinN. F., Irwin M. J., Ferguson A. M. N., Babul A., BernardE. J., Chapman S. C., Collins M., Fardal M., MackeyA. D., Navarro J., Pe˜narrubia J., Rich R. M., Tanvir N.,Widrow L., 2014, ApJ, 780..128IJordi C., Gebran M., Carrasco J. M., de Bruijne J., VossH., Fabricius C., Knude J., Vallenari A., Kohley R., MoraA., 2010, A&A, 523, A48Johnston K. V., Bullock J. S., Sharma S., Font A., Robert-son B. E., Leitner S. N., 2008, ApJ, 689..936JJuri´c M., et al., 2008, ApJ, 673, 864Kaiser N., et al., 2010, in Stepp L.M., Filmozzi R., HallH.I., eds, SPIE Conf. Ser., 7733, 14Kafle P. R., Sharma S., Lewis G. F., Bland-Hawthorn J.,F., 2013, MNRAS, 430..2973K c (cid:13) , 1–19 aking Mocks of the Stellar Halo Kazantzidis S., Kravtsov A. V., Zentner A. R., Allgood B.,Nagai D., Moore B., 2004, ApJ, 611, L73Keller S. C., Murphy S., Prior S., Da Costa G., SchmidtB., 2008, ApJ, 678, 851Kennicutt Jr. R. C., 1983, ApJ, 272, 54Libeskind N. I., Knebe A., Hoffman Y., Gottl¨ober S., YepesG., 2011, MNRAS, 418..336LMiceli A., Rest A., Stubbs C. W., Hawley S. L., Cook K. H.,Magnier E. A., Krisciunas K., Bowell E., Koehn B., 2008,ApJ, 678..865MMcConnachie A. W., Irwin M. J., Ibata R. A., DubinskiJ., Widrow L. M., Martin N. F., Cˆot´e P., Dotter A. L., etal., 2009, Nature, 461..66MMorrison H. L., Flynn C., Freeman K. C., 1990, AJ, 100,1191Morrison H. L., Mateo M., Olszewski E. W., Harding P.,Dohm-Palmer R. C., Freeman K. C., Norris J. E., MoritaM., 2000, AJ, 119, 2254Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS,283L..72NNavarro J. F., et al.,2010, MNRAS, 402,21Pontzen A. & Governato F., 2013, MNRAS, 430..121PPasetto S., Chiosi C., Kawata D., 2012, A&A, 545, A14Pillepich, A., Vogelsberger, M., Deason, A., et al. 2014,arXiv:1406.1174Pe˜narrubia J., Benson A. J., Walker M. G., Gilmore G.,McConnachie A. W., Mayer L., 2010, MNRAS, 406, 1290Prusti T., 2012, Astronomische Nachrichten, 333, 453Robertson B., Bullock J. S., Font A. S., Johnston K. V.,Hernquist L., 2005, ApJ, 632..872RSawala T., Frenk C. S., Crain R. A., Jenkins A., Schaye J.,Theuns T., Zavala J., 2013, MNRAS, 431, 1366Sch¨onrich R., Binney J., Dehnen W., 2010, MNRAS, 403,1829Sch¨onrich, R., Asplund, M., & Casagrande, L. 2014, ApJ,786, 7Sesar B., et al.,2010, ApJ, 708, 717Sesar B., Juri´c M., Ivezi´c ˇZ., 2011, ApJ, 731, 4Sharma S., Bland-Hawthorn J., Johnston K. V., Binney J.,2011, ApJ, 730, 3Sharma S., Steinmetz M., 2006, MNRAS, 373, 1293Spergel D. N., et al., 2003, ApJS, 148, 175Springel V., et al., 2008a, MNRAS, 391, 1685Springel V., et al., 2008b, Nature, 456, 73Starkenburg E., et al., 2009, ApJ, 698, 567Tissera P. B., White S. D. M., Scannapieco C., 2012, MN-RAS, 420, 255Tumlinson, J. 2010, ApJ, 708, 139Tissera P. B., White S. D. M., Scannapieco C., 2012, MN-RAS, 420..255TTissera P. B., Scannapieco C., Beers T. C., Carollo D.,2013, MNRAS, 432..3391TVogelsberger M., et al., 2014, Nature, 509, 177Walker M. G., Mateo M., Olszewski E. W., Gnedin O. Y.,Wang X., Sen B., Woodroofe M., 2007, ApJ, 667, L53Watkins L. L., et al., 2009, MNRAS, 398, 1757Xue X.-X., et al., 2014, ApJ, 784, 170Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ,137, 4377York D. G., et al., 2000, AJ, 120, 1579Zolotov A., Willman B., Brooks A. M., Governato F., HoggD. W., Shen S., Wadsley J., 2010, ApJ, 721, 738 Mart´ınez-Delgado, D., Gabany, R. J., Crawford, K., et al.2010, AJ, 140, 962 c (cid:13)000