Kinematics of the Milky Way from the Statistical Analysis of the Gaia Data Release 3

1. Introduction

Our galaxy, the Milky Way (MW), is a unique laboratory for gravity research and for understanding the formation and evolution of galaxies. In recent years, the Gaia Space Observatory has acquired a huge amount of precise astrometric, photometric and spectroscopic data on stars in the MW [23].The analysis of these data has been the subject of many thousands of publications.

The Gaia DR3 catalogue provides very rich input data for creating a kinematic map of the MW, much more detailed than previously possible. In general, it encompasses various structures on different scales, from orbiting of small gravitationally bound systems, binaries and multiple-bound systems, to the streaming motions of stellar fields in spiral arms [3,5] with various turbulences and fluctuations, to the collective rotation of the whole galactic disk with the galactic halo.

The nature of the rotation suggests the presence of dark matter, which generates a substantial part of the MW mass, see [1,8,15,16,18,28,30,37,39,48] and the overview [42]. Recent studies of MW kinematics have shown accurate results on the MW rotation represented by the rotation curve (RC) defined as the dependence of the orbital velocity on the radius. These curves are the basic input for models accounting for the presence of dark matter. This is why these two topics often appear together in the papers.

Along with gravity, the formation and evolution of the stars themselves are also controlled by the forces of the microworld (strong, electroweak - unified electromagnetic+weak) based on a well-verified Standard model. The nature and origin of dark matter at the microscopic level have not yet been explained.

Very important topics concern the detailed mapping of many kinematic substructures beyond axial symmetry [20,24,31,32,34,36,38,44,46,47]. An overview of the integrated, structural and kinematic parameters of the Galaxy is given in [10].

The first goal of our study is to use new statistical methods to determine the local velocity of the Sun and to analyse the collective orbital motion of stars in the galactic plane and beyond. Because these methods do not require knowledge of radial velocities, we can access larger samples of stars [4,35]. The second goal is to obtain the 3D axisymmetric approximation of the kinematic MW image generated by averaging local asymmetric substructures. This image can be useful in determining the scale of local asymmetries or as input for validating axisymmetric models.

In Section 2 we describe our methodology and define the basic concepts and quantities we will work with. First, we define transformation relations between galactic and Galactocentric reference frames. The first system is used to acquire and present Gaia data, while the second is suitable for simulation and data interpretation. Then the definition of the axisymmetric Monte-Carlo model in the Galactocentric reference frame follows. The model is based on a triple (partially asymmetric) Gaussian distribution, which depends on the distance from the galactic plane and is defined by six free parameters. Next, we define the sky sectors to be analysed and the cuts for selecting accurate data. In Section 3 we present our methods together with the results obtained from velocity distributions in the different sectors of the sky. Analysis of these distributions provides precise results concerning the local and orbital velocity of the Sun (Section 3.1), the different RC representations (Section 3.2) and the setting of the free parameters of the simulation model (Section 3.3). This is followed by a comparison of the simulation with all relevant distributions (Section 3.4). The obtained results and their agreement with the simulation model are discussed in more detail in Section 4. The discussion includes a comparison of the results obtained with other available data.

2. Methodology and definitions

2.1. Reference Frames

Positions of sources in Gaia data are represented in angular galactic coordinates: longitude (l) and latitude (b). Using parallax, we can also define the distance r of the source from the Sun. For our analysis also representation in Galactocentric coordinates will be important. The relation between both reference frames is illustrated in Figure 1.

The axes of the galactic and Galactocentric reference frames are defined as

$$\begin{array}{cc}\hfill x& =rcoslcosb,y=rsinlcosb,z=rsinb,\hfill \\ \hfill X& =Rcos\Phi cos\Theta ,Y=Rsin\Phi cos\Theta ,Z=Rsin\Theta .\hfill \end{array}$$

(1)

For clarity, we will use the following convention throughout the paper: positions, velocities, and their coordinates are denoted in lowercase (uppercase) in the galactic (Galactocentric) reference frame. So, the direction x points to the centre of the galaxy and the direction y is the direction of the MW rotation. The corresponding coordinates are related as

$$X=x-R_{\odot };Y=y;Z=z+Z_{\odot }.$$

(2)

Radius $R_{\odot }$ represents the distance of the Sun from the galactic center. The values obtained in the measurements are in the interval $7.1-8.92$ kpc [10]. For our analysis, we used the recent value $R_{\odot }\approx 8.178$ kpc [26]. The small parameter $Z_{\odot }\approx 0.0208$ kpc [7] defines the Sun position above the galactic plane.

Then, the transformation between spherical coordinates $\left(r,l,b\right)\to \left(R,\Phi ,\Theta \right)$ of both the frames is defined by equations:

$$\begin{array}{c}\hfill R=\sqrt{r^2+R_{\odot }^2+Z_{\odot }^2-2R_{\odot }x+2Z_{\odot }z},\Theta =arcsin{\displaystyle \frac{z+Z_{\odot }}{R}},\\ \hfill \Phi \left(x>R_{\odot }\right)=arcsin{\displaystyle \frac{y}{Rcos\Theta }},\Phi \left(x\le R_{\odot }\right)=\pi -arcsin{\displaystyle \frac{y}{Rcos\Theta }}\end{array}$$

(3)

and inversely $\left(R,\Phi ,\Theta \right)\to$$\left(r,l,b\right)$

$$\begin{array}{c}\hfill r=\sqrt{R^2+R_{\odot }^2+Z_{\odot }^2+2R_{\odot }X-2Z_{\odot }Z},b=arcsin{\displaystyle \frac{Z-Z_{\odot }}{r}}\\ \hfill l\left(X>R_{\odot }\right)=\pi -arcsin{\displaystyle \frac{Y}{rcosb}},l\left(X\le R_{\odot }\right)=arcsin{\displaystyle \frac{Y}{rcosb}},\end{array}$$

(4)

where $x,y,z,X,Y,Z$ are defined by (1). We will need these transformations to model and simulate the motion of the stars in Section 2.2 and Section 3.4. The spherical coordinates $\left(R,\Phi ,\Theta \right)$ are simply related to the cylindrical coordinates $\left(R_c,\Phi ,Z\right)$:

$$R_c=Rcos\Theta ;Z=Rsin\Theta .$$

(5)

In the following discussion, we neglect the $Z_{\odot }$ parameter and explain the validity of this approximation in Section 3.3. Next, we will need vectors forming the local orthonormal bases in both reference frames:

$$\begin{array}{c}\hfill {\mathbf{N}}_{\Phi }=\left(-sin\Phi ,cos\Phi ,0\right),{\mathbf{N}}_{\Theta }=\left(-cos\Phi sin\Theta ,-sin\Phi sin\Theta ,cos\Theta \right),\\ \hfill {\mathbf{N}}_R=\left(cos\Phi cos\Theta ,sin\Phi cos\Theta ,sin\Theta \right),\end{array}$$

(6)

$$\begin{array}{c}\hfill {\mathbf{n}}_l=\left(-sinl,cosl,0\right),{\mathbf{n}}_b=\left(-coslsinb,-sinlsinb,cosb\right),\\ \hfill {\mathbf{n}}_r=\left(-coslcosb,-sinlcosb,sinb\right)\end{array}$$

(7)

where the vectors ${\mathbf{N}}_{\alpha }$, and ${\mathbf{n}}_{\alpha }$ define directions of increasing coordinates $\alpha =\Phi ,\Theta ,R,l,b,r$. For example $-{\mathbf{N}}_{\Phi }$ is direction of MW rotation at any point.

The velocity of a star at the point $\mathbf{R}$ of Galactocentric frame can be split into two components:

$$\mathbf{V}\left(\mathbf{R}\right)={\mathbf{V}}_G\left(\mathbf{R}\right)+\Delta \mathbf{V}\left(\mathbf{R}\right);{\mathbf{V}}_G\left(\mathbf{R}\right)=-V_G(R,Z){\mathbf{N}}_{\Phi },$$

(8)

where $V_G(R,Z)>0$ is the average (over $\Phi$) rotation velocity at the radius R and distance Z from the galactic plane. $\Delta \mathbf{V}$ is the deviation from the average ${\mathbf{V}}_G\left(\mathbf{R}\right)$, which represents the velocity of the local rest frame and can be approximated by the average velocity of the stars in some neighborhood of $\mathbf{R}$,

$${\mathbf{V}}_G\left(\mathbf{R}\right)=〈\mathbf{V}\left(\mathbf{R}\right)〉;〈\Delta \mathbf{V}\left(\mathbf{R}\right)〉=0.$$

(9)

The quality of the approximation may depend on the choice of sources and the size and shape of the defined neighborhood. In Section 3.1 this issue will be discussed in more detail for ${\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)$, which is the velocity of the local standard rest frame (LSR) [40].

The proper motion of the stars in Gaia data is represented by the vector

$${\mu }_{ICRS}=\left({\mu }_{\alpha }^{*},{\mu }_{\delta }\right);{\mu }_{\alpha }^{*}\equiv {\mu }_{\alpha }cos\delta ,$$

(10)

whose components are angular velocities in directions of the right ascension and declination in the ICRS. For our analysis, we will prefer the representation of proper motion in the galactic reference frame

$${\mu }_{gal}=\left({\mu }_l^{*},{\mu }_b\right),{\mu }_l^{*}\equiv {\mu }_{l}cosb,$$

(11)

which is obtained by the transformation given in [25]. Then we will work with transverse 2D velocity defined as

$${\mathbf{v}}_{gal}=\left(v_l,v_b\right)=r\left({\mu }_l^{*},{\mu }_b\right),$$

(12)

where $v_l,v_b$ are velocity components in directions of increasing latitude $\left(l\right)$ and longitude $\left(b\right),$ and the distance r is obtained from the parallax

$$r\left[\mathrm{kpc}\right]={\displaystyle \frac{1}{p\left[\mathrm{mas}\right]-\mathrm{ZP}}};\mathrm{ZP}=-0.017\mathrm{mas},$$

(13)

where ZP is a global zero point taken from [22]. The ZP correction reduces velocity only by $\lesssim 2\%$ in our selected data. The analysis and discussion of distance extraction from parallax is thoroughly done in [6,19,21].

In Section 3.1 we will study the dependence of the mean values $〈v_l〉,〈v_b〉$ on the distance and direction from the Sun. From these curves, we will determine the local velocity of the Sun. In the rest of Section 3 we will work with slightly modified velocities $〈v_l〉,〈v_b〉$ and corresponding dispersions that are related to the LSR frame. In this reference frame, the dependence on the Sun’s local motion is eliminated.

2.2. Simulation Model of Stellar Velocities

The velocity distributions will be compared with a simple probabilistic Monte-Carlo model. The model generates the velocity of a star in the Galactocentric reference frame (Figure 1)

$${\mathbf{V}}^{gen}=\left(V_{\Phi }-V_0\left(R\right)\right){\mathbf{N}}_{\Phi }+V_{\Theta }{\mathbf{N}}_{\Theta }+V_R{\mathbf{N}}_R,$$

(14)

where $V_{\Phi },V_{\Theta }$ and $V_R$ are its components generated in the local reference frame defined by the basis (6). Velocity $V_0\left(R\right)>0$ is defined as an average of orbital velocity at the galactic plane and radius R

$$V_0\left(R\right)=V_G(R,0),$$

(15)

which is also our definition of the RC. This definition is based on direct measurements of orbital velocities in selected sectors of the MW disk, so the results obtained may differ from a global RC calculated from Jeans modelling [29] assuming an axisymmetric gravitational potential of the MW. Our definition reflects the collective orbital velocity rather than the velocity of a single star or a test particle [18].

As we shall see, a very good agreement with the observed curves and distributions provides the simulation model based on the multinormal distribution

$$P(V_{\Phi },V_{\Theta },V_R;R,Z)\sim exp\left(-{\displaystyle \frac{V_{\Phi }^2}{2{\sigma }_{\Phi }^2}}-{\displaystyle \frac{V_{\Theta }^2}{2{\sigma }_{\Theta }^2}}-{\displaystyle \frac{V_R^2}{2{\sigma }_R^2}}\right),$$

(16)

with partial asymmetry - different ${\sigma }_{\Phi }$ for positive and negative $V_{\Phi }$. The dependence on R and Z is absorbed in the standard deviations ${\sigma }_{\alpha }(R,Z)$. This important dependence will be discussed in more detail in Section 3.3. For comparison with the data, the velocities simulated in the Galactocentric system will be transformed to the galactic system as described in Section 3.4.

2.3. Data Set

Similarly as in our previous study [50,51,52], the field of stars is broken down into a mosaic of small square cells that represent a statistical input for our analysis (Figure 2).

Here the reason is rather technical, the cells allow us to exclude regions with a very high or inhomogeneous stellar density. The data sectors of the sky used for analysis are defined in Table 1.

Table 1. MW sectors used for analysis. Units: b,l [deg].

b∖l	$〈-5,5〉$	$〈85,95〉$	$〈175,185〉$	$〈265,275〉$	$〈45,135〉$	$〈225,315〉$
$〈-5,5〉$	A	B	C	D	${\mathrm{Q}}_2$	${\mathrm{Q}}_4$
			b∖l	$〈0,360〉$
			$〈-90,-80〉$	E
			$〈80,90〉$	F
	b∖l	$〈-45,45〉$	$〈45,135〉$	$〈135,225〉$	$〈225,315〉$
	$〈-60,-45〉$	${\mathrm{Q}}_{1S}$	${\mathrm{Q}}_{2S}$	${\mathrm{Q}}_{3S}$	${\mathrm{Q}}_{4S}$
	$〈45,60〉$	${\mathrm{Q}}_{1N}$	${\mathrm{Q}}_{2N}$	${\mathrm{Q}}_{3N}$	${\mathrm{Q}}_{4N}$

Table 1. MW sectors used for analysis. Units: b,l [deg].

b∖l	$〈-5,5〉$	$〈85,95〉$	$〈175,185〉$	$〈265,275〉$	$〈45,135〉$	$〈225,315〉$
$〈-5,5〉$	A	B	C	D	${\mathrm{Q}}_2$	${\mathrm{Q}}_4$
			b∖l	$〈0,360〉$
			$〈-90,-80〉$	E
			$〈80,90〉$	F
	b∖l	$〈-45,45〉$	$〈45,135〉$	$〈135,225〉$	$〈225,315〉$
	$〈-60,-45〉$	${\mathrm{Q}}_{1S}$	${\mathrm{Q}}_{2S}$	${\mathrm{Q}}_{3S}$	${\mathrm{Q}}_{4S}$
	$〈45,60〉$	${\mathrm{Q}}_{1N}$	${\mathrm{Q}}_{2N}$	${\mathrm{Q}}_{3N}$	${\mathrm{Q}}_{4N}$

The narrow sectors A-F are defined by basic directions (from the Sun, see Figure 1): towards the centre of the galaxy (A), along the direction of galaxy rotation (B), and their opposites (C,D). The perpendiculars to the galactic disk define sectors E and F. The analysis is extended by other sectors: Q_α - in the galactic disk and ${\mathrm{Q}}_{\alpha \beta }$ - outside the disk. Dimensions of the square cells are defined in Table 2 and depend on the latitude of the sector. As a result, we have approximately the same average number of stars per cell in all the sectors. Only cells with a limited number M stars are accepted in the analysis, from which we exclude stars with negative parallax. We accept:

$$M\le 50,p>0.$$

(17)

After this selection, we accept only sources that meet the condition

$${\displaystyle \frac{\Delta v_{gal}}{v_{gal}}}\equiv \eta \le 0.3;v_{gal}=\left|{\mathbf{v}}_{gal}\right|.$$

(18)

In Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8, from which we determine our resulting velocities in Table 3, we use a more strict cut: $\eta \le 0.1$. This ratio can be estimated using (12), (13) and parallax and proper motion errors in the Gaia data as

$$\begin{array}{c}\hfill {\displaystyle \frac{\Delta v_{gal}}{v_{gal}}}\approx \sqrt{{\left({\displaystyle \frac{\Delta p}{p}}\right)}^2+{\displaystyle \frac{{\left({\mu }_l^{*}\Delta {\mu }_l^{*}\right)}^2+{\left({\mu }_b\Delta {\mu }_b\right)}^2}{{\mu }^4}}};\mu =\sqrt{{\mu }_l^{*2}+{\mu }_b^2}.\end{array}$$

(19)

The resulting numbers of sources in the respective sectors are shown in the same table. Most of our calculations focus on mean values, which means that the resulting errors are much smaller than the errors of individual entries.

Table 2. Dimensions of the cells and numbers of stars in MW sectors after cuts (17), (18).

	A	B	C	D	E	F	${\mathrm{Q}}_2$	${\mathrm{Q}}_4$
$\Delta [deg]$	0.02	0.02	0.02	0.02	0.1	0.1	0.02	0.02
${\mathrm{n}}_{cut}$	332623	1862224	809392	1464150	300037	269708	11561084	10753108
	${\mathrm{Q}}_{1S}$	${\mathrm{Q}}_{2S}$	${\mathrm{Q}}_{3S}$	${\mathrm{Q}}_{4S}$	${\mathrm{Q}}_{1N}$	${\mathrm{Q}}_{2N}$	${\mathrm{Q}}_{3N}$	${\mathrm{Q}}_{4N}$
$\Delta [deg]$	0.05	0.05	0.05	0.05	0.05	0.05	0.05	0.05
${\mathrm{n}}_{cut}$	1578844	1040017	854355	1378185	1281482	1236598
832720	918059

Table 2. Dimensions of the cells and numbers of stars in MW sectors after cuts (17), (18).

	A	B	C	D	E	F	${\mathrm{Q}}_2$	${\mathrm{Q}}_4$
$\Delta [deg]$	0.02	0.02	0.02	0.02	0.1	0.1	0.02	0.02
${\mathrm{n}}_{cut}$	332623	1862224	809392	1464150	300037	269708	11561084	10753108
	${\mathrm{Q}}_{1S}$	${\mathrm{Q}}_{2S}$	${\mathrm{Q}}_{3S}$	${\mathrm{Q}}_{4S}$	${\mathrm{Q}}_{1N}$	${\mathrm{Q}}_{2N}$	${\mathrm{Q}}_{3N}$	${\mathrm{Q}}_{4N}$
$\Delta [deg]$	0.05	0.05	0.05	0.05	0.05	0.05	0.05	0.05
${\mathrm{n}}_{cut}$	1578844	1040017	854355	1378185	1281482	1236598
832720	918059

3. Methods and Results

In Figure 3 we show the distribution of star distances in the data sectors A-F defined above. The distance of most of them is $r\lesssim 3-6$ kpc. Dependencies of mean velocities $〈v_l〉,〈v_b〉,〈v_{gal}〉$ and related standard deviations on distance r are shown in the figures that follow. What information can be extracted from them?

Figure 3. Distribution of distances in the sectors A-F. Unit: r[kpc]. Binning: 0.024kpc.

Figure 3. Distribution of distances in the sectors A-F. Unit: r[kpc]. Binning: 0.024kpc.

3.1. Velocity of the Sun

According to (8) and (15) the Sun’s velocity can be decomposed as

$${\mathbf{V}}_{\odot }={\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)+\Delta {\mathbf{V}}_{\odot };{\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)=-V_0\left(R_{\odot }\right){\mathbf{N}}_{\Phi },$$

(20)

where ${\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)$ is the LSR velocity and $\Delta {\mathbf{V}}_{\odot }$ is the Sun’s velocity in the LSR frame. How can these velocities be determined? In the galactic reference frame the star velocity $\mathbf{v}\left(\mathbf{r}\right)$ is given as

$$\mathbf{v}\left(\mathbf{r}\right)=\mathbf{V}\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)-{\mathbf{V}}_{\odot }={\mathbf{V}}_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)+\Delta \mathbf{V}\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)-{\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)-\Delta {\mathbf{V}}_{\odot },$$

(21)

which implies

$$〈〈\mathbf{v}\left(\mathbf{r}\right)〉〉=〈〈\mathbf{V}\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)〉〉-{\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)-\Delta {\mathbf{V}}_{\odot }.$$

(22)

The region of averaging $〈〈...〉〉$$\left(\mathrm{over}\mathbf{r}\right)$ should be reasonably wide to suppress the influence of local fluctuations on the obtained average, then the first two terms cancel and for the local Sun’s velocity we obtain

$$\Delta {\mathbf{V}}_{\odot }=-〈〈\mathbf{v}\left(\mathbf{r}\right)〉〉.$$

(23)

For its calculation, we use the projections of galactic velocities (12) in sectors A-D shown in Figure 4 together with the relations

$$v_l\left(\mathbf{r}\right)=\mathbf{v}\left(\mathbf{r}\right).{\mathbf{n}}_l,v_b\left(\mathbf{r}\right)=\mathbf{v}\left(\mathbf{r}\right).{\mathbf{n}}_b,$$

(24)

where the vectors of local orthonormal basis ${\mathbf{n}}_{\alpha }$ are defined in (7). In the sectors considered, we get:

Figure 4. Dependence of mean velocity ($v_l$ - blue, $v_b$ - green, $v_{gal}$ - red) on distance r in the galactic reference frame in sectors A-F. Cut $\eta \le 0.1$. Units: r[kpc], v[km/s].

Figure 4. Dependence of mean velocity ($v_l$ - blue, $v_b$ - green, $v_{gal}$ - red) on distance r in the galactic reference frame in sectors A-F. Cut $\eta \le 0.1$. Units: r[kpc], v[km/s].

1) $\Delta V_{z\odot }$ in the sectors A-D Since in these sectors we have ${\mathbf{n}}_b=\left(0,0,1\right),$ so we can identify $v_z=v_b$. The mean values $〈v_b〉$ depending on distance r are for individual sectors shown in Figure 4. The velocity $\Delta V_{z\odot }$ is the average

$$\Delta V_{z\odot }\left(r_{max}\right)=-{〈〈v_z〉〉}_{\Omega },$$

(25)

where $\Omega$ means the region of averaging, which are the stars in the sectors A-D and radius $r<$$r_{max}$. The resulting curve is shown in Figure 5.

Figure 5. Velocities $\Delta V_{y\odot }$ (blue) and $\Delta V_{z\odot }$ (green) as the functions of ${\mathrm{r}}_{max}$(left) with corresponding statistical errors (right). Cut $\eta \le 0.1$. Units: ${\mathrm{r}}_{max}$[kpc], $\Delta V_{\alpha \odot }$[km/s].

Figure 5. Velocities $\Delta V_{y\odot }$ (blue) and $\Delta V_{z\odot }$ (green) as the functions of ${\mathrm{r}}_{max}$(left) with corresponding statistical errors (right). Cut $\eta \le 0.1$. Units: ${\mathrm{r}}_{max}$[kpc], $\Delta V_{\alpha \odot }$[km/s].

2) $\Delta V_{y\odot }$ in the sectors A and C In these sectors we have ${\mathbf{n}}_l=\left(0,\pm 1,0\right)$, so we identify $v_y=+v_l$ in the A and $v_y=-v_l$ in the C sector. The corresponding curves $〈v_l〉$ are in panels A,C (Figure 4) and curve

$$\Delta V_{y\odot }\left(r_{max}\right)=-{\displaystyle \frac{1}{2}}\left({〈〈v_y〉〉}_A+{〈〈v_y〉〉}_C\right)$$

(26)

is in Figure 5. The region of averaging are sectors A and C up to the distance $r_{max}.$

3) $\Delta V_{x\odot }$ and $V_0\left(R_0\right)$ in the sectors B and D Here we have ${\mathbf{n}}_l=\left(\mp 1,0,0\right)$ so we can identify $v_x=-v_l$ in the sector B and $v_x=+v_l$ in the D sector. In panels B,D in Figure 4 we show corresponding $〈v_l〉$ curves. For small r we have $\Delta V_{x\odot }=-〈v_x〉$, which implies for B,D sectors

$$\Delta V_{x\odot }=\pm 〈v_l〉.$$

(27)

but what is the reason for the steep slope of $〈v_l〉$ as r increases? In the considered sectors, the mean value $〈v_l〉$ includes a significant contribution from the collective rotation velocity ${\mathbf{V}}_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)$ proportional to r, as explained below (Figure 9 and Eq.(43)). Thus, in B,D sectors we have

$$〈v_l^B〉=\Delta V_{x\odot }-{\displaystyle \frac{r}{R}}{V}_0\left(R\right),〈v_l^D〉=-\Delta V_{x\odot }-{\displaystyle \frac{r}{R}}{V}_0\left(R\right);R=\sqrt{r^2+R_{\odot }^2}$$

(28)

and correspondingly

$$〈v_x^B〉=-\Delta V_{x\odot }+{\displaystyle \frac{r}{R}}{V}_0\left(R\right),〈v_x^D〉=-\Delta V_{x\odot }-{\displaystyle \frac{r}{R}}{V}_0\left(R\right).$$

(29)

The two curves $〈v_l〉$ are shown in Figure 6 together with the curve produced by fitting the free parameters $\Delta V_{x\odot }$ and $V_0$ in the range $0<r<3$ kpc. Assuming that rotation velocity is constant in these sectors, $V_0\left(R\right)\approx V_0\left(R_{\odot }\right)$, we obtain very good agreement of the fit to the data. In the next we abbreviate $V_0\left(R_{\odot }\right)$ by $V_0$. The obtained velocities are:

$$\Delta V_{x\odot }^B=9.43,V_0^B=228.05;\Delta V_{x\odot }^D=10.66,V_0^D=228.83.$$

(30)

Admittedly, there is a weak dependency [18,39]

$$V_0\left(R\right)=V_0-\kappa \left(R-R_{\odot }\right),$$

(31)

however the range $0<r<3$ kpc means that $8.178<R<8.711$, which is too small a range for reliable fit involving $\kappa$.

Figure 6. Curves (29): data (blue) and fit (red). Cut $\eta \le 0.1$. Units: r[kpc], $v_x$[km/s].

Figure 6. Curves (29): data (blue) and fit (red). Cut $\eta \le 0.1$. Units: r[kpc], $v_x$[km/s].

In the left panel of Figure 5 we observe nearly constant curves, while the right panel shows dependence of corresponding statistical errors. The range of distances $\left(0-3\mathrm{kpc}\right)$ involves the dominant part of analyzed stars in sectors A-D. To estimate the Sun’s local velocity $\Delta V_{y\odot },\Delta V_{z\odot }$, we take the values in the middle of the range $\left(1.5\mathrm{kpc}\right)$. The velocities $\Delta V_{x\odot }$ and $V_0$ are obtained as averages of the values in (30).

However, for consistent comparison with other analyses, the velocities $\Delta V_{y\odot }$ and $V_0$ can require further correction based on the calculation of asymmetric drift (see [14], Sec.11):

$$V_a=V_c-V_0,$$

(32)

where $V_c$ is the MW circular velocity at $R_{\odot }$ and $V_0$ is a mean orbital velocity of the stars in neigborhood of the Sun. The Sun orbital velocity $V_{\odot \Phi }$ in the Galactocentric frame can be decomposed alternatively as

$$V_{\odot \Phi }=V_0+\Delta V_{y\odot }=V_c+\Delta V_{y\odot }^{corr},$$

(33)

so

$$\Delta V_{y\odot }^{corr}=\Delta V_{y\odot }-V_a,V_c=V_0+V_a.$$

(34)

The study [49] estimates the asymmetric drift $V_a$ around the Sun’s position to be about 6 km/s. A very similar value can be extracted from Figure 3 in the paper [40], where the parameter $〈v_R^2〉$ is replaced by our parameter ${\sigma }_{R0}^2=1089$ km²s⁻², which is calculated in Section 3.3, (Table 4). In the standard notation, we can identify

$$U_{\odot }=\Delta V_{x\odot },V_{\odot }=\Delta V_{y\odot }-V_a,W_{\odot }=\Delta V_{z\odot },V_c=V_0+V_a,$$

(35)

where $\Delta V_{x\odot }$ and $V_0$ are the averages of the corresponding values in (30). We have:

$$\Delta {\mathbf{V}}_{\odot }=\left(\Delta V_{x\odot },\Delta V_{y\odot },\Delta V_{z\odot }\right)=\left(10.05,22.65,7.59\right).$$

(36)

The associated velocities $U_{\odot },V_{\odot },W_{\odot },V_c$ and $V_0$ are shown in Table 3.

Table 3. Local velocity of the Sun, mean rotation velocity $V_0$ and circular velocity $V_c$ along with the results of other analyses. Units: [km/s].

$U_{\odot }$	$V_{\odot }$	$W_{\odot }$	$V_0$	$V_c$	${\Delta }_{syst}$	$V_c+V_{\odot }$	ref.
$10.05\pm 0.05$	$16.65\pm 0.06$	$7.59\pm 0.02$	228	234	$\left(0.7,0.6,0.1,4.0,4.0\right)$	$251.0\pm 5$	this work
$11.1_{-0.75}^{+0.69}$	$12.24\pm 0.47$	$7.{25}_{-0.36}^{+0.37}$			$\left(1,2,0.5,-\right)$		[40]
$9.12\pm 0.10$	$20.80\pm 0.10$	$7.66\pm 0.08$	$230\pm 12$			$250.8$	[11]
$9.58\pm 2.39$	$10.52\pm 1.96$	$7.01\pm 1.67$					[43]
$10.00\pm 0.36$	$5.25\pm 0.62$	$7.17\pm 0.38$					[17]
	$14.6\pm 5$			$240\pm 8$		$255.2\pm 5.1$	[39]
	$26\pm 3$			$218\pm 6$		${242}_{-3}^{+10}$	[12]
$7.7\pm 0.9$	$12.4\pm 0.7$			$236\pm 3$		$248.4$	[32]
	$12.1\pm 7.6$			$240\pm 6$		$252.1$	[28]
$\approx 11.1$		$\approx 7.8$		$229\pm 0.2$	$\left(-,-,-,\approx 2\%-5\%\right)$	$\approx 245.8$	[18]
	$24\pm 1$				$\pm 2\left(V_c\right)\pm 5\left(\mathrm{large}\mathrm{scale}\right)$		[13]
				$238\pm 9$		$250\pm 9$	[41]
				$233.6\pm 2.8$			[37]
			$\approx 230$				[20]
				$224\pm 13$			[33]
				$217\pm 6$			[48]
						$247.4\pm 1.4$	[26]
						$253\pm 6$	[27]

Table 3. Local velocity of the Sun, mean rotation velocity $V_0$ and circular velocity $V_c$ along with the results of other analyses. Units: [km/s].

$U_{\odot }$	$V_{\odot }$	$W_{\odot }$	$V_0$	$V_c$	${\Delta }_{syst}$	$V_c+V_{\odot }$	ref.
$10.05\pm 0.05$	$16.65\pm 0.06$	$7.59\pm 0.02$	228	234	$\left(0.7,0.6,0.1,4.0,4.0\right)$	$251.0\pm 5$	this work
$11.1_{-0.75}^{+0.69}$	$12.24\pm 0.47$	$7.{25}_{-0.36}^{+0.37}$			$\left(1,2,0.5,-\right)$		[40]
$9.12\pm 0.10$	$20.80\pm 0.10$	$7.66\pm 0.08$	$230\pm 12$			$250.8$	[11]
$9.58\pm 2.39$	$10.52\pm 1.96$	$7.01\pm 1.67$					[43]
$10.00\pm 0.36$	$5.25\pm 0.62$	$7.17\pm 0.38$					[17]
	$14.6\pm 5$			$240\pm 8$		$255.2\pm 5.1$	[39]
	$26\pm 3$			$218\pm 6$		${242}_{-3}^{+10}$	[12]
$7.7\pm 0.9$	$12.4\pm 0.7$			$236\pm 3$		$248.4$	[32]
	$12.1\pm 7.6$			$240\pm 6$		$252.1$	[28]
$\approx 11.1$		$\approx 7.8$		$229\pm 0.2$	$\left(-,-,-,\approx 2\%-5\%\right)$	$\approx 245.8$	[18]
	$24\pm 1$				$\pm 2\left(V_c\right)\pm 5\left(\mathrm{large}\mathrm{scale}\right)$		[13]
				$238\pm 9$		$250\pm 9$	[41]
				$233.6\pm 2.8$			[37]
			$\approx 230$				[20]
				$224\pm 13$			[33]
				$217\pm 6$			[48]
						$247.4\pm 1.4$	[26]
						$253\pm 6$	[27]

Their systematic uncertainties are estimated as follows. For velocities $\Delta V_{y\odot },\Delta V_{z\odot }$ are defined by the range of variability of the curves in the left panel in interval $0.5<r_{max}<3$ kpc. For velocities $\Delta V_{x\odot }$ and $V_0$ systematic uncertainty is estimated from a set of fits in different intervals of distances $\left(0<r<2,3,4,5\mathrm{kpc}\right)$. The relatively small systematic uncertainty in $V_0$ reflects the relatively small local variations.

3.2. Rotation Curves

From now, we will substitute galactic velocity $\mathbf{v}\left(\mathbf{r}\right)$ in (21) by $\mathbf{v}\left(\mathbf{r}\right)$$\to \mathbf{v}\left(\mathbf{r}\right)+\Delta {\mathbf{V}}_{\odot }$, which is the velocity related to the LSR:

$$\begin{array}{cc}\hfill \mathbf{v}\left(\mathbf{r}\right)& ={\mathbf{V}}_G\left(\mathbf{R}\right)+\Delta \mathbf{v}\left(\mathbf{R}\right)-{\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill 〈\mathbf{v}\left(\mathbf{r}\right)〉& ={\mathbf{V}}_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)-{\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right).\hfill \end{array}$$

(38)

In this frame, the input data (12) are modified with the use of our $\Delta {\mathbf{V}}_{\odot }$ as

$$\begin{array}{cc}\hfill v_l& \to v_l+\Delta {\mathbf{V}}_{\odot }.{\mathbf{n}}_l,v_b\to v_b+\Delta {\mathbf{V}}_{\odot }.{\mathbf{n}}_b,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill v_{gal}& \to \sqrt{{\left(v_l+\Delta {\mathbf{V}}_{\odot }.{\mathbf{n}}_l\right)}^2+{\left(v_b+\Delta {\mathbf{V}}_{\odot }.{\mathbf{n}}_b\right)}^2}.\hfill \end{array}$$

(40)

After this substitution, Figure 4 is replaced by the top panel in Figure 7.

Figure 7. Dependence of mean velocity ($v_l$ - blue, $v_b$ - green, $v_{gal}$ - red) on distance r in the local rest frame at ${\mathbf{R}}_{\odot }$ in sectors A-F. For panel DATA we used cut $\eta \le 0.1$. For panel SIMULATION see Section 3.3 and Section 3.4. Units: r[kpc], v[km/s].

Figure 7. Dependence of mean velocity ($v_l$ - blue, $v_b$ - green, $v_{gal}$ - red) on distance r in the local rest frame at ${\mathbf{R}}_{\odot }$ in sectors A-F. For panel DATA we used cut $\eta \le 0.1$. For panel SIMULATION see Section 3.3 and Section 3.4. Units: r[kpc], v[km/s].

The new $v_l,v_b,v_{gal}$ refer to the LSR whose velocity is

$${\mathbf{V}}_{LSR}={\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)=\left(0,V_0,0\right).$$

(41)

The combination of the new panels A and C, which represents the RC is shown in Figure 8a.

Figure 8. Velocity curves $V=V_0\left(R\right)-V_0\left(R_{\odot }\right)$ in sectors A+C (panel (a)) and $V=V_0\left(R\right)$ in sectors B and D (panels (b) and (c)). For panel (a) we used cut $\eta \le 0.1$. Units: R[kpc], v[km/s].

Figure 8. Velocity curves $V=V_0\left(R\right)-V_0\left(R_{\odot }\right)$ in sectors A+C (panel (a)) and $V=V_0\left(R\right)$ in sectors B and D (panels (b) and (c)). For panel (a) we used cut $\eta \le 0.1$. Units: R[kpc], v[km/s].

Figure 9. Rotation of MW as seen in sectors B and D. Here the symbols ${\mathbf{V}}_{\mathbf{R}}$ and ${\mathbf{V}}_{\mathbf{S}}$ stand for ${\mathbf{V}}_G\left(\mathbf{R}\right)$ and ${\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right).$

Figure 9. Rotation of MW as seen in sectors B and D. Here the symbols ${\mathbf{V}}_{\mathbf{R}}$ and ${\mathbf{V}}_{\mathbf{S}}$ stand for ${\mathbf{V}}_G\left(\mathbf{R}\right)$ and ${\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right).$

Another representation of the RC is obtained from panels B and D. For $\left|Z\right|\approx 0$ and $r>0,$ the term $\mathbf{w}={\mathbf{V}}_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)-{\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)$ in Eq.(38) and its transverse projection $〈w_l〉$ are calculated as suggested in Figure 9 from two similar orthogonal triangles with angle $\alpha .$ We obtain

$$〈w_l〉={\displaystyle \frac{r}{R}}{V}_0\left(R\right);R=\sqrt{r^2+R_{\odot }^2}.$$

(42)

Since $w_l=-v_l$, we get

$$\begin{array}{cc}\hfill 〈v_l〉& =-{\displaystyle \frac{r}{R}}{V}_0\left(R\right),\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill V_0\left(R\right)& =-{\displaystyle \frac{R}{r}}〈v_l〉.\hfill \end{array}$$

(44)

The corresponding RCs are shown in panels b,c of Figure 8. These curves are a different representation of the $〈v_l〉$ curves in panels B and D of Figure 7.

Relation (43) holds only for sectors B and D, where $l\approx \pm \pi /2$ and $Z\approx 0$ (or $b\approx 0$). In the general case, the geometry is more complicated. With the use of (24) and (38) we have

$$〈v_l〉=\left({\mathbf{V}}_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)-{\mathbf{V}}_G\left({\mathbf{R}}_{\odot }\right)\right).{\mathbf{n}}_l.$$

(45)

This relation can be modified as

$$〈v_l〉=-\left(V_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right){\mathbf{N}}_{\Phi }\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)-V_0{\mathbf{N}}_{\Phi }\left({\mathbf{R}}_{\odot }\right)\right).{\mathbf{n}}_l.$$

(46)

With the use of relations (6),(3) and (1, vector ${\mathbf{N}}_{\Phi }$ can be expressed in galactic coordinates $r,l,b$ and after a few further modifications, we get

$$〈v_l〉=\gamma V_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)-V_{0}cosl$$

(47)

and if $\gamma \ne 0$, then

$$V_G\left({\mathbf{R}}_{\odot }+\mathbf{r}\right)={\displaystyle \frac{〈v_l〉+V_{0}cosl}{\gamma }};\gamma ={\displaystyle \frac{R_{\odot }cosl-rcosb}{\sqrt{R_{\odot }^2+r^2{cos}^{2}b-2R_{\odot }rcoslcosb}}}.$$

(48)

One can check that in sectors B and D, where $b\approx 0$, $l\approx \pm \pi /2$ and where we have assumed $V_G\approx V_0\left(R\right)$, relation (47) reduces to (43). Relation (48) allows us to analyze $V_G\left(\mathbf{R}\right)$ not only in narrow sectors B and D but also in the wider regions, which can provide higher statistics with smaller errors. This relation is not suitable for the reconstruction of $V_G$ in the vicinity of singularity $\gamma \approx 0$ $\left(R_{\odot }cosl\approx rcosb\right)$. In Figure 10

(blue points) we show RCs obtained with the use of (48) in the sectors Q₂ and Q₄. In the analyzed area we observe irregular fluctuations in the rotation velocity: $\Delta V_G/V_G\approx 5\%$. The relation allows us to calculate rotation velocity not only in the galactic plane, but also outside the plane. In Figure 11

we show the velocity curves $V_G\left(\left|Z\right|\right)$ calculated in sectors ${\mathrm{Q}}_{1S}-$${\mathrm{Q}}_{4S}$ and ${\mathrm{Q}}_{1N}-$${\mathrm{Q}}_{4N}$.

3.3. Six Parameters of the MW Collective Rotation

Panels E and F at the top of Figure 7 provide further important information. We observe $〈v_l〉\approx 0$ and $〈v_b〉\approx 0$, as expected in both narrow cones pointing perpendicularly from the galactic plane, where positive and negative $v_l,v_b$ are equally abundant. On the other hand, the value $〈v_{gal}〉$ increases with distance from the plane. This increase occurs in the galactic reference frame reflecting the slowing of collective rotation in the Galactocentric frame, which increases with $r\approx \left|Z\right|$. The same effect is seen even more clearly in Figure 11. Important information is obtained from Figure 12,

where dependencies of standard deviations are shown. The increasing standard deviations in panels E and F suggest a less collective, but more disorderly motion of high velocities away from the galactic plane.

In the distribution (16), we assume in a first approximation:

$${\sigma }_{\alpha }={\sigma }_{\alpha 0}+{\sigma }_{\alpha 1}\left|Z\right|;\alpha =R,\Theta ,\Phi ,$$

(49)

where ${\sigma }_{\alpha 0}$ and ${\sigma }_{\alpha 1}$ are the free parameters. We proceed as follows in their setting:

i) From the data panels A-D in Figure 12 where $Z\approx 0$, we estimate

$$\begin{array}{cc}\hfill {\sigma }_{\Theta 0}& \approx 〈{\sigma }_b〉\mathrm{in}\sec \mathrm{tors}\mathrm{A}-\mathrm{D},\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\sigma }_{R0}& \approx 〈{\sigma }_l〉\mathrm{in}\sec \mathrm{tors}\mathrm{B},\mathrm{D}\mathrm{and}\mathrm{at}\mathrm{small}r/R_{\odot },\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill {\sigma }_{\Phi 0}& \approx 〈{\sigma }_l〉\mathrm{in}\sec \mathrm{tors}\mathrm{A},\mathrm{C}.\hfill \end{array}$$

(52)

In directions other than A-D, the relations between $\sigma$’s in the galactic and Galactocentric frames are more complex. The data panels E and F (where $r\approx \left|Z\right|$), show in the region of the peaks ($r\lesssim 2.5$ kpc, see Figure 3, panels E,F) a linear increase of the corresponding mixture of $\sigma$’s with $\left|Z\right|$.

ii) To obtain the parameters ${\sigma }_{\alpha 1}$, we analyzed the respective distributions in all the sectors listed in Table 1. We found that the optimal shape of the distribution $V_{\Phi }$ in (16) is asymmetric, having different ${\sigma }_{\Phi }^{\pm }$ for the two opposite orientations, where $+(-)$ means in (against) the rotation direction:

$${\sigma }_{\Phi }^{+}={\sigma }_{\Phi 0},{\sigma }_{\Phi }^{-}={\sigma }_{\Phi 0}+{\sigma }_{\Phi 1}^{-}\left|Z\right|.$$

(53)

Thus, only ${\sigma }_{\Phi }^{-}$ depends on $\left|Z\right|$. This asymmetry also reflects the effective deceleration $〈\Delta V_{\Phi }〉$ of the collective rotation for larger $\left|Z\right|$ as mentioned above. We have

$$\begin{array}{c}\hfill 〈\Delta V_{\Phi }〉={\displaystyle \frac{1}{N}}\left({\int }_{-\infty }^0{V}_{\Phi }exp\left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{V_{\Phi }}{{\sigma }_{\Phi }^{-}}}\right)}^2\right)d{V}_{\Phi }+{\int }_0^{\infty }{V}_{\Phi }exp\left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{V_{\Phi }}{{\sigma }_{\Phi }^{+}}}\right)}^2\right)d{V}_{\Phi }\right);\\ \hfill N={\int }_{-\infty }^{0}exp\left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{V_{\Phi }}{{\sigma }_{\Phi }^{-}}}\right)}^2\right)d{V}_{\Phi }+{\int }_0^{\infty }exp\left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{V_{\Phi }}{{\sigma }_{\Phi }^{+}}}\right)}^2\right)d{V}_{\Phi }.\end{array}$$

(54)

After integration we get

$$〈\Delta V_{\Phi }〉=\sqrt{{\displaystyle \frac{2}{\pi }}}\left({\sigma }_{\Phi }^{+}-{\sigma }_{\Phi }^{-}\right)=-\sqrt{{\displaystyle \frac{2}{\pi }}}{\sigma }_{\Phi 1}^{-}\left|Z\right|.$$

(55)

The parameter ${\sigma }_{\Phi 1}^{-}$ is obtained by the fit from Figure 11, which suggests dependence

$$\begin{array}{c}\hfill V_G\left(\left|Z\right|\right)=V_0+〈\Delta V_{\Phi }〉\approx V_0-\varkappa \left|Z\right|;\varkappa =\sqrt{{\displaystyle \frac{2}{\pi }}}{\sigma }_{\Phi 1}^{-}\doteq 33.5\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{kpc}}^{-1}.\end{array}$$

(56)

For the remaining two parameters, the analysis showed that ${\sigma }_{\Theta 1}$ $\approx {\sigma }_{R1}$ is a good approximation. We denote ${\sigma }_1$ $\equiv {\sigma }_{R1}={\sigma }_{\Theta 1}$, this last free parameter was set up from the tuning of the slope in simulation panels E and F in Figure 12. A list of the six resulting parameters controlling simulation (16) is given in Table 4.

Table 4. Monte-Carlo simulation model parameters and corresponding parameters from other analyses. Velocity $V_0$ is taken from Table 3.

${\sigma }_{\Theta 0}\left[km{s}^{-1}\right]$	${\sigma }_{\Phi 0}\left[km{s}^{-1}\right]$	${\sigma }_{R0}\left[km{s}^{-1}\right]$	${\sigma }_1\left[km{s}^{-1}kp{c}^{-1}\right]$	${\sigma }_{\Phi 1}^{-}\left[km{s}^{-1}kp{c}^{-1}\right]$	$V_0\left[km{s}^{-1}\right]$	ref.
17	24	33	20	42	228	this work
$18.03\pm 0.03$	$24.35\pm 0.04$	$36.81\pm 0.07$	x	x	x	[2]
11	20	31	x	x	x	[45]

Table 4. Monte-Carlo simulation model parameters and corresponding parameters from other analyses. Velocity $V_0$ is taken from Table 3.

${\sigma }_{\Theta 0}\left[km{s}^{-1}\right]$	${\sigma }_{\Phi 0}\left[km{s}^{-1}\right]$	${\sigma }_{R0}\left[km{s}^{-1}\right]$	${\sigma }_1\left[km{s}^{-1}kp{c}^{-1}\right]$	${\sigma }_{\Phi 1}^{-}\left[km{s}^{-1}kp{c}^{-1}\right]$	$V_0\left[km{s}^{-1}\right]$	ref.
17	24	33	20	42	228	this work
$18.03\pm 0.03$	$24.35\pm 0.04$	$36.81\pm 0.07$	x	x	x	[2]
11	20	31	x	x	x	[45]

The relation (56) shows that the MW rotation at $Z_{\odot }\approx 0.0208$ kpc is $0.7$ km/s lower than at $Z=0$, which is significantly less than the total error in determining $V_0$. Therefore, we have neglected $Z_{\odot }$ in our calculation. The Wolfram Mathematica code of the generator is available on the website https://www.fzu.cz/~piska/Catalogue/genkinFEB24.nb. Figure 13

shows the examples of the distribution of simulated velocities $V_{\Phi },V_{\Theta }$ and $V_R$. The simulated distributions $V_{\Theta }$ and $V_R$ are symmetric for any $\left|Z\right|$, distribution $V_{\Phi }$ is asymmetric for $\left|Z\right|>0$.

3.4. Comparison of Simulation Model with Data

The comparison is done as follows:

1) Position $\left(r,l,b\right)$ of the source in the galactic reference frame is with the use of Eqs.3 transformed to the Galactocentric spherical frame $\left(R,\Phi ,\Theta \right)$. For this position, the velocity ${\mathbf{V}}^{gen}\left(\mathbf{R}\right)$ defined by (14) is generated according to distribution (16) with the parameters from Table 4.

2) Using (41) and (7) we get the corresponding LSR transverse velocities:

$$v_{\alpha }^{gen}\left(\mathbf{r}\right)=\left({\mathbf{V}}^{gen}\left(\mathbf{R}\right)-{\mathbf{V}}_{LSR}\right).{\mathbf{n}}_{\alpha };\alpha =l,b.$$

(57)

So for any source defined by input $\left(r,l,b,v_l,v_b\right)$ we generate a vector $\left(r,l,b,v_l^{gen},v_b^{gen}\right)$ and then create the desired distributions from both.

3) These distributions obtained from the input data and simulations will now be compared. The comparison in Figure 7 shows perfect agreement between data and simulation for panels B,D,E,F. In panels A and C, the simulation model does not reproduce local kinematic substructures generating deviations $\Delta v\lesssim 10$ km/s. The presence of these substructures shows the precision with which we work. At the same time, the fluctuations are not noticeable in the other panels because their velocity scale is coarser.

The simulation in panel C indicates that the velocity ${〈v_l\left(r\right)〉}_S$ increases with r, despite the constant parameter $V_0$. This small effect is because we are working inside the angle $b=\pm 5$ deg, which means a slight linear increase in average $\left|Z\right|$ and correspondingly some deceleration with r. So, a corresponding correction would therefore be necessary to evaluate the RC in this sector more accurately. We have checked that for smaller angles $\left|b\right|$ this effect disappears. At the same time, we do not observe a similar effect in sector A. The reason is that in a very dense field of this sector our cuts $M\le 50$ and $\eta \le 0.1$ accept only a narrow sector of the data: $\left|b\right|\lesssim 1$deg.

Figure 12 shows the corresponding velocity dispersion dependencies. In the upper panels A-D we again observe fluctuations that are not present in the simulation. In panels E,F with a coarser scale the fluctuations are not noticeable. Averaged data in panels A-D were together with panels E,F used to determine parameters ${\sigma }_{\alpha }$ in Table 4. After averaging the fluctuations the simulation model fits the data well.

Relatively small velocity fluctuations ($\Delta V_G\lesssim 10$ km/s, $\Delta V_G/V_G\lesssim 0.05$) also appear in the RC in Figure 10. The slightly decreasing simulated RC is due to the shape of the sectors ${\mathrm{Q}}_2$ and ${\mathrm{Q}}_4$, where a larger R correlates with a larger average $\left|Z\right|$, implying smaller $V_G$. Figure 11 shows that the simulation of decreasing $V_G\left(\left|Z\right|\right)$ controlled by the fitted parameter ${\sigma }_{\Phi 1}^{-}$ in Eq.(56) agrees well with the data.

The good agreement of the simulations with the data is confirmed by other results.

Figure 14. Distributions of $v_l$ and $v_b$ in sectors A-F: data and simulation model. Unit: v[km/s]. Binning: 1.6 km/s.

Figure 14. Distributions of $v_l$ and $v_b$ in sectors A-F: data and simulation model. Unit: v[km/s]. Binning: 1.6 km/s.

Figure 15. Histograms $\left(v_l,l\right)$ and $\left(v_b,l\right)$ in sectors E and F: data and simulation.Units: l[deg], v[km/s]. Binning $l,v_l,v_b$: $3.6$deg$,5$ km/s$,3$ km/s.

Figure 15. Histograms $\left(v_l,l\right)$ and $\left(v_b,l\right)$ in sectors E and F: data and simulation.Units: l[deg], v[km/s]. Binning $l,v_l,v_b$: $3.6$deg$,5$ km/s$,3$ km/s.

Figure 16. Asymmetry ${\sigma }_{\Phi }^{+}<{\sigma }_{\Phi }^{-}$ in the Galactocentric reference frame generates asymmetries in the galactic frame, see text and Figure 15.

Figure 16. Asymmetry ${\sigma }_{\Phi }^{+}<{\sigma }_{\Phi }^{-}$ in the Galactocentric reference frame generates asymmetries in the galactic frame, see text and Figure 15.

Figure 17. Dependence of mean velocity and its dispersion ($v_l$ - blue, $v_b$ - green, $v_{gal}$ - red) on distance z in sectors ${\mathrm{Q}}_{1N}$-${\mathrm{Q}}_{4N}$ and ${\mathrm{Q}}_{1S}$-${\mathrm{Q}}_{4S}$: data and simulation model. Units: z[kpc], v[km/s].

Figure 17. Dependence of mean velocity and its dispersion ($v_l$ - blue, $v_b$ - green, $v_{gal}$ - red) on distance z in sectors ${\mathrm{Q}}_{1N}$-${\mathrm{Q}}_{4N}$ and ${\mathrm{Q}}_{1S}$-${\mathrm{Q}}_{4S}$: data and simulation model. Units: z[kpc], v[km/s].

Figure 18. Distributions of $v_l$ and $v_b$ in sectors ${\mathrm{Q}}_{1S}-$${\mathrm{Q}}_{4S}$, ${\mathrm{Q}}_{1N}-$${\mathrm{Q}}_{4N}$: data and simulation model. Unit: v[km/s]. Binning: 1.6 km/s.

Figure 18. Distributions of $v_l$ and $v_b$ in sectors ${\mathrm{Q}}_{1S}-$${\mathrm{Q}}_{4S}$, ${\mathrm{Q}}_{1N}-$${\mathrm{Q}}_{4N}$: data and simulation model. Unit: v[km/s]. Binning: 1.6 km/s.

Table 5. Correspondence of dispersions ${\sigma }^{\pm }$ of velocity distributions of $v_l$ and $v_b$ at $l\approx 0^{\circ },{90}^{\circ },{180}^{\circ },{270}^{\circ }$ in sectors E,F with velocity dispersions in the Galactocentric reference frame ${\sigma }_{\Phi }^{\pm }$ and ${\sigma }_R.$ The origin of the asymmetries is suggested in Figure 16.

sector∖l[deg]	0	90	180	270
E,F	${\sigma }_l^{-}\approx {\sigma }_{\Phi }^{-}{\sigma }_l^{+}\approx {\sigma }_{\Phi }^{+}$	${\sigma }_l^{-}={\sigma }_l^{+}\approx {\sigma }_R$	${\sigma }_l^{-}\approx {\sigma }_{\Phi }^{+}{\sigma }_l^{+}\approx {\sigma }_{\Phi }^{-}$	${\sigma }_l^{-}={\sigma }_l^{+}\approx {\sigma }_R$
E	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{-}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{+}$	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{+}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{-}$
F	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{+}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{-}$	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{-}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{+}$

Table 5. Correspondence of dispersions ${\sigma }^{\pm }$ of velocity distributions of $v_l$ and $v_b$ at $l\approx 0^{\circ },{90}^{\circ },{180}^{\circ },{270}^{\circ }$ in sectors E,F with velocity dispersions in the Galactocentric reference frame ${\sigma }_{\Phi }^{\pm }$ and ${\sigma }_R.$ The origin of the asymmetries is suggested in Figure 16.

sector∖l[deg]	0	90	180	270
E,F	${\sigma }_l^{-}\approx {\sigma }_{\Phi }^{-}{\sigma }_l^{+}\approx {\sigma }_{\Phi }^{+}$	${\sigma }_l^{-}={\sigma }_l^{+}\approx {\sigma }_R$	${\sigma }_l^{-}\approx {\sigma }_{\Phi }^{+}{\sigma }_l^{+}\approx {\sigma }_{\Phi }^{-}$	${\sigma }_l^{-}={\sigma }_l^{+}\approx {\sigma }_R$
E	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{-}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{+}$	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{+}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{-}$
F	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{+}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{-}$	${\sigma }_b^{-}={\sigma }_b^{+}\approx {\sigma }_R$	${\sigma }_b^{-}\approx {\sigma }_{\Phi }^{-}{\sigma }_b^{+}\approx {\sigma }_{\Phi }^{+}$

Figure 14 shows distributions of $v_l$ and $v_b$ in sectors A-F along with the corresponding distributions obtained from simulations. In sectors A-D we observe a narrower distribution of $v_b$, the width of which does not depend on the sector. The distributions have roughly the same dispersion for both data and simulation, but the peak in the data is sharper. The peaks in the data also have slightly longer tails, which increase the dispersion. The simulated shape is therefore only approximate. The distributions of $v_l$ are slightly different. Note the shift in sectors B and D resulting from the decrease in $〈v_l〉$ in upper panels B and D in Figure 7. The distribution $v_l-$ DATA in sector A (where $R<R_{\odot }$) also has a sharper peak with an apparent asymmetry. This may be a manifestation of the asymmetric drift effect. So the simulated shape is only approximate here as well. In the opposite sector C (where $R>R_{\odot }$), the asymmetry is negligible. The distributions of $v_l$ in both sectors A and C are copies (up to a constant shift) of the distribution of orbital velocities in the Galactocentric frame.

The important result is shown in Figure 15. The asymmetry of histograms $\left(v_l,l\right)$ and $\left(v_b,l\right)$ in sectors E and F with the empty spaces reflects different projections of the asymmetry (53). The shape of histograms can be explained using Figure 16 and Table 5. Distributions of $v_l$ for $l\approx 0^{\circ },{90}^{\circ },{180}^{\circ },{270}^{\circ }$ are controled by the parameters ${\sigma }_l^{\pm }$ in the first row of table. Their connection with ${\sigma }_{\Phi }^{\pm }$ and ${\sigma }_R$ can be deduced from figure. At $l\approx 0^{\circ },$ the directions of $v_l$ and MW rotation are identical, so ${\sigma }_l^{\pm }\approx {\sigma }_{\Phi }^{\pm }$. But at $l={180}^{\circ }$, the two directions are opposite, so ${\sigma }_l^{\pm }\approx {\sigma }_{\Phi }^{\mp }$. At $l\approx {90}^{\circ }\left({270}^{\circ }\right)$ the situation is a little more complicated. If $\left|Z\right|/R$ or $\Theta$ are small (which is almost our case, see panels E,F in Figure 3, where r$\approx \left|Z\right|$), then the $v_l$ direction can be approximated by vector $\mathbf{R}$, so ${\sigma }_l^{\pm }\approx {\sigma }_R$. Similarly for $v_b$ distributions controlled by ${\sigma }_b^{\pm }$ in the next rows of the table. Histograms involve integrated distributions over sectors E and F. Also here, the agreement between the data and simulation model is perfect.

Correct Monte-Carlo parameter settings can be verified in wide Q-sectors outside the area of the galactic plane. Figure 17 shows $z-$dependence of mean values and dispersions of velocity distributions in these sectors. The curves together with the corresponding overall distributions of velocities in broad sectors ${\mathrm{Q}}_{1N}-$${\mathrm{Q}}_{4N}$, ${\mathrm{Q}}_{1S}-$${\mathrm{Q}}_{4S}$ in Figure 18 again confirm the perfect agreement of the simulation with data. Note in particular the projections $v_l$ in sectors ${\mathrm{Q}}_{1N}$, ${\mathrm{Q}}_{3N}$, ${\mathrm{Q}}_{1S}$, ${\mathrm{Q}}_{3S}$ and $v_b$ in sectors ${\mathrm{Q}}_{2N}$, ${\mathrm{Q}}_{4N}$, ${\mathrm{Q}}_{2S}$, ${\mathrm{Q}}_{4S}$. This is also due to the asymmetry (53) that occurs for $\left|Z\right|>0$, as illustrated by simulated $V_{\Phi }$ distribution in Figure 13.

4. Discussion and Conclusion

1) The results on local Sun’s velocity ${\mathbf{v}}_{\odot }=\left(U_{\odot },V_{\odot },W_{\odot }\right)$, the MW circular velocity $V_c$ and average orbital velocity $V_0$ can be compared with the other measurements presented in Table 3. The excellent agreement with others is obtained for the rotational velocity of the Sun $V_c+V_{\odot }=V_0+\Delta V_{y\odot }$. Within the measurement errors, this result perfectly agrees with all the others, including a very accurate measurement [26]. There is also perfect agreement with the others for $W_{\odot }$ and $V_0$ and a good agreement for $U_{\odot }$. However, differences in $V_c$ and $V_{\odot }$ from different measurements are larger, so our values agree with only some of them (within errors).

Let’s add a few remarks on our measurement of $\Delta {\mathbf{V}}_{\odot }$ and $V_0$. Both velocities are obtained independently with the use of a direct and model-independent method. The $V_{\odot }$ is measured equally as $W_{\odot }$, see Figure 5. For $r_{max}>0.5$ kpc we obtain a (nearly) constants in the band of small statistical errors. The correct determination can be verified as follows. At top half of Figure 19 we show the dependencies of $〈v\left(z\right)〉$ (magnified parts of Figure 17) in the eight ${\mathrm{Q}}_{S,N}$ sectors in the LSR reference frame. This frame defines the $\Delta {\mathbf{V}}_{\odot }$ determined from the data in sectors A-D, see Table 3. In these figures, we observe that for $z\to 0$ we have $〈v_l\left(z\right)〉,〈v_b\left(z\right)〉$$\to 0$, which means that the solar velocity related to the sectors A-D is the same as the velocity related to the nearby stars in all sectors ${\mathrm{Q}}_{S,N}$. This agreement is a simple test that our velocity $\Delta {\mathbf{V}}_{\odot }$ (Eq.(36)) is correct. In the bottom half of the figure, we have the same curves for comparison but in a galactic reference frame (Sun’s rest frame). Note the different scales on the upper and lower parts.

Figure 19. Dependence of mean velocity ($v_l$ - blue, $v_b$ - green) on distance z from the galactic plane in sectors ${\mathrm{Q}}_{1N}$-${\mathrm{Q}}_{4N}$ and ${\mathrm{Q}}_{1S}$-${\mathrm{Q}}_{4S}$ in the LSR reference frame (upper part) and the Sun’s rest frame (lower part). Units: z[kpc], v[km/s].

Figure 19. Dependence of mean velocity ($v_l$ - blue, $v_b$ - green) on distance z from the galactic plane in sectors ${\mathrm{Q}}_{1N}$-${\mathrm{Q}}_{4N}$ and ${\mathrm{Q}}_{1S}$-${\mathrm{Q}}_{4S}$ in the LSR reference frame (upper part) and the Sun’s rest frame (lower part). Units: z[kpc], v[km/s].

2) The determination of RC is based on the model-independent definition (15). In Figure 8 and Figure 10 we show RCs measured in different sectors of galactic longitudes. The curves are obtained with very high precision, so as a result, we observe local fluctuations ($\lesssim 10$ km/s) in the structure of MW rotation. These fluctuations correspond to velocity substructures and non-axisymmetric kinematic signatures mentioned in Section 1. The fluctuations do not allow us to analyze the RC slope in our limited range of R . Further, we have shown that the collective rotation velocity decreases for increasing $\left|Z\right|$, see Figure 11. A similar observation was reported in [2,47]. The slope of the curves is defined in Eq. (56) and can be compared with a prediction [9](Figure 5) of the Besançon model. Within the range of our analysis ($0-3$ kpc) the average slopes agree very well, although the shape of the curves is slightly different.

3) Except for observed local fluctuations, the analyzed kinematical distributions are very well described by a minimal MW axisymmetric model based on six free parameters in the Galactocentric reference frame. The model describes a simplified scenario in which local velocity fluctuations are averaged. The scale of averaged fluctuations increases with distance from the galactic plane and is defined by parameters of the model in Table 4. The fluctuations are most significant in the ${\mathbf{N}}_R$ direction, less in the ${\mathbf{N}}_{\Phi }$ direction and least in the ${\mathbf{N}}_{\Theta }$ direction.

Analysis and simulation of kinematical distributions in the studied region need apart of the MW parameters also other parameters related to our laboratory: its local 3D velocity $\Delta {\mathbf{V}}_{\odot }$, distance from the galactic centre $R_{\odot }$ and its position $z_{\odot }$ above the galactic plane (neglected). Thus, except for $R_{\odot }$, all the remaining parameters that we obtained in the present analysis are listed in Table 3 and Table 4. For now, we ignore the slope of the RC (Eq.(31)), which has in our region a very small effect [18,39]. The model describes the MW rotation as follows.

a) The rotation is strongly collective in the galactic disk plane with relatively small Gaussian velocity fluctuations ${\sigma }_{\Phi 0},{\sigma }_{\Theta 0},{\sigma }_{R0}$ around the much greater velocity $V_0$. This can be seen in Figure 14 in panels $v_l$ for sectors A and C, and panels $v_b$ for sectors A-D. Wider and shifted distributions $v_l$ for sectors B and D are due to the rotation effect shown in Figure 9 and expressed in Eq.(43). Our first three parameters are compared with corresponding galactic thin disc parameters obtained in another study, see Table 4. The agreement in ${\sigma }_{R0}$ and ${\sigma }_{\Phi 0}$ is very good. Some disagreement in ${\sigma }_{\Theta 0}$ obtained in [45] maybe related to the fact that the Gaussian distribution shape may not fit very well here (see comment on Figure 14 in Section 3.4).

b) The parameters ${\sigma }_1$ and ${\sigma }_{\Phi 1}^{-}$ are important outside the galactic plane, where they control the increase in fluctuations with $\left|Z\right|$ as shown in sectors E,F in Figure 12 and the slope of curves in Figure 11. These figures suggest that the collective velocity decreases with increasing $\left|Z\right|$ and the directions of the trajectories are becoming more random and probably less circular. The further effect of ${\sigma }_{\Phi 1}^{-}$ is due to the asymmetry of distribution $V_{\Phi }$, which generates deceleration of collective rotation with increasing $\left|Z\right|$ according to Eq.(55). This asymmetry is also manifested very clearly in Figure 15 and Figure 18 (see the comments in the last two paragraphs of Section 3.4).

c) Our assumption that ${\sigma }_{\Phi }^{+}$ does not depend on $\left|Z\right|$ can be verified by comparing $v_l$-distributions in the sectors ${\mathrm{Q}}_{1N}$, ${\mathrm{Q}}_{3N}$, ${\mathrm{Q}}_{1S}$, ${\mathrm{Q}}_{3S}$ in Figure 18 with the corresponding distributions in sectors A and C in Figure 14. This independence means that in the analysed region

$$\left|\mathbf{V}\left(R,Z\right)\right|\lesssim V_0+{\sigma }_{\Phi 0}\approx 252\mathrm{km}/\mathrm{s};5\lesssim R\lesssim 13\mathrm{kpc},\left|Z\right|\lesssim 3\mathrm{kpc}.$$

(58)

We have shown that the 3D Monte-Carlo model fits all studied sectors of the averaged kinematic data very well. Of course, its parameters may require further optimization in more distant regions.

To conclude, the proposed statistical methods for calculating the local velocity of the Sun, the average rotation velocity $V_0$ and generally the velocity $V_G(R,Z)$ at different positions in the MW can be useful for the analysis of the current and future Gaia data releases. It is always important to be able to compare these parameters obtained by different methods and input data samples. Averaged, axisymmetric approximation of the MW kinematics represented by Monte-Carlo simulation code can be useful in validating axisymmetric dynamic models or determining the scale of local kinematical substructures out of axial symmetry.