- Split View
-
Views
-
Cite
Cite
S. Chatzopoulos, T. K. Fritz, O. Gerhard, S. Gillessen, C. Wegg, R. Genzel, O. Pfuhl, The old nuclear star cluster in the Milky Way: dynamics, mass, statistical parallax, and black hole mass, Monthly Notices of the Royal Astronomical Society, Volume 447, Issue 1, 11 February 2015, Pages 948–968, https://doi.org/10.1093/mnras/stu2452
- Share Icon Share
Abstract
We derive new constraints on the mass, rotation, orbit structure, and statistical parallax of the Galactic old nuclear star cluster and the mass of the supermassive black hole. We combine star counts and kinematic data from Fritz et al., including 2500 line-of-sight velocities and 10 000 proper motions obtained with VLT instruments. We show that the difference between the proper motion dispersions σl and σb cannot be explained by rotation, but is a consequence of the flattening of the nuclear cluster. We fit the surface density distribution of stars in the central 1000 arcsec by a superposition of a spheroidal cluster with scale ∼100 arcsec and a much larger nuclear disc component. We compute the self-consistent two-integral distribution function f(E, Lz) for this density model, and add rotation self-consistently. We find that (i) the orbit structure of the f(E, Lz) gives an excellent match to the observed velocity dispersion profiles as well as the proper motion and line-of-sight velocity histograms, including the double-peak in the vl-histograms. (ii) This requires an axial ratio near q1 = 0.7 consistent with our determination from star counts, q1 = 0.73 ± 0.04 for r < 70 arcsec. (iii) The nuclear star cluster is approximately described by an isotropic rotator model. (iv) Using the corresponding Jeans equations to fit the proper motion and line-of-sight velocity dispersions, we obtain best estimates for the nuclear star cluster mass, black hole mass, and distance M*(r < 100 arcsec) = (8.94 ± 0.31|stat ± 0.9|syst) × 106 M⊙, M• = (3.86 ± 0.14|stat ± 0.4|syst) × 106 M⊙, and R0 = 8.27 ± 0.09|stat ± 0.1|syst kpc, where the estimated systematic errors account for additional uncertainties in the dynamical modelling. (v) The combination of the cluster dynamics with the S-star orbits around Sgr A* strongly reduces the degeneracy between black hole mass and Galactic Centre distance present in previous S-star studies. A joint statistical analysis with the results of Gillessen et al., gives M• = (4.23 ± 0.14) × 106 M⊙ and R0 = 8.33 ± 0.11 kpc.
1 INTRODUCTION
Nuclear star clusters (NSC) are located at the centres of most spiral galaxies (Carollo et al. 1997; Böker et al. 2002). They are more luminous than globular clusters (Böker et al. 2004), have masses of the order of ∼106–107 M⊙ (Walcher et al. 2005), have complex star formation histories (Rossa et al. 2006; Seth et al. 2006), and obey scaling-relations with host galaxy properties as do central supermassive black holes (SMBH; Ferrarese et al. 2006; Wehner & Harris 2006); see Böker (2010) for a review. Many host an AGN, i.e. an SMBH (Seth et al. 2008), and the ratio of NSC to SMBH mass varies widely (Graham & Spitler 2009; Kormendy & Ho 2013).
The NSC of the Milky Way is of exceptional interest because of its proximity, about 8 kpc from the Earth. It extends up to several hundred arcseconds from the centre of the Milky Way (Sgr A*) and its mass within 1 pc is ∼106 M⊙ with ∼50 per cent uncertainty (Schödel, Merritt & Eckart 2009; Genzel, Eisenhauer & Gillessen 2010). There is strong evidence that the centre of the NSC hosts an SMBH of several million solar masses. Estimates from stellar orbits show that the SMBH mass is M• = (4.31 ± 0.36) × 106 M⊙ (Schödel et al. 2002; Ghez et al. 2008; Gillessen et al. 2009). Due to its proximity, individual stars can be resolved and number counts can be derived; however, due to the strong interstellar extinction the stars can only be observed in the infrared. A large number of proper motions (PMs) and line-of-sight (los) velocities have been measured, and analysed with spherical models to attempt to constrain the NSC dynamics and mass (Genzel et al. 1996, 2000; Haller et al. 1996; Trippe et al. 2008; Schödel et al. 2009; Fritz et al. 2014).
The relaxation time of the NSC within 1 pc is tr ∼ 1010 yr (Alexander 2005; Merritt 2013), indicating that the NSC is not fully relaxed and is likely to be evolving. One would expect from theoretical models that, if relaxed, the stellar density near the SMBH should be steeply-rising and form a Bahcall & Wolf (1976) cusp. In contrast, observations by Do et al. (2009), Buchholz, Schodel & Eckart (2009), and Bartko et al. (2010) show that the distribution of old stars near the SMBH appears to have a core. Understanding the NSC dynamics may therefore give useful constraints on the mechanisms by which it formed and evolved (Merritt 2010).
In this work, we construct axisymmetric Jeans and two-integral (2I) distribution function (DF) models based on stellar number counts, PMs, and los velocities. We describe the data briefly in Section 2; for more details, the reader is referred to the companion paper of Fritz et al. (2014). In Section 3, we carry out a preliminary study of the NSC dynamics using isotropic spherical models, in view of understanding the effect of rotation on the data. In Section 4, we describe our axisymmetric models and show that they describe the kinematic properties of the NSC exceptionally well. By applying a χ2 minimization algorithm, we estimate the mass of the cluster, the SMBH mass, and the NSC distance. We discuss our results and summarize our conclusions in Section 5. The appendix contains some details on our use of the Qian et al. (1995) algorithm to calculate the 2I-DF for the fitted density model.
2 DATASET
We first give a brief description of the data set used for our dynamical analysis. These data are taken from Fritz et al. (2014) and are thoroughly examined in that paper, which should be consulted for more details. The coordinate system used is a shifted Galactic coordinate system (l*, b*) where Sgr A* is at the centre and (l*, b*) are parallel to Galactic coordinates (l, b). In the following, we always refer to the shifted coordinates but will omit the asterisks for simplicity. The data set consists of stellar number densities, PMs, and los velocities. We use the stellar number density map rather than the surface brightness map because it is less sensitive to individual bright stars and non-uniform extinction.
The stellar number density distribution is constructed from NACO high-resolution images for Rbox < 20 arcsec, in a similar way as in Schödel et al. (2010), from HST WFC3/IR data for 20 arcsec < Rbox < 66 arcsec, and from VISTA-VVV data for 66 arcsec < Rbox < 1000 arcsec.
The kinematic data include PMs for ∼10 000 stars obtained from AO assisted images. The PM stars are binned into 58 cells (Fig. 1; Fritz et al. 2014) according to distance from Sgr A* and the Galactic plane. This binning assumes that the NSC is symmetric with respect to the Galactic plane and with respect to the b-axis on the sky, consistent with axisymmetric dynamical modelling. The sizes of the bins are chosen such that all bins contain comparable numbers of stars, and the velocity dispersion gradients are resolved, i.e. vary by less than the error bars between adjacent bins.
Relative to the large velocity dispersions at the Galactic Centre (100 km s−1), measurement errors for individual stars are typically ∼10 per cent, much smaller than in typical globular cluster PM data where they can be ∼50 per cent (e.g. in Omega Cen; van de Ven et al. (2006)). Therefore, corrections for these measurement errors are very small.
We also use ∼2 500 radial velocities obtained from SINFONI integral field spectroscopy. The binning of the radial velocities is shown in Fig. 2. There are 46 rectangular outer bins as shown in Fig. 2 plus 6 small rectangular rings around the centre (not shown; see appendix E of Fritz et al. 2014). Again the outer bins are chosen such that they contain similar numbers of stars and the velocity dispersion gradients are resolved. The distribution of radial velocity stars on the sky is different from the distribution of PM stars, and it is not symmetric with respect to l = 0. Because of this and the observed rotation, the binning is different, and extends to both positive and negative longitudes. Both the PM and radial velocity binning are also used in Fritz et al. (2014) and some tests are described in that paper.
Finally, we compare our models with (but do not fit to) the kinematics derived from about 200 maser velocities at r > 100 arcsec (from Lindquist et al. 1992; Deguchi et al. 2004). As for the PM and radial velocity bins, we use the mean velocities and velocity dispersions as derived in Fritz et al. (2014).
The assumption that the NSC is symmetric with respect to the Galactic plane and the b = 0 axis is supported by the recent Spitzer/IRAC photometry (Schödel et al. 2014) and by the distribution of PMs (Fritz et al. 2014). The radial velocity data at intermediate radii instead show an apparent misalignment with respect to the Galactic plane, by ∼10°; see Feldmeier et al. (2014) and Fritz et al. (2014). We show in Section 4.2 that, even if confirmed, such a misaligned structure would have minimal impact on the results obtained here with the symmetrized analysis.
3 SPHERICAL MODELS OF THE NSC
In this section, we study the NSC using the preliminary assumption that the NSC can be described by an isotropic DF depending only on energy. We use the DF to predict the kinematical data of the cluster. Later we add rotation self-consistently to the model. The advantages of using a DF instead of the common Jeans modelling are that (i) we can always check if a DF is positive and therefore if the model is physical, and (ii) the DF provides us with all the moments of the system. For the rest of the paper, we use (r, θ, φ) for spherical and (R, φ, z) for cylindrical coordinates, with θ = 0 corresponding to the z-axis normal to the equatorial plane of the NSC.
3.1 Mass model for the NSC
3.2 Spherical model
3.3 Adding self-consistent rotation to the spherical model
We describe here the effects of adding self-consistent rotation to the spherical model, but much of this also applies to the axisymmetric case which will be discussed in Section 4. We assume that the rotation axis of the NSC is aligned with the rotation axis of the Milky Way disc. We also use a Cartesian coordinate system (x, y, z), where z is parallel to the axis of rotation as before, y is along the los, and x is along the direction of negative longitude, with the centre of the NSC located at the origin. The PM data are given in Galactic longitude l and Galactic latitude b angles, but because of the large distance to the centre, we can assume that x∥l and z∥b.
Whether a spherical system can rotate has been answered in Lynden Bell (1960). Here, we give a brief review. Rotation in a spherical or axisymmetric system can be added self-consistently by reversing the sense of rotation of some of its stars. Doing so, the system will remain in equilibrium. This is equivalent with adding to the DF a part that is odd with respect to Lz. The addition of an odd part does not affect the density (or the mass) because the integral of the odd part over velocity space is zero. The most effective way to add rotation to a spherical system is by reversing the sense of rotation of all of its counterrotating stars. This corresponds to adding f−(E, L2, Lz) = sign(Lz)f(E, L2) (Maxwell's daemon; Lynden Bell 1960) to the initially non-rotating DF, and generates a system with the maximum allowable rotation. The general case of adding rotation to a spherical system can be written f′(E, L2, Lz) = (1 + g(Lz))f(E, L2), where g(Lz) is an odd function with max |g(Lz)| < 1 to ensure positivity of the DF. We notice that the new DF is a three-integral DF. In this case, the density of the system is still rotationally invariant but f− is not.
In Fig. 5, we notice that the projected velocity dispersion in the l direction is larger than the dispersion in the b direction which was first found by Trippe et al. (2008). This is particularly apparent for distances larger than 10 arcsec. A heuristic attempt to explain this difference was made in Trippe et al. (2008) where they imposed a rotation of the form υφ(r, θ) along with their Jeans modelling, as a proxy for axisymmetric modelling. Here, we show that for a self-consistent system the difference in the projected l and b dispersions cannot be explained by just adding rotation to the cluster.
An alternative way to see this is by making a particle realization of the initial DF (e.g. Aarseth, Henon & Wielen 1974). Then we can add rotation by reversing the sign of Lz of a percentage of particles using some probability function which is equivalent to changing the signs of υx and υy of those particles. |$\overline{\upsilon _x^2}$| will not be affected by the sign change and the |$\overline{\upsilon }_x^2$| averaged over the los will be zero because for each particle at the front of the system rotating in a specific direction, there will be another particle at the rear of the system rotating in the opposite direction. In this work, we do not use particle models to avoid fluctuations due to the limited number of particles near the centre.
4 AXISYMMETRIC MODELING OF THE NSC
The reduced χ2 that corresponds to these parameter values is χ2/νSD = 0.99 for νSD = 110 d.o.f. Here, we note that there is a strong correlation between the parameters a2 and M2. The flattening of the inner component is very similar to the recent determination from Spitzer/IRAC photometry (0.71 ± 0.02; Schödel et al. 2014) but slightly more flattened than the best value given by Fritz et al. (2014), 0.80 ± 0.04. The second component is about 100 times more massive than the first, but also extends more than one order of magnitude further.
4.1 Axisymmetric Jeans modelling
In order to define our model completely, we need to determine the distance R0 and mass M* of the cluster and the black hole mass M•. To do this, we apply a χ2 minimization technique matching all three velocity dispersions in both sets of cells, using the following procedure. First, we note that the inclusion of self-consistent rotation to the model will not affect its mass. This means that for the fitting we can use |${\overline{ {\upsilon _{\rm los}^2} }^{1/2}}$| for each cell of Fig. 2. Similarly, since our model is axisymmetricm we should match to the |${ \overline{\upsilon ^2_{l,b}}^{1/2}}$| for each PM cell; the |${{\overline{\upsilon }_{l,b}}}$| terms should be and indeed are negligible. Another way to see this is that since the system is axially symmetric, the integration of |${{\overline{\upsilon }_{l,b}}}$| along the los should be zero because the integration would cancel out for positive and negative y.
With this in mind we proceed as follows, using the cluster's density parameters1 as given in (19). First, we partition the 3D space (R0, M*, M•) into a grid with resolution 20 × 20 × 20. Then for each point of the grid, we calculate the corresponding χ2 using the velocity dispersions from all cells in Figs 1 and 2, excluding the two cells at the largest radii (see Fig. 8). We compare the measured dispersions with the model values obtained from equations (24) for the centres of these cells. Then, we interpolate between the χ2 values on the grid and find the minimum of the interpolated function, i.e. the best values for (R0, M*, M•). To determine statistical errors on these quantities, we first calculate the Hessian matrix from the curvature of χ2 surface at the minimum, ∂χ2/∂pi∂pj. The statistical variances will be the diagonal elements of the inverted matrix.
First, we now look at the comparison of this model with the velocity data. Fig. 8 shows how the azimuthally averaged dispersions σl and σb compare with the measured PM dispersions. Fig. 9 shows how this best model, similarly averaged, compares with the los mean square velocity data. The maser data are also included in the plot. It is seen that the model fits the data very well, in accordance with its χ2/νJeans = 1.07 per cell. Fig. 10 shows how all three projected dispersions of the model compare. σb is slightly lower than σlos due to projection effects. The fact that all three velocity dispersion profiles in Figs 8 and 9 are fitted well by the model suggests that the assumed semi-isotropic dynamical structure is a reasonable approximation.
The model prediction in Fig. 8 is similar to fig. 11 of Trippe et al. (2008) but the interpretation is different. As shown in the previous section, the difference in projected dispersions cannot be explained by imposing rotation on the model. Here, we demonstrated how the observational finding σl > σb can be quantitatively reproduced by flattened axisymmetric models of the NSC and the surrounding nuclear disc.
Compared to the best and more flattened model, the cluster mass has increased and the black hole mass has decreased. The sum of both masses has changed only by 2 per cent and the distance only by 1 per cent. In Figs 8 and 9, we see how the projected velocity dispersions of this model compare with our best model. The main difference seen in σl comes from the different flattening of the inner component, and the smaller slope of the dispersions near the centre of the new model is because of its smaller central density slope.
4.2 Distance to the Galactic Center, mass of the star cluster, and mass of the black hole
We now consider the determination of these parameters from the NSC data in more detail. Fig 11 shows the marginalized χ2-plot for the NSC model as given in equation (19), for pairs of two parameters (R0, M•), (M•, M*), (R0, M*), as obtained from fitting the Jeans dynamical model to the velocity dispersion profiles. The figure shows contour plots for constant χ2/νJeans with 1σ, 2σ, and 3σ in the three planes for the two-dimensional distribution of the respective parameters. We notice that the distance R0 has the smallest relative error.
The best-fitting values for (R0, M*, M•) are given in equation (25); these values are our best estimates based on the NSC data alone. For the dynamical model with these parameters and the surface density parameters given in equation (19), the flattening of the inner component inferred from the surface density data is consistent with the dynamical flattening, which is largely determined by the ratio of σl/σb and the tensor virial theorem.
Statistical errors are determined from the Hessian matrix for this model. Systematic errors can arise from uncertainties in the NSC density structure, from deviations from the assumed axisymmetric 2I dynamical structure, from dust extinction within the cluster (see Section 5), and other sources. We have already illustrated the effect of varying the cluster flattening on (R0, M•, M*) with our second, rounder model. We have also tested how variations of the cluster density structure (a2, q2, M2) beyond 500 arcsec impact the best-fitting parameters, and found that these effects are smaller than those due to flattening variations.
We have additionally estimated the uncertainty introduced by the symmetrization of the data if the misalignment found by Feldmeier et al. (2014) and Fritz et al. (2014) were intrinsic to the cluster, as follows. We took all radial velocity stars and rotated each star by 10° clockwise on the sky. Then, we resorted the stars into our radial velocity grid (Fig. 2). Using the new values |${\overline{ {\upsilon _{\rm los}^2} }^{1/2}}$| obtained in the cells, we fitted Jeans models as before. The values we found for R0, M*, and M• with these tilted data differed from those in equation (25) by ΔR0 = −0.02 kpc, ΔM*(m < 100 arcsec) = −0.15 × 106 M⊙, and ΔM• = +0.02 × 106 M⊙, respectively, which are well within the statistical errors.
We note several other systematic errors which are not easily quantifiable and so are not included in these estimates, such as inhomogeneous sampling of PMs or los velocities, extinction within the NSC, and the presence of an additional component of dark stellar remnants.
Based on our best model, the mass of the star cluster within 100 arcsec converted into spherical coordinates is M*(r < 100 arcsec) = (8.94 ± 0.32|stat ± 0.9|syst) × 106 M⊙. The model's mass within the innermost pc (25 arcsec) is M*(m < 1pc) = 0.729 × 106 M⊙ in spheroidal radius, or M*(r < 1pc) = 0.89 × 106 M⊙ in spherical radius. The total mass of the inner NSC component is M1 = 6.1 × 107 M⊙. Because most of this mass is located beyond the radius where the inner component dominates the projected star counts, the precise division of the mass in the model between the NSC and the adjacent nuclear disc is dependent on the assumed slope of the outer density profile of NSC, and is therefore uncertain.
The distance and the black hole mass we found differ by 0.7 and 12 per cent, respectively, from the values R0 = 8.33 ± 0.17|stat ± 0.31|syst kpc and M• = 4.31 ± 0.36 × 106 M⊙ for R0 = 8.33 kpc, as determined by Gillessen et al. (2009) from stellar orbits around Sgr A*. Fig. 12 shows the 1σ–3σ contours of marginalized χ2 for (R0, M•) jointly from stellar orbits (Gillessen et al. 2009), for the NSC model of this paper, and for the combined modelling of both data sets. The figure shows that both analyses are mutually consistent. When marginalized over M* and the respective other parameter, the combined modelling gives, for each parameter alone, R0 = 8.33 ± 0.11 kpc and M• = 4.23 ± 0.14 × 106 M⊙. We note that these errors for R0 and M• are both dominated by the distance error from the NSC modelling. Thus, our estimated additional systematic error of 0.1 kpc for R0 in the NSC modelling translates to a similar additional error in the combined R0 measurement and, through the SMBH mass–distance relation given in Gillessen et al (2009), to an additional uncertainty ≃ 0.1 × 106 M⊙ in M•. We see that the combination of the NSC and S-star orbit data is a powerful means for decreasing the degeneracy between the SMBH mass and Galactic Centre distance in the S-star analysis.
4.3 2I-DF for the NSC
The typical shape of the VPs in the l and b directions within the area of interest (r < 100 arcsec) is shown in Fig. 14. We notice the characteristic two-peak shape of the VP along l that is caused by the near-circular orbits of the flattened system. Because the front and the back of the axisymmetric cluster contribute equally, the two peaks are mirror-symmetric, and adding rotation would not change their shapes.
The middle panels of Fig. 15 and Figs B1 and B2 in Appendix show how our best model (with parameters as given in equations 19 and 27) predicts the observed VHs for various combinations of cells. The reduced χ2 for each histogram is also provided. The prediction is very good both for the VPs in υl and υb. Specifically, for the l PMs our flattened cluster model predicts the two-peak structure of the data pointed out by several authors (Trippe et al. 2008; Schödel et al. 2009; Fritz et al. 2014). In order to calculate the VP from the model for each cell, we averaged over the VP functions for the centre of each cell weighted by the number of stars in each cell and normalized by the total number of stars in all the combined cells.
Fig. 15 compares two selected υl-VPs for our two main models with the data. The left-hand column shows how the observed VHs for corresponding cells compare to the model VPs for the less flattened model with parameters given in equations (27) and (28), the middle column compares with the same VPs from our best model with parameters given in equations (19) and (25). Clearly, the more flattened model with q1 = 0.73 fits the shape of the data much better than the more spherical model with q1 = 0.85, justifying its use in Section 4.2.
This model is based on an even DF in Lz and therefore does not yet have rotation. To include rotation, we will (in Section 4.4) add an odd part to the DF, but this will not change the even parts of the model's VPs. Therefore, we can already see whether the model is also a good match to the observed los velocities by comparing it to the even parts of the observed los VHs. This greatly simplifies the problem since we can think of rotation as independent, and can therefore adjust it to the data as a final step. Fig. B3 shows how the even parts of the VHs from the los data compare with the VPs of the 2I model. Based on the reduced χ2, the model provides a very good match. Possible systematic deviations are within the errors. The los VHs are broader than those in the l direction because the los data contain information about rotation (the broader the even part of the symmetrized los VHs, the more rotation the system possesses, and in extreme cases they would show two peaks).
4.4 Adding rotation to the axisymmetric model: Is the NSC an isotropic rotator?
As in the spherical case, to model the rotation we add an odd part in Lz to the initial even part of the DF, so that the final DF takes the form f(E, Lz) = (1 + g(Lz))f(E, Lz). We again use equation (14); this adds two additional parameters (κ, F) to the DF. Equation (16) gives the mean los velocity versus Galactic longitude. In order to constrain the parameters (κ, F) we fitted the mean los velocity from equation (16) to the los velocity data for all cells in Fig. 2. The best parameter values resulting from this 2D-fitting are κ = 2.8 ± 1.7, F = 0.85 ± 0.15, and |$\chi _{r}^2=1.25$|. Fig. B4 shows that the VPs of this rotating model compare well with the observed los VHs.
5 DISCUSSION
In this work, we presented a dynamical analysis of the Milky Way's NSC, based on ∼10 000 PMs, ∼2 700 radial velocities, and new star counts from the companion paper of Fritz et al. (2014). We showed that an excellent representation of the kinematic data can be obtained by assuming a constant mass-to-light ratio for the cluster, and modelling its dynamics with axisymmetric 2I-DFs, f(E, Lz). The DF modelling allows us to see whether the model is physical, i.e. whether the DF is positive, and to model the PM and los VHs. One open question until now has been the nature of the double-peaked VHs of the vl-velocities along Galactic longitude, and the bell-shaped VHs of vb along Galactic latitude, which cannot be fitted by Gaussians (Schödel et al. 2009). Our 2I-DF approximation of the NSC gives an excellent prediction for the observed shapes of the vl-, vb-, and vlos-VHs. The models show that the double-peaked shape of the vl-VHs is a result of the flattening of the NSC, and suggest that the cluster's dynamical structure is close to an isotropic rotator. Because both PMs and los-velocities enter the dynamical models, we can use them also to constrain the distance to the GC, the mass of the NSC, and the mass of the Galactic Centre black hole. To do this efficiently, we used the semi-isotropic Jeans equations corresponding to 2I-DFs. In this section, we discuss these issues in more detail.
5.1 The dynamical structure of the NSC
The star count map derived in Fritz et al. (2014) suggests two components in the NSC density profile, separated by an inflection point at about ∼200 arcsec ∼ 8 pc (see Fig. 7 above). To account for this, we constructed a two-component dynamical model for the star counts in which the two components are described as independent γ-models. The inner, rounder component can be considered as the proper NSC, as in Fritz et al. (2014), while the outer, much more flattened component may represent the inner parts of the nuclear stellar disc (NSD) described in Launhardt, Zylka & Mezger (2002).
The scale radius of the inner component is ∼100 arcsec, close to the radius of influence of the SMBH, rh ∼ 90 arcsec (Alexander 2005). The profile flattens inside ∼20 arcsec to a possible core (Buchholz et al. 2009; Fritz et al. 2014) but the slope of the three-dimensional density profile for the inner component is not well constrained.
The flattening for the inner NSC component inferred from star counts is q1 = 0.73 ± 0.04, very close to the value of q = 0.71 ± 0.02 found recently from Spitzer multiband photometry (Schödel et al. 2014). It is important that these determinations agree with the dynamical flattening of our best Jeans dynamical models: the dynamical flattening is robust because it is largely determined by the ratio of σb/σl and the tensor virial theorem. Because star counts, photometric, and dynamical values for the inner NSC flattening agree, this parameter can now be considered securely determined.
Assuming constant mass-to-light ratio for the NSC, we found that a 2I-DF model gives an excellent description of the PM and los velocity dispersions and VHs, in particular of the double-peaked distributions in the vl-velocities. This double-peaked structure is a direct consequence of the flattening of the star cluster; the detailed agreement of the model VPs with the observed histograms therefore confirms the value q1 = 0.73 for the inner cluster component. For an axisymmetric model rotation cannot be seen directly in the PM VHs when observed edge-on, as is the case here, but is apparent only in the los velocities. When a suitable odd part of the DF is added to include rotation, the 2I-DF model also gives a very good representation of the skewed los VHs. From the amplitude of the required rotation, we showed that the NSC can be approximately described as an isotropic rotator model, rotating slightly slower than that outside ∼30 arcsec.
Individual VHs are generally fitted by this model within the statistical errors, but on closer examination the combined vl VHs show a slightly lower peak at negative velocities, as already apparent in the global histograms of Trippe et al. (2008) and Schödel et al. (2009). Fig. 17 suggests that differential extinction of the order of ∼0.2 mag within the cluster may be responsible for this small systematic effect, by causing some stars from the back of the cluster to fall out of the sample. The dependence of mean extinction on vl independently shows that the NSC must be rotating, which could otherwise only be inferred from the los velocities. In subsequent work, we will model the effect of extinction on the inferred dynamics of the NSC. This will then also allow us to estimate better how important deviations from the 2I-dynamical structure are, i.e. whether the three-integral dynamical modelling (e.g. De Lorenzi et al. 2013) would be worthwhile.
5.2 Mass of the NSC
The dynamical model results in an estimate of the mass of the cluster from our data set. Our fiducial mass value is M*(m < 100 arcsec) = (7.73 ± 0.31|stat ± 0.8|syst) × 106 M⊙ interior to a spheroidal major axis distance m = 100 arcsec. This corresponds to an enclosed mass within a three-dimensional radius r = 100 arcsec of M*(r < 100 arcsec) = (8.94 ± 0.31|stat ± 0.9|syst) × 106 M⊙.
The fiducial mass M*(r < 100 arcsec) for the best axisymmetric model is larger than that obtained with spherical models. The constant M/L spherical model with density parameters as in Section 3, for R0 = 8.3 kpc and the same black hole mass has M*(r < 100 arcsec) = 6.6 × 106 M⊙.
There are two reasons for this difference: (i) at ∼50 arcsec where the model is well fixed by kinematic data the black hole still contributes more than half of the interior mass. In this region, flattening the cluster at constant mass leaves σl and σlos approximately constant, but decreases σb to adjust to the shape. To fit the same observed data, the NSC mass must be increased. (ii) Because of the increasing flattening with radius, the average density of the axisymmetric model decreases faster than that of the spherical density fit; thus for the same observed velocity dispersion profiles, a larger binding mass for the NSC is required.
Fig. 18 shows the enclosed stellar mass within the spheroidal radius m as in equation (26), as well as the mass within the spherical radius r. For example, the mass within 1 pc (25 arcsec) is |${M_*}(r < 1{\rm pc \sim 25\,{\rm arcsec}})=0.89 \times {10^6}{\,\mathrm{M}_{\odot }}$|. This is compatible with the spherical modelling of Schödel et al. (2009) who gave a range of 0.6–1.7 × 106 M⊙, rescaled to R0 = 8.3 kpc, with the highest mass obtained for their isotropic, constant M/L model. According to Fig. 18, at r ≃ 30 arcsec = 1.2 pc the NSC contributes already ≃ 25 per cent of the interior mass (≃ 45 per cent at r ≃ 50 arcsec = 2 pc), and beyond r ≃ 100 arcsec = 4 pc it clearly dominates.
An important point to note is that the cluster mass does not depend on the net rotation of the cluster but only on its flattening. This is because to add rotation self-consistently to the model we need to add an odd part to the DF which does not affect the density or the PM dispersions σl and σb.
Our NSC mass model can be described as a superposition of a moderately flattened nuclear cluster embedded in a highly flattened nuclear disc. The cumulative mass distributions of the two components are shown in the middle panel of Fig. 18. The NSD starts to dominate at about 800 arcsec which is in good agreement with the value found by Launhardt et al. (2002).
Approximate local axial ratios for the combined density and for the total potential including the central black hole are shown in the lower panel of Fig. 18. Here we approximate the axial ratio of the density at radius R by solving the equation ρ(R, 0) = ρ(0, z) for z and writing qρ = z/R, and similarly for |$q_\Psi$|. The density axial ratio qρ(R) shows a strong decrease between the regions dominated by the inner and outer model components. The equipotentials are everywhere less flattened. At the centre, |$q_\Psi = 1$| because of the black hole; the minimum value is not yet reached at 1000 arcsec. Therefore, we can define the NSC proper as the inner component of this model, similar to Fritz et al. (2014).
The total mass of the inner component, M1 = 6.1 × 107 M⊙ (Section 4.2), is well determined within similar relative errors as M*(m < 100 arcsec). However, identifying M1 with the total mass of the Galactic NSC at the centre of the nuclear disc has considerable uncertainties: because the outer NSD component dominates the surface density outside 100–200 arcsec, the NSC density profile slope at large radii is uncertain, and therefore the part of the mass outside ∼200 arcsec (∼64 per cent of the total) is also uncertain. A minimal estimate for the mass of the inner NSC component is its mass within 200 arcsec up to where it dominates the star counts. This gives MNSC > 2 × 107 M⊙.
Finally, we use our inferred dynamical cluster mass to update the K-band mass-to-light ratio of the NSC. The best-determined mass is within 100 arcsec. Comparing our M*(r < 100 arcsec) = (8.94 ± 0.31|stat ± 0.9|syst) × 106 M⊙ with the K-band luminosity of the old stars derived in Fritz et al. (2014), L100 arcsec = (12.12 ± 2.58) × 106 L⊙, Ks, we obtain M/LKs = (0.76 ± 0.18) M⊙/L⊙, Ks. The error is dominated by the uncertainty in the luminosity (21 per cent, compared to a total 10 per cent in mass from adding statistical and systematic errors in quadrature). The inferred range is consistent with values expected for mostly old, solar metallicity populations with normal IMF (e.g. Courteau et al. 2014; Fritz et al. 2014).
5.3 Evolution of the NSC
After ∼10 half-mass relaxation times trh a dense NSC will eventually evolve to form a Bahcall–Wolf cusp with slope γ = 7/4 (Merritt 2013); for rotating dense star clusters around black holes, this was studied by Fiestas & Spurzem (2010). The minimum allowable inner slope for a spherical system with a black hole to have a positive DF is γ = 0.5. From the data, it appears that the Galactic NSC instead has a core (Buchholz et al. 2009; Fritz et al. 2014), with the number density possibly even decreasing very close to the centre (r < 0.3 pc). This is far from the expected Bahcall–Wolf cusp, indicating that the NSC is not fully relaxed. It is consistent with the relaxation time of the NSC being of order 10 Gyr everywhere in the cluster (Merritt 2013).
From Fig. 16 we see that the rotational properties of the Milky Way's NSC are close to those of an isotropic rotator. Fiestas et al. (2012) found that relaxation in rotating clusters causes a slow (∼3trh) evolution of the rotation profile. Kim et al. (2008) found that it also drives the velocity dispersions towards isotropy; in their initially already nearly isotropic models, this happens in ∼4trh. On a similar time-scale the cluster becomes rounder (Einsel & Spurzem 1999). Comparing with the NSC relaxation time suggests that these processes are too slow to greatly modify the dynamical structure of the NSC, and thus that its properties were probably largely set up at the time of its formation.
The rotation-supported structure of the NSC could be due to the rotation of the gas from which its stars formed, but it could also be explained if the NSC formed from merging of globular clusters. In the latter model, if the black hole is already present, the NSC density and rotation after completion of the merging phase reflects the distribution of disrupted material in the potential of the black hole (e.g. Antonini et al. 2012). Subsequently, relaxation would lead to shrinking of the core by a factor of ∼2 in 10 Gyr towards a value similar to that observed (Merritt 2010). In the simulations of Antonini et al. (2012), the final relaxed model has an inner slope of γ = 0.45, not far from our models (note that in flattened semi-isotropic models the minimum allowed slope for the density is also 0.5; Qian et al. 1995). Their cluster also evolved towards a more spherical shape, however, starting from a configuration with much less rotation and flattening than we inferred here for the present Milky Way NSC. Similar models with a net rotation in the initial distribution of globular clusters could lead to a final dynamical structure more similar to the Milky Way NSC.
5.4 Distance to the Galactic Centre
From our large PM and los velocity data sets, we obtained a new estimate for the statistical parallax distance to the NSC using the axisymmetric Jeans modelling based on the cluster's inferred dynamical structure. From matching our best dynamical model to the PM and los velocity dispersions within approximately |l|, |b| < 50 arcsec, we found R0 = 8.27 ± 0.09|stat ± 0.1|syst kpc. The statistical error is very small, reflecting the large number of fitted dispersion points. The systematic modelling error was estimated from uncertainties in the density structure of the NSC, as discussed in Section 4.2.
Our new distance determination is much more accurate than that of Do et al. (2013) based on anisotropic spherical Jeans models of the NSC, |$R_0=8.92_{-0.58}^{+0.58}$| kpc, but is consistent within their large errors. We believe this is mostly due to the much larger radial range we modelled, which leaves less freedom in the dynamical structure of the model.
The new value for R0 is in the range R0 = 8.33 ± 0.35 kpc found by Gillessen et al. (2009) from analysing stellar orbits around Sgr A*. A joint statistical analysis of the NSC data with the orbit results of Gillessen et al. (2009) gives a new best value and error R0 = 8.33 ± 0.11 kpc (Fig. 12, Section 4.2). Our estimated systematic error of 0.1 kpc for R0 in the NSC modelling translates to a similar additional uncertainty in this combined R0 measurement.
Measurements of R0 prior to 2010 were reviewed by Genzel et al. (2010). Their weighted average of direct measurements is R0 = 8.23 ± 0.20 ± 0.19 kpc, where the first error is the variance of the weighted mean and the second the unbiased weighted sample variance. Two recent measurements give R0 = 8.33 ± 0.05|stat ± 0.14|syst kpc from RR Lyrae stars (Dekany et al. 2013) and R0 = 8.34 ± 0.14 kpc from fitting axially symmetric disc models to trigonometric parallaxes of star forming regions (Reid et al. 2014). These measurements are consistent with each other and with our distance value from the statistical parallax of the NSC, with or without including the results from stellar orbits around Sgr A*, and the total errors of all three measurements are similar, ∼2 per cent.
5.5 Mass of the Galactic SMBH
Given a dynamical model, it is possible to constrain the mass of the central black hole from 3D stellar kinematics of the NSC alone. With the axisymmetric Jeans modelling we found M• = (3.86 ± 0.14|stat ± 0.4|syst) × 106 M⊙, where the systematic modelling error is estimated from the difference between models with different inner cluster flattening as discussed in Section 4.2. Within errors this result is in agreement with the black hole mass determined from stellar orbits around Sgr A* (Gillessen et al. 2009).
Our data set for the NSC is the largest analysed so far, and the axisymmetric dynamical model is the most accurate to date; it compares well with the various PM and los VHs. None the less, future improvements may be possible if the uncertainties in the star density distribution and kinematics within 20 arcsec can be reduced, the effects of dust are incorporated, and possible deviations from the assumed 2I-axisymmetric dynamical structure are taken into account.
Several similar analyses have been previously made using the spherical isotropic or anisotropic modelling. Trippe et al. (2008) used the isotropic spherical Jeans modelling for PMs and radial velocities in 1 arcsec < R < 100 arcsec; their best estimate is M• ∼ 1.2 × 106 M⊙, much lower than the value found from stellar orbits. Schödel et al. (2009) constructed isotropic and anisotropic spherical broken power-law models, resulting in a black hole mass of |${M_ \bullet } = 3.6_{-0.4}^{+0.2}\times {10^6}{\,\mathrm{M}_{\odot } }$|. However, Fritz et al. (2014) find M• ∼ 2.27 ± 0.25 × 106 M⊙, also using a power-law tracer density. They argue that the main reason for the difference to Schödel et al. (2009) is because their velocity dispersion data for R > 15 arcsec are more accurate, and their sample is better cleaned for young stars in the central R < 2.5 arcsec. Assuming an isotropic spherical model with constant M/L, Fritz et al. (2014) find M• ∼ 4.35 ± 0.12 × 106 M⊙. Do et al. (2013) used 3D stellar kinematics within only the central 0.5 pc of the NSC. Applying the spherical Jeans modelling, they obtained |${M_\bullet } = 5.76_{-1.26}^{+1.76}\times {10^6}{\,\mathrm{M}_ {\odot } }$| which is consistent with that derived from stellar orbits inside 1 arcsec, within the large errors. However, in their modelling they used a very small density slope for the NSC, of γ = 0.05, which does not correspond to a positive DF for their quasi-isotropic model.
Based on this work and our own models in Section 4, the black hole mass inferred from NSC dynamics is larger for constant M/L models than for power-law models, and it increases with the flattening of the cluster density distribution.
The conceptually best method to determine the black hole mass is from stellar orbits close to the black hole (Schödel et al. 2002; Ghez et al. 2008; Gillessen et al. 2009), as it requires only the assumption of Keplerian orbits and is therefore least susceptible to systematic errors. Gillessen et al. (2009) find that the largest uncertainty in the value obtained for M• is due to the uncertainty in R0, and that M• scales as |$M_\bullet \propto R_0^{2.19}$|. Therefore using our improved statistical parallax for the NSC also leads to a more accurate determination of the black hole mass. A joint statistical analysis of the axisymmetric NSC modelling together with the orbit modelling of Gillessen et al. (2009) gives a new best value and error for the black hole mass, M• = (4.26 ± 0.14) × 106 M⊙ (see Fig. 12, Section 4.2). An additional systematic error of 0.1 kpc for R0 in the NSC modelling, through the BH mass–distance relation given in Gillessen et al (2009), translates to an additional uncertainty ≃0.1 × 106 M⊙ in M•.
Combining this result with the mass modelling of the NSC, we can give a revised value for the black hole influence radius rinfl, using a common definition of rinfl as the radius where the interior mass M(<r) of the NSC equals twice the black hole mass (Merritt 2013). Comparing the interior mass profile in Fig. 18 as determined by the dynamical measurement with M• = 4.26 × 106 M⊙, we obtain rinfl ≃ 94 arcsec = 3.8 pc.
The Milky Way is one of some 10 galaxies for which both the masses of the black hole and of the NSC have been estimated (Kormendy & Ho 2013). From these it is known that the ratio of both masses varies widely. Based on the results above we estimate the Milky Way mass ratio M•/MNSC = 0.12 ± 0.04, with the error dominated by the uncertainty in the total NSC mass.
6 CONCLUSIONS
Our results can be summarized as follows.
The density distribution of old stars in the central 1000 arcsec in the Galactic Centre can be well approximated as the superposition of a spheroidal NSC with a scalelength of ∼100 arcsec and a much larger nuclear disc (NSD) component.
The difference between the PM dispersions σl and σb cannot be explained by rotation alone, but is a consequence of the flattening of the NSC. The dynamically inferred axial ratio for the inner component is consistent with the axial ratio inferred from the star counts which for our two-component model is q1 = 0.73 ± 0.04.
The orbit structure of an axisymmetric 2I-DF f(E, Lz) gives an excellent match to the observed double-peak in the vl-PM VHs, as well as to the shapes of the vertical vb-PM histograms. Our model also compares well with the symmetrized (even) los VHs.
The rotation seen in the los velocities can be modelled by adding an odd part of the DF, and this shows that the dynamical structure of the NSC is close to an isotropic rotator model.
Fitting PMs and los dispersions to the model determines the NSC mass within 100 arcsec, the mass of the SMBH, and the distance to the NSC. From the star cluster data alone, we find M*(r < 100 arcsec) = (8.94 ± 0.31|stat ± 0.9|syst) × 106 M⊙, M• = (3.86 ± 0.14|stat ± 0.4|syst) × 106 M⊙, and R0 = 8.27 ± 0.09|stat ± 0.1|syst kpc, where the estimated systematic errors account for additional uncertainties in the dynamical modelling. The fiducial mass of the NSC is larger than in previous spherical models. The total mass of the NSC is significantly more uncertain due to the surrounding nuclear disc; we estimate MNSC = (2–6) × 107 M⊙. The mass of the black hole determined with this approach is consistent with results from stellar orbits around Sgr A*. The Galactic Centre distance agrees well with recent accurate determinations from RR Lyrae stars and masers in the Galactic disc, and has similarly small errors.
Combining our modelling results with the stellar orbit analysis of Gillessen et al. (2009), we find M• = (4.23 ± 0.14) × 106 M⊙ and R0 = 8.33 ± 0.11 kpc. Because of the better constrained distance, the accuracy of the black hole mass is improved as well. Combining with the parameters of the cluster, the black hole radius of influence is 3.8 pc (= 94 arcsec) and the ratio of black hole to cluster mass is estimated to be 0.12 ± 0.04.
It is computationally too expensive to simultaneously also minimize χ2 over the density parameters.
REFERENCES
APPENDIX A: TWO-INTEGRAL DISTRIBUTIONS FUNCTIONS
In this part we give implementation instructions for the 2I-DF algorithm of HQ. We will try to focus on the important parts of the algorithm and also on the tests that one has to make to ensure that the implementation works correctly. Our implementation is based on Qian et al. (1995) and made with mathematica (Wolfram Research, Inc. 2011). For the theory, the reader should consider the original HQ paper.
Fig. A2 shows the shape of the DF for η = 0.5 for the potential we use in the fourth section of the paper for one value of h, using the aforementioned procedure. We notice that for large energies fluctuations of the DF appear. In order to solve this, we introduce a minor improvement of the procedure, by generalizing the h value of the contour to an energy-dependent function h = h(E). The h(E) could be a simple step function that takes four or five different values. For our model, the h(E) function is a decreasing function of E. This means that the minor axis of the ellipse should decrease as the E increases to avoid such fluctuations. In general, we can write h = h(E, Lz) so that the contour depends both on E and Lz.
Once we implement the algorithm it is necessary to test it. Our first test is to check that the lower half of the integration path in Fig. A1 is the complex conjugate of the upper half. Probably, the next most straightforward test is against the spherical case. It is possible to use the HQ algorithm to calculate a DF for spherical system. This DF should be equal to that obtained from Eddington's formula for the same parameters. After calculating our 2I-DF we compare its low-order moments with those of the Jeans modelling.