Galaxy Kinematics with VIRUS-P: The Dark Matter Halo of M87

2.5 MB·English

by  Jeremy D. Murphy

Draftversion January12,2011
PreprinttypesetusingLATEXstyleemulateapjv.08/22/09
GALAXY KINEMATICS WITH VIRUS-P: THE DARK MATTER HALO OF M87
Jeremy D. Murphy, Karl Gebhardt and Joshua J. Adams
DepartmentofAstronomy,UniversityofTexas atAustin,1UniversityStationC1400, Austin,TX78712, USA
Draft version January 12, 2011
ABSTRACT
Wepresent2-DstellarkinematicsofM87outtoR=238′′ takenwiththeintegralfieldspectrograph
1 VIRUS-P. We run a large set of axisymmetric, orbit-based dynamical models and find clear evidence
1 for a massive dark matter halo. While a logarithmic parameterization for the dark matter halo is
0 preferred, we do not constrain the dark matter scale radius for an NFW profile and therefore cannot
2 ruleitout. Ourbest-fitlogarithmicmodelsreturnanencloseddarkmatterfractionof17.2+5.0%within
an aonndeegffloecbtuivlaerrcalduisutser(Rkein∼=em1a0t0i′c′)d,aritsaincgovtoer4in9g.4+−1478..528′%′ wiRthin525R4′e′.cEomxipstlientge SthAeUkRinOeNmadtaitc−ac5(o.R0ve≤rag1e3′′t)o,
10 J Rtuont=iavle4res7nec.klopOsceud(r∼mb5aessRst-eofi)ft.5lA.o7tg+−at10rh..i39itsh×mra1idc0i1ad2lyMdnias⊙mtaminccaaeklimtnhgoedMle≤ol8gs7arroeit≤nthuemronifcatdhsaetrekmllaohrsatlmomacasossms-tipvoer-ilsiggeahslta8x5ria.e3ts+−ioi22n..o54t%fh9eo.1flo+tc0ha.e2l
−0.2
] (V-band), a dark halo circular velocity of 800+75 km s−1, and a dark halo scale radius of 36+7 kpc.
O −25 −3
The stellar M/L, assuming an NFW dark halo, is well constrained to 8.20+0.05 (V-band). The stars
C −0.10
in M87 are found to be radially anisotropic out to R ∼=0.5 Re, then isotropic or slightly tangentially
h. anisotropicto ourlaststellardatapointatR= 2.4Re where the anisotropyofthe starsandglobular
p clusters arein excellentagreement. The globularclusters then become radially anisotropicin the last
- two modeling bins at R = 3.4 Re and R = 4.8 Re. As one of the most massive galaxies in the local
o universe,constraintsonboththemassdistributionofM87andanisotropyofitskinematiccomponents
r
strongly informs our theories of early-type galaxy formation and evolution in dense environments.
t
s
Subject headings: galaxies: elliptical and lenticular, cD; galaxies: individual (M87, NGC4486); galax-
a
[ ies: kinematics and dynamics

1
v 1. INTRODUCTION ing the distribution of mass on the scales of galaxy clus-
7 ters,understandingthe extentandshapeofgalaxy-sized
Dark matter is a central component of our cur-
5 dark matter halos has met with mixed success.
rent theory of large scale structure formation. Al-
9 Observationally, the study of dark matter halos in
though the nature of dark matter is unknown, sig-
1 spiral galaxies has outpaced that of ellipticals. This
nificant support for this cosmological paradigm comes
1. from well-motivated physical arguments (Gunn & Gott is largely due to the presence of extended HI discs
0 1972; Press & Schechter 1974; White & Rees 1978; found in spiral galaxies which provide a clean dy-
1 Fillmore & Goldreich 1984) and the remarkable agree- namical tracer to several effective radii (Rubin et al.
1 mentbetweenN-bodysimulationsofthegrowthofstruc- 1980; van Albada & Sancisi 1986; Jimenez et al. 2003).
: Analysis of the circular velocity curves of spiral
v ture (Frenk et al. 1985; Davis et al. 1985; Navarro et al.
galaxies provides some of the strongest evidence for
i 1995;Springel et al.2005)andobservationsofthedistri-
X the existence of dark matter on galaxy scales (see
butionofgalaxiesinthelocaluniverse(Davis et al.1982;
Sofue & Rubin 2001, for a review). Lacking the
r Colless et al. 2001).
a extended HI discs seen in spiral galaxies, progress
With the increase in computational power seen over
towards constraining the extent and distribution of
the past 30 years, the spatial resolution of numerical
dark matter in elliptical galaxies has proven a greater
simulations has improved to the point where individual
challenge. Despite this complication, evidence from
galaxiesarewellresolvedandtheirdarkmatterhaloscan
gravitational lensing (Keeton 2001; Mandelbaum et al.
be studied in detail (Moore et al. 1998a; Ghigna et al.
2006; Sand et al. 2008; Carrasco et al. 2010), X-ray gas
2000; Springel et al. 2008; Boylan-Kolchinet al. 2009).
profiles (Humphrey et al. 2006; Churazov et al. 2008;
From the study of both cosmological and galaxy scale
Das et al. 2010), planetary nebulae (PNe) and globular
simulations, different parameterizations for a universal
cluster(GC)kinematics(Cˆot´e et al.2001;Douglas et al.
dark matter density profile have emerged. Einasto in-
2007; Coccato et al. 2009) and integrated light stel-
troducedanearlyparameterization(Einasto1965,1968)
lar kinematics (Bender et al. 1994; Emsellem et al.
based on the S´ersic profile for the light distribution in
2004; Cappellari et al. 2006; Thomas et al. 2007;
galaxies(Sersic1968). Otherdarkmatter profileparam-
Weijmans et al. 2009; Forestell 2009) has shown that
eterizations have followed (Dubinski & Carlberg 1991;
elliptical galaxies are typically dark matter dominated
Navarro et al.1997;Moore et al.1998b). Whileeachpa-
rameterizationhasfoundsomelevelofsuccessatdescrib- beyond R ∼ 1.5 Re. However, not all galaxies studied
show definitive evidence for the existence of dark
Electronicaddress: [email protected] matter (Gerhard et al. 2001; Romanowsky et al. 2003;
2 Murphy, Gebhardt & Adams
Moni Bidin et al. 2010) and the best choice of dark halo sionmeasurementsextendingtoR 0.7Re andthepho-
∼
parametrization remains elusive. These open questions tometry of Young et al. (1978) to calculate the mass-to-
leave key components of our theories of the growth of light ratio (M/L) as a function of radius and estimate
structure, galaxy formation and evolution largely in the enclosed mass. Since that time, other mass estimates of
dark. M87usingX-raygas(Fabricant & Gorenstein1983;Tsai
Comparison between the results of various mass es- 1996; Matsushita et al. 2002; Das et al. 2010) and GC
timation methods return agreement for certain systems kinematics (Huchra & Brodie 1987; Mould et al. 1987;
and disagreement for others. Coccato et al. (2009) find Merritt & Tremblay 1993) have been made. A compar-
good agreement between integrated stellar light absorp- ison of these values to the mass estimate made in this
tion line kinematics and PNe data for a sample of 16 work is given in 5.2.
§
early-type galaxies. Yet in other systems the agree- The outline of the paper is as follows. In 2 we give
§
ment is poor. In an analysis of NGC1407, the cen- thedetailsofthedatasetsusedinourdynamicalmodel-
tral elliptical galaxy in a nearby evolved galaxy group, ing, with specifics on the VIRUS-P instrument given in
Romanowsky et al. (2009) find a discrepancy between 2.4. An overviewofthe data reduction steps is givenin
§
themassprofiledeterminedfromGCkinematicsandthe 3, with the complete details provided in the Appendix.
§
profiledeterminedbyX-raygas. Forthebrightestcluster 3.1 explains the extraction of the line-of-sight velocity
§
galaxyinAbell3827Carrasco et al.(2010)determinean dispersion profile and 3.3 provides details of the selec-
§
enclosedmass viastronglensing that is 10 higherthan tion of template stars and their application. In 4 we
× §
the mass determined from X-ray measurements. Mass explainthe orbit-baseddynamicalmodels. In 5 we give
§
discrepancies extend to tracers other than X-ray gas. the results of our dynamicalmodeling, with a discussion
StellarkinematicsofNGC821fromForestell & Gebhardt of our enclosed mass estimates and a comparison of the
(2010) and NGC3379 from Weijmans et al. (2009) dis- logarithmic and NFW halos found in 5.1 and 5.2. We
§ §
agreewiththe PNemeasurementsofRomanowsky et al. explore possible systematics in 5.4.
§
(2003). WeassumeadistancetoM87of17.9Mpc,correspond-
Each of these methods for estimating mass brings ing to a scale of 86.5 pc arcsec−1.
its own set of advantages, assumptions and limitations.
2. DATA
Mass estimates based on X-ray gas have the advantage
of very extended coverage, providing spatial overlap be- We make use of 3 sets of kinematic data to dynami-
tween the other methods. Yet X-ray gas analysis is lim- cally model M87. At large radii (140′′ R 540′′) we
≤ ≤
ited to massive galaxies and commonly assumes hydro- useglobularclusterkinematics(Cˆot´e et al.2001). Stellar
static equilibrium of the gas. Strong lensing mass esti- kinematics from the SAURON data set (Emsellem et al.
mates avoidthis potentialpitfall asitmakesno assump- 2004) are used within the central 13′′. New stellar kine-
tions regarding the energy distribution of the material matics, taken with VIRUS-P (Hill et al. 2008b), cover
withinthelens. However,lensingislimitedinitsflexibil- 4′′ R 238′′ and add substantially to the two-
≤ ≤
ity,astheregionsoftheuniverseavailableforexploration dimensional spatial coverage of the galaxy. We provide
aredictatedbythefixedgeometryofthelensandsource. detailsofthestellarsurfacebrightnessandglobularclus-
Velocitydispersionmeasurementsfromintegratedstellar ter data in 2.1. The SAURON stellar kinematics are
§
lightareeffectivelyavailableforalllocalsystems,butre- discussed in 2.2. In 2.3 we describe the observations
§ §
quireaparameterizationofthedarkhaloandinvolveas- made with VIRUS-P. 2.4 gives details of the VIRUS-P
§
sumptions about the degree of triaxiality of the system. spectrographand 2.5 explains the data collection.
§
There is also the challenge of getting stellar kinemat-
2.1. Photometry and Globular Cluster Kinematics
ics at large radii where the dark halo comes to domi-
nate the mass. PNe and GCs have an advantage here as The application of both the stellar surface
they typically extend to large radii, yet whether these brightness profile and globular cluster data follow
tracers follow the same dynamical history, and there- Gebhardt & Thomas (2009) (hereafter GT09). The V-
fore probe the same formation history as the stars, is band photometry comes from Kormendy et al. (2009),
not clear for all systems. A natural approach is to com- which is a combination of HST data from Lauer et al.
binevariousdatasetsandmethods inordertoapplythe (1992) and various ground-based observations. This
strengths of one method to overcome the shortcomings photometry extends from 0.02′′ to 2200′′. As the
of another. Treu & Koopmans (2004) and Bolton et al. dynamical modeling requires the stellar surface density,
(2008) take this approachto good success by using both thesurfacebrightnessprofileisdeprojectedfollowingthe
lensing and stellar kinematics to break the well known method of Magorrian (1999). Our globular cluster sur-
mass-anisotropy degeneracy (Dejonghe & Merritt 1992; face density profile comes from McLaughlin (1999) and
Gerhard 1993). isdeprojectedviaanonparametricsphericalinversionas
We focus here on the dark matter distribution in the described in Gebhardt et al. (1996). The globular clus-
giant elliptical galaxy M87, the second-rank galaxy in ter velocities are reported in Cohen & Ryzhov (1997),
the Virgo Cluster. M87 has been extensively studied Cohen (2000) and Hanes et al. (2001) and compiled
and a number of groups have made estimates of the ex- in Cˆot´e et al. (2001). We employ the same cuts to
tent of M87’s mass profile with a variety of methods. remove foreground and background contamination as
Empirical formulas, based on the virial theorem and describedin Cˆot´e et al.(2001). Thesecuts leaveus with
measurements of the central stellar velocity dispersion, 278 globular cluster velocities which we divide into 11
returned some of the earliest mass estimates for M87 modeling bins. A line-of-sight velocity dispersion profile
(Poveda 1961; Brandt & Roosen 1969; Nieto & Monnet (LOSVD) is then determined from all globular clusters
1984). Sargent et al. (1978) used stellar velocity disper- within one modeling bin as described in GT09.
The Dark Matter Halo of M87 3
2.2. SAURON Stellar Kinematics are laid out in an hexagonal array, similar to Densepak
(Barden & Wade 1988), with a one-third fill factor and
The SAURON data set provides two-dimensional spa-
tialcoverageofM87outtonearly40′′ withsuperiorspa- a large (107′′ 107′′) field of view. The large fibers
×
and field of view make VIRUS-P an extremely efficient
tial resolution to VIRUS-P. We therefore use SAURON
spectrographfor observingextended, low surface bright-
kinematicsinthecentralregionofM87. Oncethe sizeof
ness objects such as the faint outer halos of elliptical
themodelingbinsmakestheSAURONspatialresolution
irrelevant (R 8′′) the VIRUS-P data is used. We elect galaxies. Gimbal-mounted directly to the barrel of the
≥ telescope,VIRUS-Pmaintainsaconstantgravityvector.
tousebothSAURONandVIRUS-Pkinematicsbetween
8′′ R 13′′ as described in 4.2. Extensive analysis of the fiber-to-fiber wavelength
≤ ≤ § solution and fiber spatial PSF has been conducted and
The publicly available SAURON kinematics are
shows negligible evolution over a night. To quantify
parametrized by the first 4 coefficients of a Gauss-
the evolution, the location of the centers of the fibers
Hermite polynomial expansion. As our dynamical mod-
from the twilight flats taken at the start and end of the
elingfits thefullLOSVDratherthanitsmomentswere-
night are compared and found to deviate 0.1 pixels
constructtheLOSVDviaMonteCarlosimulationsbased
≤
for all nights. The wavelength solution determined
on the errors provided by SAURON. The details of this
from the arc lamps taken at dusk and dawn are also
reconstruction can be found in GT09.
compared. Typical residuals of the wavelength solution
2.3. VIRUS-P Stellar Kinematics to known arc lines show an rms scatter of 0.05 ˚A
∼
for frames taken at the same time of night. This value
The VIRUS-P data were taken during three separate
of rms scatter does not increase when arcs from both
observing runs over 10 partial nights in January 2008,
dusk and dawn are combined. The one exception occurs
February 2008 and February 2009. VIRUS-P has no with large temperature swings ( 10◦ C). Thermal
dedicated sky fibers. Therefore, sky nods are necessary ≥
contraction or expansion of the input and output ends
and constitute approximately one-third of our observ-
of the fiber bundle can lead to a change in position and
ing time. All our VIRUS-P data for M87 were acquired
stresspatternonindividualfibers. Localizedpressureon
throughacadenceof20minute scienceexposuresbrack-
a fiber can lead to focal ratio degradation (Craig et al.
eted by 5 minute sky nods. We note that while not hav-
1988; Schmoll et al. 2003) resulting in changes to the
ingdedicatedskyfiberspresentsissueswithdetermining
fiberpositionandspatialPSFoveranightandincreased
thecorrectlevelofskysubtraction,samplingtheskywith
RMSscatterinthe wavelengthsolutionresiduals. These
all 246 fibers allows us to better match the PSF varia-
effects are subtle, yet can degrade the quality of our
tionfromfiber-to-fiber while not adding substantiallyto flat-fielding. Therefore,ifatemperaturechange 10◦ C
our photon noise. A discussion of both the advantages ≥
is seen over a night, the data is split into two groups
and drawbacks of sky nods, and the details of our sky
and reduced using the calibration frames taken at
subtraction method are given in A.2.
the closest temperature. We found this approach was
The VIRUS-P data forM87 consistsof5 pointings ex-
necessary for two nights in our January 2009 observing
tending to 238.0′′ (20.6 kpc). The pointing placements
run. However, even when a steep temperature gradient
are shown in Figure 4. Exposure times and radial dis-
is seen, wavelength and flat-field calibration frames
tances for each pointing are givenin Table 1. Ten of the
are necessary only at the start and end of a night’s
51 science exposures were taken under marginal condi-
observing.
tions and withheld from the final data set as they de-
The medianspectralresolutionforthis VIRUS-Pdata
graded our signal-to-noise (S/N). The exposure times
is 4.75 ˚A FWHM as determined from Gaussian fits to
quoted in Table 1 include only the data that went into
strong emission lines in the arc lamp frames. This reso-
the final spectra and subsequent modeling.
lutioncorrespondstoaninstrumentaldispersion(sigma)
2.4. The VIRUS-P Instrument of 150kms−1 at4060˚Aand 112kms−1 at5400˚A.
∼ ∼
VIRUS-P was refocused between our January/February
The Visible Integral-field Replicable Unit
2008 and February 2009 observing runs which led to a
Spectrograph- Prototype (VIRUS-P), currently de-
non-trivialchange (∆FWHM 0.5˚A) in the instrumen-
ployed on the Harlan J. Smith 2.7 m telescope at ≃
talresolution. Aswefrequentlycombinethespectrafrom
McDonald Observatory (Hill et al. 2008b), is a pro-
differentfibersanddifferentnights,thechangeininstru-
totype for the VIRUS instrument (Hill et al. 2006).
mental resolution is taken into account when extracting
VIRUS is a replicated, fiber fed spectrograph currently
thestellarLOSVDs. Thedetailsofhowdifferencesinin-
under development for the Hobby Eberly Telescope
strumental resolution are handled can be found in 3.3.
Dark Energy eXperiment (HETDEX) (Hill et al. §
The assumption of a Gaussian spectral PSF for
2008a). Originally designed to conduct a Lyman-alpha
VIRUS-P proves to be a good one. To quantify this,
emitter survey (Adams et al. 2010; Blanc et al. 2010a),
we fit Gauss-Hermite coefficients to 4 bright lines in our
the VIRUS-P spectrograph is proving an excellent
mercury-cadmium arc lamp frames for all 246 fibers.
stand-alone instrument for a wide range of scien-
Over the 4 spectral lines and all fibers the median H3
tific problems (Adams et al. 2009; Blanc et al. 2009;
coefficient is 0.003 0.013. The median H4 coefficient
Yoachim et al. 2009; Blanc et al. 2010b; Yoachim et al. ±
is 0.0003 0.0117. Any non-Gaussian line behavior is
2010). VIRUS-P is a gimbal-mounted integral field unit ±
further mitigated by the high dispersion of M87, which
spectrographcomposed of 246 optical fibers each with a
4.1′′ on-sky diameter. The CCD is a 2048 2048 back- puts us well above the instrumental resolution.
×
illuminated Fairchild 3041 detector. The wavelength
range for these observations is 3545–5845 ˚A. The fibers 2.5. Data Collection
4 Murphy, Gebhardt & Adams
Calibrationframes, taken at the start and end of each velocitybinsandthenconvolvedwithasetof12template
observing night, consist of a set of twilight frames, mer- starstakenfromtheIndo-UStemplatelibrary. Selection
curyandcadmiumarclampframesandbiasframes. The ofthetemplatestarsisdiscussedin 3.3. Thecontinuum
§
twilight frames are used for both flat-fielding and deter- is divided out of both galaxyand template spectra prior
mining the position and shape of each fiber profile. The to fitting. The fitting routine works by allowing both
arc lamp frames are used for the wavelength solution the weights givento each of the 29 velocity bins and the
and determination of the instrumental resolution. (see weightsgiventoeachtemplatestartovary. Aparameter
A.1 for more details). The remainder of an observing to allow for an adjustment to the overall continuum of
night involves a sequence of 5 minute sky nods and 20 the template stars is also allowed to float. Minimization
minute science frames. The sky nods were taken 30′ off of the residuals of the fit of the convolved stellar tem-
the galaxy center in a region of sky with minimal field plate spectra to the galaxy spectra is used to determine
starsandcontinuum sourcesandwhere the galaxyhas a the best LOSVD for that given spatial bin and spectral
surfacebrightnessofµb 26.5(Kormendy et al.2009). region.
∼
Whilethispositionstillincludesintraclusterlightknown One of the great advantages VIRUS-P provides in the
to extend across much of the core of the Virgo Cluster extractionoftheLOSVDandsubsequenterrorestimates
(Mihos et al. 2005), the contribution to the total flux is isitswidewavelengthrange( 2200˚A).Thewidewave-
very low. length coverageallows us to d∼etermine the best LOSVD
in five different wavelength regions. Of the 5 spectral
3. DATA REDUCTION OVERVIEW regions sampled (Table 2), 4 of the spectral regions are
used in the final modeling. The Ca H + K spectral re-
We providea briefoverviewofthe data reductionpro-
cess here, up through extraction of the kinematics. The gion (3650–4150 ˚A) proves difficult to fit and exhibits a
extensive details can be found in the Appendix. large systematic offset in all of the first 4 moments of
The primary data reduction steps are completed with the LOSVD from the other 4 spectralregions,likely due
Vaccine, an in-house data reduction pipeline developed to issues with the continuum division. This region is
for VIRUS-P data. The reduction steps are as follows. therefore not included in the determination of the final
All of the science, sky and calibration frames are over- LOSVD and error estimate.1 The final LOSVD is cre-
scan subtracted. A master bias is created by combining atedbytakingthe averageofthe4LOSVDs withineach
all the overscan-subtracted bias frames taken during an ofthe29velocitybins. Figure6showstwoofthefinal88
observingrun. The arcs and twilight flats are then com- LOSVDs,witherrors,forabinatR=24′′andR=174′′.
bined using the biweight estimator (Beers et al. 1990). Overplottedin these figures are the LOSVDs from the 4
A 4th order polynomial is fit to the peaks of each of spectral regions used to generate the final LOSVDs.
A smaller systematic offset was observed for the Mg b
the 246 fibers for each night. We refer to this as the
spectral region (see Figures 8 and 9). Yet unlike the Ca
fiber trace. This polynomial fit is then used on each
H + K offset, which stems from the difficulty in deter-
science and sky frame to extract the spectra, fiber by
mining the placementof blue continuum, we believe this
fiber, within a 5 pixel wide aperture centered around
offset is inherent to the Mg b spectral region and there-
the trace of the fiber. The wavelength solution is de-
fore elect to include it in our final LOSVDs and subse-
termined for each fiber, and for each night, based on a
4th order polynomial fit to the centers of known mer- quent modeling. This decision was made as a trade-off
betweenthe 10%offsetinvelocitydispersionseenwith
cury and cadmium arc lamp lines. The twilight flats are
thisspectral∼region,andthemitigatingeffectsa4thspec-
normalized to remove the solar spectra. These normal-
tral region has on the statistics of the final LOSVD and
ized flats are then used to flatten the science and sky
uncertainty estimates. We note also that by including
data. Once the frames are flattened, the neighboring
the Mg b spectral region, our claim of a massive dark
sky frames are appropriately scaled, combined, and sub-
matterhaloisstrengthenedasthedirectionfortheMgb
tractedfromthescienceframes. Cosmicraysarelocated
offset is towards lower velocity dispersions. We pick up
and masked from each 20 minute science frame. For the
this discussion in 3.5.
dynamicalmodeling,thegalaxyisdividedintoaseriesof
§
line-of-sight radial and angular spatial bins. Therefore,
fibers that fall within a spatial bin are combined. This 3.2. Uncertainty Estimates
stepleavesuswithindividualspectrafor88differentspa-
Errorestimatesforthebest-fitLOSVDforeachspatial
tial bins. Of these 88 spectra, the 8 central spectra are
bin are determined in two ways. The first is made via
withheldfromthedynamicalmodeling,asthe SAURON
MonteCarlosimulationswhilethesecondisanempirical
data have superior spatial resolution in the central re-
method that makes use of the wide wavelengthcoverage
gion. The next step before the data is ready to model is
ofVIRUS-P.Then,foreachvelocitybinineachLOSVD,
the determination of the line-of-sight velocity dispersion
thelargestoftheuncertaintiesistakenastheuncertainty
profile, described below.
for that velocity bin. Both methods are described here.
ThefirsterrorestimateismadebyaMonteCarloboot-
3.1. Extraction of the LOSVD
strap method for each of the 4 spectral regions used in
Our method for determination of the line-of-
sight velocity dispersion profile (LOSVD) follows 1 Since the completion of the dynamical modeling, the contin-
Gebhardt et al. (2000) and Pinkney et al. (2003). We uumnormalizationissueexperiencedwiththeCaH+Kregionhas
give an overview of the method here. beensolved. However,thisregionisnotincludedinthedynamical
models as the cost of re-running 1000’s of models is prohibitive.
Tobegin,aninitialguessforanon-parametricLOSVD
We note Figure 8 where the Ca H + K region is included in the
of the stars is made. This LOSVD is distributed into 29 analysisofthesystematicoffsetseenintheMgbspectralregion.
The Dark Matter Halo of M87 5
thefinalLOSVD.Thebest-fitconvolvedLOSVDandset As focal ratio degradation is the dominant characteris-
of weightedtemplate stars providethe starting point for tic of an opticalfiber impacting instrumental resolution,
100MonteCarlorealizations. Eachrealizationinvolvesa differences inthe instrumentalresolutionacrossthe spa-
randomly chosen flux value, drawn from a Gaussiandis- tial direction of the chip are due to optical effects after
tribution, for each wavelength. The mean of the Gaus- the light exits the fiber. As resolution changes stem-
sian distribution is the flux from the best-fit convolved ming fromopticaleffects shouldbe continuous,aboxcar
template spectra, and the standard deviation is set as smoothing of the instrumental resolution values is justi-
the mean of the pixel noise values for that spatial bin as fied. Differences in the calculated instrumental resolu-
determined in the Vaccine reductions. A new LOSVD tion fromnight to night overan observing runare 1%
∼
is determined for all 100 realizations and provides a dis- and so one instrumental resolution map is made for an
tribution of values for all 29 velocity bins in the best-fit entire observing run. The worst instrumental resolution
LOSVD. The errorestimate is the standard deviationof over our data set is 5.0 ˚A FWHM at 4060 ˚A and 4.4 ˚A
the 100 realizations within each of the 29 velocity bins. FWHM at 5673 ˚A. Once an estimate of the instrumen-
This Monte Carlo simulation is run on all 4 spectral re- tal resolution for every fiber and for each observing run
gions and returns 4 error estimates for each of the 29 is made, the instrumental resolutionfor each fiber going
velocity bins in each of the 88 spatial bins. into a spatial modeling bin receives a normalized weight
The second method for estimating the uncertainty is based on the number of exposures going into the final
made by calculating the standard deviation of the 4 spectra. This approach gives more weight to fibers that
LOSVDs within each of the 29 velocity bins. This error providemoreweighttothefinalspectrawhileaccounting
estimate,combinedwiththe4fromtheMonteCarlosim- for differences in instrumental resolution between fibers
ulations, gives us 5 estimates of the uncertainty within and observing runs. Due to M87’s high velocity disper-
each of the 29 velocity bins of the LOSVD. The largest sion,thisstepamountstoanegligiblechangeinthefinal
uncertainty at each of these steps is taken as the final LOSVD.
uncertainty used in the dynamical modeling. We note Initially, we explored using template stars taken with
thatboththe MonteCarloandempiricalmethodforde- VIRUS-P to avoid the complications of convolving the
terminingtheuncertaintyreturnsimilarresults,withthe template spectra with the instrumental resolution. The
empirical method typically being larger. resultsachievedbythismethodprovedlessrobustfortwo
primary reasons. First, the S/N of the Indo-US spectra
3.3. Stellar Template Library
is very high. While it is certainly possible to reach this
The template stars used in the extraction of the S/N with VIRUS-P, there are observing time costs to
LOSVD come from the Indo-US spectral library consider. As using template stars taken with the instru-
(Valdes et al. 2004). The 12 stars in our final template ment is effective only if we are able to fully sample the
library (Table 3) were chosen from an initial list of 40 instrumentalresolutionacrosstheCCD,manyexposures
stars selected to cover a range in stellar type and metal- on the same template star are necessary. The second
licity. These12starswereselectedfromtheinitiallistas limitation is the variety of stellar types available during
they returned the lowest residuals when fit to the spec- an observing run. Although some variety in stellar type
tra while still maintaining a good range in stellar type. and metallicity is achievable, significant observing time
As the resolution of the template stars does not match would be lost in attempting to build up a sufficiently
the instrumental resolution of VIRUS-P, we must con- diverse stellar library.
volve the template stars with the instrumental resolu-
3.4. Moments of the LOSVDs
tionofVIRUS-P.Theinstrumentalresolutionvariesboth
between fibers and, as the instrument was refocused in InFigure7weplotthefirst4Gauss-Hermitemoments
April2009,betweenobservingruns. Afurthercomplica- of the LOSVDs from each of our 88 spatial bins. The
tionisthatspectrafromseveralfibersareoftencombined colored diamonds indicate the angular position on the
to reach the desired S/N. For overlapping pointings this galaxy,withblackalongthemajoraxisfollowedbyblue,
combination can involve spectra from opposite ends of green, orange and red falling along the minor axis. For
the CCD where the instrumental resolution can be dif- visual clarity, error bars are plotted only for data along
ferentbyasmuchas0.7˚AFWHM.ForagalaxylikeM87, the major axis. The error bars along the other axes
with velocity dispersions around 300 km s−1, the error are of comparable size. The vertical dashed lines indi-
introduced by ignoring this difference is small ( 2%). cate where the SAURON kinematics are used over the
∼
Asimple solution,particularlygivenM87’shighvelocity VIRUS-P data in the dynamical modeling. Overplot-
dispersion, would be to convolve all the spectra to the ted with open diamonds are moments from the best-fit
lowest instrumental resolution. However, as we are in- logarithmic model at each spatial bin, after averaging
terested in developing data reductionmethods to accept overthe angularbins. To minimize visualconfusion, the
allof the galaxiesin our sample, we avoiddegradingour model fits have not been plotted in the central region.
resolution to the lowest value in the following manner.
3.5. Systematics in Stellar Kinematics
The instrumental resolution is calculated from Gaus-
sian fits to 8 unblended arc lines from the arc lamp cal- We have found a systematic offset between our mea-
ibration frames taken each night. As the instrumental surement of velocity dispersion when compared to the
resolution values are noisy, particularly at weaker spec- SAURON data set. The offset is localized around the
tral lines, a small, smoothing boxcar (5 fibers wide) is Mg b lines. Figure 9 plots the velocity dispersion mea-
run along the spatial direction. Measurements of the fo- sured for the combined VIRUS-P wavelength regions
cal ratio degradation of the VIRUS-P fibers show mini- usedin the dynamicalmodeling (redcircles)and the ve-
malfiber-to-fibervariation( 2%)(Murphy et al.2008). locity dispersion calculated from just the Mg b region
≤
6 Murphy, Gebhardt & Adams
(green diamonds). Also plotted are the SAURON veloc- marily by fits to the Mg b lines, but rather springs from
ity dispersions for M87 (black squares). The SAURON issue in fitting that spectral region as a whole. A sys-
spectral range is 4810–5310 ˚A and shows a similar off- tematic study of various kinematic fitting methods over
set to the VIRUS-P Mg b spectral region. To highlight different spectral regions would be highly illuminating.
this difference we have plotted, in Figure 8, the VIRUS-
P spectrafor 5spectralregions,alongwiththe template 4. DYNAMICAL MODELS
fits(red)andcalculatedvelocitydispersionforeach. For
this particular spatial bin at R = 24.1′′ the velocity dis- We employ axisymmetric orbit-based dynamical mod-
eling based on the idea first presented in Schwarzschild
persion determined from the Mg b region is lower than
the mean of the other 4 regions by 30 km s−1. This (1979). The specific details of our axisymmetric mod-
∼ eling can be found in Gebhardt et al. (2000, 2003),
offset is not an outlier as can be seen in Figure 9. To
Thomas et al. (2004, 2005) and Siopis et al. (2009).
place a number on this offset we note that the average
The models have been shown accurate to 15%
velocity dispersion of all the VIRUS-P data points be-
∼
tween 7′′ R 36′′ is 301.8 km s−1 when all 4 spectral for recovery of the dark matter halo parameters
≤ ≤ (Thomas et al. 2005) and stellar M/L (Siopis et al.
regions used in the dynamical modeling are included as
2009). Several other groups have developed their
described in 3.1. The average velocity dispersion when
§ own modeling based on Schwarzschild’s orbit-based
using just the VIRUS-P Mg b region over the same spa-
tial range drops to 281.8 km s−1. Over the same spatial method. Dressler & Richstone (1988) and Rix et al.
region(7′′ to 36′′) the averageSAURONvelocity disper- (1997) developed an orbit-based dynamical model-
sion is 287.0 km s−1. ing code for spherical systems. van der Marel et al.
(1998),Cretton et al.(1999),Gebhardt et al.(2000)and
The cause for this offset is unknown and we do not
Verolme & de Zeeuw(2002)generalizedtoaxisymmetric
attempt a detailed analysis of the offset here. Consid-
systems and van den Bosch et al. (2008) has developed
ering the good agreement between the SAURON and
a triaxial code. Now a number of groups have employed
VIRUS-P results for the Mg b spectral range, and the
Schwarzschild’s orbit-based method for black hole mass
different methods used by both data reduction pipelines
determination(Cretton et al.1999;Verolme & de Zeeuw
to extract stellar kinematics, the offset is likely intrin-
2002;Cappellari et al.2002;Gebhardt & Thomas2009),
sic to this spectral region. The issues surrounding
stellar orbital structure and dark matter content
the Mg b spectral region for determination of the ve-
(Cretton et al. 2000; Gebhardt et al. 2003; Copin et al.
locity dispersion of elliptical galaxies, and the correla-
2004; Krajnovi´c et al. 2005; Cappellari et al. 2006;
tions with both galaxy luminosity and velocity disper-
Thomas et al. 2007; Forestell & Gebhardt 2010).
sionarewellknown(Terlevich et al.1981;Dressler et al.
We give an outline of the modeling procedure here.
1987; Worthey et al. 1992; Kuntschner et al. 2001).
First, the galaxy’s surface brightness is deprojected into
Barth et al. (2002) compare the velocity dispersion val-
a three-dimensional luminosity density. An edge-on in-
ues measured from the Mg b and Ca triplet spectral re-
clination is assumed and so the deprojection is unique.
gions for a sample of 33 local galaxies. They find that
Next, a trial gravitationalpotential is determined based
theMgbregionismoresensitivetochangesinthefitting
on the three-dimensional light distribution and an ini-
procedurethantheCatripletregion,andexhibitsanoff-
tial guess for the stellar M/L, central black hole mass,
set in velocity dispersion for 48% of the galaxies in their
and the dark matter halo parameters. Our orbit library
sample,yetwithroughlyequalnumbersofgalaxiesshow-
is the same as used in GT09. The galaxy models ex-
ing higher velocity dispersion values from either one or
tend to 2000′′ over 28 radial, 5 angular, and 15 velocity
the other spectral region. Barth et al. also compare the
bins. The gravitational potential and force are calcu-
velocitydispersionsoftheirMgbregioncalculatedwhen
lated on a grid that is 5 times finer than the grid used
both including and excluding the 5150–5210 ˚A spectral
to compare to the data. On average, 25000 orbits are
window. For 32 of their 33 galaxies they find lower ve-
run in the trial gravitational potential. A superposition
locity dispersion values when this region is suppressed
of these orbits is created that is both constrainedby the
from the fitting, with a clear trend towards a larger off-
lightdensityprofileandisabestmatchtothekinematic
set with higher galaxy velocity dispersion. We have ex-
data. The superposition is accomplished by giving each
ploredthistrendbysuppressingasimilarspectralregion
orbit a weight as determined by maximizing the func-
(5150–5220 ˚A) from our fitting and find similar results;
tion Sˆ = S - αχ2. Here, S is an approximation to the
velocity dispersion values calculated from the VIRUS-P
Boltzmann entropy, χ2 is the sum of squared residuals
spectra where the Mg b lines are withheld from the fit-
between the model and data LOSVDs (Eqs. 5 and 6),
ting are systematically lower than when these lines are
andαisasmoothingparameter. See Siopis et al.(2009)
included. However, the magnitude of our offset is small
( 3 km s−1) and is 10% of the offset seen by Barth for a detailed description of both the creation of the or-
∼ ∼ bit library and determination of the orbit weights. The
et al. This discrepancy in the magnitude of the absolute
stepsabovearethenrepeatedforadifferentmodel, each
offset value is likely due to differences in the two kine-
with a different stellar M/L, dark halo circular velocity
matic extraction routines used. Interestingly, it is in the
and scale radius.
opposite directionasnaivelyexpectedfromacomparison
Three types of models are run. First, we ran a set
of the SAURON and VIRUS-P Mg b regions as exclud-
of dynamical models with no dark matter halo. As the
ingtheMgblinesleadstolowervelocitydispersions,not
only free parameter is the stellar M/L, only 100 models
higher ones. This suggests that the driving force behind
areneeded to fully explore the parameterspace. For the
the overall offset between the Mg b spectral region and
cored logarithmic halo (Eq. 3) we ran 6500 dynamical
the other 4 VIRUS-P spectral regions is not driven pri-
models. Over 8500 models are run assuming an NFW
The Dark Matter Halo of M87 7
dark halo profile (Eq 4). For each model a distinct set the central model bins (R . 8′′) we elect to use just
of orbital weights is used and takes approximately 1.5 SAURONdataforitssuperiorspatialcoverage. Between
hours of cpu to run. We use the Lonestar computer at 8′′ R 16′′ we use both SAURON and VIRUS-P
≤ ≤
the Texas Advanced Computing Center (TACC) at The data. We do not combine these data, but rather send in
University of Texas, Austin to complete all our dynami- twoLOSVDsindependentlyintothedynamicalmodeling
cal modeling. routines. AtotalofNstars =25+80LOSVDs(SAURON
L
+ VIRUS-P) are used in the modeling, with each stellar
4.1. Model Assumptions LOSVD, stars, sampled by N =15 velocity bins.
vel
We calculate three types of dynamical models, each To deteLrmine the best-fit model, a χ2 minimization is
assumingadifferentmassdistribution. First,weconsider run in each trial potential. The χ2 is calculated as
a mass model for M87 with no dark matter halo. The
mass distribution (ρ) for these models takes the form NsLtarsNvel stars model(ν) 2
χ2 = Lij −Lij (5)
stars ∆ stars
ρ(r)=Υν(r)+M•δ(r) (1) i=1 j=1 Lij !
X X
where Υ is the stellar M/L, ν is the three-dimensional Here, model(ν) is the ith model LOSVD in the jth ve-
lightdensityandM• istheblackholemass. Astheblack locity Lbiijn. The orbit model is forced to reproduce ν,
hole is better constrained from GT09 we set our black
the stellar density, to machine precision. The residuals
hole mass to 6.4×109 M⊙ for allourdynamicalmodels. between the model and actual set of 105 LOSVDs are
Gebhardt et al. (2010, submitted) has refined the black
minimized for a single model potential, yielding a single
thhoilsesmmaaslsl ecshtaimngaeteisofwMith8i7ntoou6r.6u(n±ce0r.4ta)i×nti1e0s9aMnd⊙,doyeest χ2stars value.
As the globular clusters can have a different orbital
not change our results.
structure than the stars, they are treated as a separate
Both the second and third sets of dynamical models
kinematiccomponent. The GCs arehandledina similar
include a parameterization for a dark matter halo. The
fashion as the stars, with the difference that we employ
mass distribution then becomes a sum over each of the
a deprojected number density for the GCs rather than
mass terms as follows
thestellarluminositydensityasforstars. Boththestars
and GCs are then treated as massless test particles that
ρ(r)=Υν(r)+M•δ(r)+ρDM(r) (2) orbitinapotentialestablishedbytheassumedBHmass,
where the first two terms are the same as in equation stellar M/L and dark halo parameters. The weighted
1, and ρ (r) is the dark matter density term. Two orbit superposition in each trial potential is determined
DM
different parameterizations for the dark matter halo are byminimizingasimilarequationasforthestars,namely,
explored. The first is a logarithmic dark matter halo
with a density profile as given by NGLCNvel GC model(ν) 2
χ2 = Lij −Lij (6)
ρDM(r)∝vc2(2rc2rc2++r2r)22 (3) GC Xi=1 Xj=1 ∆LGijC !
where GC are the NGC = 11 globular cluster LOSVDs
L
The logarithmic halo features a flat central density core built uLp from individual GC velocities as described in
of size rc and an asymptotically constant circular veloc- 3.1. As with the stellar density, ν, the GC number
ity,vc. Theseconddarkmatterdensityparameterization §density is reproduced to machine precision.
is a Navarro-Frenk-White (NFW) profile (Navarro et al.
1996) as given by 4.3. χ2 Analysis
1 A χ2 analysis is used to determine both the best-fit
ρDM(r,rs)∝ (r/rs)(1+r/rs)2 (4) modeling parameters and their uncertainties. We can
rule out a model with no dark matter with high con-
The NFW halo diverges like r−1 towards the center and fidence. The best-fit no dark matter model returns a
drops as r−3 with radius. The concentration (c), scale stellar M/L = 11.4. However,the χ2 minimum for this
v
radius (rs) and the virial radius (rv) are related via c = model is 4898, which is a ∆χ2 3571 increase over ei-
rv/rs. Both dark matter parameterizations are included ther of the best-fit models inclu≥ding dark matter. We
in the modeling as described in Thomas et al. (2005). do not discuss these models further. The best-fit mod-
els for both the logarithmic and NFW halos returns χ2
4.2. Modeling the Stars and Globular Clusters minima of 1299.4 and 1310.1 respectively. The ∆χ2 of
The spatial grids used for the modeling are the same 10.7 between the two dark matter parameterizations is
asinGT09. Thespatialbinning is splitintoNr =28ra- statisticallysignificantwhencomparingthedifferentcon-
dialand Nθ =5 angularbins. Where the modelbins are straints we get on the stellar M/L. However, we do not
larger than the SAURON bins, we re-bin the SAURON getaconstraintoneitheroftheNFWdarkhaloparame-
data. The re-binning is accomplished by taking the ters, concentration and scale radius. This is clearly seen
average of all the LOSVDs falling within one model in the lower, right panel of Figure 10 where no clear χ2
bin, weighted by their uncertainties. This complication minimum for scale radius is seen out to 350 kpc. As our
doesn’t arise with the VIRUS-P stellar data as we sim- kinematicdatadoesnotextendbeyond50kpcweshould
plycombineallthespectrafromallfibersthatfallwithin not expect to get a constraint much beyond this radial
a given model bin before extraction of the LOSVD. For distance. As we do not constrain the NFW dark halo,
8 Murphy, Gebhardt & Adams
we focus here onthe logarithmichalo resultsfor ourdis- Further evidence for equilibrium between the large radii
cussion of the χ2 analysis, and refer to the NFW results starsandGCs is seeninthe excellentagreementin their
where appropriate. velocity anisotropy (see 5.3 and Figure 15).
§
To select the best-fit dynamical model we analyze the ThedegreetowhichtheGCshelptoconstrainthedark
χ2 values returned from each model run. The χ2 values matter profile can be seen in another light in Figure 13
plotted in Figure 10 are the additive combination of the where we plot the enclosed dark matter fraction for the
χ2 values of both stars and GCs, namely, χ2 = χ2 + logarithmic dark matter halo. The solid red line shows
stars
χ2 . Figure 10 plots these χ2 values against the three the dark matter fraction when including both GC and
GC
model parameters for both the logarithmic and NFW stellar kinematics in the analysis. The dashed blue line
dark halos. Each point gives the χ2 value from a single comes from an analysis of the stars only. It is clear that
dynamicalmodel. Thelogarithmicdarkhaloparameters kinematics at large radii are essential to a robust deter-
areplottedontheleft. Thesolidredlineisacubicspline mination of the dark matter fraction at all radii beyond
fit to the lowest χ2 values along the parameter space. the central 0.3 Re.
The dashed blue line shows the χ2 minima coming from The ∆χ2 = 1 range gives us the 68% confidence
justthe stars. Forplotting purposes,anadditive shiftof bands for each of the three parameters. For the loga-
41.5 has been givento the dashed blue line. As the shift rithmic dark halo, the best-fit stellar M/L is 9.1+0.2 (V
−0.2
is additive, the relative χ2 values are preserved. band). The best-fit dark matter halo circular velocity
On the right side in Figure 10 we plot the χ2 values is 800+75 km s−1, and dark matter halo scale radius is
for the NFW models. We do not get a constrainton the −25
radius of 36+7 kpc. The NFW dark halo, while not con-
NFW dark halo scale radius and concentration param- −3
eter. This is evident in the lower-right panel of Figure strained, still gives a robust estimate of the stellar M/L
10 where the χ2 minimum runs unconstrained to rs val- of 8.20+−00..0150. The difference in these stellar M/L values
ues, well beyond the extent of our kinematic coverage. isdrivenentirelybythe shape ofthe assumeddarkhalo,
As the NFW concentration parameter is related to the andthatthe dynamicalmodels workby constrainingto-
scale radius as c = rv/rs, we also do not constrain this tal enclosed mass. As the NFW halo allows for a higher
parameter. central concentration of mass, and the stellar M/L is
A total of 105 stellar LOSVDs and 11 GC LOSVDs assumed constant as a function of radius, mass can be
are used in the dynamical modeling. Of the 105 stellar taken up by the cuspier NFW profile, thus lowering the
LOSVDs, 25 are determined from the 4 SAURON mo- M/L of these models.
ments which provides 25 4 = 100 parameters. The 80
VIRUS-P LOSVDs used×in the modeling are fit to 15 5. DISCUSSION
velocity bins, giving 80 15 = 1200 more parameters. The parameters of the dark halo from this paper are
×
The11GC LOSVDs areconstructedfrom4parameters, different than the ones presented in GT09 which is also
giving another 11 4 = 44 parameters which totals to based on a stellar dynamical analysis. GT09 also fit a
×
1344 for each dynamical model. The best-fit dynamical cored logarithmic dark matter halo yet find a circular
modelforalogarithmichalohada χ2 =1299.4,givinga velocitythat is 10%loweranda scaleradius thatis 60%
reducedχ2 valueof0.97. Theχ2 minimumfortheNFW lowerthan these results. The reasonfor the difference is
halo was 1310.1,which gives a similar reduced χ2 value. due to the datasets; the data presented here have sub-
The constraints on stellar M/L come predominately stantiallyimprovedkinematiccoverageforthestars. The
from the stars, yet the GCs help to constrain the high stellar kinematics of GT09 end at 33′′ whereas our cov-
M/L end as can be seen in the top-left panel of Figure erage extends to nearly 240′′. The GC data is identical
10. This result is not surprising. The GC kinematics betweenthetwopapers. Thelargegapinkinematicspa-
constrain the total enclosed mass in the outer model- tial coverage in GT09 between 33′′ and 140′′ leads to
ingbins;theirkinematicsstronglyinfluencetheresulting generally poor constraints on the dark matter halo pa-
darkmatterhalomass. AsweassumeaconstantM/Lfor rameters. The new VIRUS-P data closes this gapand is
thestars,massnotaccountedforinthedarkmatterhalo therefore more robust.
must get accounted for in the stars and drive the M/L
to higher values. Therefore, kinematics that constrain 5.1. Enclosed Mass
the dark halo will also constrain regions of the model-
The best-fit dark matter halo parameters for a cored
ing where that mass would otherwise wind up, namely
logarithmicprofile returns 800+75 kms−1for circular ve-
higher values for the stellar M/L. −25
In the lower-left panel of Figure 10 we see a different locity, and 36+−73 kpc for the scale radius. In terms of
influence of the GC kinematics stemming from their ex- enclosed mass, M87’s dark matter halo is one of the
tended spatial coverage. Constraints on higher rs values largestevermeasuredforanindividualgalaxy. Figure12
come from the GC kinematics, which extend to 47 kpc. plotsenclosedmassforourbest-fitlogarithmicandNFW
This result is expected, as the stellar kinematics do not models. The black and red lines, with uncertainty, plot
extendouttothedarkhaloscaleradiusandcantherefore total enclosed mass for the logarithmic and NFW mod-
not influence the modeling. Clearly the GC kinematics els respectively. The inclusion of a 6.4 109 M⊙ black
×
areimportantforconstrainingthedarkhaloparameters, hole keeps the total enclosed mass from reaching zero at
andthe goodagreementbetweenthebest-fitstars+GC R=0. Theuncertaintiesarethemin/maxvaluesforthe
model and the stars-only model, where the kinematics 42 = 16 dynamical models that explore the parameter
overlap,isreassuringsinceitimpliesthatbothlargeradii limits of our 68% confidence bands. For the uncertainty
starsandglobularclustersareindynamicalequilibrium. in the black hole mass we use the 0.5 109 M⊙ values
± ×
from GT09. The stellar component is plotted in green
The Dark Matter Halo of M87 9
(dot-dash)withuncertaintieswithin the thickness ofthe and GCs shows good agreement over the range
lines. The dark matter profiles are plotted in gray (long 4 kpc R 20 kpc, yet diverges elsewhere (see Figure
≤ ≤
dash). The yellow vertical line shows the extent of our 14). At both larger and smaller radii the mass profile
kinematic data. The comparison between the enclosed from X-rays is lower than that determined by the stars
mass model from the best-fit logarithmic and NFW ha- andGCs. AtR=3kpctheX-rayestimatedmassisdown
los shows good agreement to the end of our stellar kine- by50%andatR=2kpcthedisagreementis 70%. A
∼
matics. The NFW enclosed mass profile then begins to similar discrepancy is seen at larger radii. The enclosed
diverge to lower total enclosed mass to the end of our massfromX-raysatR=47kpc,ourfurthestdatapoint,
kinematiccoverage. Wediscusshowourresultscompare islowerby50%thanourbest-fitvalue. Thisdifferenceis
toothermassestimatesforM87inthe followingsection. similar to the one seen in NGC4649 (Shen & Gebhardt
2010). One possible explanation for this discrepancy
comes from allowing for a turbulent component in the
5.2. Comparison to Other Mass Estimates
X-ray gas. A 50% decrease in enclosed mass can be ex-
At larger radii Doherty et al. (2009) have measured plained by a 30% non-gravitational component in the
kinematics of PNe for M87. They find a dark halo con- gas. This am∼ount of of difference is similar to the the-
sistent with the one presented here inside of 500′′, al- oreticalexpectation of Brighenti & Mathews (2001) and
though since their radial range is 400′′ to 2500′′ there has been seen in similar systems (Churazov et al. 2010).
is not much spatial overlap with our current data set. More analysis on a wider set of galaxies is necessary to
At around 600′′ they find that the mass density begins fully understand the source of these differences.
to decrease strongly, leading to a truncation of M87’s InTable4andFigure14wecompareenclosedmasses-
dark halo. At R = 1500′′, their outermost radial bin, timatesfromtheliteraturetothiswork. Ourlogarithmic
the PNe dispersion they measure is 78 25 km s−1. For andNFW halo mass profiles are plotted as in Figure 12.
the spatial overlap between our work±and theirs (400′′ Eachcolored symbol in Figure 14 indicates the methods
to 540′′), where we are now comparing globular clusters employedtodeterminetheenclosedmass. Ingeneral,we
and PNe, the kinematics disagree. Possible reasons for find a more massive dark halo for both our logarithmic
the disagreement are that the GC kinematics in this re- andNFWparameterizations,althoughourenclosedmass
gionarepoorly measuredor that the GCs arenot indy- valuesatvariousradialpositionsarenotconsistentlythe
namical equilibrium (e.g. from a recent merger event). highest reported in the literature and appear consistent
Both Cohen (2000) and Cˆot´e et al. (2001) find the GC with the scatter of the data seen in Figure 14.
population around M87 shows both chemical and kine-
matic evidence for two distinct populations of GCs. An- 5.3. Stellar Anisotropy
otherpossibilityisthatthePNemeasurementsarebiased Themechanismsbywhichmassaccumulationoccursin
in some way. Doherty et al. (2009) exclude 3 of 8 PNe galaxies leave their mark on the distribution function of
for their R = 800′′ bin as intracluster planetary nebula the stars (Lynden-Bell 1967; Valluri et al. 2007). There-
and not tracing the potential of M87. Including these 3 fore, mapping the anisotropy of both the stars and GCs
PNe raises their measured dispersion from 139 km s−1 canaddress questions surroundinggalaxyformationhis-
to 247 km s−1. Certainly a comparison to either GC toryandevolution. Ourorbit-baseddynamicalmodeling
or stellar kinematics at this radialposition would be en- return the stellar orbital structure, which we summarize
lightening. in Figure 15. Plotted is the average velocity anisotropy
Wu & Tremaine (2006) estimate the enclosed mass over the 20 angular modeling bins of both the stars and
of M87 at 32 kpc (35.1 kpc at our assumed distance) GCs. Theuncertaintiesarecalculatedinthesamewayas
to be 2.4( 0.6) 1012 M⊙ using GC kinematics and described in Figure 12 and the text. Within R 0.5 Re
assuming s±pheric×al symmetry. Our mass estimate of thestarsshowradialanisotropy,thenbecomemi≃ldlytan-
3.64+−00..8675 × 1012 M⊙ (logarithmic halo) at this radial gentially anisotropic to the last stellar data point. The
position falls within the uncertainties, yet with an off- excellentagreementbetweenthestarsandGCsinthisre-
set of 34%. Romanowsky & Kochanek (2001), us- gionisindicativeofdynamicalequilibriumbetweenthese
∼
ing stellar kinematics from van der Marel (1994) and twocomponents. Althoughwedonotconductadetailed
Sembach & Tonry (1996), and GC kinematics from sev- analysis of the anisotropy of M87 here, comparison of
eralsources,deriveanenclosedmassprofileforM87that anisotropymapstoN-bodysimulationscanbehighlyin-
shows a similar offset towards lower total mass over the formative. An example of such an analysis can be found
range1Re R 5Re. Within1Retheirmodelsdiverge in Hoffman et al. (2010) where the dynamical modeling
≤ ≤
to 50%lowertotalmass. Thisdiscrepancymaybedue ofNGC4365by van den Bosch et al. (2008) is compared
∼
tothestellarkinematicsusedoverthisradialrange. The N-body simulations.
stellarkinematicsofSembach&Tonryexhibitasystem-
5.4. A Discussion of Systematic Uncertainties
atic offset from other data sets for which Romanowsky
& Kochanek make a correction. The offset between our Given the high S/N of our data, we pay particular
enclosedmass andtheirs within 1 Re may be due to this attention to quantifying systematic uncertainties, since
effect or due to the different modeling assumptions, as they might be important for the reported uncertainties.
Romanowskyet al. assume sphericalsymmetry for their As we are using ∆χ2 to determine the parameter values
modeling. The discrepancy may also come about due to and uncertainties, if we do not have proper uncertainty
theincreaseinspatialcoveragetheVIRUS-Pdataaffords estimates for the kinematics we will bias our final mod-
over their long-slit spectroscopy. elingresults. Therearethreeinternalconsistencychecks
Comparison of the X-ray mass determination from that demonstrate that our uncertainties are properly es-
Das et al. (2010) to our mass profile from stars timated.
10 Murphy, Gebhardt & Adams
First, we estimate LOSVDs and Gauss-Hermite pa- parameterizations with the same data sets and model-
rametersfrom4differentwavelengthregions. Comparing ing methods to determine which, if any, is favored. This
the standard deviation across the four regions to the in- requires us to push the collection of stellar kinematics
dividual uncertainties from the Monte Carlo simulations to ever larger radii. The amount of observation time
provides a consistency check. We find that, in general, needed to reach to 2.4 Re with VIRUS-P was not sub-
these two uncertainties estimates are consistent. The stantial,andstellarkinematics to3 and4 Re areachiev-
large wavelength range of VIRUS-P provides this very able. Thesedatawouldallowforbothbetter constraints
important estimate, which includes both statistical and on the various dark matter halo parameters and a com-
systematic effects. parison with the other dynamical tracers (i.e. GCs and
Second,thereducedχ2forthebest-fitdynamicalmod- PNe). Asmuchofourcurrentunderstandingofthedark
els is nearunity. The χ2 is measuredfromthe LOSVDs, matter halos around elliptical galaxies depends on GC
andwecanseetheagreementintheplotofobservedand and PNe kinematics, a robust comparison between each
modeled moments (Figure 7). The deviations between tracer is needed to explore systematics.
the data and the modeling moments are consistent with A second avenue of exploration comes from the
thestateduncertainties. Whilethisconsistencydoesnot information contained in the stellar chemical abun-
directly show that systematic effects are not an issue, it dances available through a Lick index analysis.
is an indirect confirmation. Graves & Schiavon (2008) provide a publicly available
Third,whencomparingkinematicsfromdatasetstaken tool that is well-suited for this work. How elliptical
at different times we find consistent results within the galaxies formed and whether their stars were formed in
stated uncertainties. With the spatial overlap of our situoraccretedovertime requiresbothadynamicaland
pointings #3 and #4 we are able to compare the re- chemical analysis (Graves & Faber 2010). The chemi-
sulting kinematics from four of our spatial bins when cal composition of GCs at large radii have been studied
taken a year apart. We have compared the first four (Cohen 2000; Cˆot´e et al. 2001), and a detailed compar-
Gauss-Hermite moments, calculated from our extracted ison of both the kinematics and chemistry of both GCs
LOSVD,andfindthattheyareallconsistentwithintheir and stars at the same radial position should prove im-
stateduncertainties. Thesethreeinternalchecksdemon- mensely fruitful.
strate control of the measured uncertainties. Finally, work towards a more complete and uniform
Next, we discuss the two areas where systematic ef- sample of massive elliptical galaxies, both 1st and 2nd
fects may be important: sky subtraction and template rank galaxies, and equally massive field ellipticals (e.g.
mismatch. Inordertodeterminehowthelevelofskysub- NGC1600)is needed to explore the influence of environ-
tractionaffectsourextractedkinematicsandsubsequent ment on dark matter halos. Several groups have made
modeling,weexplorebothoverandundersubtractionof significant progress towards this end, yet the data sets
each 20 minute science frame. A range of sky subtrac- that involve both 2D spatial coverageat both small and
tion levels are created and taken through all subsequent large need to be expanded.
data reductions. A total of 25 different sky subtractions
are made on each science frame, over a range of 12.5% 6. ACKNOWLEDGMENTS
±
when compared to equal exposure times. The details of
The authors would like to thank The Cynthia and
these reductions are given in the Sky Subtraction sec-
George Mitchell Foundation for their generous support
tion in the Appendix. We then compare the calculated
whichmadethefabricationofVIRUS-Ppossible. J.D.M.
velocity and velocity dispersion values, taken from the
wouldliketothankthemembersofhisresearchcommit-
best-fit LOSVD. This comparison, over all 88 spectra,
tee: Michele Cappellari, Gary J. Hill, John Kormendy,
shows no systematic offsets in velocity or velocity dis-
Phillip J. MacQueen, and Milo˘s Milosavljevi´c for all the
persion for either over or under subtraction of the night
helpgiven. J.D.M.alsothanksGuillermoBlanc,RossE.
sky. The associated random errors for this full range
Falcon, and Remco van den Bosch for their insight and
of sky subtractions is within our quoted uncertainty for
many fruitful discussions. The authors also thank Dave
both velocity and velocity dispersion.
Doss, Kevin Meyer, Brian Roman, John Kuehne and all
In order to explore possible systematics due to our
of the staff at McDonald Observatory who helped im-
use of the Indo-US spectral library, we select the same
mensely with the commissioning of VIRUS-P and the
set of template stars from the Miles spectral library
successful collection of this data. This research is par-
(Sa´nchez-Bl´azquez et al. 2006) and extract kinematics
tiallybasedondatafromthepubliclyavailableSAURON
for all of our spectra. The two libraries agree very well,
Archive,whichaidedsubstantially in the modeling. The
with deviations between the libraries of 2.5km s−1,
authors thank Payel Das and her collaborators for their
∼
well within our quoted uncertainties for velocity disper-
willingness to share their X-ray data and analysis on
sion. In the case of velocity, there is a slight offset
M87. This work would not have been possible without
( 7km s−1) which is due to the lack of a velocity zero-
the resources of the Texas Advanced Computing Cen-
∼
point between the two libraries. Both of these checks
ter at the University of Texas, Austin where all of the
indicate that our systematics are under control.
dynamical modeling was run.
5.5. Next Steps
This work points the way to several other areas of in-
quiry. We have exploredtwo different parameterizations
for a dark matter halo, yet others exist and there is no
reason to dismiss any of them. A natural next step is
to rigorously explore several different dark matter halo