A new kinematic model of the Galaxy: analysis of the stellar velocity field from Gaia Data Release 3

This work presents the results of a kinematic analysis of the Galaxy that uses a new model as applied to the newest available Gaia data. We carry out the Taylor decomposition of the velocity field up to second order for 18 million high luminosity stars (i.e. OBAF-type stars, giants and subgiants) from the Gaia DR3 data. We determine the components of mean stellar velocities and their first and second partial derivatives (relative to cylindrical coordinates) for more than 28 thousand points in the plane of our Galaxy. We estimate Oort's constants $A$, $B$, $C$, and $K$ and other kinematics parameters and map them as a function of Galactocentric coordinates. The values found confirm the results of our previous works and are in excellent agreement with those obtained by other authors in the Solar neighbourhood. In addition, the introduction of second order partial derivatives of the stellar velocity field allows us to determine the values of the vertical gradient of the Galaxy azimuthal, radial and vertical velocities. Also, we determine the mean of the Galaxy rotation curve for Galactocentric distances from 4 kpc to 18 kpc by averaging Galactic azimuths in the range -30$^\circ$<$\theta$<+30$^\circ$ about the direction Galactic Centre-Sun-Galactic anticentre. Maps of the velocity components and of their partial derivatives with respect to coordinates within 10 kpc of the Sun reveal complex substructures, which provide clear evidence of non-axisymmetric features of the Galaxy. Finally, we show evidence of differences in the Northern and Southern hemispheres stellar velocity fields.


INTRODUCTION
The Gaia DR3 catalogue (Gaia Collaboration et al. (2016a, 2020, 2022)) has provided new data for studying the kinematics of stars in the Milky Way.The presence of the spatial positions and velocities of 33 million stars in this release makes it possible to obtain new information about stellar kinematics in our Galaxy.Thus, the large amount of data and the new level of its precision require us to develop new approaches to its analysis in the context of kinematic studies of the Galaxy.
Recently, some new methods are often used to analyse stellar kinematics (Poggio et al. 2018;Eilers et al. 2019;Antoja et al. 2018;Gaia Collaboration et al. 2022;Nelson and Widrow 2022).Nelson and Widrow (2022) constructed a non-parametric, smooth and differentiable model of the underlying velocity field using a sparse Gaussian process algorithm based on induction points.They ★ Contact e-mail: akhmetovvs@gmail.com(VSA)estimated the Oort constants , ,  and  and presented the maps of the velocity field and the divergence within 2 kpc of the Sun.Gaia Collaboration et al. (2022) selected various stellar populations to explore and identify non-axisymmetric features in the Milky disc and mapped the spiral structures associated with star formation in the scale range 4-5 kpc from the Sun.They used the sample of the red giant branch to obtain the velocity field maps of the Milky Way as far as 8 kpc from the Sun.Poggio et al. (2018) showed that the large-scale kinematics has a clear signature of the warp of the Milky Way appearing as a gradient of 5-6 km s −1 in the vertical velocities between 8 and 14 kpc of the Galactic radius.
Using the Ogorodnikov-Milne (O-M) model Fedorov et al. (2021) obtained the kinematic parameters of red giants and subgiants from the Gaia EDR3 catalogue along the direction of the Galactic center -the Sun -the Galactic anticentre.Fedorov et al. (2023) mapped the kinematic parameters in the plane of our Galaxy using the following approach.They solved the equations for the scalar stellar velocity field to find the velocity components of the centroid and their first order partial derivatives with respect to the Galactocentric cylindrical coordinate system.Dmytrenko et al. (2023) analysed the O-M deformation velocity tensor to determine the vertices of selected stellar samples.
In addition to classical physical models (e.g., the Oort-Lindblad, Ogorodnikov-Milne), other mathematical models are often used to analyze the stellar kinematics of the Galaxy.Assuming that stellar proper motions and radial velocities are the components of the velocity vector field, it is reasonable to use the methods of decomposing the corresponding stellar velocity field into a set of vector spherical harmonics (VSH).This approach allows the detection of all systematic components contained in the studied velocity field.The comparison of the decomposition coefficients with the model parameters shows whether the models are complete, and allows to reveal any significant systematics that are not included in the models.This approach is very effective, as previous studies have shown (Vityazev & Tsvetkov 2005;Makarov & Murphy 2007;Mignard & Klioner 2012;Velichko et al. 2020).Nevertheless, the main problem of this approach is to find the physical sense of the decomposition coefficients and to relate them with the kinematics features of our Galaxy.
In this paper, we obtain a detailed kinematics in the Galactic plane based on the Gaia DR3 data, which are provided with spatial positions and velocities.We derive the kinematic parameters by analysing with a new method the high luminosity stars, whose velocities are considered relative to the Galactic center.We exclude the velocity of the Sun  ⊙ relative to the Galactic center from the relative velocities of all stars.Consequently, this allowed us to analyze the field of stellar velocities in a Galactocentric cylindrical coordinate system.

DATA
We select Gaia DR3 sources for which the six dimensional phase-space coordinates can be computed, i.e. all stars with fiveparameters astrometric solution (sky positions, parallax and proper motions) and line-of-sight (radial) velocity (about 33 million stars in Gaia DR3).We retain only objects that have i) Renormalized Unit Weight Error ruwe < 1.4 in order to discard the sources with problematic astrometric solutions, astrometric binaries, and other anomalous cases (Lindegren 2018), ii) the relative parallax error /  > 5 (i.e. the inverse-parallax distances better than 20% ), iii) and three Gaia photometric bands ,   and   .We estimate the distance to the stars as 1/, where  is the parallax.As it was noticed in Antoja et al. (2018), the differences between distances, determined with help of a Bayesian method and of the inverse of the parallax, are between −2% and 0.6% for 90% stars with /  > 5.Such a result is expected for small relative errors in parallax.
We also exclude stars that are at large distances from the mid-plane of the Galaxy || > 1 kpc and also unbound stars or high-velocity stars with Galactocentric velocity   = √︃  2  +  2  +  2  > 500 km s −1 (Marchetti et al. 2022).Furthermore, for our kinematic analysis we select only high luminosity stars with   < 4 and exclude faint stars that are dominant near the Sun and are incomplete at large distances.The spatial distribution of our stellar tracers is more uniform and covers a wide range of Galactocentric distance range, i.e., 0 kpc <  < 18 kpc.

COORDINATE SYSTEMS AND MEAN VELOCITY DISTRIBUTION
It is convenient to use the Cartesian coordinate system (, , ) for the basic kinematic equations of various models.The origin of the system coincides with the Solar System barycenter.The  axis points to the Galactic center, the  axis coincides with the direction of the Galactic rotation, while  is perpendicular to the Galactic plane and complements the right-handed Cartesian coordinate system (see Fig. 2).The given system is called the local (heliocentric) rectangular Galactic coordinate system.The transformation of the positions and velocity components from the spherical coordinate system (, , ,   ,   ,   ) to the rectangular Cartesian one (, , ,   ,   ,   ) is where   is the radial velocity,   =     cos ,   =     are the components of stellar proper motion, and  = 4.74057 is the transformation factor from mas yr −1 to km s −1 kpc −1 .The distance  is calculated from the Gaia DR3 parallax as 1/.Formulas (2) allow us to obtain the velocity field of the stars relative to the Sun in the rectangular coordinate system.We calculated the kinematic parameters using both uncorrected Gaia DR3 data and also considering the zero point correction for parallax and proper motion (Gaia Collaboration et al. 2021;Cantat-Gaudin and Brandt 2021).It turns out that for points with a heliocentric distance less than 5 kpc, the behaviour of the kinematic parameters does not change practically, except for a small "compression" in the direction of the Sun, which does not affect the behaviour of the parameters in any way.Some differences are possible at large distances from the Sun, more than 5-6 kpc, where the errors in estimating the astrometric parameters are large.In this work we have chosen to use uncorrected Gaia DR3 data, since the main goal of the work is to demonstrate the capabilities of the presented method on Gaia DR3 data, which will be corrected and improved several times in systematic and random accuracy in the next data releases.
We calculate the velocity components of each star that determine the velocity field relative to the motionless Galactocentric coordinate system as: where   ,   ,   are the velocity components of each star relative to the Solar system barycenter and the values of the solar velocity components (  ,   ,   ) ⊙ = (9.5, 250.7, 8.56) km s −1 are relative to the Galactic center (Reid & Brunthaler. 2020;GRAVITY Collaboration et al. 2021;Drimmel&Poggio 2018).
We estimate the mean 3D velocity distribution by sampling our stellar catalogue with cubic cells of 100pc×100pc×100pc.In each cell, the mean positions (, , ) and spatial velocities (  ,   ,   ) are derived by averaging the positions and velocities of the stars.These mean values are also called "fictitious stars" in the following sections.The total number of cells (fictitious stars) is about 270 thousand.To perform the averaging, a minimum number of 3 stars is required.
Fig. 3 shows the stellar density in (x,y), (y,z) and (x,z) planes for the cell size of 100pc×100pc×100pc in the Cartesian heliocentric system (Fig. 2).Fig. 3 clearly shows the incompleteness due to the strong interstellar extinction towards low Galactic latitudes.In particular, the lack of stars close to the Galactic plane is apparent in the inner disc ( > 4 kpc) and towards  >4-5 kpc from the Sun.
The conversion from the Cartesian heliocentric system to the Galactocentric cylindrical coordinate system also requires knowledge of the distance of the Sun to the Galactic center  ⊙ = 8.249 kpc, and height of the Sun from the Galactic plane of  ⊙ = 0.0208 kpc (Bennett, M., & Bovy, J. 2019;GRAVITY Collaboration et al. 2021).We have applied the transformation procedure to the Galactic cylindrical coordinate system to each cell (fictitious star): =  , .

METHODS FOR DETERMINATION OF KINEMATIC PARAMETERS
The Ogorodnikov-Milne (O-M) model is based on Helmholtz's theorem for stellar systems: within a small region surrounding any point in a stellar system, the velocities of the centroids can be modelled as the superposition of a rigid translation and rotation velocity plus a deformation component velocity of the region (Ogorodnikov 1965).Ogorodnikov (1932) proposed this model and Milne (1935) generalised it to the three-dimensional case.The stellar velocity field, which fills a given spatial region, is treated as a continuous deformable medium, which allows us to study stellar motion using the methods of continuum mechanics.
In the Cartesian Galactic coordinate system (  ,  = 1, 2, 3), the general form of the expansion of the velocity field V(r) in the vicinity of the centroid located at r 0 is: where "0" denotes the derivatives at  0 = ( 0 ,  0 ,  0 ).In order to analyze the stellar velocities   inside the sphere with radius  ≪  0 (distance to the Galactic Centre) in the O-M model, we can restrict the expansion of the series to the first order terms: Here and 9) are components of antisymmetric  − and symmetric  + tensors respectively: and ,  = 1, 2, 3.The antisymmetric tensor  −  =   is also called the local rotation tensor, since it is equivalent to the axial vector ( 1 ,  2 ,  3 ), where  1 =  − 32 ,  2 =  − 13 ,  3 =  − 21 .The second-rank symmetric tensor V  is called the local deformation velocity tensor.Taking into account (10) we can rewrite the expression (7) in the following form: The elements of the matrices  − and  + are partial derivatives with respect to the coordinates of the projections ( 1 ,  2 ,  3 ) of the velocity vector onto the rectangular coordinate axes.These elements are referred to as the kinematic parameters.The relations between the kinematic parameters of the O-M model and the partial derivatives of the first order velocity components with respect to Galactocentric coordinates (, , ) are as follows: where   has the opposite direction with respect to the circular rotation velocity of the Galaxy  rot .We can also determine the Oort constants , ,  and  which measure the azimuthal shear, the vorticity, the radial shear, and the divergence of the velocity field in the Galactic plane respectively.The Oort constants in Galactocentric cylindrical coordinates (, ) are given by (Chandrasekhar 1942;Ogorodnikov 1965): Formulas ( 12) -( 20) do not allow us to determine the parameters and    from combinations of the O-M model parameters.It is therefore necessary to introduce some reasonable assumptions about these derivatives.
The O-M model has been used by many authors to find the kinematic parameters of the Galaxy at one point -near the Sun.Since the Gaia DR3 catalogue provides a set of spatial positions and velocities of the stars that extend well beyond the solar neighborhood, we can perform similar kinematic studies at any point of the Galaxy as shown in Fedorov et al. (2021Fedorov et al. ( , 2023)).
First of all, we create a rectangular grid specified in the Galactic plane where  = 0.The grid nodes are separated from each other along the coordinates  and  by a distance of 100 pc in the range −10 kpc and +10 kpc for  and  axes of the heliocentric coordinate system.The total number of nodes in the Galactic plane is 40 thousand.For each node of the grid, we select a spherical region with a radius of 1.0 kpc.Each sphere around a given node contains a set of  = 1, 2, 3.... cells (fictitious stars) that are located at distances not exceeding 1 kpc, i.e. each -th fictitious star satisfies the condition: The optimal size of the spheres was chosen after numerous calculations using different radii.As a result, we obtained a radius of 1 kpc, which is a compromise between a large number (accuracy) of fictitious stars and the high resolution (detail) of the parameter maps.When the radius of the spheres is small, the obtained kinematic parameters are sensitive to local variations in the stellar velocity field.Since we want to study the global behaviour of the kinematic parameters, we use a radius of 1 kpc, which is optimal for these scales.
The kinematic parameters were calculated for each node with  and  coordinates and located at a distance of 100pc×100pc in the Galactic plane.Thus, the overlap of our spherical samples provides for the continuity and differentiability of the velocity field, which makes it possible to study global kinematic parameters instead of studying them as discrete values at local points.Only spheres containing  ≥ 500 cells (fictitious stars) have been used in this investigation (Fig. 4).However, we trust the results of calculations for those points on the map for which the number of cells (fictitious stars) is more than 2000 which corresponds to a heliocentric distance of 5-6 kpc (Fig. 4).
We thus form 28139 spherical regions (centroids) with coordinates  and  in the Galactic plane ( = 0), each containing 500 <  < 3981 fictitious stars.For each spherical region we compute the cylindrical Galactocentric velocity components and their first and second order partial derivatives using the position and velocity of the fictitious stars from the solution of the equations by means of the least squares method (LSM): These formulas are obtained by expanding the cylindrical Galactocentric velocity components   (, , ),   (, , ),   (, , ) into a Taylor series limited by the second order terms of all fictitious stars with coordinates   ,   and   satisfying condition 27.The expansions were made in the vicinity of the points with coordinates (   ,    , 0).These points coincide with the nodes of the rectangular grid (, ).Condition equations ( 28), ( 29), (30) were formed for each -th fictitious star, and the system of equations was solved by least squares in 30 unknowns: three velocity components   ,   ,   , their nine partial derivatives of the first order and eighteen partial derivatives of the second order with respect to coordinates for each point with coordinate , .The system of equations was compiled and solved independently for each spherical region (centroid) with  and  coordinates.Thus, in the Galactic plane we obtained 28139 solutions of the equations system ( 28), ( 29), ( 30) and each solution containing 30 kinematic parameter estimates and their errors.The found kinematic parameters describe the stellar velocity field at any given point with  and  coordinates.This approach allows us to represent the behaviour of the calculated kinematic parameters from Galactic coordinates (, ) or (, ).
We have calculated the kinematic parameters in two ways: the first method is the use of "fictitious stars" (as presented in this paper) and the second method is the use of individual stars in equations( 28), ( 29), (30).There are practically no differences between the results of the two methods, but the differences in the speed of calculations are significant.For this reason, we show the results of the faster computational method (using "fictitious stars"), which is particularly important for future Gaia releases that will contain more data.

ANALYSIS OF THE RESULTS
In this Section we present the distribution of the kinematic parameters estimated from Eqs. 28-30 in the form of colour maps.The kinematic parameters derived by the O-M model are in good agreement over the entire coordinate range with the stellar velocity field estimated from a Taylor series limited by the first order in the cylindrical Galactocentric coordinate system, as also found by Fedorov et al. (2023).Furthermore, our approach provides additional kinematic parameters, including those that cannot be determined with the O-M model.On the colour maps, the coordinate plane (, ) coincides with the Galactic plane.The kinematic parameters in each point of this plane are calculated in the Galactocentric cylindrical coordinate system ,  and  = 0.The Sun is represented as ⊙ with coordinates (, ) = (0, 0) and ( ⊙ , ) = (8.249kpc, 180 • ) .
The mean value of the random errors of the kinematic parameters is 1.5-2 km s −1 kpc −1 up to heliocentric distances of about 8 kpc for the first order derivatives and 5-6 kpc for those of the second order; there is no correlation between the behaviour of the parameters and their errors.This gives us reason to consider these parameters as real.The actual distribution of the random errors of the kinematic parameters are shown in Appendix A (Fig. A1, A2, A3 and A4).

Oort constants 𝐴, 𝐵, 𝐶 and 𝐾
In figure 5 we show the values of the Oort constants , ,  and  calculated from the equations ( 21), ( 22), ( 23) and ( 24) in the whole Galactic plane for different Galactocentric distances and angles, and their change in different directions.Our estimations of the Oort constants are listed in Tab.A2 and are in good agreement with those made by Siebert et al. (2011), Bobylev &Bajkova (2017), Bovy (2017), Vityazev et al. (2018), Wang et al. (2021).They calculated the Oort constants only at one point using stellar samples at distances of 0.5 or 1.5 kpc near the Sun.Similarly, a reasonable agreement with the values of these parameters in the Galactic plane was found in Nelson and Widrow (2022).
It results from Figure 5 that the Oort constants  and  deviate significantly from zero, reaching the values of ±(4-5) km s −1 kpc −1 as a function of the Galactic coordinates  and .The dependence of the Oort constants  and  on the distance  and angle  indicates a non-axisymmetric Galactic potential.It is also evident from figure 5 that the Oort constants ,  and  show significant variations with .

Components of the stellar velocity field
Figure 6 shows the ordered (top) and the random (bottom) stellar velocity components   ,   ,   calculated using equations ( 28), ( 29) and (30).The stellar velocity dispersions are estimated by the standard deviation    ,    ,    resulting from the LSM solution of the spherical regions with the radius of 1.0 kpc, as described in Section 4. Our velocity distributions (  ,   ,   ) as well as their standard deviations are in good agreement with previous kinematics studies (Gaia Collaboration et al. 2022;Queiroz et al. 2021).
It is seen from the right panels of figure 6 that the values of the vertical stellar motions   and their standard deviation    are the smallest from all the components.The vertical velocity component   practically does not change in the range of Galactocentric distances from 1 to 10 kpc and has a value close to zero.The vertical velocity increases with the Galactocentric radius for distances larger than 10 kpc showing the kinematic signature of the Galactic  21), ( 22), ( 23) and ( 24) by data solution from the equation ( 28) and ( 29 28), ( 29) and ( 30).Bottom: Their standard deviation    ,    ,    obtained by solving the equation system together with LSM.al. 2020).Conversely, our velocity components   (   ,    , 0),   (   ,    , 0),   (   ,    , 0) and the corresponding velocity dispersions result from an analytical approximation of the stellar velocity field to the Galactic plane ( = 0), whose second order term, i.e.  2 / 2 , is able to model the vertical velocity gradient.This is an important improvement, since the number of stars with heliocentric distances larger than 5-6 kpc decreases drastically towards the Galactic plane due to the strong extinction, so that our estimate of the velocity field at  = 0 is dominated by stars located above the Galactic plane.

The radial gradient of the stellar velocity field
In Fig. 7 we see the maps of the radial gradients of all three velocity components   /,   / and   /.The values of these parameters are in good agreement with those obtained in (Siebert et al. 2011;Nelson and Widrow 2022) and also in works (Fedorov et al. 2021(Fedorov et al. , 2023)), where these parameters were obtained in the framework of a linear model for the velocity field decomposition.
The radial gradient of the Galactic expansion velocity   / is interpreted as the contraction (expansion) of the stellar system along the Galaxy radius.It may be the result of dynamical perturbations generated by the substructures: Galactic bar, spiral arms and warp.The value of this kinematic parameter in the Galactic plane varies within ±10 km s −1 kpc −1 .The detailed analysis of the kinematic parameter   / provides a way to determine the parameters of the spiral structure, as recently shown by Denyshchenko et al. (2024).
The kinematic parameter   / is a radial gradient of the circular velocity component, which shows the value of the slope of the Galaxy rotation curve at a given point.As can be seen from Fig. 7, the value of this parameter varies between −30 km s −1 kpc −1 near the Galaxy center to 10 km s −1 kpc −1 at the Galactocentric distances of 10 kpc and more.Note that   has the opposite sign with respect to the Galactic rotation  rot (see Eqs. ( 13) and ( 25)), hence the negative value of the derivative   / represents an increase in  rot and vice versa.Figure 7 shows a rapid increase in the rotation velocity from the Galaxy center up to  ≃ 5 kpc.A slow decrease in the Galaxy rotation velocity is seen from 7.5 kpc, corresponding to a positive value of   /.The slope   / appears in good agreement with the map of the azimuthal velocity   (fig.6) and with the rotation velocity curve (fig.10).
In the right panel of figure 7 the value of the radial gradient of the vertical velocity   / is close to zero in the Solar neighbourhood and less than 2 km s −1 kpc −1 in other regions.This means that the first order derivatives of the vertical velocity component along the Galactic radius and azimuth directions are very small values, as we can see in the right panels of Figs 7 and 8.This corresponds to a smooth and slow increase of the vertical velocity with increasing Galactocentric distance.Only at a Galactocentric distance greater than 12 kpc the value of   / is equal to +3 km s −1 kpc −1 , which is a signature of the Galactic warp.

The azimuthal gradient of the stellar velocity field
Using our method, for the first time we have obtained the estimates of the azimuth gradients of the velocity components.In the assumption of an axisymmetric potential, all velocity gradients along the changing angle  must be zero.However, as can be seen from Figure 8, the dependencies of   / on the left and   / in the middle panel are non-zero.
The kinematic parameter (1/)  / is the azimuthal gra-  dient of the radial velocity component.As can be seen from Fig. 8 (left), this gradient is a function of coordinates and it reaches values of ±(3-5) km s −1 kpc −1 .The estimation of this parameter is very important to determine the Galaxy rotation velocity, as well as to find the correct value for the Oort constants , ,  and , as we can see from the equations ( 21), ( 22), ( 23) and ( 24).
Since the parameter (1/)  / cannot be directly estimated by the O-M model, it is conveniently put to zero in the frame of the axisymmetric rotation model.In Fig. 10 we can see the difference between the Galaxy rotation curves obtained in the framework of the O-M model (blue dots -equation ( 26)), and with the model that takes into account (1/)  / (green dots) as shown in equation (25).This difference in the Galactic rotation velocity reaches ±10 km s −1 in the range of 6-12 kpc in Galactocentric distances.At other distances, if we do not consider the parameter (1/)  /, the Galaxy rotation velocity component has a significant distortion reaching several tens of km s −1 .We can conclude that not taking into account the parameters (1/)  / introduces a significant distortion in the determination of the Galaxy rotation curve.
The middle panel of Fig. 8 shows the azimuth gradient of the circular velocity component (1/)  /.As mentioned above, this parameter cannot be determined from the O-M model without additional assumptions.On the contrary, the approach we have implemented allows us to determine (1/)  / from equation (29) without additional conditions.
In figure 8, we observe a smooth change of the kinematic parameter (1/)  / as a function of the  angle.Thus, at distances greater than 8 kpc from the Galactic center and for  ≃ 210 • , the values of this parameter are about −5 km s −1 kpc −1 and increase in the direction of the Galactic rotation up to a value of +5 km s −1 kpc −1 for  ≃ 150 • .Along the line between the Sun and Galactic anticenter, the value of (1/)  / is close to zero.An asymmetric azimuthal gradient of the circular velocity with opposite slope for the line with  > 215 • and  < 215 • is well seen in the inner disc ( < ∼ 4 kpc).This is an additional kinematic signature of the Galactic bar already evidenced in Fig. 6.
We also notice that the middle panel of figure 8 shows an opposite velocity gradient (1/)  /, both a positive and a negative for  > +7 kpc and  < −7 kpc, respectively.This result derives from the higher rotation velocity,  rot ≡ −  > ∼ 250 km s −1 , estimated in the outer regions (| | > 7 kpc) with respect to the mean velocity  rot ∼ 220 km s −1 in the Solar neighbourhood (see fig. 6).The high value of the rotation velocity in the outer regions is likely overestimated because of the incompleteness of our stellar catalog close to the Galactic plane.This means that the azimuthal velocity (  )   at  = 0 may result biased because the best fit of the second order Taylor expansion (29 is mainly constrained by stars outside the Galactic plane. As described earlier, the azimuthal gradient of the vertical velocity (1/)  / is practically zero in the entire plane of the Galaxy.This parameter is associated with the kinematic evidence of the Galactic warp.Miyamoto & Zhu (1998) and Zi Zhu (2000) determined the following value (1/)  / = 3.8±1.1 km s −1 kpc −1 based on the proper motions of O-B stars in the HIP-PARCOS catalogue.At the same time, Mignard (2000) also found no confirmation of the Galaxy warp based on the HIPPARCOS data.Vityazev & Tsvetkov (2012), applying the Vector Spherical Harmonics (VSH) and Zonal Vector Spherical Harmonics (ZVSH) methods based on the stellar proper motions of the Tycho2 and UCAC3 (Zacharias et al. 2010) catalogues, obtained statistically insignificant parameter values (1/)  / both for the whole sphere and separately for the Northern and Southern hemispheres.From our analysis, the vertical velocity does not show dependence on the Galaxy coordinates  and  as we can see in Fig. s 7, 8 and 12 (right panels).

The vertical gradient of the stellar velocity field
In this section we consider the first order vertical gradients of the velocity components   /,   /, and   / illustrated in Fig. 9. Vityazev & Tsvetkov (2012) showed that the kinematic parameters   / and   /, determined separately for the Northern and Southern hemispheres using Zonal Vector Spherical Harmonics (ZVSH), have different signs but they are close in absolute values.However, if we consider the whole sphere, the values of these parameters are close to zero.In this work, we use a symmetric distribution of stars in both the Northern and Southern Galactic hemispheres.These parameters describe the difference between the vertical gradients in the Northern and Southern hemispheres -antisymmetric component of the vertical gradient.We estimate the values of the symmetric vertical gradient component using second order derivatives with respect to  (see details in the next section).
The values of the kinematic parameter   / vary within the range ±10 km s −1 kpc −1 and show a clear dependence on the coordinates  and  (Fig. 9, left).The values of   / are positive and show a slight increase with Galactocentric distance in the region with  > 160 • .This means that for this region of the Galaxy the absolute value of the vertical gradient of the Galactic expansion is larger in the Northern hemisphere than in the Southern.This conclusion can also be made from the results of Vityazev & Tsvetkov (2012), Nelson and Widrow (2022) and Wang et al. (2023).For example, Vityazev & Tsvetkov (2012) used the data of UCAC3 catalogue from (Zacharias et al. 2010).They obtained the values of   / in the region of the Sun +22.2 ± 0.5 km s −1 kpc −1 and −10.8 ± 0.5km s −1 kpc −1 for the Northern and Southern hemispheres, respectively.In the Galaxy region with  < 160 • , the kinematic parameter   / has a negative value and also shows an increase in absolute value with an increase in .
Fig. 9 (middle) shows the non-symmetry with respect to the Galactic plane of the vertical gradients of the Galaxy rotational velocity   / in the Northern and Southern hemispheres.The first proof of slowing down of the Galaxy rotation velocity with increasing distance from the Galactic middle plane was given by Majewski (1993); Girard et al. (2006), Makarov & Murphy (2007) with estimation of the corresponding parameter as +20 km s −1 kpc −1 .Vityazev & Tsvetkov (2012), using data from the UCAC3 catalogue, determined with high accuracy the values of the vertical velocity gradient of Galactic rotation of   / = +26.9± 0.6 km s −1 kpc −1 and −52.1 ± 0.6 km s −1 kpc −1 in the Northern and Southern hemispheres, respectively.It was made separately via the Zonal-VSH method.Later, Velichko et al. (2020) found similar results using the stellar proper motions of the catalogues: Gaia DR2 (Gaia Collaboration et al. 2018) and PMA (Akhmetov et al. 2017).The authors of these papers did not interpret the differences in the absolute values of   / for each hemisphere; they averaged this value and found its lower and upper estimations.Fig. 9 (middle) shows the difference in the absolute values of the vertical gradient of the Galaxy rotation velocity in the Northern and Southern hemi-spheres.In the region close to the Sun, the value of this parameter is about +7 km s −1 kpc −1 , which means that the slowdown velocity of the Galaxy rotation in the Northern hemisphere is less than in the Southern one.We can make such a conclusion from the analysis of (Nelson and Widrow 2022), their figures 4 and 6.We can see the same behavior of the asymmetry of the vertical gradient at Galactocentric distances from 0 to 8.5 kpc and beyond 12 kpc for almost all  angles.Possibly, the asymmetry of the vertical gradient of the Galaxy rotation in the region with large Galactic distances beyond 12 kpc can be also explained by the presence of the Galactic warp.
Fig. 9 (right) shows the dependence of the difference (nonsymmetry with respect to the Galactic plane) of the vertical gradient of the vertical stellar velocity between the Northern and Southern hemispheres.As can be clearly seen in this figure, this kinematic parameter is close to zero near the Sun.The value of   / begin to deviate significantly from zero only at distances greater than 10 kpc and its behavior looks like the parameter   / but with an opposite sign.We can also see an increase in the absolute value of this parameter with  and a sign change at  = 160 • .
If we do not consider the local variations of the parameter   / in Fig. 9 (middle), this parameter also shows a significant change near the line with  = 160 • .Thus, on this line we see a sign change in all kinematic parameters   /,   /, and   /, showing an asymmetry of the vertical gradients of the velocity components relative to the Galaxy middle plane (differences in the Northern and Southern hemispheres).Possibly, these vertical perturbations of the stellar velocity field are due to the influence of the Sagittarius galaxy and/or the Large Magellanic Cloud (Garavito-Camargo et al. 2019) and (Vasiliev et al. 2021).

Azimuthal stellar velocity as a function of radii: the Galactic rotation curve
Figure 10 shows the Galaxy rotation curves obtained by different methods from the azimuthally averaged values of the rotation velocity component in the range  from 150 • to 210 • .The black dots correspond to the value of the Galaxy rotation curve obtained by averaging the stellar velocity component   for all "fictitious stars" without using any kinematic model.The (blue points) represent the Galaxy rotation curve obtained from the Oort parameters  and  of equation ( 26), by assuming an axisymmetric rotation and   / = 0.The (green dots) show the Galactic rotation curve from equation (25), which takes into account the non-axisymmetric rotation   / from equation ( 28). Figure 10 also shows the Galaxy rotation curve    = −  derived from equation (29) (red dots).We note that these values are systematically larger in the range of Galactocentric distances from 4 to 9 kpc than those calculated by other methods.This is due to the vertical velocity gradient that is modeled by the quadratic term  2   / 2 , while the other methods can not correct the systematics due to the vertical gradient on the circular velocity, as the first order expansion of the circular stellar velocity along  is insufficient.The values of the Galactic rotation velocity and their errors are given in Tables A3 and A4.Thus, we argue that our model described by equation ( 29) (red dots) is able to determine the mean Galactic rotation velocity in the Galactic plane ( = 0) most accurately, because it takes into account both the influence of the vertical gradient and the non-axisymmetric rotation.Figure 10.Black dots correspond to the value of the Galaxy rotation curve obtained by averaging the stellar velocity component   for all "fictitious stars" without using any kinematics model.The Galaxy rotation curves calculated using our method from equation (29) (red dots), using O-M model with the Oort parameters  and  of the formula (26) (blue dots).We also present the rotation curve taking into account   / using eq.( 25) (green dots).
The values of the Galaxy rotation curves are obtained by averaging the corresponding data in the range of angles  from 150 • to 210 • (±30 • about the direction of the Galactic Centre -the Sun -the Galactic anticentre).We compare the rotation curves derived in this work with circular velocity curve   () derived by Eilers et al. (2019) (brown dots).

SECOND ORDER DIRECTIONAL DERIVATIVE OF THE STELLAR VELOCITY FIELD
In this paper, we present for the first time the results of the expansion of the stellar velocity field into the Taylor series up to the second order for Gaia data.As we show in the Appendix A of Fig. A3 and A4, the accuracy of the second derivatives is somewhat lower than that of the first, but all these parameters are estimated quite reliably at the level of ±0.5 -4 km s −1 kpc −2 at the heliocentric distance up to 5-6 kpc.We give below a short description and physical interpretation of only some of the most significant second derivatives.

Second order derivatives of the stellar velocity field in the radial direction
In our opinion, the most interesting are the second derivatives of the stellar velocity field along the Galactic radius  2   / 2 ,  2   / 2 , and  2   / 2 shown in figure 11.
As it can be seen in the left and middle panels of figure 11, these second order derivatives for radial   and azimuthal   components of the stellar velocity field show significant changes within ±10 km s −1 kpc −2 .Some patterns in the form of alternating ring-shaped structures are also clearly visible in the figure 11.We see alternating blue and red annular structures reflecting the local maxima and minima of the Galactic expansion velocity   and rotation velocity   on the scale 1-2 kpc, corresponding to the spatial resolution of our velocity maps.In our case, this behaviour is an analog of the Nyquist frequency.
Similar ring-like patterns are also present in the vertical velocity component  2   / 2 (Fig. 11, right panel).However, their amplitudes are significantly smaller than those observed for the other two velocity components, reaching values no higher than ±3 km s −1 kpc −2 .
The second order derivatives of the stellar velocity field in the radial direction show a wave-like dependence on Galactic distance .All these ring-like signatures evidenced in the second order derivatives of the stellar velocity field in the radial direction are probably related to kinematic substructures, such as the spiral arms and bulge, and/or by the propagation of density, bending and breathing waves.

Second order derivatives of the stellar velocity field along the azimuthal angle
Figure 12 shows the second order derivatives along the azimuthal angle.As can be seen from the left and middle panels of Fig. 12, the parameters (1/ 2 ) 2   / 2 and (1/ 2 ) 2   / 2 show a weak change that do not exceed ±3 km s −1 kpc −2 and do not have a clear dependence on Galactic coordinates.The parameter (1/ 2 ) 2   / 2 has an insignificant value practically in the whole range of Galactic coordinates (the right panel of Fig. 12).Thus, in the middle plane of the Galaxy, no clearly discernible wavelike behavior of the stellar velocity field is observed.Insignificant changes of the parameters (1/ 2 ) 2   / 2 and (1/ 2 ) 2   / 2 characterize the local change of stellar velocities and do not describe the global processes in the Galaxy.

Second order derivatives of the stellar velocity field in the vertical direction
The main value of the vertical gradient of stellar velocities results from the symmetry in absolute value and the difference in sign of the vertical gradient of stellar velocity fields in the Northern and Southern hemispheres.As Vityazev & Tsvetkov (2012) have shown, using a symmetric sample of stars relative to the Galactic middle plane, it is impossible to estimate the value of the vertical gradient within the framework of the O-M linear model or the VSH method.
Our models account for the possibility that not only the sign of the vertical gradient (symmetric component) but also its absolute value can change from the Northern to the Southern hemisphere -the antisymmetric component of the vertical gradient.Fig. 13 shows the values of the parameters  2   / 2 ,  2   / 2 and  2   / 2 corresponding to the symmetric component of the vertical gradient of the stellar velocity fields.As mentioned above, we have already considered the antisymmetric component   / and   / in Fig. 9.
The middle panel of Fig. 13 shows the map of the second order derivatives of the Galactic rotation velocity in the vertical direction.The value of  2   / 2 reaches 50 km s −1 kpc −2 in the Solar neighbourhood, that corresponds to an asymmetric drift of about −25 km s −1 at | | = 1 kpc.Such value is pretty consistent with the vertical gradient  rot = −20 ± 6 km s −1 kpc −1 and −22 km s −1 kpc −1 derived respectively by Makarov & Murphy (2007) and Levine et al. (2008), who analysed stellar samples close to the Galactic plane (| | < 1 kpc), as in our case.
In the range of Galactocentric distances  between 4 and 9 kpc, the value of the second order derivatives of the Galactic rotation velocity  2   / 2 is practically independent of the  angle and reaches a maximum value of about 50 km s −1 kpc −2 .At the Galactocentric distances  between 9 and 10 kpc, the value of the  2   / 2 quickly drops to 10 -20 km s −1 kpc −2 for all values of the  angle.Also, in the outer region with  > 10 kpc, we see a smooth decrease almost to zero at a distance of about 15 kpc.Such difference is due to the superposition of the thin disk with the slowly rotating thick disk in the inner disc, combined with the flaring of the thin disc in the outer regions (Mackereth et al. 2017;Beraldo e Silva et al. 2020).These results are consistent with previous studies carried out by Nelson and Widrow (2022, Figs. 4 and 6), Wang et al. (2023, Figs. 4-5), andGaia Collaboration et al. (2023, Fig. E.1).
The left panel of Figure 13 shows the Galactic map of the parameter  2   / 2 which, according to the physical meaning, is the second order derivative of the Galactic expansion velocity in the vertical direction.In the inner region of influence of the Galactic bar ( < 5 kpc), we see the parameter values equal to +30 km s −1 kpc −2 for  > 175 • , corresponding to an increase of the Galactic expansion velocity   with increasing height | | from the Galactic plane.In the region with an azimuth <175 • , the value of the parameter  2   / 2 has the opposite sign −30 km s −1 kpc −2 , indicating a decrease in the Galactic expansion velocity with increasing in | |.We also observe such behavior for  in the range between 5 and 10 kpc, but the value is 2-3 times smaller than in the region of the Galactic bar.
We note that, for  < 5 kpc, the effect of the parameter  2   / 2 is opposite to that of the Galactic expansion velocity   shown in the left top panel of fig.6, where   > ∼ + 20 km s −1 and < ∼ − 20 km s −1 for  < 185 • and > 185 • , respectively.This means that the radial stellar streams in the inner disc and bulge rapidly decrease with increasing | |.Thus, we point out that the second derivative  2   / 2 can provide important constraints on the dynamical models of the Galactic bar.
In Fig. 13 (right panel) the parameter  2   / 2 is the second order derivatives of the vertical stellar velocity in the vertical direction, which is positive in almost the whole Galactic plane and varies within the range from +3 to +10 km s −1 kpc −2 .Only in a small region with Galactic radii  between 6 and 8 kpc and angles  between 145 • and 240 • ,  2   / 2 has negative sign and changes in the range from −2 to −7 km s −1 kpc −2 .We can conclude that, in most of the Galaxy, the vertical velocity   slightly increases with growing distance  from the Galactic plane.As we can see from Fig. 13, the stars move on average from the Southern to the Northern hemisphere, slowing down slightly as they pass through the plane of the galaxy.The Sagittarius galaxy and/or the Large Magellanic Cloud may be responsible for these vertical motions.This will require a detailed study and dynamical modeling.

The mixed second order derivatives of the stellar velocity field
As we can see in figure 14, almost all mixed second order derivatives of the stellar velocity field have no obvious systematic behavior.The values of these parameters are close to zero within ±(3 − 4) km s −1 kpc −2 .Only two parameters  2   /  and  2   /  at the Galactocentric distance of about 12 kpc and azimuth angles  from 160 • to 200 • show a significant behavior and their values reach +10 km s −1 kpc −2 and −5 km s −1 kpc −2 , respectively.We can interpret the parameter  2   /  and  2   /  as gradients of the parameters   / and   / along the Galactic radius .
As can be seen from Fig. 9, only parameters   / and   / show a fast change in the region of Galactocentric distances =12 kpc.Apparently, this is due to the influence of the Galaxy warp into the stellar velocity field, as shown in (Poggio et al. 2018(Poggio et al. , 2020)).

SUMMARY AND CONCLUSIONS
In this paper, we have performed a kinematic analysis of the stellar velocity field for a sample of high luminosity stars from the Gaia DR3 catalogue covering a large area of the Galactic plane, up to 10 kpc from the Sun.However, we have restricted our interpretation to a heliocentric distance of up to 6 kpc, where our results appear to be most reliable.
We used a new statistical method based on the decomposition of the stellar velocity fields into a Taylor series up to the second order in the Galactocentric cylindrical coordinate system.We have shown that the kinematic parameters obtained by this method are in good agreement with the results obtained by means of the Ogorodnikov-Milne model, Oort-Lindblad, and Gaussian process models (Nelson and Widrow 2022).
Our method also provides additional kinematic parameters, including the gradients of the radial 1     and the azimuthal 1     velocity along the angle , that cannot be derived by the Ogorodnikov-Milne and Oort-Lindblad models.These two parameters are assumed to be zero by the axisymmetric Galactic models, so that significant non-zero values are a signature of non-axisymmetric structures.
We have computed the Oort constants , , , , and, for the first time, the second order partial derivatives of the stellar velocity field with respect to the Galactic coordinates (, , ).
The values of the kinematic parameters and the Oort constants for the Solar neighbourhood are listed in tables A1 and A2, respectively.
We have estimated the values of the vertical gradient of the Galactic rotation, expansion and vertical velocity using the symmetric distribution of the stars with respect to the Galactic plane.We have shown that the stellar velocity fields have different behavior in the Northern and Southern hemispheres and are in good agreement with the results of Vityazev & Tsvetkov (2012) obtained using ZVSH.
We have also determined by different methods the Galactic rotation curve for Galactocentric distances from 4 kpc to 18 kpc by averaging Galactic azimuths in the range -30 • <  < +30 • around the direction Galactic Centre -Sun -Galactic anticentre (fig.10).
In conclusion, the presented method provides a detailed description of the spatial stellar velocity field of the Milky Way and is able to reveal features related to non-axisymmetric substructures (e.g.spiral arms, Galactic bar, warp), as well as bending and breathing waves produced by satellite galaxies or other dynamical perturbations.Therefore, a detailed analysis of the new kinematic parameters in combination with dynamical modeling will significantly improve our understanding of the structure and evolution of our Galaxy.
The derived kinematic parameters and their errors can be made available to interested readers on personal request by e-mail: akhmetovvs@gmail.com.

APPENDIX A: THE DISTRIBUTION OF RANDOM ERRORS OF THE KINEMATIC PARAMETERS
In the estimation of the kinematic parameters, individual errors of the measured astrometric parameters were not taken into account.The standard errors of the model kinematic parameters were estimated using the diagonal elements of the inverse matrix of the normal equations and the unit weight error.We drew contour lines with a step size of 0.5 on the maps representing the random errors of the kinematic parameter estimates in the Galactic middle plane derived by solving the LSM the equations ( 28), ( 29) and (30). ⊙ =8.249 kpc and shown by the point ⊙.A1.The kinematic parameters for a stellar system close to the Sun with a radius of 1 kpc from equations ( 28), ( 29) and (30) .The velocity components for the stellar system are in km s −1 ; the first order velocity derivatives are in km s −1 kpc −1 and the second order velocity derivatives are in km s  A2.The Oort's constants in km s −1 kpc −1 for a stellar system close to the Sun with a radius of 1 kpc from equations ( 21),( 22),( 23) and ( 24) MNRAS 000, 1-?? (2022)

Figure 1 .
Figure 1.Selection of high luminosity stars, such as OBAF-type stars, subgiants and giants on the colour diagram:   −   -absolute magnitude   .We use stars above the red horizontal line with a value of   = 4.

Figure 3 .
Figure 3. Galactic maps of one-cell layer in Cartesian coordinates color coded according to the stellar count in the cell.The left shows the maps in Cartesian coordinates (x,y) on the galactic plane.The right top and bottom show the maps in Cartesian coordinates (x,z) and (y,z) in the y=0 and x=0 planes, respectively.The distance to the Sun  ⊙ = 8.249 kpc is marked by the point ⊙ here and in all the following figures.

Figure 4 .
Figure 4.The count of the cells (fictitious stars) in the selected spherical regions with radius of 1.0 kpc that were used for calculating the kinematic parameters.

Figure 7 .
Figure 7.The radial gradient of the velocity components (from left to right, respectively,    ,    ,    ) as a function of the Galactic coordinates.

Figure 8 .
Figure 8.The azimuth gradient of the velocity component (from left to right, respectively, 1

Figure 9 .
Figure 9.The vertical gradient of the velocity component (from left to right, respectively,    ,    ,    ) as a function of the Galactic coordinates.The figures show the value of the difference between the vertical gradient of the stellar velocity field in the Northern and Southern hemispheres of the Galaxy.

Figure 11 .
Figure 11.The second order derivatives of the stellar velocity field in the radial direction.The wavelike behavior of the stellar velocity field along the radial direction.

Figure 12 .
Figure 12.The second order derivatives of the stellar velocity field along of the azimuthal angle.The wavelike behavior of the stellar velocity field along the azimuth direction.

Figure 13 .
Figure 13.The second order derivatives of the stellar velocity field in the vertical direction.The symmetry part of the vertical gradient of the components   ,   ,   stellar velocity field with respect to the Galactic middle plane.

Figure 14 .
Figure 14.The mixed second order derivatives of the stellar velocity field as a function of the Galactic coordinates.

Figure A1 .
Figure A1.The random errors of the velocity components   ,   ,   as a function of the Galactic coordinates.

Figure A2 .
Figure A2.The random errors of the first order velocity derivatives estimates as a function of the Galactic coordinates.

Figure A4 .
Figure A4.The random errors of the mixed second order velocity derivatives estimates as a function of the Galactic coordinates.

Table A3 .
The values of the Galaxy rotation velocity are obtained by averaging the corresponding data in the range of angles  from 150 • to 210 • (±30 • about the direction of the Galactic Centre -the Sun -the Galactic anticentre) in km s −1 as functions of the Galactocentric coordinates  in kpc.