Unravelling the origins of S0 galaxies using maximum likelihood analysis of planetary nebulae kinematics

0.72 MB·English

by  A. Cortesi

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed27January2011 (MNLATEXstylefilev2.2)
Unravelling the origins of S0 galaxies using maximum
likelihood analysis of planetary nebulae kinematics
1
1 A. Cortesi1,2⋆, M. R. Merrifield1, M. Arnaboldi2, O. Gerhard3, I. Martinez-Valpuesta3,
0
K. Saha3, L. Coccato3, S. Bamford1, N. R. Napolitano4, P. Das3, N. G. Douglas5,
2
A. J. Romanowsky7, K. Kuijken6, M. Capaccioli8 and K. C. Freeman9
n
a 1Universityof Nottingham, School of Physics and Astronomy, UniversityPark, NG7 2RD Nottingham, UK
J 2European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany
6 3Max-Planck-Institut fu¨r Extraterrestrische Physik, Giessenbachstrasse, 85741 Garching, Germany
2 4Istituto Nazionale di Astrofisica, Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Naples, Italy
5KapteynAstronomical Institute, Universityof Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
] 6Leiden Observatory,Leiden University,PO Box 9513, 2300 RA Leiden, The Netherlands
A 7UCO/Lick Observatory, Universityof California, Santa Cruz, CA 95064, USA
G 8Dipartimento di Fisica, Universita` “Federico II”, Naples, Italy
9Research School of Astronomy and Astrophysics, Australian National University,Canberra, Australia
.
h
p
- ;
o
r
t
s
ABSTRACT
a
[ To investigate the origins of S0 galaxies, we present a new method of analyzing their
stellar kinematics from discrete tracers such as planetary nebulae. This method in-

1 volvesbinningthedataintheradialdirectionsoastoextractthemostgeneralpossible
v
non-parametrickinematicprofiles,andusingamaximumlikelihoodfitwithineachbin
2
in order to make full use of the information in the discrete kinematic tracers. Both
9
disk and spheroidkinematic components arefitted, with a two-dimensionaldecompo-
0
sitionof imagingdata used to attribute to eachtracera probabilityof membershipin
5
. the separate components. Likelihood clipping also allows us to identify objects whose
1
properties are not consistent with the adopted model, rendering the technique robust
0
against contaminants and able to identify additional kinematic features.
1
ThemethodisfirsttestedonanN-bodysimulatedgalaxytoassesspossiblesources
1
ofsystematicerrorassociatedwiththestructuralandkinematicdecomposition,which
:
v arefoundtobesmall.ItisthenappliedtotheS0systemNGC1023,forwhichaplan-
i etary nebula catalogue has alreadybeen releasedand analyzedby Noordermeer et al.
X
(2008). The correct inclusion of the spheroidal component allows us to show that,
r
a contraryto previousclaims,the stellar kinematicsofthis galaxyareindistinguishable
fromthoseofanormalspiralgalaxy,indicatingthatitmayhaveevolveddirectlyfrom
such a system via gas stripping or secular evolution. The method also successfully
identifies a population of outliers whose kinematics are different from those of the
main galaxy;these objects can be identified with a stellar stream associated with the
companion galaxy NGC 1023A.
Key words: galaxies: elliptical and lenticular – galaxies: evolution – galaxies: in-
dividual: NGC 1023 – galaxies: individual: NGC 1023A – galaxies: kinematics and
dynamics.
1 INTRODUCTION disk components associated with spirals. They thus poten-
tially represent a key element in attempts to understand
Lenticular or S0 galaxies lie between ellipticals and spirals
the relationship between the main types of galaxy, but at
on theHubblesequence,since they havethe featureless old
themomentthereisnoconsensusastotheendoftheHub-
stellarpopulationsofellipticalsystems,butalsocontainthe
blesequencetowhichtheyaremostcloselyrelated.Clearly,
someprocesshasconsumedtheirgasandshutofftheirstar
formation. However, this termination could bethe result of
⋆ email:[email protected] galaxy mergers, much in the manner that star formation is
2 Cortesi et al.
believed to be quenched in elliptical galaxies, or it could therunrelated backgroundobjects, orliewithin thesystem
arise from some muchgentler process, suchasram pressure but donot match theexpected kinematics.
stripping, which simply removes thegas from normal spiral Inthispaper,wethereforepresentasomewhatmoreso-
galaxies. phisticatedanalysisdesignedtocircumventtheseshortcom-
ings.Whilestillbinningthedataintheradialdirection,soas
Perhaps the best way of discriminating between these
toextractageneralnon-parametricviewofthegalaxy’sdy-
scenarios is offered by the stellar dynamics of S0 galaxy
namics,weapplyamaximumlikelihoodanalysiswithineach
disks. If they are the product of relatively gentle gas-
radialbin,soastoextractthemaximumamountofinforma-
stripping processes, one would expect the stellar dynam-
tion from the azimuthal variations in kinematic properties.
ics to be unaffected, and should be identical to those of
We also model both disk and spheroid in the system, and
the progenitor spiral system, with rotation dominating rel-
useafulltwo-dimensionalphotometricdecompositionofthe
atively modest amounts of random motion in the disk
galaxy,whichallowsustoallocatetoeachPNaprobability
(Arag´on-Salamanca2008;Williams et al.2010).If,however,
as to thecomponent to which it belongs.
theS0saretheproductofminormergers,withamassratio
The method we have developed to construct such a
higher than 10:1, one would expect the merging process to
model is described in Section 2, and in Section 3 we test
beimprintedinthestellardynamics,withrandomvelocities
the technique on a simulated galaxy for which we have
as high or even higher than the rotational velocities at all
priorknowledgeofthestructuralcomponentsandkinemat-
radii (Bournaud et al. 2005).
ics. In Section 4 we make a first real application of the
One difficulty in using such a diagnostic is that it re- methodto theNGC1023 catalogue, in which weshowthat
quires accurate stellar kinematics measured in the outer this likelihood fitting can reproduce the strange results in
parts of the galaxy where the disk dominates, and at Noordermeer et al. (2008) if we also neglect the spheroidal
lowsurface brightness,makingconventionalabsorption-line component, but that the more complete model presented
spectroscopy difficult. However, planetary nebulae (PNe) here results in a picture in which NGC 1023 has a normal
provide an excellent probe of stellar kinematics in this rotationally-supported disk. We also illustrate the power of
regime: they are easily detected and can have their ve- likelihoodfittingbyshowinghowoutlierpointscanbeiden-
locities measured from their characteristic strong emission tifiedinthecatalogue,andassociatedwithaminoron-going
lines, and they have proved to be an unbiased tracer of merger. The conclusions are presented in Section 5.
the global stellar population (Ciardullo et al. 1989, 2002;
Napolitano et al. 2001; Coccato et al. 2009). We have de-
signedandbuiltaspecial-purposeinstrument,thePlanetary
2 KINEMATIC LIKELIHOOD FITTING
Nebula Spectrograph [PN.S; Douglas et al. (2002)], specifi-
callytoexploitthistracerofstellarkinematicsatlowsurface In this section we give a brief introduction to likelihood
brightnesses. Although originally intended to study ellipti- analysis, before going on to the specifics of fitting partic-
cal galaxies (Romanowsky et al. 2003; Coccato et al. 2009; ular galaxy models. In Section 2.1, we present the case of
Napolitano et al.2009),PN.Shasalreadyprovedveryeffec- a pure disk galaxy, as it provides both a simple example of
tive in exploring the disks of spiral galaxies (Merrett et al. the technique and a useful point of comparison to previous
2006), so an obvious next step was to use it to probe the analyses that have assumed that the disk is the dominant
kinematics of S0systems. component.InSection 2.2,weproceed tothemorerealistic
situation in which we model both disk and spheroidal com-
As a pilot study, we observed the relatively nearby S0
ponent.InSection2.3wediscusshowtheindividualdiscrete
galaxy NGC 1023 with PN.S. The observations and the re-
kinematictracerscanbeassigned,atleastinaprobabilistic
sulting catalogue of positions and line-of-sight velocities of
sense, to either the disk or the spheroid using photometric
PNearedescribedinNoordermeer et al.(2008).Theanaly-
data.
sisinthatpaperconcludedthatNGC1023hasverypeculiar
In general, given a set of N independent velocity mea-
kinematics, giving theappearance of anormal rotationally-
supporteddiskgalaxyatsmallradii,buthavingentirelyran- surements, vi = (v1,...,vN), drawn from a probability den-
dom velocities at large radii, inconsistent with either of the sityfunctionF(vi;θ~),whereθ~=(V,σ)isasetofparameters
whose valueis unknown,thelikelihood function is given by
expected scenarios. However, the relatively simple dynami-
cal analysis that led to these conclusions had some signifi- = F(vi;θ~). (1)
cant shortcomings. First, it assumed that the light at large L Yi
radii was dominated by the disk component, and neglected
ThevaluesofV andσthatmaximize arethebestestima-
any contribution from the spheroid. As we will see, this as- L
tors for the true values of these parameters. Moreover the
sumptioncanleadtosizablesystematicerrorsintheinferred
method of maximum likelihood coincides with the method
disk kinematics. Second, it calculated dynamical quantities
of the least squares in the special case of a set of N Gaus-
such as velocity dispersions by binning the data both az-
sian distributed independent measurements, in which case
imuthallyandradiallyinthegalaxy.Althoughsuchbinning
the likelihood function is directly related to the usual χ2
provides an adequate signal from a sparsely-sampled veloc-
statistic,
ity field and a non-parametric measure of the local kine-
matics,aswewillshowbelowthekinematicpropertiesvary ∆χ2(~θ)= 2∆ln (~θ). (2)
− L
continuously and quite rapidly with azimuth, so such bin-
Thus,inthiscaseitisstraightforwardtodeterminethecon-
ningaveragesawaysignificantamountofinformationthatis
fidenceregion around thebest estimator:
present in the raw data. The binning process also makes it
difficulttoidentifyanycontaminationfromPNethatareei- ln (~θ)>ln max ∆ln , (3)
L L − L
The origins of S0s from PN kinematics 3
wherevaluesof∆χ2or 2∆ln ,correspondingtoadesired galaxy, liketheones we are focusing on, it does not project
− L
confidence limit, for joint estimation of m parameters, are much into the line-of-sight velocity dispersion, so the value
readily available in tabulated form. In this paper we use ofσ isdominatedbytheothertwocomponents.Ascanbe
los
a ∆χ2 that corresponds to a 1σ coverage probability. In seen from the above equation, σ2 varies sinusoidally such
los
particular, ∆χ2 = 5.39 for m = 4, and ∆χ2 = 4.11, for that its value is set by σR on the minor axis and σφ on
m=3. the major axis. Thus, by fitting the variation in σ with
los
Finally, once the best-fit parameters have been found, azimuthonecandeterminethevaluesofbothofthesecom-
onecan calculatethecontributionofeachdatapointtothe ponentsof the galaxy’s intrinsic velocity dispersion.
likelihood given this optimum model. To quantify this con- If we adopt the simplest possible model in which the
tribution in terms of whether the kinematic tracer in ques- line-of-sight velocity distribution isGaussian ateverypoint
tionisconsistentwithbeingdrawnfrom thebest-fitmodel, withameanvelocityofV andadispersionofσ ,wenow
los los
one can generate a large number of individual velocities havetherequisiteform fortheprobabilitydensityfunction,
drawnfromavelocitydistributionmatchingthemodel,and
Wseeeactanwthhautspiedrecnetnitfyileanofytohbejedcisttsrwibhuotsioenvethloecidtaietsalpieoionuttlsiiedse. F(vi;V,σr,σφ)∝exp(cid:20)−2[v[σil−os(Vσlors,(σVφ))]]22(cid:21). (6)
this confidence interval, which are most likely interlopers.
By rejecting these data points and iterating on the fit, we ThevaluesV,σr,σφthatmaximizeLarethebestestimators
for the true values of these kinematics parameters within
can render this process robust against a small amount of
each bin.
contamination,inastraightforwardgeneralization ofthe3σ
clipping procedure (Merrett et al. 2003) used in determin-
ing membership of velocity distributions that are assumed
tobeGaussian. Wedefinethelikelihood clipping probability
2.2 Likelihood analysis for a disk + spheroid
threshold asthevalueat which wecutthedistribution,and
model
we reject all thePNe beyond thislimit.
In principle, such likelihood fitting could be applied to Inadditiontothediskcomponent,mostsystemsalsousually
a full dynamical model across an entire galaxy, at all radii contain a spheroidal component comprising either a central
andazimuths.However,inthisanalysisweareinterestedin bulgeoranextendedhalo(orboth).Indeed,oneofthetra-
leaving as much freedom as possible in the resulting radial ditional defining features of S0 galaxies is that they often
profileofdynamicalproperties,soastoexplorethefullrange have rather prominent bulges. A kinematic model without
ofpossibilitiespredictedbydifferentmodelsofS0formation. suchacomponentmaythereforebesignificantlyincomplete;
Wethereforeadoptahybridapproachinwhichthedataare worse,theamountofbulgelightvarieswithazimutharound
binnedintoellipticalannulimatchedtotheinclinationofthe thegalaxy,sinceclosetotheminoraxisinahighly-inclined
disk, to extract points over a range of intrinsic radii. This systemthebulgeisusuallythedominantcomponent,soone
binning allows each section of the disk to be modeled as might expect the variation in velocity dispersion with az-
an independent non-parametric data point. However in the imuth, used above to extract the different components of
azimuthaldirectionaroundeachbinweuseafulllikelihood disk velocity dispersion, to besystematically distorted.
analysisthataccuratelymodelstheexpectedvariationinthe One slight complication in adding in such a spheroidal
line-of-sight velocity distribution with azimuth, as derived componentisthatitdoesnothavethesamesymmetryprop-
from very general dynamical considerations. ertiesasthediskcomponent.Thus,theellipticalannulithat
containdatafromasmallrangeinradiiinthediskdoesnot
correspond to thesame rangein radii in themore spherical
2.1 Likelihood analysis for a disk model
spheroidal component, but samples a larger range of radii
Consider a general rotating disk model for a galaxy. Part both due to the difference in shape and the effects of in-
of the line-of-sight velocity of each object within it, vi, is tegration along the line of sight. However, the variation in
the projection of the galaxy’s mean rotation velocity V, at kinematicswithradiusinsuchahotcomponentisgenerally
an azimuthal angle φ within the galaxy, corrected for the rather slow, so the greater averaging in radius imposed by
inclination of the galaxy to the line of sight i, and for the the choice of elliptical annuli should not have a major im-
systemic velocity of the galaxy itself Vsys: pact. Further, the validity of this assumption of slow vari-
ation with radius can be tested a posteriori by seeing how
Vlos=Vsys+V sin(i)cos(φ). (4) the inferred kinematics of the spheroid change from bin to
Superimposedonthenetrotationalvelocityaretheran- bin.
dom motions of the individual stars, which can be quanti- Accordingly, we adopt the simplest possible model for
fiedbytheirvelocitydispersionindifferentdirections.Foran the kinematics of the spheroidal component, in which the
axisymmetricdisk,thesecomponentsaremostnaturallyex- line-of-sight velocity distribution is assumed to be a Gaus-
pressedincylindricalpolarcoordinatesalignedwiththeaxis sian with zero mean velocity (relative to the galaxy). We
of symmetry,(R,φ,z), and theline-of-sight velocity disper- havethusassumedthatanyrotationinthespheroidalcom-
sion,σlos,ismadeupfromaprojectionofthesecomponents ponent is negligible, although, as we shall see below, this
such that assumptioncanalsoberelaxed.Thevelocitydispersion,σs,
2 2 2 2 2 2 2 2 2 is left as a free parameter to be modeled by the data. The
σ =σ sin isin φ+σ sin icos φ+σ cos i (5)
los r φ z only other parameter that we need to specify is the proba-
Generallyspeaking,indiskgalaxies,σz isthesmallestcom- bility that each individual tracer object at its observed lo-
ponent of the velocity dispersion, and in a nearly edge-on cation belongs to thespheroidal component, fi, so that the
4 Cortesi et al.
full probability density function can bewritten 3 APPLICATION TO A SIMULATED GALAXY
F(vi;V,σr,σφ,σs)∝ σfis exp(cid:20)−2vσi2s2(cid:21) Amulsaotdaieolnfi.grasTltahxetyessitomboutfalaintthieoidns fwmroaemsthpoardos,veilwdf-eecdopnbesyirsftKoernamntaNtkh-beSaofihdtayt(sopimria--
+1−fi exp (vi−Vlos)2 , (7) vate communication). This simulated galaxy is constructed
σlos (cid:20)− 2σl2os (cid:21) using a nearly self-consistent bulge-disk-dark halo model
(Kuijken & Dubinski1995; Widrow et al. 2008) to mimic a
where the velocity dispersions in the denominators of the typical lenticular galaxy. The value of Toomre Q is high
term in front of each Gaussian ensure that this function enough to prevent strong 2-armed spirals from forming in
is properly normalized. From this velocity distribution one thedisk.Thespheroidfollows aSersicprofilewithanindex
can then construct a likelihood function for a particular of 3.5, while the disk is exponential with a vertical sech2
set of model parameters using Equation 1. The values of profile. The dominant dark matter halo specifies the sys-
V,σR,σφ,σs that maximize then provide the estimators tem’s rotation curve. We chose to “observe” this galaxy at
L
for these parameters in the best-fitmodel. an inclination of 60 degrees. For the following tests, we set
thelikelihood clipping probability threshold to 2.1σ.
3.1 Testing the likelihood method using a priori
2.3 Spheroid–disk decomposition knowledge of PN positions
The only parameters that we havenot yet specified are the In a simulated galaxy, we have access to inside informa-
fi values, which define the probability that a given object tion which allows us to test specific aspects of this mod-
is in the spheroid. We cannot solve for these quantities di- eling process. In particular, the “stars” are tagged accord-
rectly in thelikelihood maximization since, asnoted above, ing to whether they are members of the disk or spheroid,
therelativecontributionsofspheroidanddiskvarywithaz- so we can assign the objects to specific components with
imuthwithin a single bin,so each kinematicdatapoint has certainty, rather than using the probabilistic approach de-
its own individual value. Fortunately, we have additional scribed above. Accordingly, we have carried out the likeli-
information that has not yet been used. In particular, we hoodanalysiswithouttheprobabilisticdecomposition,with
can apply a two-dimensional galaxy fitting routine such as the results presented in Figure 1. The upper panel shows
GALFIT(Peng et al. 2002) to imaging datain order tode- the results obtained with a generous catalogue of 437 ob-
compose the starlight into spheroid and disk components. jects, while the lower panel shows how well we can do with
The fit to the light distribution then provides a direct esti- only 136 kinematic tracers. Although not strictly valid for
mate of the fraction of the starlight from each component kinematically-hot S0 galaxies, an asymmetric drift correc-
at thelocation of any stellar tracer. tion hasbeenapplied toconvert thederivedrotation veloc-
This approach for determining the decomposition into ity into a circular velocity for direct comparison with the
bulgeanddiskcomponentshasthegreatadvantagethatthe simulation’s known functional form, using theformula
broad-band imaging data offer a much less sparse sampling
of the stellar spatial distribution than the PNe, providing V2=V2+σ2 σ2 1 rh + dlnσr2 , (8)
an intrinsically more accurate answer. Furthermore, there c φ− r(cid:18) − rD dlnr (cid:19)
are significant selection effects in detecting PNe, as they
are harder to find against the bright background of the in- where Vc is the circular velocity, rh is the median radius in
ner parts of a galaxy, so their spatial arrangement is not each bin,and rD is the disk scale length (Kormendy1984).
an unbiasedrepresentation of thestellar distribution.Note, The final gradient term of the equation has been estimated
however, that this bias is a purely spatial one, in that all assumingthatσ2followsanexponentialprofileandperform-
r
line-of-sight velocitiesareequallydetectable,sotheirusein ing a linear fit between lnσ2 and theradius.
r
thekinematic analysis is not in any way compromised. Whiletheerrorsare,asexpected,largerforthesmaller
The only slight subtlety in applying such analysis to sample, there are no systematic differences between them,
PNe is that their number per unit stellar luminosity has and both do a good job of recovering the simulation’s cir-
been shown to vary systematically with the colour of the cular velocity. It is interesting to see that σR is moderately
population (Ciardullo et al. 1991), with alower PN density but systematically too low for the small sample-size, and
perunitgalaxyluminosityforredderobjects(Buzzoni et al. remains low at large radii even for the larger sample. Since
2006).Ifthereisadifferenceincolourbetweenthediskand the main remaining assumption in this model is that the
spheroidal components, then one has to apply a correction velocity distributions of the individual components are in-
in order to convert the fraction of spheroidal light at any trinsicallyGaussian,itseemslikelythatthismodestsystem-
given point into theprobability that a PN detected at that atic error occurs due to the breakdown of this assumption.
pointbelongstothespheroid.Inpractice,suchcolourterms However, we note that σ and v are both quite accurately
φ
can be straightforwardly determined by performing the de- recovered,sowecanestimatewithsomeconfidencethebal-
composition on images taken in different bands, and using ancebetweenrandomandrotationalmotion,σ /v,whichis
φ
the resulting colours of the different components to correct the main physical quantity that we are seeking to use as a
theprobability. diagnostic in this analysis of the origins of S0 galaxies.
The origins of S0s from PN kinematics 5
Figure 2. Results from the isophotal analysis of the simulated
galaxy,blackfilledtriangles,andthegalaxymodelsobtainedwith
GALFITleavingtheSersicindexnasafreeparameter,redfilled
pentagons, imposingSGALFIT =3.5,greenopenpentagons and
imposingSGALFIT =4,bluestarredsymbols.Thetoprowshows
the fitted position angles, the middle row the ellipticity and the
bottom rowthesurfacebrightness.
ing simulated broad-band image. The principal parameters
returned by this fitting process are the spheroid and disk
magnitudes, their respective scale lengths, the Sersic index
ofthespheroid,andtheflatteningandposition angleofthe
components.
Interestingly, if the parameters are all left free, then
the spheroid is found to have a best-fit Sersic index of 2.8,
significantlysmallerthantheknownvalueforthissimulation
of 3.5. The value of the effective radius is also found to be
systematically too small. However, an almost equally-good
fit is found if we fix the Sersic index to 3.5. In fact, this
fit is by some measures superior: Figure 2 shows the values
Figure 1. Derived circular velocity and the components of the ofposition angle,ellipticity andsurfacebrightnessobtained
velocitydispersionversusradiusforthemodelgalaxyinthedisk by fitting elliptical isophotes to both the simulated galaxy
and in the bulge, for 437 particles, upper panel, and 136 par- image and the models in which the Sersic index is either
ticles, bottom panel.The filled symbols are from the maximum- fixedorleft free. Whilebothmodels reproducetheposition
likelihoodanalysis,withverticalerrorbarsindicatinguncertainty angle and surface brightness equally well, the ellipticity is
and horizontal error bars showing the extent of each radial bin
clearlybetterfittedbythemodelwiththeSersicindexfixed
(with the point plotted at the median radius for a PN in that
at the right value. This conflicting information underlines
bin). The dashed lines represent the model circular velocity, ve-
the complexity of non-linear model fitting, and illustrates
locitydispersioninthediskandinthespheroid.
its basic limitations.
Of course, for a real galaxy, we would not know the
“right” value for the Sersic index, so would not be able
3.2 Testing the effect of the photometric
to choose between these models. Such systematic errors in
spheroid–disk decomposition
GALFITfitting therefore pose a potential limitation to the
We now introduce the remaining aspect of the full model- effectivenessofthemodelingproceduresetoutinthispaper,
fitting process, the spheroid–disk decomposition, which al- iftheresultingkinematicsturnouttodependsensitivelyon
lows us to assign a probabilistic membership of each kine- the spheroid–disk decomposition. To assess the impact of
matic tracer to each component, since, as discussed above, such effects, we carried out the maximum likelihood mod-
in a real galaxy we do not have the luxury of this knowl- eling using both of these decompositions. We also checked
edge a priori. We constructed an image of the simulated theeffectofkinematictracersamplesizebysimulatingcat-
galaxy and convolved it with a suitable Gaussian point- alogues of 437 and 136 PNe. The resulting kinematics are
spreadfunction.WethenusedGALFITtomodeltheresult- presented in Figure 3. The good news is that the results
6 Cortesi et al.
Figure 5. Derived mean rotation speed and the components of
the velocity dispersionversus radius for NGC 1023, for the case
ofadisk-onlymodel.Thefilledsymbolsarefromthemaximum-
likelihoodanalysis,withverticalerrorbarsindicatinguncertainty
and horizontal error bars showing the extent of each radial bin
(with the point plotted at the median location for a PN in that
bin). The open symbols reproduce the results of the tilted-ring
analysisbyNoordermeeretal.(2008).
Figure 4. Positions and velocities of the PNe detected in
NGC 1023. Colors bar indicates velocities. The elliptical annuli 4.1 The disk model
showtheradialbinningofthedata.Thecircledpointsarethose
rejectedinthepurediskmodel,withthelikelihoodclippingprob- The application of the above approach, in which the data
ability threshold increased to 1.65σ. Squares show PNe associ- are binned radially (into the annuli indicated in Figure 4)
atedwiththecompaniongalaxyNGC1023A(Noordermeeretal. butmodeledazimuthallybylikelihoodanalysis,ispresented
2008). in Figure 5.For thisanalysis and thefollowing sections, we
have estimated the galaxy’s inclination from the ellipticity
of thedisk component derivedfrom thetwo-dimensional fit
are largely insensitive to systematic errors arising from the tothephotometry(seeSection4.2);assumingthatthedisk
disk – spheroid decomposition. For the smaller catalogue, is intrinsically axisymmetric, we infer a value of i = 74.3
the extra uncertainty that arises from the decomposition degrees. Likelihood clipping has been applied such that
processmeanthattheerrorsinthederivedkinematicquan- PNe with a probability of less than 2.1σ of being drawn
tities become quite large, but there is no evidence for any from the disk model are rejected. As is clear from this fig-
systematic error in the process. Thus, it appears that this ure, the analysis largely reproduces the peculiar result of
maximum likelihood procedure is quitereliable, and robust Noordermeer et al. (2008), with an inferred rotation speed
against the most likely sources of systematic error in the thatfalls rapidly outside300 arcseconds, accompanied bya
spheroid–bulge decomposition. rise in velocity dispersion. As previously noted, such kine-
matics were not predicted by any of the simpler scenarios
for S0 formation, therefore motivating this further investi-
gation.
4 APPLICATION TO NGC 1023
One immediate clue is offered by the only significant
Havingdevelopedthismethodologyforextractingkinematic differencebetween theresultsofthetwoanalyses presented
properties of disk galaxies from individual stellar tracers, in Figure 5. Specifically,while theradial dispersion exceeds
and tested its reliability, we can now apply it to the ex- the tangential component in the Noordermeer et al. (2008)
isting PNe data for the S0 galaxy NGC 1023. The cat- analysis, the opposite is the case in the maximum likeli-
alogue of 183 PNe positions and velocities, published by hoodfit.Infact,theorderingofdispersionswasfixedinthe
Noordermeer et al. (2008), are presented in Figure 4. Note firstanalysis, as theirratio was set bytheepicyclicapprox-
the “hole” in the middle of the galaxy resulting from the imation, which forces σr > σφ unless the rotation curve is
difficulty of detecting PNe against the bright stellar conti- very rapidly rising (Binney & Tremaine 1987). In the cur-
nuum at these radii, as discussed above. Application of the rentanalysis,wehaveleftbothcomponentsasfreeparame-
new method to these data has the benefit that they have ters,and it is telling usthat wethenfindσφ >σr, which is
already been studied using a more conventional tilted-ring not what the physically-motivated epicyclic approximation
model(Noordermeer et al.2008),withwhichourresultscan would produce, suggesting that there is some basic flaw in
be compared, plus the somewhat peculiar kinematics ap- themodel.
parentlyfoundinthissystemwarrantfurtherinvestigation. As noted above, the principal missing element in this
SincethepreviousanalysisbyNoordermeer et al.(2008)ne- modelisthespheroidalcomponent,whichisconsistentwith
glected any contribution to the kinematics from the bulge, the enhanced value of σ . Specifically, a spheroid will con-
φ
webeginbyconsideringthedisk-onlymodeltomakeadirect tribute a relatively small number of PNe with a large dis-
comparison. persion but zero mean velocity. Close to the minor axis of
The origins of S0s from PN kinematics 7
Figure 3. Derived circular velocity and the components of the velocity dispersion versus radius for the model galaxy in the disk and
in the bulge. The left panels show the result when the spheroid’s Sersic index is fixed at its true value of 3.5, while the right panels
show theresultswhenitisleftas afreeparameter. Upperpanels areforalargercatalogue of437kinematictracers;lower panelsfora
smallersampleof136.Thefilledsymbolsarefromthemaximum-likelihoodanalysis,withverticalerrorbarsindicatinguncertaintyand
horizontalerrorbarsshowingtheextent ofeachradialbin(withthepointplottedatthemedianlocationforakinematictracerinthat
bin).
the galaxy, the disk PNe will also have zero mean line-of- remove the contaminants from the fit. Figure 6 shows the
sight velocity, whereas on the major axis their distribution successofsuchaprocess,inwhichamoreaggressivethresh-
willbeoffsetinvelocityduetorotation.Combiningthetwo oldresultedin34PNebeingrejected.Thecomponentsofthe
distributions with different mean velocities will result in a velocity dispersions are now ordered in the manner physi-
larger incorrectly-inferred dispersion than combining them cally expected for a disk population. Further, the bizarre
wheretheirmean velocities arethesame. Sincethetangen- behavior of the kinematics has entirely vanished, with the
tialcomponentprojectsmostlyintotheline-of-sightcloseto rotationvelocitynowremainingapproximatelyconstantout
the major axis, its derived value would be more enhanced tolarge radii,andthedispersion remaininglow,just asone
by this contamination than that of the radial component, would expect for a normal disk population.
consistent with theresults in Figure 5.
Asfurtherevidencethatthecauseofthecontamination
Averysimpletestofwhethersuchcontaminationcould isthespheroid,Figure4highlightsthelocationsofthePNe
be responsible for the unphysical results can be made by rejected in this iterative clipping. The PNe appear spread
increasing the severity of the likelihood clipping, to try to throughoutthegalaxy,andnotflattenedintotheaspectra-
8 Cortesi et al.
Figure6.AsforFigure5,butwiththelikelihoodclippingprob-
abilitythresholdincreasedto1.65σ.
Figure8.Mapshowingtheprobability,f,thataPNatanygiven
location is drawn from the spheroidal component of NGC 1023.
tio of the disk, indicating that they are probably not the ThepositionsofthedetectedPNearesuperimposed.Boldsquares
resultofapoordiskmodel.Thereis,however,someindica- identifyNGC1023APNe.Thef =0.5contourisshowninblack.
tionofanasymmetryinrejectedPNebetweenthetwosides
of the galaxy, which we return toin Section 4.4.
Clearly we need to model a spheroidal component as
the best-fit model disk has a scale length of 67 arcseconds
well as the disk, but it is also heartening to see that the
fortheB-banddataand66arcsecondsfortheR-band;sim-
likelihood rejection method does a respectable job of deal-
ilarly, the spheroid has an effective radius of 24 arcseconds
ing with even such a high level of contamination, in which
in the B-band and 21 arcseconds in the R-band. There is
almost 20% of the PNe do not seem to be drawn from the
thusnoevidenceforcolorgradientswithineachcomponent.
assumedmodel.Thisresultprovidessomeconfidenceinthe
Thislackofcolorvariationbetweenandwithincomponents
robustness of the adopted procedure.
makes the translation from stellar continuum properties to
the probability that a PN belongs to the spheroid or disk
4.2 Spheroid – disk decomposition particularly simple in this case, as the probability is just
theratio ofthecomponent tototallight at each point.The
The first step toward modeling the spheroidal compo-
bulge to total light ratio is 0.31 in the R-band and 0.35 in
nent is to decompose continuum images into disk and
theB-band.Thesevaluescomparewellwithestimatesfrom
spheroidalcomponents,whichwehavedoneusingGALFIT
the cruder one-dimensional bulge–disk light decomposition
(Peng et al. 2002) to fit an exponential disk and a de Vau-
performedin(Noordermeer et al.2008),wherethebulge-to-
couleurs law spheroid to deep images of NGC 1023. Fig-
total light ratio was found to be 0.36. The two-dimensional
ure7 showstheresult ofthisprocess carried out on bothB
map of the fraction of light in the spheroidal component,
and R band images. The companion galaxy NGC 1023A is
f, as a function of position on the sky is presented in Fig-
mildly apparent even in theraw images, but shows up very
ure 8. As discussed in Section 2.3, the value of f clearly
clearly in the residuals when the model is subtracted. The
varies strongly with azimuth, with the spheroid being the
residual image also shows evidence for the bar at the cen-
dominant component at all radii close to the minor axis,
tre of this galaxy noted by Debattista et al. (2002); since
underliningthenecessityofthisfulltwo-dimensionaldecom-
this faint feature is relatively localized at small radii, and
position.
weareprimarilyinterestedinthebalancebetweendiskand
spheroid light at large radii, we do not attempt to model it
anyfurther.Wecouldalsohaveconsideredothermodelsfor
4.3 The disk + spheroid model
the spheroid, such as incorporating distinct bulge and halo
components,orpickingamoregeneral Sersicprofilefor the Having calculated the division between spheroid and disk
lightdistribution.However,aswehaveseeninSection3,the light, we can now carry out the likelihood analysis for
resultsarerelativelyinsensitivetosuchsubtleties.Sincethe NGC 1023 incorporating this extra component, as set out
primarygoalisjusttomodeltoareasonableapproximation in Section 2.2. the resulting best-fit kinematic parameters
thefractionofthelightinthedifferentcomponentsatdiffer- asafunctionof radiusareshown in Figure9.Wenowstart
ent positions, there is little benefit in trying to distinguish to see thecharacteristic properties of a normal disk galaxy.
between what are likely to be near-degenerate alternative Rotationdominatesrandommotionsinthediskatallradii,
fits. and the mean rotation speed stays approximately constant
Fromthemodelsinthetwobands,wecancalculatecol- outtothelastpointsshown.Inthisplot,wehavealsofilled
ors for the individual components, and in this case we find in the missing stellar rotation velocity at small radii from
B R = 1.62 for the disk and 1.61 for the spheroid. Thus conventional absorption-line data, and it is clear that the
−
there appears to be essentially no color difference between twotechniquesagreeextremelywellwheretheyoverlap.The
the components. We also find no significant difference be- spheroiddispersionprofileisalsoverywell-behaved,showing
tweenthescale-lengthsofthecomponentsinthetwobands: littleifanyvariationwithradius(justifyingthecompromise
The origins of S0s from PN kinematics 9
Figure 7. GALFIT analysis of NGC 1023. The top row is for a R-band image; the bottom row is for an B-band image. The data are
obtainedwiththeIsaacNewtonTelescopeon1995December 25,usingtheprimefocuscamera.The300secondexposurescover afield
ofviewof10.2′×10.2′ (Noordermeeretal.2008).Thefirstcolumnshowstheoriginaldata,thesecondcolumnpresentsthemodel,and
thethirdcolumnshowstheresultofsubtractingthemodelfromthedata.Theresidualimagehasbeennormalizedbytheoriginaldata,
sothatthevaluesrangefrom0foraperfectfit,to1or−1,wherethedatahasnotbeenfitted atall.
in radial binning discussed in Section 2.2). It is also no- in this galaxy, it seems reasonable to reduce the number of
table that the approximately-constant spheroid dispersion freeparametersbyinvokingtheepicyclicapproximation.In
isconsistentwithavalueofV/√2,whichiswhatonewould particular, for such a cold disk system with a flat rotation
obtain for thesimplest possible modelof asingular isother- curve, we expect σφ/σr = 1/√2. With this additional con-
mal sphere potential, in which circular speed and velocity straint, we obtain thekinematic parameters plotted in Fig-
dispersion are related in this way. This consistency is reas- ure 10. The error bars on all parameters are duly reduced,
suring,asthefittingprocedureinnoway imposedit onthe and we now find that the kinematics of NGC 1023 are ex-
results. actly what one would expect for avery normal disk galaxy,
As one further enhancement, one might also expect with V/σφ 4.1 throughout the disk, similar to what one
≃
some rotational velocity, Vs, in the spheroidal component, finds in the stellar component of a spiral galaxy. Thus, we
which might affect the results. Indeed, there is known to seemtohavefoundtheexplanationforthestrangekinemat-
be a strong correlation between ellipticity and rotational icsoriginally inferred from thesedataforNGC1023, andit
speed in low-luminosity spheroidal systems like this bulge is now revealed to be most likely a spiral galaxy that has
(seeFigure4.14inBinney & Tremaine1987).TheGALFIT simply been stripped of its gas.
modeling shows that the spheroid in NGC 1023 has an el-
lipticity of 0.25, which translates into a predicted value
∼ 4.4 PNe objects rejections: the footprint of an
of Vs/σs 0.5. We have experimented with including rota-
∼ ongoing merger
tion at thislevelin thespheroidal component bymodifying
the first term in Equation 7, but find that at this level the However, the story does not quite finish there. Even with
inclusion of rotation makes no significant difference to the the full disk+spheroid kinematic model, 17 PNe are still
remaining kinematic parameters. likelihood-clippedatathresholdprobabilityof2.1σ.InFig-
Theonly slight disappointmentin thefittingprocess is ure 11, we show the locations of the PNe in both velocity
that theratio between thecomponentsof disk dispersion is and right ascension, with the rejected objects highlighted.
somewhatlesswelldefinedthanwasthecaseinthedisk-only SinceNGC1023liesataposition angleverycloseto90de-
modelfit,whichpresumablyreflectstheimpactoftheextra grees,thespatialcoordinateisessentiallythedistancealong
spheroidaldispersionfreeparameterinthismodel.However, themajoraxis,sothisisaconventionalposition–velocitydi-
nowthatweclearlyhaveawell-behavednormaldisksystem agram,withtheusualantisymmetricsignatureofarotating
10 Cortesi et al.
Figure 11.Plotofradial velocity versusrightascension forthe
PNe in NGC 1023. Filled symbols show the objects rejected by
thelikelihoodclippinginthefulldisk+spheroidmodel.Asterisks
showthePNeattributedtoNGC1023A.
almost all to form a continuousstream that passes through
this companion galaxy, as one might expect if the systems
are tidally interacting Capaccioli et al. (1986). It therefore
Figure 9. Derived mean rotation speed, the components of the
seems likely that these PNe lie in the tidal debris from this
disk velocity dispersion, and the spheroid dispersion versus ra-
companion as it is stripped in an on-going minor merger.
dius for NGC 1023. The filled symbols are from the maximum-
likelihoodanalysis,withverticalerrorbarsindicatinguncertainty
and horizontal error bars showing the extent of each radial bin
(with the point plotted at the median location for a PN in that
5 CONCLUSIONS
bin. The open symbols show rotation velocities derived from
absorption-linedatabyDebattista etal.(2002). Wehavepresentedanewmethodforanalyzingthekinemat-
icsofdiskgalaxiesasderivedfrom individualstellartracers
suchasPNe.Thishybridtechniquebinsdataradiallyinthe
galaxy to maintain the maximum flexibility in the inferred
kinematics, but uses a likelihood analysis within each bin
so as to derive the maximum amount of information from
the discrete data points. In addition, we use photometric
decomposition of continuum images to assign a probability
to each kinematic tracer of belonging to either the disk or
thespheroidal component of thegalaxy being modeled.
Applicationofthemethodtosimulateddatashowsthat
itreproducesmostoftheintrinsicdynamicsofagalaxyeven
when the number of discrete kinematic data points is rela-
tively modest. Application of this technique to NGC 1023
has offered an explanation for the strange kinematics pre-
viously inferred for this system. These peculiar properties
can be entirely attributed to uncorrected contamination by
the galaxy’s spheroidal component, and when this element
is properly modeled, the galaxy is revealed to have stellar
kinematics entirely consistent with a normal spiral system,
with a cold rotationally-supported disk and a hot spheroid
with the expected velocity dispersion. This result favours a
model in which NGC 1023 formed from a spiral via a sim-
Figure10.AsforFigure9,butwiththeratioofdiskdispersions ple gas stripping process or secular evolution, rather than
components constrainedbytheepicyclicapproximation. through a more disruptiveminor merger.
Onebenefitofthelikelihoodanalysismadeclearbythis
applicationisthatitispossibleto“likelihoodclip”,toiden-
diskevidentinthePNe.TherejectedPNe,however,donot tify objects that do not seem to be drawn from the model,
show suchantisymmetry,with thevast majority located on in a reasonably robust and objective manner. Comparison
one side of the galaxy. Such an asymmetric arrangement is of Figures 10 and 6shows that such clipping can deal quite
clearly notconsistent witherrorsarisingfrom apoorchoice wellwitheventheratherextremecasewhereonlythediskis
for any axisymmetric element in the model, or from mis- includedinthemodel.Whereamorerealisticdisk+spheroid
identified unrelated background objects. model is adopted, it seems to do a good job of identifying
A clue to their origin comes from considering theloca- featureslikestellarstreams,addingtoourunderstandingof
tion of NGC 1023A in this plot, as the rejected PNe seem thedynamical properties of such systems.