Quantifying torque from the Milky Way bar using Gaia DR2 Free

ABSTRACT

1 INTRODUCTION

Gaia satellite data releases allow us to construct quite detailed models for the Milky Way (MW) stellar density distribution and its kinematics. The latest Data Release 2 (Lindegren et al. 2018) gives us an excellent opportunity to explore the solar neighbourhood (SN) and somewhat more distant regions. In the present paper, we calculate the gravitational acceleration of the MW using the Gaia DR2 data, in an ellipsoidal region within a distance of 0.5 kpc from the Sun in the Galactic plane.

Modelling the MW is very different from modelling other disc galaxies since we make observations from within the MW. Although our location within the MW can make modelling easier, (e.g. individual stars are resolved) it can also add complications to it, e.g. dust attenuation and selection function can have a strong influence on modelling. For example, it was only at the beginning of the 1980s that the first direct hints that the MW may be a barred spiral galaxy came to light (Matsumoto et al. 1982). This was possible because of near-IR observations. Due to dust attenuation and our position inside the MW, it was difficult to draw such a conclusion on the morphology of MW before that.

On the other hand, we are at a great advantage because of the wealth of observational data available for the MW, unmatched and unavailable for other galaxies. For instance, axisymmetric models developed by Piffl et al. (2014), McMillan (2017), Binney & Wong (2017) use H i and CO velocities, maser data, Sgr A* proper motions, the globular cluster system, the velocity distribution in the SN, SDSS star counts in different colours, RAVE data, detailed MW satellite data and N-body simulation data. Additional constraints on the mass distribution were derived from cold stellar streams (Bovy et al. 2016) and Gaia DR2 proper motions of globular clusters (Watkins et al. 2019). However, the assumption of axisymmetry in mass distribution models is only a first approximation.

The existence of the central bar of the MW was first confirmed by Weiland et al. (1994), by analysing asymmetries in the near-IR surface brightness distribution of the central bulge from the COBE/DIRBE data. This was further confirmed by correcting the data for extinction (Dwek et al. 1995; Binney, Gerhard & Spergel 1997).

There are currently two contrasting scenarios – a fast rotating bar and a slow rotating bar. In the first case (e.g. Binney et al. 1997; Bissantz, Englmaier & Gerhard 2003; Monari et al. 2017) the bar is rotating with pattern speed |$\Omega _\mathrm{p} = 50\!-\!70\,\mathrm{km\, s^{-1}\, kpc^{-1}}$|⁠; in the second case (Wegg & Gerhard 2013; Wegg, Gerhard & Portail 2015; Portail et al. 2015; Dias et al. 2019) the calculated pattern speed is |$25\!-\!30\,\mathrm{km\, s^{-1}\, kpc^{-1}}$|⁠. Intermediate pattern speed values were derived by Li et al. (2016), Portail et al. (2017), Pérez-Villegas et al. (2017), Sanders, Smith & Evans (2019), Bovy et al. (2019) as |$\Omega _\mathrm{p} = 35\!-\!40\,\mathrm{km\, s^{-1}\, kpc^{-1}}$|⁠. These calculated pattern speeds vary by about two times and as a result their corotation radii and outer Lindblad radii vary quite significantly. Both these scenarios agree that the angle between the major axis of the bar and the line connecting the Sun and the Galactic centre (GC) is about |$20\!-\!30\, \mathrm{deg}$|⁠.

According to the axisymmetric models, the stars in orbits are phase-mixed. According to the non-axisymmetric models, stellar orbits are somewhat perturbed and phases of stars in orbits may not be completely mixed (Dehnen 2000; Fux 2001; Monari et al. 2017; Binney 2018; Trick, Coronado & Rix 2019) and thus orbital structure is more complicated (i.e. there are resonances). For example, using the Gaia DR2 data, Ramos, Antoja & Figueras (2018) found that, in the case of the MW, some orbit phases are mixed. Similar arcs and ridges were also found by Antoja et al. (2018), Kawata et al. (2018), Trick et al. (2019). Gravitational potential disturbances due to the bar may have caused deviations of stars from their initial orbits in the case of several cold stellar streams (Hattori, Erkal & Sanders 2016; Pearson, Price-Whelan & Johnston 2017; Banik & Bovy 2019) that originate from small stellar systems. The torque from the bar is not the only reason (see e.g. Kipper et al. 2019b). These disturbances can create observed gaps in stream surface density distributions.

Unfortunately, the structural parameters of the bar and its contribution to the gravitational acceleration are still rather poorly constrained. Thus, it is important to know the gravitational acceleration distribution in the Galactic plane and also to study how this allows one to constrain the bar properties. The Gaia satellite data provide an excellent opportunity to do this.

In the present paper we calculate all three acceleration components in the SN. We use the orbital arc method, developed in Kipper, Tempel & Tenjes (2019a). The method and its specific implementation details are described in Section 2. The method is used on the Gaia DR2 data. We use two different versions of the data, from the StarHorse project (Anders et al. 2019) and from the Schönrich catalogue (Schönrich, McMillan & Eyer 2019). The selection of the data used is described in Section 3. We demonstrate in Section 4 that the derived acceleration components cannot be explained within an axisymmetric model. The final section is dedicated to the summary and discussion.

We denote (x, y, z) as Galactocentric rectangular coordinates and (R, θ, z) as corresponding cylindrical coordinates, where θ = 0 corresponds to the opposite direction from the Sun. Transformations of sky coordinates, proper motions and radial velocities to Galactocentric coordinates and velocities were carried out using the astropy package (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018). For the solar velocity, we used the values |$(\mathrm{U}_\odot , \mathrm{V}_\odot , \mathrm{W}_\odot) = (11.1,12.24,7.25)\, \mathrm{km}\, \mathrm{s}^{-1}$|⁠, and V_{g, ⊙} = V_{c, ⊙} + V_⊙, with the circular velocity |$V_{c,\odot } = 240\, \mathrm{km}\, \mathrm{s}^{-1}$| (López-Corredoira & Sylos Labini 2019).

2 METHOD AND IMPLEMENTATION

2.1 Orbital arc method

In this section, we provide a brief overview of the orbital arc method, which we have implemented in this paper to compute the gravitational acceleration, mass and torque estimates of the Galactic bar. For a detailed and thorough description of the formulation and tests of the model, please see Kipper et al. (2019a). We will refer to this particular method as the orbital arc method, since its most crucial step is the reconstruction of stellar orbits to accurately obtain the acceleration in the Milky Way, using the phase space information of stars. This has already been successfully applied to a simulation in Kipper et al. (2019a) and for the observational data in a simplified form (Kipper, Tempel & Tenjes 2018). Here we apply it for the Gaia DR2 data. The orbital arc method has five important steps.

• Step 2 – Orbital arc reconstruction: Using the initial conditions, which is the phase space information from the data, and the acceleration field from the previous step, we solve the equations of motion to obtain stellar orbits for each star. A schematic of the reconstructed orbit arcs is represented as coloured arcs in Fig. 1.

• Step 3 – Randomizing star position: The core of the proposed method lies in the oPDF, according to which the time of observing a star is random. This means that we can reposition a star along its orbit by picking a random time from a uniform distribution of time. By picking infinitely many times from this distribution, we reach a continuous distribution of the star along its orbit (a similar description to the procedure can be found in Han et al. 2016). Following this, we reposition every star in its orbit and get a new distribution of stars in the region. This is relevant in the last step where we will compare the old and new distributions.

• Step 4 – Phase space density: In order to compute the probability of finding a star in its orbit, we need to first specify a small segment of the star’s orbit. For this, we construct Voronoi tessellations by considering small subsets of data, such that in each Voronoi cell there are similar numbers of stars; this reduces the Poisson error. The Voronoi cells are shown in the schematic diagram in Fig. 1. To compute the probability of finding a star in its orbit, we use the Voronoi cells as the orbital segments. For example, for star 1 in Fig. 1, the time spent by the star in each Voronoi cell along its orbit is given as Δt_i. So, the probability of finding that star in the ith Voronoi cell is the time spent in that Voronoi cell, Δt_i, divided by the time spent in the entire region, which is, Δt_i/(t_exit − t_enter) as seen in Fig. 1. Eventually, a combined probability is calculated for each Voronoi cell, which is the sum of probabilities of all stars in each cell.

• Step 5 – Comparing phase space density distributions: In this final step we compare the phase space distribution of the original data and the phase space distribution of the newly positioned stars. The phase space distribution comparison is done statistically by computing the likelihood. If the likelihood is not maximum then the entire process is repeated with new acceleration field parameters.

Figure 1.

An illustration of the region where the orbital arc method is applied. The central point represents the Sun and the circular region up to a distance of 0.5 kpc from the Sun is the region used in this paper. The black triangle points towards the Galactic centre. The coloured arcs for star 1 and star 2 represent the reconstructed orbits for these two stars. Orbital arcs are reconstructed for all stars in this region. The grey cells are the Voronoi cells, each of which contains a similar number of stars. The time interval, Δt_i, is the time a star spends in the ith Voronoi cell. The times at which the star enters and exits the region are t_enter and t_exit, respectively.

Open in new tabDownload slide

Eventually, the orbital arc method will give the acceleration field corresponding to the maximum likelihood, which is the field that describes the true underlying acceleration of the MW. The distribution of likelihoods gives the statistical uncertainty.

Since the level of accuracy relies on the available data, we need the phase space coordinates of a sufficiently large number of stars. Hence, Gaia DR2 is aptly suited for the study.

2.2 Implementation: the smoothing kernel

In order to compare the phase space distributions of the original data and the repositioned stars, we have used Voronoi cells to get a smooth phase space density. This is achieved by computing the time stars spend in each Voronoi cell, as described in step 4 in Section 2.1. In Fig. 1, Δt represents the time a star spends in a cell. The ratio of the time spent by a star in a Voronoi cell, Δt_i, to the time that it spends in the entire region (see t_enter and t_exit in Fig. 1) gives the phase space density of this model.

The shapes and sizes of the regions in which we intend to calculate the accelerations are mainly motivated by the quality of 6D data. The shapes of these regions and the Voronoi cells used to smooth phase space¹ should be complementary to each other. For example, if the available data are of a spherical region, then a rectangular grid or cell is not the most optimal. Therefore, the best possible grid should coarsely follow the distribution of data. One of the best ways to achieve this is by Voronoi tessellations, and we have hence used this method for the paper. However, in principle, any similar grid can be used. We used a random small subset of the data of about ∼100 stars to obtain the Voronoi cells.

Each grid-cell is described by two numbers: the closest tessellation centre in ordinary space and the closest tessellation centre in velocity space. These two indices are required to avoid combining velocity and distance data into a single quantity, because this kind of combination produces an additional free parameter that we wish to avoid. For example, by using 100 data points to tessellate into a grid, we will have |$100^2 = 10\, 000$| independent cells, which is usually sufficient to describe the phase space distribution of about |$\approx 420\, 000$| stars (i.e. 42 stars per cell). For the current study, we have selected 100 cells for each of position and velocity space, unless noted otherwise.

2.3 Implementation: flux limitedness

Flux-limited observational data are a natural constraint in large surveys. There are two common approaches to overcome this: a) to construct a volume-limited sample and discard some data, or b) to use all the data and add a weight to each point.

Most dynamical modelling methods are constructed based on the assumption that we are able to observe everything, i.e. the volume-limited approach. Some specifics of the present modelling allow us to use the advantage of increased amounts of data of the flux-limited selection, while essentially using the method constructed for the volume-limited approach. This approach is described further in this section.

2.4 Requirements for data

An integral part of the method is orbit calculation. This has two ingredients: the proposed acceleration function and initial conditions for the orbits. As an analytical expression the first one is infinitely precise for each likelihood evaluation. The second one is as precise as the data allow. Imprecisions in the data are amplified by the orbit integration, i.e. Δx ∼ Δx₀ + tΔv, where Δ denotes uncertainty for positions and velocities respectively. This shows that the uncertainties accumulate with time; hence the position of a star is unknown in some cone. Due to uncertainties (especially heteroscedastic ones) in the Gaia data combined with smoothing phase space, we may reconstruct imprecise orbits, which will introduce a bias in the acceleration. The simplest way to avoid these problems is to use maximally precise data.

The second requirement is to have a sufficient amount of data. This is needed to describe the phase space density sufficiently precisely. Assuming that the data are very precise, the only source of uncertainty is the Poisson noise from the sampling.

3 OBSERVATIONAL DATA

3.1 Construction of the data sample

Six-dimensional high-quality phase space coordinates in the SN are now available from the Gaia satellite Data Release 2 for a significant number of stars. At present there are three catalogues available based on the Gaia measurements and including estimated star distances: the Gaia Collaboration catalogue (Lindegren et al. 2018), the StarHorse project catalogue, SH (Anders et al. 2019) and the Schönrich catalogue, Sc (Schönrich et al. 2019). There is a known issue concerning the zero point of parallaxes from the Gaia Collaboration, which is overcome in the latter two catalogues. Therefore we selected these two catalogues as our main sources of input data and calculated our results for both of them separately. To calculate gravitational acceleration, we need to know mass density gradients. Although the main source of density gradients results from the smooth density distribution of the MW, selection effects can produce artificial gradients. The two dominant ingredients for this kind of selection effects are Malmquist bias (covered in the previous section) and dust attenuation. To suppress the effects from dust attenuation, we use 2MASS (Skrutskie et al. 2006) catalogue J-band magnitudes where extinction is negligible.

The cross-match between the Anders et al. (2019) SH and 2MASS catalogues gave 6964 515 entries; between the Schönrich et al. (2019) Sc and 2MASS catalogues there were 6519 209 matches. We constrained the input magnitudes in such a way that the Gaia G-band completeness (being affected by dust attenuation) has substantially less effect than our selection based on the J passband (being nearly attenuation free), i.e. P(G > G_lim|J < J_lim) ≪ 1. The apparent-magnitude data within 0.5 kpc from the Sun for our selected sample are shown in Fig. 2. A strong correlation between the G- and J-band magnitudes catches the eye. The J-band limit J_lim was fixed to a value where the distribution of brighter stars in the G band ends before reaching the Gaia spectroscopic completeness limit G_lim. This is shown as a green line in the left-hand panel and the corresponding probability density distribution p(G|J < J_lim)dG in the right-hand panel. The fraction of G magnitudes crossing G_lim is 8 × 10⁻⁴ for the adopted J_lim = 10.25; hence we conclude that our sample is almost independent of the Gaia completeness limit and dust attenuation.

Figure 2.

Distribution of apparent magnitudes of all stars within 0.5 kpc from the Sun. The G-band magnitudes are from the Gaia data and the J-band magnitudes are from the 2MASS survey. The red horizontal line and the green vertical line depict the spectroscopic completeness limit and the limiting magnitude J_lim, respectively, of our main sample. The right-hand panel shows the distribution of all stars from Gaia. The black line shows all the stars and the green line shows only those with magnitudes brighter than J_lim. The distribution of our sample of stars drops before reaching the Gaia completeness limit. Only a fraction of 0.0008 stars have a higher G-band magnitude; therefore we choose our sample based on 2MASS photometry.

Open in new tabDownload slide

The smooth acceleration distribution of the MW is taken as an input in modelling process and it does not include local potential wells of stellar clusters. Hence, we cannot describe the motion of stars within clusters and must exclude these cluster stars from our sample. We excluded all stars that appeared to be cluster members in catalogues by Cantat-Gaudin et al. (2018) or the Gaia Collaboration et al. (2018). In total, 993 stars or about 0.2 per cent of stars from the final selection were excluded.

3.2 Selection of the region

In the paper where we presented the method and tested it on simulation data (Kipper et al. 2019a), we aimed to use rather small regions in order to have a simple analytical form for acceleration vector components. In the current paper, we selected a larger region to suppress Poisson noise and to increase the region size in the radial direction to have a stronger basis to also estimate the first derivative of the radial acceleration. This changes our approach somewhat: instead of using a simple form for accelerations, we now try to model the underlying acceleration field with a well-motivated analytical form.

Thus, due to available data, instead of using several small regions as we did in Kipper et al. (2019a), we selected one larger region, as shown in the schematic in Fig. 1. Our main aim was to recover the acceleration field in the plane of the Galaxy; hence we constructed a region where accelerations in the MW plane have a longer time to act on stars. In the vertical direction, density gradients are much steeper and one may expect that the corresponding accelerations may have also more complex forms. To avoid using more complex accelerations in vertical directions, we selected a thin region.

4 RESULTS

We calculated gravitational acceleration in the region around the SN as described in Section 3.2 using the method and its implementation explained in Section 2. In order to decipher various aspects of the acceleration field (e.g. deviations from axisymmetry), we used different functional forms to describe the underlying acceleration.

4.1 Calculated acceleration components

In our first attempt to model the acceleration in the region we did not specify a design-based form of an overall gravitational potential of the MW. Instead we assumed that any acceleration form can be approximated with their Taylor expansions equations (7)–(9) and we fit the coefficients of this acceleration (A_x, A_{x, x}, A_{x, y}, A_{x, z}, A_y, A_{y, x}, A_{y, y}, A_{y, z}, A_z, A_{z, x}, A_{z, y}, A_{z, z}). This way of modelling is powerful because it allows us to not only model the acceleration in a tiny region, but in principle the entire MW if we can get the overall gravitational potential. We fit a total of 12 free parameters, A_i and A_{i, j} for the flux-limited samples of stars within the selected region (see Section 3.2) for both the Sc and SH catalogues of Gaia DR2 (for more details see Section 3.1).

We used 100 random points to describe the grid; hence, there are ≈42 stars per grid bin. The fitting was done with the multinest code (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2013) using 500 live points. To include the randomness caused by the gridding, we ran the code eight times and averaged the posterior distribution of different runs. All the modelling was done in this way, unless noted otherwise. The priors of the Bayesian modelling were chosen to be of uniform distribution with the limiting values provided in Table A1. In the table we give the posterior distribution of each parameter ζ with five quantiles positioned at P(ζ) = {0.023, 0.159, 0.500, 0.841, 0.977}.

Previous studies have shown that the Sun is not located precisely at the centre of the Galactic plane, but is about |$25\, {\rm pc}$| away from it (Bland-Hawthorn & Gerhard 2016). Thus, there must be an acceleration component in the vertical direction, as confirmed in e.g. Kipper et al. (2018) based on dynamics. In Table A1 our estimation of the vertical acceleration A_z and its gradient in the z-direction A_{z, z} are given. Using these two values and by making a linear approximation at distances close to the plane, we deduce that we are located at |$\mathrm{z}_\odot \approx A_z/A_{z,z} = \lbrace -111^{+33}_{-76} ({\rm SH}),\, -117^{+58}_{-66} ({\rm Sc})\rbrace$| pc from the vertical coordinate value defined as having zero vertical acceleration in contrast to the symmetry-defined centre.

Figure 3.

The figure shows the average matter density in the solar neighbourhood and is compared with the results from Bienaymé et al. (2014), McKee, Parravano & Hollenbach (2015), Kipper et al. (2018). These results do not match very well because they use different datasets and different assumptions of the underlying acceleration. The high uncertainty in the calculated results is due to the optimization of the selected acceleration form to determine accelerations in the Galactic plane. Note that this is not the vertical component, which is usually used to determine the overall matter density.

Open in new tabDownload slide

4.2 Deviations from a simple axisymmetric MW model

In case of a stationary axisymmetric MW, the acceleration component along the direction of Galactic rotation a_y(Δx, Δy = 0, Δz) = 0 and equipotential curves are concentric circles. The median values of acceleration computed within the selected region using the coordinates (x, y, z) = (− 8.3, 0, 0) kpc as the centre are 306 and |$284\, {\rm km^2\, s^{-2}\, kpc^{-1}}$| for the SH and Sc catalogues respectively. Results of these calculations along with 1σ and 2σ limits are shown in Fig. 4 and the used priors are given in Table A2. They are designated as ‘Sc, flux’ and ‘SH, flux’. None of the |$27\, 748$| posterior samples from multinest show negative A_y values. Thus, the results are not consistent with axisymmetry.

Figure 4.

The central panel shows the correlation between the centre of acceleration, R_⊙, and the component of the acceleration vector, A_y, as described in Section 4.2. The top panel shows the distribution of R_⊙ for an axisymmetric fit. The brown point at 8.3 kpc shows the distance of the Sun to the centre of our Galaxy. This indicates that the curvature of the isopotential lines is very likely less than 8.3 kpc. The right-hand panel shows the distribution of the A_y component of the acceleration vector at the region centre. The green lines are for the overall posterior distribution and the red lines are for the subset where |$\mathrm{R}_\odot \approx 8.3\, {\rm kpc}$|⁠. This figure highlights the necessity to include the non-axisymmetric component to fit the acceleration, since the default for the MW at (R_⊙ = 8.3, A_y = 0) does not account for what is observed.

Open in new tabDownload slide

Based on the assumption that equipotential curves are concentric circles, we derived the radius of this circle. The radii are |$3.4\, {\rm kpc}$| for the SH catalogue and |$3.2\, {\rm kpc}$| for the Sc catalogue. Most of the posterior distribution had values lower than 8.3 kpc, i.e. |$P(\mathrm{R}_\odot \,\gt \,8.3\, {\rm kpc}) = \lbrace 0.048 ({\rm Sc}),\, 0.11 ({\rm SH}) \rbrace$|⁠. Hence, the ‘acceleration centre’ is most likely closer to us than the GC and equipotential curves have higher curvature than one would expect for the distance to the Galactic centre. Thus we conclude that axisymmetric potential distribution is not valid at SN, and interpret it as an argument to support the presence of a rather massive central bar.

As already explained in Section 5.1, we calculated the acceleration components assuming axisymmetry, by selecting the components a_R, a_z to be in the form of equations (4)–(6). During the fitting the solar distance R_⊙ was also taken as a free parameter. Taking the posterior in these fits close to |$\mathrm{R}_\odot =8.3\, \mathrm{kpc}$|⁠, we found the radial acceleration to be |$-6190^{+70}_{-160}\, {\rm km^2\, s^{-2}\, kpc^{-1}}$|⁠, which corresponds to the circular velocity 227 km s⁻¹. The combined estimate of the observed circular velocity at 8.3 kpc is somewhat larger, being |$238\pm 15\,{\rm km\, s^{-1}}$| (Bland-Hawthorn & Gerhard 2016), but it is consistent with the calculated result within error.

4.3 Deriving the properties of the bar

When using accelerations in the form of equations (16)–(18), substantial correlations exist in the modelled posterior samples (e.g. between A_r and A_{r, bar}). The largest correlation coefficient was found to be 1.0 between A_r and A_{r, bar} (see Fig. 5). We also found a strong correlation (correlation coefficient 0.85) between bar acceleration/total acceleration and its focal length as shown in Fig. 6. The rest of the correlations were much weaker, and the only other noticeable correlation was between A_R and A_{R, bar} and the radial acceleration derivative A_{R, R} (correlation coefficient 0.26).

Figure 5.

The relation between acceleration from the axisymmetric component and from the bar component. The contours show 1σ and 2σ confidence intervals for the Anders et al. (2019) (SH) and Schönrich et al. (2019) (Sc) datasets. The strong correlation between them shows degeneracy of accelerations in the functional form of equations (16)–(18).

Open in new tabDownload slide

Figure 6.

The top panel depicts the degeneracy between the length of the bar and the acceleration from it. The bottom panel shows the fraction of bar acceleration. The degeneracy can be broken when additional information such as bar length is used. We used the bar length value from Wegg et al. (2015). The contours show 1σ and 2σ confidence intervals for the Anders et al. (2019) (SH) and Schönrich et al. (2019) (Sc) datasets.

Open in new tabDownload slide

Currently we have only used the SN region to disentangle the bar parameters; hence, we were able to determine only some of the degenerate values. We were also able to determine that the sum of the radial acceleration due to the bar and axisymmetric components is constant, as seen in Fig. 5. Hence, if we know one then we can easily estimate the other.

5 DISCUSSION AND CONCLUSION

5.1 Validity of the sample construction

To test how well the approach described in Section 2.3 is able to cope with flux-limited data, we applied our model to two sets: a flux-limited sample and a volume-limited sample. The volume-limited approach was tested with simulation data in Kipper et al. (2019a) and was found to be consistent. Since the results for the flux-limited sample agreed well with those of the volume-limited sample, we are confident that our model is very suitable for the current case.

The acceleration components used for this test are in the axisymmetric form of equations (4)–(5). The geometry of the region was biaxial ellipsoid in the form of equation (12), but its selection and modelling had some differences: the sample was limited by the absolute J-magnitude value of 1.2^m, yielding |$54\, 819$| stars. The grid was constructed by using 70 random sample points, and fitting was done eight times to include the uncertainty from the gridding randomness.

The results of the test calculations for volume- and flux-limited selections from the Sc catalogue are given in Table A2, where calculated acceleration components are given with labels ‘Sc, vol’ and ‘Sc, flux’. The most interesting acceleration components are radial acceleration a_R and vertical acceleration a_z (see their main parameters A_R and A_z). For the volume-limited sample |$A_R = -6128 \pm 199\, {\rm km^2\,s^{-2}\, kpc^{-1}}$| and |$A_z = 183 \pm 106\, {\rm km^2\,s^{-2}\, kpc^{-1}}$|⁠; for the flux-limited sample |$A_R = -6181 \pm 82\, {\rm km^2\,s^{-2}\, kpc^{-1}}$| and |$A_z = 203 \pm 48\, {\rm km^2\,s^{-2}\, kpc^{-1}}$|⁠. The smaller errors in the flux-limited sample are due to the larger data sample. The results are clearly consistent and we may conclude that our approach to cope from here on with the flux-limited sample, described in Section 2.3, is valid.

5.2 Time dependence of acceleration due to the bar

During calculations of stellar orbits (see selected analytical forms for acceleration) we assume that accelerations do not have an explicit time dependence. However, it is known that about a quarter of galaxies contain a more or less prominent bar (Cheung et al. 2013); in the case of the MW a central bar was introduced by de Vaucouleurs (1964). A rotating bar would violate this assumption of our modelling.

To test how much a bar would influence our results, we used the same simulation (the barred one) from Garbari, Read & Lake (2011) as was used in Kipper et al. (2019a). We selected a region close to the solar radius, with an angle between the major axis of the bar and direction to the centre of the region ≈30° and fitted the acceleration components (7)–(9) (Taylor expansion) with our model. We expect that including the tangential component of the acceleration due to the bar that is changing in time would give us a somewhat wrong acceleration direction. We found that the acceleration vector was directed away by 3.27 ± 2.23° from the Galactic centre. The corresponding true angle calculated directly from the simulation gravitational potential was 2.21°. The difference between the true and calculated values is smaller than the 1σ error of the calculated value. Hence the effect of the time dependence of the bar in this case is not so significant.

5.3 Influence of uncertainties in input data

The input to this modelling does not include uncertainties. The resulting uncertainties are statistical in nature and include only sampling errors seen from the likelihood equation (7) of Kipper et al. (2019a). In order to see how observational uncertainties influence our results, we randomized phase space coordinates of stars according to their uncertainties, reconstructed the selection sample as described in Section 3, and remodelled the selected region. To account for the randomness in this process, we modelled the SN 47 times and combined the corresponding posterior distributions. The results of the calculations are shown in Table A2 with the label ‘Sc, rnd’ after the variable name. Comparing calculated accelerations with labels ‘Sc, rnd’ and ‘Sc, flux’, it is seen that randomization had very little effect on the results. We conclude that uncertainties can be ignored for this selection.

Another source of error can be due to the gridding approach employed in this study (see Section 2.2). Since there is randomization, we must include the noise caused by it. We rerun each modelling eight times to include the source of noise. All of these runs had similar posterior distributions; hence we are certain that gridding does not introduce large artificial uncertainties. To include the gridding uncertainties, we combined the posterior distributions of randomized grid runs.

5.4 Conclusions

In this paper, we have applied the orbital arc method (Kipper et al. 2019a) to Gaia DR2 and modelled the acceleration along the plane of the Galactic disc. We approximated the acceleration in the solar neighbourhood with various functional forms and came to the following conclusions:

• There are very few systematic biases between the Gaia DR2 datasets compiled by Schönrich et al. (2019) and Anders et al. (2019). Both the datasets give consistent results.

• The distribution of axisymmetric gravitational acceleration does not account for the observed acceleration for the standard distance of the Sun from the Galactic centre |$\mathrm{R}_\odot \approx 8.3\, {\rm kpc}$|⁠. The curvature of the isopotential lines is smaller than the standard R_⊙, implying that there is a component of the Galactic bar causing this acceleration.

• The acceleration vector in the solar neighbourhood is not directed towards the centre of the Galaxy. There is a significant component of the acceleration directed away from the Galactic centre. We propose that this is caused by the massive central bar. Based on our model, we calculate the torque to be |$\sim 2400\, {\rm km^2\, s^{-2}}$| per solar mass.

• Based on the assumption that isopotential surfaces of the bar are confocal ellipses, we estimate that about a third of the total acceleration in the solar neighbourhood is caused by the bar. In this computation we use the estimate of the length of the bar of Wegg et al. (2015).

• Finally, using our model, we estimated the mass of the bar as |$(1.6\pm 0.3)\times 10^{10}\, \mathrm{M}_\odot$|⁠, using the density distribution parameters from Wegg et al. (2015).

ACKNOWLEDGEMENTS

We thank the referee for helpful comments and suggestions. We thank the StarHorse core team (F. Anders, A. Queiroz, B. Santiago, A. Kalathyan, C. Chiappini) for providing their data, and G. Monari for helpful comments about the paper. This work was supported by institutional research funding IUT26-2, IUT40-2 and PUTJD907 of the Estonian Ministry of Education and Research. We acknowledge the support by the Centre of Excellence ‘Dark Side of the Universe’ (TK133) and by the grant MOBTP86 financed by the European Union through the European Regional Development Fund. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

Footnotes

REFERENCES

APPENDIX A: TABLES