Table Of Content

DRAFTVERSIONJANUARY20,2011 PreprinttypesetusingLATEXstyleemulateapjv.11/10/09 GALAXIA:ACODETOGENERATEASYNTHETICSURVEYOFTHEMILKYWAY SANJIBSHARMA,JOSSBLAND-HAWTHORN1 SydneyInstituteforAstronomy,SchoolofPhysics,UniversityofSydney,NSW2006,Australia KATHRYNVJOHNSTON DepartmentofAstronomy,ColumbiaUniversity,NewYork,NY-10027 1 1 JAMESBINNEY 0 RudolfPeierlsCentreforTheoreticalPhysics,1KebleRd,Oxford,OX13NP,UK 2 DraftversionJanuary20,2011 n a ABSTRACT J We present here a fast code for creating a synthetic survey of the Milky Way. Given one or more color- 8 magnitudebounds,asurveysizeandgeometry,thecodereturnsacatalogofstarsinaccordancewithagiven 1 modelof the MilkyWay. The modelcan be specified by a set of density distributionsor as an N-bodyreal- ization. We providefastand efficientalgorithmsforsamplingboth typesof models. As comparedto earlier ] samplingschemeswhichgeneratestarsatspecifiedlocationsalongalineofsight,ourschemecangeneratea A continuousandsmoothdistributionofstarsoveranygivenvolume. Thecodeisquitegeneralandflexibleand G canacceptinputintheformofa starformationrate,agemetallicityrelation,agevelocitydispersionrelation . andanalyticdensitydistributionfunctions. Theoreticalisochronesarethenusedtogenerateacatalogofstars h and support is available for a wide range of photometric bands. As a concrete example we implement the p Besançon Milky Way model for the disc. For the stellar halo we employ the simulated stellar halo N-body - o modelsofBullock&Johnston(2005). InordertosampleN-bodymodels,wepresentaschemethatdisperses r thestarsspawnedbyanN-bodyparticle,insuchawaythatthephasespacedensityofthespawnedstarsiscon- t s sistentwiththatoftheN-bodyparticles.Thecodeisideallysuitedtogeneratingsyntheticdatasetsthatmimic a near future wide area surveyssuch as GAIA, LSST and HERMES. As an applicationwe study the prospect [ ofidentifyingstructuresinthestellarhalowithasimulatedGAIAsurvey. Weplantomakethecodepublicly availableathttp://galaxia.sourceforge.net. 1 v Subjectheadings:Galaxy:stellarcontent–structure–methods:dataanalysis–numerical 1 6 5 1. INTRODUCTION to create a Galaxy model that is constrained by observa- 3 Generatingasyntheticcatalogofstarsinaccordancewitha tions. The earliest such attempt was by Bahcall&Soneira . givenmodelofgalaxyformationhasanumberofuses. First, (1980a,b, 1984) where they assumed an exponential disc 1 with magnitude dependent scale heights. An evolutionary 0 it helps to interpret the observational data. Secondly, it can model using population synthesis techniques was presented 1 beusedtotestthetheoriesuponwhichthemodelsarebased. byRobin&Creze(1986). Givenastarformationrate(SFR) 1 Moreoversyntheticcatalogscanbeusedto testthe capabili- andaninitialmassfunction(IMF),onecalculatestheresult- : tiesofdifferentinstruments,checkforsystematicsanddevice v ing stellar populations using theoretical evolutionary tracks. strategiesto reducemeasurementerrors. Thisis well under- i Localobservationswere then used to constrain the SFR and X stood by the architects of galaxy redshift surveys who rely IMF. Bienaymeetal. (1987) later introduceddynamical self heavily on ΛCDM simulations to remove artifacts imposed r consistency to constrain the disc scale height. The present a bytheobservingstrategy(Collessetal.2001). stateoftheartisdescribedinRobinetal.(2003)andisknown Giventhewidespreaduseofsyntheticcatalogs,aneedfor as the Besançon model. Here the disc is constructed from faster andaccuratemethodsto generatestellar syntheticcat- a set of isothermal populations that are assumed to be in alogshasrecentlyarisenduetotheadventoflargescalesur- equilibrium. Analyticfunctionsfordensitydistributions,the veysinastronomy,e.g.,futuresurveyslikeLSSTandGAIA age/metallicity relation and the IMF are provided for each haveplanstomeasureover1billionstars. Inordertogener- population.Asimilarschemeisalsousedbythephotometric ateasyntheticcatalog,onefirstneedstohaveamodelofthe codeTRILEGAL(Girardietal.2005). MilkyWay. Whilewe arefarfroma dynamicallyconsistent In spite of its popularity, the current Besançon model has model,aworkingframeworkisfundamentaltoprogress. In- important shortcomings. A web interface exists to generate evitably,thiswillrequireapproximationsorassumptionsthat syntheticcatalogsfromthe modelbutithas limitedapplica- maynotbemutuallyconsistent. Cosmologistsalreadyaccept bility for generating wide area surveys. Discrete step sizes suchcompromiseswhentheyrelatethe observedgalaxiesto forradial,andangularcoordinatesneedtobespecifiedbythe the dark-mattertest particlesthatemergefromcosmological userandresultsmightdifferdependinguponthechosenstep simulations. size. Thescaleheightandthevelocitydispersionofthedisc Therehavebeenvariousattemptsoverthepastfewdecades areinrealityafunctionofagebut,duetocomputationalcom- 1LeverhulmeVisitingProfessor,andMertonCollegeFellow,University plexity,thediscismodeledasafinitesetofisothermaldiscs ofOxford,OX13RH,UK. of different ages. Increasing the number of discs enhances 2 thesmoothnessofthemodelbutatthepriceofcomputational priori.Thiscanbeexpressedingeneralas cost(Girardietal.2005). f =f (r,v,τ,Z,m), (1) Inadditiontothedisccomponents,onealsoneedsamodel j j of the stellar halo. Under the hierarchical structure forma- j being the label of the component. An accurate form of tion paradigm, a significant fraction of the stellar halo is Equation (1) that describes all the properties of the Galaxy thoughttohavebeenproducedbyaccretioneventsandsigna- andisselfconsistentisstillanopenquestion. However,over turesoftheseshouldbevisibleassubstructuresinthestellar thepastfewdecadesconsiderableprogresshasbeenmadeto halo. Missions like GAIA, LSST and PanSTARRS are be- arrive at a working model using a few simple assumptions ing plannedthatwill enableus to detectsubstructuresin the (Robin&Creze 1986; Bienaymeetal. 1987; Haywoodetal. stellarhalo. 1997a,b; Girardietal. 2005; Robinetal. 2003). Our analyt- AsmoothanalyticstellarhaloasintheBesançonmodelis ical frameworkis based upon these models and we describe inadequatefortestingschemesofsubstructuredetection.Fur- thisbelow. thermore, such a halo does not accommodate known struc- For a given galactic component, let the stars be formed tures like the Sagittarius dwarf stream which may consti- at a rate of Ψ(τ) and the mass distribution of stars ξ(m,τ) tute a large fraction of the present halo (Ibataetal. 1995; (IMF)beaparameterizedfunctionofageτ only.Also,letthe Chouetal. 2010). Substructures have complex shapes and presentdayspatialdistributionofstars, f (r,τ)beafunc- pos hencetomodelthemwecannotusetheapproachofanalytic tionageonly. Finally,assumingthatthevelocitydistribution density distributions as discussed earlier. However, N-body f (v,r,τ)andthemetallicitydistribution,f (Z,r,τ)area vel Z models are ideally suited for this task. Brownetal. (2005) functionofageandpositiononlyonearrivesatamodelofthe attemptedtocombineasmoothgalaxymodelwithsomesim- form(labeljomittedforbrevity) ulated N-body models of disrupting satellites, but the stel- larhalowasnotsimulatedina propercosmologicalcontext. f =Ψ(τ)ξ(m,τ)f (r,τ)f (v,r,τ)f (Z,r,τ) (2) UsinghybridN-bodytechniques,Bullock&Johnston(2005) j m pos vel Z have produced high resolution N body models of the stellar h i Note, the functions on the right hand side can potentially halo that are simulated within a cosmological context; see take different forms for different galactic components, e.g., alsoCooperetal.(2010b)andDeLucia&Helmi(2008)fora the thin disc, the thick disc and so on. The IMF here is similarapproach.Thesecanbeusedtomakeaccuratepredic- normalized such that mmaxξ(m,τ)dm = 1 and m = tΛioCnDsMofpthaerasduibgsmtr.ucHtuorweseviner,thaesshteigllhalrighhatleodabnydBalrsoowtnesetttahle. mmaxmξ(m)dm is themmmeinanstellar mass. We nowhdisicuss mmin R (2005)thereareseveralunresolvedissuesrelatedtosampling thefunctionalformsofthemetallicityandthevelocitydistri- ofanN-bodymodelandthishaspreventedtheirwidespread bRution. use. Thedistributionf ismodeledasalog-normaldistribution, Z Theaimofthispaperistopresentfastandaccuratemeth- ods to convertanalytic and N-body models of a galaxy into fZ = 1 e−(logZ−logZ¯(τ))/(2σl2ogZ¯(τ)), (3) asyntheticcatalogofstars. Thiswouldrelievetheburdenof σlogZ¯(τ)√2π generating catalogs from modelers on one hand and on the the mean and dispersion of which is given by an age- other hand would allow the testing of models generated by dependent function. The mean metallicity as function of different groups. This should also facilitate rapid testing of age,Z¯(τ), ispopularlyknownastheage-metallicityrelation newmodels. (AMR).Ingeneralaspatialdependencemightalsobeadded Wepresentanewschemeforsamplingtheanalyticalmod- toEquation(3),e.g.,dlog(Z¯)/dR=constant. elsthatenablesusto generatecontinuousvaluesofthevari- The velocity distributionf can be modeled as a triaxial ableslike positionand age ofstars. Instead of a set of discs vel Gaussian, at specified ages, our methodology allows us to generate a disc that is continuous in age. As a concrete example, we 1 v2 v2 usetheBesançonanalyticalmodelforthedisc. Tomodelthe fvel=σ σ σ (2π)3/2exp −2σ2(Rτ,R) − 2σ2(τz,R) × disc kinematics more accurately, we employ the Shu (1969) R φ z (cid:18) R z (cid:19) distributionfunctionthatdescribesthenon-circularmotionin (v v (R) v (τ,R))2 the plane of the disc. For the stellar halo we use the simu- exp φ− circ − ad (4) lated N-bodymodels of Bullock&Johnston (2005) that can − 2σφ2(τ,R) ! reproducethesubstructureinthehalo.Weshowaschemefor where,R,φ,zarethecylindricalcoordinates,v (R)thecir- samplingtheN-bodyparticlessuchthatthesampledstarspre- circ cularvelocityasafunctionofcylindricalradiusRandv the servetheunderlyingphasespacedensityofN-bodyparticles. ad asymmetricdrift, whichis givenbythe Stromberg’srelation (equation4.228inBinney&Tremaine2008) 2. METHODS v (τ,R)= σR2 ad 2.1. Analyticframeworkformodelingthegalaxy 2vc × Wefirstdescribetheanalyticframeworkthatisusedbyus dlnρ + dlnσφ2 +1 σφ2 +1 σz2 (5) formodelingtheGalaxy. ThestellarcontentoftheGalaxyis dlnR dlnR − σR2 − σR2 ! modeled as a set of distinct components, e.g., the thin disc, the thick disc, the stellar halo and the bulge. The distribu- Alternatively,onecanmodelthe distributionofv usingthe φ tionfunction,i.e.,thenumberdensityofstarsasafunctionof Shu(1969)distributionfunction(Scho¨nrich&Binney2009), position (r), velocity(v), age (τ), metallicity (Z), and mass thisisdescribedinAppendixA. ThedispersionsoftheR,z (m)ofstarsforeachcomponentisassumedtobespecifieda andφcomponentsofvelocityincreaseasafunctionagedue 3 tosecularheatinginthediscandmoreoverthereisalsoradial densityregionsthatrequirerejectionofa lotofstars. More- dependence,the dispersionis higherclose to the center. We over,thecomputationaltimedoesnotscalewithoutputsam- modelthese effects as in Binney (2010) using the following plesize, sotogenerateevenasmallspecificsampleofstars, functionalform e.g.,starslyingwithinagivencolormagnitudelimitsorstars lyinginsomespecificregionofspace,onehastogenerateall σthin (R,τ)=σ τ +τmin β Σthin(R) q(6,) thestarsinthegalaxy. R,φ,z R0,φ0,z0 τ +τ Σthin We now discuss an efficient technique to sample a spec- (cid:18) max min(cid:19) (cid:18) ⊙ (cid:19) ified multidimensional distribution. We first divide the en- where Σ is the surface density of the disc (age dependence tire domain into bins that extend perpendicular to the time fromAumer&Binney2009). axis. Thenwesubdividethesebinsintosmallersub-domains, The modelas givenby Equation(2) has some limitations. calledleafnodes,withaspatialoct-tree(eachsubnodehaving For example, Equation (2) by itself is not dynamically self 1/8ththe volumeof its parentnode). Havingsubdividedthe consistentandthishastobeimposedexternally. Foragiven system,thevonNeumannrejectiontechniqueisnowapplied vcirc(R),Ψ(τ) andσz(τ) the verticalstructureofthediscin individually to each of the nodes. This has two immediate thesolarneighborhoodcanbemadedynamicallyselfconsis- advantages: tentbyusingtheapproachofBienaymeetal.(1987)(seealso Just&Jahreiß 2010). This leads to a constraint of the form Dependinguponthegivengeometryofthesurvey,one • ǫ = ǫ(τ), ǫbeingtheellipticityofahomoeoiddisc(orscale can check if the boundariesof the node intersect with heighth =h (τ)foradoubleexponentialdisc). thatofthesurveyand,ifnot, onecanskipthegenera- z z The second limitation of Equation (2) is that it does not tionofstars. have an explicit dependenceon the birth radii of a star. Re- cently semi-analytic models of the Galaxy that treat chem- Foracceptingandrejectingstarsitisnowpossibletoset ical evolution with radial mixing have been proposed by • themaximumdensitygmax(r,τ,l),lbeingavectorrep- Scho¨nrich&Binney(2009)wherethepropertiesofastarare resentingthelengthofthesidesofthenode,adaptively linked to their birth radii. A possible way to accommodate for each node. Since the variationin stellar densityin suchmodelsistointroduceinEquation(2)anadditionalpa- any node will be limited, g/gmax will never be very rameter R , the radius at which a star is born in cylindrical small,withtheconsequencethattherejectionsampling i coordinates.Themodeldistributionfunctionisthengivenby methodwillnotbeinefficient. Ψ(τ,R ) Westillneedtodecideonanoptimumtruncationcriterion i f= ξ(m,τ) tostopthesplittingofnodes. Anidealtruncationcriterionis m × h i tohavetheleastdensityvariationwithinanodebutthiswould fpos(r,Ri,τ)fvel(v,r,Ri,τ)fZ(Z,Ri,τ), (7) resultinalotofnodeswithnegligiblenumberofstarsbutan f (r,R ,τ) being the present day spatial distribution of otherwisestrongnumberdensitygradient. Instead,weadopt pos i a criterion that strives to achieve a fixed number of stars in starsthatwerebornτ yearsagoatradiusR . i a node similar to the one used in adaptive mesh refinement schemes. ForatotalofN stars,wesettheresolutionlimit 2.2. Samplingananalyticmodel:adaptivevonNeumann tot of a node as N /106, a choice that effectively generates rejectiontechnique tot about 106 nodes. However, calculating the number of stars Havingspecifiedtheanalyticalframeworkwenowdiscuss inanodeinvolvesanintegralofthenumberdensityoverthe theschemetosamplestarsfromamodelspecifiedwithinthis volumeof thenode, whichcan becomputationallyintensive framework.Inthepresentpaperwerestrictourselvestomod- whendoneforalargenumberofnodes. Hence,forthenode elsofthe formofEquation(2)buttheschemewe discussis splittingpurpose,wecomputethenumberofstarsas equallyapplicabletoextensionsoftheformofEquation(7). For a model given by Equation (2), if r and τ are given, Nn′ode =gmax(r,τ,l)l1l2l3l4 (9) sampling the metallicity distribution and the velocity distri- whereg (r,τ,l)is an analyticdensityfunctionrepresent- butionistrivial. Similarly,samplingthemassofastarisalso max ing the maximum density within a node and l ,l ,l ,l are straightforwardasitisa one-dimensionaldistribution. Once 1 2 3 4 thelengthofthesidesofthenode. Thefinalnumberofstars theage,metallicityandmassofastarareknown,asynthetic in a node, N , is calculated by numerically integrating libraryofisochronescanthenbeusedto generatestellar pa- node thenumberdensityoverthenodewithincreasingspatialand rameterslikecolor,magnitude,gravityandsoon. Sothekey temporalresolutiontillthenumberconvergeswithanuncer- taskistosamplethepositionandagecoordinatesfromaspec- taintylessthan0.1√N . Theaboveconvergencecriterion ifiedanalyticalfunction node roughlycorrespondsto1/10thoftheexpectedPoissonnoise. Ψ(τ)f (r,τ) Tosummarize,ourschemeforgeneratingasyntheticcata- g(r,τ)= pos , (8) logof starsusing theadaptivevonNeumannsamplingalgo- m h i rithmisasfollows. Foreachnode,wefirstdetermineNnode, whichismultidimensional. the number of stars that need to be generated. Next, to as- Asimpleapproachtosamplepointsfromananalyticalmul- signpositionandagetostars,werandomlyselectrandτ and tidimensionaldistribution as given by Equation (8) is to use then accept or reject them using the von Neumann scheme. thevonNeumannrejectiontechnique—generateastarhav- Next, the mass m for a star is sampled from the IMF ξ(m) ing random coordinates r,τ; then generate a random num- andZfromf (Z,r,τ). Thenearestisochronecorresponding Z ber x between 0 and g and finally acceptthe star if x < to(Z,τ)isselectedandthestellarparameterscorresponding max g(r,τ). Thisnaiveapproachis computationallyveryexpen- to m are determined. The velocities are assigned according sive. Alotofcomputationaleffortisspenttosamplethelow tothedistributionfunctionf (v,r,τ),whichisgenerallyin vel 4 theformofanellipsoidorcanbesubjectedtomoreadvanced fixedstepsize(astepsizeof0.1magindistancemodulus). treatment,anissuethatwerevisitlater. Finally,weusease- lection functionto decide if a star entersthe final catalog of 2.4. SamplinganN-bodymodel the survey. The processof spawningstars isrepeatedforall An N-body model consists of a finite set of particles dis- nodestogeneratethefullcatalog. tributedinphasespace(r,v).LetfN-body (r,v)bethephase During the star spawning process we do not take into ac- phase-space space distribution function describing this. An N-body par- count the effect of photometric errors or extinction. This is ticle in general represents a collection of stars rather than becauseeachobservationalsurveywillhavedifferentphoto- metricerrorlaws,anditisnotpossibletoformulateageneral a single star. Let mN-body be the mass of the stellar pop- scheme for this. Also, there is no perfect 3d model for ex- ulation corresponding to the N-body particle and this rep- tinction. Although,wedoprovideadefaultextinctionmodel resents the mass of stars that were ever formed (summed (seeSection2.6),butausermightbeinterestedintryingout over all ages). Also let us assume the age, metallicity and other models too. However, photometric errors and extinc- mass distribution to be given by Ψ(τ), fZ(Z,τ), and ξ(m) tion can be easily handledas a postprocessingstep and this respectively. The distributions being normalized such that istheapproachthatwefollow. Toaccomplishthis,first,one Ψ(τ)fZ(Z,τ)ξ(m)dZdτdm = 1. The distribution func- hastogenerateacatalogwithcolormagnituderangesslightly tionofstarsisthengivenby R largerthan the desired rangesand then for each star add the f∗ =fN-body (r,v) expectedextinctionandsubsequentlythephotometricerrors. phase-space × Finally, using these transformed magnitudes one can select Ψ(τ)f (Z,τ)ξ(m)mN-body/ m (11) Z thestarslyingwithinthedesiredcolormagnitudelimits. h i where the simulations provide the spatial and kinematic 2.3. Speedingupthestarspawningprocess behavior represented by the phase space density term fN-body (r,v). Theprocessofspawningstarsfromanode,asdiscussedin phase-space theprevioussection,canbefurtheroptimizedandwediscuss Foran N-bodymodel, the factthat the numberof N-body thisbelow. Tobeingwith,thedistributionofstellarmassesis particlesisfiniteleadstothefollowingproblem—ifthenum- such thatthere are a lot of low mass stars butveryfew high ber density of visible stars at a given point in space is more massstars. Lowmassstarsarealso lowinluminosityhence thanthenumberdensityofN-bodyparticlesatthesamepoint for a magnitude limited survey only nearby low-mass stars oneneedstooversampletheparticles. Theconditionforthis arevisible. Ontheotherhand,thevolumeexploredincreases isgivenby asthecubeofthedistance. Henceformostofthenodesthe mmax low mass stars do not enter the final catalog. We use this dτdZ f∗dm>fN-body (r,v), phase-space informationtooptimizeourscheme—foreachnodewefirst Z Zmmin(r) determine the lowest stellar mass mnode that can generate a orinotherwords min visible star and then exclusively generate stars having mass above this limit. The number of stars from a node that fall mN-body ξ(>m (r)) >1, (12) min intoasurveyisthengivenby m h i Nvniosde(mnmoidne,mmax)= mmaxξ(m)dm g(r,τ)d3rdτ mahmeilnio(rc)enbteriicngditshteanmceinoimfru.mTmypaiscsalolyf asusctahrathsaittuiastivoinsiabrlieseast Zmnmoidne Znode whenrissmall,i.e.,inregionsclosetothesun.Henceunless =Nnodeξ(>mnode) (10) onesimulatestheN-bodymodelwiththenumberofparticles min equalto the numberof stars in the modelgalaxy,one hasto We use stochastic rounding to convert Nnode to an integer live with oversampling. The question now is how are we to vis —i.e., ifthe fractionalpartofNnode isless thana Poisson- assignpositionsandvelocitiestotheoversampledparticles? vis distributedrandomnumberwitharangebetween0and1,we Thepositionsofoversampledparticlesshouldbeassigned incrementtheintegralpartby1. in such a way that they preservethe underlyingphase space We now describe the procedure to calculate mnode. If density. Forthiswefirstreviewapopularschemeofdensity min r is the heliocentric distance of a node, having sides of estimation known as the kernel density estimation, which is length l , then calculating mnode requires computing the alsostatisticallyquiterobust.Inthisschemethenumberden- minimumnodme ass of a star that wmilinl be visible at a distance sity at a given point x for an N-body system of particles is r √3l foreachisochroneandthencomputingtheglobal givenby node m−inimum, i.e., over isochrones of all ages and metallicities. f(x)= 1 K rj (13) Doing this individually for each node might be computa- hn h tionally expensive. However, it is possible to optimize this. Xj j (cid:18) j(cid:19) First, for a given galactic component, we identify the set of where K(u) is a kernel function such that its integral over isochronesthathaveanon-zeroprobabilityofbeingsampled all space is unity, r is the distance of the j-th particle from j and compute mnode over these set of isochrones only. Sec- x, and n the dimensionality of the space. In a multivariate min ondly, we note that for a given set of isochrones mnode is setting,r isgivenby min j amonotonicallyincreasingfunctionofheliocentricdistance. Since,wecanchooseanyvalueofmnodethatissmallerthan r = (x x )TΣ−1(x x ) (14) min j − j j − j thesmallestvisiblemass, weaccessthenodesinasequence q sortedbytheirminimumheliocentricdistance,r √3l , whereΣ−1isthematrixrepresentingthedistancemetricand node j and recompute mnode whenever the distance ch−anges by a isnormalizedsuchthatdet(Σ ) = 1. Thesmoothinglength min j 5 TABLE1 GEOMETRYOFSTELLARCOMPONENTS.THEFORMULASUSEDAREFROMROBINETAL.(2003).NOTE,(R,θ,z)ARETHECOORDINATESINTHE GALACTOCENTRICCYLINDRICALCOORDINATESYSTEMANDa2=R2+ z−zwarp (FORTHETHINDISC). kflareǫ(τ) Component Age(Gyr) densitylawρ(r,τ) ThinDisc ≤0.15 kρflcaΨre(ǫτ(τ)){exp(−(a/hR+)2)−exp(−(a/hR−)2)} where: hR+=5000pc,hR−=3000pc IMF-ξ(m) m−1.6form<1M⊙andξ(m) m−3.0form>1M⊙ ∝ ∝ ThinDisc 0.15–10 ρcΨ(τ) exp( (0.52+ a2 )) exp( (0.52+ a2 )) kflareǫ(τ){ − h2R+ − − h2R− } where: hR+=2530pc,hR−=1320pc, IMF-ξ(m) m−1.6form<1M⊙andξ(m) m−3.0form>1M⊙ ∝ ∝ Thickdisc 11 if|z|≤xl, ρcδ(τ−11)exp(−R−hRR⊙)×(1− xl×(21./+hxzl/hz) ×z2) if|z|>xl, ρcδ(τ−11)exp(−R−hRR⊙)× e1x+p(xxll//2hhzz)exp(−|hzz|) where: hR=2500pc,hz=800pc,xl=400pc IMF-ξ(m) m−0.5 ∝ Spheroid 14 ρcδ(τ−14)(cid:16)MaxR(⊙ac,a)(cid:17)nH where: a2=R2+ z2, ǫ2 aIMcF=-ξ5(0m0)pc,ǫm=−00..654,nH=−2.77 ∝ Bulge 10 ifpx2+y2<Rc, ρcδ(τ−10)exp(−0.5rs2) ifpx2+y2>Rc, ρcδ(τ−10)exp(−0.5rs2)×exp(−0.5(√x2+0.y52−Rc)2) where: r2= [( x )2+( y )2]2+( z )4, s q x0 y0 z0 Rc=2.54,x0=1.59,y0=z0=0.424,α=78.9◦,β=3.5◦,γ=91.3◦ IMF-ξ(m) m−2.35 ∝ ISM ρcexp(−R−hRR⊙)×exp(−|hzz|) where: hR=4500pc,hz=140pc Darkhalo ρc (1.+(a/Rc)2) where: Rc=2697pcandρc=0.1079 hofaparticleisdeterminedasthedistancetothek-thnearest BiDschemewiththenumberofnearestneighborsksetto64, neighborfromthelocationoftheparticle.Thephysicalinter- whichwefindisadequateforreproducingasmoothdistribu- pretationofthedensityschemeisthatifeachpointisassumed tion of stars in the phase space. Too small a value leads to tobeahyper-ellipsoidalballwithmassdistributedaccording lumpinessaroundthe N-bodyparticlesandtoolargea value tothekernelfunctionKandhavingacovarianceofh2Σ,then erodesthepositionalaccuracyofthedata.Weassumethema- theunderlyingsmoothdensityfieldissimplytheresultofsu- trixΣtobediagonalandusethecubiccellschemeofEnBiD, perposition of such balls. This suggests that if the spawned i.eΣ =Σ =Σ andΣ =Σ =Σ (dimensions1,2 11 22 33 44 55 66 starsarealsodistributedaroundtheparticleinasimilarfash- and3representingpositioncoordinateswhilethedimensions ion, i.e., distributed according to the kernel function K and 4,5and6thevelocitycoordinates).TheEnBiDalgorithmwas havingacovarianceish2Σ,thentheywillessentiallybesam- applied individuallyto each of the disrupted satellites in the plingthedensityfieldoftheN-bodysystem. simulations. Thequestionthatremainsishowtodeterminetheoptimum In order to assess the usefulness of EnBiD for our partic- covariancematrixΣ? Thesimplestschemeistouseaglobal ular application, we plot in Figure 1 (lower right panel) the metric based on the variance of the data along each dimen- distribution of the velocity scale factor Σ /Σ , which 11 44 sion. A globalmetric, however,isinsufficientforaccurately is the ratio of position to velocity scale and has units of calculating the phase space densities. What is needed is a kpc/kms−1, as estimated by EnBiD fopr the N-body parti- locallyadaptivemetricthatchangeswiththelocalconfigura- cles from one of the simulations of a disrupted satellites by tion of the data at each point in space (Sharma&Steinmetz Bullock&Johnston(2005) (see Section 3.4). It can be seen 2006; Sharma&Johnston 2009). Such a metric can be cal- that the scale factor varies by about 3 orders of magnitude, culatedusingthepubliclyavailablemultidimensionaldensity which demonstratesthe inappropriatenessof a globalmetric estimationcodeEnBiD by Sharma&Steinmetz (2006). En- andtheneedforalocallyadaptivemetricscheme. Thehuge BiD, which is an improvement over a scheme proposed by variationisduetothefactthatduringphasemixingthephase Ascasibar&Binney(2005),makesuseofabinaryspacepar- space density of particles is approximately conserved. This titioningschemealongwith the use ofthe informationtheo- meansthathighdensityregionsinpositionspacewill corre- retic concept of Shannon entropy to handle the issue of de- spondtolowdensityregionsinvelocityspaceandviceversa. termining the optimal locally adaptive distance metric Σ−1 ThiscanbeeasilyseeninthepanelsofFigure1,wherescat- j ina multi-dimensionalspace. Forouruse, we adopttheEn- ter plots of particles in position and velocity space has been 6 FIG.1.—Distributionofthevelocityscalefactorinunitsof kpc/kms−1 foroneofthedisruptedsatellitesfromBJ05simulationsanditsdependence onlocationinphasespace. Thebottomleftpanelshowsthedistributionof velocityscalefactorwhiletheotherpanelsshowthescatterplotofparticles inpositionandvelocityspace. Theparticleslyingintheinner(r<40)and outerregions(r >80)ofthehaloarecoloredredandgreenrespectively(r beingtheradialdistancefromthecenterofthestellarhalo). Highdensity regionsinpositionspacearelowdensityregionsinvelocityspaceandvice versa. FIG.2.— Positionandvelocityscatterplotsofadisruptedsatellitesampled shown. The highdensityregionsin positionspace (r < 40) fromanN-bodysimulationofBullock&Johnston(2005).Theupperpanels showtheoriginalNbodydatawhiletherestofthemareforsampleddata. arecoloredredwhilethelowdensityregionsinpositionspace Thesamplingfractionfsampleislabeledonthepanels.Thepanelsinsecond (r > 80) are colored green. The mean value of scale fac- rowuseanimprovedschemethatdispersesthespawnedstarsinphasespace onlywhenmorethanonestarisspawnedbyanN-bodyparticle. Note,the tor Σ /Σ was found to be 0.057 and 2.9 for the red h 11 44i stellarmassassociatedwithanN-bodyparticleisnotnecessarilyconstant. andgreenregions(seelowerrightpanelfordistributions).As a chepck we computed the scale factor using variance of the particles (σ2+σ2+σ2)/(σ2 +σ2 +σ2 ). The values orequaltoonestarneedstobespawnedfromtheN-bodypar- x y z vx vy vz ticle,onecaneliminatethenoisebysimplyassigningtheco- obtainedqwere 0.037 and 2.3 for red and green particles re- ordinatesoftheparticledirectlytothespawnedstar. Butitis spectively, which are in agreement with results reported by notpossibletoknowaprioriastohowmanystarsfromapar- EnBiD. ticlewillenterfinalcatalog.Toovercomethisweimplementa Summarizing, the scheme for sampling N-body distribu- schemewhichmakessurethat,amongthestarsspawnedbya tions is as follows. Select a particle i, determine if the dis- particleifoneofthestarscanbereassignedthecoordinatesof persing volume, defined by the kernel function K and the itsparentparticlewithoutitbeingexcludedfromthefinalcat- smoothinglengthh , intersects with the surveyvolume, and i alogwedoso. Theeffectofthisimprovedsamplingscheme calculate the number of visible stars it spawns, i.e., Ni = vis is shown in Figure 2 which shows the position (left panels) mN-bodyξ(> m (r))/ m . We use stochastic rounding as min and velocity (right panels) scatter plot of a disrupted satel- h i discussedearliertoconvertNviis toaninteger. Next,foreach litesystemfromanN-bodysimulationofBullock&Johnston spawned star we generate τ, Z and m according to the dis- (2005). The panels in the top row show the N-body parti- tributionfunctionsΨτ,fZ(Z,τ)andξ(m)anddispersethem cles in the simulation, while the other rows show the stars in phase space using the scheme discussed above. Finally, sampledfromit. Thesamplingfractionf ,i.e.,average sample weuseaselectionfunctiontodecideifastarentersthefinal numberofspawnedstarsperN-bodyparticle,isalsolabeled catalogofthesurvey. on each of the panels. For illustrative purpose we assume the numberof stars spawnedper unitmass of the particle to 2.4.1. ImprovingthepositionalaccuracywhilesamplinganN-body be constant but in a realistic simulation it will depend upon model the heliocentric distance of the particle. First, we compare Dispersing the stars spawned by an N-body particle in the first row with the second row that shows the results of phase space has its side effects. Thisaddsnoise to the data, the improved scheme for the case of f = 1. The dis- sample theeffectofwhichistoreducethepositionalaccuracy.When tribution of position coordinates is nearly identical for both aparticleisoversampledonecannotavoidthenoisebutwhen rowswhichiswhatwe expectwhenf = 1. On closer sample a particle is under-sampledor equally-sampledi.e, less than scrutiny,insomeregionsalossofpositionalaccuracycanbe 7 seen.Alsosomeregionsapparentlyseemtobeundersampled. mass stars using the bolometric correction (BC) tables from This departurefrom the ideal behavior is because the stellar http://stev.oapd.inaf.it/dustyAGB07/,which mass associated with an N-body particle obtained from the isadatabaseofBCtablesfordifferentphotometricbandspro- simulationsthatweuse,isnotnecessarilyconstant. Sosome videdbythePadovagroup(Marigoetal.2008;Girardietal. N-bodyparticlescanspawnmorethanonestarwhichsubse- 2002,2004). quentlyneedstobedispersedinphasespace. Next,wecom- Note, in the present version of the code we do not model paretheresultsoftheimprovedscheme(secondrow)withthe the white dwarfs. White dwarfs are quite faint (M > 10) V naivescheme(thirdrow)thatdispersesallstarsforthesame and rare comparedto other type of stars and hence for most valueoff =1. Itcanbeclearlyseenthattheimproved applicationsthey do notdominatethe star counts. However, sample schemehasbetterpositionalaccuracy. Forreference,thelast inapplicationswhereoneexpectstofindwhitedwarfscaution row shows the case of f = 4. We revisit this issue in shouldbeexercisedwheninterpretingtheresultsofGalaxia. sample Section 5.2 and Section 5.3 and discuss in greater detail the consequencesofourscatteringschemeforrealapplications. 2.6. Extinction Moststellarobservationssufferfromextinctionbydustand 2.5. StellarIsochrones this needs to be taken into account when comparing simu- Theoretical stellar isochrones are one of the main ingre- latedcatalogswithstellarsurveys. Althoughthetotalextinc- dients of our scheme. These are used to assign parameters tion along a given line of sight has been accurately mapped like luminosity, effective temperature, magnitude and color (Schlegeletal.1998),thesamecannotbesaidforthe3Ddis- to a star of a given age, metallicity and mass. We employ tribution of dust. Drimmeletal. (2003) provide a sophisti- Padova isochrones (Marigoetal. 2008; Marigo&Girardi catedmodelforthe3Ddustdistribution,andotheralternatives 2007; Girardietal. 2000; Bertellietal. 1994), which are includeMarshalletal.(2006). Sinceacomprehensiveextinc- availableforawidevarietyofphotometricsystems2 andun- tionmodelisbeyondthescopeofthispaper,weprovidehere likeotherpopularisochronesalsocovertheasymptoticgiant onlyasimpleextinctioncorrectionscheme. Note,extinction branchandtheredclumpphasesofstellarevolution. isbasicallyapostprocessingstepandonceanextinction-free Wenowdescribetheschemeforinterpolatingacrossdiffer- cataloghasbeengeneratedusingGalaxia, onecan subjectit ent isochrones. An isochrone, for a given age and metallic- totheextinctionmodelofone’schoice. ity, is a table of stellar parameters, e.g., magnitudes, effec- Theextinctionschemeweuseisarefinementofthemethod tive temperature and gravity, as a function of stellar mass. proposedbyBland-Hawthornetal.(2010b).The3Ddistribu- Hence in a given isochrone, stellar parameters for any in- tionofthedustismodeledasadoubleexponentialdisc, termediate mass can be easily computed by linear interpo- lation. However, interpolating across isochrones of differ- ρ0 R R⊙ z zwarp ρ (R,z)= exp − exp | − | . entagesandmetallicitiesinvolvestheuseofequivalentevo- Dust k − h × − k h flare (cid:18) R (cid:19) (cid:18) flare z (cid:19) lutionary points, which are not always provided with the (16) isochrones. Hence,forthesakeofflexibilityaswellascom- putational ease we adopt a simpler scheme. To begin with, wherezwarp(R,φ) andkflare(R) describethe warpandflare a grid of 182 34 isochrones spanning an age interval of ofthediscandaredescribedinSection3.2. Thevaluesofthe 6.6 < logt/yr×< 10.22(step size of ∆log(t) = 0.02) and parameterscontrollingthewarpandflarewereassumedtobe metallicity interval of 10−4 < Z < 3 10−2 (mean step sameasthatofthestellardisc.Next,weusetheSchlegeletal. sizeof∆log(Z) 0.072)werechosen.×Next,forastarofa (1998)E(B V)dustmapstoevaluatethebestfitparameters givenageandmet∼allicity,thenearestisochroneinthegridwas ofthedustm−odelandobtainedhR =4200pc,hz =88pcand identified. Subsequently,thisisochronewasusedtocompute ρ0 = 0.54mag/kpc. The fit parameters were obtained by thestellarparametersbylinearinterpolationinmass. Where minimizing needed, the metallicity Z was convertedto [Fe/H] using the Emap Emodel 2 relation χ2 = i − i (17) Emap [Fe/H]=log(Z/Z⊙); whereZ⊙ =0.019 (15) X(cid:18) i (cid:19) whereEmodel = ∞ρ (l ,b ,r)dristheextinctioninthe Intheabovementionedscheme,tosimulatethecolormag- i 0 Dust i i nitude diagrams accurately, the adopted grid should have a cellidefinedinthe(l,b)space. Thefitwasperformedinthe fineenoughresolution;andouradoptedresolutionwasfound range 15◦ < bR< 15◦ assuming R⊙ = 8.0 which defines − tobeadequateforthispurpose.However,forgreateraccuracy thedistancescale. ifdesiredonecanincreasetheresolutionofthegrid. Themodeleddustdiscisthenusedtogeneratea3Dextinc- The Padova isochrones are limited to stellar masses tionmapingalacticcoordinatesl,bandr, byintegratingthe greater than 0.15 M⊙. For lower masses extending up dustdensityalongdifferentlinesof sight. Theangularreso- to the hydrogen mass burning limit (0.07 M⊙ < m < lutionofthe3Dmapswasassumedtobe1024×1024while 0.15M⊙)weusetheisochronesfromChabrieretal.(2000). intheradialdirections512logarithmicallyspacedbinsinthe Since, these isochrones are only available for solar metal- range0.01<r <30kpcwereused. Valuesforanyarbitrary licity, we used them for all metallicities. Also, the positionwereobtainedfromthemapsbylinearinterpolation Chabrieretal. (2000) isochrones are only available for in(r,l,b). Alternately,onecanalsoassumeadirectionalde- Johnson-Cousin bands, hence, in order to able to sup- pendentnormalizationconstantρ0 suchthatthe totalextinc- port a variety of photometric bands, we choose to com- tionatinfinityalongagivenlineofsightequalsthevaluein pute absolute magnitudes for the above mentioned low theSchlegeletal.(1998)maps. Forthisweusehighresolu- tionSchlegeletal.(1998)mapshavinganangularresolution 2 The isochrones were downloaded from of 4096 4096 in l and b. We employ this as our default http://stev.oapd.inaf.it/cgi-bin/cmd scheme.× 8 TABLE2 TABLE3 AGEANDMETALLICITYDISTRIBUTION(MEANANDDISPERSION)OF VELOCITYELLIPSOIDOFSTELLARCOMPONENTS.NOTE,(R,φ,z)ARE GALACTICCOMPONENTS.THEVALUESSHOWNAREFROM THECOORDINATESINTHEGALACTOCENTRICCYLINDRICAL ROBINETAL.(2003) COORDINATESYSTEM. Age(Gyr) [Fe/H] σ[Fe/H] d[Fe/H]/dR σR0 σφ0 σz0 q β Thindisc 0-0.15 -0.01 0.12 kms/s kms/s kms/s 0.15-1 -0.03 0.12 Thindisc 50 32.3 21 0.33 0.33 1-2 -0.03 0.10 Thickdisc 67 51 42 0.33 0.33 2-3 -0.01 0.11 -0.07 Spheroid 141 75 75 0 0 3-5 -0.07 0.18 Bulge 110 110 100 0 0 5-7 -0.14 0.17 7-10 -0.37 0.20 for the thin and thick disc reproduce the Besançon results. Thickdisc 11 -0.78 0.30 0 Note, the radial dependence of velocity dispersions is mod- StellarHalo 14 -1.78 0.50 0 eled here using Equation (6) (values of q and β being from Bulge 10 0.00 0.40 0 Scho¨nrich&Binney 2009), rather than the Besançon model Note, Schlegel extinction maps are known to overesti- thatassumesdlnσR2/dR = −0.2kpc−1 andzeroderivative mate the extinction in the Galactic plane (regions where for other components. Note, the circular velocity profile is A > 0.5 mag), by about 30% (Arce&Goodman 1999; computedfromtheBesançonmassmodelandatthelocation CaVmbre´syetal. 2005). While Arce&Goodman (1999) oftheSunithasavalueof226.84kms−1. report results in a strip of few degrees in length only, Cambre´syetal. (2005) report mean values in galactic anti- 3.2. Warpandflare center hemisphere only. Since a proper recalibrated map is Thethinandthickdiscsareassumedtohaveawarpanda not yet available it is difficult to take these effects into ac- flare that is modeledfollowing the prescriptionof R03. As- count. In addition to directlyaffecting the total value of ex- suminggalactocentriccylindricalcoordinates(R,φ,z), stars tinctionalongaparticularlineofsight, theoverestimationis with radius R > R are displaced perpendicular to the warp also expected to affect the determinationof our scale height planebyanamount forthedustdisc,whichinturnwillaffectthevariationofex- tinctionwithdistance.Ifweassumetheoverestimationfactor zwarp(R,φ)=γwarpMin(Rwarp,R Rwarp)cos(φ φmax) − − to be a monontonicfunctionof extinction we expectthe ac- (19) tual scale height of the dust disc to be larger than what we whereφmax is the directionin which the warp is maximum. findhere. Hence,weadvisecautiontobeexercisedwhenus- For flaring, stars with R > Rflare have their scale heights ingourextinctionschemeinthegalacticplane. increasedbyafactor k (R)=1+γ Min(R ,R R ) (20) 3. AWORKINGMODELOFTHEGALACTICCOMPONENTS flare flare flare − flare Themainingredientsofourgalacticmodelasdescribedby Fortheparametersweadoptthesamevaluesasthatusedby Equation(2)arethestarformationrate(SFR),theagevelocity R03, namely φmax = 90.0, γwarp = 0.18kpc−1, Rflare = relation(AVR),theIMFandthedensityprofiles.Wenowdis- 1.12R⊙andγflare =0.0054kpc−1. cusstheadoptedfunctionalformsforeachofthesefunctions so as to arriveata workingmodelof the Galaxy. Instead of 3.3. Bulge doingadetailedmodelingwehereimplementthewellknown IthasbeenwellestablishedthattheMilkyWayhostsabar- and well tested Besançon model. The density distributions shapedbulge(Blitz 1993). With thearrivaloftheCOBE in- and the ages and metallicity of various galactic components frared maps, it has been possible to construct 3d models to are given in Table 1 and Table 2, these are the same as the fit the data. Various models using different triaxial analytic oneusedbyRobinetal.(2003)(R03hereafter).Thethindisc functionswerepresentedbyDweketal.(1995),ofwhichthe spansan ageintervalof 0 10Gyr. Onthe otherhand, the G2 model (boxy triaxial Gaussian type functions) provided − thick disc, the bulge and the halo are all assumed to have a the best fit to the data. The G2 bar is used by the Besançon constantage. modelandisadoptedbyusinourinitialdemonstrationofthe Galaxiamodel. TheG2densitydistributionhasacoreatthe 3.1. Thinandthickdisc centerdefinedbyradiusR andthe distributioniselongated c IntheBesançonmodel,theaxisratioǫandvelocitydisper- alongoneaxiswhichdefinesthemajoraxis. Theorientation sionsσR,σφandσz,forthethindisc,aretabulatedasafunc- ofthe majoraxisis definedbythreeanglesα,β andγ. The tionofageτ (velocityellipsoidbeingtakenfromGomezetal. values for these parametersare taken from R03, where they (1997).Hereweparameterizethemasfunctions.Theaxisra- obtained the values by fitting the model to the near infrared tioǫisparameterizedas starcountsoftheDENISsurvey. Next, we describe the kinematic properties of our ǫ(τ)=Min 0.0791,0.104 τ +τmin 0.5 , (18) model. Self-consistent dynamical models have been pre- τ +τ sented that use either the Schwarzschild technique Zhao (cid:18) max min(cid:19) ! (1996) orare extractedfromN-bodysimulations(Fux1999; and for velocity dispersions we use Equation (6). To match Athanassoula&Misiriotis 2002), but these are beyond the the Besançon values, the velocity dispersion for the thin scope of the current paper. Instead, we attempt to create a disc is assumed to saturate at 0.862σ and a value of simple working modelthat roughlysatisfies the existing ob- R0,φ0,z0 τ = 10.0andτ = 0.1is used. Thevaluesof param- servationaldataandwouldbeusefulforstudyingthesystem- max min eters used are summarized in Table 3. The adopted values aticeffectsthatcloudtheinterpretationofobservationdata. 9 Anumberofearlierstudieshaveclaimedthatthebulgeisin dark matter particles 3. The stellar masses depend upon the solidbodyrotationbyfittingastraightlinetotheplotofmean adopted star formation time scale and the later is chosen by radialvelocityasafunctionoflongitude(alsoreferredtoasa requiring that the simulated satellite stellar distributions re- rotationcurve),e.g.,byMenzies(1990)usingMiravariables producethestructuralpropertiesofLocalGroupdwarfgalax- andbyIzumiuraetal.(1995)usingSiOmaserstars. Aslope ies. Thethreemainmodelparametersofanaccretingsatellite ofabout10kms−1deg−1hasgenerallybeenreported,which arethetimesinceaccretion,t ,itsluminosity,L ,andthe acc sat translatestoapatternspeedofΩ=71.62kms−1kpc−1,as- circularity of its orbit, defined as ǫ = J/J (J being the circ suming8kpcasthedistancetothegalacticcenter. Recently, angularmomentumoftheorbitandJ theangularmomen- circ Howardetal. (2009) have also measured the rotation curve tumofacircularorbithavingsameenergy). Thedistribution using bulge M-giants in two strips at b = 4 and b = 8. ofthesethreeparametersdescribestheaccretionhistoryofa − − Theyreportthetworotationcurvestobenearlyidentical,sug- halo. To study the sensitivity of the properties of structures gestingthatthebulgeisrotatingcylindrically. Astraightline in the stellar halo to accretion history, additionally a set of fit(byeye)totheirrotationcurvealsoseemstosuggestavalue six artificial stellar halo models (referred to as non-ΛCDM closeto10kms−1deg−1. Hence,forsimplicityweadoptthe halos)were generatedbyJohnstonetal. (2008). These have value of Ω = 71.62kms−1 kpc−1 for the pattern speed of accretion events that are predominantly(i)radial(ǫ < 0.2), thebulge. Note,Howardetal.(2008)alsoreportthatthero- (ii)circular(ǫ > 0.7), (iii)old(t > 11Gyr), (iv)young acc tationcurvefortheb= 4strip,showsevidenceofflattening (tacc <8Gyr),(v)high-luminosity(Lsat >107L⊙),and(vi) beyond l >4. − low-luminosity(Lsat <107L⊙). | | Having specified the rotation we next review the velocity dispersions. Specifically,weattempttofixthevaluesforthe 3.5. Analyticsmoothstellarhalo velocity dispersions σR,σφ and σz of our model, expressed Forsomeapplicationsitisalsonecessarytohaveasmooth in cylindricalcoordinateswith respect to the galactic center. stellar halo. We do this by assuming an oblate power law A value of σr = 110.0kms−1 has typically been reported distributionwithellipticityǫ = 0.76andpowerlawindexof for the line of sight velocity dispersion in Baade’s window n = 2.44assuggestedbyR03.Thevaluesofǫandn can H H (Terndrupetal. 1995; Binney&Merrifield 1998) and this is bealte−redifdesired,e.g.,recentresultsbyJuric´ etal.(2008) what we assume for σR. Recent results by Howardetal. usingSDSSpredictaflatterhalo(ǫ=0.64)andsteepervalue (2009,2008)alsoconfirmthisbutshowavariationwithboth forthe powerlaw index(n = 2.77). In accordancewith H latitude and longitude– at b = 4, σr is found to increase theresultsofBondetal.(2010),−theprincipalaxesoftheve- from 70kms−1 to 110kms−1 w−ith l varyingfrom 10◦ to locityellipsoidareassumedtobealignedwithasphericalco- 0◦,whileatb= 8,σr isfoundtobe|re|lativelyconstantwith ordinatesystemandthevelocitydispersionsareassumedtobe avalueof70km−s−1. σ = 141kms−1 andσ = σ = 75kms−1. Thespheroid r θ φ Using proper motions, velocity dispersions along l and b isassumedtohavenonetrotation. havealsobeenmeasuredandananisotropyofσ /σ = 1.15 l b hasbeenreported(Rattenburyetal. 2007; Vieiraetal. 2007; 3.6. Initialmassfunctionandnormalizations Spaenhaueretal. 1992). The anisotropy could both be due Having specified the density functions, we now provide tointrinsicanisotropyaswellasduetotherotationalbroad- the the initial mass function ξ and the normalization con- ening as demonstrated by Zhaoetal. (1996). Since it is not stants like ρ and the star formation rate Ψ(τ). As in the c possibletoexactlydeducethevalueofσφfromσlweassume Besançon model, the IMF was assumed to be of a power for simplicity σφ = σR and choose the value of σz so as to law form, ξ(m) = m−α, the values of α are tabulated in reproducetheobservedvalueofanisotropyratioinBaades’s Table 1. The mass limits of the IMF were defined to be in window. For a choice of σz = 100kms−1, the bulge stars the range 0.07 < m < 100. Next, the density normal- lyinginthefield(l,b)=(1◦, 4◦)andwithheliocentricdis- ization constants were chosen to reproduce the local mass − tancebetween7 9kpc,werefoundtohaveσl/σb =1.17. density of current stars ρ for the various galactic compo- − 0 nents in the solar neighborhood (Table 2 in R03). This re- sulted in a star formation rate of Ψ = 2.85M⊙/yr and ρc 3.4. Simulatedstellarhalo of 1.55 109 and 1.31 107 M⊙kpc−3 for the thick disc × × andthestellarhalorespectively.However,forthethindiscin To make theoretical predictions of structures in the the regime 7 < τ < 10Gyr the SFR had to be lowered by stellar halo, we use the eleven stellar halo models of 20%inordertomatchthestarcountsintheBesançonmodel. Bullock&Johnston(2005),whichweresimulatedwithinthe This is because in this regime the isochrones employed by context of the ΛCDM cosmologicalparadigm. These simu- us differ from those in the Besançon model. For the bulge laastiNon-sbofdoyllopwarttihceleascycsrteetmiosnoonftoinadigvaildauxaylwsahtoelsleitedsismc,obduellgeed, ρdecn=sit3y.4o9f×131.7069Mst⊙arksppcc−−33winasascecloercdteadncwehwicihthgtihveesBaescaennçtroanl and halo potential is represented by time dependent analyt- model. ical functions. Semi-analytical prescriptions are used to as- For our analysis, we assume the Sun to be located at signastar formationhistoryto eachsatellite andaleakyac- a radial distance of 8 kpc from the center of the galaxy creting box, chemicalenrichmentmodelis used to calculate (Eisenhaueretal. 2003; Reid 1993) Since our adopted ra- themetallicityasafunctionofageforthestellarpopulations dial distance of Sun is different from the Besançon model (Robertsonetal.2005;Fontetal.2006). Thedarkmatterha- (8.5kpc), we need to recomputethe SFR. The re-calibrated los are assumed to follow an NFW density profile whereas, valueofSFRforouradoptedlocationoftheSunturnsoutto thestellarmatterisassumedtofollowaKingsprofileembed- ded within the dark matter halos. To simulate Kings profile 3Weactuallyusethemasslesstestparticleswhichhave10timesthephase within dark matter profiles the dark matter particles are as- spaceresolutionoftheoriginaldarkmatterparticle(seeBullock&Johnston signedastellarmasswhichisafunctionoftheenergyofthe (2005)forfurtherdetails). 10 FIG.4.—Comparisonofcolor-magnitudediagramobtainedfromtheHip- parcoscatalog(leftpanel)withthatofsimulationfromGalaxia(rightpanel). TheplotsshowstarsinthesolarneighborhooddefinedbyV < 8anddis- tancer<100pc. FIG.3.—StarcountpredictionsofGalaxiainvariousdirectionscompared withthoseoftheBesançonmodel.ForthedirectionalongtheNorthpolethe contributionsfordifferentgalacticcomponentsareshownseparately. be 2.37M⊙/yr. Also, to match the star counts for the stel- larhaloatlargedistances(alongthenorthgalacticpole)with those of the Besançon model, we had to increase the local massdensityofthestellarhaloby10%. 4. TESTSANDRESULTS Inthissectionwetestthecodebycomparingitspredictions withvariousobservationalconstraints. 4.1. ComparisonwithBesançon ThestarcountpredictionsoftheBesançonmodelhavebeen testedagainstavarietyofobservations,andsinceGalaxiauses the same disc model, it suffices for most situations to show thatGalaxiareproducestheBesançonmodel.Inordertokeep thetestsimpledustextinctionwasneglected. Totestthecor- respondenceofourresultswiththatofBesançon, weplotin Figure3thestarcountsasafunctionVbandmagnitudealong threedirections;thenorthgalacticpole,(l,b)= (90◦,10.0◦) and(l,b)= (270◦,10.0◦). Thelasttwodirectionswerecho- sentoillustratetheeffectofwarpingwhichresultsinabifur- cationinthestarcountsatfaintermagnitudes.Overallwefind goodagreementwiththecurvesfromGalaxiacorresponding closelytotheBesançoncurves. Forthedirectionalongthegalacticpole,thestarcountsfor each of the galactic components are also shown separately. It can be seen that the thick disc and the stellar halo are in perfectagreementbutthethindiscshowssomesubtlediffer- FIG.5.— Comparison of distributions of radial distance, absolute mag- nitudeandcolorforstarsobtainedfromtheHipparcoscatalogwiththatof ences. FortheGalaxiathindisc,thereisaslightexcessofin- simulationfromGalaxia. termediatemagnitudestarsandashortageofextremelyfaint magnitudestars. Thissameeffectcanbeseenintheothertwo stars, accurate distances can be obtained, which enables the directionsalsobutisshiftedtoevenfaintermagnitudes.These constructionofthecolor-absolutemagnitudediagramofstars. discrepanciesareduetothedifferencesintheisochronesem- TheHipparcosmission(ESA1997)producedanastrometric ployedbyusascomparedtothoseusedbyR03. Specifically data base of 117955 stars down to a limiting magnitude of forthethindisc in themassrange0.15 < m < 0.6, theab- V 12.4. For stars with V < 9, the data have a median solutemagnitudeofstarspredictedbythePadovaisochrones pre∼cisionofabout1maswhichisidealforanalyzingthesolar werefoundtobebrighterthanthatoftheisochronesusedby neighborhood. ThecompletenesslimitoftheHipparcoscat- theBesançonmodel. alogisestimatedtobebetween7.3-9maginVband. Hence we use the following criterion to extract a volume complete 4.2. ComparisonwiththeHipparcosdata samplefromtheHipparcoscatalog(vanLeeuwen2007,new The predictionsof stellar propertiesin the solar neighbor- reduction)—V < 8andparallaxπ > 10mas. Binarystars hood is an essential test for any galactic model. For these wereexcludedfromtheanalysis. TheHipparcosmagnitude,