Table Of Content

clumpy: a code for γ-ray signals from dark matter structures AldeéCharbonniera,Ce´lineCombetb,c,DavidMaurina,b,c,d aLaboratoiredePhysiqueNucleáireetHautesEnergies,CNRS-IN2P3/Universite´sParisVIetParisVII,4placeJussieu,Tour33,75252ParisCedex05,France bDept.ofPhysicsandAstronomy,UniversityofLeicester,Leicester,LE17RH,UK cLaboratoiredePhysiqueSubatomiqueetdeCosmologie,CNRS/IN2P3/INPG/Universite´JosephFourierGrenoble1,53avenuedesMartyrs,38026Grenoble, France dInstitutd’AstrophysiquedeParis,UMR7095CNRS,Universite´PierreetMarieCurie,98bisbdArago,75014Paris,France 2 1Abstract 0 We presentthefirstpubliccodeforsemi-analyticalcalculationoftheγ-rayfluxastrophysicalJ-factorfromdarkmatterannihila- 2 tion/decayintheGalaxy,includingdarkmattersubstructures. Thecoreofthecodeisthecalculationofthelineofsightintegralof n thedarkmatterdensitysquared(forannihilations)ordensity(fordecayingdarkmatter). Thecodecanbeusedinthreemodes: i) a JtodrawskymapsfromtheGalacticsmoothcomponentand/orthesubstructurecontributions,ii)tocalculatethefluxfromaspecific halo(thatisnottheGalactichalo,e.g. dwarfspheroidalgalaxies)oriii)toperformsimplestatisticaloperationsfromalistofal- 3 2lowedDMprofilesforagivenobject.ExtragalacticcontributionsandothertracersofDMannihilation(e.g.positrons,antiprotons) willbeincludedinasecondrelease. ] EKeywords: Cosmicrays,Cosmology,DarkMatter,Indirectdetection,Gamma-rays H . h PROGRAMSUMMARY 1. Introduction p ProgramTitle:CLUMPY - oProgramminglanguage: C/C++ rComputer:PCandMac t Despite the several astrophysical evidences pointing at the sOperatingsystem:UNIX(Linux),MacOSX existenceofdarkmatter(flatrotationcurvesofgalaxies,gravi- aRAM: depends on the requested size of skymaps ( 40 Mb for a [ ∼ tationallensing,“bulletcluster”,etc.),itsnaturestillevadesus. 500 500map) × Thissearchhasbecomeoneofthemajortopicinbothparticle 1Keywords: darkmatter,indirectdetection,gamma-rays physicsandastrophysicsandistackledusingeitherdirectorin- vClassification: 1.1,1.9 8Externalroutines/libraries:CERNROOT1,Doxygen2(optional) directdetectionmethods.Fortheformer,thehopeistodirectly 2Natureofproblem: Calculationofγ-raysignalfromdarkmatteran- witness the interaction of a dark matter particle with a detec- 7nihilation(resp. decay). Thisinvolvesaparticlephysicstermandan tor. Theindirectapproachesaimatmeasuringtheendproducts 4 astrophysicalone.Thefocushereisonthelatter. ofdarkmatterannihilation/decay(e.g.,γ-rays,positrons,anti- . 1Solution method: Integrationof theDMdensity squared (resp. den- protons). Thedetectionofγ-rayswassoonrecognisedtobea 0sity)alongalineofsight. ThecodeisoptimisedtodealwiththeDM promisingchannel[1,2]. Ifsuchasignalisyettobemeasured 2densitypeaksencounteredalongthelineofsight(DMsubstructures). bythe existingγ-raysobservatories(FERMI,HESS, MAGIC, 1Asemi-analyticalapproach(calibratedonN-bodysimulations)isused : VERITAS), the prospect may significantly improve with the vforthespatialandmassdistributionsofthedarkmattersubstructures forthcomingnext-generationinstruments,suchasCTA[3]. iintheGalaxy. XRestrictions: Some generic dark matter annihilation spectra arepro- Having a modelling tool that can compute the expected γ- rvidedbutarenotincludedinthecalculationsofarasitisassumedthat ray flux from dark matter annihilation/decay in a wide range a theparticlephysicsisindependentoftheastrophysicsoftheproblem. of astrophysical configurations and making it available to the Runningtime:thisishighlydependentoftheDMprofilesconsidered, community should prove very useful. This is the motivation therequestedprecisionǫandintegrationangleαint: for developing the code presented in this paper. This paper about60mnfora5 5 maptowardstheGalacticcentre,with highlightsclumpy’smainfeaturesbutamorethoroughdescrip- ◦ ◦ • × α =0.01 ,NFWdarkmatterprofilesandǫ=10 2; tioncanbefoundinthedocumentationcomingwiththesource int ◦ − code. The paper is organised as follows. Section 2 presents about2hforthesameset-uptowardstheanti-centre; • the formulation of the generic γ-ray signal calculation that is 0.1to10DMmodelspersecond,dependingonintegrationangle performed by clumpy, with special emphasis on the contribu- • andDMprofile. tionfromDMclumps. Thecodeisdescribedinsect.3,which presentsthe structure,parametersandfunctionsthatare atthe coreofclumpy. Section4dealswithtwoexamplesofthesci- Emailaddress:[email protected](Ce´lineCombet) 1http://root.cern.ch encethatcanbeperformedwithclumpy. Weconcludeanddis- 2http://www.doxygen.org cussthefuturedevelopmentsofthecodeinsect.5. PreprintsubmittedtoComputerPhysicsCommunications January24,2012 2. Calculating the γ-ray flux from dark matter annihila- squaredalongthelineofsight(l,ψ,θ),intheobservationalcone tion/decay ofsolidangle∆Ω(seeAppendix Aandgeometry.hfordefi- nitionsandthegeometryused) The γ-ray flux dΦ /dE from dark matter annihilat- γ γ ing/decaying particles is expressed as the product of a parti- ∆Ω J(ψ,θ,∆Ω)= ρ2(l(ψ,θ)) dldΩ , (3) cle physicstermbyanastrophysicalcontribution. Foragiven Z Z 0 l.o.s experimentwithaspatialresolutionα correspondingtoanin- int whereρ(l(ψ,θ))isthelocalDMdensityatdistancelfromEarth tegrationsolidangle∆Ω = 2π(1 cos α ), atenergyE and − int γ inthedirectionψ,θ (longitudeandlatitudeinGalacticcoordi- pointinginthedirection(ψ,θ),thefluxiswrittenas: nates). dΦ dΦPP Thisintegralisthecoreoftheclumpycodeandcanberather γ γ (Eγ,ψ,θ,∆Ω)= (Eγ) J(ψ,θ,∆Ω). (1) complicatedtoevaluate.Thedifficultycomespartiallyfromthe dE dE × γ γ steepnessofsomeDMprofilesconsideredbutalso,andmainly, 2.1. Annihilationanddecay:generalstatements from the existence of DM clumps, the contribution of which must be added to the smooth Galactic DM background. Both 2.1.1. Annihilation:theparticlephysicsterm componentsaredealtwithhereafter. Itcanbegenericallyexpressedas ddΦEγPP(Eγ)≡ 41π hσ2amnn2vi · ddNEγf Bf, (2) 2.1γ.3-r.ayDseccaanyianlsgodbaerkprmodauttceerdfromthedecayofDM.Formally, γ χ Xf γ thefluxisstillgivenbyEq.(1)buttheastrophysicalandparticle physicstermstakedifferentforms: wherem is themassoftheDM particle, σ istheannihila- χ ann tion cross section, σ v the annihilation rate averaged over h ann i The particle physics term now corresponds to the decay theDMvelocitydistribution,B isthebranchingratiointothe • f spectrum and no parametrisation of that term has been finalstate f anddNγf/dEγ isthephotonyieldperannihilation. included in the code so far. The parametrisation from Awidelyacceptedvaluefor σ v estimatedfrompresentday ann Bertoneetal. [7]couldbeused. h i thermal relic density is σ v = 3 10 26 cm3s 1 [4]3. The ann − − h i · photonspectrumisdependentontheannihilationchannels.We Becausewearenowconsideringadecayreaction,theas- • restrictourselvestothepopularMinimalSupersymmetricStan- trophysics term is the integration of the density (rather dard Model (MSSM). In this framework, the neutralino is the thatdensitysquaredasinEq.(3))alongthelineofsight, lighteststablesupersymmetricparticleandisthefavouredDM namely candidate. Neutralinoannihilationcanproduceγ-raysinthreeways: i) ∆Ω J(ψ,θ,∆Ω)= ρ(l(ψ,θ)) dldΩ . (4) DM annihilates directly into two photons or Z0γ, giving rise Z Z 0 l.o.s to monochromatic lines, ii) or it can annihilate into primary products,thehadronizationanddecayofwhichproduceaγ-ray This quantity can be straightforwardly com- continuum,iii)finally,ifchargedparticlesareproduced,inter- puted in CLUMPY by setting the flag parameter nal bremsstrahlung(IB) can also contribute to the continuum. gDM IS ANNIHIL OR DECAYtofalse(seetable3). The lines are very model-dependentand given the uncertain- ties on the particle physicsmodels, we did not include any in 2.2. Darkmatterdistribution clumpy. Theuserishoweverfreetoimplementtheirownfunc- Todate,thereisnoconsensustowhattheGalacticDMpro- tionifneedbe. Thereexistsintheliteratureseveralparametri- file, ρ(r), should be. There is some dynamical evidence that sationsfortheγ-raycontinuumandwehaveimplementedthat theGalacticDMhalomightbetriaxial[8]andnumericalsim- of Bergstro¨m et al. [4], and Tasitsiomi and Olinto [5]. The ulations, such as the Aquarius[9] or the Via Lactea runs [10] IBcontributionsisalsoverymodeldependent,andweprovide alsofindnon-sphericalhalos. Whenincluded,thebaryonstend onebenchmarkmodelfromBringmannetal. [6]. Wereferthe howevertomakethehalosmorespherical[11]. Inthisfirstver- readerto 3.3.7(spectra.h)formoredetailsonthesespectra. sion of clumpy we simplifythe problem(asoften done)using § sphericallysymmetricDMhalos,butallowingfortri-axialityis 2.1.2. Annihilation:theastrophysicalcontribution being considered for future releases. The total averaged den- Theγ-raysbeingproducedfromtheannihilationofpairsof sityprofilesoftheDMhalosmeasuredinthesesimulationsare darkmatter particles, the γ-rayflux is proportionalto the DM generallyfittedbyaformoftheZhaoprofile[12]: densitysquared.ThesecondtermofEq.(1),termedastrophys- ρ ical contribution,correspondsto the integrationof the density ρ (r)= s , (5) tot (r/r )γ [1+(r/r )α](β γ)/α s s − whereγ,αandβarerespectivelytheinner,transitionandouter 3Notethatavelocity dependent crosssectionshouldbeincluded inafu- slope of the profile, ρ and r a scale density and radius. The turerelease ofthecodeinordertoaccount fortheSommerfeldeffectinthe s s Navarro et al. [13], Moore et al. [14], and Diemand et al. calculation. Thiseffectresultsinmuchlargercrosssectionsforsomeresonant neutralinomassesandthereforeboostthesignalfortheseDMcandidates. [15] profiles fall in this category. Some other profiles having 2 alog-varyinginnerslope[16,17]arealsoincludedinclumpy. butions are assumed independent4. If N is the total number tot The reader is referred to 3.3.4 (profiles.h) for a detailed of clumps in a DM host halo, then the overall distribution of § description. clumpsiswrittenas: SomeoftheseDMprofilesaresteepenoughintheinnerre- d2N d (r) d (M) gionstoleadtoasingularityoftheγ-rayluminosityatthecen- = N PV PM , (8) tot treofthehalo. However,acut-offradiusr naturallyappears, dVdM dV dM cut within which the DM density saturates due to the balance be- where the spatial and mass distribution are probabilities, re- tween the annihilation rate and the gravitational infalling rate spectivelynormalisedas: ofDMparticles. Thesaturationdensityreads[18] dPV(r)dV =1 and Mmax dPM(M)dM =1.(9) ρsat =3.1018(cid:18)100mGχeV(cid:19)× 10−26σcmv3s−1!M⊙ kpc−3. (6) ZV dV ZMmin dM h i AnalysisofN-bodysimulationshaveshownthemassdistribu- Plugged in any profile parametrisation, this saturation density tiontovaryasasimplepowerlaw defines the cut-off radius below which the annihilation rate d (M) is constant. This is a very crude description, but this is not PM M−αM withαM 1.9. (10) dM ∝ ≈ important as this cut-off matters only for the steepest profiles (γ 1.5)whichare disfavouredby currentnumericalsimula- Given the mass and spatial distribution, the total number of tion≥s. clumpsNtotcanbedeterminedintwoways.Ifthemassfraction The hierarchicalscenario of galaxy formationin a Λ-CDM f ofclumpsandthemassofthehostMtot haloareknown, cosmologyis characterisedby a highdegreeof clumpinessof fM theDMdistribution(supportedbyN-bodysimulations)sothat Ntot = tot , M thetotalaverageddensitycomputedfromthesesimulationscor- h 1cli respondsto with M1cl theaveragemassofoneclump. Anotherapproach h i is to know the number of clumps N in a given mass interval 0 ρtot(r)=ρsm(r)+hρsubs(r)i (7) [M1,M2]andtocalculatethetotalnumberofclumpsas whereρsm(r)isthe“true”smoothcomponentand ρsub(r) rep- N = N0 . creosmenptosntehnetsavareeradgiescduesnsesidtybefrloowm.the substructuresh. Thesei two tot RMM12 dPdMM(M) Clumpyusesthefirstapproachwhendealingwithsubstructures in individual host halos (that are not the Milky Way) and the 2.2.1. Theclumps secondonefortheGalactichaloclumps(numericalsimulations Whethersubstructuresinsubhalosarescaled-downversions have found about 100 clumps with masses above 108 M in ofsubstructureinmainhalosremainsanopenquestion[9]. In MilkyWay-likehalos). clumpyallowstoselectanyofthep⊙ro- thisreleaseofclumpy, weonlyconsideronelevelofsubstruc- filesgivenin 3.3.4(profiles.h)fortheclumpspatialdistri- turewithinthehalounderscrutiny(Galactichaloorindividual § bution, independentlyof the choice made for the total density halo such as a dwarf spheroidalgalaxy),and the propertiesof profile. these sub-halos can be independentlychosen from that of the parenthalo(seebelow). 2.2.2. Thesmoothcomponent From the spatial distribution d (r)/dV one can define the Clump individualproperties. Inthe code,allclumpshavethe PV averageclumpdensity ρ usedinEq.(7)as: same inner DM distribution, which can be any of the profiles h cli discussedin 3.3.4.Themassoftheclumpgenerallysufficesto d (r) § ρ (r) = fM PV . (11) determineallitsproperties,i.e.,itssizeRcl ,andonceaninner h subs i tot dV vir density profile is chosen, its scale radius r and scale density s Foragiventotaldensityprofileandclumpspatialdistribution, ρ . This is done using the so-called concentration parameter s Eq.(7)thengivesthesmoothdensityprofileas c . Parametrisations of the mass-concentration relation have vir been established from numericalsimulationsfor isolated field ρ (r)=ρ (r) ρ (r) . (12) sm tot subs halos[19,20]buttheconcentrationgenerallypresentsasignifi- −h i cantscatter[21]. Theserelationshavebeenimplementedinthe Using this approach, rather than providing independently a code, regardlessof that scatter or of the environmentand for- smoothandaclumpprofile,ensures(byconstruction)thatthe mationhistoryofthehalosthatcanalsoaffectthevalueofthe totaldensityprofilesfollowstheonemeasuredinN-bodysim- concentration. We refer the reader to 3.3.6 (clumps.h) and ulations.ThishasbeendiscussedinPierietal. [23]. § theDoxygendocumentationformoredetails. 4KeepinmindthatsuchanassumptionfailsneartheGalacticcentre,where Spatialandmassdistributionoftheclumps. FollowingLavalle smallclumpsareexpectedtobedisruptedbytidalforces.Inthisdirection,the et al. [22], sphericalsymmetryis assumedand the two distri- signalfromthesmoothcontributiondominatesanyway. 3 2.3. CalculatingtheannihilationJ-factor where (M)istheintrinsicluminosityofaclumpgivenby L Theastrophysicalcontributiontotheannihilationfluxisfor- (M) (ρ )2dV . (19) cl mallygivenbyEq.(3). TheJ-factorismoreexplicitlywritten L ≡Z Vcl as: Theaverageofthedistancetoanypowerniswrittenas J = ∆Ω lmax 1 ρ + ρi 2l2dldΩ. (13) ln = ∆Ω lmaxlndPVl2dldΩ. (20) Z0 Zlmin l2  sm Xi cl CombhiniingZth0ese,ZiltmicnanbedVshownthattheaveragecontribution whereρ isgivenbyEq.(12)andthesecondtermcorresponds sm fromtheclumpsisexpressedas to the sum of the inner densities squared of all clumps i con- tainedwiththevolumeelement. Thelattershouldnotbecon- J = N ∆Ω lmax dPVdldΩ Mmax (M)dPMdM ,(21) subs tot fusedwith ρ (r) . Threetermsarisefromthisequation: h i Z Z dV Z L dM h subs i 0 lmin Mmin Inthesamespirit,inthecontinuumlimitthecross-productbe- ∆Ω lmax J ρ2 dldΩ, (14) comes sm ≡Z Z sm 0 lmin ∆Ω lmax J =2 ρ ρ dldΩ. (22) cross prod sm subs ∆Ω lmax 2 h − i Z0 Zlmin h i Jsubs ≡Z Z  ρicl dldΩ, (15) The importance of the cross product can easily be identified 0 lmin Xi  wtohtaelnacvoenrasgideerdienngstihtyatptrhoeficlelu.mInptshpaatticaalsdei,strρibuti(orn) fo=llofwρst(hre) andthecrossproduct h subs i tot andρ (r) = (1 f)ρ (r),leadingdirectlyto J /J = sm tot cross prod tot ∆Ω lmax 2f(1 f). Fo−r a typical clump mass fraction f− = 10%, Jcross−prod ≡2Z0 Zlmin ρsmXi ρicldldΩ. (16) tJhis c−rosρs2pdrloddΩuc.tTahmisosuhnotwsstoth1a8t%wheonf athneavreerfaegreendcceluJm-fapcdtoer- tot ≡ tot scriptiRonis usedto integratethesignalalongthe lineof sight, 2.4. Statisticalconsiderations this cross product is not negligible and should be taken into Taking α = 1.9 as a canonical value for the slope of account. In Appendix B, it is however shown that on an in- M the mass distribution, with 100 clumps above 108 M as de- dividual basis, the cross product of one clump in an underly- termined from N-body simulations, and a minimal the⊙oretical ing smoothGalactic componentcan besafely discarded. This mass of 10 6 M for the smallest ones, the total number of willapplytotheclumpsthatneedtobestatisticaldrawninour − clumpsintheMi⊙lkyWayhalo(i.e.,for∆Ω = 4π)isestimated model(seebelow). tobe 1014,whichiscomputationallyprohibitive. Evenifwe ∼ 2.4.2. Varianceoftheclumpcontributionσ2 are interestedin a5◦ 5◦ skymaponly(i.e., ∆Ω 6 10−3), clumps thenumberofclumps×is 5 1010. Thisistoolar∼gea×number Alongagivenlineofsightandforagivenintegrationangle ∼ × α , the average flux between l and l given by Eq. (21) torealiseastatisticaldrawingaccordingtothed /dV distri- int min max V P canalsobeexpressedas butionandthennumericallyintegratetheJcomponentforeach oftheseclumps. However,theirhighnumbermeansthatasta- J = N J , (23) subs cl 1cl tisticalaveragecanbesafelyperformed,aslongasthevariance h i h ih i is small. The clumpy code computes the γ-ray flux from DM where hNcli is the average number of clumps in this volume (given the distance interval and mass range considered), and clumpsbyeitherrandomlydrawingclumpsgivenuser-defined J theaveragefluxofoneclump(seealso[24]). Usingthe spatialandmassdistributionsandinnerDMprofiles,orbycal- h 1cli point-likeapproximation,aclumpofmassMatadistancelhas culating the mean contribution whenever the variance of this aflux contributionissmallerthansomeuser-definedvalue. (M) J = L , (24) 1cl l2 2.4.1. Averagedquantitiesfromthecontinuumlimit sothat,recallingthatthemassoftheclumpanditslocationare Giventheclumpmassandspatialdistributionsgivenabove, independentvariables, it is possible to define average values for the primary random variables of the problem (a clump mass and position) and the 1 J = . (25) one deriving from them. The mean mass and the mean lumi- h 1cli hLi*l2+ nosityarerespectivelydefinedby Thevarianceonthefluxofasingleclumpisthendefinedas M = Mmax MdPMdM , (17) σ2 = J2 J 2 = 2 1 J 2 . (26) h i Z dM 1cl h 1cli−h 1cli hL i*l4+−h 1cli Mmin Thevarianceofapopulationof N clumpsfollowsas cl = Mmax (M)dPMdM . (18) σ2 = N σ2 . h i (27) hLi ZMmin L dM cl h cli 1cl 4 103 102 aTnatbil-ece1n:trEex(acmriptelreioonflRcrEitJacnludm<ptsh5e%as)sfoocriaate2◦n×2u◦msbkeyrmoafpclwumithpsantoindtreagwrattoiownaradngthlee 101 (instrument)ofαint=0.1◦.Actually,agivenlcritgives<ncl>numberofclumps: 100 nMclatossdrdaewcaisdeobta#inceldumfropmsaPolisso(nkdpicst)ributio<n#ocflmumeapnsv>alue<#nctlo>.draw 10-1 Ψ = 180o log Mcl (Galaxy) kcprict (5 5 ) (5 5 ) %) 10 M ◦× ◦ ◦× ◦ RE (Jcl 1100--23 Ψ = 90o −−456:::(cid:16)−−345⊙(cid:17) 543...645××111000111234 111...000××111000−−333 211...485××111000−−543 000 10-4 − − × × − × − 3: 2 7.0 1011 7.3 10 3 1.2 10 3 0 − − 1100--65 ααiinntt == 00..01o1 o Ψ = 0o −−−21::−−01 18..18×××11001100 45..08×××1100−−12 27..86××10−2 03 10-7 1 - 10 M 0:1 1.4×109 2.0 3.7×10+1 45 sun 1:2 1.7 108 6.7 1.0 10+2 117 10-08.01 1 100 2:3 2.2×107 1.6 101 8.7×10+1 99 l (kpc) × × × min 3:4 2.8 106 3.3 101 3.6 10+1 44 × × × 4:5 3.5 105 5.8 101 1.0 10+1 7 Figure 1: Relative errorbetween theaverage andtrue J value from clumps 5:6 4.4×104 9.5×101 2.2× 1 whenthesignalisintegratedbetweenlminandlmax,accordingtoEq.(28),for 6:7 5.6×103 1.5×102 4.1 10 1 0 αMinctl =∈0[.10,11◦0](dMas⊙h.edTlhinee)s.igTnahlreiesdinirteecgtriaotnesdionvtehreαGinatla=cti0c.1p◦lan(seoalirdelcionne)sido-r 7:8 7.0××102 2.5××102 6.6××10−−2 0 ered: towards theGalactic centre (blue), towards theanti-centre (black) and 8:9 8.8 101 2.9 102 9.5 10−3 0 × × × towardstheGalacticEast(red). Forthelatterthegreendottedlineillustrates 9:10 1.1 101 2.9 102 1.3 10 3 0 − × × × howtogetthecriticaldistance(here2kpc)belowwhichclumps(inthismass range)mustbedrawniftheuserrequiresanaccuracyof5%. Alldistributions (smooth, clumpdistribution andinnerDMclumpprofile)aretakenasNFW Threshold mass mthresh for the sub-clumps in an “far-away” profiles. halo. The second case where the criterion is used is for the sub-clumpsinahalolocatedawayfromus,e.g.,adSphgalaxy. 2.4.3. Criterionforusingtheaverageddescription Inthatcase, thedistanced totheclumpsisknownastheyare For a given integration domain, there is a threshold mass alldistributedwithintheirparenthalo’sradiusandwecande- above which the clumps are not numerous enough to be determine in a similar fashion the threshold mass above which scribed by an averaged description. Conversely, for a given clumpsmustbedrawn. mass decade there is a critical distance below which the aver- ageddescriptionfailsandwhereclumpsshouldbestatistically 3. Descriptionofthecode drawn. The averageddescriptioncan be safely used, as long as the Clumpy is written in C/C++ and relies of the CERN ROOT relative error RE made on the clump flux with respect to the librarythatispubliclyavailablefromhttp://root.cern.ch. totalflux(i.e. includingthesmoothcontribution) Themainfeaturesofthecodearedescribedhereafteralongwith brief descriptionsof its most importantfunctions. The code’s √N σ REJclumps = N J cl+J1cl (28) structure is standard, with separate directories for the various clh 1cli smooth pieces of code: declarations are in include/*.h, sources in issmallerthansomeuser-definedprescription. src/*.cc,compiled libraries, objects and executablesare re- spectivelyinthelib/,obj/,andbin/directories. Critical distance l for the clumpsin the Galaxy. TheRE is crit plottedinFig.1asafunctionofthelowerlimitoftheintegra- 3.1. Codeexecutables tionl . Ascanbeseen,foranyrequiredprecision(e.g.,5%in min The code has one executable that can take three main op- Fig.1,taggedbythedottedgreenline),onecaninferthecorre- tionsdependingwhethertheuseri)isinterestedinperforming spondingdistance l below which the relative errorbecomes crit agenericanalysisontheGalactichalo,ii)isfocusingonalist largerthantherequest. ofspecifichalos(e.g. dSphgalaxies),oriii)wishestoperform In practise, it is convenient to perform the calculation in astatisticalanalysisonasingleobject.Thesemainthreeusages termsofmassdecades:foreachdecade,wefind(bydichotomy) the critical distance below which clumps must be drawn5. A ofthecodeallaccommodateseveralcases,thedetailsofwhich lower limit of 10 3 kpc is set on l to avoid divergencesfor aregivenbelow. − crit a clumpsitting at the observer’slocation. An exampleforthe 3.1.1. Galacticmode:./bin/clumpy -g[option] finalnumberofclumpstodrawbelowl isprovidedinTab.1. crit This calculates the J-factor for the Galactic halo, including TheflowchartgiveninFig.2summarisesthevariousstepsim- plementedinclumpytoobtaintheskymap(see,e.g.,Fig.3). boththesmoothandtheclumpcomponent.Fortheclumpcom- ponent,themeanfluxorastochasticrealisationcanbechosen. Note that a list of specific halos can be added to the calcula- 5ThismassdecadecorrespondstoatypicalclumpradiusRcl andwealso tion. Several options allow to select the quantity to be calcu- checkthatthepoint-likehypothesisusedtoderivedthecontinuumlimitisap- propriate(lcrit≫Rcl). lated/plotted from a text interface. A mandatory input is the 5 READ IN USER−DEFINED PARAMETERS CLUMPS DEFINE FIELD OF VIEW SMOOTH COMPUTE L_CRIT(M) GIVEN CLUMP DISTRIB For a given mass decade, need to determine the distance l_crit below which the analytical clump estimation fails and clumps need to be drawn from dP/dV COMPUTE NUMBER OF CLUMPS TO COMPUTE ANALYTICAL CONTRIBUTION COMPUTE SMOOTH CONTRIBUTION DRAW IN THE FIELD OF VIEW OF THE CLUMP "CONTINUUM" e FROM 0 to L_CRIT(M) d a c e d DRAW CLUMPS s s a m ADD ALL CONTRIBUTIONS PIXEL BY PIXEL h IDENTIFY CLUMPS BELONGING c TO THE F.O.V a e r o F COMPUTE CLUMP FLUX AND STORE IN MAP PIXELS STORE ALL COMPONENTS IN FILES AND/OR DISPLAY WITH ROOT Figure2:Flowchartofthedifferentstepsleadingtothecreationofaγ-rayskymap. Table2:FormatoftheASCIIfiletodescribehalosprocessedby./bin/clumpy dsph -h[option].’Type’shouldtakeoneofthevaluesofgENUM TYPEHALOES, namelyDSPH,GALAXYorCLUSTER.’long’,’lat’and’d’arethegalacticlongitude,latitude,anddistancetothehalo,zitsredshiftandRviritsvirialradius.Theprofile parametersareρs,rs,’prof.’ (seeprofile.h)andtheshapeparameters. Inparticular,if’prof’=kZHAO,see 3.2),theprofilecorrespondstothe(α,β,γ)profile § (inwhichcasethenext-to-lastthreecolumnsareused).IfankEINASTOisusedinstead,onlythefirstshapeparameterisconsidered. Name Type long. lat. d z R ρ r prof. α β γ vir s s - - deg deg kpc - kpc M kpc 3 kpc - - - - − ⊙ parameter file clumpy params.txtdescribed in 3.2 and ta- ’option = -g4’ is the same as ’-g3’ but including a § • ble3. Theoptionsarethefollowing: listofspecificobjects,suchasdSphgalaxies. ’option = -g1’ plots the total, smooth and averaged ’option = -g5’plotsa 2D skymapof the J-factor, us- • • clump density in 1D, given the user-definedtotal density ingtheaveragedclumpdescription.where andthenumberofclumpsinagivenmassrange. ’option = -g6’issimilar ’-g5’,butaddsspecificha- • ’option = -g2’ plots J + J + J 6, re- lostotheplots(e.g.dSph),toidentifytheircontrasttothe sm subs cross prod • spectivelydefinedbyEq.(14),Ehq.(2i1)ahndEq−.(22i)asa background. functionoftheintegrationangleα . int ’option = -g7’ and -g8 are respectively the same as • ’-g5’ and ’-g6’ but with the clumps randomly drawn ’option = -g3’isthesameas’-g2’butasafunction • fromtheiruser-definedmassandspatialdistributions. of the angulardistanceto the Galacticcentre, fora given integrationangle. 3.1.2. Halolistmode: ./bin/clumpy -h[option] Thisusageofclumpyperformsoperationsthatarebasically 6DoesnotexistfordecayingDM. the identical to those performed in the ’Galactic’ mode (see 6 3.1.1), but for a list of halos that are not the Galactic halo. enum gENUM GAMMASPECT kSUSY BUB98, § • { This mode was initially implemented for the study of dwarf kSUSY TO02, kSUSY BBE08 fortheγ-rayspectrum; } spheroidalgalaxies[25], butcan be extendedto anyDM halo enum gENUM PROFILE kZHAO, kEINASTO, (e.g., galaxyclusters). If onlythe smoothtotaldensity profile • { kEINASTO N for the DM spatial distributions (smooth is considered, the halo profiles are read from a simple ASCII } and clumpy components). The corresponding shape list of profiles along with their parameters, as exemplified in Tab.2. Ifsub-clumpsareconsidered,theclumpyparameterfile parameters(e.g. (α,β,γ)forkZHAO,αforkEINASTOand n forkEINASTO N) shouldbeprovidedbythe user in the clumpy params.txtmustalsobeprovided(seetable3). inputparameterfile(seebelow). Again, furtheroptionsareaccessed forthismodebymeans ofasimpletext-userinterface: Thisenumerationsallowtheusertousetheirexplicitdenomina- tionkXXX,ratherthanthecorrespondingintegernumber.These ’-h1’to’-h3’areidenticalto’-g1’to’-g3’buttrans- • parametersaregatheredintheASCIIfileclumpy params.txt posedtoanyhalothatisnottheGalacticDMhalo; which is formatted as “Parameter name”, “Units”, “Value”. ’-h4’ to ’-h5’ create 2D skymaps of the halo under The function load parameters(file name) defined in • scrutiny using respectively an averaged description or a params.h reads in the input file and assigns the appropriate statisticalrealisationforthesub-clumpswithinthishalo; valuesto the parameters. The user canchangeany“Value”in thisfile7. 3.1.3. ’Statistical’mode:./bin/clumpy -s[option] Confidence levels on the astrophysical factor J, for e.g. 3.3. clumpylibraries dSphs,canalsobeestimatedfromanASCIIstatisticalfile.The Thefunctionsofclumpyarelocatedineleven“library”files, formatforthelatterisquitesimilartothehalolistfile(seedoc- accordingly to their field of action. We briefly describe here- umentation) as it requiresthe profile parameters, but it differs after the most important functions of these libraries. We are in the sense that for a single halo, thousands of profiles ’se- omitting in this document the more minor functions, the de- lected’byastatisticalanalysismustbeprovided(theirrespec- scriptionofwhichcanbefoundinthedocumentationprovided tive χ2 value isnecessary, althoughnotused forall CL calcu- withthecode. Also,wedonotdetailherewhatargumentsare lations). Such a statistical analysis can be the Markov Chain passedinthesefunctionsasitcanbeeasily retrievedfromthe MonteCarlo(MCMC),asdiscussedinCharbonnieretal. [25]; documentation. Walkeretal. [26]. Themostimportantoptionsinthestatistical modeare: 3.3.1. integr.h This library contains standard integration algorithms. The ’-s2’ plots the probability density functions (PDF) and • Simpsonintegratorwithdoublingstepisthemostwidelyused correlationsoftheparametersgiveninthestatisticalfile; throughoutthecodetoensurethattheuser-definedprecisionis ’-s4’locatesthe x%confidenceintervalsonthePDFof reached. Depending on the function to integrate, either linear • eachparameters; orlogarithmicsteppingcanbeused. ’-s5’computes,foreachradiusfromthedSphcentre,the • 3.3.2. geometry.h medianandconfidencelevelsoftheDMdensityρ(r); In geometry.hare gatheredfunctionsperforminggeomet- ’-s6’isthesameas’-s5’butfor J(α ),i.e. asafunc- ricaltasks(changeofcoordinates,etc.). Theyarenotessential int • tionoftheintegrationangle; forthepresentationofthecodeandgenerallyself-explanatory. We refer the reader to the extensive Doxygen documentation ’-s7’issameas’-s5’butfor J(θ),i.e. asafunctionof andtoAppendix Aformoredetails. • theangletothecentreofthehalo. 3.3.3. inlines.handmisc.h 3.2. User-definedparameters These contain some definitions (unit conversion and inline Formostoftheoptionsusedabove,auser-defineparameter functions)andmiscellaneousformattingfunctionsrespectively. fileisread. Table3givesthelistofthemainparameters. Each parameteris a single word formattedas g (forglobalparame- 3.3.4. profiles.h ter)+physicskeyword+parametername.Forinstance,COSMO Thislibrarycontainsallthe densityprofilesavailable in the refersto cosmology-relatedparameters, TYPEdefinesthe type codeandoperationsonthose.Amongthefunctionsinthe“pro- ofhalounderscrutiny(dSPhgalaxy,galaxyorgalaxyclusters) files” section of the code, some act in a very straightforward andGALreferstotheMilkyWayhaloparameters.Notealsothat way (e.g. squaring the density) and are not included in the the parameters ending in FLAG CVIRMVIR, FLAG SPECTRUM, present document. The reader is, once again, referred to the andFLAG PROFILEshouldtakevaluesfromthefollowingenu- Doxygendocumentationforafulldescription. merations: enum gENUM CVIRMVIR kB01, kENS01, kJS00 for 7Some of these parameters are optional and used only for specific • { } theclumpconcentration; modes[options] 7 Table3:clumpyuser-definedinputparameterfile(clumpy params.txt). Name Definition gCOSMO RHO0 C Criticaldensityoftheuniverse[M kpc 3] − gCOSMO OMEGA0 M Present-daypressure-lessmatterde⊙nsity gCOSMO OMEGA0 LAMBDA Present-daydarkenergydensity gDM GAMMARAY FLAG SPECTRUM γ-rayspectrumflag[gENUM GAMMASPECT] gDM MMIN SUBS Min.massofaDM(sub-)clump[M ] gDM MMAXFACTOR SUBS Definesmax.massofclumpinhost⊙haloasM =Factor M max × host gDM RHOSAT Saturationdensity[M kpc 3] − ⊙ gTYPE CLUMPS FLAG CVIRMVIR ClumpconcentrationflagintheparentTYPEhalo[gENUM CVIRMVIR] † gTYPE CLUMPS FLAG PROFILE Clumpinnerprofileflag[gENUM PROFILE] gTYPE CLUMPS SHAPE PARAMS[0-2] Shapeparametersforsub-clumpsinnerprofile gTYPE DPDM SLOPE Slopesubclumpmassfunction gTYPE DPDV FLAG PROFILE Spatialdistributionsubclumps[gENUM PROFILE] gTYPE DPDV RSCALE Scaleradiusforthesubclumpdistributionintheparenthalo[kpc] gTYPE DPDV SHAPE PARAMS[0-2] Shapeparametersfortheclumpspatialdistribution gTYPE SUBS MASSFRACTION Massfractionoftheparenthaloinsub-clumps gGAL CLUMPS FLAG CVIRMVIR Concentrationflag[gENUM CVIRMVIR] gGAL CLUMPS FLAG PROFILE Clumpinnerprofileflag[gENUM PROFILE] gGAL CLUMPS SHAPE PARAMS[0-2] ShapeparameterfortheGalacticclumpsinnerprofile gGAL DPDM SLOPE Slopetheclumpmassfunction gGAL DPDV FLAG PROFILE Clumpnumberdistributionflag gGAL DPDV RS Scaleradiusforclumps[kpc] gGAL DPDV SHAPE PARAMS[0-2] ThreeshapeparametersfortheGalacticclumpspatialdistribution gGAL SUBS M1 ReferencemassforgGAL SUBS N INM1M2[M ] gGAL SUBS M2 ReferencemassforgGAL SUBS N INM1M2[M⊙] gGAL SUBS N INM1M2 #ofclumpsin[M ,M ] ⊙ 1 2 gGAL RHOSOL LocalDMdensity[GeVcm 3] − gGAL RSOL DistanceSun–Galacticcentre[kpc] gGAL RVIR VirialradiusoftheGalaxy[kpc] gGAL TOT FLAG PROFILE TotalDMprofileoftheMilkyWay[gENUM PROFILE] gGAL TOT RSCALE ScaleradiusforDMhalooftheMilkyWay[kpc] gGAL TOT SHAPE PARAMS[0-2] ShapeparameterfortheGalactictotaldensityprofile gSIMU ALPHAINT DEG Integrationangle[deg] gSIMU EPS Defaultprecisionfornumericalintegrations gSIMU IS ANNIHIL OR DECAY ForannihilatingordecayingDM gSIMU SEED Seedofrandomnumbergenerator †TYPEcorrespondstothetypeofhalounderscrutiny(nottheGalactichalo):DSPH,GALAXYorCLUSTERarethevalues allowedbygENUM TYPEHALOESinthisversion. 8 rhosimplyreturnsthevalueofthedensitybycallingany Table4:StandardDMdensityprofilesparameters,followingEq.(29). • ofthefunctionsabove,dependingonthevalueofthepro- α β γ Reference fileflaggiveninclumpy params.txt. ISO 2 2 0 - NFW97 1 3 1 Navarroetal.[13] Fortwodensityprofilesρ andρ ,rho mixreturnseither M98 1.5 3 1.5 Mooreetal.[14] • 1 2 thedifferenceρ ρ ortheproductρ ρ . Thisisneeded DMS04 1 3 1.2 Diemandetal.[15] 1− 2 1 2 whencomputingEq.(12)forthesmoothdensityprofileor Eq.(22)forthecrossproduct. rho ZHAOreturnstheDMdensityofauser-definedprofile • havingthegenericfunctionalformgiveninZhao[12] 3.3.5. integr los.h This is the core of the clumpy code as it contains the rou- ρ ρ(r)ZHAO = s . (29) tinesperformingtheintegrationalongthelineofsightdirection (r/rs)γ [1+(r/rs)α](β−γ)/α (θ,ψ)andoverthechosenangularresolutionαint.Givenafunc- tiong(l,α,β;θ,ψ)where(l,α,β)aresphericalcoordinatesinthe Theisothermal,Navarroetal. [13],Mooreetal. [14],and (θ,ψ) direction(see Fig. 6), the integrationalongthe observa- Diemandetal. [15]profilesareallfromtheZhaofamily. tioncone,froml tol isgenericallywrittenas: 1 2 TheyaredefinedusingEq.(29)withthesetofparameters oftable4. Theuserisrequiredtospecifytheseparameters 2π αint l2 I(θ,ψ)= dβ sinαdα g(l,α,β;θ,ψ)l2dl.(32) intheASCIIinputfile. Z Z Z 0 0 l1 rho EINASTOandrho EINASTO N:ithasbeenarguedthat integrand l: returns the value of the function • • aprofilewithalogarithmicinnerslopegivesabetterfitto f(l,α,β;θ,ψ) g(l,α,β;θ,ψ)l2. A switch is includedto ≡ the simulated data. This so-called Einasto is generically integratedifferentfunctionsofthedensityprofileρ. Note writtenas thathereρcanrefertoanyprofileordistribution,beitthe density or the spatial distribution of clumps dP /dV de- V ρ(r) exp( Arα) . pendingonthenormalisationused. Thefunctionsthatcan ∝ − beintegratedalongthelineofsightare: In Navarroet al. [16] and Springelet al. [9] the Einasto – f(l,α,β;θ,ψ)=ρwhencalculatingtheclumpsaver- profileisdefinedas agedcontribution J ofEq.(21). subs h i – f(l,α,β;θ,ψ) = ρ2 in orderto evaluateJ such as in α 2 r ρEINASTO(r)=ρ exp 1 , (30) Eq.(14). −2 −α" r−2! − #! – f(l,α,β;θ,ψ) = l2ρ is used to calculate the mass fractionornumberfractionofclumpsintheintegra- where r is the radius at which the profile’s slope 2 tionvolume; dlnρ/dl−nr = 2, and ρ the corresponding density. 2 Both Navarro et−al. [16]−and Springel et al. [9] find – f(l,α,β;θ,ψ) = l3ρ is needed when estimating the α 0.17tobeagoodfittothesmoothcomponentofsim- meandistanceoftheclumpsasinEq.(20). ∼ ulatedhalos. Springeletal. [9]furthernotesthatα 0.68 – f(l,α,β;θ,ψ) = l4ρ is needed when estimating the ∼ allowsagooddescriptionofthespatialdistributionofsub- varianceofdistanceoftheclumps. structuresdP/dV. – f(l,α,β;θ,ψ) = l 2ρ must be evaluated to evaluate − Merrittetal. [17]useanEinastor1/n profile,slightlydif- thevarianceonJasinEq.(26). ferentfromtheaboveandgivenby Note that if decaying rather than annihilating DM is be- r 1/n ingconsidered,i.e. gDM IS ANNIHIL OR DECAYisfalse ρEINASTO N(r)=ρ exp d 1 . (31) insteadoftrue,anyρ2isreplacedbyρintheabovefunc- e − n· re! −  tions8. A variation of this function, called integrand lrel, is where n = 6 and d = 3n 1/3+0.0079/n. The scale n − usedwheneverthesteepnessoftheprofilerequiresit(see radius r in this relation correspond to the radius within e Doxygendocumentationformoredetails). whichhalfofthehalomassisenclosed. fn beta alpha:computesthethirdintegrandinEq.(32), The normalisation density of all these profiles (ρ , ρ • • or ρ ) can be done by i) either requesting a givensvalu−e2 namely the integration over l. The functions to inte- e grateroughlybehaveaspowerlaws,withavaryingslope, for the density at a given point or ii) by requesting a given mass within a given radius. The functions set par0 given rhoref and set par0 given mref 8Thismakessomeofthesefunctionsredundantfordecay,buttheyarestill achievethenormalisationusingbothapproaches. requiredforannihilation. 9 and log-step Simpson integration schemes are used (see frac nsubs in m1m2returnsthefractionofclumpswith • integr.h). Some profiles can be very steep and several a mass in the range [m ,m ], given the total number of 1 2 testsandtricksareneededtoensurethenumericalstabil- clumps. ity/optimisationoftheintegration. Thereaderisagainre- ferred to the documentationfor an explicitdescriptionof jcrossprod continuum returns the value of the cross- • theintegrationstrategy. productintegrated along the l.o.s as defined by Eq. (22). Thisisasimplecalltolos integral mix. fn beta: computes the integration over the angle α in • Eq. (32). A simple linear Simpson integration scheme is jcrossprod interpolationinterpolatesthevalueof J used. • in the pixelsof a skymap, providedthat J haspreviously beenestimatedinsomelocationsofthemap. Thismakes los integral:itisthefinalpartofthetripleintegration, • use ofROOT interpolationfunctionbased on Delaunay’s i.e. the integration over the angle β. This is again per- triangles and allows to speed-up the calculation signifi- formed quite easily with a linear Simpson scheme. This cantlywhentheskymaphasa(too)largenumberofpixels. functioniscalledwheneverasingledensityprofileistobe dealtwith. j smooth returns the J-factor when a single density • los integral mix: as seen previously, it may be nec- profile is used (e.g. ρtot). This is a simple call to • essary to integrate a difference or a product of two los integral. densities along the line of sight. A flag parameter j smooth mix is the same as above but deals with a switch rho is introduced to tell rho mix to return ei- • combination of densities, e.g. when computing Eq. (14) ther ρ ρ (switch rho = 0) or ρ ρ (switch rho = 1 2 1 2 − where ρ = ρ ρ . This is a simple call to 1) (see profiles.h). Thisfunctionthenbehavesidenti- sm tot cl − h i los integral mix. cally to los integral and uses the same flags with re- specttothequantitytointegrate. j sub continuum computes and returns the averaged J • generatedby the sub-clumpsin a halo, in a given line of 3.3.6. clumps.h sight. Allthefunctionsrelatedtotheclumpsmassandspatialdistri- butionsaregatheredinthisfile. Here,weonlygiveanoverview lum singlehaloreturnsthe totalintrinsicluminosityof ofthemostimportantones. • aDMhalo,givenitsinnerdensityprofileasinEq.(19). add halo in map: when using the code to generate a • skymap where clumps are randomly drawn from their mass singlehalo returns the total mass of a DM halo, • mass and spatial distribution, this addsthe J contribution givenitsinnerdensityprofile. of the drawnclumpto theappropriatepixels. Theclump isonlyaddedtoapixelifJcl(pixel)>ǫ×Jsm(pix0)where 3.3.7. spectra.h pix0 is the pixel containing the centre of the clump and Thislibrarycontainsthedifferentparametrisationsoftheav- ǫ is the requested integration precision gSIMU EPS. This erageγ-rayspectra, avoids the clump from being integrated “too far” which would only be time consuming without changing the re- dN dNf γ B (33) sult. dE ≡ dE f γ Xf γ dpdv XXXandset normprob dpdv: theformerfunction • intheparticlephysicstermofEq.(2). Asmentionedin 2.1.1, returnsthespatialdistributionoftheclumpsasdescribed § we have included the Bergstro¨m et al.[4] and Tasitsiomi & in 2.2.1 and the second calculates the appropriate nor- § Olinto [5] parametrisations for the γ-ray continuum and the malisation to make it a probability according to Eq. (9). Bringmannetal. [6]expressionoftheinternalbremsstrahlung. The functiondpdv XXX is written in several versions de- Defining x E /m , the three parametrisations are given by pending of the reference frame and coordinate system in γ χ ≡ theircorrespondingfunctionsbelow. whichtheclumpsaredrawn(cartesianXXX=xyz,spherical centredonthe galacticcentreXXX=rthetapsi,spherical dNdE BERGSTROMreturns centredonEarthXXX=lthetapsi,orcentredonEarth,in • thecone∆Ω,inthedirection(ψ,θ)XXX=lalphabeta). dN 1 exp( 7.8x) γ =0.73 − . dpdm and set normprob dpdm: this is the equivalentof dE · m · x1.5 γ χ • the above but for the mass distribution of the clumps as givenbyEq.(10). dNdE TASITSIOMIcomputestheγ-rayspectrumas • frac nsubs in foi returns the number fraction of • dN 1 10 5 5 5 clumpsinagivenfieldofintegrationw.r.t. thetotalnum- γ = x0.5 x 0.5+ x 1.5 . − − berofclumpsinthehosthalo dEγ mχ · 3 − 4 − 2 12 ! 10