On the Run: mapping the escape speed across the Galaxy with SDSS

1.3 MB·

by  Angus A. Williams

Mon.Not.R.Astron.Soc.000,1–??(2016) Printed9January2017 (MNLATEXstylefilev2.2)
On the Run: mapping the escape speed across the Galaxy
with SDSS
Angus A. Williams1, Vasily Belokurov1, Andrew R. Casey1 & N. Wyn Evans1
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
7
AcceptedReceived;inoriginalform
1
0
2
ABSTRACT

n WemeasurethevariationoftheescapespeedoftheGalaxyacrossarangeof 40kpc
a inGalactocentricradius.Thelocalescapespeedisfoundtobe521+46kms−1,∼ingood
J −30
agreement with other studies. We find that this has already fallen to 379+34kms−1
5 −28
at a radius of 50kpc. Through measuring the escape speed and its variation, we ob-

tain constraints on the Galactic potential as a whole. In particular, the gradient in
]
A the escape speed with radius suggests that the total mass contained within 50kpc is
29+7 1010M , implying a relatively light dark halo for the Milky Way. Our method
G −5× (cid:12)
represents a novel way of estimating the mass of the Galaxy, and has very different
.
h systematics to more commonly used models of tracers, which are more sensitive to
p the central parts of the halo velocity distributions. Using our inference on the escape
- speed,wetheninvestigatetheorbitsofhigh–speedMilkyWaydwarfgalaxies.Foreach
o
dwarf we consider, we predict small pericenter radii and large orbital eccentricities.
r
t This naturally explains the large observed ellipticities of two of the dwarfs, which are
s
a likely to have been heavily disrupted as they passed through pericenter.
[
Key words: Galaxy: halo – galaxies: kinematics and dynamics – dark matter

1
v
4
4
4 1 INTRODUCTION authors used data from the RAVE spectroscopic survey,
1 which provided an abundance of information on the kine-
The fastest moving stars have long been a subject of fas-
0 matics of local stars. P14, who defined the escape speed as
cination and speculation. Vibert Douglas (1956) recounts a
. the minimum speed needed to reach three virial radii, con-
1 discussion between Kopff and Eddington on horse-racing,
cluded that it was 533+54kms−1 (90% confidence), albeit
0 an enthusiasm of the latter. Kopff stated he was not in- −41
7 on the basis of a small sample.
terested, because of course one horse will always run faster
1
than another. “But why”, retorted Eddington. “When one There are some obvious drawbacks. First, there is no
:
v starmovesfaster,youareveryinterested!”Thefastestmov- guaranteethatthehighvelocitytailofthedistributionfunc-
i ing stars are intriguing both because of the processes that tionisactuallyoccupiedallthewayuptotheescapespeed.
X
acceleratedthemandbecausetheyareprobesoftheirenvi- This may mean that the velocity of the fastest moving star
r
ronment. is an under-estimate of the true escape speed. Second, the
a
The local escape speed v is a measure of the depth method is sensitive to interlopers or contaminants, which
esc
ofthepotentialwellatthesolarpositionΦ(R ).Therefore, maybeunrepresentativeofasmooth,relaxedstellarpopula-
(cid:12)
thevelocitiesofthespeedieststarspassingthroughthesolar tion.ThisincludesstarsintheprocessofleavingtheGalaxy,
neighbourhood can in principle provide information on the such as hypervelocity stars ejected by interaction of bina-
enclosed mass. There were early investigations on the local ries with black holes (e.g., Brown 2015; Boubert & Evans
escape speed by Caldwell & Ostriker (1981) and Alexan- 2016). Stars may also be unbound from the Milky Way but
der (1982), who suggested values ∼450kms−1. The analy- nonetheless bound to the Local Group. Although no such
sis methods were extended and systematized by Leonard & stars are known, the phenomenon is familiar to us through
Tremaine (1990), who emphasised the importance of mod- the intergalactic stars identified in nearby clusters like For-
ellingtheshapeofthetailofthevelocitydistribution.They nax(Theuns&Warren1997).Third,andperhapsmostawk-
provided a number of physically-inspired models and con- wardly, the spatial distribution of the highest energy stars
cluded that v lay between 450 and 650kms−1 with 90% issetbythestochasticpatternsofcosmicaccretion.There-
esc
confidence. The recent spectroscopic surveys of the Galaxy fore,theconceptofasmoothvelocitydistributionmaybea
have stimulated a surge of activity, notably by Smith et al. fictionatthehighestenergies,evenatlocationsintheinner
(2007)andPiffletal.(2014)(hereafterS07andP14).These galaxy.
(cid:13)c 2016RAS
2 Williams, Belokurov, Casey & Evans
These disadvantages are offset by a number of assets. Thismodelhasbeenappliedtodataseveraltimes,most
First,onlythehighvelocitytailofthedistributionfunction notably by S07 and P14 using data from the RAVE survey
(DF) need be modelled, as opposed to its entirety. This is (Kordopatisetal.2013).Bothstudiesusedsmallsamplesof
clearly a massive simplification that sidesteps much of the stars (<100) close to the sun and found v ∼530kms−1.
esc
complexityinbuildingDFs.Second,themethodprovidesa We seek to extend their work by constraining the escape
nice counterpoint to other ways of measuring the potential speedoftheGalaxyatavarietyoflocations.Wedothisby
(orequivalentlythemass)oftheGalaxy.Forexample,meth- parameterisingtheescapespeedasafunctionofpositionx,
ods based on the Jeans equations use the first and second so that
moments of the distribution, and so are controlled by the (cid:40)
mainbulkratherthanthetail.And,third,futureprospects p(v |x)= C·(vesc(x)−|v|||)k+1 if vmin (cid:54)v|| <vesc(x),
||
arebright,withhugenewdatasetsofradialvelocitiesfrom 0 otherwise,
spectroscopic surveys and of proper motions from the Gaia
(3)
astrometric satellite becoming available.
All previous work has focussed on measuring the es- whereC isalocation-dependentnormalisationfactor,given
cape speed locally. This is because samples of high veloc- by
ity stars have been small (sometimes minute, for example k+2
S07 used just 16 stars, whilst P14 relied on 86 stars) and C = . (4)
(v (x)−v )k+2
esc min
concentrated in the solar neighbourhood. In this paper, we
By far the largest source of uncertainty in our analysis is
presentthefirstmeasurementsoftheescapespeedthrough-
in the distance to each star. Consequently, we consider the
out the Galaxy using a variety of tracers – main-sequence
uncertainty in the radial velocity, longitude and latitude to
turn-offstars(MSTOs),bluehorizontalbranchstars(BHBs)
be negligible. Our likelihood function should therefore be
andK-giants–extractedfromtheSloanDigitalSkySurvey.
the probability of a radial velocity, given Galactic coordi-
Although the MSTOs are located at heliocentric distances
natesandanimperfectinferenceonthedistancetothestar.
within ∼ 3 kpc, the BHBs and K-giants in our sample ex-
Writingx=((cid:96),b,s),where((cid:96),b)areGalacticlongitudeand
tend out to Galactoctocentric radii of ∼ 50 kpc, providing
latitude, respectively, and s is the measured line of sight
much greater reach.
distance to the star, we then have
Section 2 describes the likelihood function and mod-
els for the tail of the velocity distribution, primarily fol- (cid:90)
p(v |(cid:96),b,s)= p(v |(cid:96),b,s(cid:48))p(s(cid:48)|s)ds(cid:48). (5)
lowing the formalism established by Leonard & Tremaine || ||
(1990). We discuss our tracers and distance estimators in
p(s(cid:48)|s) is the probability that s(cid:48) is the true distance to the
Section 3. This includes the cuts used to build the samples
stargivenourimperfectinferences.Finally,wealsoinclude
of MSTOs, BHBs and K giants, as well as to extract the
a Gaussian outlier model for possibly unbound stars
high velocity stars. The choice of potential is given in Sec-
tion 4 and enables us to broaden the discussion of escape A −v 2
p (v )= √ exp || , (6)
speedintoenclosedmassandcircularspeed.Priorsandnu- out || 2πσ2 2σ2
merical implementation are described in Section 5, whilst
Section 6 presents our results, and discusses their implica- wherewefixσ=1000kms−1,andAisthenormalisationof
tionsforthreefastmovingGalacticsatellitegalaxies(Bootes theGaussianovertheinterval[vmin,∞].Wethenintroduce
III, Triangulum II and Hercules). a free parameter f for the fraction of outliers, so that the
overall likelihood function is
p (v |(cid:96),b,s)=(1−f)p(v |(cid:96),b,s)+fp (v ). (7)
tot || || out ||
2 METHOD
ThisEquationrepresentsthelikelihoodthatwewillusefor
Leonard & Tremaine (1990) proposed that the velocity dis- the remainder of the paper, while making specific choices
tribution of high-speed stars is a power law of the form for vesc(x). In practice, we compute the RHS of Equation
(5) using Monte-Carlo integration (e.g. Evans et al. 2016;
(cid:40)
p(v)∝ (vesc−v)k if vmin (cid:54)v<vesc, (1) Bowden et al. 2016), so that
0 otherwise.
N
1 (cid:88)
p(v |(cid:96),b,s)(cid:39) p(v |(cid:96),b,s ), (8)
vesc is the escape speed and vmin is a cut-off, such that || N || n
n=1
p(v < v ) begins to deviate from a power law. Since the
min
model depends only on the speed v, the full velocity dis- where each of the sn is drawn from p(s(cid:48)|s).
tributionfunctionisimplicitlyisotropic.Asradialvelocities
aremeasuredfarmorepreciselythantransversevelocities,it
isusefultomarginalisetheabovemodeloverpropermotion. 3 SAMPLE SELECTION
This gives
In order to measure the variation in the escape speed with
(cid:40)
(v −|v |)k+1 if v (cid:54)v <v , position in the Galaxy, we require a set of halo tracers that
p(v)∝ esc || min || esc (2)
spans a sufficiently large volume, and has measured radial
0 otherwise,
velocities and distances. To obtain such a sample, we use
wherev isGalactocentricradialvelocity.Sinceweareonly the SDSS 9th data release (Ahn et al. 2012). We choose to
||
interested in the absolute value of v in this work, we shall utilise three distinct sets of tracers: main sequence turnoff
||
now refer to |v | as v for brevity. (MSTO), K-giant and blue horizontal branch (BHB) stars.
|| ||
(cid:13)c 2016RAS,MNRAS000,1–??
Fast moving stars in SDSS 3
550
MSTO K-giant BHB
500 N =1573 N =343 N =44
450
1ms−400
k
/350
v||
300
250
200
0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50
r/kpc
Figure 1. The distribution in Galactocentric radius and v of each of our tracer samples. The MSTO stars probe a smaller range in
||
radius,butarenumerous,whereastheK-giantsprobeamuchlargerdistancerange,butarefewerinnumber.TheBHBsprobeasimilar
range to the K-giants, but form a far smaller sample. The black points in the MSTO panel highlight the 11 objects that are discussed
inSection6.
MSTO stars are numerous, and are mostly observed at dis- with a mother sample of 22071 MSTO tracers. In order to
tances ∼ 3kpc. These stars allow us to constrain the local drawsamplesfromp(s(cid:48)|s)(Equation5),wetakethequoted
escape speed, and its variation relatively close to the sun. uncertainties on g, r, i and [Fe/H] from the SDSS pipelines
K-giants and BHBs are fewer in number, but are bright andassumethattheyarenormallydistributedanduncorre-
and have been observed at distances ∼ 50kpc from the lated.WethendrawMonte-Carlosamplesfromeachofthese
sun, pushing the spatial extent of our sample to a range distributions and compute the distance for each sample.
of ∼ 40kpc in Galactocentric radius. Before selecting our
high-speed sample of stars, we first constructed a ‘mother
3.2 BHB sample
sample’foreachclassoftracerusingaseriesofcuts.Inaddi-
tiontothecutsdescribedbelow,wealsoremovedstarswith To obtain a clean sample of BHBs, we first select in the
latitudes |b|<20◦ in order to remove possible disc contam- de-reddened colour-colour box
inants, and stars with radii r>50kpc. We compute Galac-
−0.25<(g−r) <0,
tocentricradiibyassumingasolarradiusR =8.5kpc.The 0
(cid:12)
full SQL queries used are given in Appendix A. 0.9<(u−g)0 <1.4, (12)
which was used by Deason et al. (2011). There will still
be significant contamination from high surface-gravity blue
3.1 MSTO sample
stragglerswithinthisbox,sowethenmakethespectroscopic
To extract the mother sample of MSTO stars, we start by cuts
selecting in the de-reddened colour-magnitude box
−2<[Fe/H]<−1,
0.2<(g−r)0 <0.6, 3<logg<3.5, (13)
14.5<r/mag<20, (9) 8300<T /K<9300,
eff
where the r-band extinction A < 0.3. We then make cuts giving us 1039 BHBs in total. We estimate the absolute g-
r
on spectroscopic parameters, so that band magnitudes, and hence distances, of the BHBs using
the relation derived by Deason et al. (2011)
4500<T /K<8000,
eff
3.5<logg<4, Mg =0.434−0.169(g−r)0+2.319(g−r)20
−4<[Fe/H]<−0.9. (10) +20.449(g−r)3+94.617(g−r)4.
0 0
This gives us a sample of metal-poor MSTO stars. We also (14)
makecutstoensurehigh-qualityphotometryandradialve- We then use the same Monte-Carlo strategy as for the
locity measurements for the sample. To compute distances MSTO stars to sample from the distance uncertainties.
tothesestars,weestimatetheabsolutemagnitudeinther-
band using the prescription derived by Ivezi´c et al. (2008).
3.3 K-giant sample
With x=(g−i) , we have
0
δM =4.5−1.11[Fe/H]−0.18[Fe/H]2, Rather than gathering a sample of K-giants directly from
r
SDSS, as we did for our other two tracer samples, we use
M =−5.06+14.32x−12.97x2+6.127x3
r0 the catalog constructed by Xue et al. (2014), taken from
−1.267x4+0.0967x5, SEGUE. The SEGUE survey (Yanny et al. 2009) explicitly
targeted K-giants in three different categories: “I-colour K-
M =M +δM . (11)
r r0 r giants”, “red K-giants” and “proper-motion K-giants”. All
Using these formulae, we compute a point estimate of the three categories satisfy the photometric constraints
distancetoeachstarandremoveobjectswithimpliedhelio-
0.5<(g−r) <1.3,
centric distances > 15 kpc in order to remove objects with 0
0.5<(u−g) <2.5, (15)
spuriousabsolutemagnitudeestimates.Thisfinallyleavesus 0
(cid:13)c 2016RAS,MNRAS000,1–??
4 Williams, Belokurov, Casey & Evans
andhavemeasuredpropermotions<11masyr−1.Eachtar- distributionsofeachofourfinaltracersamplesinther−v
||
getcategorythenpopulatesdistinctregionsofthe(u−g)− planeisshowninFigure1.Wenotethatoursampleof1960
−(g−r) plane. The precise details of the targetting crite- starsrepresentsafactorof>100increasefromthe16stars
ria are given in Yanny et al. (2009). On top of the SEGUE used by S07 and a factor >20 compared to the 86 used by
selection, Xue et al. then make further cuts such that P14.
logg<3.5,
E(B−V)<0.25mag.
4 CHOICE OF v (x)
esc
Besides our cuts on latitude and Galactocentric radius, we Our primary model for the escape speed is a spherically–
make one further cut on this sample to ensure a clean halo symmetric power law model (SPL)
population, namely
(cid:18) r (cid:19)−α/2
[Fe/H]<−0.9. (16) vesc(r)=vesc(R(cid:12)) R , (17)
(cid:12)
After applying our latitude and metallicity cuts, there re- wherer isGalactocentricradiusand0(cid:54)α(cid:54)1,onphysical
main 5160 K-giants in our mother sample. Xue et al. com- grounds.Wenormalisethemodeltotheescapespeedatthe
pute posterior samples for the distance moduli, µ, of the solarradius.ThegravitationalpotentialΦandescapespeed
K-giants. Their table provides the maximum a posteriori are related by
valueofµforeachstar,percentilesoftheposteriorsamples, √
and a summary of the uncertainty, δµ. In order to compute vesc = −2Φ, (18)
distance samples for the K-giants, we simply assume that meaningthattheSPLmodelcorrespondstoagravitational
theuncertaintyinthedistancemodulusisGaussian,witha potential of the form
standard deviation of δµ.
v (R )2 (cid:18) r (cid:19)−α
Φ(r)=− circ (cid:12) . (19)
α R
(cid:12)
3.4 High-speed sample
Thus,unlikeotherstudiesthatmeasureonlythelocalescape
Givenourthreemothersamples,wenowseektoselectmem- speed, we can translate our measurements of vesc(R(cid:12)) and
bers of the high-speed tail of the velocity distribution. In α into inference of the circular speed of the potential, and
order to do this, we must make a decision for the value of the mass enclosed within spherical shells.
v in Equation (3). S07 did not have access to measure- We will also report on results when two other models
min
mentsofstellarparametersfortheirsample,andadopteda are used. The first is a simple generalisation of the SPL, a
very high threshold of v = 300kms−1 in order to avoid power law in elliptical radius (EPL)
min
contaminationfromtheGalacticdisc.P14usedpropermo- (cid:18) r (cid:19)−β/2
tionmeasurementstoremovestarswithdisc-likesignatures v (r )=v (R ) q , (20)
esc q esc (cid:12) R
ofrotation(althoughinpracticethiscutoperatesverysimi- (cid:12)
larlytoametallicitythreshold,seetheirFigure9),andthen withr =(cid:112)R2+z2/q2,whereRandz aretheusualcylin-
q
introduced a cut of vmin =200kms−1. Our cuts in latitude drical coordinates, and 0(cid:54)q(cid:54)∞. If the escape speed falls
andmetallicityservetoremovedisccontaminants,whichis off more rapidly with height above the Galactic plane than
muchlessofaconcernwhenusingmoredistantSDSStrac- it does in the radial direction, then q<1 (oblate), whereas
ers than it is for samples of RAVE stars. Since we expect larger escape speeds high above the plane suggest q > 1
very little disc contamination, our choice of vmin now de- (prolate). This model corresponds to a gravitational poten-
pends only on the range of velocities for which we believe tial of the same form as Equation (19), but with elliptical
the power law model in Equation (3) is valid. radius r replacing spherical radius r.
q
P14 performed a detailed analysis of the simulation The final model we consider is the truncated flat rota-
suite of Scannapieco et al. (2009) in order to make an edu- tion curve (TF) model (see Gibbons et al. 2014; Wilkinson
cated choice for vmin, and found that the distribution func- &Evans1999).Weprimarilyinvestigatethismodelinorder
tion of v|| in the simulations does not significantly deviate to facilitate a direct comparison with the mass inference of
from a power-law when v|| > 150kms−1. In our work, we Gibbonsetal.(2014).Themodelhasaflatrotationcurveof
must consider the fact that our cut must be appropriate amplitudev intheinnerparts,whichthenstartstodecline
0
across a range of locations. Since the vast majority of our as a power-law γ/2 around a scale radius r
s
sample are further from the Galactic center than the sun,
and physical reasoning suggests that the escape speed can- v2(r)= v02rsγ . (21)
notincreaseasafunctionofradius,thenacutof200kms−1 0 (rs2+r2)γ/2
shouldguaranteethatthepower-lawmodelofEquation(3)
Theformoftheescapespeedforthismodelcannotbewrit-
is appropriate at the locations of all the stars in our study.
ten in terms of elementary functions, so we omit it here.
Consequently, we set v =200kms−1. Note that this cut
min
isappliedtoGalactocentricradialvelocities,andsothemo-
tionofthesunmustberemovedbeforehand.Forthis,weas-
5 PRIORS AND NUMERICAL
sumealocalstandardofrestv =240kms−1 andasolar
LSR IMPLEMENTATION
peculiar motion (U ,V ,W ) = (11.1,12.24,7.25) kms−1
(cid:12) (cid:12) (cid:12)
(Scho¨nrich et al. 2010). Once this cut is applied, there re- Having constructed our likelihood function and gathered
main 1573 MSTO stars, 343 K-giants and 44 BHBs. The our data, we now need to choose explicit priors on each of
(cid:13)c 2016RAS,MNRAS000,1–??
Fast moving stars in SDSS 5
our model parameters. We split these parameters into two where µ = 15kpc and σ = 7kpc. We adopted this prior
groups: global parameters, which are the model parameters afterexperimentationwithlessinformativepriors,wherewe
thatareindependentofourchoiceforv (x),andpotential foundthatourmodel,incombinationwiththesedata,does
esc
parameters, which concern our explicit model of the escape notconstrainr (i.e.,theposteriordistributiononr always
s s
speed. resembles the prior). The power law index γ can take any
value between 0 and 1, and hence we adopt the same prior
as Equation (22).
5.1 Global priors
Weallowadifferentpowerlawslopeinthevelocitydistribu-
tion of each tracer k = (kMSTO,kK−giant,kBHB), and we fit 5.3 Sampling method
fortheoutlierfraction,f.Thesefourparametersconstitute
Armedwithourpriors,wecannowperformaBayesiananal-
our global model parameters. We anticipate that the three
ysis. Our full parameter space is Θ=(k,f,θ), which has 6
values of k should be similar, since all three sets of tracers
(SPL) or 7 (EPL and TF) dimensions. Bayes’ theorem for
belong to the halo. On the other hand, the MSTO sample
our problem can be written in the following way
spansadifferentradialrangetotheBHBandK-giantsam-
ples, and so any differences in our inferred k values could
p(Θ|data)=
suggest a radial variation in k.
S07 and P14 both used cosmological simulations in or- p(θ)p(k,f) (cid:81)3 N(cid:81)j p (vj |(cid:96)j,bj,sj,k ,f,θ)
der to place informative priors on the value of k. S07 used tot ||,i i i i j
j=1i=1
,
a flat prior over the range 2.7 < k < 4.7 and P14 adjusted p(data)
this range to 2.3<k<3.7. These priors were necessary for (26)
bothanalysesbecauseofthelimitednumberofstarsinboth
the S07 and P14 samples. Given our much larger sample of wheretheindexj referstothetracertype(MSTO,K-giant
stars, we opt for a much less informative prior on each k, or BHB). Since we do not seek to compare the evidence for
which is flat in the range 0 < k < 10. This loose prior will our three models, we do not explicitly compute the denom-
allowustocriticallyassessthesimilaritybetweencosmolog- inator of Bayes’ theorem. In order to constrain the model
icalsimulationsandtheGalaxy.Ourprioronf isflatinthe parameters, we use an MCMC approach.
permissible range 0(cid:54)f (cid:54)1. Inordertocomputethelikelihood,weevaluatetheinte-
gralinEquation(5)using200Monte-Carlosamplesfromthe
distance uncertainties of each star. We draw these samples
5.2 Potential priors once, and then use the same set for each posterior evalua-
tion,thusavoidingrandomnoiseintheposterior(McMillan
Wenowdetailthepriorsusedontheparametersoneachof
& Binney 2013). We then use the emcee code (Foreman-
thethreeformsforv (x).TheSPLmodelhastwoparam-
esc Mackeyetal.2013)toMonte-Carlosampletheposteriordis-
eters:θ=(v (R ),α).IftheGalaxyistopossessaflator
esc (cid:12) tribution.emceeisapythonimplementationoftheaffine–
declining rotation curve, as we would expect from physical
invariant ensemble sampling approach suggested by Good-
reasoning, then
man & Weare (2010), where an ensemble of N ‘walkers’
0(cid:54)α(cid:54)1. (22) is used, and the proposal distribution for a given walker is
basedonthecurrentpositionsofN/2oftheremainingwalk-
Besidesthisbasicphysicalconstraint,wedonotincludeany
ers. Here, we choose N =80. We initialise the ensemble by
more information about the value of α, so our prior is flat
randomly sampling from our prior distributions, and then
over the above range. Since v (R ) is a positive definite
esc (cid:12) evolveeachwalkerfor5000steps.Wetheninspectthetrace
scale parameter, we adopt a scale invariant Jeffreys prior
plots in each dimension in order to prune our samples for
p(v (R ))∝1/v (R ). (23) burn-in, which typically requires ∼200 steps. To assess the
esc (cid:12) esc (cid:12)
convergence of our chains, we first compute the acceptance
TheEPLmodelisidenticaltotheSPLmodelotherthanthe fraction a . For all the analyses presented in this paper,
f
inclusion of the axis ratio, q, so that θ = (vesc(R(cid:12)), β, q). 0.3 < af < 0.5 for the entire ensemble of walkers. In addi-
The priors we use for β and vesc(R(cid:12)) are the same as in tion,wecomputetheintegratedautocorrelationtimeτf for
Equations (22) and (23). For q, we follow Bowden et al. each of our chains, which is the number of posterior eval-
(2016) and use a prior uations required to produce independent samples. We find
1 τf (cid:39) 50 in each dimension, so that our chains are run for
p(q)∝ , (24) ∼100autocorrelationtimes,whichprovidesuswith∼8000
1+q2
independent samples from the posterior.
whichplacesequalweightonoblate(0(cid:54)q<1)andprolate
(q>1) axis ratios.
Our final model, the TF model, has three parameters:
θ = (v , r , γ). For the velocity normalisation v , we use a
0 s 0 6 RESULTS
Jeffreys prior, limited so that 150 < v /kms−1 < 400. For
0
the scale length, rs, we use a prior that mirrors the result Figure 2 shows the one and two dimensional projections of
found by Gibbons et al. (2014), so that the posterior samples of the six SPL model parameters, as
well as the inferred median posterior values and uncertain-
1
logp(rs)∝−2(rs−µ)2/σ2, (25) ties based on the ±34% credible intervals. All of the model
(cid:13)c 2016RAS,MNRAS000,1–??
6 Williams, Belokurov, Casey & Evans
4.09+1.01
0.68
−
3.90+1.02
0.74
− SPL
nt 7.5
a
gi
− 5
K
k
2.5 6.57+1.75
10 −1.52
B 7.5
H
B 5
k
2.5
2.04+2.03
1.22
−
18
f∗12
N
6
0 521.26+45.79
1s− 800 −30.23
km 700
/)(cid:12)600
R( 500
vesc 0.37+00..0099
−
.75
α .5
.25
2 4 6 8 2.5 5 7.5 2.5 5 7.5 10 0 6 12 18 500 600 700 800 .25 .5 .75
kMSTO kK−giant kBHB N∗f vesc(R(cid:12))/kms−1 α
Figure2.Oneandtwo–dimensionalprojectionsofourMCMCsamplesfortheSPLfit.The68%and94%credibleintervalsareshownin
the2Dprojections,andthemedianparametervaluesanduncertaintiescomputedusingthe±34%credibleintervalsofthe1Dprojections
areshownabovethe1Dhistograms.Wemultiplytheoutlierfractionf bythenumberofstarsinoursample,N∗=1960inordertomake
the inferred value easier to interpret. All of our parameters are well constrained, except kBHB, which is unsurprising given that this is
thetracersamplewiththesmallestnumberofstarsbyanorderofmagnitude.NotabledegeneraciesarebetweenkMSTO andkK−giant,
whichareconstrainedtobeclosetoequal, andbetweenvesc(R(cid:12))andkMSTO andkK−giant.Largervaluesofksuggestlargervaluesof
vesc(R(cid:12)),asexpected.Seetextfordiscussion.
parametersarewellconstrained,savefork ,whichisun- 3showstherunofv withradiusimpliedbyourinference,
BHB esc
surprising given that BHBs are by far the least numerous with associated 68% and 94% credible intervals, and the
tracer in our sample. steep drop in the escape speed is clear. For perspective, we
also show the distribution of the mother samples of each
Ourresultsimplyalocalescapespeedof521+46kms−1,
−30 tracergroupinther−v planeinFigure4.TheMilkyWay
which is in good agreement with S07 and P14. We infer a ||
loosens its grip on its inhabitants significantly: our model
powerlawindexα=0.37±0.09,suggestingthattheescape
predicts that the local escape speed is 521+46kms−1, and
speed is falling rapidly as a function of radius. The middle −30
by 50kpc this has dropped to 379+34kms−1.
panel of Figure 1 is prophetic of this, because the edge of −28
theK-giantdistributioninther−v planeissteep.Figure A priori, it is unclear what the value of k should be.
||
(cid:13)c 2016RAS,MNRAS000,1–??
Fast moving stars in SDSS 7
samplesizesinpreviousstudies,k hasneverbeenmeasured
750
from data on the Milky Way. Given our significantly larger
700 sample of stars, we are able to do this for the first time.
The two tracer samples containing the most stars, MSTO
650
and K-giants, both favour k (cid:39) 4±1, which is in comfort-
P14
600 able agreement with simulations. These results also suggest
1ms−550 that k is not a strong function of position, given the rather
/k different radial ranges probed by the MSTO and K-giant
v(r)esc500 swaemakpeler,s.aTndhefaivnofuerresnacesloignhtklyfohrigthheerBvaHluBe.saSm07plpeoiisntmsuocuht
450
that this is to be expected for small sample sizes. Nonethe-
400 less,theinferenceonk isnotinsignificanttensionwith
BHB
350 the hypothesis that k is constant. Our results vindicate the
choice of prior by S07, while the range used by P14 is a
300
10 20 30 40 50 touch on the low side.
r/kpc Figure2showsastrongdegeneracybetweenk and
MSTO
v (R ),whichcanbeencodedbytheempiricalcovariance
esc (cid:12)
matrix of the samples
Figure 3. Our inference on the escape speed as a function of (cid:20) 0.84 37kms−1 (cid:21)
Galactocentric radius. The median posterior result is shown as Cov(k , v (R ))= . (27)
MSTO esc (cid:12) 37kms−1 1713km2s−2
a dark blue line, and the 68% (94%) credible interval is a
dark (light) blue band. The result using RAVE data from P14 Thisistobeexpected,andisthereasonwhyanarrowprior
and the associated 90% credible interval is also shown, and is onk wasnecessaryinpreviouswork.Figure1ofP14nicely
in good agreement with our inference. We measure a significant demonstrates the appearance of this degeneracy for vary-
gradient in the escape speed, such that it has already fallen by
ingsamplesizes.Fortunately,oursampleislargeenoughto
∼100kms−1 byaradiusof30kpc.
locatethemaximumalongthedegeneracy.Thesamedegen-
eracy is seen between k and the local escape speed,
K−giant
though it is broader. Note that this explains why our sta-
600 tistical uncertainty on the local escape speed is larger than
thatofofP14,whofound533+54kms−1 at90%confidence,
−41
400 comparedtoour90%credibleintervalof521+88kms−1.Our
−45
larger 95th percentile of 690kms−1 is a consequence of the
200 degeneracy between k and vesc(R(cid:12)): the 95th percentile of
the posterior on k is 6, which is considerably larger
1− MSTO
ms 0 than the upper end of P14’s prior.
/k The inferred outlier fraction is very small, f (cid:39) 0.001,
v|| but non-zero. This suggests that there are one or two out-
200
− liersinoursample.InspectionofFigure1suggestsoneclear
candidate: there is an MSTO star at r (cid:39) 10kpc, shown
−400 as a black point, with a measured line of sight velocity of
518.2kms−1,whichismorethan100kms−1largerthanany
600 otherstaratacomparableradiusinoursample.Otherwise,
− 10 20 30 40 50
r/kpc therearenoobviousoutliersthroughvisualinspection.Asa
checkofthisintuition,wecalculatedtheoutlierprobability
of each star in our sample as
Figure4.Thedistributionofourmothersamplesofstarsinthe f¯p (v )
r−v plane,withhorizontaldashedlinesatv =±200kms−1, p(outlier|v ,(cid:96),b,s)= out || ,
|| || || f¯p (v )+(1−f¯)p(v |(cid:96),b,s,θ¯,k¯)
our cut in radial velocity. The coloured bands are our inference out || ||
on the escape speed as a function of radius. The ‘spur’ at nega- (28)
tiveradialvelocitiesisfromK-giantsbelongingtotheSagittarius
using the model parameters obtained by taking the median
stream. Note that the contamination in our high speed sample
valuesofeachoftheonedimensionalmarginalisedposterior
from these stars is negligible, since the maximum velocity that
the stream centroid reaches is ∼ 150kms−1 (Belokurov et al. distributions, Θ¯. The largest outlier probability is >0.999,
2014)withadispersionof∼20kms−1. andbelongstotheobjectidentifiedvisuallyinFigure1.Oth-
erwise, the largest outlier probability is < 0.01, and so we
conclude that this object is the only probable outlier in the
Leonard&Tremaine(1990)pointoutthatviolentrelaxation sample.Havingidentifiedthisoutlier,wevisuallyinspected
would lead to k = 3/2, whereas collisional relaxation gives its spectrum and image data from SDSS. From the image
k=1(Spitzer&Shapiro1972).S07furthershowedthatthe data it is clear that this object is a galaxy, and has been
PlummerandHernquistspheres(Binney&Tremaine2008) misclassified by the spectroscopic pipeline of SDSS. Having
have k = 2.5 and k = 3.5, respectively. The simulations found a galaxy contaminant in our sample, we added a fur-
analysedbyS07andP14bothsuggestk(cid:39)3.Clearly,there therconstrainttoourSQLquerythatalloftheMSTOtar-
is a relatively large range of possible values. Due to small getsshouldbemorphologicallyclassifiedasstars(aswellas
(cid:13)c 2016RAS,MNRAS000,1–??
8 Williams, Belokurov, Casey & Evans
spectroscopically,asinouroriginalquery)inordertocheck
for any other similar contaminants. The high speed sample 521.49+46.29
30.65
isreducedinsizeby11objectswhenthisstricterconstraint −
is applied: all 11 of these targets, including the galaxy out-
lier previously discussed, are shown as black points in Fig-
EPL
ure 1. We visually inspected the spectra and image data
for all of these objects, and concluded that, other than the
0.35+0.09
original outlier that we identified, they are likely partially 0.09
−
blended stars. Although this might mean that the photom- .75
etry for some of these objects may be unrepresentative to
.5
somedegree,visualinspectionoftheirspectrasuggeststhat β
theylikelyhavereliablestellarparametersandradialveloci- .25
ties.Forthisreason,weoptednottoremovethemfromour
0 1.03+0.63
sample, but we did re-run our analysis with these objects 4 −0.32
removed to check for inconsistencies. We found that none
3
of our conclusions change when these targets are removed
from the sample: the comparison of the posterior distribu- q 2
tions with and without these objects is shown in Appendix 1
B. Later, when we compare our model to the data, we re-
move the bona-fide galaxy contaminant from the sample, 500 600 700 800 0 .25 .5 .75 1 2 3 4
q
but retain the other 10 objects. vesc(R )/kms−1 β
(cid:12)
Figures 5 and 6 show the results of the EPL and TF
analyses. We removed the global parameters from the fig-
ures because the results were essentially identical to those
Figure 5.PosteriordistributionfortheEPLmodel.Theresults
foundfortheSPLmodel,withoutanyinterestingadditional
areverysimilartotheSPLmodel,andtheextraparameterq is
degeneracies.Theadditionofahaloaxisratio,q,leadstothe foundtobe∼1,suggestingthattheGalacticpotentialisprobably
conclusionthattheescapespeedfallsatthesamerateinall spherical,thoughouruncertaintiesarequitelarge.
directionsfromtheGalacticcentre,withq=1.03+0.63.This
−0.32
inturnimpliesthattheGalacticpotentialislikelyspherical,
although our uncertainties are large: the data are compati-
219.97+39.99
blewithq=0.7−1.6.Sinceq correspondstotheflattening 30.75
−
of the potential, this corresponds to a relatively wide range
offlatteningsinthedarkmatterhalo.Hopefully,withmore
data, the method will provide a more useful constraint on TF
thesymmetriesoftheGalacticpotential.Theothertwopa-
rameters of the EPL model are the same as those in the 16.20+6.63
SPL model, and the inferred values and uncertainties are 45 −6.57
indistinguishable.
Our inference on the TF model suggests an inner rota- pc 30
k
tioncurvewithamplitudev0 =220+−4301kms−1,whichstarts r/s 15
tofallasγ/2=0.22+0.07 overascaler =16+7kms−1.The
−0.06 s −7
posterior on the scale radius is no different from our prior, 0.43+0.13
0.11
and so is uninteresting. However, having assumed similar 1 −
inference on the scale radius as Gibbons et al. (2014), the .75
amplitude and power law index are entirely consistent with
their result, which we will discuss further in Section 6.1. γ .5
When translated into an escape speed profile, the inference .25
appears virtually identical to Figure 3. This suggests that 160 240 320 400 15 30 45 .25 .5 .75 1
oteurrisiantfieornesncoef tihsereelsactaivpeelyspreoebdu.st against different parame- v0/kms−1 rs/kpc γ
Figure 6. Posterior for the TF model when the prior based on
the inference of Gibbons et al. (2014) is used for rs. We infer a
rotation curve amplitude of 220kms−1 and a power law decline
ofγ(cid:39)0.4,bothofwhichareentirelyconsistentwiththeresults
ofGibbonsetal.(2014).
6.1 The mass and rotation curve of the Milky
Way
Having mapped the escape speed across the Galaxy, we are
now able to convert this measurement into a mass profile
M(r)androtationcurvev (r)fortheGalaxy.Foraspheri-
c
(cid:13)c 2016RAS,MNRAS000,1–??
Fast moving stars in SDSS 9
cally symmetric escape speed profile, we have
600
M(r) = −r2vescdvesc, (29) BooIII Herc
G dr 400
(cid:114)
dv
v (r) = −rv esc. (30)
c esc dr 200
1−
Givenourmodelofvesc(r),andourinferenceonitsparam- ms
eters, we can compute posterior distributions on M(r) and /k 0
vc(r) using these formulae. For example, the local circular √3v||
speed that we measure is 200
−
v (R )=223+40kms−1. (31)
c (cid:12) −34
400
Note that our inferences on the mass profile and rotation − Tuc2 Gru1
curverestheavilyontheassumptionthatthespeeddistribu- TriII
√ 600
tionofstarsintheGalaxytruncatesat −2Φ.Ifthestarsdo − 50 100 150 200 250
r/kpc
notfillouttothethisvalue,thenouranalysiswillunderes-
timatethedepthofthepotentialwell,whichwillleadtoun-
derestimatesinthe inferred massprofileandrotation curve
oftheGalaxy.P14showedthattheinferredhalovirialmass Figure 7. The√distribution of dwarf galaxies around the Milky
Wayinther− 3v plane.Thebluebandsareourinferenceon
increased by 20% if the escape speed was instead defined ||
(cid:112) the escape speed as a function of radius. Radial velocities have
as −2(Φ(r)−Φ(r )), with r = 3R ∼ 600kpc. √
max max vir beenmultipliedby 3,asisdoneintheliteraturetoaccountfor
Since we do not attempt to track the Galaxy’s mass out
unknown tangential velocities. If the true speeds of Triangulum
to such large radii, the possible bias incurred by effectively √
II,TucanaII,Grus1,BootesIIIandHerculesarecloseto 3v ,
||
setting rmax = ∞ will be significantly smaller than 20%, thentheyarelikelytobeunbound.
and so we henceforth assume that the velocity distribution
√
reaches −2Φ. The value of v (R ) that we obtain while
c (cid:12)
making this assumption is pleasingly aligned with a mul- jects,andsotheiruncertaintyislarge.Theirpreferredmass
titude of other methods, and provides us with confidence of M(50kpc) = 54+−236×1010M(cid:12), like the other two distri-
that systematic uncertainty caused by these considerations bution function approaches, is significantly larger than our
is unimportant relative to our statistical uncertainties. result, although the large asymmetric uncertainty removes
It is worth noting that our method clearly possesses any possible tension.
very different systematic uncertainties when compared to The final study to which we compare is that of Gib-
morecommonapproachesintheliterature.Mostdynamical bonsetal.(2014,G14),whomodelledthedisruptionofthe
models of halo tracers, like distribution function and Jeans Sagittarius stream. They exploited the fact that the apoc-
analyses,aremostsensitivetothecentralpartsoftheveloc- entric precession of the stream should be sensitive to the
itydistributions.ThisisparticularlytrueofJeansanalyses, details of the gravitational potential. Their inference pro-
whichgenerallyonlymodelthefirstandsecondmomentsof duces very similar results to our work, with M(50kpc) =
thevelocitydistributions.Distributionfunctionssatisfythe 29±5×1010M(cid:12). We can make this comparison even more
full collisionless Boltzmann equation, and therefore the en- explicit because we have also estimated the parameters of
tire infinite hierarchy of Jeans equations, but this generally themodeltheyusedintheiranalysis(TF).Whenweusethe
comesatthecostoflargesystematicuncertaintiesthatarise TFmodel,themassenclosedisM(50kpc)=33+−86×1010M(cid:12).
from the chosen form of the model (Wang et al. 2015). Our The two analyses, though very different in detail, produce
approach moves the focus to the tail of the velocity distri- near identical results.
bution,andisthereforecomplimentarytootherapproaches.
Figure8showsthemassandcircularspeedprofilesim-
6.2 The orbits of Milky Way dwarf galaxies
plied by the SPL model, along with associated 68% and
94% credible regions. Our model predicts M(50kpc) = Figure 7 shows the√distribution of known Milky Way dwarf
29.8+−65..92×1010M(cid:12). For reference, we have also plotted the galaxiesinther− 3v|| plane.Iti√stypicalintheliterature
results from a selection of other studies. Xue et al. (2008, to multiply the radial velocity by 3 as a crude way of ac-
X08), Deason et al. (2012, D12) and Williams & Evans countingforunknowntangentialvelocities.Weseethatmost
(2015, WE15) all used samples of halo BHBs taken from ofthedwarfsareenvelopedbytheescapespeedcurves,with
SDSS.D12andWE15applieddistributionfunctionmodels a similar shape to the r−v distribution of stars (Figure
||
tothedata,andinfersystematicallyhighermassesthanwe 4).However,someofthedwarfsseemlikelyunboundbased
dohere,withM(50kpc)(cid:39)45×1010M .Bothareconsistent on our estimate of the escape speed. On the other hand,
(cid:12)
with the 94% credible interval of our inference, but there is ΛCDM simulations predict that 99.9% of subhalos should
ahintthatthereisadiscrepancybetweendistributionfunc- beboundtotheirhosts(Boylan-Kolchinetal.2013).Arec-
tion methods and the present approach. X08, on the other onciliation of these two statements is to conclude that the
√
hand, compared SDSS BHBs to cosmological simulations, 3v approximation for the total speed of these dwarfs is
||
and their result is comfortably in agreement with ours. likely unrealistic in these cases. Given our inference on the
W99, like D12 and WE15, used a distribution function escape speed, the assumption that these objects are bound
approach,butappliedtheirmethodtoglobularclustersand allowsustoplaceconstraintsontheirorbits.Theredpoints
dwarf galaxies. Their sample was small, with only 27 ob- in Figure 7 are Bootes III (Grillmair 2009), Triangulum II
(cid:13)c 2016RAS,MNRAS000,1–??
10 Williams, Belokurov, Casey & Evans
(Laevens et al. 2015) and Hercules (Belokurov et al. 2007), Given a set of parameters in our posterior samples for the
all of which are likely associated with the Milky Way and SPL model, we first draw a tangential velocity from Equa-
havelargeradialvelocities.Theyarealsoatradiiwhereour tion(34).Fork,weusek ,andwenormalisethespeed
MSTO
inferenceontheescapespeedistrustworthy(althoughHer- distributionbetweenv andv (r).Then,wesolveEquation
r esc
cules is a somewhat marginal case). The two yellow points (35)forr ,r and(cid:15)andstoretheresult.Thisprocessis
peri apo
areTucana2andGrus1(Koposovetal.2015;Bechtoletal. repeated for every set of model parameters in our posterior
2015), which are probably dwarfs of the Magellanic clouds samples. The resulting histograms in r , r and (cid:15) are
peri apo
(Jethwaetal.2016).Wedonotseektoconstraintheorbital then faithful representations of p(X|v ,r).
r
propertiesofthesedwarfsbecauseoftheirmorecomplexor- The results of this procedure are shown in Figure 9.
bital histories. For the three Milky Way dwarfs, we would Allthreedwarfsareexpectedtobeonveryeccentricorbits.
liketocharacterisetheirorbitsthroughtheirpericenterradii Bootes III and Triangulum II have (cid:15)(cid:39)0.95, while Hercules
r , apocenter radii r and orbital eccentricities (cid:15). has (cid:15)(cid:39)0.8, though with a somewhat broader distribution.
peri apo
Inordertocomputeposteriordistributionsonr ,r Thisisalignedwithourintuition:iftheradialvelocityalone
peri apo
and (cid:15), we first write: is relatively close to the escape speed at the radius of the
(cid:90) dwarf,thenthetangentialvelocitycannotbelargeandhence
p(X|vr,r)= p(X|vr,r,Θ)p(Θ|data)dΘ, (32) the orbit must be eccentric. As a consequence, the dwarfs
havelargeapocenters,r ∼100−300kpcandconsiderably
apo
whereX =(rperi,rapo,(cid:15))andp(Θ|data)istheposteriordis- smallerpericentersrperi ∼10kpc(althoughtheposterioron
tributiononourmodeloftheescapespeed.Wehavewritten the pericenter of Hercules is significantly less peaked than
v|| (cid:39)vr,wherevr isthevelocityawayfromthecentreofthe for Bootes III and Triangulum II).
Galaxy, because these objects are distant. We then have Ku¨pper et al. (2016) argued that the observed elliptic-
(cid:90) ity of Hercules (e (cid:39) 0.7) and its large radial velocity are
p(X|vr,r,Θ)= p(X|vr,vT,r,Θ)p(vT |vr,r,Θ)dvT, suggestivethatithas‘exploded’asaconsequenceofitslast
pericenter passage. Using N-body simulations, they arrived
(cid:90)
= δ[X−f (v ,v ,r,Θ)] p(v |v ,r,Θ)dv , at an estimate of the orbital eccentricity (cid:15)(cid:39)0.95, which is
X r T T r T
larger than our value but not in significant tension with it.
(33)
Bootes III has a similar morphology (e (cid:39) 0.5) in keeping
wherev isthetransversevelocity.Inasphericalpotential, with the picture that these satellites are on orbits that will
T
the two velocities (v ,v ) and radius r completely specify cause them to disrupt into streams after comparatively few
T r
the orbit, and hence r , r and (cid:15). This gives rise to orbital periods. Both Bootes II and Hercules have positive
peri apo
the delta–function in Equation (33), where f represents radial velocities: they are travelling away from the Galac-
X
the relation between (v ,v ,r) and X. Finally, we need to ticcentre,suggestingthattheyhaveundergoneatleastone
r T
specify p(v |v ,r,Θ). To do this, we assume that the high pericenter passage. Triangulum II, on the other hand, has
T r
speed dwarfs follow the same distribution function as the a negative radial velocity and a relatively small ellipticity
high speed stars, given by Equation (1) (e (cid:39) 0.2). Therefore, it is plausible that Triangulum II is
on first infall and is about to undergo a large amount of
(cid:18) (cid:113) (cid:19)k
p(v |v ,r,Θ)∝ v (r)− v2 +v2 . (34) disruption on its first pericenter passage.
T r esc T r
Wemakethisassumptionforsimplicity,althoughErkaletal.
6.3 Model performance and systematics
(2016) note that the velocity distributions of subhaloes in
theVLIIsimulations(Diemandetal.2008)arenotthesame We now seek to assess how well our model fits the data by
asthoseofthedarkmatterparticles,andhavevelocitydis- performing posterior predictive checks. Since our model is
persionsofsize160−200kms−1.Measuredvelocitydisper- onlygenerativeinradialvelocities,andnotinthepositions
sionsforMilkyWayhalostarstendtobesomewhatsmaller of stars, we should compare the distribution in v of the
thanthis,at∼100−150kms−1 (Evansetal.2016).Hence, datawithourmodel.Thereisasubtletyinhowthis|m| ustbe
for a given radial velocity, larger tangential velocities could done,however.Thevelocitydistributionofhigh-speedstars
bemoreprobablethantheestimateofEquation(34),which is position dependent in the model, owing to the spatial
in turn means that the orbit of the dwarf may be more ec- variation of the escape speed. The fastest star at a given
centric than our simple calculation suggests. radius is likely to be travelling slower than its counterparts
Since fX is not analytic, we must solve for rperi, rapo atsmallerradii,wheretheescapespeedislarger.Thiseffect
and (cid:15) numerically. We do this by solving the equation means that we must take into account the number of stars
L2 that have been observed at each radius in the Galaxy for
E− −Φ(r)=0, (35) our comparison between model and data to be meaningful.
2r2
We therefore write
where Φ(r) is the gravitational potential implied by our
(cid:82)
p(v |r,Θ)p(r)dr
tmhoedteoltoafltahneguelsacrapmeosmpeenedtu,mE.iTshthiseeoqrubaittaiolnenhearsgytwaondroLotiss: p(v|||Θ) = (cid:82) p(v|||||r,Θ)p(r)drdv
rperi and rapo. (cid:15) is then given by = (cid:82) (cid:82)(vesc(r)−v||)k+1p(r)dr , (37)
r −r (k+2)−1(v (r)−v )k+2p(r)dr
(cid:15)= apo peri. (36) esc min
r +r
apo peri where p(r) is the probability of observing a star at radius
In practice, we sample p(X|v ,r) in the following way. r. We approximate p(r) by binning the mother sample of
r
(cid:13)c 2016RAS,MNRAS000,1–??