Application of the viscosity-expansion method to a rotating thin fluid disk bound by central gravity

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The 2D rotation of a thin fluid disk with a porous sink around the center is studied within the Navier–Stokes and Euler equations. The viscosity (ν)-expansion method is applied to the viscous fluid bound to the central mass via gravity. The Navier–Stokes equations yield various types of rotation curve, including a flat one, depending on the choice of the pressure function that is not determined within the fluid dynamics. Stationary flow is achieved through the balance of the pressure gradient, gravity, and the centrifugal force. These features of the stationary flow survive in the inviscid limit. The stability of the inviscid flow is examined by the Euler equations for the perturbations. At large distances, the real part of eigenfrequencies (EFs) are dominantly positive and decreasing with distance for flat and rising rotation curves, meaning that the spiral pattern of the perturbations is trailing. One real increasing EF exists for the decaying rotation curve, for which the spiral pattern is leading. Complex frequencies always emerge when the disk has m-fold rotational symmetry with m ≥ 2. The shape of the perturbed rotation curve has azimuthal as well as temporal dependences. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .


Introduction
Spiral galaxies are vortices in which most of the masses due to visible matter are concentrated in the bright central core and the remaining visible masses form the arms in the disk.The mass distribution of this type is believed to be unable to account for the observed motions of stars in the disk of the spiral galaxy, whose azimuthal velocity distribution does not decay with the distance from the center but remains almost constant or even increases.Such an anomaly in the rotation curve was first noted by Babcock [1] for the Andromeda nebula, and has been continuingly confirmed by the improved observations to be common in many spiral galaxies [2][3][4][5][6].
The simplest way to avoid the above difficulty in understanding the observed rotation curve within the framework of the Newtonian gravity theory is to assume the existence of dark matter, which is supposed not to emit radiation while being a few times as massive as the visible matter in a galaxy.The idea of dark matter seems compatible with the assertion of the Virgo project [7,8] that ubiquity of dark matter enables the large-scale structure in the present universe to develop within the age of the universe.This also reminds us of the "dunkle Materie" suggested by Zwicky [9] in relation to the large velocity dispersion of galaxies in the Coma cluster.
Explaining the observed galaxy rotation curves without invoking dark matter has been attempted by modifying the Newtonian gravity theory in the low-acceleration regime [10,11] or by extending the Einstein gravity theory so as to incorporate new fields other than the metric tensor [12,13].Although PTEP 2015, 073J01 K. Takahashi these works may be successful in achieving their primary purpose, in this paper, we will focus our attention on the traditional Newtonian dynamics viewed from a somewhat novel standpoint.
Spiral galaxies have winding arms of interstellar gases and stars that together form a disk that is relatively thin compared to its radial size.The metamorphosis from a three-dimensionally extended object to a quasi-2D disk may be driven by the rotation under gravity.The initially rotating dark matter will also eventually form a disk as long as it interacts with each other and with visible masses via gravity.If the dark matter and the interstellar gas are cold enough, the motion of the system will be viewed as almost free fall accompanied with a small energy loss due to occasional inelastic collisions of constituents.Then the motion of the galaxy matter including dark matter as a whole will be amenable to the dynamics of a nearly inviscid fluid.
In an astronomical system of large scale, the event of the encounter of two objects may be rare.Therefore, in studying astronomical clusters of gravitating objects, it is usual to exploit the collisionless Boltzmann equation for the distribution function in the phase space or the Jeans equation in the configuration space (see, e.g., Ref. [5]).The fluid dynamical counterpart for the latter is the Euler equation for the inviscid fluid.This was the fundamental idea in the pioneering works on the formation of arms in a spiral galaxy by Lindblad [14,15], Toomre [16], Lin and Shu [17], and Goldreich and Lynden-Bell [18].On the grounds of hydrodynamics, they showed that the instability induced by gravitation sets forward the formation of a spiral-armed structure of a fluid disk via density waves [17,18].
At first glance, the Reynolds number of a galaxy composed of dark and visible matter is large enough to validate the application of the Euler equation to a galaxy.The flows of fluid with a large Reynolds number are, however, not always described by the Euler equation, even if one lets the Reynolds number be infinity.(Recall, e.g., that large Reynolds numbers and finite viscosity lead to the emergence of turbulences through the cascade of energies.Namely, fluid dynamics is singular at zero viscosity, as von Neumann [19] summarized compactly.)During the cosmologically long term, the astronomical system will experience various processes that bring about energy dissipations.The momentum transport under a non-uniform velocity field v will yield viscosity that is expressed for the Newtonian fluid by the Laplacian term ∇2 v in the Navier-Stokes equations.Gradient in collisions with finite cross section, σ , will be another cause of small but finite viscosity.Assuming that it is proportional to σ , a dimensional analysis together with the rotational symmetry suggests, e.g., σ ∇ • v∇ 2 v, σ ∇ω 2 (ω ≡ ∇ × v is the vorticity) etc. for the variation rate of the velocity field.In the case in which the quadrupole moment in the mass distribution is significant, the gravitational wave will also cause energy dissipation, which will be of a still higher power of v. Whatever the mechanism is, it is likely that dissipation will govern the final state of the galaxy as long as the dynamics is singular at zero viscosity.That the majority of spiral galaxies have not leading but trailing spirals seems to support this point of view [20].Various processes of dissipation as the origin of viscosity, albeit "small", will be active during the course of the evolution of an astrophysical fluid.
How can small but non-vanishing viscosity give rise to distinctly observable phenomena?For vortical flows, the final patterns of the swirling are determined by the component of the velocity field that vanishes in inviscid limit [21].Considering this effect of viscosity, the astronomical fluid may pertinently be treated as a viscous fluid.As was mentioned above, various physically different processes can give rise to viscosity.Among them, the momentum transport is the simplest one.Thus we are led to apply the Navier-Stokes equations to a rotating disk of a viscous fluid bound by gravity and observe the pattern of rotation of the disk in order to observe the residual effect of energy dissipation.
The flow that is to be considered in this paper has an infinitesimally small viscosity but is not described by the Euler equation.We may call such a flow a non-Eulerian inviscid flow (NEIF).The idea of the NEIF is substantiated by subjecting the physical quantities in fluid dynamics to the Maclaurin expansion in the kinematic viscosity, ν, provided that all of the physical quantities are analytic functions of ν.This presumption immediately leads to the ν-expansion method in fluid dynamics.Although careful treatments of the Navier-Stokes equations are generally required because of their singular nature [22,23], infinite kinds of new vortex solutions, including the missing link between the Burgers vortex [24] and the Sullivan vortex [25], have been found by this method [21].
As was noted in Ref. [21], the significance of the ν-expansion lies in that it relates the finite-velocity component of a steady and axisymmetric vortex directly to another component that vanishes at zero viscosity.This is the consequence of retaining the shear stress term in the Navier-Stokes equation that is neglected in the Euler equation or the Jeans equation and is one of the manifestations of the idea that some characteristics at finite viscosity are preserved during the history of letting the viscosity approach zero.Considering their long histories, the spiral galaxies, when viewed as fluid, would also have been influenced by non-vanishing viscosity during their formation processes.
The existence of spiral galaxies suggests that the rotating uniform fluid disk may be unstable.Elucidating the origin and nature of the instability together with its interrelation to the spiral structures is one of the most important problems of galaxy dynamics.In addition to the momentum transfer mentioned above, the angular momentum transport will play a role in the evolution of the fluid disk [21,22] (for a review, see Ref. [5]).Theoretically, the normal mode expansions and/or Wentzel-Kramers-Brillouin (WKB) approximation within the thin-disk model are prevalent [16][17][18][26][27][28][29][30][31] (for application to Saturn's ring system, see, e.g., Ref. [32]) but are accompanied by an embarrassing problem of consistency between ansatz and outcome.This problem may be avoided by a stricter specification of such model contents as mass and velocity distributions and so on [33], with a cost of obscuring the dynamical basis of the presumed simplifications.
The purpose of this paper is to explore within the ν-expansion scheme a thin swirling NEIF bound by gravity exerted by a central mass.Being distinct from the standard method mentioned above, our method will be expected to provide in some cases exact results for the stability of the system.In the next section, we briefly review the ν-expansion method.In Sect.3, the method is applied to a thin gravitationally bound disk and a family of new solutions for vortex motions and mass distributions is presented.In Sect.4, the stability of the vortex found in Sect. 3 is investigated by the linear perturbation method.In Sect.5, the patterns formed by perturbations are discussed.Section 6 is devoted to a summary and remarks.

Navier-Stokes equations for the large Reynolds number
The model consists of a single fluid of very small viscosity, reflecting the fact that the dark matter, the main ingredient of the galaxy, may be almost inviscid.The basic equations for the velocity field v of a fluid with density ρ and kinematic viscosity ν are given by ) There exists another way to take the limit of zero viscosity.Let us divide both sides of (2.1a) by ν and then let ν approach zero: Here the velocity field and the pressure are supposed to be functions of ν.This equation tells us that v at ν = 0 on the r.h.s. is related to quantities of O(ν) in the braces on the l.h.s. of (2.2).If (2.2) leads to equations with nontrivial solutions, they will describe the NEIF whose pattern is governed by the stress term.This is in fact the case for the well known Burgers [24] and Sullivan [25] vortices with the free parameter in those solutions being appropriately scaled.For other vortex solutions with a similar property, see Ref. [21].In the following sections, we will see how (2.2) is realized in the rotation of a fluid disk.
together with the mass-conservation equation All of the physical quantities except ν are functions of the 3D coordinate r = (r, θ, z) as well as the time t.
PTEP 2015, 073J01 K. Takahashi Suppose that the flow and the forces are steady and axisymmetric.In this case, the terms of time and angle derivatives are dropped and we have where . Now, in the ν-expansion method, physical quantities are supposed to be expressed by regular functions of ν.In the present problem, v r and v z are supposed to be odd in ν, while the remaining quantities are even.Regarding the density and gravity as independent of ν, we subject v and p to the expansions and substitute (3.3a)-(3.3c)for (3.2a)-(3.2d).By matching the coefficients of the power series in ν, we have equations for the expansion coefficients in (3.3).
The assumption of the expansions (3.3) is based on the fact that Eqs.(3.2) for the steady and axisymmetric flow are invariant under the "ν-inversion": The positivity of ν in the real world implies that the velocity gradient acts to decrease itself or that the energy dissipates.In order to maintain the steadiness of the motion, the energy lost is supplied by the presence of the pressure gradient and the external force (or by heat transfer, if any).The negativity of ν implies the opposite, i.e., that the velocity gradient acts to increase itself or that the energy concentrates.In order to maintain the steady flow, the increment of the energy must be consumed somehow.This can be done by inverting the motion and flowing against the pressure gradient and the external force.Since the system is axisymmetric, the energy-consuming motion is achieved by inverting the signs of v r and v z in (3.2).This is the meaning of the ν-inversion invariance.
The system we are considering is a thin and axisymmetric disk with an inner edge of radius r 1 > 0. There may exist a porous sink or a porous source at r = r 1 .That the fluid disk is thin and axisymmetric is embodied in the simplest way by letting ρ(r) have the factor δ(z), Dirac's delta function, as where r is a distance in the disk plane from the center of the disk.In the following, we scrutinize the resulting equations in the order of the power of ν.

0th power of ν
Representing the lowest-order balance of forces and the pressure gradients in (3.2a) and (3.2c), we have Equation (3.6a) together with (3.5) implies that p n (r), n = 0, 2, . .., also has the factor δ(z): Noting that f Gz given by (2.1b) and (2.1c) is odd in z and since the system is confined on the plane z = 0, (3.6b) merely states that 0 = 0. Equivalently, the same identity is derived for (3.7) by integrating both sides of (3.6b) over a small interval of z containing z = 0. Equation (3.6a) is used to determine p 0 (r ) as once ρ(r ), v θ 0 (r ), and f Gr (r ) are known.

1st power of ν
The relevant equations are (3.1b) and (3.1d), which lead to Note that, in this (and the following) order, the external force is not involved in the equations, which is an aspect peculiar to the ν-expansion method.The equations have been reduced to those of the standard fluid dynamics.Since the fluid is confined in the disk at z = 0, it seems to be natural to posit that Note, however, that the velocity field at the region where matter does not exist is not a physical observable in fluid dynamics.In other words, the field at vanishing density is allowed to take any value.Thus (3.10a) can be generalized to This latter condition is fulfilled in the axisymmetric vortex solutions found by Burgers [24], Sullivan [25], and Takahashi [21] and is anticipated to preserve the possibility of natural continuation to the 3D flows.Either of (3.10a) and (3.10b) will work in our scheme.With the delta function in (3.7) being a mathematical idealization of the thin but finite thickness of real disks, the condition (3.10b) may be more appealing.However, for simplicity, we first put forward the argument with the simpler condition (3.10a).We will return to (3.10b) later.Then (3.9b) leads to v r 1 is also a function of r only.Note that the above relation is independent of ν.Next, we solve (3.9a) with (3.10a) by keeping the derivative term with respect to z in the Laplacian, because the observable is not the velocity field itself but, e.g., the one multiplied by the density and there is no a priori reason to omit the z-variable as far as (3.9a) is concerned.Separating the variables by and substituting (3.12) for (3.9a), we have where K 1 is a real constant.The prime stands for the derivative with respect to the variable of each function.Equation (3.13b) immediately yields If it were not for the K 1 term, (3.13a) would be solved as This is nothing but the solution that is shared in the steady axisymmetric vortex solutions of the Navier-Stokes equation with no boundary [21,24,25].These vortices have v r 1 (r ) that is negative and linear in r at sufficiently large r , so that f behaves as 1/r for large r .In the case of non-vanishing K 1 , f behaves in an entirely different way, as we shall see below.
Obviously, the azimuthal velocity component v θ is one of the most fundamental quantities in understanding the rotating disk.The important fact that we must notice is that the nature of v θ is governed by v r , which vanishes in the inviscid limit, as (3.13a) together with (3.15) for a special case shows.On the other hand, v r is a self-determining quantity when viscosity is nonzero (see (3.17)).Thus, we first have to know v r in order to elicit the mechanism by which the pattern of the rotation is determined.

2nd power of ν
Extracting the ν 2 terms in (3.1a) and (3.1c), we have In deriving (3.16b), use has been made of the fact that p 2 is a function of r only.v z1 = 0 is the solution of (3.16b).
If the second-order coefficients v θ 2 and p 2 identically vanish, (3.16a) reduces to 7/20 where c is an arbitrary real constant.Physically interesting solutions to (3.18) exist for c > 0. In this case, the solution behaves as v r 1 takes negative values, which means that the flow swirls in toward the origin when ν is not zero.Together with (3.11), (3.19) implies that ρ ∼ 1/r 2 near the origin and ∼1/r at infinity.The density function is singular at r = 0.This is why we assume a sink (porous surface) at r = r 1 .The mass M(r ) within the radius r varies as ∼ ln r near the origin and ∼r at large distances.
The above result for the behaviour of ρ at r ∼ 0 is a direct consequence of setting v z1 (r) ≡ 0 in the whole space.This constraint can be relaxed by adopting the condition (3.16b) instead and by assuming the form v z1 (r) = zχ(r ).Together with (3.9b), we have ) Owing to (3.19), (3.20b) has two independent solutions that are finite in the whole region of r .They behave near r = 0 as where b is a constant.Numerical calculations show that, depending on the boundary condition, χ can be positive or negative in the region of r that we are interested in.This means that, when the viscosity is nonzero, the velocity field outside the disk directs toward or away from the disk plane.Now, (3.20a) yields ρ ∼ r −2+2b/c or r −2 e −b(ln r ) 2 /2 for r ∼ 0. (3.21b) That the behaviours of ρ and χ are correlated shows that the vestige of the effect of nonzero viscosity remains in the spatial distribution of density in the inviscid limit.For r → ∞, χ decays as exp − (2c) 1/2 r .Accordingly, we again have ρ ∝ 1/r and M(r ) ∝ r for large r , irrespective of the profiles of the rotation curve f .Equations (3.13a) and (3.18) are simultaneously solved by a numerical method.Discarding the solution in which v θ 0 diverges at r = 0, the results are depicted in Fig. 1 for some values of K 1 /c.For all cases, v θ 0 increases linearly near r = 0, i.e., a uniform rotation.For not too small K 1 /c, v θ 0 does not decay at large distances but exhibits exponential growth ∼ exp ((2c In an intermediate region, it maintains moderate values.If K 1 is negative, f exhibits an exponential decay at long distances.The positivity of K 1 seems crucial for the rotation curve of the disk to behave in a way consistent with the observations of spiral galaxies. We have so far considered the case that p 2 (r ) identically vanishes.Since p n≥2 (r ) are not uniquely determined within the ν-expansion scheme that exploits the Navier-Stokes equations and the continuity equation only, it is worthwhile to try another choice for p n (r ), which corresponds to a different boundary condition.For this purpose, let us recall that the Navier-Stokes equation allows 3D nonsingular axisymmetric vortex solutions for an incompressible fluid under the choice ∂ r p 2 (r ) ∼ r at infinity [21].In the present problem of a compressible fluid, the simplest choice that generalizes the where K is constant.Then (3.16a) becomes Equations (3.20) still hold for v z1 .Choosing the sign of K to be positive, this gives the solution for v r 1 that corresponds to the Burgers vortex [24]: Different solutions corresponding to other known vortex solutions will be possible [21,25], but, for simplicity, we restrict ourselves to the solution (3.24).Equation (3.24) implies that ρ(r ) and M(r ) are given by m/r 2 and m ln r , respectively, with m > 0. Of course, this makes sense as long as M(r ) is positive.Now, (3.13a) becomes Apart from the normalization, f is a function of √ K r and K 1 /K , and asymptotically behaves as r K 1 /K −1 .Some solutions that are linear near the origin were numerically sought and the results are depicted in Fig. 2. We see that the flat rotation curve together with concave and convex ones are reproduced in accordance with the relative sizes of K 1 and K .
The equations for the third or higher order of ν are derived in a similar way.However, these equations are trivially satisfied inasmuch as the coefficient functions of second or higher order except for p 2 (r ) are all zero.Thus the equations that we have dealt with so far are complete in the ν-expansion scheme.The velocity of the NEIF obtained by ν → 0 is represented by v θ 0 (r, z = 0)[= f (r )] only.
In this section, we saw that there exists a class of solutions that exhibit flat or almost flat rotation curves.This was possible because v θ 0 was oscillatory in the z-direction with the period l z = 2π/K 1/2 1 (see (3.14)).The velocity field outside the disk has by itself no physical significance in the razor-thin-disk model considered here.However, δ(z) in ρ(r) will always be extended, e.g., as exp −z 2 /ε /(π ε) 1/2 , and thereby the velocity at z = 0 comes to bear a meaning in a compressible fluid.On the other hand, the characteristic length within the disk may be given by l r = 2π/K 1/2 .The profile of v θ 0 with l r ∼ l z resembles the observed rotation curves of spiral galaxies, as is shown in Fig. 2.   Now that v r 1 and v θ 0 are determined, the density ρ and the pressure p 0 can be calculated by (3.8) and (3.20) once the strength of the central gravity f Gr is given.By parametrizing f Gr by

9/20
we here set α = 2 in order to gain a gross measure on the behaviours of ρ and p 0 .We will later discuss the physical role of the parameter α (see Eq. (4.5c)).The differential equations are solved in 0.3 ≤ r ≤ 12 by choosing boundary conditions that result in an asymptotically vanishing pressure gradient, although this is not an absolute requisite.Some results of calculations are shown in Fig. 3 for three χ .The density decreases monotonically with r in the whole region except for a very vicinity of r = 0.3, while the pressure first decreases to a minimum near r = 1.0, then increases and approaches a constant value.(This minimum point shifts inward as α decreases, although this is not shown here.)This implies that the pressure is not determined solely by the density.In other words, the fluid is not barotropic.This hollow in p 0 becomes shallower when χ ini = χ(0.3) is negative and approaches a certain value around −0.12.However, the asymptotic value of p 0 cannot be arbitrarily small because we are neglecting the disk contribution in (2.1b) while assuming the flat rotation curve.
In the Jeans equations, the pressure is given by the isotropic velocity dispersion v 2 .Its characteristic value in the disk may be grossly estimated by p 0 / ρ for large r .In our model, its possible PTEP 2015, 073J01 K. Takahashi minimum value will be found by setting the minimum of p 0 to be zero.The values of p 0 defined in this way at, e.g., r = 6 are 0.0015, 0.04, and 0.26 for the solid, dotted, and dashed curves in Fig. 3, respectively.In this case, we have v ∼ ( p 0 /ρ) 1/2 | r =6 ∼ 6.1, 3.7, and 10 for the above values of p 0 .These values are very large as compared to the asymptotic value of the azimuthal velocity v θ 0 ∼ 1.6 of the flat rotation curve (see Fig. 2) and the resultant p 0 still seem insufficiently small.Incorporating the self-gravity due to the disk fluid, the second term in (2.1b), may be necessary to construct the disk rotation with sufficiently small velocity dispersions.

Linear perturbation and stability of the disk
The r -dependences of the azimuthal velocity shown in Figs. 1 and 2 do not exhibit power law decays with r characteristic of the Kepler motion and therefore are in gross agreement with the observation of the galaxy rotations.However, the rotation curves of the majority of real galaxies have more or less oscillating structures [4].On the other hand, it is known that the rotating axisymmetric fluid disks are generically unstable under small perturbations [14][15][16][17][18] and will come to acquire nontrivial radial as well as azimuthal structures through excitation of the unstable modes.Thus the structures of galaxies will be elucidated by scrutinizing the stability of the rotating disk.In this section, we will examine the linear stability of the rotating axisymmetric fluid disk without recourse to the WKB approximation.

Perturbed Euler equation
Hereafter, for ρ, we simply write ρ, which will cause no confusion with the 3D density ρ(r).In the inviscid limit, the velocity field is given by v = (0, r , 0), (4.1) where ≡ v θ 0 /r with v θ 0 being the solution of the Navier-Stokes equation found in the previous section.We are interested in the stability of (4.1) under the linear perturbations which are governed by the Euler equation.We allow for the r -dependence of ω.
In the linear perturbation analysis, the Fourier expansions of perturbations of similar form are usually assumed.A common factor exp[i(mθ − ωt)] is involved in (4.2), too, but this does not mean that the expressions in (4.2) are Fourier expansion approximations in the usual sense.This is because i) ω may be a function of r , and ii) the amplitudes A, B, C, and D may generally be functions of t as well as of r .With these presumptions, we expect that the above expressions provide exact solutions within the linear perturbation analysis.Notice that, although the double Fourier expansion a(k, ω)e ikr−iωt dkdω of perturbations is assumed in the prevalent WKB method, the factor exp[−iω(r )t] with a non-constant function ω(r ) is generally not subjected to the above expansion.
Substituting (4.2) for (3.1a), (3.1b), and (3.1d) and keeping the terms linear in the perturbations, we have PTEP 2015, 073J01 K. Takahashi where the prime stands for a derivative in r .When ω is not zero, terms linear in t are present in (4.3a) and (4.3c).This is the reason that some of the amplitudes A, B, C, and D are also expected to have t-dependences.In fact, it is easy to see that Eqs.(4.3) have consistent solutions when only the amplitude D is linearly dependent on t as and A, B, C, D 0 , and D 1 are time-independent.Algebraic but cumbersome treatments of (4.3) and (4.4) are needed to our end.We here summarize the case of the flat rotation curve in (4.1) that gives a typical result.
i) The angular frequency ω obeys the equation where α in (4.5c) is a positive constant for Newtonian central gravity, to which case we restrict ourselves.Four EFs are thus exactly known in the sense that the algebraic equation (4.5a) can be exactly solved by Ferrari's method.See the appendix for a brief explanation of the derivation of (4.5a).ii) Four EFs exist, as are depicted in Fig. 4 for the flat rotation curve.Generally, two ω are real and the other two are complex.An exponentially growing unstable mode exists and drives the formation of the density wave [15,17,18].iii) ω depends on r .The winding of spirals of perturbations takes place simply because |Reω|t monotonically increases with time.iv) There exist points where the conditions of the corotation resonance, ω = m , and the Lindblad resonances, ω = ω p ≡ m + κ or ω = ω m ≡ m − κ, where κ ≡ (4 ( + r /2)) 1/2 is the epicycle frequency, are satisfied.Singularities associated with these resonances are inherent in the system of Eqs.(4.3).The EF equation (4.5a) does not suffer from these singularities and its roots depicted in Fig. 4 are all finite.v) Except for the mode denoted as I in Fig. 4(a), the stability nature discontinuously changes at points A and B, where bifurcation and coalescence of ω are taking place, respectively.Consequently, there exists a region of r that we shall call the polynomial-instability section (PIS) where all ω are real but D grows linearly in time.(This linearity will in turn cause the linearity of the perturbation of gravitational potential [9].) vi) The PIS shifts inward (outward) as the central gravity is weakened (strengthened).vii) The real part of ω in the outer region of the disk is positive and decreasing with r , thereby revealing a trailing profile for the pattern of perturbations.This fact holds not only for stable but also for unstable modes.

The cases of non-flat rotation curves
The rotation curves of galaxies are not necessarily flat.Let us examine the effect of deviation from the flatness specific to iv)-vii) above.As examples, we adopt two rotation curves depicted in Fig. 2, excluding the flat one.in units of 10 4 ly•(10 5 ms −1 ) 2 .The solutions to (4.3a) are shown in Fig. 5.One of the remarkable features observed in Fig. 5 is that the case of K 1 /K = 1.2 resembles the case of the flat rotation curve.That is, four Reω are dominantly positive and decreasing in the outer region of r and the PIS exists in the calculated region.This implies that the spatial patterns generated by the perturbations are trailing for t > 0.
The other remarkable feature is that, for K 1 /K = 0.8 or for the rotation curve decreasing at large distances, the PIS does not exist in the calculated region (it is in the further region from the origin, although not shown here).This means that the stability nature of each of the four modes is not altered in the whole region of r in Fig. 5.We also notice that there exists one neutral mode, whose ω is increasing in r > 1.When this mode is excited, the leading spiral pattern will emerge.

Wave patterns formed by the perturbations
In this section, the spatial pattern and the propagation of the perturbations are scrutinized for the case of the flat rotation curve (K 1 /K = 1 in Fig. 2).As for ω, we take the sequence II -IV-VI in Fig. 4, which has a negative imaginary part, and the system is stable, from which the behaviour for the unstable mode may be easily inferred.

Phase velocity
The r -and θ -components of the phase velocity of the nonzero mode are given respectively by ) (5.1b) v r > 0 for our choice of EF and the perturbation propagates radially outward with a speed that slows down as time elapses.The azimuthal component v θ has positive signs, i.e., the perturbation propagates in the same direction as v θ .The perturbation pattern is mostly inside the disk rotation.
The spatial pattern of the perturbation at fixed t is determined by where −tReω is the shape function of the spirals.Since Reω/m is positive and a decreasing function of r , the spatial patterns given by (5. in Fig. 6.The spirals also differentially rotate with the angular velocity Reω/m in the same direction as the disk rotation. Because of the factor t in the r.h.s. of (5.2), tight winding-ups of spirals result in as time elapses.In addition, discontinuous changes of curvature in the spiral arms are observed for the cases of α = 10 and 20.These are due to the PIS.For α = 2 (weak gravity), the PIS lies in the (high-density) region r < 1 and does not appear in the figure.Similarly, if α is rendered large enough (strong gravity), the PIS will be expelled to the (low-density) region r > 10 and will not appear in the figure .For the stable mode, the spiral pattern of matter clumps dies off in the decay time of the order 1/|Imω|.

Amplitudes
The relative magnitudes of the amplitudes in (4.2) are determined via (4.3) and (4.4) as ) where a has been defined by (4.5c).At the same time, a first-order differential equation of, e.g., A is obtained as

.4)
D 0 /A is undetermined, which is ascribed to the arbitrariness of the initial density perturbation.One natural choice is D 0 = 0, i.e., the variation in the density distribution takes place after the velocity field is varied.When the PIS exists, ω is not defined at its edges.This problem may be partially PTEP coped with by rewriting (5.4) as where the last term on the r.h.s. of (5.4) has been dropped.It is easier to solve (5.5) numerically.
The results for the flat rotation curve are shown in Fig. 7 for m = 2.All of the amplitudes have a peak around r = 3-5 and approach zero near the origin and at infinity.The snags of D 1 appearing in Figs.7(b) and (c) are due to the discontinuities of ω at the edges of the PIS and are concomitant to the strong gravity, as the insets in the figures reveal.These singularities in amplitudes will give rise to further instability due to nonlinearity that the linear perturbation fails to treat.
In (b) and (c), appreciable negative pressure gradients ∂ r p = aρv 2 θ /r must exist in the inner part of the disk (see (3.8)).These cases are not favourable for the following reason.Let ρ g , p g , and d g be the density, the pressure, and the thickness of a real characteristic spiral galaxy.Regarding our ρ as the projected density on the 3D galaxy, these are qualitatively related to our ρ and p by ρ ∼ ρ g d g and p ∼ p g d g .Dividing both sides of (3.8) by d g and differentiating with respect to r yields ∂ r p g ∼ aρ g v 2 θ /r .By substituting our galaxy's average density 10 −19 kg/m 3 for ρ g , we have |∂ r p g | ∼ 10 −9 Pa • (10 4 ly) −1 .One possible origin of the negative pressure gradient may be radiation from the central mass that we have assumed.Together with the graphs shown in the insets of Fig. 7, as an effective temperature at around r ∼ 10 4 ly we have ∼10 2 K, which seems too high compared to reality.We thus discard the models of α = 10 and 20.
In the case of α = 2 shown in Fig. 7(a), a(r ) > 0 (positive pressure gradient) and all of the amplitudes are regular in r > r 1 , where r 1 is the assumed inner radius of the disk.The regular amplitudes in this region of r are actually brought about for α < 2.12 and the PIS lies in the region r < r 1 (see Fig. 4).δv r , δv θ , and δp develop with the characteristic time τ c = 1/|Imω(r )|.Numerical calculations show that τ c takes values from 1.1 (i.e., 3.3 × 10 7 yr) at r ∼ 1.3 to 6.3 (1.9 × 10 8 yr) at r ∼ 10.On the other hand, the matter rotates with azimuthal velocities of about 1.5 (150 km/s) at these radii (Fig. 2, K 1 /K = 1), which correspond to 5.4 (i.e., 1.6 × 10 8 yr) to 42 (1.3 × 10 9 yr) of the rotation period.The perturbations develop about one order of magnitude faster than the fluid's circular motion.
The perturbed azimuthal velocity is given by where the plus and minus signs correspond to unstable and stable modes, respectively.A plot of v θ (r, θ, t) vs. r gives the rotation curve of the disk.It is dependent on θ and t.Some examples for the flat background rotation curve (K 1 /K = 1) of stable m = 2 mode are shown in Fig. 8.The rotation curves are almost flat as long as the perturbation is small.With the lapse of time, the velocity converges to the background velocity with smaller rates at larger distances.For the unstable mode, the pattern of the rotation curves temporarily varies from right to left in Fig. 8.The velocity dispersions in galaxies have been observed and are possibly ascribed to non-circular motions of stars [34].Azimuthal angle dependences of circular motions due to perturbations may be far more difficult to observe.

Summary and remarks
The Navier-Stokes equations together with the ν-expansion method reveal that several types of rotation curves are possible for a thin disk of a steady flow surrounding a mass.The fundamental profiles of the axisymmetric and steady velocity field, together with its relation to the density, are determined in an independent way to the external force (in our case, gravity).As ν becomes small, the radial velocity approaches zero, while its ratio to the viscosity is finite and maintains a close relation to the rotation velocity.The types of rotation include, as asymptotic radial behaviours, exponential increases (Fig. 1), constant, power-law increases, and power-law decreases (Fig. 2).Near the symmetry axis, the rotation is uniform.
In the proximity of the symmetry axis, namely r → 0, the density profile shows two possible patterns.One is the power law behaviour 1/r 2−δ and the other is a rapid convergence to zero.We showed that these behaviour are correlated to the form of v r (r) and v z (r).As for the halo of real galaxies, although a universality of density distribution has been proposed [35] based on the cold dark matter model, observations [36,37] and other simulative studies [35,38] have revealed that the dark matter density profile of galaxies also exhibits variety to some extent.Simon et al. [36] also reported the detection of radial motions in some galaxies.In view of our fluid model, the vertical component will be non-vanishing, too, if finiteness of viscosity makes some contribution to these motions.Further exploration of the correlation between the density profile and the distribution of velocity components may be intriguing.
The pressure p 0 turned out to peak around the central part of the disk.Recalling that the pressure is translated to the velocity dispersion v 2 in the Jeans equations, this profile of p 0 implies that v 2 is larger near the center of the stellar disk.If an older star's velocity dispersion tends to be large, as

Fig. 4 .
Fig. 4. Solutions of (4.5a) for the flat rotation curve with K 1 /K = 1 in Fig. 2 and m = 2. (a) Reω vs. r .(b) Imω vs. r .Parameter values are α = 10, K = 10 ly −2 .In (a), = v θ /r , ω p = m + κ, and ω m = m − κ are also plotted by dotted, dashed, and dot-dashed curves, respectively (for κ, see text).Six branches are drawn by solid curves and are labeled by Roman numerals I-VI.Branches I, III, and V are real, while II and IV are complex.The section between A and B in (b) is the PIS.

Fig. 6 .
Fig. 6.Temporal evolutions of the spatial phase pattern of m = 2 spirals for the model parameter α = 2, 10, and 20 as observed in an inertial frame.All are undergoing anticlockwise rotations.The radius range is 1 ≤ r ≤ 10.Snapshots at t = 1, 5, and 10 are shown.The fluid also flows anticlockwise.The unit of time is 3 × 10 7 yr.

Fig. 7 .
Fig. 7. Amplitudes of perturbation A, iB, iC, and D 1 in 1 ≤ r ≤ 10 for α = (a)2, (b)10, (c)20.Other parameter values are same as those in Fig. 4. Note the difference in the scales of the abscissa.Within a figure of a given α, only the relative values of four amplitudes have meanings.The insets are a(r ), whose negativity means that the gravity overwhelms the centrifugal force.