Galaxy structure from multiple tracers: II. M87 from parsec to megaparsec scales

2.2 MB·

by  L. J. Oldham

Mon.Not.R.Astron.Soc.000,000–000(2015) PrintedJanuary8,2016 (MNLATEXstylefilev2.2)
Galaxy structure from multiple tracers: II. M87 from
parsec to megaparsec scales
L. J. Oldham(cid:63), M. W. Auger
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
6
1
0 January8,2016
2

n
ABSTRACT
a
J
Following a number of conflicting studies of M87’s mass profile, we undertake a
6
dynamicalanalysisofmultipletracerpopulationstoconstrainitsmassoveralargera-

] dius range. We combine stellar kinematics in the central regions with the dynamics of
A 612globularclustersoutto200kpcandsatellitegalaxiesextendingtoscalescompara-
G ble with the virial radius. Using a spherical Jeans analysis, we are able to disentangle
the mass contributions from the dark and baryonic components and set constraints
.
h on the structure of each. Assuming isotropy, we explore four different models for the
p darkmatterhaloandfindthatacentrally-coreddarkmatterdistributionispreferred.
-
o We infer a stellar mass-to-light ratio Υ(cid:63),v =6.9 0.1 – consistent with a Salpeter-like
±
r IMF–andacoreradiusrc =67 20kpc.Wethenintroduceanisotropyandfindthat,
t while the halo remains clearly co±red, the radial stellar anisotropy has a strong impact
s
a on both Υ and the core’s radius; here we find Υ =3.50+0.32 – consistent with a
(cid:63),v (cid:63),v −0.36
[ Chabrier-likeIMF–andrc =19.00+−88..3384 kpc.Thusthepresenceofacoreatthecentre
1 of the dark halo is robust against anisotropy assumptions, while the stellar mass and
v core size are not. We are able to reconcile previously discrepant studies by showing
3 that modelling the globular cluster data alone leads to the very different inference of
2 a super-NFW cusp, thus highlighting the value of multiple-population modelling, and
3 we point to the possible role of M87’s AGN and the cluster environment in forming
1
the central dark matter core.
0
. Keywords: galaxies:ellipticalandlenticular,cD–galaxies:kinematicsanddynamics
1
– galaxies: haloes – galaxies: individual: M87 – galaxies: structure
0
6
1
:
v
i 1 INTRODUCTION the central regions via adiabatic contraction (e.g. Blumen-
X
thal et al. 1986; Gnedin et al. 2004). The current observa-
The ΛCDM paradigm of structure formation has been very
r tionalpicturereflectsthiscomplexity,withthehaloesofan
a successful in describing the Universe on large scales, but
increasing number of systems being found to favour cored
there remains some tension regarding its predictions about
or only weakly cuspy central profiles. For instance, the re-
galaxy structure. For instance, one main prediction of cold,
centlocalsurveysTHINGSandLITTLETHINGS(Hunter
collisionlessgravitationalcollapseistheformationofacen-
et al. 2007) found a large fraction of dwarf field galaxies
tral cusp in the density profile of the dark matter (DM)
to have DM density profiles that go as ρ ∝ r−0.4 within
haloes that envelope galaxies, with ρ ∼ r−γ and γ = 1
DM the central kiloparsec, and studies of low-surface-brightness
at small radii (Navarro et al. 1997; Navarro et al. 2010).
galaxies also point to relatively flat central profiles with a
However, real haloes also contain baryons, and the imprint
large scatter (e.g. de Naray & Spekkens 2011).
of baryonic physics on the DM distribution could be signif-
icant.Forinstance,feedbackfromsupernovaeandAGN,as While a great deal of progress has been made in con-
wellasdynamicalfrictionfrominfallingsatellites,couldlead strainingthehalostructureoflow-surface-brightnessgalax-
to some degree of heating and expansion (e.g. Mashchenko iesanddwarfspheroids,whereDMdominatesoverthebary-
etal.2006;Laporteetal.2012;Governatoetal.2012;Vellis- onicmass,thesituationismuchmorecomplicatedfortheir
cigetal.2014),thushollowingouttheDM,whilethecooling massive elliptical counterparts. Here, our ignorance about
and condensation of baryons could increase the density in thestellarinitialmassfunction(IMF)introducesadegener-
acybetweendarkandluminousmatter,whichmakesithard
to constrain the behaviour of the DM in the inner regions;
(cid:63) E-mail:[email protected] equally, the task of probing the gravitational potential at
(cid:13)c 2015RAS
2 L. J. Oldham and M. W. Auger
large radii is made hard by the fact that their outskirts are 2 DATA
notoriously faint. However, one way of significantly allevi-
To constrain M87’s density profile across a wide radius
atingthesedegeneraciesistousemultipledynamicaltracer
range, we use multiple dynamical tracers, combining stel-
populations, spanning a range of galactocentric radii (e.g.
lar kinematics in the central regions with GC dynamics at
Schuberth et al. 2010; Walker & Pen˜arrubia 2011; Napoli-
large radii and satellite galaxies on cluster scales. We can
tano et al. 2014). Indeed, massive elliptical galaxies are of-
thensolvetheJeansequationforeachpopulationseparately,
ten home to large populations of planetary nebulae (PNe),
provided that the underlying density distribution is known.
globular clusters (GCs) and even, in the case of brightest
We therefore take data from a number of sources, as sum-
clustergalaxies(BCGs),satellitegalaxies,andeachofthese
marised in Table 1.
populations,withitsownsignaturespatialdistributionand
kinematic profile, can be used as an independent probe of
thegravitationalpotential.ThepoolofETGsforwhichsuch 2.1 Stars
an analysis has been carried out is currently too small for
anymeaningfulconclusionstobeextracted,thoughthesug- Thedeprojectedstellarsurfacebrightnessprofilecomesinto
gestion also appears to be that the NFW profile may not our analysis at two points: first, the use of the stars as dy-
provideagooduniversalfit:studiesofBCGsbySandetal. namical tracers in the Jeans equation requires us to know
(2004), Newman et al. (2011) and Newman et al. (2013) their 3D density distribution; second, our goal is to model
have found evidence for sub-NFW density profiles in clus- M87’s mass as the sum of dark and luminous components,
ters,whileahandfulofstudiesoffieldellipticals(e.g.Grillo andthelatterissimplytheproductoftheintegrated3Dlu-
2012;Sonnenfeldetal.2015),havefoundsuper-NFWdensi- minositydensitywithsomeconstantmass-to-lightratio,Υ(cid:63),
ties.Clearly,alotremainstobeunderstoodhere,fromboth whichcanbeinferredfromthedata.Weusetheradialpro-
observational and theoretical perspectives. file for M87 presented in Kormendy et al. (2009), in which
20setsofobservationsacrossarangeofradiiandfiltersys-
ThemassiveellipticalM87,locatedatthecentreofthe tems were synthesised into a single profile in the V band.
Virgo cluster, is an ideal subject for continuing such stud- As M87 is known to have a very extended cD envelope in
ies, as it has an enormous GC population (estimated as ∼ additiontoastellarcore(eg.Chakrabarty2007),itssurface
12,000,e.g.McLaughlinetal.1994;Tamuraetal.2006;Old- brightnessprofilecannotbeaccuratelymodelledbyaS´ersic
ham & Auger 2016), making it one of the richest GC hosts profile at both small and large radii. We therefore chose to
in the local Universe. Further, the sample of this popula- model it using a more flexible Nuker profile, according to
tion for which we have high-resolution kinematic data has the following relation
been greatly expanded in recent years through the wide-
field study of Strader et al. (2011). It is also understood to
(cid:16)r (cid:17)−ζ(cid:16) (cid:104)r (cid:105)α(cid:17)ζ−η
beaslowrotatorandnearlyspherical(e.g.Cappellarietal. I(R)=I 1+ α (1)
0 r r
2006),indicatingthatitissuitableformassmodellingunder b b
theassumptionofsphericalsymmetry(thoughseealsoEm- with amplitude I , break radius r = 1.05 kpc, inner
0 b
sellemetal.2014,forrecentevidencefornon-axisymmetry). slope ζ = 0.186, outer slope η = 1.88 and break softening
However,tworecentstudiesofM87’smassdistribution,both α=1.27.Whilethishasthedisadvantageofhavingnoana-
primarily based on Strader’s GC dataset, are markedly in- lyticdeprojectionornormalisation,itismuchmoreflexible
consistent. The first of these, Agnello et al. (2014), divided than a S´ersic profile as it allows both the inner and outer
theGCsampleintothreeindependentpopulationsandused slopes greater freedom. Assuming spherical symmetry, we
a virial analysis to infer a very cuspy central density profile deproject this profile to give the 3D density shown in Fig-
(γ =1.57),whileZhuetal.(2014)combinedSAURONcen- ure 1.
tralstellarkinematicswiththeStraderGCdata(alongwith Thekinematicsoftheinner33(cid:48)(cid:48)×41(cid:48)(cid:48) ofM87havebeen
anadditionalGCsamplefromHectospec)andmodelledthe observed with the integral-field unit (IFU) SAURON, and
density profile using a logarithmic potential, thus imposing a catalogue of the first four moments of the Gauss-Hermite
a core, which they infer to be r =42±10 kpc. Their total expansion of the line-of-sight velocity distribution is avail-
s
inferred stellar masses also differ by almost a factor of two. able online1. We use the velocity dispersions, which were
AsM87isoneofthenearestandmostwell-observedBCGs, obtained from the spectra using a direct pixel fitting rou-
itseemsunsatisfactorythatitsmassdistributionshouldstill tine(e.g.vanderMarel1994).Thespectrawereadaptively
besopoorlyconstrained,andclearlythereremainsworkto binnedtoensureasignal-to-noiseofatleast60perspectral
bedone.Theaimofthispaperisthereforetouseasynthe- resolution element; uncertainties are generally less than ∼
sisofGC,satelliteandstellarkinematicdatainconjunction 20 kms−1, with a mean uncertainty of ∼9kms−1.
with flexible mass models in order to infer a density profile M87 is known to host a supermassive black hole
which is free to be cuspy or cored as the dynamics dictate. (SMBH)ofmass∼6.6×109M (Gebhardt&Thomas2009;
(cid:12)
Gebhardt et al. 2011), and this should make a significant
The paper is organised as follows: in Section 2, we in-
contribution to stellar velocity dispersions at the smallest
troduce the tracer populations used in our analysis and the
radii.TheSAURONdatasetextendsrightdowntothecen-
associated datasets. In Section 3.1 we describe our mass tre, though its resolution of 1(cid:48)(cid:48) is too low to be able to set
modelandJeansanalysis,theresultsofwhicharepresented
constraints on the SMBH mass. As we are mainly inter-
in Section 4. We discuss the implications of our findings in
ested in the behaviour of the DM and stellar components,
Section5,anduseSection6tosummariseourmainconclu-
sions. Throughout this work, we assume a distance to M87
DL =16.5 Mpc. 1 http://www.strw.leidenuniv.nl/sauron/
(cid:13)c 2015RAS,MNRAS000,000–000
M87 from parsec to megaparsec scales 3
radialcoverage datatype instruments sources
stars 1.6×10−3 -210kpc photometry multiple Kormendyetal.(2009)
0-2.5kpc kinematics SAURON/WilliamHerschelTelescope Emsellemetal.(2004)
0-0.17kpc kinematics NIFS/GeminiTelescope Gebhardtetal.(2011)
GCs 1.3-240kpc photometry MegaPrime/CFHT Oldham&Auger(2016)
2-200kpc kinematics multiple Straderetal.(2011)
satellite 35-1000kpc kinematics multiple(mostlySDSS) Kimetal.(2014)
galaxies Blakesleeetal.(2009)
Table 1.Varioussourcesofphotometricandspectroscopicdataforthedifferenttracerpopulationsusedinthedynamicalanalysis.
10−2 10−2
10−3 10−4
10−4 10−6
p(R)dR1100−−65 p(r)dr1100−−108
10−7 stars 10−12 stars
redGCs redGCs
10−8 10−14
blueGCs blueGCs
10−190−1 100 101 102 103 10−1160−1 100 101 102 103
R/kpc r/kpc
Figure 1.Left:NormalisedsurfacebrightnessprofilesforthestarlightandtheredandblueGCs,scaledarbitrarily.Thestellarsurface
brightness was obtained from a fit to the V-band profile presented in Kormendy et al. (2009), while the GC distributions come from
themodellinginOldham&Auger(2016).Right:Normalised3Dluminositydensityprofilesforthesamethreepopulations,deprojected
assumingsphericalsymmetry.
oneoptiontodealwiththiswouldbetosimplyexcludethe andcolourdistributionsandtherelativefractionsofthetwo
apertures within the central ∼ 3(cid:48)(cid:48) from our analysis. How- GCpopulations.ThesurfacedensityforeachGCcomponent
ever, though the contribution of the SMBH to the enclosed was modelled as a S´ersic profile
mass becomes sub-dominant beyond this approximate ra-
dius, it still contributes non-negligibly at larger radii and
(cid:34) (cid:35)
maybecovariantwiththestellarmass.Wethereforeinclude N(R)=N exp −k (cid:16) R (cid:17)n1 (2)
the SMBH in our mass model and constrain it using high- 0 n Re
resolutionkinematicsofthecentral2(cid:48)(cid:48),asobservedwiththe
IFU NIFS on the Gemini North Telescope (Gebhardt et al. withkn =2n−0.324,andwiththeprofilesnormalisedsuch
2011).Asexplainedinthatpaper,thesedatawereobtained that (cid:82)0RmaxN(R)dR=1, where Rmax =240 kpc is the ra-
from spectra which used laser adaptive optics corrections, diusoftheoutermostGCinthecatalogue.Theg−r,g−i
and have a resolution of 0.08(cid:48)(cid:48) and a signal-to-noise gener- colours and the GC luminosity function for each GC pop-
ally greater than 50. The velocity moments are provided in ulation were modelled as Gaussians, with radius-dependent
radius and position angle bins, though our simplifying as- colour gradients in all but the latter. The 3D deprojected
sumption of spherical symmetry allows us to combine bins profiles are also shown in Figure 1.
azimuthally. Straderetal.(2011)presentsaspectroscopiccatalogue
of737GCcandidatesaroundM87atradiifrom2kpcto200
kpc (plus one object at 800 kpc, which we exclude because
itliesoutsidetheregionoverwhichourmodelfromthepho-
2.2 Globular clusters
tometryisstrictlyvalid).Thecataloguecombinesnewmea-
WeusethecolourandradialprofilesfortheGCpopulations surements for 451 GCs – obtained using Keck/DEIMOS,
ofOldham&Auger(2016),asshowninFigure1.Thedetails Keck/LRIS and MMT/Hectospec – with literature data,
oftheinferenceareexplainedfullyinthatpaper,butinbrief, and provides a ‘classification’ of objects as blue GCs, red
archivalCFHT/MegaPrimeimagesintheugriz bandswere GCs, ultra-compact dwarfs and transient/unknown along
used to compile a sample of 17620 GC candidates, selected withSDSS-bandg−r andg−icolours.Wecross-correlate
according to their colours, magnitudes and apparent sizes, ourphotometriccataloguewiththisspectroscopiccatalogue,
andtheresultingcataloguewasmodelledtoinfertheradial selectingonlytheobjectsclassifiedinStraderetal.(2011)as
(cid:13)c 2015RAS,MNRAS000,000–000
4 L. J. Oldham and M. W. Auger
1000 16
900
800
12.8 15
600
600
14 300
400
12.6
200 0
dec/deg12.4 0200 1v/kms− dec/deg1123 −300 1v/kms−
− 600
−
12.2 400 11
−
900
600 −
−
10
12.0 800 −1200
−
1000 9
187.2 187.4 187.6 187.8 188.0 188.2 188.4 − 182 184 186 188 190 192
ra/deg ra/deg
Figure 2. Maps of the tracer populations, coloured by velocity with respect to M87. Left: kinematic GC sample. Note that this
does not follow the distribution of the underlying population, which is characterised independently using the photometric catalogue of
Oldham&Auger(2015,accepted).Right:satellitegalaxysample.M87itselfisplottedasanavystar.Notethatwedonotexpectthis
spectroscopically-selectedsampletobecomplete.
GCs,toobtainasampleof612GCswithcompleteluminos- satellite sample is most likely highly incomplete, we do not
ity,spatial,colourandkinematicinformation.Intheanalysis have access to this quantity. Most importantly, the sample
that follows, we choose to use the Oldham & Auger (2016) of satellites we use has been selected spectroscopically, and
photometry over that provided in Strader et al. (2011), for thisimposesanon-trivialselectionfunctionwhichmayalter
consistency with our GC colour distributions. the spatial distribution from the true underlying one, lead-
ing us to draw incorrect conclusions from a Jeans analysis.
Instead,wechoosetousethesesatellitestogiveanestimate
2.3 Satellite galaxies of the total mass at large radii and so give a further con-
straintinourinferenceonthemassprofile.Wedothisusing
The Extended Virgo Cluster Catalogue (EVCC, Kim et al.
thevirialmassestimatordevelopedinWatkinsetal.(2010)
2014)providestheredshiftsandpositionsontheskyof1589
andappliedtoagalaxygroupinDeasonetal.(2013),which
galaxiesinafootprintof725deg2centredonM87,extending
is designed to be robust against simple approximations to
to3.5timestheclustervirialradius.Theredshiftsarecom-
the true distributions. We use the mass estimator given by
piledfromtheSDSSDR7releaseandtheNASAExtragalac-
ticDatabase(NED),andeachobjectisclassifiedaseithera
certainclustermember,apossiblememberorabackground M(R<R )= C(cid:10)v2 rµ(cid:11), (3)
out G los
source based on morphological and spectroscopic criteria.
As it is important that our sample only contains galaxies where
moving in M87’s halo potential, we selected only those ob-
jects classed as certain members according to both criteria, µ+ν−2β
and we further cross-correlated these with the catalogue of C = I ro1u−tµ (4)
µ,ν
Blakesleeetal.(2009),whichusedsurfacebrightnessfluctu-
ationstocalculatedistancemoduli.Toavoidcontamination and
fromtheWcloud,aslightlymoredistantcomponentofthe
Virgoclusterat acharacteristic distance of23Mpc,weim- π0.5Γ(µ +1)
posed a distance cut of 20 Mpc; further, to separate the A Iµ,ν = 4Γ(µ2+ 5) [µ+3−β(µ+2)]. (5)
cloud (centred on M87) from the B cloud(centredon M49) 2 2
wealsoimposedadeclinationanglecutof9.5degreesanda Here β is the anisotropy parameter
radius cut of 1 Mpc. The spatial and velocity distributions
of the resulting sample are shown in Figure 2: it comprises σ2
60 galaxies, with radii relative to M87 ranging from 35 kpc β =1− t (6)
σ2
to 1 Mpc. r
We could use this tracer population in the same way forradialandtangentialvelocitydispersionsσ2 andσ2,and
r t
asthestarsandtheGCsandrequireitsvelocitydispersion µ and ν are the slopes of the potential and tracer density
profile to satisfy the Jeans equation so as to obtain a fur- respectively, both assumed to be scale-free. As mentioned
ther probe of the mass at scales comparable to the virial previously, our satellite sample is likely to be incomplete,
radius. However, as noted earlier, the Jeans equation re- andthismeanswecannotinferβ,µandν directlyfromthe
quiresustoknowthetracerdensitydistributionand,asour data;insteadwecalibratetheirvaluesusingsimulations.Fol-
(cid:13)c 2015RAS,MNRAS000,000–000
M87 from parsec to megaparsec scales 5
14
ρ(r)=ρ (r)+ρ (r)+ρ (r). (8)
(cid:63) DM BH
12
The stellar density profile is obtained by deprojecting the
stellarsurfacebrightnessandscalingbythestellarmass-to-
10
light ratio Υ , which we assume to be constant with radius
(cid:63)
andincludeasafreeparameterinthemodel.Toallowflex-
8
ibility in the DM density profile and to explore the impact
of changing the halo model on the mass inference, we carry
6
out our analysis using four different halo profiles. The first
is a standard NFW profile,
4
ρ (cid:16) r (cid:17)−1(cid:16) r (cid:17)−2
ρ (r)= 0 1+ , (9)
2 DM 4π r r
s s
with scale radius r and density scale ρ , and cumulative
0 s 0
0.25 0.20 0.15 0.10 0.05 0.00 0.05 0.10 0.15 0.20 mass
− − − − log−10(Mtrue/Mest)
M (r)=ρ r3(cid:16)ln(cid:0)1+ r (cid:1)− r (cid:17). (10)
Figure3.Theratioofthetruehalomasstothatestimatedfrom DM 0 s r r+r
s s
thevirialmassestimator,ascalibratedfromthesimulations.The
Following Zhu et al. (2014), we also use a cored isothermal
spreadinestimatesreflectsthefactthattheparametersβ,µand
ν doinfactvarybetweenthedifferentgroups,andarenottruly profile with a logarithmic potential (the LOG model)
global as we have assumed. However, the scatter is small at 0.1
dexandthemedianoffsetisnegligible,indicatingtherobustness ρ r2 3r2+r2
oftheestimator. ρ (r)= 0 s s . (11)
DM 4π (r2+r2)2
s
Thismodelisinherentlycored,butincontrasttotheNFW
lowingDeasonetal.(2013),weusethez=0halocatalogue has a ρ ∼ r−2 dependence at large radii: for a given core,
ofthefirstMultiDarksimulation,desribedindetailinPrada then, this profile allows more dark matter to be placed at
etal.(2012).ThisusestheWMAP5cosmologyandcontains large distances from the galaxy centre. Its cumulative mass
about 8.6 billion particles per Gpc/h3: the halo finder uses is given by
the bound density maximum technique described in Klypin
& Holtzman (1997). We identify all haloes with more than
ρ r2r3
30 subhaloes and treat each subhalo as a distinct satellite M (r)= 0 s . (12)
DM r2+r2
galaxy.Wethenusethesubhalovelocitiesandpositionsand s
theparenthalomassprofilestoinfertheposteriordistribu- We also use a generalised NFW (gNFW),
tionsonβ,µandν thatbestdescribetheglobalproperties
of the population, and use these to generate a posterior on
ρ (cid:16) r (cid:17)−γ(cid:16) r (cid:17)γ−3
the total mass M(R < Rout) of our cluster, whose median ρDM(r)= 4π0 r 1+ r , (13)
s s
and standard deviation we use to re-select haloes from the
simulationanditeratetheprocedureuntilourtheinference withfreeparametersγ,ρ0andrs,whereγistheinnerslope.
onµ,β andν converges.Wereportmedianvaluesβ =0.30, This has the advantage of leaving the profile free to choose
µ=0.15, ν =1.9, and find a cluster mass between cusps (γ (cid:62) 1) and cores (γ = 0) in the centre,
while still becoming NFW-like at large radii in a way that
is consistent with both simulations and observations. This
log(M(R<985kpc)/M(cid:12))=14.11±0.19. (7) gives a cumulative mass profile that can be related to the
Gauss hypergeometric function F as
Totestthecalibration,weusethemedianvaluesofβ,µand 2 1
ν toapplythevirialmassestimatortothesimulatedparent
haloes, and the resulting comparison is shown in Figure 3. M (r)= ρ0rs3(cid:16) r (cid:17)ω F [ω,ω;1+ω;− r ] (14)
Encouragingly, the masses are consistent, with a negligible DM ω rs 2 1 rs
median offset and a scatter of 0.1 dex, thus illustrating the where ω = 3−γ. The final profile we use is a cored gener-
robustnessoftheestimator.TheapplicabilityofthistoM87 alised NFW (cgNFW)
is then dependent on the assumption that the properties of
thesimulatedgalaxysatellitepopulationsarerepresentative
(cid:16)r+r (cid:17)−γ (cid:16) r (cid:17)γ−3
of real galaxies. ρ (r)=ρ c 1+ . (15)
DM 0 r r
s s
inwhichr isnowthescaleradiusofthecore.Thisisamore
c
general case of the gNFW, and, in addition to allowing the
3 MODELLING data to choose between cusps (r = 0, γ = 1) and cores
c
(r > 0, or r = 0 and γ = 0), it has additional flexibility
3.1 Mass models c c
at intermediate radii. We carry out the integration of the
Wemodelthegalaxydensityprofileasthesumofluminous cumulative mass for this profile numerically.
and DM components plus a SMBH: The BH is simply a point mass at the origin,
(cid:13)c 2015RAS,MNRAS000,000–000
6 L. J. Oldham and M. W. Auger
ρBH(r)= M4πBrH2 δ(r), (16) 21GI(R)σlos(R)2 =(cid:90) ∞ √lrσ2r2−rdRr2 −R2(cid:90) ∞ r√βrlσ2r2−drR2
R R
(20)
givingaconstantterminthecumulativemassdistribution.
which, for certain choices of anisotropy parameterisations,
As the SAURON dataset far outweighs the GC, NIFS
can be written in the simple form
andsatellitedatasetsintermsofsize,weregulariseitscon-
tributiontothelikelihoodcalculationbyadditionallyfitting
for a ‘noise’ parameter ∆ . This can be interpreted as ac- (cid:90) ∞ dr
counting for scatter in thσe(cid:63)2data not included in the uncer- I(R)σl2os(R)=2G lKβM r (21)
R
tainties or, alternatively, as modifying the relative weight
with the kernel K dependent on the particular anisotropy
giventothesedata,andisaddedinquadraturetothemea- β
model. Initially, we assume all orbits to be isotropic; later
sured uncertainties on the velocity dispersion.
we consider the effect of anisotropy on our inference by
Ouroverallmodelthereforehasanumberoffreeparam-
modelling the stellar population as following a radially-
eters dependent on the halo model in question and our as-
dependentanistropyprofileandeachGCpopulationwitha
sumptionsabouttheanisotropy.Intheisotropiccase,which
constant,non-zeroanisotropy.Forthestars,weuseascaled
wetreatfirst,thenumberofparametersvariesbetweenfive
Osipkov-Merritt (Osipkov 1979; Merritt 1985) profile
and seven. In common for all halo models are the mass-to-
lightratioΥ ,thenormalisationoftheDMhalolog(ρ ),the
(cid:63) 0
scaleradiusrs,theBHmasslog(MBH)andthenoiseinthe r2
β(r)=β ; (22)
SAURON data, ∆σ(cid:63)2. The gNFW and cgNFW halo models ∞r2+ra2
thenaddthefreeparameters(γ),(γ,r )respectively.When
c
which is centrally isotropic and becomes radially
the orbital anisotropy of each tracer population is allowed
anisotropic at large radii, tending to the asymptotic
to vary, this adds three parameters, as will be discussed in
anisotropy β . The kernel for the scaled Osipkov-Merritt
Section 3.2. In our notation in the following sections, we ∞
profile is presented in the Appendix; the kernel for the
use the gNFW parameter set whenever we write out the
constant-anisotropycaseisgiveninMamon&L(cid:32)okas(2005),
model parameters explicitly, but this should be understood
and we refer the reader there for further details.
as standing in for any of the halo models.
Given a prescription for the anisotropy, the surface
brightness and luminosity profiles of Figure 1 and the mass
model, then, we can calculate line-of-sight velocity disper-
3.2 Jeans modelling sionsσlos fromtheJeansequationandsoinfertheposterior
probability distribution of the model parameters, based on
Given the stellar velocity dispersion and the GC and satel-
the data in hand.
litevelocities,wewanttoinfertheposteriorprobabilitydis-
tribution on M87’s density profile. For the satellite galax-
ies,thisinvolvesadirectcomparisonofthemasscalculated
from our virial estimator and that obtained by integrating 3.3 Statistical analysis
the proposed density profile. For the stars and GCs, on the
WeusethemethodsofBayesiananalysistoinfertheposte-
otherhand,wecanrelatetheobservedvelocitiesandveloc-
rior probability distribution of our density profile parame-
itydispersionstothedensityprofileviatheJeansequation.
ters,giventhedata.Bayestheoremstatesthattheposterior
Assuming spherical symmetry and dynamical equilibrium,
distribution is proportional to the product of the likelihood
the Jeans equation has the simple form
function of the data given the model and the priors on the
model, and so the task here is to construct sensible likeli-
d β(r) GM(r) hoodfunctionsforeachdataset.Thetotallikelihoodisthen,
(lσ2)+2 lσ2 =l(r) (17)
dr r r r r2 inturn,theproductofthese,aseachconsitutesanindepen-
dent set of measurements.
where l(r) is the luminosity density of the tracer, σr(r) the First, as we have velocity dispersion measurements for
radialvelocitydispersionandβ(r)istheanisotropyparam- the stars, the likelihood of observing a particular velocity
eter defined in Equation 6. This is a first-order differential dispersion at radius R is assumed to be Gaussian, with a
equation with general solution standard deviation equal to the uncertainty. Thus the kth
stellarvelocitydispersionmeasurementgivesacontribution
1 (cid:90) ∞ GM(s) to the likelihood:
l(r)σ2(r)= f(s)l(s) ds (18)
r f(r) s2
r
(cid:32) (cid:33)2
where lnL (σ ,R |M(cid:126))=−0.5 σk2−σm2 −0.5ln(2πδ2 )
(cid:63),k k k δσ2 σk2
k
(cid:104)(cid:90) r ds(cid:105) (23)
f(r)=f(ri)exp ri 2β(s) s (19) fpoarraumnecteerrtsaiMn(cid:126)ty=δσ(kΥ2 a,nMd mo,d∆el p,ρred,γic,trion). Aσm2s,thaendobmseorvdae-l
(cid:63) BH σ(cid:63)2 0 s
(e.g.vanderMarel1994;Mamon&L(cid:32)okas2005).Projecting tions in each aperture are independent, the total log likeli-
this along the line of sight gives σ2 (R) as hood of observing the ensemble is just the sum:
los
(cid:13)c 2015RAS,MNRAS000,000–000
M87 from parsec to megaparsec scales 7
(cid:88) (cid:88) (cid:88)
lnL(cid:63) = lnL(cid:63),k. (24) lnL= lnL(cid:63),k+ lnLGC,k+lnLMsat. (27)
k k k
Note that, for measurement uncertainty ∆2 , the regu- We explore the parameter space using the ensemble-
σk2 sampling code emcee (Foreman-Mackey et al. 2013).
larisation of the SAURON data gives a total uncertainty
δ2 = ∆2 +∆2 . For the NIFS data, on the other hand,
σk2 σk2 σ(cid:63)2
δ2 =∆2 .
σ2 σ2
k k 4 RESULTS
The virial mass estimate from the satellites can be
treatedinasimilarway,thoughhereweonlyhaveonemea- Initially, we treat all orbits as isotropic and investigate the
surement: uncertaintyintroducedtothemassinferencebythechoiceof
halomodel.Theinferenceforthefourdifferenthalomodels
(cid:32) (cid:33)2 are presented in Figures 4 - 7, and are summarised in Ta-
log(M )−log(M ) (cid:16) (cid:17)
lnL =−0.5 sat mod −0.5ln 2πδM2 ble 2 along with the associated virial parameters. We find
sat δMsat sat the results of the gNFW, cgNFW and LOG models to be
(25) very similar: the gNFW and cgNFW haloes are both cen-
forvirialmassestimateMsat,modelmassMmodandthelog- trallycoredwithbreakradii∼80kpc,whiletheLOGmodel
arithm of the uncertainty from the mass estimator, δMsat. is inherently cored and has a slightly smaller scale radius
Thisisadirectcomparisonwithouttheneedforrecourseto r =48±5kpc(note,however,thatthesearenot coreradii,
s
the Jeans equation. and that the scale radius is defined differently in the latter
For the GCs we can assign probabilities, based on the case). All three models also agree closely on the high stel-
colourandpositioninformationprovidedinthephotometric larmass-to-lightratioΥ ∼6.9andthevirialparameters,
(cid:63),v
catalogue,ofeachGCbelongingtoeithertheredortheblue withlog(M /M )∼14.2andR ∼1.3Mpc.Indeed,the
vir (cid:12) vir
population, although we are not able to classify them with best cgNFW model approximately recovers the gNFW so-
certainty.Asthetwopopulationsareassumedtobedynam- lution, with an intermediate slope γ ∼3 and both the core
ically decoupled, we do not expect their velocity dispersion andbreakradiicomparabletothegNFWscaleradius.Like
profilestobethesame.Incontrasttotheothertracerpop- the gNFW, this represents a solution in which the halo is
ulations,then,wecalculatethelikelihoodofobservingaGC coredcentrallyandbecomesNFW-likeatlargeradii;unlike
withaparticularvelocityundertheassumptionthattheve- thegNFW,itallowsmoreflexibilityintheintermediatere-
locity distribution of each GC population can be described gions, and the fact that this intermediate slope is found to
byaGaussianwithastandarddeviationgivenbytheveloc- be close to the NFW value γ ∼ 3 suggests that the NFW
ity dispersion, such that profile is an adequate description at intermediate radii as
well as large radii.
(cid:32) (cid:33) On the other hand, the NFW model predicts a lower
(v −v )2 (cid:16) (cid:17)
lnLGC,k =−0.5 δkvk2+gσam2l −0.5ln 2π(δvk2+σm2 ) mthaesso-tthoe-rligmhotdraeltsioinΥt(cid:63)e,Vrm=s6o.f6i0t±s 0v.i0ri5al,tphaoruagmhetiterdso.eTshmeadtcifh-
(26) ference at small radii arises because of the NFW profile’s
wherev andδv arethemeasuredline-of-sightvelocityand hard-wired cusp, which places more DM in the central re-
k k
velocityuncertaintyofthekthGC,v =1284kms−1isthe gionsattheexpenseofstellarmass.InFigure8weplotthe
gal
heliocentric velocity of M87 (Cappellari et. al. 2011) and massprofilesinferredineachcase,alongwiththeassociated
the line-of-sight velocity dispersion at the location of each rotation curves in order to better highlight the differences
globular cluster, σ , is modelled separately for the red and between the four models. As can be seen, deviations start
m
the blue populations. toariseatlargerradii,wheretheconstraintsfromourdata
ToaccountfortheuncertaintyinassigningeachGCto areweakest.Thesearemainlyduetothedifferentstructural
either the red or the blue population, for each set of model featuresofthemodels:forinstance,theLOGprofiledecays
parameters M(cid:126) we draw 1000 Monte Carlo samples which asr−2atlargeradiiwhereastheotherthreegoasr−3,which
stochastically explore the population distribution based on allowstheformertoplacemoremassatlargeradii.Equally,
the colour, magnitude and position information for the in- theNFWprofilegoesasr−1 atsmallradiiwhereastheoth-
dividualGCs.Wethenmarginaliseoverthesamplestogive ersareeitherintrinsicallycoredorhaveflexibleinnerslopes
afinalcontributiontothelikelihood.WhilesomeGCshave (which are found to favour cores here), and this causes the
either very high or very low probabilities of belonging to NFW profile to start deviating from the others at smaller
one of the GC populations, with small uncertainty, there radii.
also exists a significant fraction with comparable probabil- Following this, we test the robustness of our inference
ities of belonging to either: for these objects, it would not on the halo structure against model assumptions by relax-
be meaningful to simply assign them to one population or ing the isotropy condition and introducing some element of
the other. A further advantage of this stochastic sampling anisotropy into each of the tracer populations. Guided by
is that it allows us to explore different combinations of red cosmological N-body simulations, (e.g. Diaferio 1999; Col´ın
and blue GCs. etal.2000;Wojtaketal.2005),whichfindorbitstochange
Aseachtracerpopulationisindependent,thefinallog- smoothly from being radially anisotropic in the galaxy out-
likelihoodofanysetofmodelparametersisthesumoverall skirts to being isotropic in the centre – a result of hierar-
contributions: chical formation in which BCGs form by the accretion of
(cid:13)c 2015RAS,MNRAS000,000–000
8 L. J. Oldham and M. W. Auger
halo
model M(cid:63)/L log(ρDM) rs γ δσ(cid:63)2 rc log(Mvir) Rvir
NFW 6.6±0.1 6.6±0.1 448±75 – 12.1±0.2 – 14.39+1.11 1620+460
−0.53 −360
LOG 6.9±0.1 7.7±0.1 48±5 – 12.0±0.2 – 14.21+0.60 1410+240
−0.37 −200
gNFW 6.9±0.1 8.3±0.1 79±10 <0.14 12±0.1 – 14.13+0.76 1320+270
−0.44 −240
6.9±0.1 5.6+1.4 273+430 2.54+0.33 12.0±0.2 63+11 14.16+0.67 1350+240
cgNFW* −1.2 −180 −1.18 −14 −0.40 −200
Table 2. Final inference on the parameters in the different β = 0 models. We report the median values of our inferred posterior
distributions,alongwiththe16thand84thpercentilesasameasureofouruncertainty.Fortheinnerslopeγ inthegNFWandcgNFW
models, the 95 % confidence value is given. All quantities are measured in units of solar mass, solar luminosity, kilometers per second
andkiloparsecs.*NotethatthecgNFWposteriorisbimodalduetodegeneraciesinherentintheprofile:seethepanelofFigure7.
infalling satellites which become phase-mixed by the time ourmodelwhichcouldaffectourinferenceontheBHmass.
they reach the galaxy centre – we choose to model the stel- For instance, our anisotropy profile enforces isotropy in the
lar anisotropy using the scaled Osipkov-Merritt profile of centre,which,whileaphysicallysoundassumption,restricts
Equation 22. As the GC data are comparatively sparse, we therangeofparameterspaceexploredbytheBHmass.We
treat each GC population as having some constant, non- arealsoonlymakinguseofthesecond-ordermomentsofthe
zeroanisotropy.Thisincreasesourparameterspacebyfour, line-of-sightvelocityprofile:whilethisshouldnotcauseany
as we are now also fitting for the scale radius of the stel- systematic bias in the BH mass inference, we may be able
lar scaled Osipkov-Merrit profile r , the asymptotic stellar to obtain tighter constraints by incorporating higher-order
a
anisotropy β , and the anisotropies of the two GC popula- moments.
∞
tons, β and β . Again, we carry out the inference for Wecomparethegoodnessoffitofthefourmodelsusing
red blue
thefourhalomodels,andpresenttheposteriorsinFigures9 a reduced chi-squared criterion, and find that the cgNFW
- 12 and summarise the results in Table 3. profile provides the best description of the dynamics – not-
ing that the cgNFW, gNFW and NFW are nested models.
Interestingly,wenowfindthesomesignificantdeviation
We therefore present the best cgNFW mass profile in Fig-
among the halo models. Firstly, the LOG profile recovers a
ure 13 and the residuals on the SAURON and NIFS data
modelforthestarsthatisquiteisotropic,withβ >0.99at
∞ in Figure 14. The mass profile appears to be generally con-
the95%confidencelevelandr =9.29±0.05kpc;giventhat
a sistentwithearlierwork;theresidualsarenoticeablybetter
theSAURONdataonlyextendoutto∼2kpc,thisimplies
than in the isotropic case, indicating the importance of the
anisotropiesβ <0.05wherewehavedata(thoughourdata
anisotropy in the stellar dynamics.
are in projection). This is presumably a result of the hard-
We also note that the assumption that the MultiDark
wiredcore.Theremainingthreemodelsagree,withinuncer-
simulationsareabletoaccuratelyreproducethedynamicsof
tainties,onthemassprofile–though,asmightbeexpected,
real galaxies may introduce some additional systematic un-
not so closely as when isotropy is enforced. In these cases,
certainty into the mass measurement made using the virial
we find that the mass-to-light ratio drops significantly and
estimatorinSection2.3.Wethereforetestthesensitivityof
thatthestarsbecomeradiallyanisotropicwithinscaleradii
ourinferencetothisuncertaintybycarryingouttheanalysis
r < 4 kpc. This covariance between the central mass and
a with twice the calculated uncertainty, and we find that this
anisotropy is a direct manifestation of the mass-anisotropy
hasanegligibleeffectonourinferreddensitymodel.Essen-
degeneracy: in the isotropic models, the assumed lack of
tially,thelikelihoodisdominatedbycontributionsfromthe
anyradialanisotropydrovethecentralmasstohighvalues,
starsandtheGCs,asthedatasetsforthesepopulationsare
whereas here we are able to go some way towards break-
far more extensive; the satellite mass estimate rejects mod-
ing the degeneracy through the use of multiple populations
elswithdiscrepantextrapolationsofthetotalmasstolarge
tracing the same gravitational potential. Note also that the
radii, but most of the models consistent with the smaller-
introductionofradialanisotropy,coupledwiththedecrease
radius data are still well-behaved at 1 Mpc.
instellarmass,nowallowstheDMhalotobecomelesscored
inthecentre:wefindagNFWinnerslopewhichisstilldis-
tinctly sub-NFW, but less extremely so, with γ < 0.81 at
the 95 % confidence level, and the cgNFW model finds a 5 DISCUSSION
smaller central core (r = 19.00+8.38 kpc) and a sub-NFW
c −8.34 5.1 Resolving previous discrepancies
intermediate slope γ =2.39+0.39 before transitioning to the
−0.43
r−3regimebeyondthescaleradiusr =412.1+123.0kpc.For Part of the impetus for this study was to resolve the dis-
s −143.0
thesethreeprofiles,bothGCpopulationsarefoundtohave crepancies between recent models of M87’s mass structure.
amildradialanisotropies.WefindaslightlyhigherBHmass The two recent analyses of Zhu et al. (2014) and Agnello
thanpreviousestimates(e.g.Gebhardtetal.2011),though et al. (2014), both of which relied heavily on the Strader et
it is consistent within uncertainties. However, we note that al. (2011) GC kinematics, disagreed on both the structure
this is the first time the BH mass has been inferred jointly of the DM halo and the total stellar mass of the system,
with the halo and stellar mass parameters using the high- with Agnello et al. (2014) finding a stellar mass-to-light ra-
spatial-resolution NIFS dataset; as we might anticipate a tio Υ ∼ 4.5 and a super-NFW cusp γ = 1.6 and Zhu
(cid:63),V
strongcovariancebetweentheseparameters–especiallybe- et al. (2014) inferring Υ = 6.0±0.3 (which converts to
(cid:63),I
tweenthestellarandBHmasses–aslightchangeisnottoo a V-band measurement of ∼ 7.5) under the assumption of
surprising.Ontheotherhand,therearestillassumptionsin a cored halo. There are a number of possible reasons for
(cid:13)c 2015RAS,MNRAS000,000–000
M87 from parsec to megaparsec scales 9
NFW model
isotropic
6.9
6.8
6.7
)0 6.6
ρ
(
g 6.5
o
l
6.4
6.3
6.2
700
600
rs 500
400
300
200
12.5
2σ? 12.0
∆
11.5
6.4 6.5 6.6 6.7 6.8 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9200 300 400 500 600 700 11.5 12.0 12.5
ϒ?,v log(ρ0) rs ∆σ?2
Figure 4.InferenceontheNFWmodelparameters,assumingisotropy.Thismodelfavoursalowermass-to-lightratiothantheothers,
duetothelargeramountofDMthatthecuspnecessarilyputsatsmallradii.AsinTables2and3,allquantitiesaremeasuredinunits
ofsolarmass,solarluminosity,kilometerspersecondandkiloparsecs.
M(cid:63)/L log(ρDM) rs γ M1B0−H9× δσ(cid:63)2 ra β∞ βred βblue rc log(Mvir) Rvir
NFW 3.05+−00..2244 7.51−+00..1132 128.4+−2274..81 – 7.04+−00..22 9.34+−00..2109 2.78+−11..2360 0.69+0.00.404 0.16+−00..2211 0.34+−00..1157 – 14.00+−00..3372 1190+−214200
LOG 6.48+−00..0022 8.25−+00..2139 22.1−+65..77 – 9.5−+00..2223 10.2+−00..22 9.29+−00..0055 >0.99 −0.17+−00..2231 0.08+−00..2239 – 14.03+−00..4249 1200+−215900
gNFW 3.33+−00..3485 8.64−+00..2345 40.3+−186.8.4 <0.81 7.2+−00..345 9.3+−00..22 1.06+−00..0170 0.67+−00..0055 0.13+−00..2210 0.23+−00..1177 – 13.74+−00..4332 990+−217400
cgNFW 3.50+−00..3326 4.76−+00..5450 412.1−+112433..00 2.39−+00..3493 7.22+−00..3440 9.33+−00..2210 1.06+−00..0085 0.67+−00..0036 0.13+−00..2210 0.23+−00..1168 19.00+−88..3384 13.87+−00..4325 1090+−212500
Table 3. Final inference on the parameters for the anisotropic models. We report the maximum-likelihhood values of our inferred
posteriordistributions,alongwiththe16thand84thpercentilesasameasureofouruncertainty.Allquantitiesaremeasuredinunitsof
solarmass,solarluminosity,kilometerspersecondandkiloparsecs.
(cid:13)c 2015RAS,MNRAS000,000–000
10 L. J. Oldham and M. W. Auger
LOG model
isotropic
7.9
7.8
)
0
(ρ 7.7
g
o
l
7.6
7.5
65
60
55
rs 50
45
40
35
12.5
2σ? 12.0
∆
11.5
6.80 6.85 6.90 6.95 7.00 7.05 7.10 7.5 7.6 7.7 7.8 7.9 35 40 45 50 55 60 65 11.5 12.0 12.5
ϒ?,v log(ρ0) rs ∆σ?2
Figure5.InferenceontheLOGmodelparameters,assumingisotropy.Allquantitiesaremeasuredinunitsofsolarmass,solarluminosity,
kilometerspersecondandkiloparsecs.
this: first, while both studies overlapped in the majority of Zhuetal.(2014)treatedalltheGCsasasingletracerpop-
their GC data, each used a different mass model, and it is ulation, but used the same SAURON data as in this study
possible that these may have unnecessarily constrained the to constrain the stellar mass-to-light ratio. Thus it is pos-
inference on the mass. Specifically, Zhu et al. (2014) chose sible that Agnello et al. (2014) lacked the data coverage in
tousealogarithmicpotentialmodelforthedarkmatter,so these central regions that would have permitted a reliable
enforcingacore,whereasAgnelloetal.(2014)usedapower distinction between cusps and cores.
law for the dark halo, thus requiring a constant slope at all
radii.Neitheroftheseallowsthehalomuchflexibilityinthe By exploring multiple mass models and data combina-
central regions. Further, while Agnello et al. (2014) relied tions, we are able to reproduce the results of both stud-
solely on the GCs, separating them into three independent ies. First, excluding the stellar kinematics and carrying out
populationsbasedontheirvelocities,positionsandcolours, the inference using just the satellites and GCs in a way
more similar to that of Agnello et al. (2014), we infer an
(cid:13)c 2015RAS,MNRAS000,000–000