PreprinttypesetusingLATEXstyleemulateapjv.05/04/06 DARK MATTER SUBSTRUCTURE AND GAMMA-RAY ANNIHILATION IN THE MILKY WAY HALO Ju¨rg Diemand1,2, Michael Kuhlen1,3, & Piero Madau1,4 ABSTRACT We present initial results from “Via Lactea”, the highest resolution simulation to date of Galactic CDM substructure. It follows the formation of a Milky Way-size halo with M = 1.8×1012 M halo ⊙ in a WMAP 3-year cosmology, using 234 million particles. Over 10,000 subhalos can be identified at z=0: Their cumulative mass function is well-fit by N(> M ) = 0.0064(M /M )−1 down sub sub halo to M = 4×106 M . The total mass fraction in subhalos is 5.3%, while the fraction of surface sub ⊙ mass density in substructure within a projected distance of 10 kpc from the halo center is 0.3%. Because of the significant contribution from the smallest resolved subhalos, these fractions have not converged yet. Sub-substructure is apparent in all the larger satellites, and a few dark matter lumps areresolvedeveninthesolarvicinity. Thenumberofdarksatelliteswithpeakcircularvelocitiesabove 10 kms−1 (5 kms−1) is 124 (812): of these, 5 (26) are found within 0.1r , a region that appeared vir practicallysmoothinprevioussimulations. The neutralino self-annihilationγ-rayemissionfromdark matter clumps is approximately constant per subhalo mass decade. Therefore, while in our run the contribution of substructure to the γ-ray luminosity of the Galactic halo amounts to only 40% of the total spherically-averagedsmooth signal, we expect this fraction to grow significantly as resolution is increasedfurther. An all-sky map of the expected annihilation γ-ray flux reaching a fiducial observer at8kpcfromtheGalacticcentershowsthatatthecurrentresolutionasmallnumberofsubhalosstart to be bright enough to be visible against the background from the smooth density field surrounding the observer. Subject headings: cosmology: theory–darkmatter–galaxies: dwarfs–galaxies: formation–galaxies: halos – methods: numerical 1. INTRODUCTION parablecircularvelocity(Mooreetal.1999;Klypinetal. 1999) has been the subject of many recent studies. So- In structure formation scenarios dominated by cold lutions involving feedback mechanisms that make star dark matter (CDM), the halos of galaxies and galaxy formation in small halos very inefficient offer a possi- clusters are assembled via the hierarchical merging and ble way out (e.g. Bullock et al. 2000; Somerville 2002; accretion of smaller progenitors. This process causes Kravtsov et al. 2004; Moore et al. 2006). Other models structures to virialize to a new equilibrium by redis- have attempted to solve the apparent small-scale prob- tributing energy among the collisionless mass compo- lems of CDM at a more fundamental level, i.e. by re- nents. Early low resolution numerical simulations and ducing small-scale power (e.g. Kamionkowski & Liddle simple analytical models found that the end products of 2000). Even if most dark matter satellites have no op- this “bottom-up” scenario were smooth, triaxial halos. tically luminous counterparts, the substructure popula- In recent years, higher resolution cosmological N-body tionmaybedetectable viafluxratioanomaliesinstrong simulations have modified this picture: the merging of gravitational lenses (e.g. Metcalf & Madau 2001; Chiba progenitors is not always complete, and the cores of ac- 2002),throughitseffectsonstellarstreams(e.g. Ibataet creted halos often survive as gravitationally bound sub- al. 2002; Mayer et al. 2002), or possibly via γ-rays from halos orbiting within a larger host system. CDM halos darkmatter annihilationin their cores(Bergstromet al. arenotsmooth,theyhaveawealthofsubstructureonall 1999;Calcaneo-Roldan&Moore2000;Stoehretal.2003; resolved mass scales (e.g. Moore et al. 1999; Klypin et Colafrancesco et al. 2006; Diemand et al. 2006). al. 1999;Moore et al. 2001;Stoehr et al. 2003; Diemand The possibility of observing the fingerprints of the et al. 2004b;Reed et al. 2005;Gao et al. 2005;Diemand small-scale structure of CDM hinges on the ability of et al. 2006). subhalostosurvivethe hierarchicalclustering processas The amount and spatial distribution of subhalos substructurewithinthehost. Thisisinturnparticularly around their host provide unique information and clues sensitive to resolution issues, as subhalos with numeri- on the galaxy assembly process and the nature of the callysoftened centraldensities are moreeasily disrupted dark matter. The mismatch between the dozen or so bytidalforces. Inthispaperwepresentinitialresultson dwarf satellite galaxies observed around the Milky Way halo substructure from a new N-body cosmological sim- andthepredictedlargenumberofCDMsubhalosofcom- ulation of unprecedented dynamic range: dubbed “Via 1Department of Astronomy & Astrophysics, University of Lactea”,itresolvesthe galaxyformingregionofa Milky California,SantaCruz,CA95064. Way-size halo at z = 0 with over 200 million particles. 2HubbleFellow. This is an order of magnitude more than used in previ- 3School of Natural Sciences, Institute for Advanced Study, oussimulations. Therunwascompletedin320,000CPU EinsteinDrive,Princeton,NJ08540. 4Max-Planck-Institut fuer Astrophysik, Karl-Schwarzschild- hoursonNASA’sProjectColumbiasupercomputer,cur- Strasse1,85740GarchingbeiMuenchen, Germany. rently one of the fastest machines available. Electronic address: diemand@ucolick.org, mqk@ias.edu, pmadau@ucolick.org. 2. SIMULATIONS 2 Dark matter substructure in the MW Fig. 1.—Top: logarithmicslopeofthedensityprofileofourViaLactearun,asafunctionofradius. Densitieswerecomputedin50radial logarithmicbins,andthelocalslopewasdeterminedbyafinitedifferenceapproximationusingoneneighboringbinoneachside. Thethin lineshowstheslopeofthebest-fitNFWprofilewithscaleradiusof24.6kpc. Theverticaldottedlineindicatestheestimatedconvergence radius: local densities (but not necessarily the logarithmic slopes) should be correct to within 10% outside of this radius. Bottom: the residualsinpercentbetween thedensityprofileandthebest-fitNFWprofile,asafunctionofradius. The simulation was performed with PKDGRAV ical parameters from the WMAP three-year data re- (Stadel 2001; Wadsley, Stadel & Quinn 2004) and em- lease (Spergel et al. 2006): Ω = 0.238, Ω = 0.762, M Λ ployedmultiple mass particle gridinitial conditions gen- H = 73 kms−1Mpc−1, n = 0.951, and σ = 0.74. The 0 8 erated with the GRAFICS package (Bertschinger 2001). lower values of σ and n compared to 1-year WMAP 8 The high resolution region was sampled with 234 mil- results have the effect of delaying structure formation lionparticles particles of mass 2.1×104M and evolved andreducingsmall-scalepower. Thesimulationwascen- ⊙ with a force resolution of 90 pc. It was embedded teredonanisolatedhalothathadnomajormergerafter within a periodic box of comoving size 90 Mpc, which z = 1.7, making it a suitable host for a Milky Way-like was sampled at lower resolution to account for the large disk galaxy. In this work we focus on results obtained scale tidal forces. We adopted the best-fit cosmolog- atz=0,theformationhistoryoftheViaLacteahaloand Diemand, Kuhlen, & Madau 3 zi zf ǫ Np mp solved. 48.4 0.0 90pc 2.34 108 2.09 104M⊙ 48.4 0.0 378pc 6.50×107 5.64×105M⊙ 3.1. Subhalo identification × × We have constructed substructure catalogs using Note. —Initialandfinalredshiftszi andzf,(spline)softening the phase-space friends-of-friends algorithm 6DFOF de- length ǫ, total number Np and mass mp of dark matter particles scribedinDiemandetal.(2006). Weadoptaspacelink- for the Via Lactea simulation (top) and a lower resolution run (bottom). Force softening lengths ǫ are constant in physical units inglengthof0.28physicalkpc(0.033ofthemeanparticle backtoz=9andconstant incomovingunitsbefore. separationatz =0),avelocityspacelinkinglengthof16 physical kms−1 (0.1 of the mean 1D velocity dispersion itssubstructureswillbepresentedinasubsequentpaper. of the host halo), and keep all groups with at least 16 The host halo mass at z =0 is M =1.77×1012 M halo ⊙ members. Around their centers, spherical density pro- withinaradiusofr =389 kpc(Wedefine”r ”asthe vir vir files (using 50 logarithmic bins) were constructed: the radius within which the enclosed average density is 200 bins reach out to r in the case of isolated halos or out times the mean matter value). We used adaptive time- vir to the center of the nearest larger halo in the case of stepsasshortas(13.7/2×105)Gyr=68,500yrsatisfying subhalos, i.e. all structures are allowed to contain their the condition ∆t<0.2 ǫ/a, where ǫ is the force soften- own,smallersubstructures. The resultingcircularveloc- ing length and a is theplocal acceleration. This criterion ityprofileswerethenfittedwiththesumofanNFWpro- isinsufficienttoensureconvergenceatsmallradii,where fileandaconstantdensitybackground. Theradialrange the dynamical time becomes too short and runs using used in these fits starts at 3 force softening lengths, i. e. the ǫ/a criterion start to underestimate matter densi- 0.27 kpc in the Via Lactea run, and ends either at the ties.pIn comparable runs this problem was observed to largestradialbinstoredduringatruntimeorfurtherinif togive10%toolowlocaldensitiesat6×105ρcrit(z =0), thecircularvelocitygrowslinearly(within10%)overfive see Figure 1 of Diemand et al. 2005. This density is consecutive bins. Figure 4 shows nine example subhalo reached at 0.0034rvir in the Via Lactea run, therefore circular velocity profile fits. A subhalo tidal radius rt we set rconv = 0.0034rvir. The local densities (but not wasdefinedasρsub(rt)=2×ρBG,whereρBG isthelocal necessarily the logarithmic slopes) should be correct to matter density of the host halo. This is the tidal radius within 10% beyond 0.0034rvir, but detailed convergence ofanisothermal(ρsub ∝r−2)satelliteonacircularorbit tests using galaxy scale halos would be needed to con- within an isothermal host. We define the subhalo mass firm this. To resolve the density profile further in than to be all the mass within r , i.e. we did not perform an t 0.0034rvir, adaptive time-steps that scale with dynam- unbinding procedure. The speed of 6DFOF and its par- ical timescale should be adopted instead (Zemp et al. allel implementation into PKDGRAV allowed us to use 2006). TheforcecalculationsinthetreealgorithmPKD- it at run-time with negligible computational cost (com- GRAV includes multipole moments up to hexa-decapole pared to the gravity calculations) on all 200 snapshots toreachhighforceaccuracy. Tocheckfornumericalcon- that were stored (separated by 68.5 Myr). vergence of substructure properties we resimulated the The6DFOFsubhalolistwascheckedinvariousregions same halo (plus a few neighboring galaxyhalos) with 27 oftheVia LacteahaloagainstalistobtainedwithSKID timesmoremassiveparticles. Theparametersofthetwo (Stadel 2001). The results compare very well down to simulations are summarized in Table 2. about 40 particles per halo except for the very largest Figure 1 (upper panel) shows the logarithmic slope of subhalos, where the SKID algorithm tends to underesti- the spherically-averaged density profile, dlogρ/dlogr, mate the masses and is unable to properly account for for our Via Lactea halo, plotted versus radius. The the presence of resolved sub-substructure. Detailed evo- fitting formula proposed by Navarro, Frenk, & White lutionarytracksofsubhalopropertiesanddensityprofiles (1997) with scale radius of 24.6 kpc provides a reason- will be published in a subsequent paper. able approximation to the density profile down to our 3.2. Subhalo abundance convergence radius r . Within the region of conver- conv gence, deviations from the best-fit NFW matter density In Figure 5 we present the cumulative mass func- are typically less than 10%. From 0.025rvir down to tion for our entire subhalo populationw ithin rvir as rconv ViaLacteaisactuallydenserthanpredictedbythe well as for the subpopulation within the inner 0.1rvir. NFW formula. Near rconv the density approaches the The former is well-fit by a power-law N(> Msub) = NFWvalueagainwhilethelogarithmicslopeisshallower 0.0064(M /M )−1 above M = 4 × 106 M , as sub halo sub ⊙ (−1.0 at rconv) than that the NFW fit. While we have shown by the solid line (formally, the best-fit slope checked that these conclusions are robust to variations in the mass range M = 4 × 106 − 4 × 109 M is sub ⊙ inthe binning usedtoconstructthe densityprofile,they −0.97 ± 0.03). This slope implies an equal mass per might be affected by numerical flattening near our con- subhalomassdecade,makingthetotalmassinsubstruc- vergence radius. We plan to run new simulations with ture quite dependent on resolution. Note how numeri- improved time-stepping (Zemp et al. 2006) to address cal resolution effects start to flatten the distribution al- this issue. readybelowmassescorrespondingtoabout200particles (i.e. below M = 4 × 106 M and M = 108 M 3. SUBSTRUCTUREINGALAXYHALOS for the Via Lacsuteba and the lowe⊙r resolutiosnubrun, respec⊙- The wealth of substructure resolvedin our Via Lactea tively). The fraction of the host halo mass in subha- run is clearly see in Figure 2. About 10,000 surviving los is f = 5.3%(0.87%) within r (0.1r ). The sub vir vir satellites, some sub-subhalos (a few examples are shown contribution from subhalos below 108 M is significant: ⊙ in Fig. 3), and even a couple of dark matter clumps at f (M < 108 M ) = 2.3%(0.15%). The total sub- sub sub ⊙ less than 8 kpc from the Galactic Center are now re- structure mass fractions in the lower resolution run are 4 Dark matter substructure in the MW Fig. 2.— Projected dark matter density-squared map of our simulated MilkyWay-size halo (“Via Lactea”) at the present epoch. The imagecoversanareaof800 600kpc,andtheprojectiongoesthrougha600kpc-deepcuboidcontainingatotal of110millionparticles. Thelogarithmiccolorscalec×overs20decades indensity-square. Diemand, Kuhlen, & Madau 5 Fig. 3.—Projecteddarkmatterdensity-squaremapofthefourmostmassivesubhaloswithinthesimulatedMilkyWayhostatthepresent epoch. Sub-substructure is clearly visible. Only dark matter particles within the tidal radius rt are used for the projections. Clockwise (fr2o.4m to1p09lefMt:⊙(,M14su.7bk,rptc,,r6V.m1akxp)c=).(9T.h8e×m10ea9nMs⊙ub,4h0a.l1okdpecn,s7i.t6iekspcw)i,t(h3i.n7×th1e0t9idMal⊙r,a3d3i.u4skp(icn,4u.0niktpsco),f(t3h.0e×co1s0m9icMb⊙a,c2k8g.r0okupncd,4d.9arkkpcm),atatnedr densi×ty)are1002, 654, 904, and4950, respectively. Thesevalues arerelatedtothelocal matterdensityofthehost(72, 46,59and397in thesameunits),andcorrelateonlyweaklywiththesubhalodistancefromtheGalacticcenter (345,374,280and185kpc). 6 Dark matter substructure in the MW R=239 kpc, R=2.6 kpc, M=2.3e+06 M , r = 20 r R=150 kpc, R=2.3 kpc, M=8.2e+06 M , r =176 r R= 30 kpc, R=0.6 kpc, M=4.3e+06 M , r =3929 r t t SUN BG crit t t SUN BG crit 16 t t SUN BG crit 3.5 7 14 3 6 12 m/s]2.5 m/s] 5 m/s] 10 [kvcirc1.52 [kvcirc 34 [kvcirc 68 1 2 4 0.5 1 2 00 2 4 6 8 10 12 14 00 1 2 3 4 5 6 7 8 9 00 0.5 1 1.5 2 2.5 3 3.5 4 4.5 r [kpc] r [kpc] r [kpc] R= 15 kpc, R=0.4 kpc, M=5.8e+06 M , r =19892 r R=132 kpc, R=2.6 kpc, M=2.0e+07 M , r =154 r R=238 kpc, R=11.6 kpc, M=3.2e+08 M , r = 27 r t t SUN BG crit t t SUN BG crit t t SUN BG crit 22 10 18 20 18 8 16 16 14 m/s] 14 m/s] 6 m/s] 12 [kvcirc11802 [kvcirc 4 [kvcirc 180 6 6 4 2 4 2 2 00 0.5 1 1.5 2 2.5 00 2 4 6 8 10 12 00 10 20 30 40 50 60 r [kpc] r [kpc] r [kpc] R=318 kpc, R=14.4 kpc, M=6.3e+08 M , r = 33 r R=279 kpc, R=14.6 kpc, M=8.7e+08 M , r = 31 r R=211 kpc, R=14.4 kpc, M=1.3e+09 M , r = 51 r t t SUN BG crit t t SUN BG crit 30 t t SUN BG crit 20 14 18 25 12 16 m/s] 10 m/s] 1124 m/s] 20 [kvcirc 68 [kvcirc 180 [kvcirc 15 10 4 6 4 5 2 2 00 2 4 6 8 10 12 14 16 18 20 00 10 20 30 40 50 60 00 10 20 30 40 50 60 70 r [kpc] r [kpc] r [kpc] Fig. 4.— Examples of subhalo circular velocity profiles (crosses). The fitting functions (dashed lines) are combinations of an NFW subhalo(dottedlines)andalinearcontributionformaconstantbackgrounddensity(solidlines). Statedaboveeachpanelarethesubhalos distanceformthegalaxycenter, tidalradius,tidalmassandthelocalbackgrounddensity. substantially smaller, f = 2.6%(0.11%). This is in nitude larger than the dozen or so dwarf satellite galax- sub part because the low resolution realization happens to ies of comparable circular velocity observed around the contain slightly fewer massive satellites at z = 0, and MilkyWay(Mooreetal.1999;Klypinetal.1999),awell in part because of the reduced number of resolved small known mismatch often referred to as the “missing satel- subhalos. lites problem”. The WMAP 3-year cosmologicalparam- Figure 6 shows the cumulative peak circular velocity eters do not alleviate the discrepancy significantly, and (V ≡ max{ GM ( V ) = (1/48)(V )/V )−3 above nounced. Indeed, we find a new “missing inner satellites max max max,host 5 kms−1, as shown by the dotted line. This is in good problem”within0.1rvir =39kpc,aregionthatappeared agreement with the mean abundance obtained from a practically smooth in previous lower resolution simula- largesample ofhosthalosby Reedetal.(2005)after ac- tions. In this inner region we now resolve 5 (26) dark counting for the effect of the different cosmology.5 Note satellites with Vmax > 10 kms−1 (Vmax > 5 kms−1), how numerical effects start to flatten the velocity func- compared to only one known Milky Way dwarf galaxy, tions below 5 kms−1 and already below 15 kms−1 for Sagittariusatadistanceof24±2kpc(Mateo1998). Even the lower resolution run. when the suspected stripped cores of dwarf spheroidals fromtherecentcompilationofvandenBergh (2006)are included,thenumbersofobservedinnersatellitesremain 3.3. A ”missing inner satellites” problem? well below the predicted number of substructures in the The number of dark matter subhalos with peak circu- inner halo. lar velocity larger than 10 kms−1 is 124 (compared to Recent Local Group models (Kravtsov et al. 2004; 74 in the lower resolutionrun). This is an order of mag- Moore et al. 2006) were based on simulations like the lowerresolutionrunusedinthiswork. Theoverallabun- 5Reedetal.(2005)simulationswerenormalizedusingσ8=1.0. dance of subhalos does increase with resolution in the Decreasing σ8 to 0.74 (the value used here) lowers the subhalo relevant V range, but not dramatically, i.e. the to- velocity function by about a factor two (see Zentner & Bullock max 2003andtheerratumtoDiemandetal.2004b). talabundance of satellite galaxiesinthese models seems Diemand, Kuhlen, & Madau 7 ) 110034 b u Ms > ( 102 N 10 1 106 107 108 109 1010 M [M ] sub Fig. 5.— Cumulative subhalo mass functions within rvir (upper curves) and 0.1rvir (lower curves). Solid crosses: Via Lactea run. Dashed crosses: lowerresolutionrun. Solid line: power-lawfit,N(>Msub)=0.0064 (Msub/Mhalo)−1. × not to be very affected by their lower numerical resolu- imagedquasars(Mao&Schneider1998;Metcalf&Zhao tion. But what would these two models predict about 2002). The discrepancy between the predicted and ob- satellite galaxies within the newly resolved inner sub- servedfluxratiosiscommonlyreferredtoasthe“anoma- halos? The Kravtsov et al. (2004) model allows galaxy lous flux ratio problem”, and dark matter substructure formation only in subhalos that were relatively massive within the lens halo is one of the leading interpreta- (>109 M )beforetheyfellintothegravitationalpoten- tions of such anomaly (e.g. Metcalf & Madau 2001; ⊙ tial of the host. In the Moore et al. (2006) scenario only Chiba2002;Dalal& Kochanek2002;Bradacetal.2004; the earliest forminghalosabovethe atomiccoolingmass Mao et al. 2004; Amara et al. 2006; Maccio et al. 2006). at z > 12 become proto-galactic building blocks (many Ratherthanthetotalmassfraction,lensingobservations ofthesesystemsremainwellbelow109 M atalltimes). are sensitive to the mass fraction in substructure pro- ⊙ Byz =0theyhavebuiltupahaloofstars,globularclus- jected through a cylinder of radius 5-10 kpc around the ters and a few surviving dwarfs. We followed the assem- lens center. Figure 7 depicts the fraction of the host bly history of the inner subhalos backwards in time and halo mass within a sphere of radius r that is bound found that only two of them lie above the (time depen- up in substructure, f (< r), as well as the substruc- sub dent) minimum mass from Kravtsov et al. (2004). The ture mass fraction in cylindrical projection of radius R, same two systems are also the only ones that form early f ( ( N 10 1 1 2 3 4 5 6 7 8 910 20 30 40 50 V [km/s] max 102 ) x a m V > ( 10 N 1 1 2 3 4 5 6 7 8 910 20 30 40 50 V [km/s] max Fig. 6.— Cumulative peak circular velocity function for all subhalos within rvir (upper panel) and for the subpopulation within the inner 0.1rvir (lower panel). Solid lines without error bars: ViaLactea run. Dashed lines: lower resolution run. Solid line: power-law fit N(>Vmax)=(1/48)×(Vmax/Vmax,host)−3. Solidlineswitherrorbars: observednumberofdwarfgalaxysatellitesaroundtheMilkyWay (asinMateo(1998);Klypinetal.(1999),plustherecentlydiscoveredUrsaMajordwarf(Willmanetal. 2005;Kleynaetal.2005)),within rvir andwithin0.1rvir (onlySagittarius). TheverticalerrorbarsshowthePoissonian√N scatter. Diemand, Kuhlen, & Madau 9 1100--11 )) 1100--22 RR << (( BB UU SS f f 1100--33 ,, )) rr << (( BB UU SS ff 1100--44 --55 1100 1100 110022 rr [[kkppcc]] ,, RR [[kkppcc]] Fig. 7.—Cumulativemassfractionofsubstructure,shownspherically-averagedasafunctionofradiusr(dashed lines)andincylindrical projectionasafunctionofprojectedradiusR(solidlines). Averagesover1000randomprojectionsareshown,andonlysubtructureswithin rvir areincluded. Thick lines: allsubhalos. Thin lines: subhaloswithMsub<(107 M⊙,108 M⊙,109 M⊙)frombottom totop. 10 Dark matter substructure in the MW we measure, dN/dM ∝ M−2, extended all the way tion indicates that more than 500 dark matter particle sub sub down to 104 M , for example, then the mass fraction per subhaloare neededto produce a signalthat is of the ⊙ in subhalos with 104 M < M < 4×106 M would right order of magnitude. Below this scale (i.e. below ⊙ sub ⊙ be comparable to that in all higher mass satellites with about107 M⊙inViaLacteaand3×108 M⊙forthelower 4×106

DARK MATTER SUBSTRUCTURE AND GAMMA-RAY ANNIHILATION IN THE MILKY WAY HALO ...

Similar DARK MATTER SUBSTRUCTURE AND GAMMA-RAY ANNIHILATION IN THE MILKY WAY HALO ...