Modeling of Spiral Structure in a Multi-Component Milky Way-Like Galaxy

: Using recent observational data, we construct a set of multi-component equilibrium models of the disk of a Milky Way-like galaxy. The disk dynamics are studied using collisionless-gaseous numerical simulations, based on the joined integration of the equations of motion for the collision-less particles using direct integration of gravitational interaction and the gaseous SPH-particles. We ﬁnd that after approximately one Gyr, a prominent central bar is formed having a semi-axis length of about three kpc, together with a multi-armed spiral pattern represented by a superposition of m = 2-, 3-, and 4-armed spirals. The spiral structure and the bar exist for at least 3 Gyr in our simulations. The existence of the Milky Way bar imposes limitations on the density distributions in the subsystems of the Milky Way galaxy. We ﬁnd that a bar does not form if the radial scale length of the density distribution in the disk exceeds 2.6 kpc. As expected, the bar formation is also suppressed by a compact massive stellar bulge. We also demonstrate that the maximum value in the rotation curve of the disk of the Milky Way galaxy, as found in its central regions, is explained by non-circular motion due to the presence of a bar and its orientation relative to an observer.


Introduction
Attempts to understand the phenomenon of spiral structure in galaxies have a long history [1][2][3][4][5]. There is a general consensus that the spiral structure is a manifestation of density perturbations propagating in a multi-component stellar-gaseous disk. However, there is not agreement as to the mechanism for spiral generation, one that can successfully explain the rich morphological variety observed in galaxies. To describe the observed patterns, a few mechanisms have been suggested. Some researchers treat the spiral structure as the long-lived global modes that last in galactic disks for tens of galactic rotations [1,3,6]. Others consider the spiral structure to be a transient phenomenon, so that the spiral pattern changes many times during a galaxy's lifetime [7][8][9][10][11]. The tidal encounter of a small mass companion with a disk galaxy can also be a mechanism for the formation of spiral structure under certain conditions [12]. Khoperskov et al. [13] showed that a nonaxisymmetric halo can generate a large-scale spiral density wave in an embedded stellar disk. These authors showed that the pattern is observed during many revolutions, even if the disk is gravitationally stable.
Recent observations have provided a significant step in understanding the nature of the spiral structure in galaxies. Using stellar-cluster catalogs for three nearby spiral galaxies (NGC 1566, M51a, and NGC 628) from the Legacy ExtraGalactic UV Survey, Shabani et al. [14] measured the gradients in stellar ages across the spiral arms within these galaxies. The authors ρ * (r, z) = ρ * 0 (r) sech 2 where h * is the vertical scale height of the disk and * 0 is the central volume density. Deviations from (2) are determined by large number of factors, such as an inhomogeneity of disk velocity dispersion in vertical direction, multi-component nature of the disk that reflects its chemical and dynamical evolution, etc. An additional factor causing deviation from (2) is a dark matter halo, as seen from the equilibrium Jeans equation in vertical direction [46,47]: where G is the gravitational constant, ρ(z) is the volume mass density, ρ 0 = ρ(z = 0), σ is the surface density, z 0 = ∞ 0 (ρ(z)/ρ 0 ) dz, V c is the circular velocity in the plane z = 0, ρ h is the local halo density, and c z is the velocity dispersion of stars in the vertical direction. The solution (2) corresponds to the case c z = const, so that the expression in curly brackets is equal to zero. The initial equilibrium of the disk in our simulations was built using the iteration procedure that is described in [47], which was modified for presence of a thin gaseous component.
Juric et al. [48] find that the three-dimensional distribution of old stellar population of Milky Way can be represented by two subsystems: the thin disk with a radial scale length about 2.5 kpc and a vertical scale height of 300 pc and the thick disk with a scale length and height of 3.6 kpc and 900 pc, respectively.
The vertical scale height of the density distribution of the Milky Way disk is another important parameter influencing the gravitational stability of the system and, thus, in determining the morphology of the spiral pattern. To explore this, in our numerical simulations we vary the vertical scale height of the collisionless disk from 200 to 380 parsec.
The observational data show that the Milky Way scale height increases with radius [49]. Nonetheless, we build our equilibrium disk models with constant vertical scale heights. Figure 1 shows the observed rotational curve of the Milky Way disk based on various observations. There is fairly good agreement between the different determinations of the rotation curve in the range for r (6)(7)(8)(9)(10)(11)(12)(13)(14) kpc. However, the curves differ considerably in the inner regions of the Galaxy. The circular velocity V c and rotation velocity of the gas disk V g are shown: curve 1 is the circular velocity from WISE, 2MASS, and Gaia [50]; 2 is the result of compiling a large amount of data from various spectroscopic and trigonometric measurements from radio to optical wavelengths [51]; 3 is the reconstructed circular velocity using Gaia Collaboration data and simulations that are based on the Jeans equation [52]; 4-maser observations from [53]; gaseous rotation curves inside the solar circle [54]: atomic hydrogen (curve 5), molecular hydrogen (curve 6); the maser rotation curve from Reid et al. [55] is shown by line 7.

Rotation Curves
New kinematical data that are based on maser observations provide independent and promising data on the kinematical properties of the cold component of the Milky Way disk [53,55,56]. However, despite the fact that maser observations directly determine the parallaxes of sources, one has to keep in mind that the reconstruction of the Galactic rotation curve involves a number of assumptions, such as, for example, the model of the Galactic potential, which is necessary to reconstruct the rotation curve [57].
The Galactic rotation curve has a local maximum in the inner disk region, as observed in atomic and molecular hydrogen lines [58][59][60]. The maximum of the rotational velocity is probably associated with the non-circular motion of gas, as caused by the bar, in the central regions of the Galactic disk.
The angular velocity of the stellar disk is usually lower than that of the gaseous disk. Estimates of the rotational velocity of stars in the solar neighborhood vary from 200 km s −1 [61] to 237 km s −1 [62]. A recent estimate by [63] gives the value of the rotational velocity of the old stellar disk (G-and K-dwarfs) as V = 205 ± 20 km s −1 .
Outside the stellar disk, the rotation curve of the Milky Way slowly decreases with a radius reaching about 160 km s −1 at a distance of 100 kpc from the center of the Galaxy [64], which apparently points to the presence of the extended and massive halo of the Galaxy.    [38] taking into account a decrease by 16% for RGB-stars, 9-velocity dispersion within bulge region [70], 10-as based on Gaia-ESO Survey Guiglion et al. [71], 11 and 12 for giant stars from Gaia in layer −200 pc ≤ z ≤ 200 pc for negative and positive azimuths, respectively [71], 13 is the central velocity dispersion of bulge [72], 14 is the velocity dispersion from Gaia data analysis [50].

Velocity Dispersion of Stars
In general, the observational data can be approximated by the function c r = 91 [km s −1 ]/ exp(r/7.4 [kpc]). The velocity dispersion profile that was obtained by Tiede & Terndrup [73] gives good agreement with the data in the disk's central regions of the Galaxy, but in the solar neighborhood with R = 8 kpc the estimate by Tiede & Terndrup [73] gives the value 17 km s −1 , which is too low, and contradicts observations. The velocity dispersion in the central regions of the Milky Way disk is about c r (110-130) km s −1 [65,70,72] and is likely associated with the bulge stellar populations.
Using the Gaia-ESO spectroscopic survey, [71] determined the velocity dispersions (c r , c ϕ , c z ) for the thin and thick disk separately, within 6.5 kpc ≤ r ≤ 9 kpc. They found, for the thin disk, a velocity dispersion of c r = 33 km s −1 (See points 10 in Figure 2). Piffl et al. [40] used a sample of 200,000 red giant stars from the Gaia-ESO spectroscopic survey and estimated the radial velocity dispersion for the thin disk within 1.5 kpc of the plane to be c r = 34 km s −1 [40]. It is known that the presence of young stellar populations decreases the effective velocity dispersion in a stellar disk. A detailed study of this effect was conducted by Rafikov [74], who discussed the axisymmetric stability criterion of a disk consisting of a few stellar components and a gaseous component. We take this into account by using a multiplicative factor, c r → c r /1.1, based on a statistically averaged estimate of the velocity dispersion for the multi-component stellar disk (giants, main sequence stars, white and brown dwarfs): where σ l and c l are the surface density and dispersion velocities for the l-th component, respectively. Using data of [74] for nine types of stellar populations, we obtain the above mentioned decrease in the radial velocity dispersion of the stars. It is known that, for a collisionless disk in equilibrium, the condition is satisfied (κ is the epicyclic frequency, Ω is the angular speed). We use this condition to determine the initial dependence of the azimuthal velocity dispersion along the radial direction. Our numerical simulations show that this condition is satisfied during the disk evolution. Different authors provide values of 0.45 − 0.6 for the ratio c z /c r in the solar neighborhood. A recent estimate by [63] of the velocity dispersions c z = 16 km s −1 , c r = 35 km s −1 , gives a value for c z /c r of 0.46, which is what we used to build our equilibrium models. The ratio c z /c r = 16 km s −1 /35 km s −1 is equal to 0.46, based on SCUSS and SDSS observational data for G and K dwarfs of the thin disk [63]. Estimates by Katz et al. [25] and Adibekyan et al. [75] give c z /c r = 0.5 for this ratio. We choose a value of 0.5 in our simulations for the ratio c z /c r . We note that the kinematical properties of the stellar disk change with time during disk evolution, due to disk heating Hayden et al. [76].

Stellar Bulge and Bar
A number of observational studies demonstrate that the Milky Way has a boxy/peanutshaped bulge. Recent evidence points to the presence of an additional spheroidal component [77]. The structural properties and origin of the bulge are the subject of intense study. In this paper, we neglect the possible evolution of the bulge and consider it to be a fixed structure of the Galaxy contributing to the total gravitational potential. Estimates of the mass of the Milky Way bulge vary over a rather wide range, (0.9-2.0) × 10 10 M [78,79], which, in part, is a result of different assumptions regarding the stellar density distribution of the bulge. Estimates of the Galactic bulge's mass are based on photometric and kinematical data. A recent estimate for the bulge, based on photometric data [78], gives the total mass of stars and remnants in the bulge as 2.0 ± 0.3 × 10 10 M , which is consistent with the estimate by Portail et al. [79], who derive a total dynamical mass for the bulge equal of 1.85 ± 0.05 × 10 10 M . The kinematical estimates of the bulge's mass are based on the assumption that the maximum of the rotational velocity of gas close to the galactic center at 0.5 kpc is caused only by the galactic bulge. This gives a total bulge mass of approximately M b = 0.99 × 10 10 M [80]. We vary the total mass of the bulge in our simulations within (0.7-1.6)×10 10 M .
The presence of a prominent stellar bar with a large semi-axis of about (3.1-3.5) kpc and with an axial ratio of about 1:0.4 is another peculiarity of the Milky Way galaxy [77]. The major semi-axis of the bar deviates by approximately 23 • from the direction to the Galactic center [81]. A successful dynamical model of the Milky Way must agree with the available observational data and reproduce the observed parameters of the bar.

Gas Distribution
Gas plays a significant, possibly crucial, role in forming and sustaining the spiral structure in galaxies. Goldreich & Lynden-Bell [82] pointed to the importance of gas in the formation of spiral structure: S0 galaxies are topographically similar to normal spirals, but they have no gas, no dust, and no spiral arms. This suggests that stellar dynamics is not solely responsible for arm formation.
The total amount of gas within 30 kpc of the Milky Way disk is approoximately 8 × 10 9 M [83]. In our simulations, we vary the total amount of the disk's gaseous component from 3.5 × 10 9 M to 6.5 × 10 9 M within the optical radius of the disk, R opt = 4r d = 9 kpc, s listed in Table 1. As a result, the observationally limited ratio of masses of gas and stars in the Milky Way disk varies within µ g = M g /M d = (0.13-0.17). In our simulations, the gas surface density in the solar neighborhood was kept within the observed estimates of (H 2 , H I and H I I) σ g = 13.7 ± 1.6 M pc −2 [42]. The observational data for distant galaxies indicate a significant role of gas accretion, minor and major mergers, which makes the stellar-gas disk a non-conservative system, exerting a strong influence on the kinematics and morphology of galaxies [84].
The current computational possibilities allow following the dynamics of Milky Way-like galaxy for more than 10 Gyr. Recent cosmological simulations undertaken by a few groups for following the formation of galaxies in isolated Milky Way-mass are impressive haloes [85][86][87][88]. These simulations follow many aspects of galaxy formation, including black hole accretion and its feedback, feedback from massive stars, stellar and chemical evolution, metallicitydependent cooling, star formation, and an influence of magnetic fields, e.g., mock images of the Milky Way-like galaxy in FIRE-2 simulations demonstrate, at z =0, a developed multi-armed spiral pattern. Our simulations use a more simplified model when gas is treated as a quasiisothermal one, but we base the simulations of disk dynamics on the current observational knowledge of Milky Way stellar and gaseous rotation curves, stellar, and gaseous velocity dispersions, the density distributions of gaseous and stellar components in the Milky Way disk, as well as on the knowledge of size and density distribution in the stellar bulge of the Milky Way galaxy.

Dark Matter Halo
The equilibrium of the galactic disks depends on the dark mass distribution in halos [89,90]. The rotation of the Galactic disk is determined by the common gravity of the disk itself, the bulge, and of the massive dark matter halo. Once the rotation curve is specified, as well as the mass and density distributions of the disk and bulge, one can determine the mass of the dark-matter spherical halo. Within the optical radius of the Milky Way disk, the halo mass exceeds that of the disk (see Table 1). The Galactic halo reaches a mass of approximately (5.5-10.3)×10 10 M within R opt . In numerical models, with the value of R opt differing depending on the adopted radial scale length of the disk r d (See Table 1). Based on kinematical data for 200,000 giant stars, [40] estimated the mass of the dark-matter halo within the solar radius to be M h (r ≤ R ) = (6 ± 0.9) × 10 10 M . However, the galactic halo extends to much larger distances than the optical radius of the stellar disk (for a review see [89,91] and the references there).

Equations and Numerical Algorithm
We treat the dynamics of the 3D stellar-gaseous galactic disk self-consistently. The gaseous subsystem is modeled by N g Smoothed Particle Hydrodynamics (SPH) particles. The dynamics of the collisionless (stellar) disk is modeled using N * particles and a direct particle-particle integration scheme. We set the numbers of stellar and gaseous particles equal, N g = N * = N/2, where N = N g + N * is the total number of particles used in the numerical simulations.
Because the star formation is not taken into account in our model, the relative number of the gaseous and stellar particles does not influence the results of simulations. Therefore, based on the convenience of modeling, we chose N * = N g .
The equations of motion of gravitationally interacting SPH and collisionless disk particles in the external gravitational fields of an axisymmetric isothermal halo f h and bulge f b are as follows: where t is the time, ∇ is the Hamilton operator, the radius-vector r i (t) specifies the position of the i-th particle in space, while i and p i are the mass density and gas pressure of the i-th particle, respectively. We use the analytical halo model described in the works [92,93]. The gravitational interaction between the i-th and j-th particles is given by the equation: where m j and r j are the mass and radius-vector of the j-th particle, respectively, and δ is the gravitational softening length, preventing unrealistic accelerations during close encounters of the particles [94,95]. In our simulations, the number of particles is N = (2 19 -2 21 ), and the softening length is taken to be δ 40 pc. With the maximum resolution N * = N g = 2 20 in simulations listed in Table A1, masses of the stellar and the gaseous particles are within m * ∼ 2-3 × 10 4 M , m g ∼ 2.4-6.5 × 10 3 M .
In order to describe the dynamics of the gaseous disk, one must supplement the equations of motion with the equation of specific internal energy conservation, e i , and with the equation of state of the gas e = e(p, ): where v i is the velocity vector of the i-th particle and γ is the adiabatic index. We choose an equation of state for the gas close to an isothermal one with γ = 1.05 to decrease the effect of gas heating by shock fronts during the evolution of the gaseous disk. Such an approach is employed as the simplest cooling model in various studies of astrophysical gas dynamics [96][97][98][99]. A more realistic consideration of thermal processes requires a cooling function that depends on temperature, density, and gas metallicity [100][101][102][103]. A correct description of gas cooling is possible for multicomponent gas models that take chemical reactions and gas metallicity into account, as is done in the modeling of the molecular clouds in the disks of spiral galaxies [104][105][106]. An accurate accounting of gas cooling should also include radiation transfer (see, e.g., [101]), and we leave this for further consideration.
In accordance with the SPH-approach [107], the density of gas that is associated with the i-th gas particle, the equation of motion (6), and the energy conservation Equation (8) can be written in the following form (details can be found in [108]): Here, W is the Monaghan smoothing kernel [107] and W p is the smoothing kernel that is used for the approximation of pressure forces [108,109], and h ij = 0.5 (h i + h j ) is the effective smoothing length, where the smoothing length for each particle depends on its mass and is the symmetric approximation of the pressure forces, and ν a ij is the artificial viscosity [108]. We have, for the artificial viscosity, ν a ij : The empirical constants α, β, and η determine the value of the artificial viscosity. In our simulations ee used the values α = 0.5, β = 1, and η = 0.1.
If the smoothing kernel W (Monaghan cubic spline [107]) is used to calculate the pressure gradient, a nonphysical numerical clustering of particles will occur in high-pressure regions [111]. This is caused by the weakening of the interaction between the particles in the vicinity of 0 < ξ < 2 3 and lim i-th and j-th particles. To eliminate this and improve the stability of the numerical algorithm, we use the smoothing kernel W p taken from [109]: From Equation (13), it follows that lim For the numerical integration of the differential Equations (11) and (12), we use the predictor-corrector scheme of second-order accuracy (the so-called leapfrog method). We use the direct particle-particle algorithm to calculate the gravitational forces. The leapfrog method allows for us to simulate the dynamics of the disk systems, even in the cartesian coordinate system. Khrapov et al. [112] have shown that use of the leapfrog method in double-precision simulations conserves the total energy, momentum, and angular momentum of the equilibrium system with N = 2 20 particles with an accuracy of 10 −5 , 10 −15 and 10 −13 , respectively. For single-precision simulations, the accuracy of conservation of the above-mentioned quantities is equal to 10 −3 , 10 −2 and 10 −3 , respectively [see 112].
Details of the realization of the predictor-corrector method are described by [108,112]. Here, we outline the coordination procedure for simulations of gaseous and stellar disks. For the gaseous disk, the integration time step ∆t g (t) is limited by the stability condition of the SPH-algorithm [see e.g. 108], while, for the stellar disk, the time step is limited by the condition of applicability of Newtonian gravity, namely the time step in the collisionless simulations should be greater than the time of light propagation in the region of the simulations ∆t * ≥ ∆t crit . Here, ∆t crit is the propagation time of light within the region. If we choose ∆t * = ∆t g , the angular momentum conservation is satisfied with an accuracy of 10 −13 . However, in this case, the condition of applicability of Newtonian gravity fails, and the total integration time of the problem increases by factors of tens to hundreds. Therefore, we choose the value ∆t * = ∆t crit 2 × 10 5 years in our simulations. The calculation of the gravitational interaction between all of the particles was carried out once for a time interval (t, t + ∆t * ) using the expression (7). For the gaseous disk, the number of time steps is large (n g 1) for each time interval (t, t + ∆t * ), in order that the gravitational force vector is constant during each time interval, which leads to an accuracy of conservation of angular momentum of 10 −2 . The following correction procedure for the velocities of gaseous particles at the last time step ∆t g (t n g ) allows us to increase the accuracy of angular momentum conservation to 10 −8 : where v x , v y , and v z are the components of the velocity vector of gaseous particles in the Cartesian coordinate system (x, y, z), Here, (F R , F ϕ , F z ) are the components of the total gravitational force in the cylindrical system of coordinates (R, ϕ, z). For multiple GPUs, the details of a parallel OpenMP-CUDA implementation of SPH and N-body numerical algorithms are presented in the following papers [108,112].

Stability Criteria
Toomre [113] derived a stability criterion that determines the growth of small-scale spiral perturbations. The criterion is: where σ * is the surface density of stars, κ is the epicyclic frequency of the disk, and c r is the velocity dispersion of the stellar disk in the direction along the disk radius. For a gaseous disk, the criterion is: where c s is the sound speed of gas and σ g is the surface density of the gaseous disk. Disks of real galaxies are multi-component, and the formulation of a local stability criterion for such systems is not an easy task. A few different local stability criteria have been suggested. We use the two-component criterion, as suggested by [114]: We note that gravitating disks are unstable with respect to large-scale spiral perturbations (global modes), even when the local Q-stability parameter is greater than unity, which, in particular, is the case in our numerical simulations. The nonlinear stage of the large-scale instability that is discussed in Sections 4.2 and 4.3 re-distributes the surface density and velocity dispersion of the disk.

Simulations
It is convenient to choose the following units for the numerical modeling: With such a definition, the dimensionless gravitational constant and dimensionless mass of the stellar disk for any model will be equal to one for arbitrary values of the parameters k m and k r that were taken from the Table 1. We assume that the optical radius of the disk is equal to four radial scale lengths R opt = 4r d = r , where r d is the radial scale length of the thin stellar disk, meaning that the dimensionless value of the optical radius of the disk is equal to the one as well (R opt = 1). For example, for models 706 and 707, we have k m = 3.72, k r = 0.9 which gives m = 3.72 × 10 10 M , r = 9 kpc, v 133.7 km s −1 and t 63.2 Myr.
Henceforth, we will use the dimensionless units unless the dimension of the quantity is explicitly indicated.  When building equilibrium models of the Milky Way galaxy, it is important that the model parameters are consistent with observational data. Below are listed the specific parameters and equilibrium profiles that must be matched with the available observational data for the Milky Way.

1.
Rotation curve of the cold gaseous component V g (r).

2.
Rotation curve of the stellar disk V * (r).

3.
Velocity dispersion of the stars in the solar neighborhood c r together with its radial dependence c r (r).

4.
Surface density of the thin stellar disk in the solar neighborhood σ * .

5.
Radial scale length of the density distribution of the thin stellar disk r d . 6.
Surface density of the gaseous disk in the solar neighborhood σ g together with the surface density dependence on radius σ g (r). 7.
Vertical scale height of the thin stellar disk h.

8.
Size of the stellar bar r bar . 9.
Radial scale length of the bulge density distribution b and its mass M b .
We have constructed over one hundred numerical models of the Milky Way galaxy using different sets of equilibrium parameters. Table 1 lists the models satisfying the observational properties of the Milky Way. To discuss the results of the numerical simulations, we choose model 706 (see Table 1), which agrees with the available observational data for the Milky Way galaxy.
Equilibrium particle distributions were constructed using the force balance equations. For the collisionless (stellar) particles, the Jeans equations were solved along the radial and vertical coordinates. For the gas particles, the hydrostatic-balance equations using a polytropic equation of state were solved in the z and r -directions. Because of the assumptions used to constrain the equilibrium state, and due to the finite number of particles modeling both the collisionless and the gaseous disks, the equilibrium solution is approximate and it contains so-called numerical noise. Thus, the evolution of a stellar-gas system begins from a quasistationary equilibrium that includes practically a full spectrum of initial perturbations. This instability leads to the growth of the most unstable modes that form the observed spiral structure.

Spiral Structure
For the models considered, Figure 3 shows the typical radial dependence of the Toomre Qparameter (Q T (r)). There is a moderate increase of the Q-parameter in the region r d < r < 4r d with a broad minimum at the value Q T ∑ = (1-1.5). As is known from previous studies, in this case we have conditions for a spiral global instability and the formation of spiral structure in the disk, including a central bar. Figure 4 shows the typical temporal evolution of a stellar gaseous disk with equilibrium parameters in agreement with observational data. The formation of an open spiral/bar in the central regions of the disk occurs rather quickly after approximately 1.5 disk revolutions and, after approximately two disk revolutions, the bar is formed. Outside the stellar bar, a complex spiral pattern grows. Even in the stellar disk, with its more regular spirals, a superposition of two-and three-armed spiral patterns with different amplitudes grows and reaches the saturation stage. The gaseous disk exhibits an even more complicated spiral structure due to nonlinear interactions of the unstable modes, which leads to the branching of spirals, the appearance of rows, and that of a multi-level spiral pattern at the periphery of the disk.
Fourier analysis of the surface density perturbations provides a means for visualizing the growing spiral and bar-like structures. It is convenient to characterize the growth of perturbations by the time dependence of the global Fourier amplitudes for different azimuthal numbers m [94,115]: where I = √ −1, ϕ j is the azimuthal angle (radians) of the j-particle, and the summation is taken for the particles inside a radial layer r k ≤ r j ≤ r k+1 , N k is the total number of particles in this layer. Figures 5 and 6 show the time dependence of the amplitudes of the m = (2-4) harmonics for different regions of the stellar and gaseous disks. In the central regions of the stellar disk r < 2r d (see the left and middle panels in figure. 5), the dominant mode is m = 2 harmonics, related to the bar. In addition to the m = 2 mode, the disk is susceptible to the instability of m = 3 and m = 4 spirals, which are more noticeable at the periphery of the disk. The growth rates of perturbations in the central regions of the disk are noticeably larger when compared to that in the periphery of the disk. The temporal behavior of perturbations in the gaseous disk is qualitatively similar to that of the stellar disk. As one can see from Figure 6, in the central regions of the gaseous disk, the mode m = 2 also prevails, while, outside the disk's central regions, the perturbations are represented by a combination of m = 2, m = 3, and m = 4 spiral modes.
The number of spiral arms in the Milky Way galaxy, as well as their pitch angle, has been the subject of studies for a long time without consensus among the researchers as to the results. From observations of the stellar and the gaseous disks, estimates of the pitch angle of the Milky Way spiral pattern range from 5 • to 25 • [116][117][118]. Unsuccessful attempts to understand the properties of the global spiral pattern in the Milky Way disk have led to an approach where segments of the spiral arms are treated separately [26,118,119]. The disagreement between various determinations of the properties of the Milky Way spiral arms that are based on observational data (e.g., masers, ionized. neutral and molecular hydrogen gas, and 2MASS sources) that we suggest is a manifestation of the complex and non-stationary spiral structure of the Milky Way disk. We believe that this structure is a superposition of nonlinear spiral patterns with different azimuthal wavenumbers, angular speeds of the patterns, and amplitudes.    Here we have paid special attention to the reconstruction of the maximum of the Milky Way disk's rotational velocity in its central region r 1 kpc. Figure 7 illustrates the typical spatial distributions of the gas and of stars in models 401, 405, 455, and 474, where the rotational velocity in the disk has an internal maximum in the region r 1 kpc due to the bulge-disk mass distribution. We find that, due to a massive and centrally concentrated bulge, which is necessary for reconstructing the internal peak on the velocity rotation curve, the central bar either does not form or has a short major axis and is short-lived. Even if the bar appears at the initial stages of the disk evolution in a relatively cold stellar disk, it is destroyed at later stages. A strong gradient in the gravitational potential in the central regions of the disk has a well-known [120,121] destructive role on bar formation that is caused by the presence of a compact massive bulge. An additional factor suppressing bar instability is the increase of gas density in the central regions of the disk in the process of disk evolution. These arguments lead us to choose models in which the central peak of the disk's rotation curve is not explained by the presence in the disk of a massive and compact bulge. Instead, the central peak of the disk's rotation curve is formed by the process of nonlinear evolution of spiral/bar perturbations. Examples of these models are 600, 602, 610, 701 shown in Figure 7 and 706 in Figure 4. These models develop a pronounced stellar bar with gas kinematics that are similar to that observed for the Milky Way bar.
Outside the central regions of the stellar disk, models with a bar demonstrate a mixture of two-, three-, and four-armed spiral patterns of lower amplitude. The gaseous disks show a more complicated behavior, with the appearance of a non-stationary spiral pattern.
We note that models with a relatively low stellar velocity dispersion (e.g., model 701 in Table 1) develop unrealistically high-amplitude spiral perturbations. Actually, the development of a strong spiral perturbations leads to a heating of the stellar disk, and the velocity dispersion of the stars increases to approximately 40 km s −1 (see Figure 11) in a short period of time. Thus, the observed amplitudes of spiral perturbations in the disk of the Milky Way impose additional restrictions on the velocity dispersion of stars in the solar neighborhood. Models with a velocity dispersion c r = (29-33) km s −1 support a long-lived stellar bar with a semi-axis of ∼3 kpc, together with a spiral pattern that has observable amplitudes. Edge-on images of the numerical models demonstrate features that are similar to those observed in the central regions of the Milky Way, such as the so-called X-structure (see the edge-on images of models 401, 455, 602, 610, and 701 in Figure 7. The radial scale length of the stellar disk is of critical importance in the development of a bar in the disk's central regions. The models with a radial scale length of the stellar disk larger than r d = 2.6 kpc, such as model 730 with r d = 2.6 kpc or model 720 with r d = 3 kpc, lead to a relatively low surface density disk in the central region σ 0 700 M pc −2 . That, in turn, requires the presence of a massive dark matter halo to explain the observed disk rotation curve. Such models do not allow the formation of a bar-mode in the disk (see, e.g., model 730 in Figure 7). Besides that, for model 730, the values of the parameters Q T * and Q Tg are (15)(16)(17)(18)(19)(20)% larger within the optical radius when compared to models 706 and 707, which also leads to the stabilisation of the bar-mode. Attempts to explain the observed central maximum of the rotational velocity of the disk in these modes (curves 1, 2 in Figure 1) were unsuccessful. We then considered a set of models with a less concentrated bulge, which explained the appearance of the maximum of the disk's rotational velocity as a result of bar appearance and the nonlinear evolution of perturbations.

Density Profile and Rotation Curve
Here, we discuss the behavior of model 707 (see Table 1), which satisfies all of the limitations imposed by the observational data for the Milky Way galaxy, as discussed above. Figure 8 shows the decomposition of the rotation curve for this model.
As one can see from Figure 8, the nonlinear stage of instability leads to a steepening of the gas's rotation curve in the disk's central regions. However, note that the amplitude of the central peak in the gaseous rotation curve is smaller than that observed for the velocity peak in the rotation curve of the Milky Way gaseous disk. Figure 10 shows the azimuthal density profiles of the stellar σ * (ϕ) and of the gaseous disks σ g (ϕ) at three fixed radii r = 0.15 (inside the bar), r = 0.51 (at two exponential scale lengths of the disk), and at the periphery of the stellar disk close to the Sun's orbit r = 0.89. Dashed black lines correspond to the averaged values of the density over the azimuthal angle ϕ. Azimuthal density distributions demonstrate the features of the spiral disturbances. In the central zone, the stellar density distribution shows a pronounced two-arm mode. As for the gaseous component, in addition to the main mode, it shows the presence of two other modes at lower amplitudes. At two radial scale lengths of the disk (r 2r d ), the stellar system demonstrates a four-armed spiral structure. At the periphery of the stellar disk, the amplitudes of the density perturbations decrease considerably and remain within 10-20 percent. Gas shows more complicated behavior, which can be described as a superposition of two-, three-, and four-armed spiral patterns (see Figure 6).
The radial scale length of the stellar density distribution is one of the key parameters that determines the mass of the stellar disk and its stability properties. The estimated radial scale length of the Milky Way's stellar density distribution is about 2.25 kpc [36] and we use that value in our simulations. In the set of models with r d = 3 kpc and r d = 2.6 kpc, the surface density of the stellar disk decreases in the central regions to σ 0 = 481 M /pc 2 and σ 0 = 725 M /pc 2 (corresponding to models 720 and 730). In order to satisfy the observed rotation curves in these models, one has to increase the mass of the dark matter halo relative to the mass of the disk to µ 5. This, as a result, suppresses the formation of the central bar in the Milky Way galaxy. In other words, with the observed surface density of the solar neighborhood, the dynamical models do not allow us to reproduce the stellar bar in the central regions of the disk if the exponential scale length of the surface density distribution is more than 2.5 kpc.
We find that the dynamical models that conform not only to the observed rotation curve of the galactic disk, but also to the surface densities of the stellar and gaseous components in the solar neighborhood, together with their observed velocity dispersions, give a relatively small total mass of the disk M d = M * + M g 4.5 × 10 10 M within 2R opt = 18 kpc. Therefore, dark matter dominates within 2R opt = 18 kpc, composing about 74 percent of the total mass of the Milky Way galaxy.

Stellar Velocity Dispersion
As an initial condition, we choose a radial dependence of the velocity dispersion that corresponds to the observational data for the Milky Way (curve 6 on Figure 2). At the galactocentric radius of the Sun, this choice gives the value c r (30-32) km s −1 . The tangential component of the velocity dispersion in the collisionless disk c ϕ is prescribed according to Equation (5).
During the evolution of the perturbations, the velocity dispersion of the stars evolves. Figure 11 illustrates the temporal evolution of the velocity dispersion for a number of models.
As one can see, the typical behavior of the velocity dispersion is its slow evolution during approximately one billion years. At roughly this time, the formation of spiral structure and of a bar takes place, and it is accompanied by essentially no heating of the collisionless disk. At later times, some increase of the velocity dispersion of the disk is observed. Figure 11. Dependence on time of the velocity dispersion c r (t) for different models.
During the development of perturbations, the velocity dispersions deviate from those of the axisymmetric case. Figure 12 shows the dependence of the velocity dispersion c r as a function of the azimuthal angle ϕ at three different radii of the disk for model 707 at time t = 9. As one sees, large amplitudes of the density perturbations in the central regions of the disk are accompanied by strong variations in the velocity dispersion c r . The variations of the velocity dispersion are less prominent at the periphery of the disk.

Disk Kinematics in Central Regions
Because of presence of the bar, the kinematical properties of the stellar and gaseous disks inside r 3 kpc demonstrate some specific features. Figure 13 shows the rotation curves of the disk projected along lines representing different viewing orientations with respect to the semimajor axis of the bar. As one can see, the rotation curve of the disk along the semi-minor-axis of the bar has a local maximum inside 1 kpc from the center of galaxy, similar to the observed kinematical properties of the Milky Way disk within the central kiloparsec (see Figure 1). The estimates of the angular velocity and mass of the Milky Way bar [81] provide values of Ω bar < 50 km s −1 kpc −1 and M bar 2 × 10 9 M , respectively.
The conditions for the formation of a bar are rather sensitive to the mass of the dark matter halo µ and that of the bulge µ b , relative to the mass of the thin disk. An important parameter is the radial scale length of the disk density distribution together with the scale length of the mass distribution of the bulge, as mentioned above. In the models with µ 2, the bar mode does not develop due to the suppression of bar instability by the massive halo. The classical stability criterion of the global bar mode [122] requires that the ratio of the kinetic energy of disk rotation E rot to the total gravitational energy |E Φ | does not exceed some critical value (E rot /|E Φ |) crit 0.14. Thus, galaxies with a dark halo are more stable with respect to the formation of the bar mode. The authors [122] made a fundamental conclusion, that stabilization of the bar mode requires the value of the parameter µ within (1-2). Criteria of this kind have been tested using better N-body models and they show that the exact value of µ depends on a number of parameters of the galactic subsystems that determine the properties of the bulge, the gaseous disk, the radial profile of the dark matter, etc. [22,27,81,89,123]. Similarly, in the models with a moderate ratio of masses of the halo and of the disk µ 1.6, but with a massive and compact bulge b/R opt 0.014, the bar is not reproduced in the numerical simulations due to scattering of the orbits of the collisionless particles in the central regions where the gravitational potential has a large gradient. Thus, the presence of a bar in the Milky Way disk puts rather strict limitations on its equilibrium properties.

Discussion
Athanassola [124] showed that a bar grows faster in galactic models with a live halo as compared to models with a rigid one. Athanassoula attributed this effect to the destabilizing influence of the resonant halo particles. The total number of particles in Athanassoula's numerical simulations was of order of 10 6 with the mass of individual particle being about 2.5 × 10 5 M . This leads, apart from the resonant effects, to the generation of high amplitude noise in the disk component. Polyachenko et al. [123] returned to the problem of bar formation using one to two orders of magnitudes larger numbers of particles. They found that, while a live halo decreases the time of growth for the bar instability from 500 to 250 Myr, the pattern speed and other parameters of the bar remain approximately the same. Because we are interested in the late stages of disk evolution, the possible influence of a live halo on the growth rate of the bar instability is insignificant for our purpose.
A possible influence of the "live" halo on the dynamics of the disks can be obscured in the numerical simulations by the particle resolution when massive dark matter "particles" generate a high amplitude noise in the disk.
Long et al. [125] and Collier et al. [126] showed that spinning dark live halos suppresses the growth of stellar bars. Recent studies based on GAIA DR2 observational data [127] demonstrate that the Milky Way halo has close to zero rotational velocity with average value of 1 ± 4 km s −1 .
Chemin discussed the dependence of the observed rotation curve on bar orientation [128]. He used the results of high-resolution hydrodynamical simulations by Renaud et al. [129] and demonstrated that the rotation curve of the disk depends on bar orientation relative to an observer. We confirm this result in a more realistic model using combined, direct N-body-gas dynamical simulations.
Deg et al. [130] presented a new version of the GlactICS code, and tested it simulating the dynamics of a Milky Way-like galaxy. The authors used a two-component collisionlessgaseous model, and included bulge and dark matter halo in their simulations. Similar to our results, they find that, after approximately 1 Gyr, a bar formation with semi-axis of about 3 kpc occurs in the disk, and a multi-armed spiral pattern, which is more prominent in the gaseous component, is formed outside the bar region. Strictly speaking, to model the dynamics of a particular galaxy, one must guess its initial equilibrium, one that fits the final stage of observed equilibrium properties for the galactic disk. However, the disk evolution does not noticeably alter the initial density distribution of the massive stellar disk, and up to 0.5 Gyr does not affect considerably the density distribution of the gaseous disk, thus justifying our approach, as seen in Figure 9. Therefore, the set of unstable spiral modes is well-determined by the current state of a slowly evolving disk. The succcessful modeling of spiral structure in nearby spiral galaxies NGC 1566 [131] and NGC 5247 [93], which was done using the observed equilibrium properties of these systems, further justifies the validity of our approach.

Swing Amplification vs. Global Modes
In the review paper by Dobbs and Baba [7] the authors state "We argue that, with the possible exception of barred galaxies, spiral arms are transient, recurrent, and initiated by swing amplified instabilities in the disc". The swing amplification mechanism of spiral formation was suggested by A. Toomre (see, e.g., [132]) and it remains in many studies the governing mechanism for explaining the origin of spiral structure. However, there are a number of unanswered questions related to this approach. The observed amplitudes of the spirals is one of the oldest of issues with the swing amplification theory. Swing amplification can amplify the initial perturbations by about a hundred times ( [132]), which is not sufficient for explaining the observed amplitudes of the spiral arms in galaxies. Similar to other hydrodynamical or plasma systems, galactic disks oscillate, and they have their own global modes that can be unstable. Why are these global modes not considered in the explanation of spiral structure in galaxies? What happens to these global modes in an explanation involving the swing amplification mechanism? Why are global modes observed in barred galaxies described as the possible exceptionato the recurrent and swing amplified instabilities? Differential rotation and gravity also work in barred galaxies, so the swing amplified noise should be seen in these galaxies too.
It is assumed that the origin of the spiral pattern recently discovered in the galaxy M51a is caused by its satellite galaxy, NGC 5195. However, the closest approach between these two systems was about 900 Myr ago, according to Wahde and Donner [133]. Is the observed two-armed global spiral structure attributed to an interaction that occurred about one Gyr ago? Why would swing amplification, operating on much shorted time scales, be suppressed in this galaxy?
Sellwood & Carlberg concluded in their recent paper [134] "We argue that apparently shearing transient spirals in simulations result from the superposition of two or more steadily rotating patterns, each of which is best accounted for as a normal mode of the non-smooth disc." We completely agree with this statement, and will also cite here Vadim Antonov, the author of classical works in stellar dynamics (Antonov's theorems in Binney and Tremaine's Galactic Dynamics), who said: "The galaxies are similar to copper trampets". We can add here that some galaxies, like NGC 1566, "play a pure spiral tone", while others "play an accord" of spiral patterns.
In light of this, a number of questions arise: -What is the mechanism of spiral saturation at the nonlinear stage of instability? Laughlin et al. [135] concluded that nonlinear self-interaction of a growing mode is responsible for mode saturation, but further study is needed. -Why does the presence of gas in the stellar-gaseous gravitating disk extend the lifetime of the spiral pattern? Goldreich & Lynden-Bell [82] pointed out, back in the sixties, the importance of gas to the formation of spiral structure: S0 galaxies are topographically similar to normal spirals but they have no gas, no dust, and no spiral arms. The importance of gas in sustaining the spiral structure was also confirmed in numerical simulations [93]. In light of new observational data, the spiral structure in galaxies continues to pose questions to theory.

Summary
To model the dynamics of a Milky Way-like collisionless-gaseous disk, we use, for the first time in such a study, direct integration in calculating the gravitational potential, contrary to the approximate tree-code or particle-mesh codes that were used before. It is not clear whether an approximate treatment of the gravity can correctly reproduce the dynamics of the multi-component disks over cosmological time intervals. An agreement of simulation results in the hydrodynamical, collisionless, and the collisionless-hydrodynamical approaches [93] provides validity to our model. As for particular results, we demonstrate that: -An axisymmetric two-component gravitating disk with parameters close to those of the azimuthally averaged Milky Way galaxy-in terms of observed rotation curve, velocity dispersion profile, and masses of the bulge and of the stellar and gaseous disk components-is unstable towards near-exponentially growing spirals having numbers of arms m = (2)(3)(4).
-At the nonlinear stage of instability, the spirals saturate at few tens of percent in the central regions of the disk, to a few percent at the disk's periphery.
-At the nonlinear stage, a prominent bar is formed in the central regions of the disk with a large semi-axis of about 3 kpc. -Outside the bar region, a complex spiral structure, represented by a superposition of two-, three-, and four-armed spiral patterns, rotating with different angular velocities, is formed. The spiral structure of the Milky Way galaxy is interpreted in some papers as single spiral pattern with a fixed number of arms and spiral pitch angle, and one rotation resonance with fixed position. We show that this is not the case for the Milky Way-like disk. -We demonstrate that the peak in the rotation curve of the disk of the Milky Way galaxy, which is located in its central regions, is a result of the non-circular motions caused by the bar which develops in the disk.
-We confirm that the presence of a massive and centrally concentrated bulge prevents the formation of a bar, in agreement with [122].
-We also show that the presence of gas of about ten percent the disk's mass extends the lifetime of the spiral structure to a few Gyr as compared to what is found in purely collisionless models.
This conclusion though should be validated in further studies with a multi-phase gas model using proper feedback, which is beyond the modelling capabiities of this study.
Author Contributions: Conceptualization, methodology: V.K. and A.K.; Writing-original draft preparation, V.K. and A.K.; Writing-review & editing, all authors; software and validation, S.K.; numerical simulations, S.K. and A.K.; visualization, S.K. and A.K.; supervision and project administration, V.K. All authors have read and agreed to the published version of the manuscript.

Funding:
The development of software for modeling galaxies and numerical experiments were funded by the Ministry of Science and Higher Education of the Russian Federation (the government task no. 0633-2020-0003, S.K. and A.K.).

Institutional Review Board Statement: please add
Informed Consent Statement: please add Data Availability Statement: please add Acknowledgments: Authors thank T. Girard and the referees for valuable comments. Numerical simulations were carried by using the equipment of the shared research facilities of "Supercomputer Center of Volgograd State University". V.K. acknowledges financial support by Southern Federal University, 2020 (Ministry of Science and Higher Education of Russian Federation).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

COBE / DIRBE
Cosmic Background Explorer / Diffuse Infrared Background Experiment SPH Smoothed-particle hydrodynamics OpenMP-CUDA Open Multi-Processing -Compute Unified Device Architecture VLBI Very Long Baseline Interferometry