TY - JOUR
AU - Tenjes,, Peeter
AB - ABSTRACT We present a method to calculate gravitational potential gradients within regions containing few 10s of thousands stars with known phase space coordinates. The central idea of the method is to calculate orbital arcs for each star within a given region for a certain parametrized potential (gravitational acceleration) and to assume that position of each star on its orbital arc is a random variable with a uniform probability density in time. Thereafter, by combining individual probability densities of stars it is possible to calculate the overall probability density distribution and likelihood for a given region as a function of gravitational acceleration parameters. The likelihood has a maximum if the calculated probability distribution and the observed distribution of stars in phase space are consistent. This allows us to constrain gravitational accelerations and potential gradient values. The method assumes that phases of stars are mixed within the regions where stellar orbits are calculated. We tested the method for 12 small rectangular regions within simulated disc galaxy from Gaia Wiki. Tests show that even with a rather simple acceleration form the calculated accelerations in galactic plane coincide with their true values from simulation about 5 per cent, misalignment between the calculated and true acceleration vector directions is less than 1 deg (median values). The model can be used with the Milky Way Gaia complete solution data. methods: data analysis, stars: kinematics and dynamics, Galaxy: kinematics and dynamics 1 INTRODUCTION One of the most important functions determining kinematics and dynamics of galaxies is the gradient of the gravitational potential. When calculating gravitational potential and/or their derivatives, it is usually assumed that a galaxy is in a stationary state and its gravitational potential has some symmetry. Starting from the surface-brightness distribution a galaxy is assumed to consist of one or several axisymmetric or spherical components, the luminosity density distribution of which being described with a sufficiently simple analytical formula. Then, assuming that a similar formula is valid also for the mass density distribution it is possible to calculate gravitational potential derivatives by solving, for example, Jeans equations (Binney, Davies & Illingworth 1990; van der Marel, Binney & Davies 1990; Emsellem, Monnet & Bacon 1994; Cappellari 2008). Free parameters in the density distribution formulae can be derived by comparing calculated model predictions with the spectral and kinematic data of the galaxy. These kind of models can be rather simple (e.g. two-component models) and are thus suitable to model a large number of galaxies (see McGaugh, Lelli & Schombert 2016; Kalinova et al. 2017; Barone et al. 2018; Li et al. 2018). More detailed but also more time-consuming methods are models using Schwarzschild orbit based method (Schwarzschild 1979; Thomas et al. 2004; Valluri, Merritt & Emsellem 2004; Cappellari et al. 2006, 2007; van de Ven, de Zeeuw & van den Bosch 2008; Kowalczyk, Łokas & Valluri 2017) or action-angle torus-mapping technique (McGill & Binney 1990; Copin, Zhao & de Zeeuw 2000; Bovy 2014; Binney & McMillan 2015, 2016). These models can be rather flexible, that is, with a sufficient number of components it is possible to describe a galaxy with nearly arbitrary density distributions. These models take into account three integrals of motion and allow to describe galactic components with three-axial ellipsoidal density distributions and even bars and spiral arms (Binney 2018). Both methods assume global stationarity of a galaxy as complete orbits are calculated there. Unfortunately, corresponding calculations are more time consuming and the number of galaxies modelled with these methods is limited.1 Modelling of individual galaxies via N-body simulations (e.g. made to measure) is usually free of symmetry and stationarity assumptions but these models are even more time consuming and their spatial and mass resolution is yet moderate (Syer & Tremaine 1996; de Lorenzi et al. 2007; Long & Mao 2010; Zhu et al. 2014). Several galactic main components (stellar disc, bulge, stellar halo) often seem to be symmetric. But knowing the importance of interactions in galaxy evolution it is not always clear how good symmetry assumptions really are. As an example, for our Milky Way (MW) galaxy, it is not clear how important is the central bar while modelling the overall gravitational potential (Binney 2018). An excellent opportunity to study gravitational potential of the MW is given by the Gaia data (Gaia Collaboration et al. 2018). To calculate accelerations directly for individual stars from Gaia data requires very precise full six-dimensional phase space information, which is not available in Gaia. However, there is a possibility to calculate acceleration components from Gaia data within a sufficiently small regions containing many (several thousands) stars. In this paper, we propose a method to do this and we test the method by using simulated data. Application of the proposed method assumes that Gaia complete solution data are available, i.e one has sufficiently precise three-dimensional coordinates and velocities for a large number of stars. By selecting a region around a point in a galaxy containing several 10s of thousands stars and calculating orbits of these stars it is possible to calculate gravitational potential derivatives at a given point rather precisely. This paper continues our previous work (Kipper, Tempel & Tenjes 2018) where only Solar neighbourhood in one dimension was studied. With Gaia complete solution data the entire MW covered by sufficiently precise Gaia data can be analysed. The structure of the paper is following. In Section 2, we describe our proposed method to calculate the gravitational acceleration within discrete localized regions. In Section 3, we apply the method to simulated data taken from Gaia Challenge project. And finally, in Section 4 and Section 5 we present our results and conclude the paper, accordingly. 2 METHOD In this section, we present a method allowing to calculate gravitational acceleration vector components in a small bounded region D that contains sufficient number of stars with known phase space coordinates. These individual regions should be taken sufficiently small in order to make an assumption that gravitational potential derivatives within these region can be approximated with a simple analytical form. No assumptions about symmetries of the overall mass distribution are necessary. The whole gravitational acceleration field of a galaxy can be calculated by selecting many regions of this kind all over the galaxy. We assume that phases of stars are mixed only within these individual regions D, that is, they are locally mixed. This is different from the assumption of completely mixed phases. Let us consider first an orbit (a set of similar valued integrals of motion) extending all over the galaxy and containing sufficiently large number of stars. Let us describe positions of a star within this orbit as phases τi on the orbit. If the phases are completely mixed the distribution of τi is uniform, or equivalently, the phase distances between stars on the orbit are taken from exponential distribution with the mean value of Δτ. If the phases are not completely mixed but only moderately mixed (corresponding to a somewhat non-stationary galaxy) then Δτ slightly changes along the orbit. In this case we may introduce a concept of locally mixed phases, i.e. globally Δτ values can be different, but within a sufficiently small arc of the orbit Δτ can be approximated as Δτ ≃ const. For different arcs (or different regions D) this constant is different. Thus, the method allows that for a whole galaxy the phases of stars in a given orbit are not completely mixed (e.g. near to a bar). However, by selecting a sufficiently small region D, one can approximate the phases to be mixed locally. Knowing the phase space coordinates of a star i in a particular moment t0 (xi(t0), |$v$|i(t0)) and accelerations a(xi) caused by an external potential it is possible to integrate the equations of motion for the star i within a region D \begin{eqnarray*} \frac{\mathrm{d}{{\boldsymbol x}}_i}{\mathrm{d}t} = {\boldsymbol v}_i(t) \end{eqnarray*} (1) \begin{eqnarray*} \frac{\text{d}{\boldsymbol v}_i}{\mathrm{d}t} = {\bf {\boldsymbol a}} ({{\boldsymbol x}}_i). \end{eqnarray*} (2) Here t is time and also a parameter that determines the position of a star in an orbit. Initial conditions xi(t0) and |$v$|i(t0) can be taken from Gaia observations. For a given acceleration a(xi), a solution (xi(t), |$v$|i(t)) can be derived. Next we use q as a shorter designation of all six phase space coordinates. For all stars orbits are calculated within a region D and we assume that within this region acceleration components in equation (2) are only functions of coordinates or can be approximated within D as such. Let |$p_i ({{\boldsymbol q}}_i(t))\, \text{d}t$| be the probability to find a star i at its orbital point qi(t) within an interval determined by t and t + dt. For all stars these probabilities can be calculated by solving equations of motion (1) and (2), and then calculating how much time a star spends at these phase space coordinate intervals. The probability to find a star near q is an integral of corresponding probability density over the time star spends at corresponding orbital arc ∫pi(qi) dt. If we consider the whole region D this probability should be normalized to one \begin{eqnarray*} \int \limits _{{\bf q}_i\in D} p_i({{\boldsymbol q}}_i(t))\,\text{d}t = 1. \end{eqnarray*} (3) In case of a stellar system, where the phases of stellar orbits are randomized, the probability to find a star at a specific phase on its orbit does not depend on time \begin{eqnarray*} \frac{\text{d}p_i({{\boldsymbol q}}_i (t))}{\text{d}t} = 0. \end{eqnarray*} (4) Thus, taking into account equations (3) and (4) we deduce that the probability density in (3) can be written in form \begin{eqnarray*} p_i({{\boldsymbol q}}_i (t)) = \text{const} = \frac{1}{t_\text{max}-t_\text{min}}. \end{eqnarray*} (5) Here coordinates and velocities are taken along the orbit of a star i and tmin and tmax indicate the times when this star enters and leaves the region D. As we assumed that phases of stars are mixed only within a region D, for different regions the constant in equation (5) can be different. Combining individual stellar probability distribution functions of stars, we can calculate the overall probability distribution function at all points q in D \begin{eqnarray*} p({{\boldsymbol q}}) = \frac{1}{U} \sum \limits _{i}\int p_i({{\boldsymbol q}}_i(t))\, K({{\boldsymbol q}}_i(t), {{\boldsymbol q}})\, \text{d}t. \end{eqnarray*} (6) Here U is a normalizing constant, summation is over the orbits of all stars and integral is taken over the time a star spends in D. The kernel function K is needed to convert the orbital arcs into contiguous six-dimensional function by summing over all orbits. A simple example of the kernel function is a top-hat function describing whether an orbital arc qi is near to a particular point q. For a given gravitational acceleration form, we maximize the likelihood function (or likelihood multiplied by a prior when using a Bayesian analysis) to find the best parameters for the acceleration. Likelihood is constructed using observations and probability densities from equation (6) as follows \begin{eqnarray*} \mathcal {L} = \prod \limits _k p ({{\boldsymbol q}}_k) \quad \mathrm{or} \quad \log \mathcal {L} = \sum \limits _k \log p ({{\boldsymbol q}}_k). \end{eqnarray*} (7) The maximum of this function gives the values to a priori best-fitting acceleration function a(x) \begin{eqnarray*} {\boldsymbol a}({\boldsymbol x}) = {arg \, max}_{a} \mathcal {L}, \end{eqnarray*} (8) where the maximum is taken over all possible values (e.g. using a fixed parametrized form) of a. 3 TESTING OF THE METHOD ON SIMULATION To test the method, we chose two galaxy snapshots taken from the simulation (see Garbari, Read & Lake 2011) published in the Gaia Challenge wiki.2 Initial conditions of the simulation correspond to conditions giving a system close to equilibrium. Simulated galaxy consists of three components – the bulge, the disc, and the dark matter halo. One snapshot was taken at 0.049 Gyr from the start of the simulation and corresponds to a nearly axisymmetric galaxy. We denote this snapshot as a simple galaxy.3 The second snapshot was taken from 4.018 Gyr from the beginning of the simulation and is similar to a barred disc galaxy. We denote this snapshot as a barred galaxy. Surface density isophotes of these two snapshot mocks are seen in Figs 1 (a simple galaxy) and 3 (a barred galaxy). Figure 1. View largeDownload slide Calculated (red arrows) and true (green arrows) acceleration vectors in the plane of the simple snapshot of the simulated galaxy. True accelerations are normalized to constant length and simulated ones according to them. We note that some of the red and green arrows are exactly on top of each other. Simulated regions together with their names are presented as black rectangular boxes (see also Table 1). Figure 1. View largeDownload slide Calculated (red arrows) and true (green arrows) acceleration vectors in the plane of the simple snapshot of the simulated galaxy. True accelerations are normalized to constant length and simulated ones according to them. We note that some of the red and green arrows are exactly on top of each other. Simulated regions together with their names are presented as black rectangular boxes (see also Table 1). In case of the second snapshot due to the presence of the bar, accelerations in equation (2) are probably also functions of time. This snapshot is used to check sensitivity of the method to non-stationarity effects and to test the method in a more realistic situation. Regions D used for calculations may have rather different shapes (box, ellipsoid, cone etc) and sizes. To calculate p with a sufficient accuracy, these regions should contain a sufficient number of stars. In addition, as our intention is to determine acceleration of stars but since acceleration acts slowly (it takes time or needs a sufficient distance) these regions need to be large enough. Simulation results used to test the method are given in cartesian coordinates. Due to this we selected regions D as rectangular boxes with sides aligned with the coordinate axes. Relative proportions of boxes were chosen from the criteria that in all directions for fast stars (mean velocity + two standard deviations) it would take about the same time to pass through the box. Overall dimensions of regions were chosen by demanding that regions should contain sufficient number of stars. We selected 12 different regions. All regions contain |$100\, 000 \pm 20\, 000$| stars. Since we tested our results also for shot noise with smaller samples and in order to keep the tests and results consistent, all these regions were modelled with |$50\, 000$| stars. Locations of regions were chosen to include most relevant regions in the test galaxies: intermediate ‘Solar neighbourhood’ regions, outer disc regions, and regions surrounding the central bulge or bar. In vertical |$z$| direction centres of the regions were selected to lie outside of the galactic plane to test the method in various conditions. For both snapshots regions positions were the same, although due to slightly different kinematical properties their sizes were slightly different. These 12 regions superposed to the surface density isophotes for two snapshot mocks are shown in Figs 1 and 3 as rectangular boxes. Analytical form for the acceleration a(x) (or gravitational potential) within D is rather free, that means it may contain quite a lot of free parameters (e.g. Taylor series with many terms). In principle, this sets a limit to the maximum size of a region. However, these limiting conditions for gravitational potential are rather weak. In this paper, acceleration components are taken to be linear functions of the coordinates \begin{eqnarray*} a_j = a_{0,j} + b_{j}\cdot (x_j - x_{\mathrm{0},j}). \end{eqnarray*} (9) Here j indicates cartesian coordinates, aj, 0 is the acceleration at the centre of the region, x0, j and bj are linear term parameters. This form (compared to e.g. aj = a0, j + bjxj) helps to keep the model parameters aj and bj a priori maximally independent, which is advantageous during modelling, and keeps the interpretation of aj as the acceleration in the centre of the box. In case of |$z$| component it is possible to use an assumption of intrinsic symmetry of a galaxy with respect to the galactic plane and to write |$z$| component of equation (9) in a different form: a|$z$| = b|$z$| · ||$z$| − |$z$|0|. In this paper, for all components we used a general form (equation 9). In this case for every region we have six free parameters to fit. To test the accuracy of the model, true accelerations were calculated directly by using least square regression to the gravitational potential values aj = −∂Φ/∂xj at the locations of simulation particles. Regression details (e.g. whether we use all the stars from the region to calculate potential gradient or only a subsample of them, whether to add mixed terms to the aj form etc) introduce uncertainties to the true acceleration values up to few hundred km2 s−2 kpc−1. This is the maximum accuracy we may expect from our tests. For the kernel function K in equation (6) we adopted a gridding approach: in a particular region D a six-dimensional grid was constructed and for each orbit the stellar orbital probabilities were summed within a grid cell \begin{eqnarray*} K ({{\boldsymbol q}}_i, {{\boldsymbol q}}) = {1}(\lbrace {{\boldsymbol q}}_i, {{\boldsymbol q}} \rbrace \in g). \end{eqnarray*} (10) Here |$\bf {1}$| is an indicator function describing whether qi and q are in the same grid-cell g. In this case, p(q) will be a sum of top-hat functions. This approach does not widen probability density profile perpendicular to the orbit but bins, hence there will be no artificial distortion to the p. Number of grid cells for all dimensions were taken six, thus the total number of cells in a particular region is rather large (about 47 000) and in average there was one star per cell. For the fitting process, we used multinest code for the Bayesian analysis (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2013). It is a well-tested code that is aimed to converge to the solution in cases of multimodal posterior distributions. The selected prior distributions were uniform with rather wide limits |$-28\, 000\dots 28\, 000$| km2 s−2 kpc−1 for aj, 0 and |$-28\, 000\dots 28\, 000$| km2 s−2 kpc−2 for bj. We used 5000 live points during the fitting, and adopted the mean of the equal weight output posterior samples as our modelled values. 4 RESULTS 4.1 Overall recovery of the accelerations Calculated and true radial acceleration values |$a_{ R} = \sqrt{a_x^2 + a_y^2}$| and their directions for the simple galaxy mock are seen in Figs 1 and 2. Lengths of arrows that describe true values were fixed and modelled values were adjusted to those with correct ratios. It is seen that recovery of true accelerations is rather good, calculated aR values differ from true ones by only 4 per cent (median). Directions of acceleration vectors differ from true ones by only 1 deg (median). Individual fitting results are given numerically in Table 1. Figure 2. View largeDownload slide Modelled and true parameters together with error bars of accelerations [see equation (9)] for 12 regions in the simple snapshot of the simulation. The errors are marginalized standard deviations of the posterior distribution samples. It is seen that acceleration components are rather well recovered (upper panels), but correction terms for accelerations for several regions have rather large uncertainties. Figure 2. View largeDownload slide Modelled and true parameters together with error bars of accelerations [see equation (9)] for 12 regions in the simple snapshot of the simulation. The errors are marginalized standard deviations of the posterior distribution samples. It is seen that acceleration components are rather well recovered (upper panels), but correction terms for accelerations for several regions have rather large uncertainties. Table 1. General parameters of selected regions and accuracies of the calculated accelerations. The first columns give the name of the region and galactocentric cylindrical coordinates (R and |$z$|⁠) of the regions centres. The rest of the columns are divided according to the test snapshots, both having four parameters: Δ|$z$|/2 as semi-height of the modelled region, aR as planar acceleration relative error [(true value − modelled value)/true value], a|$z$| as vertical acceleration relative error, and directional difference between modelled and true acceleration vectors. Simple galaxy Barred galaxy Name R z Δ|$z$|/2 aR a|$z$| Direction Δ|$z$|/2 aR a|$z$| Direction (kpc) (kpc) (kpc) (%) (%) (deg) (kpc) (%) (%) (deg) R1 3.13 0.2 0.2 −0.13 −1.1 1.6 0.28 9.04 23.87 11.4 R2 3.47 0.34 0.23 4.3 7.84 0.4 0.25 3.07 −1.93 6.5 R3 9.07 −0.06 0.2 5.72 14.5 0.6 0.21 5.89 63.6 0.5 R4 3.1 −0.22 0.18 8.57 7.45 0.9 0.23 11.92 16.06 13.1 R5 4.62 0.1 0.18 3.98 23.24 1.1 0.25 8.64 10.68 0.5 R6 7.99 0.15 0.21 6.31 15.22 1.2 0.21 1.63 0.48 0.5 R7 5.58 −0.15 0.19 1.66 3.82 1.3 0.2 4.36 22.83 1.5 R8 6.38 −0.15 0.21 5.88 15.12 0.7 0.22 6.36 6.36 1.1 R9 3.08 −0.15 0.21 2.24 17.64 0.6 0.31 4.66 −0.16 2.7 R10 5.39 −0.3 0.23 0.78 7.08 1.3 0.29 4.88 7.06 1.2 R11 11.76 0.1 0.22 6.69 32.01 0.7 0.22 12.27 4.81 1.7 R12 11.92 0.11 0.3 0.41 12.51 0.5 0.23 −3.18 18.05 0.8 Simple galaxy Barred galaxy Name R z Δ|$z$|/2 aR a|$z$| Direction Δ|$z$|/2 aR a|$z$| Direction (kpc) (kpc) (kpc) (%) (%) (deg) (kpc) (%) (%) (deg) R1 3.13 0.2 0.2 −0.13 −1.1 1.6 0.28 9.04 23.87 11.4 R2 3.47 0.34 0.23 4.3 7.84 0.4 0.25 3.07 −1.93 6.5 R3 9.07 −0.06 0.2 5.72 14.5 0.6 0.21 5.89 63.6 0.5 R4 3.1 −0.22 0.18 8.57 7.45 0.9 0.23 11.92 16.06 13.1 R5 4.62 0.1 0.18 3.98 23.24 1.1 0.25 8.64 10.68 0.5 R6 7.99 0.15 0.21 6.31 15.22 1.2 0.21 1.63 0.48 0.5 R7 5.58 −0.15 0.19 1.66 3.82 1.3 0.2 4.36 22.83 1.5 R8 6.38 −0.15 0.21 5.88 15.12 0.7 0.22 6.36 6.36 1.1 R9 3.08 −0.15 0.21 2.24 17.64 0.6 0.31 4.66 −0.16 2.7 R10 5.39 −0.3 0.23 0.78 7.08 1.3 0.29 4.88 7.06 1.2 R11 11.76 0.1 0.22 6.69 32.01 0.7 0.22 12.27 4.81 1.7 R12 11.92 0.11 0.3 0.41 12.51 0.5 0.23 −3.18 18.05 0.8 View Large Table 1. General parameters of selected regions and accuracies of the calculated accelerations. The first columns give the name of the region and galactocentric cylindrical coordinates (R and |$z$|⁠) of the regions centres. The rest of the columns are divided according to the test snapshots, both having four parameters: Δ|$z$|/2 as semi-height of the modelled region, aR as planar acceleration relative error [(true value − modelled value)/true value], a|$z$| as vertical acceleration relative error, and directional difference between modelled and true acceleration vectors. Simple galaxy Barred galaxy Name R z Δ|$z$|/2 aR a|$z$| Direction Δ|$z$|/2 aR a|$z$| Direction (kpc) (kpc) (kpc) (%) (%) (deg) (kpc) (%) (%) (deg) R1 3.13 0.2 0.2 −0.13 −1.1 1.6 0.28 9.04 23.87 11.4 R2 3.47 0.34 0.23 4.3 7.84 0.4 0.25 3.07 −1.93 6.5 R3 9.07 −0.06 0.2 5.72 14.5 0.6 0.21 5.89 63.6 0.5 R4 3.1 −0.22 0.18 8.57 7.45 0.9 0.23 11.92 16.06 13.1 R5 4.62 0.1 0.18 3.98 23.24 1.1 0.25 8.64 10.68 0.5 R6 7.99 0.15 0.21 6.31 15.22 1.2 0.21 1.63 0.48 0.5 R7 5.58 −0.15 0.19 1.66 3.82 1.3 0.2 4.36 22.83 1.5 R8 6.38 −0.15 0.21 5.88 15.12 0.7 0.22 6.36 6.36 1.1 R9 3.08 −0.15 0.21 2.24 17.64 0.6 0.31 4.66 −0.16 2.7 R10 5.39 −0.3 0.23 0.78 7.08 1.3 0.29 4.88 7.06 1.2 R11 11.76 0.1 0.22 6.69 32.01 0.7 0.22 12.27 4.81 1.7 R12 11.92 0.11 0.3 0.41 12.51 0.5 0.23 −3.18 18.05 0.8 Simple galaxy Barred galaxy Name R z Δ|$z$|/2 aR a|$z$| Direction Δ|$z$|/2 aR a|$z$| Direction (kpc) (kpc) (kpc) (%) (%) (deg) (kpc) (%) (%) (deg) R1 3.13 0.2 0.2 −0.13 −1.1 1.6 0.28 9.04 23.87 11.4 R2 3.47 0.34 0.23 4.3 7.84 0.4 0.25 3.07 −1.93 6.5 R3 9.07 −0.06 0.2 5.72 14.5 0.6 0.21 5.89 63.6 0.5 R4 3.1 −0.22 0.18 8.57 7.45 0.9 0.23 11.92 16.06 13.1 R5 4.62 0.1 0.18 3.98 23.24 1.1 0.25 8.64 10.68 0.5 R6 7.99 0.15 0.21 6.31 15.22 1.2 0.21 1.63 0.48 0.5 R7 5.58 −0.15 0.19 1.66 3.82 1.3 0.2 4.36 22.83 1.5 R8 6.38 −0.15 0.21 5.88 15.12 0.7 0.22 6.36 6.36 1.1 R9 3.08 −0.15 0.21 2.24 17.64 0.6 0.31 4.66 −0.16 2.7 R10 5.39 −0.3 0.23 0.78 7.08 1.3 0.29 4.88 7.06 1.2 R11 11.76 0.1 0.22 6.69 32.01 0.7 0.22 12.27 4.81 1.7 R12 11.92 0.11 0.3 0.41 12.51 0.5 0.23 −3.18 18.05 0.8 View Large In vertical direction calculated accelerations a|$z$| differ from the true values by 14 per cent (median). This poorer fit is expected as in case of |$z$| direction our used form for a|$z$| [see equation (9)] is too simple to describe rather quick density changes with |$z$|⁠. In |$z$| direction near to the galactic plane it would be better to use a more complicated form for a|$z$|(x, y, |$z$|⁠) with square and even cubic terms and to use thinner regions. However, in this case the time stars spend in these regions decreases also. Increasing the number of stars should compensate this, but unfortunately, we are limited by the available simulation data. We want to point out that the purpose of this paper is to introduce the method and to make simple tests. Thus, selected form for accelerations should be handled as a first approximation only. In case of more sophisticated modelling the form for acceleration should depend on a specific application and availability and quality of the data. Overall, the coefficients (bj) in equation (9) contributed about 10 per cent to the overall acceleration values. Taking into account that within these rather small regions it is not easy to derive accelerations from velocities one may say that derived bj values are quite uncertain. In principle, increasing the number density of stars in regions should allow to decrease uncertainties of bj, but the number density of stars in the used simulation is limiting that. We can notice some correlation (correlation coefficient 0.5) between the relative error of recovered a|$z$| and aR, and between the orientation of the region and underlying rotational direction of stars (correlation coefficient 0.6). This is related with the choice of the box location and geometry and can be taken into account and improved in real applications. Real galaxies often include non-axisymmetric components and non-stationarities. To test the method for these galaxies we used the method to model a barred galaxy snapshot. Results of the model calculations are shown in Figs 3 and 4, and numerically given in Table 1. Modelling accuracies in this case for aR and a|$z$| are 5 and 9 per cent (median values), respectively. Directions of acceleration vectors were recovered with less than 2 degs (median). Figure 3. View largeDownload slide Calculated (red arrows) and true (green arrows) acceleration vectors in the plane of the barred snapshot of the simulated galaxy. True accelerations are normalized to constant length and simulated ones scaled according to them. Simulated regions together with their names are presented as black rectangular boxes (see Table 1 for numerical accuracies). Figure 3. View largeDownload slide Calculated (red arrows) and true (green arrows) acceleration vectors in the plane of the barred snapshot of the simulated galaxy. True accelerations are normalized to constant length and simulated ones scaled according to them. Simulated regions together with their names are presented as black rectangular boxes (see Table 1 for numerical accuracies). Figure 4. View largeDownload slide Modelled and true parameters together with error bars of accelerations [see equation (9)] for 12 regions of the barred snapshot. The errors are marginalized standard deviations of the posterior distribution samples. It is seen that acceleration components are rather well recovered (upper panels), but correction terms for accelerations for several regions have rather large errors. The two outliers (blue dots with their error bars) in |ay| panel are the ones lying near to the bar (R1 and R4). Figure 4. View largeDownload slide Modelled and true parameters together with error bars of accelerations [see equation (9)] for 12 regions of the barred snapshot. The errors are marginalized standard deviations of the posterior distribution samples. It is seen that acceleration components are rather well recovered (upper panels), but correction terms for accelerations for several regions have rather large errors. The two outliers (blue dots with their error bars) in |ay| panel are the ones lying near to the bar (R1 and R4). It is seen from Table 1 that three least accurate solutions are for regions R1, R4, and R11. In case of R1 and R4, the model did not reproduce directions of acceleration vectors. This is not surprising, since both regions are near to the bar (see Fig. 3). In case of the region R11 derived acceleration value aR is quite inaccurate. This is seen also for the simple galaxy snapshot. From simple snapshot results we found a mild correlation (0.6) between the position angles of the regions and recovered accuracies – for regions with best alignment with x- or y-axis the accuracy of modelled accelerations are the best. This can be interpret in a way that rotating stars move through the corners of these regions very fast and do not sufficiently influence overall acceleration value. Most of the information for modelling is from bins that are near to the centre of the box and are thus derived acceleration values are geometrically biased. 4.2 Robustness of the modelling First, we checked whether multinest fitting converges to the same posterior distribution when running the fitting process with different random seeds. We reran complete fitting process for all 12 regions in case of the barred galaxy snapshot and calculated the differences of the mean values of derived aR distributions, then divided the result by the averaged standard deviation of the two runs. We found that on average (over all regions) this quantity was 0.5, the maximum difference was 1.0. We can conclude that there are some convergence issues, but as we will see, these are smaller than shot noise effects. We conclude that this difference is small enough that it does not change overall our results and conclusions. Subsequent tests were done for the region R10 in the barred galaxy sample. We preferred the barred snapshot over simple one to have more realistic estimate on how well would the model behave in real observations. To test for robustness against shot noise, we created three subsamples for R10, all of them containing |$38\, 297$| stars. The results of modelling for each of these subsamples are seen in Fig. 5. The contours show one and two sigma posterior distribution regions. We conclude that the shot noise has larger impact on resulting uncertainties than statistical uncertainties. Figure 5. View largeDownload slide Precision of calculated accelerations in galactic plane for different sampling in case of the region R10. Narrow grey lines are curves of constant acceleration in galactic plane with labels indicating the fraction from the true value. The contours show one and two sigma confidence intervals. All the samples are independent and contain equally |$38\, 297$| stars. Figure 5. View largeDownload slide Precision of calculated accelerations in galactic plane for different sampling in case of the region R10. Narrow grey lines are curves of constant acceleration in galactic plane with labels indicating the fraction from the true value. The contours show one and two sigma confidence intervals. All the samples are independent and contain equally |$38\, 297$| stars. In calculations we used regions containing about 105 test stars (although half of them were actually used) and divided modelled regions into a grid with 66 cells per region. Next we will test how robust these choices are. To test the requirements for the number of test stars, we calculated accelerations by using different numbers of stars. For this we created from the initial number of 105 stars a random subsets of 3000, |$10\, 000$|⁠, and |$30\, 000$| stars. Box sizes remained the same. Fitting results are given in Fig. 6. We can see that for this region about 30 000 stars are needed to recover acceleration ax and ay components within about 5 per cent. Of course the minimum number of stars depends on how smooth is the gravitational potential, are there perturbations in phase density distribution, what is the size and shape of the region etc. This number of stars needed for accurate modelling should be handled as an order-of-magnitude estimate. Figure 6. View largeDownload slide Precision of calculated accelerations in galactic plane as a function of the number of test stars used in modelling in case of the region R10. Narrow grey lines are the curves of constant acceleration in galactic plane with labels indicating the fraction from the true value. The contours show one and two sigma confidence intervals. Figure 6. View largeDownload slide Precision of calculated accelerations in galactic plane as a function of the number of test stars used in modelling in case of the region R10. Narrow grey lines are the curves of constant acceleration in galactic plane with labels indicating the fraction from the true value. The contours show one and two sigma confidence intervals. Next we tested sensitivity of the modelling process to grid parameters, more precisely to the number of cells within a given region. Number of test stars remained unchanged (50 000 stars). Four different combinations of numbers for grid coordinates and velocities were used. In all combinations there is about one test star per grid cell. Results of these tests are given in Fig. 7. It is seen that there are no systematic trends how the grid is formed as long as there are enough stars to fill the grid. Figure 7. View largeDownload slide Precision of calculated accelerations in galactic plane for different grids in case of the region R10. Narrow grey lines are curves of constant acceleration in galactic plane with labels indicating the fraction from the true value. The contours show two sigma confidence intervals for each grid. Grid sizes are indicated separately for coordinate space (Nx) and for velocity space (Nv). Figure 7. View largeDownload slide Precision of calculated accelerations in galactic plane for different grids in case of the region R10. Narrow grey lines are curves of constant acceleration in galactic plane with labels indicating the fraction from the true value. The contours show two sigma confidence intervals for each grid. Grid sizes are indicated separately for coordinate space (Nx) and for velocity space (Nv). 5 SUMMARY AND DISCUSSION We developed a method to calculate acceleration vectors (gravitational potential gradients) within regions containing 10s of thousands of stars with known coordinates and velocities. The method assumes that phases of stars on orbits within a given region are mixed, i.e. the exact position of a star on its orbit is a random variable with a probability density uniform in time. By knowing all phase space coordinates for stars in a given region at a particular moment, orbits of stars within the region can be calculated for some parametrized form of acceleration. Combining these orbital arcs together we calculate the probability function to find stars at some particular phase space point. Then the likelihood can be calculated using modelled stellar probability distributions with the observed phase space coordinates of the stars. The likelihood (or posterior distribution) has a maximum if the calculated probability distribution and the observed distribution of stars are consistent. This allows to set constraints to the accelerations used in stellar orbit calculations. The method assumes that phases of stars are mixed within selected regions where orbits are calculated. In all regions orbits and accelerations are calculated independently to other regions. Thus, there is no need to assume overall stationarity of a galaxy. We tested the method with two snapshots from a mock galaxy simulation. The first one (denoted as a simple galaxy) corresponds to an axisymmetric galaxy. 12 regions were selected each containing 100 000 stars. Method allows to reproduce true radial acceleration values within 4 per cent and directions within 1 deg (median values). The second snapshot corresponds to the barred galaxy having thus perturbed kinematics. In this case, the method allows to reproduce true radial acceleration values within 5 per cent and directions within 2 deg (median values). A similar representation of potential gradient or acceleration, circular velocity, can be calculated within 3 per cent and are thus also rather accurate. One of the assumptions of the method in its present form was that acceleration in equation (2) is a function of coordinates only, i.e. we assume that there are no significant perturbations near and within the studied regions D. As the tests for barred galaxy regions R1 and R4 showed (Fig. 3), indeed, the model did not reproduce true accelerations well. Our preliminary tests gave that fitting accuracy for these regions can be improved by slightly changing equation (5) i.e. allowing mild non-uniformity. Proposed method can be used for the MW using Gaia complete solution data with known spatial coordinates and velocities. By selecting modelled regions D to be sufficiently small the gravitational potential derivatives can be calculated at the whole grid of coordinates covered by Gaia data. A minimum size of a region is constrained by the timing argument (acceleration needs time) and by the statistical argument (large number of stars are needed for an accurate likelihood evaluation). Using typical velocities (>100 |$\mathrm{km\, s^{-1}}$|⁠) and typical accelerations in the MW galaxy and demanding that velocity changes due to acceleration should exceed, for example five times the velocity measurement error one obtains that in the halo or outer disc region the minimal box size is about 0.5 kpc. This box contains also sufficient number of stars for statistics. Thus the grid in the MW where it is possible to calculate acceleration components is quite dense. For a typical star it takes around 2–5 Myr to move through a box with this minimum size. This gives the minimum time-scale how long typical stellar orbits should be integrated. Another constraint is the accuracy of the measured parallaxes: the uncertainty of parallaxes should be smaller than the grid cell size. Although this hinders the applications, the errors of parallaxes are constantly reduced by longer timebase of Gaia and increasing number of observations. We presented here only one possible usage of the method. With modifications the model can be used to calculate accelerations close to the Galactic disc plane, to study disturbances near to the bar and in other interesting cases. ACKNOWLEDGEMENTS We thank the referee for useful comments and suggestions that helped to improve the paper substantially. This work was supported by institutional research funding IUT26-2 and IUT40-2 of the Estonian Ministry of Education and Research. We acknowledge the support by the Centre of Excellence ‘Dark side of the Universe’ (TK133) and grant MOBTP86 financed by the European Union through the European Regional Development Fund. We thank the Gaia Wiki team, who made their simulation data publicly available. Footnotes 1 However, surprisingly large number of galaxies have been modelled within Califa project (Zhu et al. 2018). 2 http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php?id=tests:discs:mockreadfullnbm 3 In the Gaia wiki, this snapshot is referred as unevolved snapshot. REFERENCES Barone T. M. et al. , 2018 , ApJ , 856 , 64 https://doi.org/10.3847/1538-4357/aaaf6e Crossref Search ADS Binney J. , 2018 , MNRAS , 474 , 2706 https://doi.org/10.1093/mnras/stx2835 Binney J. , McMillan P. J. , 2015 , TM: Torus Mapper , Astrophysics Source Code Library , record (ascl:1512.014) Binney J. , McMillan P. J. , 2016 , MNRAS , 456 , 1982 https://doi.org/10.1093/mnras/stv2734 Crossref Search ADS Binney J. J. , Davies R. L. , Illingworth G. D. , 1990 , ApJ , 361 , 78 https://doi.org/10.1086/169169 Crossref Search ADS Bovy J. , 2014 , ApJ , 795 , 95 https://doi.org/10.1088/0004-637X/795/1/95 Crossref Search ADS Cappellari M. , 2008 , MNRAS , 390 , 71 https://doi.org/10.1111/j.1365-2966.2008.13754.x Crossref Search ADS Cappellari M. et al. , 2006 , MNRAS , 366 , 1126 https://doi.org/10.1111/j.1365-2966.2005.09981.x Crossref Search ADS Cappellari M. et al. , 2007 , MNRAS , 379 , 418 https://doi.org/10.1111/j.1365-2966.2007.11963.x Crossref Search ADS Copin Y. , Zhao H. S. , de Zeeuw P. T. , 2000 , MNRAS , 318 , 781 https://doi.org/10.1046/j.1365-8711.2000.03827.x Crossref Search ADS de Lorenzi F. , Debattista V. P. , Gerhard O. , Sambhus N. , 2007 , MNRAS , 376 , 71 https://doi.org/10.1111/j.1365-2966.2007.11434.x Crossref Search ADS Emsellem E. , Monnet G. , Bacon R. , 1994 , A&A , 285 , 723 Feroz F. , Hobson M. P. , 2008 , MNRAS , 384 , 449 https://doi.org/10.1111/j.1365-2966.2007.12353.x Crossref Search ADS Feroz F. , Hobson M. P. , Bridges M. , 2009 , MNRAS , 398 , 1601 https://doi.org/10.1111/j.1365-2966.2009.14548.x Crossref Search ADS Feroz F. , Hobson M. P. , Cameron E. , Pettitt A. N. , 2013 , preprint (arXiv:1306.2144) Gaia Collaboration , 2018 , A&A , 616 , A1 https://doi.org/10.1051/0004-6361/201833051 Crossref Search ADS Garbari S. , Read J. I. , Lake G. , 2011 , MNRAS , 416 , 2318 https://doi.org/10.1111/j.1365-2966.2011.19206.x Crossref Search ADS Kalinova V. , van de Ven G. , Lyubenova M. , Falcón-Barroso J. , Colombo D. , Rosolowsky E. , 2017 , MNRAS , 464 , 1903 https://doi.org/10.1093/mnras/stw2448 Crossref Search ADS Kipper R. , Tempel E. , Tenjes P. , 2018 , MNRAS , 473 , 2188 https://doi.org/10.1093/mnras/stx2441 Crossref Search ADS Kowalczyk K. , Łokas E. L. , Valluri M. , 2017 , MNRAS , 470 , 3959 https://doi.org/10.1093/mnras/stx1520 Crossref Search ADS Li H. et al. , 2018 , MNRAS , 476 , 1765 https://doi.org/10.1093/mnras/sty334 Long R. J. , Mao S. , 2010 , MNRAS , 405 , 301 https://doi.org/10.1111/j.1365-2966.2010.16438.x McGaugh S. S. , Lelli F. , Schombert J. M. , 2016 , Phys. Rev. Lett. , 117 , 201101 https://doi.org/10.1103/PhysRevLett.117.201101 Crossref Search ADS PubMed McGill C. , Binney J. , 1990 , MNRAS , 244 , 634 Schwarzschild M. , 1979 , ApJ , 232 , 236 https://doi.org/10.1086/157282 Crossref Search ADS Syer D. , Tremaine S. , 1996 , MNRAS , 282 , 223 https://doi.org/10.1093/mnras/282.1.223 Crossref Search ADS Thomas J. , Saglia R. P. , Bender R. , Thomas D. , Gebhardt K. , Magorrian J. , Richstone D. , 2004 , MNRAS , 353 , 391 https://doi.org/10.1111/j.1365-2966.2004.08072.x Crossref Search ADS Valluri M. , Merritt D. , Emsellem E. , 2004 , ApJ , 602 , 66 https://doi.org/10.1086/380896 Crossref Search ADS van de Ven G. , de Zeeuw P. T. , van den Bosch R. C. E. , 2008 , MNRAS , 385 , 614 https://doi.org/10.1111/j.1365-2966.2008.12873.x Crossref Search ADS van der Marel R. P. , Binney J. , Davies R. L. , 1990 , MNRAS , 245 , 582 Zhu L. et al. , 2014 , ApJ , 792 , 59 https://doi.org/10.1088/0004-637X/792/1/59 Crossref Search ADS Zhu L. et al. , 2018 , Nature Astron. , 2 , 233 https://doi.org/10.1038/s41550-017-0348-1 Crossref Search ADS © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - A method to calculate gravitational accelerations within discrete localized regions in the Milky Way
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty2598
DA - 2019-01-11
UR - https://www.deepdyve.com/lp/oxford-university-press/a-method-to-calculate-gravitational-accelerations-within-discrete-zjDPP6sfhg
SP - 1724
VL - 482
IS - 2
DP - DeepDyve
ER -