Quantifying torque from the Milky Way bar using Gaia DR2

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 et al., 2018). Here we apply it for the Gaia DR2 data. The orbital arc method has five important steps.

1. Step 1 – Acceleration field:

   We first select a region with a sufficiently large number of stars and known phase space coordinates. Next, we guess an acceleration field and use this to get the orbits. In the orbital arc method, the acceleration field is a free function of the model and contains free parameters. For instance, we can take advantage of the axisymmetric property of a galaxy, and choose an acceleration field described by the cylindrical coordinates:

   	$\displaystyle a_{x}$	$\displaystyle=$	$\displaystyle a_{R}\cos\theta,$		(1)
   	$\displaystyle a_{y}$	$\displaystyle=$	$\displaystyle a_{R}\sin\theta+A_{y},$		(2)
   	$\displaystyle a_{z}$	$\displaystyle=$	$\displaystyle a_{z}.$		(3)

   Components $a_{R}$ and $a_{z}$ are taken in the form of functions

   	$\displaystyle a_{z}$	$\displaystyle=$	$\displaystyle A_{z}+A_{z,z}z+A_{z,R}\Delta R+A_{z,Rz}z\Delta R,$		(4)
   	$\displaystyle a_{R}$	$\displaystyle=$	$\displaystyle A_{R}+A_{R,R}\Delta R,$		(5)
   	$\displaystyle\Delta R$	$\displaystyle=$	$\displaystyle R-R_{\odot},$		(6)

   where $A_{y}$, $A_{z}$, $A_{R}$, $A_{z,z}$, $A_{R,R}$, $A_{z,R}$, $A_{z,Rz}$ and in some cases also $R_{\odot}$ are free parameters obtained via fitting.

   If we do not assume axisymmetry, acceleration vector components are taken in the form of their first-order Taylor expansion:

   	$\displaystyle a_{x}=A_{x}+A_{x,x}\Delta x+A_{x,y}\Delta y+A_{x,z}\Delta z$		(7)
   	$\displaystyle a_{y}=A_{x}+A_{y,x}\Delta x+A_{y,y}\Delta y+A_{y,z}\Delta z$		(8)
   	$\displaystyle a_{z}=A_{x}+A_{z,x}\Delta x+A_{z,y}\Delta y+A_{z,z}\Delta z.$		(9)

   Here $\Delta x$, $\Delta y$, $\Delta z$ denote the distances from the region’s centre.

2. 1.

   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.

3. 2.

   Step 3 – Randomizing star position: The core of the proposed method lies in the oPDF (orbital Probability Density Function), according to which the time of observing a star is random. This means, 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.

4. 3.

   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 $\Delta t_{i}$. So, the probability of finding that star in the $i^{\mathrm{th}}$ Voronoi cell is the time spent in that Voronoi cell, $\Delta t_{i}$, divided by the time spent in the entire region, which is, ${\Delta t_{i}}/{\left(t_{\rm exit}-t_{\rm enter}\right)}$ 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.

5. 4.

   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.

   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.

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 Sect. 2.1. In Fig. 1, $\Delta t$ represents the time a star spends in a cell. The ratio of the time spent by a star in a Voronoi cell, $\Delta t_{i}$, to the time that it spends in the entire region (see $t_{\rm enter}$ and $t_{\rm 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¹¹1The Voronoi tessellation of the region is done in order to compare the original and the new phase space distributions. 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 $\sim 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.

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.

A volume-limited selection is one in which both the stellar distance from an observer and absolute magnitudes of its stars are constrained by a flux limit of the sample. This grants that all of the stars would remain observable if we randomize their position in the region. Our aim is to combine volume-limited selections to acquire methodology that allows us to use flux-limited data. Let us denote $m_{\rm lim}$ as the completeness limit of the flux-limited sample. Then the corresponding absolute-magnitude limit $M_{\rm lim}$ and the distance limit $d_{\rm lim}$ are related by $5\log_{10}d_{\rm lim}=m_{\rm lim}-M_{\rm lim}+5$ (at the moment we ignore the attenuation correction). It is possible to construct a volume-limited sample by selecting only stars that have $M<M_{\rm lim}$ and $d<d_{\rm lim}$. The same applies if we use an additional cut from higher absolute magnitudes, leaving $M$ in the range $M_{\rm lim}-\Delta M<M\leq M_{\rm lim}$. Following the denotations from Kipper et al. (2019a), the observed phase space density as a function of phase space coordinates (${\bf q}$) for a volume-limited sample can now be written as $p_{\rm obs}({\bf q}|M_{\rm lim},d_{\rm lim})$, and the model one as $p({\bf q}|M_{\rm lim},d_{\rm lim},\zeta)$. Here $M_{\rm lim}$ and $d_{\rm lim}$ are not independent, but tied to the absolute-magnitude definition. For the correct gravitational acceleration parameters $\zeta$, and irrespective of the absolute-magnitude limit, these distributions must match:

If this relation matches for each small volume-limited sample, then the relation must hold for the sum (or integral) of these small volume-limited samples as well:

This means that we may integrate an orbit until the apparent magnitude of the corresponding star reaches the limiting magnitude $m_{\rm lim}$ due to its changing distance, and smooth the position of the star along its orbit within that limit. In Sect. 5.1 we test the validity of this approach.

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. $\Delta x\sim\Delta x_{0}+t\Delta v$, where $\Delta$ 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.

Six-dimensional high-quality phase space coordinates in the SN are now available from 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 catalog, SH, (Anders et al., 2019) and the Schönrich catalog, 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 6 964 515 entries; between the Schönrich et al. (2019) Sc and 2MASS catalogues there were 6 519 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_{\rm lim}|J<J_{\rm lim})\ll 1$. The apparent-magnitude data within $0.5~{}{\rm 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_{\rm 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_{\rm lim}$. This is shown as a green line in the left-hand panel and the corresponding probability density distribution $p(G|J<J_{\rm lim}){\rm d}G$ on the right-hand panel. The fraction of $G$ magnitudes crossing $G_{\rm lim}$ is $8\times 10^{-4}$ for the adopted $J_{\rm lim}=10.25$; hence we conclude that our sample is almost independent of the Gaia completeness limit and dust attenuation.

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.

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.

The region size is selected to balance two previous effects: maximally small to have a simple acceleration form and maximally close to keep the observational uncertainties low, and at the same time maximally large to let acceleration act for a sufficiently long time in the model. The boundary of the selected region is described by a biaxial ellipsoid:

with $x_{\rm max}=y_{\rm max}=0.494\,{\rm kpc}$ and $z_{\rm max}=0.218\,{\rm kpc}$. Here $\Delta x$, $\Delta y$ and $\Delta z$ denote the coordinates from the centre of the region. The centre of the region is at $(x,y,z)=(-8.3,0.0,0.0)$ in kpc. The position of the Sun with respect to the region centre is $(-0.040,0.0,0.027)$ in kpc. Within this region there are $417\,727$ stars when using the Schönrich et al. (2019) catalogue (Sc), and $426\,767$ stars when using the StarHorse (SH) catalogue. Larger regions would require the use of a precise selection function, which would complicate the analysis.

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 Eqs. (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 Sect. 3.2) for both the Sc and SH catalogues of Gaia DR2 (for more details see Sect. 3.1).

We used $100$ random points to describe the grid; hence, there are $\approx 42$ stars per grid bin. The fitting was done with the Multinest code (Feroz & Hobson, 2008; Feroz et al., 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 1. In the table we give the posterior distribution of each parameter $\zeta$ with five quantiles positioned at $P(\zeta)=\{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 1 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 $z_{\odot}\approx A_{z}/A_{z,z}=\{-111^{+33}_{-76}({\rm SH}),\,-117^{+58}_{-66}({\rm Sc})\}$ pc from the vertical coordinate value defined as having zero vertical acceleration in contrast to the symmetry-defined centre.

The Poisson equation combines the gravitational potential and total mass density. By using calculated acceleration components in the Poisson equation, we can compute the average total matter density in our selected region as:

where $\Phi$ is the gravitational potential and $\rho_{\rm total}$ is the total mass density inside the region. Calculated Taylor expansion fit components for two different Gaia catalogues give us $\rho_{\rm total}=0.070^{+0.016}_{-0.016}{\rm\,M_{\odot}\,pc^{-3}}$ (SH catalogue) and $\rho_{\rm total}=0.069^{+0.016}_{-0.023}{\rm\,M_{\odot}\,pc^{-3}}$ (Sc catalogue). The full probability density distribution of the total matter density $\rho_{\rm total}$ can be seen in Fig. 3. One must bear in mind, that this total density value applies as the average in this region (extent in the $z$-direction is $2\times 0.22$ kpc). In order to describe the changes in the vertical component of the acceleration, we need a more sophisticated form of acceleration as the linear form cannot capture quick density changes along the $z$-direction. Therefore, if we use another form by assuming axisymmetry with respect to the centre, then the number of free parameters can be significantly reduced and more concrete conclusions can be made. Selected Taylor expansion might not fully grasp all the details of the vertical structure, since it is described with just one free parameter ($A_{z,z}$).

In case of a stationary axisymmetric MW, the acceleration component along the direction of Galactic rotation $a_{y}(\Delta x,\Delta y=0,\Delta 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\sigma$ and $2\sigma$ limits are shown in Fig. 4 and the used priors are given in Table 2. 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.

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(R_{\odot}~{}>~{}8.3\,{\rm kpc})~{}=~{}\{0.048({\rm Sc}),\,0.11({\rm SH})\}$. 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 Sect. 5.1, we calculated the acceleration components assuming axisymmetry, by selecting the components $a_{R},a_{z}$ to be in the form of Eqs. (4) – (6). During the fitting the solar distance $R_{\odot}$ was also taken as a free parameter. Taking the posterior in these fits close to $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^-1. 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 errors.

To calculate the total mass of the bar, we assumed that spatial density distribution of the bar has the same form as that derived by Wegg et al. (2015):

The values of the parameters $x_{0},y_{0},z_{0},\sigma_{\rm in},\sigma_{\rm out},c_{\perp},R_{\rm in},R_{\rm out}$ were also taken to be the same as those derived by Wegg et al. (2015). In this form Wegg et al. (2015) excluded symmetric parts of the Galaxy (e.g. bulge) by cut-off. By calculating tangential accelerations for this bar and fitting the calculated values with the values derived by us and referred to in the previous subsection, we derived that the mass of the bar (using the cut-off in previous equations) will be $0.41^{+0.07}_{-0.08}10^{10}\mathrm{M}_{\odot}$ for the SH catalogue and $0.40^{+0.07}_{-0.11}10^{10}\mathrm{M}_{\odot}$ for the Sc catalogue. Without the cut-off, by adopting their profile (and inferring their $M_{\rm bar}$) the overall bar mass would be $\{1.59^{+0.27}_{-0.31}({\rm SH}),\,1.55^{+0.29}_{-0.43}({\rm Sc})\}\,10^{10}\mathrm{M}_{\odot}$.

In the previous subsection we concluded that the assumption that equipotential curves are circles is clearly not valid. Therefore, we assumed that these curves are confocal ellipses in the Galactic plane and then computed the acceleration due to the bar. An analytical form for accelerations describing the potential of the bar as confocal ellipsoids was chosen to be

where $a_{R}$ and $\Delta R$ are given by Eqs. (5) and (6). The normalized vector components of the potential gradient of the confocal bar, $\Delta x^{\prime}$ and $\Delta y^{\prime}$, are described by

$L_{\rm bar}$ is the focal length of the equipotential curves and $a$ describes the size of the ellipsoid. The coordinate transformation from $(x,y)$ to $(x^{\prime},y^{\prime})$ is done by rotating the original axes by an angle of $29.5^{\circ}$ which is the position angle of the major axis of the bar (Wegg et al., 2015). Results from these calculations are given in Table 3, they contain the coefficients obtained from the fits.

When using accelerations in the form of Eqs. (16)–(18), substantial correlations exist in the modelled posterior samples (e.g. between $A_{r}$ and $A_{r,{\rm bar}}$). The largest correlation coefficient was found to be $1.0$ between $A_{r}$ and $A_{r,{\rm 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,{\rm bar}}$ and the radial acceleration derivative $A_{R,R}$ (correlation coefficient $0.26$).

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.

Based on the results from Sect. 4.2, we modelled the tangential acceleration value $A_{y}$. From this, we can calculate the $z$-component of the torque caused by the bar per solar mass using:

or

Here $R_{\rm GC}$ is our physical distance from the Galactic centre (instead of the modelled radius of curvature of equipotential curves $R_{\odot}$). The degeneracy is because a long bar closer to us and a short massive bar engender the same force, making it difficult to distinguish the two scenarios. This is seen in Fig. 6 where the acceleration value due to the bar and the fraction of acceleration caused by the bar as a function of the length of the bar are given. The degeneracy can be broken when we fix the length of the bar from independent measurements. As an example, Wegg et al. (2015) estimated that the half-length of the bar is $4.6\pm 0.3\,{\rm kpc}$. If we fix this value as $L_{\rm bar}$, it gives the fraction of acceleration caused by the bar as about a third: $0.34\pm 0.07$ in the case of the SH catalogue and $0.29\pm 0.09$ in the case of the Sc catalogue.

To test how well the approach described in Sect. 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 Eqs. (4) – (5). The geometry of the region was biaxial ellipsoid in the form of Eq. (12), but its selection and modelling had some differences: the sample was limited by the absolute $J$-magnitude value of $1.2^{\rm 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 2, 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 Sect. 2.3, is valid.

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 et al. (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 $\approx 30^{\circ}$ 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 which is changing in time would give us a somewhat wrong acceleration direction. We found that the acceleration vector was directed away by $3.27\pm 2.23^{\circ}$ from the Galactic centre. The corresponding true angle calculated directly from the simulation gravitational potential was $2.21^{\circ}$. The difference between the true and calculated values is smaller than the $1\sigma$ error of the calculated value. Hence the effect of the time dependence of the bar in this case is not so significant.

Another approach to estimate the effect of a time dependent gravitational potential is to use available data from the literature. The average angular speed of the bar is about $\sim 40\,{\rm km\,s^{-1}kpc^{-1}}$, although there is significant uncertainty in this value (see Sect. 1). The angular speed of the Sun is $\simeq 30\,{\rm km\,s^{-1}kpc^{-1}}$(Bland-Hawthorn & Gerhard, 2016). The strong similarity between these two values suggests that the potential of the bar near the Sun does not change very fast. A hypothetical star moving with a speed of $220\,{\rm km\,s^{-1}}$ passes the half box-size distance of $0.5$ kpc within $\sim 2.3$ Myr. Within this time, the angle of the bar orientation with respect to our comoving location changes only by $2.3\,\mathrm{Myr}\,\times\,(40-30)\,\mathrm{km\,s^{-1}kpc^{-1}}\approx 1.4^{\circ}$. If we assume that the angle between us and the major axis of the bar is $\theta_{0}=30^{\circ}$ and the ’tip of the bar’ is about $L=5$ kpc away from the centre of the Galaxy²²2We make an approximation that the mass of the bar is a point mass at the tip of the bar. This gives an upper limit for the bar influence. then our distance to the tip of the bar is

The force due to the bar changes within $2.3$ Myr maximally by

Hence, the force due to the bar changes by about $5$ per cent if the force due to the bar is not steeper than $\propto w^{-2}$ and the bar dominates the potential. If it does not, then the result must be multiplied by the acceleration fraction of the bar. This can be considered as a component of systematic uncertainty. While calculating the orbital arcs for stars within the selected region, the time dependence of accelerations due to the bar has quite a small effect. Therefore in the current study we ignore this effect.

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 Sect. 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 2 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 Sect. 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.

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:

1. 1.

   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.

2. 2.

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

3. 3.

   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.

4. 4.

   Based on the assumption that isopotential surfaces of the bar are confocal ellipses, we estimate that about one 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).

5. 5.

   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).