BALRoGO: Bayesian Astrometric Likelihood Recovery of Galactic Objects -- Global properties of over one hundred globular clusters with Gaia EDR3
MMNRAS , 1–13 (2021) Preprint 10 February 2021 Compiled using MNRAS L A TEX style file v3.0
BALRoGO: Bayesian Astrometric Likelihood Recovery of Galactic Objects– Global properties of over one hundred globular clusters with
Gaia
EDR3
Eduardo Vitral ★ Institut d’Astrophysique de Paris (UMR 7095: CNRS & Sorbonne Université), 98 bis Bd Arago, F-75014 Paris, France
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present BALRoGO: Bayesian Astrometric Likelihood Recovery of Galactic Objects, a public code to measure the centers,effective radii, and bulk proper motions of Milky Way globular clusters and Local Group dwarf spheroidals, whose data aremixed with Milky Way field stars. Our approach presents innovative methods such as surface density fits allowing for stronginterloper contamination and proper motion fits using a Pearson VII distribution for interlopers, instead of classic Gaussian-mixture recipes. We also use non-parametric approaches to represent the color-magnitude diagram of such stellar systems basedin their membership probabilities, previously derived from surface density and proper motion fits. The robustness of our methodis verified by comparing its results with previous estimates from the literature as well as by testing it on mock data from 𝑁 − bodysimulations. We applied BALRoGO to Gaia EDR3 data for over one hundred Milky Way globular clusters and nine Local Groupdwarf spheroidals, and we provide positions, effective radii, and bulk proper motions. Finally, we make our algorithm availableas an open source software. Key words: methods: data analysis – astrometry – proper motions – (Galaxy:) globular clusters: general – galaxies: dwarf –stars: kinematics and dynamics
The galactic neighborhood is a replete garden of important objects,such as globular clusters and dwarf spheroidal galaxies, that helpto explain the key ingredients of galactic evolution and other un-solved mysteries in astrophysics. Globular clusters (GCs) are amongthe oldest relics of the Universe, sometimes reaching ages close to13 Gyr (Marín-Franch et al. 2009), which are spherically shaped,compact collections of stars that orbit the Milky Way (MW), andwhose origin is still uncertain (Peebles & Dicke 1968; Peebles 1984;Searle & Zinn 1978; Bullock & Johnston 2005; Abadi, Navarro &Steinmetz 2006; Peñarrubia, Walker & Gilmore 2009). They haverecently been discovered to harbor multiple stellar populations (e.g.Carretta et al. 2009) and are frequently served as test laboratories forintermediate-mass black holes search (e.g. van der Marel & Ander-son 2010; Noyola, Gebhardt & Bergmann 2008, and more recentlyVitral & Mamon 2021). Dwarf spheroidal galaxies (dSphs), as well,are composed of old stellar systems and are located farther away thanGCs, but still in the Local Group. They are much fainter group ofstars and are often studied as tracers of dark matter (e.g., Battaglia,Helmi & Breddels 2013, Walker 2013 and Boldrini et al. 2020).The study of both GCs and dSphs made a major leap as data fromthe Gaia astrometric mission became available, specially with itssecond and early third release (hereafter Gaia DR2 and Gaia EDR3,respectively). This mission provided an overall astrometric coverageof stellar velocities, positions and magnitudes of more than 10 stars ★ E-mail: [email protected] in the MW and beyond and therefore allowed many deeper studiesof spherical systems such as GCs and dSphs, which can be wellseparated from MW field stars (i.e., interlopers) when combiningboth position, distance, and most importantly, proper motions. Newbulk proper motions were derived in Gaia Collaboration et al. (2018),Baumgardt et al. (2019) and Vasiliev (2019) with Gaia DR2 for morethan 150 GCs, and different methods were used to do so, althoughthey often relied on Gaussian mixtures to differentiate tracers frominterlopers. Recently, by analyzing the GC NGC 6397, Vitral &Mamon (2021) showed that the separation between GC stars andinterlopers was much more robust by assigning distribution functionsto the proper motion space which were not based on a Gaussianmixture, as usual, but rather considering a GC Gaussian componentplus a Pearson VII (Pearson 1916) distribution for the interlopers.In fact, with the arrival of Gaia data, extracting GC members,as well as dSphs became a very important topic, for which a robustapproach is required in order not to let interloper contamination biasone’s results. For example, if the selection of GC stars at its outskirtsis not well made, one could misinterpret the presence of tidal tails,and thus derive wrong conclusions on the cluster evolution.Given the importance of a separation algorithm between suchgalactic objects and interlopers, Bustos Fierro & Calderón (2019)proposed a method based on Gaussian mixtures and clustering algo-rithms, with only one step not performed by a machine and applied itto eight GCs. Their work joined previous attempts to separate clustermembers from interlopers through parametric and non-parametricmethods (e.g., Vasilevskis et al. 1965 and Galadi-Enriquez et al.1998, respectively). In this paper, we also present an algorithm to ex-tract GCs and dSphs members embedded in MW interlopers, which © a r X i v : . [ a s t r o - ph . GA ] F e b Eduardo Vitral employs a combination of Bayesian fits, giving the membership prob-ability of each star in the field, and non-parametric cuts on the color-magnitude diagram. Our method takes as an input a cone searchdelivered by Gaia and, from that data, estimates the object centerin right ascension and declination (i.e., 𝛼 and 𝛿 ), the bulk propermotions, 𝜇 𝛼, ∗ (i.e., [ d 𝛼 / d 𝑡 ] cos 𝛿 ) and 𝜇 𝛿 , and the effective radius(half projected number radius) by fitting discrete data. To validate ourmethod, we compare bulk proper motions provided in the literaturefor 13 galactic GCs using Gaia DR2, and we compare our effec-tive radii estimates for over a hundred sources with Gaia EDR3. Weuse Gaia EDR3 data from the two GCs NGC 6752 and NGC 6205(M 13) along with the Sculptor dSph galaxy as illustrative examplesin our Figures, throughout the paper.In addition, we derive effective radii, bulk proper motions andcenters for over a hundred GCs from the New Galactic Catalog(NGC) using Gaia EDR3 and make it available in Table format. Ouralgorithm developed with Python is also provided as an open sourcecode named BALRoGO: Bayesian Astrometric Likelihood Re-cover of Galactic Objects , and allows the user to adjust theinputs of our algorithm to better suit any particular source. Our algorithm is specially designed to deal with data from the Gaiaastrometric mission, and therefore we begin by acquiring all thestars in a two degrees cone search for each of the GCs and dSphsin Table C1. The first two steps of our approach deal with all thestars in the two degrees field of view, and use only positional ( 𝛼 , 𝛿 )information. As it is going to be shown, they allow for interlopers intheir fitting routines. The first step of our analysis is the estimation of the GC or dSph(hereafter, galactic object for short) center. Although BALRoGOallows for an iterative frequentist approach, this can be dangerouswhen dealing with objects that suffer from crowding issues (Arenouet al. 2018), since the amount of stars decreases when one approachesthe center too much. Therefore, as for many other steps, we opt for aBayesian fit of the galactic object center.This is done by fitting a Plummer profile (Plummer 1911) plus aconstant contribution of interlopers to the surface density of the starsinside the two degrees cone (in some cases, we selected smaller radiidue to high MW contamination, or due to parasite satellites). Suchapproach aims for nothing more than finding a center of mass of thestellar distribution, by fitting the center that best suits an sphericalshaped collection of stars. Even if the galactic object does not followprecisely a Plummer model (e.g., core-collapsed GCs), other pro-files, such as Sérsic (Sérsic 1963; Sersic 1968), Hernquist (Hernquist1990) or (Kazantzidis et al. 2004) for example also follow a circularsymmetry, and therefore should have their centers well determinedby a Plummer fit. In addition, allowing for the contribution of inter-lopers in the fit enables the user not to lose many galactic object starsand thus have a better knowledge of the radial extent of the source,as well as the radial membership probability of a star with respect tothe galactic object. We use the first ten GCs from the Messier catalog, plus NGC 5139 ( 𝜔 Cen), NGC 6397 and NGC 6752. https://gitlab.com/eduardo-vitral/balrogo Figure 1.
Geometry of the center estimation:
The original cone search, cen-tered in ( 𝛼 S , 𝛿 S ) is shown as the gray circle of radius 𝑅 max , while the newfitted center is ( 𝛼 , 𝛿 ). The angles 𝜙 and 𝜙 represent the effective circularsections (of respective radii 𝑅 and 𝑅 ) where we consider circular symme-try. The distance between the two centers is labeled as 𝑑 , and we considerthis value to approach zero for the cone searches around the centers from theSIMBAD data base. Following the normalization of Vitral & Mamon (2020, equa-tions 15–17): Σ ( 𝑅 ) = 𝑁 ∞ 𝜋𝑅 (cid:101) Σ (cid:18) 𝑅𝑅 scale (cid:19) , (1a) 𝑁 ( 𝑅 ) = 𝑁 ∞ (cid:101) 𝑁 (cid:18) 𝑅𝑅 scale (cid:19) , (1b)where 𝑁 ∞ is the projected number of tracers at infinity, one can writethe global surface density as: (cid:101) Σ ( 𝑋 ) = (cid:101) Σ sys ( 𝑋 ) + R (cid:101) 𝑁 sys , tot 𝑋 − 𝑋 , (2)where R = 𝑁 ilop , tot / 𝑁 sys , tot is the ratio between the number ofinterlopers and system (galactic object) stars, (cid:101) 𝑁 sys , tot = (cid:101) 𝑁 sys ( 𝑋 max )− (cid:101) 𝑁 sys ( 𝑋 min ) is the total number of system stars, and 𝑋 = 𝑅 / 𝑎 is thenormalized projected radius, with 𝑎 being the Plummer scale radius,or the Plummer projected, two-dimensional, half number radius. Thenormalized surface density of the analyzed system, for a Plummerprofile reads: (cid:101) Σ sys ( 𝑋 ) = (cid:0) + 𝑋 (cid:1) , (3)while the normalized projected number of stars for the same modelis: (cid:101) 𝑁 sys ( 𝑋 ) = 𝜋 ∫ 𝑋 max 𝑋 min 𝜙 ( 𝑋 ) 𝑋 (cid:101) Σ sys ( 𝑋 ) d 𝑋 , (4)where 𝜙 ( 𝑋 ) is the angle corresponding to the effective circular sec- MNRAS000
The original cone search, cen-tered in ( 𝛼 S , 𝛿 S ) is shown as the gray circle of radius 𝑅 max , while the newfitted center is ( 𝛼 , 𝛿 ). The angles 𝜙 and 𝜙 represent the effective circularsections (of respective radii 𝑅 and 𝑅 ) where we consider circular symme-try. The distance between the two centers is labeled as 𝑑 , and we considerthis value to approach zero for the cone searches around the centers from theSIMBAD data base. Following the normalization of Vitral & Mamon (2020, equa-tions 15–17): Σ ( 𝑅 ) = 𝑁 ∞ 𝜋𝑅 (cid:101) Σ (cid:18) 𝑅𝑅 scale (cid:19) , (1a) 𝑁 ( 𝑅 ) = 𝑁 ∞ (cid:101) 𝑁 (cid:18) 𝑅𝑅 scale (cid:19) , (1b)where 𝑁 ∞ is the projected number of tracers at infinity, one can writethe global surface density as: (cid:101) Σ ( 𝑋 ) = (cid:101) Σ sys ( 𝑋 ) + R (cid:101) 𝑁 sys , tot 𝑋 − 𝑋 , (2)where R = 𝑁 ilop , tot / 𝑁 sys , tot is the ratio between the number ofinterlopers and system (galactic object) stars, (cid:101) 𝑁 sys , tot = (cid:101) 𝑁 sys ( 𝑋 max )− (cid:101) 𝑁 sys ( 𝑋 min ) is the total number of system stars, and 𝑋 = 𝑅 / 𝑎 is thenormalized projected radius, with 𝑎 being the Plummer scale radius,or the Plummer projected, two-dimensional, half number radius. Thenormalized surface density of the analyzed system, for a Plummerprofile reads: (cid:101) Σ sys ( 𝑋 ) = (cid:0) + 𝑋 (cid:1) , (3)while the normalized projected number of stars for the same modelis: (cid:101) 𝑁 sys ( 𝑋 ) = 𝜋 ∫ 𝑋 max 𝑋 min 𝜙 ( 𝑋 ) 𝑋 (cid:101) Σ sys ( 𝑋 ) d 𝑋 , (4)where 𝜙 ( 𝑋 ) is the angle corresponding to the effective circular sec- MNRAS000 , 1–13 (2021) lobular clusters and dwarf spheroidal galaxies tion where we analyze the data (see Fig 1), and whose analyticalexpression for small cone apertures (i.e. 𝑅 max (cid:28) 𝜙 ( 𝑋 = 𝑅 / 𝑎 ) = 𝜋 , if 𝑅 ≤ 𝑅 max − 𝑑 , (cid:34) 𝑅 + 𝑑 − 𝑅 𝑅 𝑑 (cid:35) , if 𝑅 > 𝑅 max − 𝑑 , (5)where 𝑅 max is the maximum radius of the original cone search (usu-ally two degrees) and 𝑑 is the distance between the fitted center andthe center from the original cone search ( 𝛼 S , 𝛿 S ), set as the sourcecenter on SIMBAD by the automatic Gaia advanced query. Theprojected radius 𝑅 is defined, in spherical trigonometry, as: 𝑅 = arccos [ sin 𝛿 sin 𝛿 + cos 𝛿 cos 𝛿 cos ( 𝛼 − 𝛼 )] , (6)where ( 𝛼 , 𝛿 ) is the center of the galactic object. The likelihoodfunction is therefore written as: L = (cid:214) 𝑖 𝜙 ( 𝑋 ) 𝜋 𝑋𝑎 (cid:101) Σ ( 𝑋 ) (cid:101) 𝑁 sys , tot ( + R) . (7)We minimize − log L with the differential_evolution routinefrom the scipy.optimize and find the centers that best fit the data,later displaying them in Table C1. For simplicity, since we expectthe quantity 𝑑 to approach zero, given the reliable previous measure-ments of the centers in SIMBAD, one can assume the case where 𝜙 ( 𝑋 ) = 𝜋 , and thus derive: (cid:101) 𝑁 sys ( 𝑋 ) = 𝑋 + 𝑋 . (8)Notably, for more general cases, the BALRoGO routine allows touse the more general representation of 𝜙 ( 𝑋 ) , and thus to account for amore complicated expression for (cid:101) 𝑁 sys ( 𝑋 ) , presented in appendix A. The next step of our algorithm is to fit the surface density of thegalactic object plus interlopers once again, using a maximum like-lihood estimation (MLE), but this time with fixed centers from theprevious step. It may seem as an unnecessary step given the previousone, but it actually allows the fitting routine to better explore the otherparameters ranges like the scale radius and the interloper fraction,and thus derive more robust results. In this step, we also reduced themaximum projected radius whenever there was a parasite galacticobject on the field or when the MW contamination was too intense,disturbing the fits.In this step, we not only tested Plummer profiles, but also Sérsicmodels (Sérsic 1963; Sersic 1968), which can be useful for core-collapsed clusters such as NGC 6397, and the Kazantzidis model,which is motivated from dynamical simulations of repeated tidalencounters (Kazantzidis et al. 2004). Equations 2, 6 and 7 remainvalid for these two other models, while the normalized surface densityand projected number of stars become: (cid:101) Σ sys ( 𝑋 ) = 𝑏 𝑛 ( 𝑛 ) 𝑛 Γ ( 𝑛 ) exp (cid:104) − 𝑏 ( 𝑛 ) 𝑋 / 𝑛 (cid:105) , (9a) http://simbad.u-strasbg.fr/simbad/ (cid:101) 𝑁 sys ( 𝑋 ) = 𝛾 ( 𝑛, 𝑏 ( 𝑛 ) 𝑋 / 𝑛 ) Γ ( 𝑛 ) . (9b)for Sérsic, where 𝛾 ( 𝑎, 𝑥 ) = ∫ 𝑥 𝑡 𝑎 − 𝑒 − 𝑡 d 𝑡 is the lower incompletegamma function and 𝑋 = 𝑅 / 𝑅 e , 𝑅 e being the effective projectedradius. And for Kazantzidis: (cid:101) Σ sys ( 𝑋 ) = 𝐾 ( 𝑋 )/ , (10a) (cid:101) 𝑁 sys ( 𝑋 ) = − 𝑋 𝐾 ( 𝑋 ) . (10b)Where 𝐾 and 𝐾 are modified Bessel functions of the second kindand 𝑋 = 𝑅 / 𝑎 K , / 𝑎 K being the radius of density slope − = AIC + 𝑁 free ( + 𝑁 free ) 𝑁 data − 𝑁 free − , (11)where AIC is the original Akaike Information Criterion (Akaike1973):AIC = − L MLE + 𝑁 free , (12)Our motivation for this criterion is the same one portrayed in Vitral& Mamon (2021), section 8.2. Whenever the difference between twomodels (i.e., Δ AICc = AICc i − AICc j ) was smaller than two (i.e., noparticular preference), priority was given to Plummer, Kazantzidisand Sérsic profiles, respectively. The results of our Bayesian Plummermodel plus constant interloper contribution are displayed in Figure 2for NGC 6752, NGC 6205 (M 13) and the Sculptor dSph, respectively. The two previous steps required no data cleaning, since the filteringcould bias the results, for example, by forcing radial incompletenesstowards the center of the galactic object (mostly because of crowd-ness issues) and cutting of fainter stars at the tracer outskirts, thusdamaging the fitting of the center and other surface density parame-ters. However, the next step, which is aimed to look into the propermotion space is, in its turn, much more sensible to poor astrome-try. This is a consequence of the much more subtle measurements ofproper motions than sky coordinates, given that the Gaia astrometricmission has not yet a long timeline data base to fully trust its mainvelocity values without any filtering. Another reason why a conser-vative data cleaning of proper motions is necessary is that our propermotion fits themselves do not take into account the measurement er-rors provided by Gaia, which will be latter justified. Therefore, ourdata cleaning follows a similar approach as done in Vitral & Mamon(2021) for the GC NGC 6397, detailed in the following.
We first filtered the Gaia stars inside a maximum projected radius:we selected stars inside a cone of radius 𝑟 cut , where 𝑟 cut = 𝑅 scale . 𝑅 scale was 𝑅 e , 𝑎 , and 𝑎 K for Sérsic, Plummer and Kazantzidis, re-spectively. For some sources, this limit was changed in order toaccount for a more representative subset, but never out-passed therange between 0.03 and two degrees. MNRAS , 1–13 (2021)
Eduardo Vitral − R [degree] Σ ( R ) [ d e g r ee − ] NGC 6752
Joint fitGC aloneField stars − R [degree] Σ ( R ) [ d e g r ee − ] NGC 6205
Joint fitGC aloneField stars − − R [degree] Σ ( R ) [ d e g r ee − ] Sculptor
Joint fitdSph aloneField stars
Figure 2.
Surface density fits:
Fits for the surface density of the galactic object (NGC 6752, NGC 6205 or M 13 and the Sculptor dSph, respectively) plusa constant contribution of interlopers, according to section 2.2. The histogram shows the empirical profile, using logarithmic radial bins extending from theinnermost bin point to 2 degrees. The curves show different models: our MLE fit ( red ) of a Plummer, Sérsic and Plummer models respectively ( dashed ) plusconstant field stars surface density ( dotted ), as well as the total ( solid ) to compare with the data.
The next step was to retain only stars with 𝜖 𝜇 < 𝜎 limit , where 𝜎 limit is an error threshold and 𝜖 𝜇 is the proper motion error (semi-majoraxis of the error ellipse), defined in Lindegren et al. (2018, eq. B.2): 𝜖 𝜇 = √︂ ( 𝐶 + 𝐶 ) + √︃ ( 𝐶 − 𝐶 ) + 𝐶 , (13) 𝐶 = 𝜖 𝜇 𝛼, ∗ , (14) 𝐶 = 𝜖 𝜇 𝛼, ∗ 𝜖 𝜇 𝛿 𝜌 , (15) 𝐶 = 𝜖 𝜇 𝛿 , (16)where 𝜖 denotes the error or uncertainty and where 𝜌 is the correlationcoefficient between 𝜇 𝛼, ∗ and 𝜇 𝛿 . The error threshold was derivedby default, as the following routine:(i) First, we selected the stars inside the maximum radius men-tioned above.(ii) With this subset, we generated a two-dimensional histogramin proper motion space and assigned the two highest local peaks ofthe distribution with Python’s routine peak_local_max from theskimage.feature method, which would correspond to the tracer andinterloper clumps of stars. The binning of the histogram was chosensuch that each dimension 𝑥 𝑖 had a number of bins computed as: 𝑁 bin = (cid:22) × ( max [ 𝑥 𝑖 ] − min [ 𝑥 𝑖 ]) min [ 𝜂 . ( 𝑥 𝑖 ) − 𝜂 . ( 𝑥 𝑖 ) , 𝜂 . ( 𝑥 𝑖 ) − 𝜂 . ( 𝑥 𝑖 )] (cid:25) , (17)where 𝜂 .𝑛 ( 𝑥 𝑖 ) is the 𝑛 -th percentile of the 𝑥 𝑖 data.(iii) We naively estimated the Gaussian dispersion of each clumpby taking the distance between the peaks found above and the closestpoint where the histogram reached e − / times the value at the peak.(iv) We assigned the interloper peak as the one with great disper-sion, while the narrow dispersion clump (calculated as above) wasconsidered as the tracer velocity dispersion 𝜎 tracer .(v) The error threshold 𝜎 limit was taken as 𝜎 limit = 𝑘 𝜎 tracer ,where the default 𝑘 was set as 0.5. For many clusters this value waschanged in order to select a more representative subset. The GCs had We use the standard notation 𝜇 𝛼 ∗ = cos 𝛿 [ d 𝛼 / d 𝑡 ] , 𝜇 𝛿 = d 𝛿 / d 𝑡 . all 𝑘 values between 0.2 and two, while the dSphs had 𝑘 up to five,given their greater distances, and consequently, much higher errors. We only kept stars whose astrometric solution presented a sufficientlylow goodness of fit: √︄ 𝜒 𝜈 (cid:48) − < . { , exp [− . ( 𝐺 − . )]} , (18)where 𝜈 (cid:48) is the number of points (epochs) in the astrometric fit of agiven star and 5 is the number of free parameters of the astrometric fit(2 for the position, 1 for the parallax and two for the proper motions).Eq. (18) gives a sharper HR diagram, removing artifacts such as dou-ble stars, calibration problems, and astrometric effects from binaries.It is more optimized than the ( astrometric_excess_noise < 𝐺 (cid:46) RUWE parameter.Second, we only kept stars with good photometry:1 . + . ( 𝐺 BP − 𝐺 RP ) < 𝐸 < . + . ( 𝐺 BP − 𝐺 RP ) , (19)where 𝐸 = ( 𝐼 BP + 𝐼 RP )/ 𝐼 G is the flux excess factor. Eq. (19) per-forms an additional filter in the HR diagram, removing stars withconsiderable photometric errors in the BP and RP photometry, af-fecting mainly faint sources in crowded areas. This poor photometrybroadens the CMD, leading to more confusion with field stars.Those variables correspond to the following quantities in the GaiaDR2 and EDR3 archive: • 𝜒 : astrometric_chi2_al • 𝜈 (cid:48) : astrometric_n_good_obs_al • 𝐸 : phot_bp_rp_excess_factor • 𝐺 BP − 𝐺 RP : bp_rp After having cleaned the data set according to the previous section,we proceed to fit the proper motion joint distribution of galactic
MNRAS , 1–13 (2021) lobular clusters and dwarf spheroidal galaxies −
10 0 10 µ α, ∗ [mas yr − ] − − µ δ [ m a s y r − ] NGC 6752 −
20 0 20 µ α, ∗ [mas yr − ] − − µ δ [ m a s y r − ] NGC 6205 −
20 0 20 µ α, ∗ [mas yr − ] − − − µ δ [ m a s y r − ] Sculptor −
10 0 10 [mas yr − ] . . . P D F PMs projected on semi − minor axis − −
10 0 10 [mas yr − ] . . . . P D F PMs projected on semi − minor axis [mas yr − ] . . . . P D F PMs projected on semi − minor axis − −
10 0 10 [mas yr − ] . . . P D F PMs projected on semi − major axis −
20 0 20 [mas yr − ] . . . . P D F PMs projected on semi − major axis −
20 0 20 40 [mas yr − ] . . . . P D F PMs projected on semi − major axis Figure 3.
Proper motion fits:
We display here the results of the proper motion fits of NGC 6752, NGC 6205 (M 13) and the Sculptor dSph in the first, secondand third columns, respectively. The first row displays the entire proper motion subset color-coded by stellar counts, from light blue to dark blue. The dashedellipse displays the Pearson VII asymmetric distribution, with its semi-minor and semi-major directions as two perpendicular dotted lines, while the continuouscircle represents the galactic object (globular cluster or dwarf spheroidal) proper motion mean with a radius equals eight times its intrinsic dispersion, for bettervisualization. The second and third rows display the fits (solid red) projected on the semi-minor and semi-major axis respectively, with the data in blue. objects plus interlopers. These two components present two clumpsof data in proper motion space which allow one to assign membershipprobabilities to each clump.Baumgardt et al. (2019) and Vasiliev (2019) used Bayesian ap-proaches based on Gaussian mixture models (a Gaussian distributionfor both cluster and interloper clumps), while the original work fromGaia Collaboration et al. (2018) used an iterative method based onclustering approaches. The use of Gaussian mixtures must, however,be taken with caution, since the distribution of field stars has oftenmuch wider wings than a Gaussian distribution. This was first noticedin Vitral & Mamon (2021, see appendix B), where the authors optedfor a symmetric Pearson VII distribution for the interlopers and theresults were much more robust. This trend is confirmed in this workwhere we employ a refined version of their Bayesian method, by as- signing a symmetric Gaussian distribution for the galactic object plusa non-symmetric Pearson VII distribution for interlopers. The choiceof a non-symmetric distribution allows for a much better adjustmentof the interlopers, which improves the membership probability of thestars in the subset, while the impact on the bulk proper motions ofthe galactic object varies from cluster to cluster .The main drawback of such an approach is that in order to intro-duce a non-symmetric Pearson VII distribution, which is analyticallymore complicated than a symmetric distribution, we abandon the For example, there is no strong difference when using symmetric and non-symmetric interloper distributions for some GCs such as NGC 6397 andNGC 6121 (M4) MNRAS , 1–13 (2021)
Eduardo Vitral convolution of our set with Gaussian errors. This was not an issuefor Vitral & Mamon (2021), who used a polynomial approximationto this convolution for the symmetric case, which was simpler, inorder to avoid extra integrals. This becomes much more computa-tionally costly whenever accepting the non symmetricity of the data.We therefore try to contour that we do not use the Gaia proper mo-tion errors in our fits by trusting in our data cleaning described in theprevious section, which applied conservative cuts on proper motionerrors. We further test this assumption in section 3.4.
To construct the probability distribution function (PDF) of the ana-lyzed subset, we consider the PDF from these tracers and from theMW contaminants:PDF = 𝑓 GO PDF GO + ( − 𝑓 GO ) PDF MW , (20)where 𝑓 GO is the fraction of galactic object stars. The PDF of galacticobjects is a straightforward Gaussian:PDF GO ( 𝜇 𝛼, ∗ , 𝜇 𝛿 ) = exp (cid:16) − 𝜁 (cid:17) 𝜋 𝜎 , (21a) 𝜁 = ( 𝜇 𝛼, ∗ − 𝜇 𝛼, ∗ GO ) + ( 𝜇 𝛿 − 𝜇 𝛿, GO ) 𝜎 (21b)where 𝜎 GO is the proper motion dispersion of the galactic objectstars, and 𝜇 𝛼, ∗ GO and 𝜇 𝛿, GO are their bulk proper motions in ( 𝛼 , 𝛿 ). For the PDF of interlopers, we first shift the origin to the bulkproper motions of the interlopers clump, and then rotate the referenceframe into the main axis of the interlopers proper motion ellipsoidaldistribution: 𝜇 𝑥 = ( 𝜇 𝛼, ∗ − 𝜇 𝛼, ∗ MW ) cos 𝜃 + ( 𝜇 𝛿 − 𝜇 𝛿, MW ) sin 𝜃 , (22a) 𝜇 𝑦 = −( 𝜇 𝛼, ∗ − 𝜇 𝛼, ∗ MW ) sin 𝜃 + ( 𝜇 𝛿 − 𝜇 𝛿, MW ) cos 𝜃 , (22b)where 𝜇 𝛼, ∗ MW and 𝜇 𝛿, MW are the contaminants bulk proper mo-tions in ( 𝛼 , 𝛿 ) and 𝜃 is the angle between the original ( 𝜇 𝛼, ∗ , 𝜇 𝛿 )frame and the new one. Then, if we call the semi-major and semi-minor axis of the Pearson VII ellipsoidal distribution 𝑎 𝑥 and 𝑎 𝑦 , wecan write the interlopers PDF as:PDF MW = Γ (cid:16) − − 𝜏 (cid:17) Γ (cid:16) − − 𝜏 (cid:17) (cid:26)(cid:20) + (cid:16) 𝜇 𝑥 𝑎 𝑥 (cid:17) (cid:21) (cid:20) + (cid:16) 𝜇 𝑦 𝑎 𝑦 (cid:17) (cid:21) (cid:27) ( + 𝜏 )/ 𝜋 𝑎 𝑥 𝑎 𝑦 , (23)where Γ ( 𝑥 ) is the gamma function of 𝑥 and 𝜏 is an intrinsic slope of thedistribution (with 𝜏 < − In order to derive errors of our Bayesian estimates, such as the bulkproper motion uncertainties shown in Table C1, we used Python’snumdifftools.Hessian method to compute the Hessian matrix ofthe proper motion PDF (i.e., eq [20]). After, we assigned the un-certainties of each parameter as the square root of the respectivediagonal position of the inverted Hessian matrix.
Once one has a precise analytical description of both the surfacedensity and the proper motion distribution of the ensemble of fieldstars plus the analyzed galactic object, it is much easier to extracttracer members by means of membership probabilities. However,the incredible amount of astrophysical information from the Gaiareleases allows to go even further, by analyzing how likely it is for astar to be part of the color-magnitude diagram of the tracer. It alsoallows to spot particular groups of stars such as binaries and bluestragglers (e.g., Leonard 1989) that lie away from the CMD.The last step of our default filtering routine aims at constrainingthe region covered by the galactic object on the color-magnitudediagram (CMD). Since we do not have an analytical form to cor-rectly describe the CMD, we opt, for the first time until now, to notuse a Bayesian method, but rather a non-parametric Kernel DensityEstimation (KDE) approach. We do so with the Python methodscipy.stats.gaussian_kde, similarly to Vitral & Mamon (2021).Nevertheless, our approach has some few differences: • We first select stars with a membership probability greater than0.8 (this limit is a modifiable parameter in BALRoGO). • We used a KDE bandwidth of half the one derived by the Sil-verman’s rule (Silverman 1986). • We used the membership probabilities of each star as weightsto the KDE routine. • We selected stars inside a 3 − 𝜎 density contour in the CMD (this 𝑛 − 𝜎 limit is a modifiable parameter in BALRoGO).The CMD of NGC 6752, NGC 6205 (M 13) and the Sculptor dSphare presented in Figure 4, with the 3 − 𝜎 contour region displayed asa black thick line. Another important parameter from the Gaia releases is the parallax,which can be converted to distances. In this article, we choose notto analyze this quantity since the errors on the parallax are stillconsiderably high for most of the sources in Table C1, and thereforedo not provide very reliable estimates of distance.It is possible, however, to perform such an analysis with BAL-RoGO, whose method parallax allows for a conversion of parallaxto distances, eliminating negative parallaxes and then computing themode of the distance distribution, also by means of non-parametricKDE approaches. This is done by generating a non-parametric func-tion which is the PDF of distances (for that reason, eliminating thenegative parallaxes can help to describe the PDF with more precision,given a certain KDE bandwidth), and maximizing it with Python’sdifferential_evolution method.
For astronomers interested in generating mock data sets, BALRoGOhas a mock method equipped with many capabilities such as (1) Pro-
MNRAS000
MNRAS000 , 1–13 (2021) lobular clusters and dwarf spheroidal galaxies B − R G NGC 6752 − − − KDE s PDF B − R G NGC 6205 − − − KDE s PDF B − R G Sculptor − − − KDE s PDF Figure 4.
Color-magnitude diagram fits:
We display the color-magnitude diagram (CMD) of NGC 6752, NGC 6205 (M 13) and the Sculptor dSph color-codedby the Kernel Density Estimation (KDE) non-parametric PDF. The black line contours indicate 3 − 𝜎 confidence regions, which we used to filter out interlopersand binaries that lie away from the CMD.
14 16 18 20 G mag . . . . . . H i s t og r a m o f G m ag ( m − / Data
14 16 18 20 G mag − − (cid:15) µ α , ∗ [ m a s y r − ] .
26 ( G mag − . Data
14 16 18 20 G mag − − (cid:15) µ δ [ m a s y r − ] .
26 ( G mag − . Data
Figure 5.
Gaia EDR3 uncertainties:
On the left , we show the cumulative histogram of Gaia EDR3 𝐺 magnitudes in blue, along with the function 10 ( 𝑚 − )/ in red. In the middle , we display the Gaia EDR3 𝜇 𝛼, ∗ uncertainties blue color-coded by stellar counts, with the function 10 . ( 𝑚 − . ) in red. In the right , wedisplay the Gaia EDR3 𝜇 𝛿 uncertainties blue color-coded by stellar counts, with the function 10 . ( 𝑚 − . ) in red. The Gaia EDR3 data used for those plotsis the stack of all stars in a two degrees cone search around the nearby globular clusters NGC 6121 (M 4), NGC 5139 ( 𝜔 Cen), NGC 6397 and NGC 6752. jecting cartesian data into sky coordinates , (2) Adding field starsuniformly distributed in an spherical cap and following a Pearson VIIdistribution in proper motions space, and (3) Adding realistic GaiaEDR3 errors to velocities. We describe below these three function-alities. BALRoGO’s method cart6d_to_gaia is able to convert cartesiancoordinates into plane of sky coordinates by making use of Astropy(Astropy Collaboration et al. 2013) coordinates method, allowingthe users to decide weather they want a realistic data set with propermotion uncertainties of not. Moreover, if the users want to shift a cer-tain source from one ( 𝛼 , 𝛿 ) position in the sky to another one, this canbe promptly done by calling the method angle.transrot_source, 𝛼 , 𝛿 , 𝜇 𝛼, ∗ and 𝜇 𝛿 which takes into account the spherical symmetry of the sky projec-tion. When observing regions of the sky narrow enough (such as the twodegrees cone searches we made in this work), one can expect tofind uniformly distributed field stars in the field of view. By takinginto account spherical trigonometry, one can generate such randomdistribution by inverting the probabilities Pr { 𝑟 < 𝑅 } and Pr { 𝜃 < Θ } ,where 𝑅 and Θ are the angular distance of a tracer from the source’scenter and the angle between this tracer and the source’s center withrespect to the increasing declination axis. We have, thus: 𝑅 = arccos [( − U ) + U cos ( 𝑅 lim )] , (24a) Θ = 𝜋 U , (24b) MNRAS , 1–13 (2021)
Eduardo Vitral where 𝑅 lim is the maximum distance from the source’s center andU is a random variable uniformly distributed between zero and one.Similarly, one is able to generate random proper motions from asymmetric Pearson VII distribution ( 𝑎 𝑦 = 𝑎 𝑥 in equation [23]) byinverting the probability Pr {| 𝜇 | < M } :M = 𝑎 √︁ U /( + 𝜏 / ) − , (25)where 𝑎 and 𝜏 are the scale radius and slope from the PearsonVII distribution, and once again U is a random variable uniformlydistributed between zero and one. We distribute those proper motionmoduli azimuthally by choosing angles from a distribution such as theone from equation (24b). The precise derivation of the probabilitiesabove is presented in appendix B. The proper motion uncertainties in Gaia EDR3 present a clear depen-dence on the apparent magnitude: In Figure 5, we stacked all GaiaEDR3 stars in a two degrees cone search around the nearby GCsNGC 6121 (M 4), NGC 5139 ( 𝜔 Cen), NGC 6397 and NGC 6752,in order to show that dependence. For that reason, to generate re-alistic Gaia EDR3 proper motion uncertainties, BALRoGO firstgenerates random magnitudes, again, by inverting the probabilityPr { 𝐺 mag < 𝑚 } , which can be easily derived from the cumula-tive distribution of 𝐺 mag magnitudes from Gaia EDR3 (i.e., the phot_g_mean_mag parameter), by imposing a threshold magnitude 𝑚 lim . The respective equations for a distribution of magnitude 𝑚 are:Pr { 𝐺 mag < 𝑚 } = ( 𝑚 − 𝑚 lim )/ 𝑓 , (26a) 𝑚 = 𝑓 log U + 𝑚 lim , (26b)where U is a random variable uniformly distributed between zeroand one, and 𝑚 lim =
21 and 𝑓 =
4, according to the fits displayedin Figure 5. Once one has a random set of magnitudes from equa-tion (26b), we can again make use of the fits displayed in Figure 5 toassign: 𝜖 𝜇 𝛼, ∗ = . ( 𝑚 − . ) , (27a) 𝜖 𝜇 𝛿 = . ( 𝑚 − . ) . (27b)Those uncertainties are provided to the user, after adding Gaussianerrors to 𝜇 𝛼, ∗ and 𝜇 𝛿 , with respective standard deviation of 𝜖 𝜇 𝛼, ∗ and 𝜖 𝜇 𝛿 . In Table C1, we display the main results of our fits, including ourBayesian estimates of the center, bulk proper motion and half numberradii for over a hundred GCs and a few dSph galaxies measured byGaia EDR3. This section aims to compare our results with previousestimates and discuss the implications of our new methods. R / [arcmin] , This work R / [ a r c m i n ] , L i t e r a t u r e Vasiliev 2019Baumgardt et al . Figure 6.
Effective radii:
Comparison between the half number radii in arcminderived by BALRoGO from Gaia EDR3 data in the x-axis, and the samequantity derived by Baumgardt et al. (2019, blue triangles) and Vasiliev(2019, red squares) from textscGaia DR2 data in the y-axis. We display the 𝑥 = 𝑦 line in dashed black. Our GC centers, derived according to section 2.1, were compared tothe estimates from Goldsbury et al. (2010) for a robustness check: Themedian separation between their centers and ours is of 0 .
12 arcsec,while the median of their reported uncertainties is of 0 . 𝑑 ≈
0, inFigure 1), and gives us confidence in our Bayesian center estimation.The maximum separation between our measurements and those fromGoldsbury et al. (2010) was of 26 . . Figure 6 displays the effective (two-dimensional) radii derived fromGaia EDR3 by BALRoGO in the x-axis along with the same quantityderived from Gaia DR2 by Vasiliev (2019) and Baumgardt et al.(2019), in red and blue respectively, in the y-axis. The conversion toprojected half number radii was straightforward for both Sérsic andPlummer fits, since the scale radius used in both models was alreadythe half number radius. For the Kazantizidis model, we multiplied theKazantizidis scale radius 𝑎 K from equations 10a and 10b by 1 . MNRAS000
0, inFigure 1), and gives us confidence in our Bayesian center estimation.The maximum separation between our measurements and those fromGoldsbury et al. (2010) was of 26 . . Figure 6 displays the effective (two-dimensional) radii derived fromGaia EDR3 by BALRoGO in the x-axis along with the same quantityderived from Gaia DR2 by Vasiliev (2019) and Baumgardt et al.(2019), in red and blue respectively, in the y-axis. The conversion toprojected half number radii was straightforward for both Sérsic andPlummer fits, since the scale radius used in both models was alreadythe half number radius. For the Kazantizidis model, we multiplied theKazantizidis scale radius 𝑎 K from equations 10a and 10b by 1 . MNRAS000 , 1–13 (2021) lobular clusters and dwarf spheroidal galaxies − . . . D i ff e r e n ce s (cid:2) m a s y r − (cid:3) Helmi et al . µ α, ∗ ∆ µ δ − . . . D i ff e r e n ce s (cid:2) m a s y r − (cid:3) Baumgardt et al . µ α, ∗ ∆ µ δ − . . . D i ff e r e n ce s (cid:2) m a s y r − (cid:3) Vasiliev 2019∆ µ α, ∗ ∆ µ δ N G C N G C N G C N G C N G C N G C N G C N G C N G C N G C N G C N G C N G C − . . . D i ff e r e n ce s (cid:2) m a s y r − (cid:3) Mean of previous estimates∆ µ α, ∗ ∆ µ δ Figure 7.
Proper motion means:
Comparison of proper motion means derivedfrom Gaia DR2 data by BALRoGO and other works in the literature (GaiaCollaboration et al. 2018, Baumgardt et al. 2019 and Vasiliev 2019) forthe first ten globular clusters from the Messier catalogue plus NGC 6397,NGC 6752 and NGC 5139 ( 𝜔 Cen). We display, in the y-axis, the differences Δ 𝜇 𝛼, ∗ = 𝜇 𝛼, ∗ this work − 𝜇 𝛼, ∗ other work and Δ 𝜇 𝛿 = 𝜇 𝛿 this work − 𝜇 𝛿 other work as red squares and blue triangles, respectively, for the 13 globular clustersmentioned above, distributed along the x-axis. In the bottom plot, we displayagain this difference, but with respect to the mean 𝜇 𝛼, ∗ and 𝜇 𝛿 from GaiaCollaboration et al. (2018), Baumgardt et al. (2019) and Vasiliev (2019). Theerrors bars were calculated as explained in section 3.3, and we plot a dashedblack horizontal line at the value of zero as a reference. One can notice thatthe differences moduli lie all bellow 0.05 mas yr − . any data filtering in this first step: It avoids a forced incompletenesstowards both the cluster center and outskirts, due to crowdness andfainter stars respectively, which are associated with worse astrometricsolutions that would be filtered out in most filtering routines. Gaia Collaboration et al. (2018), Baumgardt et al. (2019) and Vasiliev(2019) measured bulk proper motions for over a hundred GCs withGaia DR2 by using different methods previously mentioned. In thiswork, we update those information with the new Gaia EDR3, which
Table 1.
Comparison of estimates on the bulk 𝜇 𝛼, ∗ and 𝜇 𝛿 from the cleanand inaccurate mock data sets (see section 3.4).Data set 𝜇 𝛼 ∗ 𝜇 𝛿 [mas yr − ] [mas yr − ]Clean 4 . ± .
011 4 . ± . . ± .
012 4 . ± . has more precise and robust measurements of proper motions givenits longer baseline. Due to this fact, the values of proper motionsmeans for nearly all GCs changed by more than their respectiveerrors published in the works mentioned above. Therefore, it wouldbe unfair to compare our results using Gaia EDR3 with their resultsobtained from Gaia DR2 modelling.In order to make such a comparison possible, we decided to per-form the same analysis, but this time using Gaia DR2, for 13 GCs: thefirst ten GCs of the Messier catalogue, plus NGC 6397, NGC 6752and NGC 5139 ( 𝜔 Cen). We show, in Figure 7, the differences be-tween the bulk proper motions ( 𝜇 𝛼, ∗ , 𝜇 𝛿 ) estimated from Gaia DR2by BALRoGO and by the three methods cited above. In the bottompanel, we display again this difference, but with respect to the mean 𝜇 𝛼, ∗ and 𝜇 𝛿 from Gaia Collaboration et al. (2018), Baumgardt et al.(2019) and Vasiliev (2019). The uncertainty of those mean valueswas calculated as 𝜖 = (cid:104) 𝜖 𝑖 (cid:105) + 𝜎 , where 𝜖 𝑖 stands for the uncer-tainties on the estimated values, and 𝜎 stands for their standarddeviation. The error bars on those plots, in their turn, was calculatedas 𝜖 = √︃ 𝜖 + 𝜖 .We observe an outstanding agreement between the measurementsusing BALRoGO and the measurements from Gaia Collaborationet al. (2018), Baumgardt et al. (2019) and specially Vasiliev (2019).When taking the mean over their values, we see even better theconsistency of our method, with all the estimates agreeing within1 − 𝜎 . This gives us confidence on the use of a non-Gaussian mixturein our Bayesian fit, along with the choice of a Pearson VII non-symmetric distribution of proper motions for interlopers, even thoughwe neglected Gaia errors after cleaning the data. As previously pointed out, our conservative data cleaning from sec-tion 2.3 should be enough to contour the fact that we do not convolveequation (20) with Gaussian errors. In order to test this assumption,Pierre Boldrini kindly provided a GC mock using the initial condi-tion 𝑁 − body generator magi (Miki & Umemura 2018). Adoptinga distribution-function-based method, it ensures that the final real-ization of the cluster is in dynamical equilibrium (Miki & Umemura2018). This GC mock was inspired by the real cluster NGC 6397, andtherefore followed a Sérsic profile with stellar mass of 1 . × 𝑀 (cid:12) ,Sérsic radius of 3 .
14 pc, Sérsic index of 3 . .
91 kpc (Vasiliev 2019) and59 . stars in both subsets (hereafterclean and inaccurate subsets, for simplicity) and ran BALRoGO’s MNRAS , 1–13 (2021) Eduardo Vitral −
20 0 µ α, ∗ [mas yr − ] − − µ δ [ m a s y r − ] Mock data − − [mas yr − ] . . . . . P D F PMs projected on semi − minor axis −
10 0 10 [mas yr − ] . . . . P D F PMs projected on semi − major axis Figure 8.
Mock data:
We display here the results of the proper motion fits of our mock data set with realistic Gaia EDR3 errors. The image on the first columndisplays the entire proper motion subset color-coded by stellar counts, from light blue to dark blue. The dashed ellipse displays the fitted Pearson VII symmetricdistribution, with its main axis directions as two perpendicular dotted lines, while the continuous circle represents the galactic object (mock globular cluster)proper motion mean with a radius equals five times its intrinsic dispersion, for better visualization. The second and third columns display the fits (solid red)projected on the semi-minor and semi-major axis respectively, with the data in blue. routine on them. Figure 8 displays the proper motion fits of theinaccurate subset, similar to Figure 3.The fitted bulk 𝜇 𝛼, ∗ and 𝜇 𝛿 from the clean and inaccurate subsets,along with their respective uncertainties are displayed in Table 1. Wecan verify that even without convolving the global proper motionPDF with Gaussian errors, the fits on the bulk proper motion fromthe data sets with and without uncertainties agree within the 1 − 𝜎 error bars. This strengthens our confidence in the proper motioncleaning routine from BALRoGO, as well as in its fits. This work confirms the tendency of field stars proper motions tofollow a Pearson VII distribution, instead of a Gaussian, and weaddress now the interpretation of this result. Such an effect couldbe explained by the fact that the proper motions delivered by Gaiaare not a direct measurement of the velocity, i.e. a measure of spacevariation per time, but rather a measurement of angular velocity,which neglects the distance of the stars.According to the central limit theorem (Laplace 1810), when oneconsiders the ensemble of independent and identically distributedrandom variables, their properly normalized sum tends toward anormal distribution regardless of the variables original distribution.In the case of GCs and dSph, this can be considered as an adequateapproximation, since these sources contain generally tens or hundredsof thousands of stars (Binney & Tremaine 2008) which are locatedat distances that can be considered barely the same for a distantobserver. This means that when converting the spatial velocities ofits stars into angular velocities (i.e., dividing by their distance andturning them into proper motions), their originally quasi-independentand identically distributed velocities remain as such, and thereforethe variation around their mean is close to a Gaussian.In the case of MW field stars however, since they have completelydifferent distances, the distribution of proper motions is drawn awayfrom an independent and identically distributed assumption. In fact,their variation around their mean depends on the distance they liefrom the observer, some of them much closer and others much fartherthan the galactic object analyzed. As a consequence, we expect tofind more outliers, i.e., stars with proper motions that deviate more strongly from their mean, and thus a higher kurtosis (wider tails) inthe distribution. That is why the Pearson VII distribution, with itswider tails, is better adapted to fit the proper motions of interlopersthan a Gaussian distribution.
We present a new algorithm aimed at measuring some of the mainstructural parameters of galactic objects such as globular clustersand dwarf spheroidal galaxies. This algorithm, named BALRoGO:Bayesian Astrometric Likelihood Recovery of Galactic Ob-jects, performs Bayesian and non-parametric fits in order to extractstars drowned in Milky Way interlopers, by accessing their member-ship probabilities (with respect to the galactic object analyzed).Our approach presents innovative points which have been pre-viously used in Vitral & Mamon (2021), but with some new im-provements. Among the new approaches used in this algorithm, wehighlight the following: • Bayesian surface density fits considering a constant contributionfrom Milky Way interlopers. • We allow the surface density modelling of globular clusters withprofiles such as Sérsic and Kazantzidis, instead of the generally usedKing models (King 1966). • The proper motion Bayesian fit is not based on a Gaussianmixture, but rather a Pearson VII distribution for interlopers, giventheir distribution wider tails. • We use the membership probabilities of the surface density andproper motion fits in order to derive a non-parametric representationof the color-magnitude diagram, and then select cluster membersinside confidence regions from this representation.Comparisons between our method and previous works such asGaia Collaboration et al. (2018), Baumgardt et al. (2019) and Vasiliev(2019) indicate strong agreement of bulk proper motions for globularclusters using Gaia DR2, and of scale radii using Gaia EDR3. Inaddition, we make available our measurements of center, bulk propermotions and scale radii for over a hundred globular clusters fromthe NGC catalog, along with a few dwarf spheroidal galaxies in
MNRAS000
MNRAS000 , 1–13 (2021) lobular clusters and dwarf spheroidal galaxies Table C1, and in text format at this footnote . We also make all codesfrom BALRoGO available as an open source software at the addresshighlighted in this footnote .The dynamics of such stellar systems is a fundamental key tounderstand some of the main aspects of galaxy evolution, as well asthe astrophysical impacts of dark matter. With new releases of theGaia astrometric mission, along with other astrometric data sets, thefuture of dynamical modeling is very enticing. For that reason, webelieve it is important to make algorithms such as ours available,which can be easily used in order to derive important parameters ofmany stellar systems. ACKNOWLEDGEMENTS
First, I want to thank Alexandre Macedo for all his efforts towardsmaking BALRoGO public, readable and accessible by anyone. I ac-knowledge the important comments of Gary Mamon, in particularfor proposing equation (4) and providing equation (5) in a more com-plicated form than presented here. I also thank Pierre Boldrini, whoprovided the globular cluster mock data presented in section 3.4 andgave insightful comments on this work.Eduardo Vitral was funded by an AMX doctoral grant from ÉcolePolytechnique.This work greatly benefited from public software packages such asNumpy (van der Walt et al. 2011), emcee (Goodman & Weare 2010),Matplotlib (Hunter 2007), as well as the Spyder Integrated Devel-opment Environment (Copyright (c) 2009- Spyder Project Contrib-utors and others).
DATA AVAILABILITY
The data that support the plots within this paper and other find-ings of this study are available from the corresponding author uponreasonable request.
REFERENCES
Abadi M. G., Navarro J. F., Steinmetz M., 2006, MNRAS, 365, 747Akaike H., 1973, Information Theory and an Extension of the MaximumLikelihood Principle. Springer New York, New York, NY, pp 199–213Arenou F., et al., 2018, A&A, 616, A17Astropy Collaboration et al., 2013, A&A, 558, A33Battaglia G., Helmi A., Breddels M., 2013, New Astron. Rev., 57, 52Baumgardt H., Hilker M., Sollima A., Bellini A., 2019, MNRAS, 482, 5138Bertin G., Varri A. L., 2008, ApJ, 689, 1005Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. PrincetonUniversity Press, Princeton, NJ, USABoldrini P., Miki Y., Wagner A. Y., Mohayaee R., Silk J., Arbey A., 2020,MNRAS, 492, 5218Bullock J. S., Johnston K. V., 2005, ApJ, 635, 931Bustos Fierro I. H., Calderón J. H., 2019, MNRAS, 488, 3024Carretta E., et al., 2009, A&A, 505, 117Gaia Collaboration et al., 2018, A&A, 616, A12Galadi-Enriquez D., Jordi C., Trullols E., 1998, VizieR Online Data Catalog,pp J/A+A/337/125Goldsbury R., Richer H. B., Anderson J., Dotter A., Sarajedini A., WoodleyK., 2010, AJ, 140, 1830 https://gitlab.com/eduardo-vitral/balrogo/-/raw/master/table_gc_dsph.dat https://gitlab.com/eduardo-vitral/balrogo Goodman J., Weare J., 2010, Communications in Applied Mathematics andComputational Science, 5, 65Hernquist L., 1990, ApJ, 356, 359Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Kazantzidis S., Mayer L., Mastropietro C., Diemand J., Stadel J., Moore B.,2004, ApJ, 608, 663King I. R., 1966, AJ, 71, 64Laplace P., 1810, Mémoires de l’Académie Royale des Sciences de Paris, 10,301Leonard P. J. T., 1989, AJ, 98, 217Lindegren L., et al., 2018, A&A, 616, A2Marín-Franch A., et al., 2009, ApJ, 694, 1498Miki Y., Umemura M., 2018, MNRAS, 475, 2269Noyola E., Gebhardt K., Bergmann M., 2008, ApJ, 676, 1008Peñarrubia J., Walker M. G., Gilmore G., 2009, MNRAS, 399, 1275Pearson K., 1916, Philosophical Transactions of the Royal Society of LondonSeries A, 216, 429Peebles P. J. E., 1984, ApJ, 277, 470Peebles P. J. E., Dicke R. H., 1968, ApJ, 154, 891Plummer H. C., 1911, MNRAS, 71, 460Searle L., Zinn R., 1978, ApJ, 225, 357Sérsic J. L., 1963, Bull. Assoc. Argentina de Astron., 6, 41Sersic J. L., 1968, Atlas de galaxias australes. Cordoba, Argentina: Observa-torio AstronomicoSilverman B. W., 1986, Density estimation for statistics and data analysisSugiura N., 1978, Communications in Statistics - Theory and Methods, 7, 13Vasilevskis S., Sanders W. L., van Altena W. F., 1965, AJ, 70, 806Vasiliev E., 2019, MNRAS, 484, 2832Vitral E., Mamon G. A., 2020, A&A, 635, A20Vitral E., Mamon G. A., 2021, Astronomy & AstrophysicsWalker M., 2013, Dark Matter in the Galactic Dwarf Spheroidal Satellites.p. 1039, doi:10.1007/978-94-007-5612-0_20van der Marel R. P., Anderson J., 2010, ApJ, 710, 1063van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in ScienceEngineering, 13, 22
APPENDIX A: PROJECTED NUMBER
We present here the solution of equation 4, for a distribution followingthe Plummer (Plummer 1911) profile, in the approximation of smallcone apertures in the sky (i.e. 𝑅 max (cid:28) 𝑅 max and 𝑑 can be seen in Figure 1). For the case where 𝑅 ≤ 𝑅 max − 𝑑 ( 𝑑 alsodefined in Figure 1), one has equation 8.However, whenever 𝑅 > 𝑅 max − 𝑑 , the indefinite integral of equa-tion 4 yields: (cid:101) 𝑁 ( 𝑥 ) = − 𝑎 arccos (cid:104) 𝑅 + 𝑑 − 𝑅 𝑅 𝑑 (cid:105) 𝜋 ( 𝑎 + 𝑅 ) + (A1) (cid:32)√︁ Ξ (cid:32) − √︁ Ξ × arctan (cid:34) ( 𝑑 − 𝑅 ) − 𝑅 ( 𝑑 + 𝑅 )( 𝑑 − 𝑅 ) √ Ξ (cid:35) +( 𝑎 + 𝑑 − 𝑅 ) arctan (cid:20) Ξ √ Ξ √ Ξ (cid:21)(cid:33)(cid:33) ÷ (cid:16) 𝜋 √︁ Ξ √︁ Ξ (cid:17) + C . where C = 𝑎 is the Plummer effectiveradius, and Ξ , Ξ and Ξ are defined as: Ξ = − 𝑑 − ( 𝑅 − 𝑅 ) + 𝑑 ( 𝑅 + 𝑅 ) (A2a) MNRAS , 1–13 (2021) Eduardo Vitral
Figure B1.
Spherical geometry:
Representation of the physical situation of asource projected in the plane of sky. Ξ = 𝑎 + ( 𝑑 − 𝑅 ) + 𝑎 ( 𝑑 + 𝑅 ) (A2b) Ξ = ( 𝑑 − 𝑅 ) + 𝑎 ( 𝑑 + 𝑅 ) − 𝑅 ( 𝑎 + 𝑑 + 𝑅 ) (A2c)This treatment can be chosen in the BALRoGO method po-sition.find_center(), by providing the argument method="mle_robust" . APPENDIX B: FIELD STARS MOCK DATA
In this section, we derive equations (24) and (25).
B1 Random positions
In order to generate 𝑛 points uniformly distributed in a spherical cap,we shall first considerate Figure B1, with the geometry of the prob-lem. From classical spherical trigonometry relations, it is straight-forward to write:cos 𝑅 = sin 𝛿 sin 𝛿 + cos 𝛿 cos 𝛿 cos ( 𝛼 − 𝛼 ) , (B1a)sin 𝜙 = cos 𝛿 sin ( 𝛼 − 𝛼 ) sin 𝑅 , (B1b)sin 𝛿 = cos 𝑅 sin 𝛿 + sin 𝑅 cos 𝛿 cos 𝜙 . (B1c)We wish to generate an uniform distribution of points in a sphericalcap of radius 𝑅 lim , so the probability of having a radius smaller than 𝑅 can be written as:Pr { 𝑟 < 𝑅 } = Surface ( 𝑅 ) Surface ( 𝑅 lim ) , (B2)or, more precisely:Pr { 𝑟 < 𝑅 } = ∫ 𝑅 ∫ 𝜋 𝜌 sin 𝜃 d 𝜃 d 𝜑 ∫ 𝑅 lim ∫ 𝜋 𝜌 sin 𝜃 d 𝜃 d 𝜑 = − cos 𝑅 − cos 𝑅 lim , (B3) where 𝜌 is the radius of the sphere, 𝜃 and 𝜑 are the longitudinaland latitudinal angles of the sphere, that in our description have theirorigin at ( 𝛼 , 𝛿 ) . Similarly, the probability that 𝜙 is smaller than anangle Θ is:Pr { 𝜙 < Θ } = Θ / 𝜋 . (B4)Therefore, in order to derive equations (24), one just needs to invertthe relations (B3) and (B4), with respect to 𝑅 and Θ , respectively.The following step, which is to convert the set of known 𝛼 , 𝛿 , 𝑅 and 𝜙 into pairs of ( 𝛼 , 𝛿 ) is done by first deriving 𝛿 with equation (B1c)and then 𝛼 with: 𝛼 = 𝛼 + arccos ( Λ ) , for Λ > ≤ 𝜙 ≤ 𝜋𝛼 + arccos (− Λ ) , for Λ < ≤ 𝜙 ≤ 𝜋𝛼 − arccos ( Λ ) , for Λ > 𝜋 < 𝜙 < 𝜋𝛼 − arccos (− Λ ) , for Λ < 𝜋 < 𝜙 < 𝜋 (B5)where we have the following correspondences: Λ = √︄ − sin 𝜙 sin 𝑅 − sin 𝛿 (B6a) Λ = cos 𝑅 − sin 𝛿 sin 𝛿 cos 𝛿 cos 𝛿 (B6b) B2 Random proper motions
For the field stars, we generated random variables that followed asymmetric Pearson VII distribution. The statistical approach to doso was to invert the probability that a field star proper motion mod-ulus is smaller than M, i.e. the cumulative distribution function ofM for a symmetric Pearson VII distribution of scale radius 𝑎 andcharacteristic slope 𝜏 :CDF ( M ) ≡ ∫ M0 𝑓 PM ( 𝜇 ) d 𝜇 = − ∫ M0 𝜏 + 𝑎 𝜇𝑎 (cid:20) + (cid:16) 𝜇𝑎 (cid:17) (cid:21) 𝜏 / d 𝜇 . (B7)where 𝑓 PM is the distribution function of proper motions moduli forthe Pearson VII symmetric distribution. The result of the integralabove is:CDF ( M ) = − (cid:104) + ( M / 𝑎 ) (cid:105) + 𝜏 / . (B8)Thus, if U is a uniform random variable with boundaries [0, 1](and thus, 1 − U ≡ U, random PM variables can be generated withequation (25).
APPENDIX C: STRUCTURAL PARAMETERS OFGLOBULAR CLUSTERS AND DWARF SPHEROIDALGALAXIES
We present the main results of our paper in Table C1, such as effectiveradii and bulk proper motions for over a hundred globular clusters
MNRAS000
MNRAS000 , 1–13 (2021) lobular clusters and dwarf spheroidal galaxies and some of the main Local Group dwarf spheroidal galaxies. Thistable can be accessed in text format at the footnote link . This paper has been typeset from a TEX/L A TEX file prepared by the author. https://gitlab.com/eduardo-vitral/balrogo/-/raw/master/table_gc_dsph.dat MNRAS , 1–13 (2021) Eduardo Vitral
Table C1.
Catalog of proper motions and other dynamical parameters of NGC globular clusters and Local Group dwarf spheroidal galaxies.Name Other ID 𝛼 𝛿 𝐷 𝜇 𝛼, ∗ 𝜇 𝛿 𝜖 𝜇 𝛼, ∗ 𝜖 𝜇 𝛿 Σ 𝑅 / 𝜖 𝑅 / [deg] [deg] [kpc] [mas yr − ] [mas yr − ] [mas yr − ] [mas yr − ] [arcmin] [arcmin](1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)NGC 104 47 Tuc 6 . − . . − .
568 0 .
003 0 .
003 S 7 .
84 0 . . − . . − .
052 0 .
009 0 .
009 P 1 .
38 0 . . − . . − .
655 0 .
005 0 .
005 P 2 .
08 0 . . − . . − .
599 0 .
007 0 .
007 S 2 .
71 0 . . − . . − .
169 0 .
008 0 .
008 S 1 .
18 0 . . . − . − .
567 0 .
017 0 .
016 S 0 .
91 0 . . − . .
998 0 .
277 0 .
005 0 .
005 S 3 .
11 0 . . − . . − .
700 0 .
005 0 .
005 S 2 .
57 0 . . − . . − .
966 0 .
002 0 .
002 S 4 .
50 0 . . − . . − .
545 0 .
011 0 .
011 S 2 .
49 0 . . . − . − .
122 0 .
014 0 .
014 P 0 .
86 0 . . − . − .
409 3 .
307 0 .
005 0 .
005 S 4 .
00 0 . . − . − .
746 1 .
777 0 .
006 0 .
006 P 2 .
11 0 . . − . − . − .
959 0 .
006 0 .
006 S 2 .
48 0 . . . − . − .
274 0 .
008 0 .
008 S 2 .
47 0 . . . − . − .
209 0 .
011 0 .
011 S 2 .
29 0 . 𝜔 Cen 201 . − . − . − .
756 0 .
003 0 .
003 S 11 .
25 0 . . . − . − .
667 0 .
004 0 .
004 P 3 .
85 0 . . − . . − .
153 0 .
007 0 .
007 S 1 .
90 0 . . . − . − .
838 0 .
008 0 .
008 P 2 .
30 0 . . − . − . − .
409 0 .
019 0 .
020 S 1 .
50 0 . . − . − . − .
081 0 .
020 0 .
019 P 0 .
63 0 . . − . − . − .
234 0 .
011 0 .
010 S 1 .
11 0 . . − . − . − .
398 0 .
008 0 .
008 S 2 .
29 0 . . . . − .
866 0 .
004 0 .
004 S 4 .
26 0 . . − . − . − .
215 0 .
007 0 .
007 S 3 .
13 0 . . − . − . − .
649 0 .
010 0 .
010 S 1 .
05 0 . . − . − . − .
561 0 .
010 0 .
009 S 1 .
86 0 . . − . − . − .
588 0 .
010 0 .
010 P 1 .
80 0 . . − . . − .
271 0 .
006 0 .
006 P 2 .
27 0 . . − . − . − .
011 0 .
004 0 .
004 S 6 .
17 0 . . − . − . − .
698 0 .
010 0 .
010 S 1 .
40 0 . . − . − . − .
618 0 .
008 0 .
008 P 1 .
63 0 . . − . − . − .
970 0 .
006 0 .
006 S 2 .
09 0 . . . − . − .
570 0 .
003 0 .
003 S 3 .
87 0 . . − . − . − .
811 0 .
007 0 .
007 S 2 .
88 0 . . . − . − .
458 0 .
018 0 .
018 S 1 .
52 0 . . − . − . − .
595 0 .
010 0 .
010 S 1 .
03 0 . . − . − . − .
615 0 .
007 0 .
007 P 3 .
60 0 . . − . − . − .
640 0 .
019 0 .
018 K 1 .
17 0 . . − . − . − .
955 0 .
010 0 .
010 S 3 .
06 0 . . − . − .
252 1 .
656 0 .
007 0 .
007 S 2 .
62 0 . . − . − . − .
015 0 .
012 0 .
012 S 1 .
00 0 . . − . − . − .
875 0 .
013 0 .
013 P 1 .
16 0 . . − . . − .
327 0 .
011 0 .
011 P 1 .
51 0 . Notes : Columns are (1): Source name; (2): Alternative ID; (3): right ascension derived according to section 2.1; (4): declination derived according to section 2.1;(5): Distances in kpc from Baumgardt et al. (2019); (6): Bulk proper motion in right ascension (i.e., [ d 𝛼 / d 𝑡 ] cos 𝛿 , in mas yr − ); (7): Bulk proper motion indeclination (i.e., d 𝛿 / d 𝑡 , in mas yr − ); (8): uncertainty of 𝜇 𝛼, ∗ , in mas yr − ; (9): uncertainty of 𝜇 𝛿 , in mas yr − ; (10): Surface density model preferred by AICc(see section 2.2). ‘P’ stands for Plummer, ‘S’ for Sérsic and ‘K’ for Kazantzidis; (11): Effective (two-dimensional) radius, in arcmin; (12): uncertainty on theeffective radius, in arcmin.MNRAS000
51 0 . Notes : Columns are (1): Source name; (2): Alternative ID; (3): right ascension derived according to section 2.1; (4): declination derived according to section 2.1;(5): Distances in kpc from Baumgardt et al. (2019); (6): Bulk proper motion in right ascension (i.e., [ d 𝛼 / d 𝑡 ] cos 𝛿 , in mas yr − ); (7): Bulk proper motion indeclination (i.e., d 𝛿 / d 𝑡 , in mas yr − ); (8): uncertainty of 𝜇 𝛼, ∗ , in mas yr − ; (9): uncertainty of 𝜇 𝛿 , in mas yr − ; (10): Surface density model preferred by AICc(see section 2.2). ‘P’ stands for Plummer, ‘S’ for Sérsic and ‘K’ for Kazantzidis; (11): Effective (two-dimensional) radius, in arcmin; (12): uncertainty on theeffective radius, in arcmin.MNRAS000 , 1–13 (2021) lobular clusters and dwarf spheroidal galaxies Table C1 – continued Name Other ID 𝛼 𝛿 𝐷 𝜇 𝛼, ∗ 𝜇 𝛿 𝜖 𝜇 𝛼, ∗ 𝜖 𝜇 𝛿 Σ 𝑅 / 𝜖 𝑅 / [deg] [deg] [kpc] [mas yr − ] [mas yr − ] [mas yr − ] [mas yr − ] [arcmin] [arcmin](1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)NGC 6304 – 258 . − . − . − .
071 0 .
016 0 .
016 S 0 .
95 0 . . − . − . − .
600 0 .
015 0 .
015 P 0 .
96 0 . . − . − . − .
013 0 .
014 0 .
014 S 0 .
78 0 . . − . − . − .
229 0 .
008 0 .
008 S 2 .
30 0 . . . − . − .
630 0 .
004 0 .
004 S 2 .
66 0 . . − . − . − .
118 0 .
011 0 .
010 S 2 .
16 0 . . − . − . − .
443 0 .
006 0 .
006 S 3 .
00 0 . . − . − . − .
561 0 .
015 0 .
015 S 1 .
13 0 . . − . − . − .
401 0 .
008 0 .
008 S 2 .
01 0 . . − . − . − .
763 0 .
003 0 .
003 S 2 .
74 0 . . − . − . − .
159 0 .
005 0 .
005 S 3 .
38 0 . . − . − . − .
223 0 .
019 0 .
019 S 0 .
90 0 . . − . − . − .
707 0 .
008 0 .
008 S 2 .
48 0 . . − . . − .
654 0 .
003 0 .
003 S 5 .
77 0 . . − . − .
762 1 .
434 0 .
016 0 .
016 S 0 .
69 0 . . − . − . − .
070 0 .
008 0 .
008 S 2 .
21 0 . . . − . − .
996 0 .
012 0 .
012 P 0 .
83 0 . . − . − . − .
969 0 .
018 0 .
018 S 1 .
07 0 . . − . − . − .
386 0 .
009 0 .
009 P 1 .
05 0 . . − . − . − .
268 0 .
009 0 .
009 P 1 .
43 0 . . − . − . − .
466 0 .
014 0 .
013 S 1 .
15 0 . . − . . − .
475 0 .
024 0 .
024 S 0 .
46 0 . . − . − . − .
659 0 .
026 0 .
026 K 0 .
57 0 . . − . − . − .
948 0 .
012 0 .
012 K 1 .
06 0 . . − . − . − .
542 0 .
008 0 .
008 S 1 .
91 0 . . − . − . − .
806 0 .
026 0 .
023 P 0 .
40 0 . . − . . − .
844 0 .
006 0 .
005 S 2 .
71 0 . . − . − . − .
617 0 .
010 0 .
010 K 4 .
57 0 . . − . . − .
422 0 .
016 0 .
016 P 1 .
29 0 . . − . − . − .
150 0 .
020 0 .
019 K 0 .
61 0 . . − . − . − .
355 0 .
010 0 .
010 P 0 .
77 0 . . − . − . − .
223 0 .
012 0 .
011 S 1 .
32 0 . . − . . − .
958 0 .
014 0 .
013 P 1 .
24 0 . . − . − . − .
945 0 .
013 0 .
013 S 1 .
70 0 . . − . − . − .
831 0 .
009 0 .
009 P 1 .
73 0 . . − . − . − .
076 0 .
012 0 .
012 S 0 .
76 0 . . − . − . − .
911 0 .
017 0 .
017 P 0 .
70 0 . . − . − . − .
257 0 .
013 0 .
012 S 0 .
93 0 . . − . . − .
595 0 .
011 0 .
010 S 8 .
02 0 . . − . . − .
727 0 .
008 0 .
008 P 1 .
67 0 . . − . . − .
454 0 .
010 0 .
010 P 1 .
98 0 . . − . − . − .
380 0 .
003 0 .
003 S 2 .
73 0 . . − . − . − .
020 0 .
016 0 .
016 S 0 .
94 0 . . − . . − .
419 0 .
006 0 .
006 S 2 .
33 0 . . . − . − .
995 0 .
014 0 .
015 P 1 .
38 0 . . − . − . − .
035 0 .
003 0 .
003 P 5 .
17 0 . . . − . − .
615 0 .
009 0 .
009 S 2 .
32 0 . . . − .
007 1 .
616 0 .
006 0 .
006 S 1 .
60 0 . , 1–13 (2021) Eduardo Vitral
Table C1 – continued Name Other ID 𝛼 𝛿 𝐷 𝜇 𝛼, ∗ 𝜇 𝛿 𝜖 𝜇 𝛼, ∗ 𝜖 𝜇 𝛿 Σ 𝑅 / 𝜖 𝑅 / [deg] [deg] [kpc] [mas yr − ] [mas yr − ] [mas yr − ] [mas yr − ] [arcmin] [arcmin](1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)NGC 6809 M 55 294 . − . − . − .
315 0 .
003 0 .
003 S 3 .
94 0 . . . − . − .
657 0 .
004 0 .
004 S 3 .
37 0 . . − . − . − .
803 0 .
011 0 .
011 P 1 .
08 0 . . . − . − .
700 0 .
011 0 .
010 P 1 .
24 0 . . − . − . − .
378 0 .
016 0 .
014 P 1 .
13 0 . . . − . − .
646 0 .
016 0 .
016 P 0 .
44 0 . . . − . − .
808 0 .
005 0 .
005 P 3 .
12 0 . . − . . − .
172 0 .
007 0 .
007 S 2 .
54 0 . . − . − . − .
293 0 .
006 0 .
006 S 2 .
16 0 . . − . . − .
315 0 .
019 0 .
019 S 1 .
02 0 . . . − . − .
036 0 .
046 0 .
045 P 12 .
45 2 . . − . .
526 0 .
114 0 .
014 0 .
014 P 8 .
71 0 . . . . − .
191 0 .
012 0 .
012 S 6 .
98 0 . . − . . − .
342 0 .
003 0 .
003 S 14 .
80 0 . . . − . − .
116 0 .
018 0 .
017 S 2 .
93 0 . . . − . − .
148 0 .
039 0 .
038 S 2 .
41 0 . . − . . − .
143 0 .
005 0 .
005 P 9 .
71 0 . . − . − . − .
001 0 .
022 0 .
023 P 21 .
39 1 . . . − .
127 0 .
054 0 .
010 0 .
010 S 13 .
97 0 .000