Astronomy&Astrophysicsmanuscriptno.GC_starcounts_Gallego (cid:13)cESO2017 January17,2017 The distribution of old stars around the Milky Way’s central black hole: I. Star counts E.Gallego-Cano1,R.Schödel1,H.Dong1,F.Nogueras-Lara1,A.T.Gallego-Calvente1,P.Amaro-Seoane2,andH. Baumgardt3 1 InstitutodeAstrofísicadeAndalucía(CSIC),GlorietadelaAstronomías/n,18008Granada,Spaine-mail:[email protected] 2 InstitutdeCiènciesdel’Espai(CSIC-IEEC)atCampusUAB,CarrerdeCanMagranss/n08193Barcelona,Spain InstituteofAppliedMathematics,AcademyofMathematicsandSystemsScience,ChineseAcademyofSciences,Beijing100190, China KavliInstituteforAstronomyandAstrophysics,Beijing100871,China 7 ZentrumfürAstronomieundAstrophysik,TUBerlin,Hardenbergstraße36,10623Berlin,Germany 1 3 SchoolofMathematicsandPhysics,UniversityofQueenslandSt.Lucia,QLD4068,Australia 0 2 Received;accepted n a ABSTRACT J 3 Context.Theexistenceofdynamicallyrelaxedstellardensitycuspsindenseclustersaroundmassiveblackholesisalong-standing 1 predictionofstellardynamics,butithassofarescapedunambiguousobservationalconfirmation. Aims.InthispaperwerevisittheproblemofinferringtheinnermoststructureoftheMilkyWay’snuclearstarclusterviastarcounts, ] toclarifywhetheritdisplaysacoreoracusparoundthecentralblackhole. A Methods.WeusejudiciouslyselectedadaptiveopticsassistedhighangularresolutionimagesobtainedwiththeNACOinstrumentat G theESOVLT.ThroughimagestackingandimprovedPSFfittingwepushthecompletenesslimitaboutonemagnitudedeeperthan inprevious,comparablework.Crowdingandextinctioncorrectionsarederivedandappliedtothesurfacedensityestimates.Known . h young,andthereforedynamicallynotrelaxedstars,areexcludedfromtheanalysis.Contrarytopreviouswork,weanalysethestellar p densityinwell-definedmagnituderangesinordertobeabletoconstrainstellarmassesandages. - Results. WefocusonstarsintheRedClump(RC),withobservedmagnitudes K ≈ 15.5,andonstarswithobservedmagnitudes o K ≈18,whichweexpecttohavesimilarmeanagesandmassesthantheformer.TheRCandbrightergiantstarsdisplayacore-like r t surfacedensityprofilewithinaprojectedradiusR<0.3pcofthecentralblackhole,inagreementwithpreviousstudies,butshowa s cusp-likesurfacedensitydistributionatlargerR.Thesurfacedensityofthefainterstarscanbedescribedwellbyasinglepower-law a atR < 2pc.Thecusp-likeprofileofthefaintstarspersistsevenifwetakeintoaccountthepossiblecontaminationofstarsinthis [ brightnessrangebyyoungpre-mainsequencestars.Thedataareinconsistentwithacore-profileforthefaintstars.Weexcludefrom 1 ouranalysistheregionR < 1”/0.04pc,wherethedensityislessthanexpectedatallbrightnessranges,becausethisregionmaybe v toocrowdedandalsohassufferedchangestoitsstellarpopulationbeyonddynamicalrelaxation.Finally,weshowthata3DNuker 6 lawprovidesaverygooddescriptionoftheclusterstructure. 1 Conclusions.WeconcludethattheobservedstellardensityattheGalacticCentre,asitcanbeinferredwithcurrentinstruments,is 8 consistentwiththeexistenceofastellarcusparoundtheMilkyWay’scentralblackhole,SagittariusA*.Thiscuspiswelldeveloped 3 inside the influence radius of about 3pc of SagittariusA* and can be described by a single three-dimensional power-law with an 0 exponent γ = 1.23±0.05. This corroborates existing conclusions from Nbody simulations performed in a companion paper. The . apparentlackofRCstarsandbrightergiantsatprojecteddistancesofR (cid:46) 0.3pc(R (cid:46) 8”)ofthemassiveblackholemayindicate 1 thatsomemechanismhasalteredtheirdistributionorintrinsicluminosity.Weestimatethenumberofpossiblymissinggiantstoa 0 fewdozen.ThisvalueagreeswellwiththetheoreticalideathatapreviouslyexistinggaseousfragmentingdiscaroundSagittariusA*, 7 precursoroftheobservedyoung,massivestars,wasresponsibleforrenderinggiantstarsinvisible. 1 : Keywords. Galaxy:center–Galaxy:kinematicsanddynamics–Galaxy:nucleus v i X r1. Introduction can test observationally the existence of a stellar cusp. Such a a stellar cusp is a prediction of stellar dynamics for the case of Thisisthefirstoneofthreepapersaddressingthedistributionof adynamicallyrelaxedclusterandhasbeenderivedandstudied stars around SagittariusA* (SgrA*), the massive black hole at by analytical, Monte Carlo and N-body methods (e.g., Bahcall thecentreoftheMilkyWay.Theyarecloselyrelated,butfocus & Wolf 1976; Lightman & Shapiro 1977; Murphy et al. 1991; ondifferentmethodsandstellarpopulations.Inthisworkweuse Baumgardt et al. 2004; Amaro-Seoane et al. 2004; Alexander the method of star counts, while in our other paper (Schödel et & Hopman 2009; Preto & Amaro-Seoane 2010). These consis- al.,referredtoasPaperIIinthefollowing),weanalysethedif- tenttheoreticalresultshave,however,notyetbeenconfirmedob- fuselightfromtheunresolvedstellarpopulation.Finally,Baum- servationally.Thereexistcurrentlyonlymeasurementsofabout gardtetal.(PaperIII)presentnewNbodysimulationsthatcon- twodozensystems,wherewecanactuallyresolvetheradiusof firmtheagreementbetweenmodellingandobservations. influence of the central black hole with several resolution ele- ThedistributionofstarsaroundSgrA*isofgreatastrophys- ments(see,e.g.,Table1inGültekinetal.2009).Thegreatdis- icalinterestbecauseitistheonlymassiveblackholewherewe Articlenumber,page1of14 A&Aproofs:manuscriptno.GC_starcounts_Gallego tanceofmostextragalacticsystemsmeansthatwecanonlystudy (Yusef-Zadehetal.2012).Hence,wehaveonlyobservedthetip the light density of hundreds to thousands, or even millions, of of the iceberg, which may not be representative for the overall, starsperresolutionelement.Sincenuclearstarclustersareenti- underlyingstructure. tieswithcomplexstellarpopulations,manyofwhichshowsigns This work continues similar efforts carried out by Genzel ofrecentstarformation(seereviewby Böker2010),measured et al. (2003) and Schödel et al. (2007). As experiments are re- lightdensitiesmayfrequentlybedominatedbyasmallnumber peated, both their accuracy and precision tend to increase be- ofbrightstars,whicharegenerallytooyoungtobedynamically cause of factors such as improvements in observational tech- relaxed.Thiscanleadtoambiguousorerroneousresults. niques, increasing know-how on data reduction and analysis, Ontheotherhand,withitsdistanceofonlyabout8kpcfrom progressintheoryandinterpretation,andanincreasingamount Earth,theGalacticCentre(GC)allowsustoresolvethestarsob- of data. The novel aspects of this work are, in particular, the servationallyonscalesofabout2milli-parsecs(mpc),assuming stackingofhighqualityimageswithalargefield-of-view(FOV) diffractionlimitedobservationsatabout2µmatan8-10mtele- toreachfaintercompletenesslimitsinthemostcrowdedregions scope.AtthecentreoftheMilkyWay,a4×106M (e.g.,Boehle near SgrA*, an extension of high angular resolution data to a (cid:12) et al. 2016; Gillessen et al. 2016) massive black hole, Sagittar- largerfieldofabout1.5(cid:48)×1.5(cid:48),improvementsindatareduction iusA*(SgrA*),liesembeddedina∼2.5×107M nuclearstar (rebinningoftheimages,removalofsystematicnoisefromde- (cid:12) cluster(Schödeletal.2014a,b).Therefore,theGCappearstobe, tector electronics), and analysis (improved PSF fitting with use inprinciple,theidealtestcasefortheexistenceofstellarcusps. ofnoisemapsandaspatiallyvariablePSF).Finally,moreexplic- Surprisingly,nounambiguousobservationalevidenceforthe itly than in previous work – and because our deeper data allow existence(ornot)ofastellarcusparoundSgrA*hasbeenpre- us to do so – we focus on clearly delimited stellar brightness sentedsofar.Whilethefirsthighangularresolutionobservations rangesinordertominimisethemixingofdifferentstellarpopu- at an 8m telescope appeared to indicate the existence of a stel- lations.Weaddnewdataonthestellardistributionatprojected lar cusp (Genzel et al. 2003; Schödel et al. 2007), it was later radii R (cid:38) 2pc from the literature (Schödel et al. 2014a; Fritz realised that the star counts within about 0.5pc of SgrA* were et al. 2016) to facilitate the interpretation of the data and the contaminated by a significant number of young, and therefore derivationofthe3Ddensitystructureofthestarsnearthemas- dynamicallyunrelaxed,stars.Whenomittingtheyoungstars,the siveblackhole. projectedstellardensityofgiantsappearsalmostflat,core-like, within a few 0.1pc of SgrA* (Buchholz et al. 2009; Do et al. 2. Datareductionandanalysis 2009; Bartko et al. 2010). Observations of the stellar surface brightness from old stars also appeared to indicate a possibly 2.1. Basicreduction core-like structure (Fritz et al. 2016). These findings led to the missing cusp problem and have given rise to a large number of We use H and Ks-band data obtained with the S27 camera theoreticalpapers,tryingtoexplainitsabsence.Thehypotheses (0.027”pixelscale)ofNACO/VLT.TheAOwaslockedonthe reach from a very long relaxation time (Merritt 2010), over the NIRbrightsupergiantGCIRS7thatislocatedabout5.5”north destructionoftheenvelopesofgiants–thusrenderingthemin- of SgrA*. The data used are summarised in Table1. All data visible–viastellarcollisions,whichcannotfullyexplaintheob- wereacquiredwithasimilar4pointditherpattern,roughlycen- servations(Daleetal.2009)orafragmentedgaseousdisc,which tered on SgrA*, with the exception of the data from 11 May can (Amaro-Seoane & Chen 2014). Other explanations involve 2011, which covered a shallow, but wider mosaic with a 4×4 theapparentstellarstructurearisingfromsubsequentepochsof dither pattern, centered on SgrA*. Preliminary data reduction starformationand/ortheaccumulationofstellarremnantsnear was standard, with sky subtraction, bad pixel removal and flat SgrA*(e.g.,Löckmannetal.2010;Aharon&Perets2015). fielding.Subsequently,asimpleshift-and-add(SSA)procedure When evaluating the observational evidence, it is, however, was applied to obtain final images. A quadratic interpolation ofutmostimportancetobeawareofitslimitations.Firstly,due with a rebinning factor of two was used because tests showed to the extreme interstellar extinction toward the Galactic Cen- thatthisimprovedthephotometryandreducedresidualsinPSF tre(e.g.,Nishiyamaetal.2009;Schödeletal.2010;Fritzetal. fitting,inparticularsincetheS27pixelscalebarelysamplesthe 2011), observations need to be performed in the near-infrared angularresolution,whichwasroughly0.06”FWHMforallim- (NIR).AsecondrequirementforobservingtheGalacticCentre ages. Along with the mean SSA images we also created noise (GC)istouseanangularresolutionashighaspossibletoover- mapsthatcontaintheerrorofthemeanofeachpixelintheSSA cometheextremesourcecrowding.Hereweusetheadaptiveop- images. tics(AO)assistedNIRcameraNACOinstalledattheESOVLT. ImagingdataatHandK areusedtobeabletoestimateextinc- S 2.2. Alignmentandstacking tionand–bythesamemeans–toexcludeforegroundstars. Because of these observational difficulties, our knowledge We treated each of the four pointings toward SgrA* indepen- about the stellar population at the GC is limited to the bright- dently to avoid problems arising from camera distortions near est few percent of stars: A few million-year-old hot post main the edges of the NACO S27 camera’s field (see Trippe et al. sequence (MS) giants and MS O/B stars (the latter being al- 2008; Schödel et al. 2009). The final images from all epochs readyatthefaintlimitofspectroscopiccapabilities),ontheone were aligned with the one from 09 Aug 2012 via a polynomial hand, and, on the other hand, giants with luminosities equal to fitofdegreeone(IDLPOLYWARPandPOLY_2Dprocedures). or higher than RC stars. In fact, the typical spectroscopic com- Theparametersofthelatterweredeterminedviaaniterativefit pleteness reaches only about K = 15.5 stars and thus only usinglistsofdetectedstarsintheimage.Theimageswerecom- S half the RC (see Do et al. 2009; Bartko et al. 2010; Do et al. binedinasimplemeanandthecorrespondingnoisemapswere 2015).StudiesofthestellarsurfacedensityattheGChavesofar quadraticallycombined.Apossibleconcerninthisstackingpro- beendominatedbyRCstarsandbrightergiants.Theonlyrecent cedure is the use of different observing epochs because of the study that focused on the light from stars fainter than this limit largepropermotionsofthestarsintheGC.Atthedistanceofthe didindeedfindacusp-likestructurewithin5”/0.2pcofSgrA* GC(hereassumedas8kpc,see,e.g.,Genzeletal.2010;Meyer Articlenumber,page2of14 E.Gallego-Cano etal.:ThedistributionofoldstarsaroundtheMilkyWay’scentralblackhole:I.Starcounts Table1.Detailsoftheimagingobservationsusedinthiswork. 2.3. Patternremoval Datea λ ∆λ Nb NDITc DITd The images from the individual epochs contained horizontal central [µm] [µm] [s] stripe patterns from the detector electronics. These horizontal 09May2010 1.66 0.33 4 64 2 stripescanbedetrimentalforsourcedetectionbecausetheymay 17May2011 2.18 0.35 4 9 2 either mask faint sources or be deblended into rows of stars 09Aug2012 2.18 0.35 8 60 1 by the PSF fitting program. It is therefore important to remove 11Sep2012 2.18 0.35 8 60 1 them. We proceeded as follows: We used the StarFinder pro- 12Sep2012 2.18 0.35 8 60 1 gramtodetectandsubtractrobustlydetectedpointsourcesfrom eachimage(conservativesettingsoftheStarFinderparameters: a UTCdateofbeginningofnight. min_correlation= 0.85 and deblend= 0) and to fit the diffuse emission (from unresolved stars or dust and gas in the inter- b Numberof(dithered)exposures stellarmedium).Thelatterwasfittedwithanangularresolution of about 0.25”, a non-critical value that just needs to be large c Number of integrations that were averaged on-line by the enough to remove the variable background due to unresolved read-outelectronics stellar emission and small enough to roughly correspond to the d Detectorintegrationtime.Thetotalintegrationtimeofeach sizeofdiffuse,uncresolvedstructuresinthemini-spiral(seePa- observationamountstoN×NDIT×DIT. perII). While fitting of the diffuse background is important in this procedure,the exact choiceof its variabilityscale is not. It caneasilybechosentobeafactor2largerorsmaller. Theresultingresidualimages,i.e.imageminusdiffuseemis- sion minus point sources, were then dominated by small-scale etal.2012),avelocityofabout40kms−1ontheplaneofthesky (ontheorderafewpixelswidth)randomandsystematicnoise. correspondstoapropermotionofonemilli-arcsecondperyear. We determined the pattern of horizontal stripes induced by the AdisplacementbyonepixeloftheNACOS27cameraperyear electronicsthroughmediansmoothingeachrowofpixelswitha thereforecorrespondstoavelocity>1000kms−1.Sinceweare medianboxwidthofabout2.7”,correspondingto200pixels(in notinterestedinhighprecisionastrometryorphotometry,thisef- therebinnedimages).Thispatternwasthensubtractedfromthe fectisthereforenegligibleforourdata,exceptpossiblyasmall SSAimages.Wecouldthusremovemostofthesystematicnoise numberofveryfastmovingstarswithin∼0.1”ofSgrA*,which withoutintroducinganysignificantbiasonthepointsourcesor isnotrelevanttotheproblemandscalesaddressedinthispaper. onthediffuseemissionbecausemoststarshadalreadybeensub- tractedandbecausethemediansmoothingboxwasafactorofa few to ten larger than the scales of the diffuse emission, of the August 2012 Deep mean image size of PSF residuals, or of faint, unresolved sources. Finally, afterhavingcleanedtheimagesofeachepoch,theywerecom- bined to a deep mean image. This last step further reduced any remaining systematics. To be conservative, we used the noise mapsderivedfromtheuncleanedimages.Figure1showsdetails of K -images to illustrate the effect of the systematic readout s noiseandtheimprovementafterremovingit. 2.4. Sourcedetection PointsourceextractionwascarriedoutwiththePSFfittingpro- 5" grammeStarFinder(Diolaitietal.2000).Sincetheimagescover areassimilartoorlargerthantheisoplanaticanglesattherespec- tivewavelengths,carewastakentodealwiththespatialvariabil- ity of the PSF. We parted the images into sub-fields of approx- imately 10.5” × 10.5” size. Subsequently, ten of the brightest, most isolated stars in each sub-field were used for an iterative extraction of a local PSF (similar to what was done in Schödel et al. 2010). Because of variable extinction and source density, notallsub-fieldscontainedPSFreferencestarsofsimilarbright- ness,whichwouldleadtoasystematicchangeofthezeropoint acrossthefield.Also,nottakingintoaccounttheextendedsee- inghalofromthelightthatcouldnotbecorrectedbytheAO,can leadtoanenhanceddetectionofspuriousfaintstarsnearbright Fig.1. Cleaningofhorizontalstripes(systematicreadoutnoise).Upper stars.AsremarkedbySchödel(2010),theseeinghaloisaffected left:DetailofAugust2012Ks-bandimage.Lowerleft:Asupperleft, inaratherminorwaybyanisoplanaticeffects.Wethereforeused butcleaned.Upperright:Detailofdeep,meancombinedKs-bandim- the brightest star in the field, GCIRS7, to estimate the seeing agewhentheinputimageshavenotbeencleaned.Lowerright:Detail halo. The local PSFs were masked beyond radii of about 0.3”, ofdeep,meancombinedKs-bandimageaftercleaningoftheinputim- ages.Thedisplayedfieldislocatedabout12.0”westand1.7”northof uptowhichtheycouldbereliablydetermined.Thentheywere SgrA*.Thecolourscaleislogarithmicandidenticalforallimages. matchedtotheseeinghalo(usingaleast-squaresfittodetermine fluxoffsetsandnormalisationfactors).Thus,wecouldcreatelo- calPSFsthatavoidedlargesystematicphotometriceffectsacross Articlenumber,page3of14 A&Aproofs:manuscriptno.GC_starcounts_Gallego thefield.Sometests(similartowhatwasdoneinSchödel2010) artificialstar,thentheartificialstarwasconsideredasdetected. showed that we could constrain the systematic photometric ef- This latter point is critical to avoid bias because the relatively fectsfromthevariabilityofthePSFtoafewpercentacrossthe highdensityofartificallyintroducedstarswouldotherwiselead field. to non-detection of real sources and thus an over-estimation of In PSF fitting we have to walk a thin line between achiev- incompleteness. ing an almost complete detection of sources while, at the same As mentioned above, to estimate the systematic errors in- time,avoidingtopickupspuriousones,whichcanarise,inpar- ducedbyeitherthenon-detectionofrealsourcesorthedetection ticular, close to bright stars or due to systematic effects from ofspurioussources,werepeatedthesourcedetectionandcom- thedetectorelectronics.WevisuallyverifiedthattakingthePSF pleteness determination for the following combinations of the seeing halos into account, along with the use of our SSA noise StarFinder parameters: min_correlation= 0.80,0.85,0.85,0.90 maps, effectively suppressed the detection of spurious sources and deblend= 0,0,1,1 (each value in the first list corresponds near bright stars (see also Schödel et al. 2013). The PSF halos to the value with the same index in the second list). In Fig.2, includeeffectssuchasdiffractionspikesandstaticspeckles.Our weshowthevaluesofcompletenessfortwoofthesecasesand empirical noise maps seemed to deal well with suppressing the for different projected distance ranges from SgrA*. The dif- detectionofspurioussourcesnearbrightstars. ferences between the different choices of parameters are gen- Finally, since there can be no absolute certainty in the reli- erally small, on the order of a few percentage points, except abilityofsourcedetection,wealsorepeatedlyanalysedtheim- for the faintest magnitudes, where the differences are some- ageswithdifferentvaluesoftheStarFinderparametersthatdom- whatmorepronounced.Also,wecanobservetheexpectedgen- inate the probability of source detection (for a fixed detection eral trend of less completeness for fainter magnitudes and in threshold, which was chosen as 3σ in all cases). These param- the more crowded areas near SgrA*. Finally, Fig.2 also shows etersaremin_correlationanddeblend.Fortheminimumcorre- thatthedifferencesofcompletenessbetweenthefourpointings lationvaluewechose0.8−0.9,alwaysmoreconservativethan are small. For all cases, we found that source detection was, at the standard value of 0.7, the default value of StarFinder. The all projected distances, at least 50% complete for magnitudes keyworddeblendcanbesettodeblendclosesources.Whilede- Ks ≤18.5. blendingcanbeveryuseful,itcanleadalsotothedetectionof asignificantnumberofspurioussourcesinacrowdedfield.We 2.7. Extinction includedmeasurementswithandwithoutsettingthiskeyword. Weusedthe H−K photometryandtheintrinsicsmallcolours s of stars at these bands to create an extinction map for the en- 2.5. Photometriccalibrationandsourceselection tire field, with the same method as applied in Schödel et al. (2010). We do not consider stars with H − K < 1.5 because Finally, the photometry was calibrated with the stars s we consider them as foreground stars. Neither do we consider IRS16C, IRS16NW, and IRS33N (apparent magnitudes starswithH−K >3.0becausetheymayeitherbebackground K = 9.93,10.14,11.20 and H = 11.90,12.03,13.24, see s s starsorintrinsicallyreddenedobjects(inanycase,theirnumber Schödel et al. 2010). The uncertainty of the zero points was a is very small, see Fig4 in Schödel et al. 2010). Median stellar few percent. We note that for the purposes in this paper we do colourswereobtainedfromtheindividualcoloursofthe20near- not require any high accuracy/high precision photometry and eststarsateachpositionandtheextinctionwasthencalculated astrometry. asinSchödeletal.(2010),assumingA ∝λ−2.2. Almost all stars in the field have intrinsic colours −0.1 ≤ Ks Ontheonehand,theextinctionmapwasusedtocorrectthe H−K ≤0.3(see,e.g.,Doetal.2009;Schödeletal.2010).The individual stellar magnitudes for differential extinction. On the meancolourduetoreddeningis H−K ≈ 2.1.Weexcludedall other hand, we applied the methodology of Chatzopoulos et al. stars with H − K < 1.5 as foreground stars. We also excluded (2015b)tocomputethestellardetectioncompletenessvariation spectroscopically identified young stars from our final star list causedbyvariableextinction:Wemodeledtheluminosityfunc- (Doetal.2009;Bartkoetal.2010).Subsequently,wecreatedan tion (LF) by taking the product of a power-law stellar LF and extinctionmap,byusingthe20starsnearesttoeachpoint.The anerrorfunctionthatrepresentsthecompletenessfunction,asin resulting map is similar, to within the uncertainties, to the one expression (2) in Chatzopoulos et al. (2015b). The approxima- presentedinSchödeletal.(2010). tion of the LF with a power-law – which ignores the presence of the RC bump – does not introduce any significant error be- causeourdataaresensitiveenoughtoreachwellbelowtheRC 2.6. Crowdingandcompleteness bump over the entire field and because including the RC bump We determined the source detection completeness in the K - wouldonlyhaveaminoreffectasshownbyChatzopoulosetal. s imagesthroughthetechniqueofinsertingandrecoveringartifi- (2015b).First,wemeasuredtheobservedKLFforeachpointing cialstars.Weusedamagnitudestepof0.5magandinsertedthe and each StarFinder parameter set (excluding a radius of about starsona0.5”×0.5”grid.Withthisrelativelywidespacingwe 5”aroundSgrA*,wheretheKLFismoreincompletebecauseof avoidedartificiallyincreasingthecrowding.Thegridwasshifted crowding). We computed the power law for each case. Finally, severaltimestofinallyprobecompletenessonadense0.1”×0.1” we used these power-laws, combined it with the measured lo- grid (as done by Schödel et al. 2007). We used the respective cal extinction and computed the corresponding local correction local PSFs (see above). Subsequently, PSF fitting was carried factorsaccordingtoequation(5)inChatzopoulosetal.(2015b), out with StarFinder and a source was considered as detected if but using the approximation of a single extinction screen, i.e. it was found within a magnitude range of 0.5mag of the input no variability of A along the line-of-sight. Since we approxi- Ks magnitudeandwithinadistanceof0.054”oftheinputposition mate the LF with a power law, the equation takes on the form (correspondingto2pixelsoftheS27cameraorroughlythean- p= L(−∆Ak)=10−γ∗∆Ak (seeChatzopoulosetal.2015b),where gularresolutionofthedata).Ifarealstarofasimilarmagnitude pisthereductionfactorforthenumberoflocallydetectedstars, was already present within this distance to the grid point of an ∆A isthedifferencebetweenthelocalextinctionandthemean k Articlenumber,page4of14 E.Gallego-Cano etal.:ThedistributionofoldstarsaroundtheMilkyWay’scentralblackhole:I.Starcounts a) b) c) Fig. 2. Completeness of the star counts in the deep NACO K images. a) Completeness in pointing 1 for different projected distance ranges s fromSgrA*,formin_correlation= 0.80anddeblend= 0.b)Asina),butformin_correlation= 0.90anddeblend= 1.c)Completenessforall four pointings and within 5” of SgrA*, for min_correlation= 0.80 and deblend= 0. The corresponding plots for other used combinations of min_correlationanddeblendlookverysimilar. extinctionoverthefield,andγisthepowerlawindexofthelu- minosityfunction.Ifthelocalextinctionislowerthanthemean extinction,then p > 1,andifthelocalextinctionishigherthan themeanextinction,then p<1. We apply the correction factor 1/p to each detected star. In Fig.3weshowthe K LFforpointing1,formin_correlation= s 0.8 and deblend= 0, along with the LF corrected for differen- tial extinction, a power law fit to the stars 11 < K < 14.5 ( s γ = 0.26±0.02,similar to the value obtained in Schödel et al. (2010)),andthecompletenessfunction(blueline),asdefinedby Chatzopoulos et al. (2015b). For the latter, we use m = 19.4 0 andσ = 0.2becausethesevaluesapproximateour K LFwell. s ThesevaluesaredifferentfromthoseusedinChatzopoulosetal. (2015b) because our data are significantly deeper than theirs. In Fig.4 the percentage reduction in observed stars versus pro- jectedradiusisrepresentedforthedetectedstarsinpointing1,for min_correlation= 0.8anddeblend= 0.Onecanseethatextinc- tion is, on average, higher at larger distances from SgrA*, but thattheeffectofdifferentialextinctiononcompletenessisrela- tivelyminor,typically<10%andatmost20%. Fig.3. TheKsLFinpointing1formin_correlation=0.8anddeblend= 2.8. Widefield 0ifshownasablackline.ThedottedredlineistheKsLFaftercorrect- ingthemagnitudeofeachstarfordifferentialextinction.Thegreenlines showsapowerlawfittothebrightstars11<K <14.5withpowerlaw s indexof0.26±0.02.Thebluelineshowstheeffectofthecompleteness The2011dataareofexcellentquality,butrelativelyshallow.On functionaccordingtoequation(2)inChatzopoulosetal.(2015b). theotherhand,thereare16pointingsthatcoverafieldofabout 1.5(cid:48)×1.5(cid:48),comparedtothesmallerfieldsofabout40”×40”cov- ered by the other NACO observations. The 2011 observations 3. ThesurfacedensityofoldstarsintheGC arethereforeideallysuitedtoextendthesensitivecentralobser- vationswithhighangularresolution,albeitsomewhatshallower, 3.1. K luminosityfunction s dataouttolargerdistances.Fig.5.Inordertodealwiththedis- tortionsoftheNACOS27camera(see,e.g.,Trippeetal.2008; WeshowtheK LF(KLF)determinedfromourdeepmosaicin s Schödeletal.2009)wealignedeachpointingoftheNACOmo- Fig.6.TheKLFscorrespondingtothefourdifferentStarFinder saic with a reference frame created from positions measured in parameter settings are shown. The KLF derived in this work is HSTWFC3observationsofthesamefield(Dongetal.,inprep.). aboutonemagnitudedeepercomparedtopreviouswork(Fig.10 WeapplyvariablePSFfittingasweexplaininsection2.4.Inthis in Schödel et al. 2007). This is a decisive advantage. When we case,min_correlation=0.8anddeblend=0wereselectedforthe wanttoprobetheexistenceofadynamicallyrelaxedstellarcusp StarFinderparameters. around SgrA*, we need to focus on stars that are at least sev- Articlenumber,page5of14 A&Aproofs:manuscriptno.GC_starcounts_Gallego other, fainter magnitude range with a high fraction of old stars (17.5≤ K ≤18.5).Also,thestarsinthesetwobrightnessranges s havesimilarmasses,whichletsusexpectasimilarsurfaceden- sitydistribution. Thegoalofthisstudyistoinvestigatetheexistenceofastel- lar cusp at the GC. This requires us to select stars old enough tohaveundergonedynamicalrelaxation.Therelaxationtimeat the GC is roughly a few Gyr (e.g., Alexander 2005, 2011). We specifically exclude all spectroscopically identified early-type, i.e. young and massive, stars from our sample (using the data of Do et al. 2013). Unfortunately, spectroscopic stellar classi- fication is limited to stars of about K ≤ 15.5 at the GC. For s fainter stars, we can only use their luminosity as a proxy for their type. Fig.16 of Schödel et al. (2007) illustrates the LF, mean masses, and old star fractions for stars of different mag- nitudes at the GC, assuming a basic continuous star formation model.Itshowsthatwecanprobeold((cid:38)1Gyr),low-massstars in the range 15 (cid:46) K ≤ 16. This is the RC, which dominates s allpreviousstardensitymeasurements.Thefractionofoldstars risesagainabove∼50%forstarsK >17.5,reachingpractically s 100%atK ≈18.Therefore,inthisstudy,wefocusonthemag- s nituderanges15 ≤ K ≤ 16,theRC,and17.5 ≤ K ≤ 18.5,i.e. s s thefainteststarsaccessiblebyourdata. Fig. 4. Relative detection frequency due to extinction versus pro- jecteddistancetoSgrA*,inpointing1formin_correlation= 0.8and deblend= 0.Thegreenlinerepresentsthemeanofthe p(%)consider- ingdetectedstarsatthesamedistancefromSgrA*.Wecanobservethat forclosedistancestoSgrA*p(%)ishigherthanforlargedistances,as weexpected,becausetheextinctionnearSgA*islower. Fig. 6. KLF for the deep K mosaic. The different colours corre- s spond to the different combinations of the values of min_correlation N anddeblend. 1 parsec E 3.2. Radialprofileofthenumberdensity Fig.5. Wide-fieldmosaicoftheobservationsfrom11thMay2011.The field-of-viewis1.5(cid:48)×1.5(cid:48).Thefieldofabout40”×40”thatcorresponds In order to analyse the surface density of stars at the Galactic tothedeepimagingdataismarkedbyawhitesquare. Center,weassumethattheunderlyingspatialdistributionofthe starsinthecentralparsecissphericallysymmetric.Althoughthe nuclearclusterisflattened,asphericalapproximationshouldbe eral Gyr old. As shown in the illustrative Figure16 in Schödel acceptable at projected radii R (cid:46) 2pc because the difference et al. (2007), the only magnitude range where this was previ- between the density profiles along the orthogonal directions of ouslypossiblewasaroundtheRC(15.25 ≤ K ≤ 16.25).Now, maximumdifferenceisonlyontheorderof10%−20%inthis s with the deeper data from our new analysis, we can probe an- region(seeSchödeletal.2014a;Fritzetal.2016). Articlenumber,page6of14 E.Gallego-Cano etal.:ThedistributionofoldstarsaroundtheMilkyWay’scentralblackhole:I.Starcounts a) b) c) Fig.8. Azimuthallyaveragedextinctionandcrowdingcorrectedsurfacedensitiesvs.distancefromSgrA*.Straightdottedlinesillustratefitswith asinglepower-law.a)Surfacedensityformin_correlation=0.85anddeblend=1.Thevalueofthepower-lawindex,Γ,forRCstars(redline)is Γ=−0.36±0.04andforfainterstars(blackline)isΓ=−0.52±0.02.b)Surfacedensityformin_correlation=0.90anddeblend=1.Thevalue oftheindexforRCstars(redline)isΓ=−0.34±0.04andforfainterstars(blackline)isΓ=−0.44±0.02.c)Surfacedensityprofilesforstarsof 17.5≤K ≤18.5forthedataobtainedwiththedifferentStarFinderparameters. s runs. For faint stellar magnitudes, we masked the regions with completeness below 30%. We tested the effect of different masks, with completeness< 30%,40%,50% respectively, and found that the results did not vary significantly. Fig.8 shows the radial surface density profiles for RC stars (K ≈ 15.5) s and fainter stars (K ≈ 18) for two different combinations of s StarFinder parameters (left panel: min_correlation= 0.85 and deblend= 1;middlepanel:min_correlation= 0.9anddeblend= 1).Ascanbeseen,theresultingsurfacedensitiesareverysim- ilar.TherightpanelofFig.8showsthesurfacenumberdensity profilesforallfourdifferenceStarFinderparametersettings.Al- though there are offsets between the densities, their shapes are verysimilar. Finally, we combined the surface density profiles obtained withthefoursettingsoftheStarFinderparameters.Meandensi- ties and standard deviations were computed and all uncertain- ties were quadratically combined. The resulting final number density profiles for the magnitude ranges 15 ≤ K ≤ 16 and s 17.5 ≤ K ≤ 18.5 are shown in Fig.9. Simple power laws s were fit to the data at 1” ≤ R ≤ 20” (0.04pc ≤ R ≤ 0.8pc), i.e. we deliberately left out the innermost data points from the most crowded region around SgrA*. The corresponding power law indices are Γ = −0.45 ± 0.03 for the faint stars and cusp Γ = −0.34 ± 0.04 for RC stars. There is a significant dip Fig.7.Optimalbinning.WeshowtheRLPasafunctionofthenumber cusp aroundR=0.2pc(5(cid:48)(cid:48))inthedensityprofileoftheRCstarsand, ofbins.ThemaximumfortheRLPisreachedfor21bins(reddotted possibly, an excess at R ≈ 7”. It can also be seen that a simple line). power law is a better approximation for the faint stars than for theRCstars. We computed the azimuthally averaged stellar surface den- sitiesinannuliaroundSgrA*.Itisimportanttochooseanum- 3.3. Starcountsbeyond20” berofbinssufficientlylargetocapturethemajorfeaturesinthe datawhileignoringfinedetailsduetorandomfluctuations.Fol- Inordertostudythestellarnumberdensityinabroaderrangeof lowingthestudiesofKnuth(2006)andWitzeletal.(2012)we distancesfromSgrA*weanalysedthelargemosaicimagefrom firstdeterminethebestbinsize.ThedependenceoftheRelative the 2011 data. We did not apply any extinction and complete- LogarithmicPosteriorProbability(RLP)onthebinnumberfor nesscorrectionsbecausetheeffectoftheextinctioncorrectionis pointing1isshowninFig.7.ThemaximumfortheRLPforthe smallandcrowdingdoenotposeanyseriousproblematR>20” starnumberisreachedfor21bins,andthebestbinsizeis1(cid:48)(cid:48).We and with the high angular resolution data used here. We did, appliedthismethodologyforallpointings,withsimilarresults. however mask all the regions occupied by the dark clouds that We computed extinction and completeness-corrected stellar can be seen show in Fig.5. Finally, the surface densities from surfacedensitiesforthestarsdetectedinthedifferentStarFinder thewidefielddatawerescaledtothecompletenessandextinc- Articlenumber,page7of14 A&Aproofs:manuscriptno.GC_starcounts_Gallego Faint stars RC stars Fig.10. Combineddeepfield(blue)pluswidefield(green)surfacedensityplots. 3.4. Theprojectedstellarsurfacedensity As Figs.9 and 10 show, the projected surface number density of the faintest stars in our sample (17.5 ≤ K ≤ 18.5) can be s describedwell(reducedχ2 = 2.8)byasinglepower-lawofthe formΣ(R) ∝ R−Γ,whereΣisthesurfacenumberdensity,Rthe projecteddistancefromSgrA*,andΓthepower-lawindex.We findΓ=0.51±0.02forthefitshowninFig.10. TheRCstarspresentadifferentpicture:Asinglepower-law providesonlyagoodfittothedataatR (cid:38) 6”/0.24pc(reduced χ2 = 3.9), and with a higher value of Γ = 0.68±0.03 than in caseofthefaintstars(seeFig.10).Thisvalueagreeswellwith what has been found by previous authors for RC and brighter stars(seeSchödeletal.2007,andreferencestherein).Also,the flatteningoftheRCsurfacenumberdensityinsideR≈8”agrees verywellwithpreviousstudies(Buchholzetal.2009;Doetal. 2009;Bartkoetal.2010).Itisinterestingtonote,however,that Γ=0.52±0.03ifweforceasingle-powerlawfit(reducedχ2 = 7.5) to the RC number density profile, in agreement with the valuefoundforthefaintstars. 4. Discussion 4.1. Influenceofthecorrectionfactors Fig.9.Azimuthallyaveragedextinctionandcrowdingcorrectedsurface densityvs.distancefromSgrA*forallStarFinderdatasetsandforthe Whendeterminingthenumberdensityprofile,wehavetoapply Ks-magnituderanges15−16(redline)and17.5−18.5(blackline).Blue correction factors to compensate the effects of variable stellar linesrepresentsimplepower-lawfits,withthecorrespondingpower-law crowdingandinterstellarextinction.Fig.11showsthemeasured indicesindicatednexttothelines. surfacedensityprofileforstarsofmagnitude17.5 ≤ K ≤ 18.5 s without any correction applied (blue), after applying the com- pleteness correction for crowding (red), and after applying the tion corrected ones from the deep images in the overlap region completeness corrections for crowding and extinction (black). fromR = 10(cid:48)(cid:48) −20(cid:48)(cid:48).Fig.10showsthecombinednumberden- As we can see, the completeness correction steepens the pro- sityplotsforstarsinthemagnituderanges16.5≤ K ≤17.5and file somewhat. The extinction correction only introduces minor s (17.5 ≤ K ≤ 18.5(left)andRCstars(right),alongwithsimple changes because the azimuthal averaging compensates most of s power-lawfits. theeffectsofdifferentialextinctionacrossthefield.Insummary, Articlenumber,page8of14 E.Gallego-Cano etal.:ThedistributionofoldstarsaroundtheMilkyWay’scentralblackhole:I.Starcounts cending branch or main sequence (MS) stars of ∼2.5M . They (cid:12) couldalsobepre-MSstarsofafewsolarmassesorless(Luetal. 2013). From what is known about the star formation history of theNSCweexpectthatthemajorityofstarsisold(∼80%ofthe NSC’smasswereformed> 5Gyrago,accordingtoPfuhletal. 2011)andthatmostofthefaintstarsinoursamplearethusold, (sub-)giants.However,sinceweknowthatastarformationevent created on the order 104M of young stars in the region about (cid:12) 0.5pcaroundSgrA*,contaminationbyasignificantamountof pre-MSobjectsispossible.Wewilldiscussthispossibilityinthe nextsubsection. 4.3. Possiblecontaminationbypre-MSstars Fig.11. Meansurfacedensityprofileforstarswith(17.5≤K ≤18.5), s afteraveragingoverthefourrunswithdifferentStarfinderparameters. Blue:Uncorrecteddata.Red:Datacorrectedforcrowding.Black:Data correctedforcrowdingandextinction. the applied correction factors, albeit necessary, do not signifi- cantlyalterourresults.Thisshowsthatourdataarerobust. Table2.Parametersusedintheestimationofthesurfacedensityprofile ofpre-MSstarsandresultingcorrectedΓforthedensityprofileofstars withmagnitudes17.5≤K ≤18.5. s ID ηa Nb Γc 1d 1.40 146 -0.41±0.03 2e 0.93 200 -0.41±0.03 3f 1.10 200 -0.39±0.03 4g 0.93 300 -0.37±0.03 Fig.12. Measuredextinctionandcrowding-correctedsurfacenumber 5h 1.10 300 -0.36±0.03 densityprofilesfor16.5 ≤ K ≤ 17.5(redline)and17.5 ≤ K ≤ 18.5 s s (green line). The black line represents the surface density profile for Notes. 17.5≤K ≤18.5aftercorrectionforthepossiblecontaminationbypre- s MSstars,usingη=0.93±0.09(Doetal.2013)andN =300,roughly a Power-lawindexofthesurface-densityprofile. estimatedbyusingtheinformationinFig.15ofLuetal.(2013).The dottedlinesaresimplepower-lawfitstothedataandtheircorrespond- b Totalnumberofpre-MSstarsatR=0.8”−12.5”. ingindicesareindicatedinthepanel. c Power-lawindexofthesurface-densityprofileoffaintstars (17.5≤ K ≤18.5)aftercorrectionforpre-MSstars, Our primary goal in this study is to determine the spatial s distributionoftheold,relaxedstellarpopulationattheGC.Care d Valueoftheη-parameterfrom(Bartkoetal.2010). mustthereforebetakentoexcludeyoungandthereforeprobably dynamicallyunrelaxedstars.Wehaveexcludedfromtheanaly- e Valueoftheη-parameterfrom(Doetal.2013). sisallspectroscopicallyidentifiedearly-typestarsfromDoetal. f,h Intermediatevalueoftheη-parametertotest. (2013). Unfortunately, the limit of spectroscopic identification of early-type stars with current instruments is K ≈ 16 in the S g Valueoftheη-parameterfrom(Doetal.2013). GC (with the exception of a few, very deep exposures of small fieldsthatreached K ≈ 17.5,seePfuhletal.2011).Asargued s intheprevioussection,thestarformationeventthattookplace afewMyragointhecentralR ≈ 0.5pcmeansthatwehaveto considerexplicitlythepossibilitythatourstellarsurfacenumber 4.2. Tracerpopulations densitiesarecontaminatedbypre-MSstars. With an approximate magnitude of K = 18, the faintest stars To do this, we assume that the density profile of the young s in our sample are consistent with being (sub-)giants on the as- starscanbedescribedbyasimplepower-law(e.g.,Bartkoetal. Articlenumber,page9of14 A&Aproofs:manuscriptno.GC_starcounts_Gallego 2010;Luetal.2013).Themainparameterstotakeintoaccount the NSC stellar population is well mixed and that its average are: η, the index of the power law and N, the total number of properties(massfunction)donotchange.Wescaledthedataof young stars in the magnitude range 17.5 ≤ K ≤ 18.5 and Fritzetal.(2016)tooursintherange0.5pc≤R≤0.9pctoour s within R ≤ 12”. We assumed different values of η and N from surfacenumberdensityforstars17.5≤ K ≤18.5,correctedfor s the literature, subtracted the resulting pre-MS number density thepossiblepresenceofpre-MSstarswiththeparametersfrom profiles and then derived the power-law indices of the profiles model ID4 in Table2. Subsequently, we subtracted a constant forthecorrectednumberdensityprofiles,assummarisedinTa- star density offset to take into account the contribution of the ble2.InFig.12weshowthesurfacedensityprofileforstarsof starsthatdonotbelongtothenuclearclusterbuttootherstruc- 17.5 ≤ K ≤ 18.5bothwithandwithoutcorrectionofcontami- tures,suchasthenucleardiskortheGalacticBulge.Thissimple s nationbypre-MSstars.Ascanbeseenfromthedifferentvalues procedureispossiblebecausethescalelengthsofthelattercom- in Table2 , taking the possible bias from pre-MS stars into ac- ponents are one to several orders of magnitude larger than the countwitharangeofreasonablevalues,canreducethemeasured half-light radius of the NSC (see, e.g., Launhardt et al. 2002; power-law index from Γ = −0.51±0.02 to Γ = −0.36±0.03. Bland-Hawthorn&Gerhard2016). This value is the most conservative one from our estimates be- We used a 3D Nuker model (Lauer et al. 1995) as given in causeitassumesthehighestcontaminationratebypre-MSstars, equation 1 of Fritz et al. (2016) and projected it onto the sky and we will adopt it here as a measure for the projected sur- to fit the surface density profiles at R ≤ 20pc. The fits were face density of old, faint stars. We note that in PaperII we find repeated several times, using different values for the subtracted Γ = −0.28±0.03 for the surface density of even fainter, constantstardensityandfortheinnermostradialdatarangethat diffuse unresolvedstars.Thetwovaluesoverlapwithintheir2σuncer- wasincludedinthefit.Wealsoperformedafitwithoutthecor- tainties. rectionforthepotentialpre-MSstars.Theresultsofourfitsare The value of Γ = −0.36 ± 0.03 lies > 10σ away from a listed in Table3. While the parameters r (break radius) and β b flatcore.Finally,wealsonotethatthebestpower-lawfittothe (outer power-law index) vary significantly, γ (inner power-law surface number density of stars in the range 16.5 ≤ K ≤ 17.5 index) shows consistent values for all fits, which is reassuring s resultsinaverysimilarΓ=0.39±0.03.Althoughthereissome becauseitisthemostimportantparameterforthiswork,where uncertaintyastowhetherstarsinthisbrightnessrangearesuit- we are interested in whether there exists of stellar cusp at the abletracersofarelaxedpopulation,itisreassuringtonotethat GC. Fit IDs 1 and 2 explore extreme values for the subtracted thisvaluecoincideswiththepre-MScorrectedΓofthefaintest constantstardensity.Weusethemeananderrorofthemeanof tracers.Sincethemagnituderange16.5≤ K ≤17.5wouldcor- theparametersoftheotherfourfitstoobtainorientativevalues s respondtomoremassive–andthisshorter-lived–pre-MSstars, fortheaveragebestNukermodel:r = 3.0±0.4pc(77±10”), b weexpectaweakerpre-MScontaminationinthisrange(seealso γ =1.29±0.02,β=2.1±0.1,andadensityatthebreakradius Fig.5in Luetal.2013). of ρ(r ) = 3900±900pc−3. The best fit according to line 4 in b Table3isshowninFig.13. TheNukerfitshowsthatthefaintstarsshowacusp-likedis- 4.4. 3Dprofile:Nukerfit tribution around SgrA*. A flat core can be excluded with high significance.Weexplicitlynotethatinthefitspresentedinthis workweomittheregionR≤1”(0.04pc).Inthisregion,thestar counts appear to drop slightly below the expected levels. How- ever,thisregionisalsothemostcrowdedregion,whichmaylead tostrongsystematicsinthestarcounts.Additionally,thestellar population in the extremely close environment of SgrA* may havebeenaltered,asisindicatedbythepresenceofthesocalled ”S-stars”, apparently B-type MS stars that appear concentrated withinR < 1”ofSgrA*andmayhavebeendepositedthereby individual scatter or capture events (see, e.g., Eisenhauer et al. 2005;Genzeletal.2010;Alexander2011). 4.5. DistributionofgiantstarsnearSgrA* In agreement with previous work, we have found an unexpect- edlyflatsurfacedensityforRCstarsandbrightergiantswithin about 0.3pc of SgrA*. This may indicate a possible deficit of giants in this region. Here we produce a 3D Nuker model fit forthegiantstarsandtrytoconstrainthenumberofpotentially Fig. 13. Nuker model fit (red line) to the surface density profile for missinggiants.Weproceededasinsection4.4.Atlargeradii,we 17.5≤K ≤18.5correctedforthecontaminationbypre-mainsequence used the data from Fritz et al. (2016) and subtracted a constant s starsintheregion0.8”-12.5”. densityof0.07starspersquarearcsecondfromthedata.Wepro- ducedfitsforRCstars,usingasaproxyfortheselectionofRC In order to convert the measured 2D profile into a 3D den- starsthattheRCbumpintheK LFiscomprisedapproximately s sitylaw,weneedtodealwithprojectioneffects,whichrequires between15≤ K ≤16,andforspectroscopicallyidentifiedgiant s ustoconstrainthesurfacedensityonscaleslargerthanwhatwe stars(K ≤ 15.5)fordifferentrangesofR.Theresultingbest-fit s could measure with NACO. For this purpose we use the data parametersaregiveninTable4. fromFritzetal.(2016),whichtheyacquiredfromobservations As can be seen from the reduced χ2 values, the quality of with NACO/VLT, WFC3/HST, and VIRCAM/VISTA. To com- thefitimprovesasweomitthecentralmostdatapoints.Overall, binethesedatawithours,wehavetoassumethat,onlargescales, the surface densities of RC stars and spectroscopically identi- Articlenumber,page10of14