Galaxia: a code to generate a synthetic survey of the Milky Way

1.6 MB·English

by  Sanjib Sharma

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
Besanc¸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 Besanc¸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 Besanc¸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-
AsmoothanalyticstellarhaloasintheBesanc¸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
usetheBesanc¸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 Besanc¸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 Besanc¸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; computedfromtheBesanc¸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 Besanc¸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 Besanc¸on
constantage.
modelandisadoptedbyusinourinitialdemonstrationofthe
Galaxiamodel. TheG2densitydistributionhasacoreatthe
3.1. Thinandthickdisc
centerdefinedbyradiusR andthe distributioniselongated
c
IntheBesanc¸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 Besanc¸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 Besanc¸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%inordertomatchthestarcountsintheBesanc¸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 Besanc¸on model. For the bulge
laastiNon-sbofdoyllopwarttihceleascycsrteetmiosnoonftoinadigvaildauxaylwsahtoelsleitedsismc,obduellgeed, ρdecn=sit3y.4o9f×131.7069Mst⊙arksppcc−−33winasascecloercdteadncwehwicihthgtihveesBaescaennc¸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 Besanc¸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
withthoseoftheBesanc¸onmodel.ForthedirectionalongtheNorthpolethe
contributionsfordifferentgalacticcomponentsareshownseparately.
be 2.37M⊙/yr. Also, to match the star counts for the stel-
larhaloatlargedistances(alongthenorthgalacticpole)with
those of the Besanc¸on model, we had to increase the local
massdensityofthestellarhaloby10%.
4. TESTSANDRESULTS
Inthissectionwetestthecodebycomparingitspredictions
withvariousobservationalconstraints.
4.1. ComparisonwithBesanc¸on
ThestarcountpredictionsoftheBesanc¸onmodelhavebeen
testedagainstavarietyofobservations,andsinceGalaxiauses
the same disc model, it suffices for most situations to show
thatGalaxiareproducestheBesanc¸onmodel.Inordertokeep
thetestsimpledustextinctionwasneglected. Totestthecor-
respondenceofourresultswiththatofBesanc¸on, weplotin
Figure3thestarcountsasafunctionVbandmagnitudealong
threedirections;thenorthgalacticpole,(l,b)= (90◦,10.0◦)
and(l,b)= (270◦,10.0◦). Thelasttwodirectionswerecho-
sentoillustratetheeffectofwarpingwhichresultsinabifur-
cationinthestarcountsatfaintermagnitudes.Overallwefind
goodagreementwiththecurvesfromGalaxiacorresponding
closelytotheBesanc¸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-
theBesanc¸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,