Triaxial orbit-based modelling of the Milky Way Nuclear Star Cluster
We construct triaxial dynamical models for the Milky Way nuclear star cluster using Schwarzschild’s orbit superposition technique. We fit the stellar kinematic maps presented in Feldmeier et al. (2014). The models are used to constrain the supermassive black hole mass M, dynamical mass-to-light ratio , and the intrinsic shape of the cluster. Our best-fitting model has M = () , = (0.90) /, and a compression of the cluster along the line-of-sight. Our results are in agreement with the direct measurement of the supermassive black hole mass using the motion of stars on Keplerian orbits. The mass-to-light ratio is consistent with stellar population studies of other galaxies in the mid-infrared. It is possible that we underestimate M and overestimate the cluster’s triaxiality due to observational effects. The spatially semi-resolved kinematic data and extinction within the nuclear star cluster bias the observations to the near side of the cluster, and may appear as a compression of the nuclear star cluster along the line-of-sight. We derive a total dynamical mass for the Milky Way nuclear star cluster of M = (2.1) 10 within a sphere with radius = 2 = 8.4 pc. The best-fitting model is tangentially anisotropic in the central = 0.5-2 pc of the nuclear star cluster, but close to isotropic at larger radii. Our triaxial models are able to recover complex kinematic substructures in the velocity map.
keywords:Galaxy: center; kinematics and dynamics.
The Milky Way nuclear star cluster is the ideal object to study the dynamics of a stellar system around a supermassive black hole. At a distance of 8 kpc it is close enough to resolve the individual stars, and measure discrete velocities in three dimensions. Modelling the stellar kinematics can constrain the mass distribution of the star cluster, and reveal the presence of a central dark massive object. In the special case of our own Galaxy, it it possible to observe Keplerian orbits of stars around a dark, point-mass-like object in the Galactic centre. These observations constrain this dark object to be a supermassive black hole with a mass of (4.1 0.6) (Ghez et al., 2008), (4.3 0.39) (Gillessen et al., 2009), or (4.02 0.20) (Boehle et al., 2016). Unfortunately, similar high-resolution observations are not yet possible in other galaxies.
Already in the 1970s the requirement of a central supermassive black hole in the Galactic centre was discussed to explain observational data (e.g. Oort, 1977). Several studies used stellar radial velocities to constrain the mass distribution in the Galactic centre (e.g. Rieke & Rieke, 1988; McGinn et al., 1989; Sellgren et al., 1990; Haller et al., 1996; Genzel et al., 1996). Also stellar proper motions were used to study the Galactic centre mass distribution (Schödel et al., 2009). Several studies combined radial velocity and proper motion data (Trippe et al., 2008; Do et al., 2013; Fritz et al., 2016). The mass distribution was derived using the spherical Jeans (1922) equations or the projected mass estimators of Bahcall & Tremaine (1981) for spherical systems. These studies found that a central dark mass of 25 10 is required to explain the observations.
Together with the increase of observational data, also the modelling became more advanced. Trippe et al. (2008) included the rotation of the nuclear star cluster in the modelling, although the rotation velocity of their data was too high (Schödel et al., 2009; Feldmeier et al., 2014). Feldmeier et al. (2014) and Chatzopoulos et al. (2015a) studied the Milky Way nuclear star cluster using axisymmetric Jeans models. Chatzopoulos et al. (2015a) showed the advantages of axisymmetric models over spherical Jeans models, which cannot explain the observed asymmetry of the velocity dispersion of proper motions parallel and perpendicular to the Galactic plane. The nuclear star cluster appears to be flattened in its light distribution (Schödel et al., 2014) as well as in the kinematics (Chatzopoulos et al., 2015a). Most studies showed that the nuclear star cluster kinematics is in agreement with isotropy (Schödel et al., 2009; Do et al., 2013; Chatzopoulos et al., 2015a), although the uncertainties are quite large (e.g. Do et al. 2013). All these models assumed a constant mass-to-light ratio for the light distribution of the cluster.
In this study we relax the assumption of axisymmetry and use triaxial orbit-based Schwarzschild (1979) models. Orbit-based models make no assumptions on the velocity anisotropy of the stellar motions, as Jeans models do. Further, the higher moments of the kinematics can also be included (Rix et al., 1997), which is important to break the degeneracy of mass and anisotropy in dynamical models.
Orbit-based models are commonly used to analyse line-of-sight velocity data of other galaxies (e.g. van der Marel et al., 1998; Gebhardt et al., 2000; Valluri et al., 2005; van den Bosch et al., 2008), and are an excellent tool to detect and measure the masses of supermassive black holes and dark matter halos. For extragalactic systems, the data are usually obtained from integrated light observations. Each data point contains the accumulated kinematics of many stars, weighted by their respective brightness. However, modelling the dynamics of integrated light data may be prone to systematic uncertainties, and bias the results of the central black hole mass. Therefore, it is interesting to test dynamical models on systems for which we know the central black hole mass from other independent measurements. The Milky Way nuclear star cluster is a good object for this kind of test. Also megamaser disc galaxies are useful to validate stellar dynamical black hole measurements. Black hole mass measurements from megamasers are very precise with uncertainties of only about 10 per cent. However, there is currently only one megamaser disc galaxy with a stellar dynamical black hole mass measurement (van den Bosch et al., 2016), NGC 4258. Different dynamical studies found either a 15 per cent lower or a 25 per cent higher black hole mass than the maser measurement (Siopis et al., 2009; Drehmer et al., 2015).
We use the triaxial orbit-based code by van den Bosch et al. (2008) to model the light distribution and line-of-sight kinematics of the Milky Way nuclear star cluster. We use the spectroscopic data cube constructed by Feldmeier et al. (2014) for the kinematic data, and derive a surface brightness distribution using Spitzer 4.5 and NACO band images. We assume a galactocentric distance of 8 kpc (Malkin, 2012) and a position angle 3140 East of North (J2000.0 coordinates, Reid & Brunthaler 2004) with respect to the Galactic plane. This paper is organised as follows: We describe the kinematic and photometric data in Section 2. The dynamical models are introduced in Section 3. Section 4 discusses the results, and Section 5 summarizes the main conclusions.
2 Description of the Data
2.1 Kinematic data
The line-of-sight velocity distribution (LOSVD) provides constraints on the dynamical structure of stellar systems. To extract this information, we used the near-infrared band spectroscopic data cube of Feldmeier et al. (2014), which has a pixel scale of 2.22 arcsec pixel. We used the data cube that was cleaned of foreground stars and bright stars. The cleaned data cube contains the light of the old red giant star population.
We fitted the LOSVD as in Feldmeier et al. (2014), i.e. on the stellar CO absorption lines (2.29022.365 ) with the IDL routine pPXF (Cappellari & Emsellem, 2004) and the high resolution spectra of Wallace & Hinkle (1996) as template stars. We applied the same spatial binning as Feldmeier et al. (2014), resulting in 175 spatial bins. Feldmeier et al. (2014) fitted only the velocity and velocity dispersion . However, we fitted in addition also higher moments of the LOSVD, in particular the Gauss-Hermite parameters and . We added noise to each of the 175 integrated light spectra in 100 Monte Carlo simulation runs and obtained a distribution for each moment of the LOSVD. The mean and standard deviation of the Monte Carlo distribution are taken as measurement and 1 uncertainty of the kinematics.
Since the Milky Way nuclear star cluster is at a distance of only 8 kpc, the spectroscopic observations are spatially semi-resolved. Bright stars can be resolved individually, and contribute a large fraction of the flux. For that reason we used the cleaned data cube of Feldmeier et al. (2014), where bright stars were excluded. However, the kinematic maps still show stochastic shot noise. As a consequence, the difference of the kinematics in adjacent bins can be higher than their uncertainties, which causes problems when we model the kinematics. The stochastic noise can be mistaken for signal, and this means the best fit will be achieved by modelling the shot noise. To prevent this, we increased our kinematic uncertainties such that the difference of the measurement in two adjacent bins () is less than the sum of their uncertainties (). We did this for the velocity , velocity dispersion , , and , and find that it is required for about 68 per cent of the kinematic data uncertainties. Additionally, we point-symmetrised the kinematics using the procedure of van den Bosch & de Zeeuw (2010). The median uncertainties or , , , and are 24.6 km s, 18.4 km s, 0.15, and 0.17. Our resulting kinematic maps are consistent with the maps of Feldmeier et al. (2014). We find rotation in the velocity map of approximately 50 km s and an increase in the velocity dispersion from about 65 km s towards =135 km s at the centre. The kinematic maps are shown on the top row of Fig. 1, the uncertainties are shown on the bottom row.
2.2 Imaging data and surface brightness distribution
The light distribution of the nuclear star cluster traces the stellar density. We require the two-dimensional light distribution of the red giant stars, which are our kinematic tracers. The extinction is high at optical wavelengths in the Galactic centre ( 30 mag, Scoville et al., 2003; Gao et al., 2013), therefore we used near- and mid-infrared images.
For the central 40.4 arcsec 40.4 arcsec (1.6 pc 1.6 pc) we used the high-resolution NACO -band mosaic of Schödel et al. (2009), which has a spatial scale of 0.027 arcsec pixel. We preferred the band over the band in order to avoid light from gas emission lines in the band (Br and He I, Paumard et al., 2004). Our kinematic tracers are cool late-type stars, but there are also more than 100 hot, young stars located in the centre of the cluster, within a projected radius = 0.5 pc (12.8 arcsec, Paumard et al. 2006). We masked out the young stars from the image with a 15 pixel radius. For the bright red supergiant IRS 7 we used a larger mask with a 30 pixel radius. Beyond the central 0.5 pc, the nuclear star cluster light is dominated by cool stars, and the contribution of young stars is negligible (Feldmeier-Krause et al., 2015).
For the large-scale light distribution, we used Spitzer IRAC images (Stolovy et al., 2006). These images were corrected for dust extinction and PAH emission by Schödel et al. (2014). We used the extinction and emission corrected 4.5 image to measure the light distribution. The image was smoothed to a scale of 5 arcsec pixel, and extends over 270 pc 200 pc. We excluded a central circle with = 0.6 pc (15.4 arcsec) to avoid contribution from ionised gas emission and young stars. In addition we masked out the young Quintuplet star cluster (Figer et al., 1999), and the dark 20-km s-cloud M-0.13-0.08 (García-Marín et al., 2011).
We used the MGE_FIT_SECTORS package (Cappellari, 2002) to derive the surface brightness distribution. The Multi-Gaussian-Expansion model (Emsellem et al., 1994) has the advantage that it can be deprojected analytically. We measured the photometry of the two images along the major axis and the minor axis. We assumed that the cluster’s major axis is aligned along the Galactic plane, as found by Schödel et al. (2014), and constant. The centre is the position of Sgr A*, which is the radio source associated with the Galactic centre supermassive black hole. We fitted a scale factor to match the photometry of the two images in the region where they overlap (1627.8 arcsec). Then we measured the photometry on each image along 12 angular sectors, and converted the NACO photometry to the Spitzer flux. Assuming four-fold symmetry, the measurements of four quadrants are averaged on elliptical annuli with constant ellipticity. Using the photometric measurements of the two images, we fitted a set of two-dimensional Gaussian functions, taking the point-spread-function (PSF) of the NACO image into account.
A comparison with the surface brightness profile of Fritz et al. (2016, their Fig. 2) showed that our profile is steeper in the central 30 arcsec. A possible reason is the small overlap region of the Spitzer and NACO images, and that the Spitzer flux could be too high at the centre. Maybe the PAH emission correction of the Spitzer image was too low. The mid-infrared dust emission is significant out to almost one arcmin. Fritz et al. (2016) used NACO - and -band images in the central = 20 arcsec. Out to 1 000 arcsec (39 pc) they used Hubble Space Telescope WFC3 data (M127 and M153 filters) and public VISTA Variables in the Via Lactea Survey images ( and bands, Saito et al. 2012). We lowered the intensities of the central Gaussians by scaling our averaged profile to the one-dimensional flux density profile of Fritz et al. (2016). As a result the amplitudes of the inner Gaussians become smaller, but the outer Gaussians (>100 arcsec 4 pc) are nearly unchanged. We list the components of the Multi-Gaussian Expansion in Table 1 and plot the surface brightness profile in Fig. 2 (upper panel). We also show the projected axial ratio as a function of radius in the lower panel of Fig. 2. Out to the central 1 pc, is increasing from 0.4 to 0.7. Schödel et al. (2014) found a mean axial ratio of 0.710.02 for the nuclear star cluster. This is in agreement with our maximum value of . However, decreases at larger radii, as the contribution from the nuclear stellar disc becomes more important and the light distribution therefore flatter.
We note that there are three main differences with the surface brightness distribution derived by Feldmeier et al. (2014): (1) We used an -band instead of a -band NACO image to avoid ionised gas emission; (2) We masked young stars in the NACO image to match the distribution of stars used as kinematic tracers; and (3) We scaled the central photometry to the flux density data of Fritz et al. (2016) to avoid a possible overestimation of the central flux when scaled to the Spitzer image. All three changes influence only the central part of the surface brightness distribution, as ionised gas emission and light from young stars are only important in the central parsec.
The surface brightness profile is deprojected to obtain the three-dimensional spatial stellar light distribution. In general, the deprojection of a surface brightness profile is non-unique (e.g. Rybicki, 1987; Franx, 1988), and our deprojection is only one possible solution. The MGE deprojection produces smooth intrinsic densities, which are in agreement with the photometric observations (Cappellari, 2002).
3 Dynamical Models of the Milky Way nuclear star cluster
3.1 Schwarzschild’s method
Orbit-based models or Schwarzschild models are a useful tool to model the dynamics of stellar systems by orbit superposition. The first step of Schwarzschild’s method is to integrate the equations of motion for a representative library of stellar orbits in a gravitational potential . Then the observables for each orbit are computed, considering projection, PSF convolution and pixel binning. The next step is to find orbital weights to combine the orbits such that they reproduce the observed data. Schwarzschild models are a powerful tool to recover the intrinsic kinematical structure and the underlying gravitational potential (Schwarzschild, 1979; van de Ven et al., 2008; van den Bosch & van de Ven, 2009). We refer the reader for further details to van den Bosch et al. (2008) for implementation and van de Ven et al. (2008) for verification of the triaxial Schwarzschild code.
We calculated orbits in the combined gravitational potential of a supermassive black hole and the star cluster , inferred from the imaging data. As we run triaxial models, there are three intrinsic shape parameters, , , and , for the cluster. The shape parameters characterise the axial ratios for the long, intermediate and short axes , , and . They are defined as , , and , where is the length of the longest axis projected on the sky. Thus, represents the compression of due to projection on the sky. Each set of axial ratios refers to a set of viewing angles (, see also van den Bosch et al. 2008). The surface brightness distribution is deprojected given the intrinsic shape parameters , and multiplied with the dynamical mass-to-light ratio to get the intrinsic stellar mass density . From Poisson’s equation one calculates the gravitational potential. We did this for different values of the black hole mass M, dynamical mass-to-light ratio , and different shape parameters. In total our model has five free parameters, M, , , , and .
Besides the considered stellar population and the supermassive black hole, there are other components within the nuclear star cluster, which we neglected. We measured a dynamical mass-to-light ratio, which combines the stellar mass-to-light ratio with other components. These components are the young stars, ionised gas, neutral gas, and dark matter. The young stars are at a distance of about 0.5 pc from the supermassive black hole. The lower limit of the total mass of young stars is 12 000 . However, the total enclosed extended mass in the same region is 10 (Oh et al., 2009; Feldmeier et al., 2014), and the mass of the supermassive black hole is 4 10 . The mass of the young stars is therefore probably negligible. The hot ionised gas has a mass of only a few 100 (Ferrière, 2012), and cannot influence the stellar dynamics significantly. The neutral gas in the circum-nuclear disc may contribute more mass, estimates range from 10 (Etxaluze et al., 2011; Requena-Torres et al., 2012) to 10 (Christopher et al., 2005), though this is probably the upper limit (Genzel et al., 2010). The circum-nuclear disc extends over a distance of about 1 pc to more than 5 pc from the centre. At 5 pc, the total enclosed mass is 10 (McGinn et al., 1989; Feldmeier et al., 2014). We decided to neglect the mass distribution of the circum-nuclear disc in our dynamical models, since it is very uncertain, and makes up only 0.1 to 10 per cent of the enclosed mass. The contribution of dark matter to the nuclear star cluster mass is also neglected. Linden (2014) showed that the fraction of dark matter in the central 100 pc of the Milky Way is about 6.6 per cent, assuming the traditional dark matter profile of Navarro et al. (1996).
The orbit library should be as general as possible and representative for the gravitational potential. We assumed that the orbits are regular and that three integrals of motion, , , and , are conserved. The orbit families consist of box orbits, which can cross the centre and have an average angular momentum of zero, and three types of tube orbits, which avoid the centre. The tube orbits are divided in short-axis-tube orbits, which have non-zero mean angular momentum around the short axis, outer and inner long-axis-tube orbits, which have non-zero mean angular momentum around the long axis. The orbit grid should sample the entire phase space. It has to be dense enough to suppress discreteness noise, but integration has to be done in a reasonable amount of computing time.
We followed van den Bosch et al. (2008) and sampled the orbit energy using a logarithmic grid in radius. Each energy is linked to the radius by calculating the potential at = (, 0, 0). We sampled = 35 energies calculated from in logarithmic steps ranging from = to = , i.e. 3.16 arcsec to 4.4 degree or 0.12 pc to 616.5 pc. We note that the outer radius is about 3.5 times the outermost Gaussian of the MGE fit. We tested lower values of the inner radius but found consistent results. For each energy, the starting point of an orbit was selected from a linear grid over 14 values each. For details on the orbit sampling we refer to van den Bosch et al. (2008). In total, we have = 6860 orbits. Each orbit was integrated over 200 periods, and sampled on 100 000 points per orbit. For each orbit we stored the intrinsic and projected properties. The projected orbits are stored in a grid, with PSF convolution and pixel size of the observed data taken into account. The velocities were stored in 183 bins between 7.4 and +7.4 . These numbers guarantee a proper sampling of the observed velocity profiles (Cretton et al., 1999).
Solving the orbital weight distribution
The model has to fit the kinematic data, the intrinsic and the projected mass distribution. The fit was done by finding a linear combination of the orbits, and solving for orbital weights . Each orbital weight corresponds to a mass on the respective orbit , and the weights are therefore non-negative. We used the non-negative least-squares (NNLS) algorithm of Lawson & Hanson (1974) which was also used by Rix et al. (1997), van der Marel et al. (1998), and Cretton et al. (1999). One of the fitting constraints is to make sure that the model is self-consistent. It is required that the orbit superposition reproduces the intrinsic and projected aperture masses within two per cent, which is the typical accuracy of the observed surface brightness (van den Bosch et al., 2008).
3.2 Constraining the input parameters
We ran 4899 models with different parameter combinations of M, , . The black hole mass M was sampled in logarithmic steps of 0.2 from 5.5 to 7.5, starting with 6.3 (i.e. M 2 ). The mass-to-light ratio was linearly sampled between 0.1 and 2.0 with steps of 0.04, with a starting value of 0.6 (in units of ). The starting model had () = (0.29, 0.84, 0.99). We sampled different combinations of () with a step size of (0.01, 0.02, 0.01), with the boundaries 0.05 < <0.29, 0.40 < <0.99, and 0.70 < 0.99. When fitting the orbital weights, we already applied the self-consistency criteria to make sure that the photometry is fitted with a high accuracy of at least two per cent for each model. The best fit of the five parameters M, , , however, was found by calculating the of the kinematic moments. After the starting model we computed different combinations of (M, , ) in the grid. We used a analysis to find regions with good fits and started new models, with a finer sampling at low . This was done in several iterations, and it was not necessary to compute each point in the (M, , ) grid.
We measured four kinematic moments in 175 spatial bins of the spectroscopic map in Sect. 2.1. In each bin we measured the velocity , velocity dispersion , and the higher Gauss-Hermite moments , and . Thus, the number of observables is the number of spatial bins times the number of kinematic moments, in our case N = 175 4 = 700. We minimised the function
where and denote the kinematic measurements, the respective measurement uncertainties, , , , and the model kinematics, to find the best-fitting model.
3.3 Modelling results
The best-fitting model
Our best-fitting parameters are M = , = 0.90, = 0.28, = 0.64, = 0.99. This corresponds to best-fitting viewing angles = 80°, = 79°, = 91°. We show the surface brightness map and the symmetrised kinematic maps in Fig. 3. The upper row are the data, the lower row are the maps of the best-fitting model. The misalignment of the kinematic rotation axis with respect to the photometry, and the perpendicular rotating substructure at 20 arcsec (0.8 pc) found by Feldmeier et al. (2014) are well reproduced in the model velocity map. We also show the data and the best-fitting model in Fig. 4, the panels denote surface brightness, , , and . The surface brightness map is reproduced within one per cent. The highest and lowest values of , , , and the lowest values of in the data have higher absolute values than in the best-fitting model, but are consistent within their uncertainties. The best fit has = 290. With M = 5 fitted parameters and N = 4175=700 observational constraints, this means = 0.42. That is less than one is partially due to the large uncertainties of the kinematics, and the fact that the kinematic measurements are correlated.
We illustrate the distribution of for the 4899 models in Fig. 5. We plot each combination of parameters. Red colours denote low , bluer, smaller symbols denote high . The black cross denotes the best-fitting model. The observed projected flattening of the surface brightness profile constrains the viewing angles and thus also the intrinsic shape parameters. In particular, a flat 1 means that the stellar system is observed along one of the principal planes (van den Bosch et al., 2008). We obtained as lowest value = 0.30, this limits the possible projections and puts an upper limit on the parameter . Likewise, the value of = 0.99 is the boundary value of the grid. The values of and denote the respective minimum values of the radially varying intrinsic axial ratios , and over the entire radial range of the photometry, i.e. the nuclear stellar disc and the embedded nuclear star cluster. The upper left panel of Fig. 5 shows that for each value of , the best-fitting is approximately 0.90. A similar behaviour is found with and . There is only a slight increase of the best-fitting with higher values of . At the same time, the best-fitting values of M do not show a strong dependence on or (second row) in the allowed parameter ranges. The intrinsic shape parameters do not influence our best fit for M, as this measurement is mostly made from the inner bins and the outer bins contribute little. The outer bins, however, contribute to the intrinsic shape fit. The supermassive black hole mass and the dynamical mass-to-light ratio are correlated. For higher values of , a lower M fits the data.
We show how depends on the different parameters in Fig. 6. The best-fitting model, which has the lowest , is marked as blue asterisk symbol. The blue lines denote the 1, 2, and 3 confidence limits, corresponding to = 5.9, 11.3, and 18.2. The red line illustrates the standard deviation of itself, i.e. = 37.3, where N = 700, and M = 5. This value was used as confidence limit by van den Bosch & van de Ven (2009). In Table 2 we list the 1 and 3 uncertainties.
Mass profile and intrinsic shape
We show the enclosed total mass as a function of the projected radius in Fig. 7, grey shaded contours are the 3 uncertainties. The mass was computed within spherical shells. We also plot the results of various other studies. Most studies assumed a spherical shape of the cluster, Feldmeier et al. (2014) and Chatzopoulos et al. (2015a) assumed axisymmetry. Some of the studies used also different Galactocentric distances, so we scaled the masses using = 8.0 kpc. Our results are in agreement with several studies in the central 100 arcsec. However, we obtain a lower enclosed mass than Trippe et al. (2008) and Oh et al. (2009), who used the Jeans equation for a spherical system to obtain the enclosed mass. At larger radii 400 arcsec (15.5 pc) beyond the reach of our kinematic data, we obtained a higher mass than Lindqvist et al. (1992), but we are in agreement at =750–1600 arcsec. Their data extend to larger radii, but their assumption of spherical symmetry does no longer hold at such large radii.
The enclosed total mass within a sphere with the radius = 8.4 pc, i.e. about two times the effective radius of the nuclear star cluster, is M = (2.10.7) 10 , here we give the 3 uncertainty.
The black hole influences the stellar kinematics only at the centre of the nuclear star cluster. Out to = 53 arcsec (2 pc), the best-fitting mass of the black hole (M = ) is higher than the enclosed stellar mass of our best-fitting model. Merritt (2004) defined the radius of influence of a black hole as the radius where the enclosed stellar mass equals two times the black hole mass. With this definition and a black hole mass of 4 , we obtain = (104) arcsec, i.e. approximately (4.0) pc. This result is in agreement with Alexander (2005), who found = 3 pc. The kinematic measurements at larger radii have little influence on the black hole mass measurement, but are important to constrain the orbital structure and dynamical mass-to-light ratio.
The shape of the nuclear star cluster is illustrated in Fig. 8. We show the intrinsic axial ratios and as a function of radius . The axial ratio is low in the centre (=0.3), increases to =0.8 at 35 arcsec, and then decreases to =0.6 at =150 arcsec. The best-fitting shape parameter, =0.28, is approximately the central value. The axial ratio is also low in the centre (0.65), increases to =0.9 at =40 arcsec, and then decreases to =0.75 at =150 arcsec. Though, the uncertainty of is rather high, and the decrease not significant. Also the best-fitting axial ratio, =0.64, is close to the central value. We plot the triaxiality , it varies between 0.45 and 0.7. The triaxiality increases from =70 arcsec to the outer radius =150 arcsec. This may be a signature of the increasing influence of the triaxial Galactic bulge at larger radii (Tsatsi et al., 2017). However, given the large uncertainty of , this increase is not significant. A constant or a decreasing are not excluded by the large uncertainties.
The best-fitting model has tangential anisotropy in the centre of the cluster. The value of the anisotropy = is negative, where is the tangential velocity dispersion and is the radial velocity dispersion. We show the anisotropy as a function of radius in Fig. 9, top panel. We plot the mean anisotropy of the models within the 1 uncertainty limit. The uncertainty of is about 0.1. The plot extends to the outer edge of the kinematic data at 150 arcsec. The cluster kinematics becomes nearly isotropic at radii >70 arcsec.
We show the angular momentum distribution of the orbits in Fig. 10. The colours denote the density of orbits passing radius with mean angular momentum (top panel) or (bottom panel). The plot of denotes rotation about the short z-axis. Orbits with 0 are contributed by short-axis-tube orbits, while long-axis-tube orbits have = 0. On the other hand, denotes rotation about the long x-axis (bottom panel), and orbits with 0 are contributed by long-axis-tube orbits. Short-axis-tube orbits have = 0. Long-axis-tube orbits are most important in the central 20–60 arcsec and at larger radii 80 arcsec. Short-axis-tube orbits, which contribute in total more mass than long-axis-tube orbits, are most important at = 60–140 arcsec. We illustrate the distribution of the stellar mass on the different orbit types also in Fig. 9 (bottom panel) as a function of radius. Most stars (>50 per cent) are on short-axis-tube orbits, i.e. they orbit the minor axis. Long-axis-tube orbits contribute about 40 per cent to the luminosity and thus stellar mass in the central 30 arcsec. They produce the perpendicular rotating substructure at 20 arcsec (0.8 pc) found by Feldmeier et al. (2014). At larger radii, long-axis-tube orbits contribute only about 30 per cent to the stellar mass. Box orbits contribute little mass in the centre (<10 per cent), but their relative mass increases towards larger radii. At = 150 arcsec (5.8 pc), they contribute 20 per cent.
4.1 Difference of the resulting black hole mass
The currently best black hole mass estimate is (4.02 0.20) , derived by Boehle et al. (2016) from Keplerian stellar orbits around the supermassive black hole. Using axisymmetric Jeans models and the same spectroscopic data as our study, Feldmeier et al. (2014) found a lower value of M = (1.7) . Our best fit using triaxial Schwarzschild models is ( ) . This measurement is consistent with the direct measurements of Boehle et al. (2016) within the 1 uncertainty limit. The result is also in agreement with the lower black hole mass of Feldmeier et al. (2014). We derived a 3 lower limit for the black hole of 0.7 , and an upper limit of 5.4 . We briefly discuss the model degeneracies, possible reasons for the different black hole mass measurements, and why our results are closer to the direct measurement than the black hole mass derived by Feldmeier et al. (2014).
Some model parameters seem to be correlated. This becomes clear when looking at Fig. 5. The best-fitting value of apparently increases with increasing dynamical mass-to-light ratio (second column of the first row). However, the value of has little effect on M, as can be seen in the second column of the second row in Fig. 5 . With a higher value of , the best-fitting M increases slightly. At larger , the -contours of M and broaden. This means that for a more oblate axisymmetric cluster with closer to one, M and are not as well constrained as with smaller values of . The black hole mass M=4.1 can be obtained with a higher value of =0.82 and =1.22, this model is within the 1 uncertainties.
The dynamical mass-to-light ratio is inversely correlated with the black hole mass (fourth column of the first row in Fig. 5). The higher , i.e. the more massive the cluster, the less massive is the black hole. This degeneracy is often obtained in dynamical models. Valluri et al. (2004) found that the degeneracy of M depends on how well the black hole’s sphere of influence is resolved. The measurement of is better constrained when the data extend to larger radii, provided that is constant over the entire field. We have several kinematic data bins within the radius of influence of the supermassive black hole, and our data extend to one effective radius. This may not be sufficient to put strong constraints on . To get agreement with the measurement of (4.02 0.20) (Boehle et al., 2016), we would require a lower value of 0.75, which is within the 1 uncertainties.
Influence of the surface brightness profile
The shape of the surface brightness profile is important to estimate the mass of the supermassive black hole. The surface brightness profile has to represent the density of the kinematic tracer. We excluded young stars and ionised gas from the surface brightness profile, as these components contribute little mass compared to the cool, old stars we used as kinematic tracers. Excluding these components results in a lower surface brightness and stellar mass in the centre compared to Feldmeier et al. (2014). The stellar mass we obtain at = 48 arcsec (1.8 pc) is 1.3 less. Our black hole mass is therefore higher, and closer to the direct measurement of M 4 . We ran the same axisymmetric models (Cappellari, 2008), using the same kinematic data as Feldmeier et al. (2014), but our surface brightness distribution from Table 1. The best fit is obtained with M = () , = 0.89 , and a constant anisotropy of = 0.3. This result is in agreement with the triaxial Schwarzschild models, and confirms that the surface brightness profile has a strong influence on the results of the black hole mass and dynamical mass-to-light ratio.
The deprojection of the surface brightness profile is non-unique and only one possible solution. We applied the MGE method, which produces smooth intrinsic densities (Cappellari, 2002). We assumed that the deprojected density is smooth, as this is what we observe in other galaxies, where the clusters are observed from various viewing angles.
van den Bosch (1997) studied how much central density can be hidden in an oblate axisymmetric galaxy without effects on the projected surface brightness. They found that the percentage of hidden density in the centre of a Staeckel potential is 0 per cent for inclination =90, and 10 per cent for 80 of the total galaxy mass. The percentage of hidden density can increase in steeper cusps. However, the effect on the dynamics is still negligible. For our models, a hidden density is equivalent to a spatially varying mass-to-light ratio.
Spatially varying mass-to-light ratio
We assumed a constant dynamical mass-to-light ratio for the Schwarzschild models. We obtained = 0.90 (1 uncertainty). The dynamical mass-to-light ratio combines the stellar mass-to-light ratio with other components, it is sensitive to the presence of gas or dark matter.
Our best-fitting value of = 0.9 is consistent with stellar-population studies. Norris et al. (2014) found =0.8–1.2 at 4.6 for stellar populations with ages >7.0 Gyr and metallicities  ranging from -1.0 dex to +0.3 dex. For populations with lower metallicity =–2.18 dex and 13 Gyr age, can be even 1.5. However, for younger stellar populations 5 Gyr, Norris et al. (2014) obtained lower values 0.6. Our measurement of is averaged over the entire field of the kinematic data. We cannot exclude that the stellar age or metallicity changes over the range of the kinematic data. Stellar population studies of the red giant population were so far confined to the central 1–2 pc (e.g. Pfuhl et al., 2011; Do et al., 2015; Feldmeier-Krause et al., 2017). The stars in this region are mostly older than 5 Gyr and rather metal-rich. Our knowledge of the stellar population at the outer region of our field is based on only a few bright stars (e.g. Blum et al., 2003; Feldmeier et al., 2014). But these stars are brighter and probably younger than our kinematic tracer stars. However, the mass-to-light ratio for old stars in the mid-infrared varies modestly with age and metallicity in comparison to the optical mass-to-light ratio (Meidt et al., 2014). Therefore we do not expect a change of by more than 0.4 within the cluster. Should vary with radius, our mass profile (Fig. 7) would have a different shape. For example, if was lower in the centre than outside, this would increase M, and there would be less mass in the stellar component.
However, the stellar mass-to-light ratio may also increase towards the central =0.5 pc, as massive stellar remnants may migrate to the centre. The mass and distribution of dark stellar remnants, i.e. stellar-mass black holes and neutron stars, in the central parsec of the nuclear star cluster is uncertain. For a top-heavy initial mass function, there could be >1 in dark remnants (Morris, 1993), though Löckmann et al. (2010) found a lower mass of about 1 10 for a canonical initial mass function.
In our models we neglected the mass of molecular gas in the circum-nuclear disc. The molecular gas may contribute 10 . The gas disc extends from 17 pc along the Galactic plane, but only to 3 pc along the minor axis (Ferrière, 2012). Thus, the molecular gas is located in the central part of our spectroscopic field, but absent in the North. If the gas contributes significantly to the cluster mass, our assumption of spatially constant would be violated, and would be higher than for a stellar component alone. When we assume the maximum gas mass of , the value of a constant decreases to about 0.85, which is within our 1 uncertainty limit.
The spatial distribution of dark matter in the Galactic centre is uncertain. A classical cuspy Navarro et al. (1996) dark matter profile results in a dark matter fraction of about 6.6 per cent in the central 100 pc (Linden, 2014). However, black hole accretion, dark matter annihilation, and scattering alter the shape of the dark matter distribution in the Galactic centre. Vasiliev & Zelnikov (2008) found that these effects produce a shallower dark matter profile in the central 2 pc than further out. The dark matter mass inferred from the classical cusp is reduced by up to 50 per cent in the central 2 pc. The contribution of dark matter to the nuclear star cluster mass should therefore be negligible. Although the dark matter distribution may be different from the luminous baryonic matter, and the dynamical mass-to-light ratio for that reason not spatially constant, the effect on the cluster mass distribution should be only minor.
4.2 Triaxial cluster shape
Our best-fitting model has axial ratios of = = 0.28 , = = 0.64 , and = = 0.99 . These axial ratios correspond to viewing angles = 80°, = 79°, and = 91°. The angle denotes the polar viewing angle, the azimuthal viewing angle, and is the misalignment angle between photometric major axis and the projected intrinsic long axis (van den Bosch et al., 2008; van den Bosch & van de Ven, 2009). For the best-fitting model the angle between the cluster’s major axis and the line-of-sight is about 79°. The cluster’s shape is illustrated in Fig. 8, the triaxiality paramter varies between 0.45 and 0.7, with an average at 0.6. An oblate axisymmetric system has = 0, a prolate axisymmetric system has = 1.
Also the Milky Way’s bulge is triaxial, the axial ratios are = 0.26 and = 0.63 (Wegg & Gerhard, 2013). The shape was derived from the density of red clump stars in the central 800 pc of the bulge. The Milky Way bulge is much larger than the nuclear star cluster, and extends out to about 2.5 kpc. At the outer edge of our data, at =150 arcsec, we obtain for the intrinsic shape parameters =0.60 and =0.75. Both and are higher than found in the bulge. However, they are also decreasing, though the decrease of is not significant. The bulge has a peanut or X-shape (Nataf et al., 2010; McWilliam & Zoccali, 2010). The angle between the bulge major axis and the line-of-sight to the Galactic centre is about 27° (Rattenbury et al., 2007; Wegg & Gerhard, 2013), while we obtained 79° for the nuclear star cluster. There are also indications for another bar within the Galactic bulge from star count data in the inner 140 pc (Alard, 2001) or even 560 pc (Nishiyama et al., 2005). Rodriguez-Fernandez & Combes (2008) derived an angle 60-75° for a thick triaxial nuclear bar with axial ratios 0.55 and 0.75.
One possible scenario for nuclear star cluster formation is that massive star clusters (10–10 ) formed in the galactic disc, migrated to the galaxy’s centre and merged (Neumayer et al., 2011; Guillard et al., 2016). Simulations of multiple star cluster mergers and of star cluster accretion on a nuclear stellar component can produce triaxial nuclear star clusters (Bekki et al., 2004; Hartmann et al., 2011; Perets & Mastrobuono-Battisti, 2014). However, so far no observational study was able to constrain the triaxial shape of nuclear star clusters in general. Hartmann et al. (2011) constrained the shape of two nuclear star clusters and found agreement with an axisymmetric shape. For the Milky Way nuclear star cluster, may be as high as 0.94 within the 3 uncertainties. Thus, a nearly axisymmetric shape is consistent with our data.
4.3 Caveats and Considerations
Regime of semi-resolved populations
We used integrated light spectroscopy to measure the stellar kinematics. This is the common approach for extragalactic systems, which have a distance of several Mpc. The measured kinematics are weighted by the respective luminosities of different stars. As the Milky Way nuclear star cluster is only 8 kpc distant, we are in the regime of semi-resolved populations. The brightest stars can be resolved individually, and these stars contribute a large fraction of the flux. In consequence, individual spatial bins can be dominated by a single star. Instead of measuring the spectrum of an ensemble of stars, one measures a spectrum in which a large percentage of the flux is contributed by one single star. This causes shot noise, and high differences between neighbouring spatial bins. We accounted for this problem by excluding the brightest stars from the spectroscopic map. This method helps to significantly reduce the intrinsic scatter of the velocity dispersion (Lützgendorf et al., 2011; Bianchini et al., 2015). We further increased the kinematic uncertainties such that the data in two neighbouring bins have consistent values within their uncertainties. This helps to prevent that the models fit only stochastic shot noise. Due to the large kinematic uncertainties, the intrinsic shape parameters , , , and the dynamical mass-to-light ratio are not very well constrained, and have large error margins.
At a distance of only 8 kpc, also the relative distances of the stars become more important. A star located on the near side of the nuclear star cluster, at a distance d = 7.9 kpc, contributes 1.05 times more flux than a star with the same absolute magnitude at the far side of the cluster, at d = 8.1 kpc. In an extragalactic system, the distance of a star at the near side and the distance of a star at the far side with respect to the observer are approximately the same, as the system is farther away. For a galaxy at d = 5 Mpc, a relative difference of 200 pc changes the flux only by a factor 1.00008. Even foreground stars that belong to the outer parts of the stellar system contribute roughly the same flux as a star with the same magnitude that is located in the galactic nucleus.
Another observational complication is interstellar extinction in the Galactic centre, which varies on arcsecond scales (Schödel et al., 2010). In particular, the field of view of the kinematic data contains the so-called 20-km s-cloud (M-0.13-0.08, e.g. García-Marín et al., 2011) in the Galactic southwest. It lies at a projected distance of about 70 arcsec (3 pc) from the centre, and probably about 5 pc in front of Sgr A* (Ferrière, 2012). This cloud blocks the light from stars of the nuclear star cluster. We cannot access the kinematics of stars behind this cloud. There is also interstellar dust within a projected distance of 20 arcsec (0.8 pc) from the centre, i.e. within the radius of influence of the black hole. This dust causes extinction within the nuclear star cluster by up to 0.8 mag (Chatzopoulos et al., 2015b). As a consequence, the two effects of dimming by distance and by extinction add up and stars that lie on the far side of the nuclear star cluster appear even more faint than the stars on the near side.
Both the semi-resolved stellar population and the inter-cluster extinction cause that our observations are biased to the near side of the nuclear star cluster. As a consequence, we measured a lower limit of the velocity dispersion. Feldmeier et al. (2014) found that the velocity dispersion in the projected radial range 6 arcsec < < 20 arcsec is smaller compared to the velocity dispersion computed from proper motion data of Schödel et al. (2009), which is based on resolved stars. For resolved stars, the velocity dispersion is not weighted by the flux of the stars. An underestimated velocity dispersion means that the black hole mass measurement is biased to lower values.
This observational bias also influences the measurements of , and . In particular, the cluster may appear compressed along the line-of-sight, and thus the value of = 0.64 may be too low. As a consequence, = 0.90 would be underestimated (see second column of the first row in Fig. 5).
Influence of figure rotation
The Galaxy rotates, and with it the nuclear star cluster. In a non-axisymmetric, rotating system, centrifugal and Coriolis forces play a role. However, figure rotation and the resulting forces were not included in our triaxial models. Figure rotation influences the stellar orbits. The prograde and retrograde tube orbits no longer fill the same volumes, while the box orbits acquire net mean angular momentum (e.g. Heisler et al., 1982; Schwarzschild, 1982; Sellwood & Wilkinson, 1993; Skokos et al., 2002). As a result, orbit-based, tumbling, triaxial models are computationally expensive. Other than an early attempt by Zhao (1996), no such models have been constructed that take into account kinematic data. It is difficult to predict how our results would change in a rotating model. The inferred orbital structure will be affected (depending on the tumbling speed of the nuclear star cluster), but our results on the mass distribution are likely to be fairly robust, as the assumption of a constant mass-to-light ratio is probably more important.
5 Summary and Outlook
We constructed for the first time triaxial orbit-based Schwarzschild models of the Milky Way nuclear star cluster. We used the spectroscopic integrated light maps of Feldmeier et al. (2014) to measure the cluster kinematics of the central 60 pc of the Milky Way. As photometry we used Spitzer 4.5 and NACO -band images, and measured a two-dimensional surface brightness distribution. We excluded young stars, avoided gas emission and dark clouds in the photometric data. Our triaxial models were based on the code by van den Bosch et al. (2008). Our best-fitting model contains a black hole of mass M = () , a dynamical mass-to-light ratio of = (0.90) /, and shape parameters = 0.28, = 0.64, and = 0.99. Our black hole mass measurement is in agreement with the direct measurement of (4.02 0.20) by Boehle et al. (2016). We obtain a total cluster mass M = (2.10.7) 10 within a spherical shell with radius = 2 = 8.4 pc. The best-fitting model is tangentially anisotropic in the central = 0.5-2 pc of the nuclear star cluster, but close to isotropic at larger radii. The model is able to recover the long-axis rotation in the central = 0.8 pc found by Feldmeier et al. (2014), and the misalignment of the kinematic rotation axis from the photometric minor axis.
There are several possible ways to extend the dynamical models in the future. One way is to include a component for the neutral gas disc inside the nuclear star cluster. If the gas mass is close to the upper limit of 10 , the dynamical mass-to-light ratio would probably decrease slightly, and in return would slightly increase the black hole mass. Modelling a spatially varying mass-to-light ratio may provide a better representation of the cluster’s intrinsic properties. Further, proper motions can be included in combination with discrete line-of-sight velocities, as shown by van de Ven et al. (2006) and van den Bosch et al. (2006) for axisymmetric Schwarzschild models. Watkins et al. (2013) extended axisymmetric Jeans models and implemented discrete kinematic data without binning. Using discrete data means that the stars are not weighted by their luminosities. This prevents the previously discussed bias towards the near side of the cluster.
GvdV acknowledges partial support from Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject A7 and A8) funded by the German Research Foundation. RS acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n. 614922. We thank Remco van den Bosch for helpful discussions about the project. We finally thank the anonymous reviewer for useful comments and suggestions.
- thanks: Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere, Chile (289.B-5010(A)).
- pubyear: 2016
- pagerange: Triaxial orbit-based modelling of the Milky Way Nuclear Star Cluster
thanks: Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere, Chile (289.B-5010(A)).–Triaxial orbit-based modelling of the Milky Way Nuclear Star Cluster thanks: Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere, Chile (289.B-5010(A)).
- Alard C., 2001, \aap, 379, L44
- Alexander T., 2005, \physrep, 419, 65
- Bahcall J. N., Tremaine S., 1981, \apj, 244, 805
- Bekki K., Couch W. J., Drinkwater M. J., Shioya Y., 2004, \apjl, 610, L13
- Bianchini P., Norris M. A., van de Ven G., Schinnerer E., 2015, \mnras, 453, 365
- Blum R. D., Ramírez S. V., Sellgren K., Olsen K., 2003, \apj, 597, 323
- Boehle A., Ghez A. M., Schödel R., Meyer L., Yelda S., Albers S., Martinez G. D., Becklin E. E. et al, 2016, \apj, 830, 17
- Cappellari M., 2002, \mnras, 333, 400
- —, 2008, \mnras, 390, 71
- Cappellari M., Emsellem E., 2004, \pasp, 116, 138
- Chatzopoulos S., Fritz T. K., Gerhard O., Gillessen S., Wegg C., Genzel R., Pfuhl O., 2015a, \mnras, 447, 948
- Chatzopoulos S., Gerhard O., Fritz T. K., Wegg C., Gillessen S., Pfuhl O., Eisenhauer F., 2015b, \mnras, 453, 939
- Christopher M. H., Scoville N. Z., Stolovy S. R., Yun M. S., 2005, \apj, 622, 346
- Cretton N., de Zeeuw P. T., van der Marel R. P., Rix H.-W., 1999, \apjs, 124, 383
- Deguchi S., Imai H., Fujii T., Glass I. S., Ita Y., Izumiura H., Kameya O., Miyazaki A. et al, 2004, \pasj, 56, 261
- Do T., Kerzendorf W., Winsor N., Støstad M., Morris M. R., Lu J. R., Ghez A. M., 2015, \apj, 809, 143
- Do T., Martinez G. D., Yelda S., Ghez A., Bullock J., Kaplinghat M., Lu J. R., Peter A. H. G. et al, 2013, \apjl, 779, L6
- Drehmer D. A., Storchi-Bergmann T., Ferrari F., Cappellari M., Riffel R. A., 2015, \mnras, 450, 128
- Emsellem E., Monnet G., Bacon R., 1994, \aap, 285, 723
- Etxaluze M., Smith H. A., Tolls V., Stark A. A., González-Alfonso E., 2011, \aj, 142, 134
- Feldmeier A., Neumayer N., Seth A., Schödel R., Lützgendorf N., de Zeeuw P. T., Kissler-Patig M., Nishiyama S. et al, 2014, \aap, 570, A2
- Feldmeier-Krause A., Kerzendorf W., Neumayer N., Schödel R., Nogueras-Lara F., Do T., de Zeeuw P. T., Kuntschner H., 2017, \mnras, 464, 194
- Feldmeier-Krause A., Neumayer N., Schödel R., Seth A., Hilker M., de Zeeuw P. T., Kuntschner H., Walcher C. J. et al, 2015, \aap, 584, A2
- Ferrière K., 2012, \aap, 540, A50
- Figer D. F., McLean I. S., Morris M., 1999, \apj, 514, 202
- Franx M., 1988, \mnras, 231, 285
- Fritz T. K., Chatzopoulos S., Gerhard O., Gillessen S., Genzel R., Pfuhl O., Tacchella S., Eisenhauer F. et al, 2016, \apj, 821, 44
- Gao J., Li A., Jiang B. W., 2013, Earth, Planets, and Space, 65, 1127
- García-Marín M., Eckart A., Weiss A., Witzel G., Bremer M., Zamaninasab M., Morris M. R., Schödel R. et al, 2011, \apj, 738, 158
- Gebhardt K., Richstone D., Kormendy J., Lauer T. R., Ajhar E. A., Bender R., Dressler A., Faber S. M. et al, 2000, \aj, 119, 1157
- Genzel R., Eisenhauer F., Gillessen S., 2010, Reviews of Modern Physics, 82, 3121
- Genzel R., Thatte N., Krabbe A., Kroker H., Tacconi-Garman L. E., 1996, \apj, 472, 153
- Ghez A. M., Salim S., Weinberg N. N., Lu J. R., Do T., Dunn J. K., Matthews K., Morris M. R. et al, 2008, \apj, 689, 1044
- Gillessen S., Eisenhauer F., Trippe S., Alexander T., Genzel R., Martins F., Ott T., 2009, \apj, 692, 1075
- Guillard N., Emsellem E., Renaud F., 2016, \mnras, 461, 3620
- Haller J. W., Rieke M. J., Rieke G. H., Tamblyn P., Close L., Melia F., 1996, \apj, 456, 194
- Hartmann M., Debattista V. P., Seth A., Cappellari M., Quinn T. R., 2011, \mnras, 418, 2697
- Heisler J., Merritt D., Schwarzschild M., 1982, \apj, 258, 490
- Jeans J. H., 1922, \mnras, 82, 122
- Lawson C. L., Hanson R. J., 1974, Solving least squares problems
- Linden T., 2014, in IAU Symposium, Vol. 303, The Galactic Center: Feeding and Feedback in a Normal Galactic Nucleus, Sjouwerman L. O., Lang C. C., Ott J., eds., pp. 403–413
- Lindqvist M., Habing H. J., Winnberg A., 1992, \aap, 259, 118
- Löckmann U., Baumgardt H., Kroupa P., 2010, \mnras, 402, 519
- Lützgendorf N., Kissler-Patig M., Noyola E., Jalali B., de Zeeuw P. T., Gebhardt K., Baumgardt H., 2011, \aap, 533, A36
- Malkin Z., 2012, ArXiv e-prints
- McGinn M. T., Sellgren K., Becklin E. E., Hall D. N. B., 1989, \apj, 338, 824
- McWilliam A., Zoccali M., 2010, \apj, 724, 1491
- Meidt S. E., Schinnerer E., van de Ven G., Zaritsky D., Peletier R., Knapen J. H., Sheth K., Regan M. et al, 2014, \apj, 788, 144
- Merritt D., 2004, Coevolution of Black Holes and Galaxies, 263
- Morris M., 1993, \apj, 408, 496
- Nataf D. M., Udalski A., Gould A., Fouqué P., Stanek K. Z., 2010, \apjl, 721, L28
- Navarro J. F., Frenk C. S., White S. D. M., 1996, \apj, 462, 563
- Neumayer N., Walcher C. J., Andersen D., Sánchez S. F., Böker T., Rix H.-W., 2011, \mnras, 413, 1875
- Nishiyama S., Nagata T., Baba D., Haba Y., Kadowaki R., Kato D., Kurita M., Nagashima C. et al, 2005, \apjl, 621, L105
- Norris M. A., Meidt S., Van de Ven G., Schinnerer E., Groves B., Querejeta M., 2014, \apj, 797, 55
- Oh S., Kim S. S., Figer D. F., 2009, Journal of Korean Astronomical Society, 42, 17
- Oort J. H., 1977, \araa, 15, 295
- Paumard T., Genzel R., Martins F., Nayakshin S., Beloborodov A. M., Levin Y., Trippe S., Eisenhauer F. et al, 2006, \apj, 643, 1011
- Paumard T., Maillard J.-P., Morris M., 2004, \aap, 426, 81
- Perets H. B., Mastrobuono-Battisti A., 2014, \apjl, 784, L44
- Pfuhl O., Fritz T. K., Zilka M., Maness H., Eisenhauer F., Genzel R., Gillessen S., Ott T. et al, 2011, \apj, 741, 108
- Rattenbury N. J., Mao S., Sumi T., Smith M. C., 2007, \mnras, 378, 1064
- Reid M. J., Brunthaler A., 2004, \apj, 616, 872
- Requena-Torres M. A., Güsten R., Weiß A., Harris A. I., Martín-Pintado J., Stutzki J., Klein B., Heyminck S. et al, 2012, \aap, 542, L21
- Rieke G. H., Rieke M. J., 1988, \apjl, 330, L33
- Rix H.-W., de Zeeuw P. T., Cretton N., van der Marel R. P., Carollo C. M., 1997, \apj, 488, 702
- Rodriguez-Fernandez N. J., Combes F., 2008, \aap, 489, 115
- Rybicki G. B., 1987, in IAU Symposium, Vol. 127, Structure and Dynamics of Elliptical Galaxies, de Zeeuw P. T., ed., p. 397
- Saito R. K., Hempel M., Minniti D., Lucas P. W., Rejkuba M., Toledo I., Gonzalez O. A., Alonso-García J. et al, 2012, \aap, 537, A107
- Schödel R., Feldmeier A., Kunneriath D., Stolovy S., Neumayer N., Amaro-Seoane P., Nishiyama S., 2014, \aap, 566, A47
- Schödel R., Merritt D., Eckart A., 2009, \aap, 502, 91
- Schödel R., Najarro F., Muzic K., Eckart A., 2010, \aap, 511, A18
- Schwarzschild M., 1979, \apj, 232, 236
- Schwarzschild M., 1982, \apj, 263, 599
- Scoville N. Z., Stolovy S. R., Rieke M., Christopher M., Yusef-Zadeh F., 2003, \apj, 594, 294
- Sellgren K., McGinn M. T., Becklin E. E., Hall D. N., 1990, \apj, 359, 112
- Sellwood J. A., Wilkinson A., 1993, Reports on Progress in Physics, 56, 173
- Siopis C., Gebhardt K., Lauer T. R., Kormendy J., Pinkney J., Richstone D., Faber S. M., Tremaine S. et al, 2009, \apj, 693, 946
- Skokos C., Patsis P. A., Athanassoula E., 2002, \mnras, 333, 847
- Stolovy S., Ramirez S., Arendt R. G., Cotera A., Yusef-Zadeh F., Law C., Gezari D., Sellgren K. et al, 2006, Journal of Physics Conference Series, 54, 176
- Trippe S., Gillessen S., Gerhard O. E., Bartko H., Fritz T. K., Maness H. L., Eisenhauer F., Martins F. et al, 2008, \aap, 492, 419
- Tsatsi A., Mastrobuono-Battisti A., van de Ven G., Perets H. B., Bianchini P., Neumayer N., 2017, \mnras, 464, 3720
- Valluri M., Ferrarese L., Merritt D., Joseph C. L., 2005, \apj, 628, 137
- Valluri M., Merritt D., Emsellem E., 2004, \apj, 602, 66
- van de Ven G., de Zeeuw P. T., van den Bosch R. C. E., 2008, \mnras, 385, 614
- van de Ven G., van den Bosch R. C. E., Verolme E. K., de Zeeuw P. T., 2006, \aap, 445, 513
- van den Bosch F. C., 1997, \mnras, 287, 543
- van den Bosch R., de Zeeuw T., Gebhardt K., Noyola E., van de Ven G., 2006, \apj, 641, 852
- van den Bosch R. C. E., de Zeeuw P. T., 2010, \mnras, 401, 1770
- van den Bosch R. C. E., Greene J. E., Braatz J. A., Constantin A., Kuo C.-Y., 2016, \apj, 819, 11
- van den Bosch R. C. E., van de Ven G., 2009, \mnras, 398, 1117
- van den Bosch R. C. E., van de Ven G., Verolme E. K., Cappellari M., de Zeeuw P. T., 2008, \mnras, 385, 647
- van der Marel R. P., Cretton N., de Zeeuw P. T., Rix H.-W., 1998, \apj, 493, 613
- Vasiliev E., Zelnikov M., 2008, \prd, 78, 083506
- Wallace L., Hinkle K., 1996, \apjs, 107, 312
- Watkins L. L., van de Ven G., den Brok M., van den Bosch R. C. E., 2013, \mnras, 436, 2598
- Wegg C., Gerhard O., 2013, \mnras, 435, 1874
- Zhao H., 1996, \mnras, 283, 149