A perturbative analysis of Quasi-Radial density waves in galactic disks
aa r X i v : . [ a s t r o - ph . GA ] S e p Revista Mexicana de Astronom´ıa y Astrof´ısica , , ?? – ?? (2009) A PERTURBATIVE ANALYSIS OF QUASI-RADIAL DENSITY WAVES INGALACTIC DISKS.
X. Hernandez and I. Puerari Received October 31, 2018; accepted October 31, 2018
RESUMENABSTRACTThe theoretical understanding of density waves in disk galaxies starts fromthe classical WKB perturbative analysis of tight-winding perturbations, the keyassumption being that the potential due to the density wave is approximately ra-dial. The above has served as a valuable guide in aiding the understanding of bothsimulated and observed galaxies, in spite of a number of caveats being present. Theobserved spiral or bar patterns in real galaxies are frequently only marginally consis-tent with the tight-winding assumption, often in fact, outright inconsistent. Herewe derive a complementary formulation to the problem, by treating quasi-radialdensity waves under simplified assumptions in the linear regime. We assume thatthe potential due to the density wave is approximately tangential, and derive thecorresponding dispersion relation of the problem. We obtain an instability criterionfor the onset of quasi-radial density waves, which allows a clear understanding of theincreased stability of the higher order modes, which appear at progressively largerradii, as often seen in real galaxies. The theory naturally yields a range of patternspeeds for these arms which appears constrained by the condition Ω p < Ω ± κ/m .For the central regions of galaxies where solid body rotation curves might apply,we find weak bars in the oscillatory regime with various pattern speeds, includingcounter rotating ones, and a prediction for Ω p to increase towards the centre, asseen in the rapidly rotating bars within bars of some numerical simulations. Wecomplement this study with detailed numerical simulations of galactic disks andcareful Fourier analysis of the emergent perturbations, which support the theorypresented. Key Words:
INSTABILITIES — STELLAR DYNAMICS — GALAXIES:KINEMATICS AND DYNAMICS — GALAXIES: STRUC-TURE — GALAXIES: SPIRAL
1. INTRODUCTIONThe development of short-wavelength, tight-winding disk wave dynamics, notably by Lin & Shu(1964), Shu (1970) and Toomre (1969), still formsthe basic analytical substrate used to understandgalactic structure observations and numerical sim-ulations of galactic dynamics. This WKB perturba-tive analysis exploits the fact that for tight-windingspirals, the long range character of the gravitationalpotential due to the density wave one is modelingdisappears, the response becomes local, and an ana- Instituto de Astronom´ıa, Universidad NacionalAut´onoma de M´exico, M´exico ([email protected]). Instituto Nacional de Astrof´ısica, Optica y Electr´onica,Santa Mar´ıa Tonantzintla, M´exico ([email protected]). lytical formulation to the problem becomes feasible.The dispersion relation that results has been usedsuccessfully over the decades to gain some under-standing of the expected physical scalings betweenpattern speed, the orbital frequency, the epicycle fre-quency and the mass density of galactic disks.However, several weak points in the applicationof the classical WKB density wave dynamics to realgalaxies have been recognised since they where orig-inally proposed. First, one would like spiral armsto be very tightly wound, making the approximationon the radial dominance of the perturbed potentialclearly valid. For typical early type spirals, spiralarms are clearly tightly wound, but often one findslate type spirals with almost radial ’spiral’ arms.1 X. HERNANDEZ AND I. PUERARIThe spiral patterns of many galaxies appear to beon the limit of applicability for this approximation,and many, galactic bars included, are outright out ofit. One of the strong points of the density wave the-ory was the substantial easing of the winding prob-lem, which can actually be used to rule out any mate-rial arms interpretation for spiral structure in galax-ies, provided one seeks a long lived structure.Here we explore a complementary, and in a cer-tain sense, perpendicular approach. We shall assumethat the density wave being treated is quasi-radial,and that it represents only a small disturbance onthe overall potential. Under this assumptions, thepotential of the disturbance will only be felt nearthe density enhancement, and will be tangential, di-rected towards the density enhancement itself, par-ticularly if one assumes no strong radial gradients inthe amplitude of the density wave itself. Treatingonly the circular motions of a collisionless compo-nent, and ignoring the vertical gradients of all quan-tities involved, lands us in the thin disk approxima-tion for stellar disks, a scenario often used in galacticdynamical works and in studies of galactic densitywaves in analytical developments, orbital structureanalysis or N-body simulations, being Hernandez &Cervantes-Sodi (2006), Jalali & Hunter (2005), Ma-ciejewski & Athanassoula (2008) and O’Neil & Du-binski (2003) examples of the above.This allows an approximate analytical treatmentof the problem which we present here in section (2).A dispersion relation for the problem appears, theconsequences of which we explore under two ide-alised galactic regimes, flat rotation curve, and solidbody rotation, with the aim of sampling the resultingphysics under the approximate dynamical conditionswhere one finds spiral arms and bars. As in the caseof the Lin-Shu formalism, we begin by treating only astellar disk, thinking of modeling the old stellar pop-ulation of a present day galaxy. The development isidealised and remains within the linear regime, forquasi-radial density waves, for a disk made up ofclose to circular orbits. While lacking the detail ofmore advanced studies (e.g. Evans & Read 1998,Polyachenko 2004, Jalali & Hunter 2005), it has theadvantage of identifying clearly a physical instabilitycriterion which captures some of the details missedby the Lin-Shu formalism, highlighting the role ofstellar diffusion. As an independent confirmation ofthe ideas presented, we also run an extensive gridof numerical models, and analyse the details of thesmall amplitude density perturbations which appear.It is encouraging that the results of these direct nu- merical simulations validate the general results of theideas developed in this study.In section (3), for the case of the flat rotationcurve regime, we derive a dispersion relation forthe problem, and a corresponding stability criterionfor the onset of quasi-radial density perturbations,which is seen to be clearly satisfied for any reason-able galactic disk, for low values of m and Q , thesymmetry number for the pattern and Toomre’s pa-rameter. The case of a solid body rotation regimeis also treated, where we show that for weak barsin the oscillatory regime of the dispersion relation,relatively constant values for Ω p appear quite natu-rally, at values which represent ’slow’, ’fast’ or evencounter-rotating bars, for natural values of m and Q . The extent of these bars being naturally limitedby the corotation condition. A detailed numericalmodel is presented in section (4), including a carefulFourier analysis of the density perturbations whichappear, yielding results in consistency with the the-ory developed. Section (5) presents our conclusions.2. PHYSICAL SETUPIn this section we develop the simplified physicalmodel for quasi-radial density waves. As opposedto the Lin-Shu formalism, designed to model tight-winding waves, we assume that the derivatives of thepotential due to the perturbation are dominated bythe angular component, rather than the radial one.We shall start from the angular moments equationwritten in cylindrical coordinates ( R, θ, z ), in a ref-erence frame which rotates with a constant angularfrequency Ω P parallel to the z axis, but without as-suming azimuthal symmetry. In spite of the multi-component nature of a galactic system, we shall be-gin by ignoring the gaseous component, in an at-tempt to model the underlying density perturba-tion represented by loosely-wound arms in the oldstellar population. We shall also assume no ver-tical or radial velocities, in modelling weak radialarms which only result in weak angular distortionsin both the density and velocity fields. Taking alsoan isotropic velocity ellipsoid with no vertex devia-tion, i.e. σ ij = δ ij σ , results in the following equation(e.g. Vorobyov & Theis 2006): ∂ ( R Σ v θ ) ∂t + 1 R ∂ ( R Σ v θ ) ∂θ + ∂ (Σ σ ) ∂θ = − Σ ∂ Φ eff ∂θ (1)In the above equation Φ ef = Φ − | ~ Ω p × ~r | / + ǫ Σ ( θ ) v θ = v ( R ) + ǫv ( θ ) σ = σ + ǫσ ( θ )Φ( R, θ ) = Φ ( R ) + ǫ Φ ( R, θ )In the above Φ is the axially symmetric gravita-tional potential, which cancels the centrifugal forcefixing ~V ( R ) = (0 , R (Ω − Ω p ) , ( R ) the cen-trifugal equilibrium orbital frequency. Φ will hencebe the potential perturbation due to the quasi-radialarms, resulting in angular perturbations Σ ( θ ), v ( θ )and σ ( θ. ) The form taken for the velocity impliesno further inertial terms are necessary. The unper-turbed state trivially satisfies eq.(1) for the standardaxisymmetric centrifugal equilibrium state.To first order in the perturbation, eq.(1) now be-comes: R Σ ∂v ∂t + R ∆Ω ∂ Σ ∂t + 2 R ∆ΩΣ ∂v ∂θ +( R ∆ Ω + σ ) ∂ Σ ∂θ + 2 σ Σ ∂σ ∂θ = − Σ ∂ Φ ∂θ (2)In the above equation we have introduced ∆Ω =(Ω − Ω p ), and have omitted making explicit the ra-dial dependence of all the above variables, which re-mains implicit.We now introduce an isothermal equation of state(e.g. Binney & Tremaine 1987) or equivalently theconservation of phase-space density Σ /σ = cte. toeliminate σ for Σ through: σ = (cid:18) σ (cid:19) Σ , (3)with which eq.(2) becomes: R Σ ∂v ∂t + R ∆Ω ∂ Σ ∂t + 2 R ∆ΩΣ ∂v ∂θ + (cid:0) R ∆ Ω + σ (cid:1) ∂ Σ ∂θ = − Σ ∂ Φ ∂θ (4)We see the assumption of eq.(3) in the numericalconstant of the fifth term in the above equation. Anyother similar equation of state would only changethis constant slightly, leaving all conclusions quali-tatively unchanged, and only modified to within asmall numerical factor of order unity. The partialangular derivative of the above equation reads: R Σ ∂ v ∂θ∂t + R ∆Ω ∂ Σ ∂θ∂t + 2 R ∆ΩΣ ∂ v ∂θ + (cid:0) R ∆ Ω + σ (cid:1) ∂ Σ ∂θ = − Σ ∂ Φ ∂θ (5)We will change the dependencies of the aboveequation on v for dependencies on Σ through theuse of the continuity equation, which reads: ∂ Σ ∂t + 1 R ∂ ( R Σ v R ) ∂R + 1 R ∂ (Σ v θ ) ∂θ = 0 . (6)Introducing the perturbation conditions withonly tangential velocities into the above continuityequation yields: ∂ Σ ∂t + Σ R ∂v ∂θ + ∆Ω ∂ Σ ∂θ = 0 . (7)From this last equation we can solve for ∂v /∂θ to construct: ∂ v ∂θ = − R ∆ΩΣ ∂ Σ ∂θ − R Σ ∂ Σ ∂θ∂t (8) ∂ v ∂θ∂t = − R ∆ΩΣ ∂ Σ ∂θ∂t − R Σ ∂ Σ ∂t (9)The last two above are now used to replace thedependence on v appearing in eq.(5) for a depen-dence on Σ , giving: − R ∂ Σ ∂t − R ∆Ω ∂ Σ ∂θ∂t + (cid:0) σ − R ∆ Ω (cid:1) ∂ Σ ∂θ = − Σ ∂ Φ ∂θ (10)Finally, we relate Σ and Φ through the firstorder Poisson equation of the problem, which underthe assumption of tangential forces due to the almostradial perturbation dominating over the other twodirections gives: 1 R ∂ Φ ∂θ = 4 πGρ , (11)allows to write eq.(10) as: − R ∂ Σ ∂t − R ∆Ω ∂ Σ ∂θ∂t + (cid:0) σ − R ∆ Ω (cid:1) ∂ Σ ∂θ = − πGρ Σ (12)In the above equation we have used h , the heightscale of the disk, as h = Σ /ρ = Σ /ρ .We now turn to well known galactic structurescalings to re-write the right hand side of eq.(12) interms of more helpful dynamical variables, specifi-cally, we shall make use of vertical virial equilibriumin the disk and Toomre’s stability criterion, X. HERNANDEZ AND I. PUERARI κσ πG Σ = Q. (13)In the above κ = 4Ω + 2Ω R ( d Ω /dR ) is theepicycle frequency as functions of the radius. Thestar formation processes in the disks of spiral galax-ies have often been thought of as constituting a selfregulated cycle. The gas in regions where Q < a cer-tain critical value finds itself in a regime where lo-cal self-gravity, appearing through ( G Σ ) in eq.(13),dominates over the combined effects of rotationalshears (or equivalently global tidal forces) and thestabilising thermal ’pressure’ effects of the product( κσ ). This leads to collapse and the triggering of starformation processes, which in turn result in signifi-cantly energetic events. The above include radiativeheating, the propagation of ionisation fronts, shockwaves and in general an efficient turbulent heatingof the gas media, raising σ locally to values result-ing in Q > a certain critical threshold. This restoresthe gravitational stability of the disk and shuts offthe star formation processes. On timescales longerthan the few ×
10 million years of massive stellar life-times, an equilibrium is expected where star forma-tion proceeds at a rate equal to that of gas turbulentdissipation, at time averaged values of Q ∼ Q crit .Examples of the above can be found in Dopita &Ryder (1994), Koeppen et al. (1995), Firmani et al.(1996) and Silk (2001).The preceding argument applies to the gaseouscomponent, not the stellar one, which is the one weare interested in modeling, particularly the old stel-lar component. However, if stars retain essentiallythe velocity dispersion values of the gas from whichthey formed, one also expects Q ∼ Q crit for the stars.Alternatively, if we consider the dynamical heatingof the stellar populations in the disk through encoun-ters with giant molecular clouds or the spiral armsthemselves, slightly larger values of Q would be ex-pected for the stars than for the gas. We shall there-fore take eq.(13) as valid throughout the disk (e.g.O’Neill & Dubinski 2003, Maciejewski & Athanas-soula 2008), without specifying any particular val-ues of Q at this point. Imposing viral equilibriumin the vertical direction in the disk (e.g. Binney &Tremaine 1987) yields: h = σ πG Σ , (14)Since h = Σ /ρ , we can replace the dependenceon Σ for one on ρ in eqs.(13) and (14), and solvefor ( σ /h ) in both eqs.(13) and (14) to obtain: πGρ = 2 (cid:18) κQ (cid:19) . (15)Equation (15) serves to re-write the right handside of eq.(12) and obtain: − R ∂ Σ ∂t − R ∆Ω ∂ Σ ∂θ∂t + (cid:0) σ − R ∆ Ω (cid:1) ∂ Σ ∂θ = − (cid:16) κQ (cid:17) (16)Equation (16) is now a second order partial dif-ferential equation for the temporal and angular vari-ations of the density of the perturbation. This isvalid at each radius, and shows dependence on boththe gravitational dynamics of the system through κ and Ω, and on the local stability and structure ofthe disk through Q and σ . Now we impose periodicsolutions in θ for the perturbation in density, withan oscillatory or growing time dependence:Σ ∝ e − i ( mθ − ωt ) , (17)with m an integer, which when introduced intoeq.(16) gives:( ω − m ∆Ω) = 53 (cid:16) mσ R (cid:17) − (cid:18) κQ (cid:19) . (18)Equation (18) is the dispersion relation for theproblem. For a density wave of the type being de-scribed to develop, we require for ω to be imaginary,i.e., 2 κQ > (cid:18) (cid:19) / mσ R (19)This last is an instability criterion which gives thevalues of σ below which a certain m-symmetry pat-tern will begin to grow, for a given rotation curve andvalue of Q. We can also obtain directly the ranges ofvalues Ω p is expected to take, considering that thislimit values will be given by the condition ω = 0 inequation (18). Since the pattern will develop onlywhen the second term in the R.H.S. of eq.(18) domi-nates over the first, we can start by ignoring the firstterm on the R.H.S. of this equation. For example, ina flat rotation curve region, with a circular velocityof 220 km/s and a σ = 30 km/s , even with a Q = 3,the first term above is only 1.7 % of the second for m = 1, one has to go to m = 4 for this first termto pass 20 % of the second. We can hence estimatethe limit values of Ω p by neglecting this first term togive:UASI-RADIAL DENSITY WAVES IN GALAXIES 5 m ∆Ω p,lim = ± √ κQ . (20)Remembering that ∆Ω = Ω − Ω p , we obtain:Ω p,lim = Ω ± √ Q ! κm . (21)Since Q is of the order of 2 √
2, the above equa-tion is a direct analytical explanation for the rangesof pattern speeds seen in numerical experiments, in-dependent of the swing amplification hypothesis.Going back to equation(19), we see that the lowervalues of m are easier to excite than the higher ones,as often found in numerical simulations, providing anexplanation for observed galaxies with strong, highpitch angle arms, being preferentially m = 2 sys-tems. Indeed, Martos et al. (2004) showed that the m = 4 optical spiral pattern in the Milky Way canform as a consistent hydrodynamical response of thegas to an underlying m = 2 pattern in the old stel-lar population, as traced by COBE-DIRBE K-bandphotometry.Remembering that κ/Q is really just a place-keeper for ( Gρ ) / , the instability criterion ofeq.(19) in the original terms reads: Rmσ > (cid:18) π (cid:19) / (cid:18) Gρ (cid:19) / (22)In this terms, the instability criterion of eq.(22)is just a comparison between the local gravitationaltimescale of the disk, the free-fall timescale, and theinter-arm diffusive crossing time. That is, the insta-bility will not develop if diffusion is such that theperturbation is blurred out before it can condensegravitationally. The lower the number of arms, theeasier it is to excite the pattern, as in high m pat-terns, it is easier for stars to diffuse from one arm toanother, blurring and erasing the proposed pattern.To conclude, we identify the dimensionless number P m = (cid:18) (cid:19) / mQσκR = (cid:18) πGρ (cid:19) / mσR (23)as the critical indicator for an n-armed quasi ra-dial density wave to develop, with the instabilitythreshold appearing at P < m and Q . 3. CONSEQUENCES FOR GALACTIC DISKSIn this section we shall trace some of the con-sequences for the onset of the instability presented,in the context of thin galactic discs. For simplicity,in this first treatment, we shall consider only twodistinct regimes for the rotation curve of a galaxy,a strictly flat rotation curve regime, of relevance tothe development of galactic arms, and a strictly solidbody rotation curve region, appropriate for the mod-eling of dynamics of bars. The above with the inten-tion of developing as clear a physical understandingof the effect being introduced as possible. A moredetailed treatment of the problem using more real-istic rotation curves and detailed numerical galacticmodels, appears in section (4).3.1 . Galactic Arms in Flat Rotation Curves We start this section with a consideration of theflat rotation curve limit of the relations derived pre-viously, as a suitable first approximation to the dy-namics of a galaxy in the region over which galac-tic arms are seen. In a flat rotation curve regime, κ = 2 / Ω, and therefore, eq.(19) becomes: V r σ > (cid:18) (cid:19) / nQ . (24)In real galactic systems, typically with 5 8, we see that in a flat rotation curve regimethe instability criterion will always be satisfied for m = 2, for example, in a Q = 1 or Q = 2 disk. For m = 4 it is still easy to satisfy the instability crite-rion for the onset of a strong perturbation, for Q < m would be hard to accommo-date. Indeed, in the detailed stability analysis ofsimple power-law disks in idealised rotation curvesof Evans & Read (1998), it is shown explicitly thatlower m modes are always the most unstable ones.Actually, from the radial dependence of P m , it can beseen that higher m patterns will be easier to exciteas one moves further out in a galactic disk, in con-sistency with the bifurcations often seen in galacticspiral arms, patters which tend to be characterisedby higher order symmetries as one moves to largerradii.It is interesting to note that since Q scales lin-early with σ , P m which scales with Qσ , is a quadraticfunction of σ . This last means that for small valuesof σ , P m < Q , with the opposite being true for largevalues of σ . By setting P m = Q we obtain4 V r mσ = (cid:18) (cid:19) / , (25) X. HERNANDEZ AND I. PUERARIfor low values of m this condition will result in avalue of σ larger than that present in normal galac-tic disks, and hence we see that typically galaxieswill lie in the region P m < Q . This is interesting,as perhaps the identification of values of Q ∼ Q stable, P unstable disks, with quasi-radial densitywaves arising, and then being wound up into tightspirals.Remembering again that ( κ/Q ) is a proxy for Gρ , we can compare equation (18) to the classi-cal WKB equivalent relation, e.g. eq.(6.61) Binney& Tremaine (2008). We see that the term includ-ing σ in equation (18) appears in replacement ofthe standard κ term. By considering a purely radialperturbation with a purely tangential force field andresulting matter flows, we have lost the contribu-tion of differential rotation, in what is essentially ananalysis at constant radius, in favour of a dynamicalpressure term through the stellar velocity dispersion,which does not appear in the classical criterion.3.2 . Central Bars in Solid Body Rotation Curves In moving to central galactic bars we now go toa solid body rotation regime. The above considera-tion leads to the use of m = 2 in what follows. Inthis regime, κ = 2Ω , with Ω a constant, and theinstability criterion of eq(19) is modified by the sub-stitution of a 6 for the 3 in the denominator of thesquare root term in the right. In terms of Ω thisreads: R Ω σ > (cid:18) (cid:19) / Q p forthese weak pulsating bars can be found from eq.(18)setting ω = 0, which under solid body rotation gives:Ω p Ω = 1 ± " (cid:16) σ R Ω (cid:17) − (cid:18) Q (cid:19) / . (27)We find a slow and a fast solution, as a func-tion of the radial coordinate, Ω , σ and Q . Fig-ure (1) now gives values of Ω p / Ω , as a function ofthe dimensionless radial coordinate υ = R Ω /σ , for3 values of Q , 1, 2 and 3, inner, middle and outer -20-15-10-5 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Ω p / Ω υ Fig. 1. Values of Ω p / Ω for stationary weak bars, asa function of the dimensionless radial coordinate υ = R Ω /σ , for Q=1, inner curves, Q=2, middle curves andQ=3, outer curves. The vertex of the curves defines theextent of the weak bars in question, always at the coro-tation radius. curves, respectively. Notice that for large values of υ , relatively constant values of Ω p appear, for val-ues of Q > 2. In going towards the centre however,the divergence towards R = 0 becomes evident, andlarge variations in Ω p appear. This last feature mightexplain the appearance of bars within bars, whichexhibit increasingly rapid pattern rotation rates ingoing towards the central regions (e.g. Heller et al.2007, Maciejewski & Athanassoula 2008 and Shen& Debattista 2008). We see that the lower branchof the curves, corresponding to the minus sign ineq.(28), can have negative values over much of itsextent. This corresponds to counter-rotating bars,where the stars are not counter rotating themselves,but the density enhancement is.Also, we see that replacing the inequality for anequality in eq.(27) we can solve for R M , the max-imum extent of a weak bar of the type being de-scribed. This last is given by the vertex of thecurves in figure(1), and as can be seen from eq.(27),will scale linearly with Q . From eq.(28) we see thatthese weak bars end at corotation, as already clearfrom eq.(21), the transition between the oscillatoryand the unstable regimes coincides with the coro-tation condition. The interpretation provided hereto eq.(22), highlighting the role of orbital diffusion,might offer a simple explanation to the commonlyfound truncation of bars at corotation in both or-bital structure studies and N-body simulations, e.g.the unstable modes of Jalali & Hunter (2005), thesimulated bars of O’Neill & Dubinski (2003) or theUASI-RADIAL DENSITY WAVES IN GALAXIES 7 Fig. 2. Disk particles distribution at T = 0. The upperpanel is 12 × LU , and the other ones are 12 × LU . outer of the double bars of Maciejewski & Athanas-soula (2008) and Shen & Debattista (2008).In summary, weak bars can develop in the os-cillatory regime of eq.(18), and in the case of solidbody rotation curves, present an almost constantnon-winding Ω p for a large range of radii, for val-ues of Q > 2, perhaps related to the ’pulsating bars’seen in some numerical experiments, e.g. Shen &Debattista (2007).4. NUMERICAL MODELS AND SIMULATIONSIn order to check the predictions of the modelpresented in the more realistic cases of fully self consistent 3D galactic models, we have run a largegrid of simulations of isolated, stable disk mod-els. The models contain three self-consistent compo-nents, corresponding to a disk, a bulge and a halo.The construction of the models and their followingevolution were performed using the package Nemo(http://carma.astro.umd.edu/nemo/).Setting up a stable three-dimensional disk forN-body experiments is a difficult task. Sev-eral strategies have been suggested including grow-ing the disk mass distribution in a self-consistenthalo/bulge model (Barnes 1988, Athanassoula, Puer-ari & Bosma 1997), or treating the spherical compo-nents halo+bulge as a static background (Sellwood& Merritt 1994, Quinn, Hernquist & Fullagar 1993,Sellwood 2011). Some authors have also directlysolved the Jeans equation for the complete disk,bulge and halo system to find the velocity disper-sions (Hernquist 1993). The more important prob-lem in the construction of self-consistent 3 compo-nent models is the fragility of a cool disk, generatinggrand design spiral structure and bars, as well as aseries of other known instabilities (bending modes,disk warping, etc).Our models were constructed using the Nemoroutine mkkd95 . The routine is based on the methoddescribed in Kuijken & Dubinski (1995). Their strat-egy is the choice of analytic forms for the distribu-tion functions (DF) of the three components. Forthe bulge, the DF is taking as a King model (King1996). For the halo, the DF is a truncated Evans(1993) model for a flattened logarithmic potential.For the disk, the more difficult component to buildup in the model, the authors used a vertical exten-sion of the planar DF discussed by Shu (1969) andKuijken & Tremaine (1992). All the DF’s are usedto calculate the spatial density of each galaxy com-ponent, and the Poisson equation ∇ Ψ( R, z ) = 4 πG [ ρ disk ( R, θ, z )+ ρ bulge ( R )+ ρ halo ( R )](28)is solved by using spherical harmonic expansionby Prendergast & Tomer (1970) with some modifi-cations (see details in Kuijken & Dubinski 1995).For the evolution of the models, we used theNemo routine gyrfalcON , which is a full-fledged N-body code using a force algorithm of complexity O ( N ) (Dehen 2000, 2002), which is about 10 timesfaster than an optimally coded tree code. In theprogram, the gravitational constant is taken equalto 1. To match this value, an appropriate normal-ization has been chosen with length unit LU = 3kpc, and time unit T U = 10 years. Using these val- X. HERNANDEZ AND I. PUERARI Fig. 3. Upper panel: tangential velocity curves at T = 0. The thick solid line is for the model with σ R, = 0 . 27. The solid, dashed and dotted lines are for σ R, = 0 . , . , and 0.57, respectively. Bottom panel:Toomre’s Q . The lines are as in the upper panel, i.e.,the thick solid line represents the cooler model, and thedotted line represents the warmer one. ues, the units of mass, velocity and volume densityare 6 × M ⊙ , 293 km s − , and 2.22 M ⊙ pc − , re-spectively, making the unit of length 3 kpc. In thissection, we plot all quantities in computational units.With gyrfalcON , all the simulations were run usinga tolerance parameter θ = 0 . 5, a softening parameter ǫ = 0 . 05 and a time step of 1 / . These values ensurea total energy and angular momentum conservationof the order of 10 − and 10 − , respectively.Our models were chosen to be Milky WayA (MW-A) like-models from Kuijken & Dubinski(1995), having masses of M D : M B : M H = 0 . 82 :0 . 42 : 5 . N D : N B : N H =160 , 000 : 40 , 000 : 400 , Fig. 4. The P number from eq.(23) as a function ofradius for 4 of our 16 simulations. Upper panels show themodels with σ R, = 0 . 37 ( z = 0 . 05, left; z = 0 . 1, right),and bottom panels show colder models with σ R, = 0 . z = 0 . 1, left; z = 0 . 15, right. The plotted curves arethe mean P from T = 0 to 200, 2Gyr. The lower curvein each panel is for m = 2 and the upper one for m = 7.The horizontal lines show limits for P between 1.0 and1.15. The upper and lower cuts of these limits and each m curve constrain the radial region outwards of whichwe expect the power of that given m to arise. time when using the fully self-consistent model with300,000 particles; only the disk scale length changesof the order of 15%. All our models start with novertex deviation, and result in average values close to σ θθ /σ RR = 0 . 75 and σ zz /σ RR = 0 . σ R, = 0 . 47 and a vertical scale length z = 0 . σ R, =0 . , . , . , . 57 and z = 0 . , . , . , . 2, andtherefore covering the space parameter from coolthin disks, to warmer/thicker ones. Here we presentresults for colder/thinner models; in warmer/thickermodels, the structures present much lower ampli-tude, and they are much more chaotic. Our fully self-consistent colder/thinner models support low ampli-fication of several modes, unlike the unstable models,UASI-RADIAL DENSITY WAVES IN GALAXIES 9 Fig. 5. Power of the m components as a function ofradius for the same simulations in figure (4). Differentlines are the means for different time intervals: solid line,0 < T < 50, dashed 50 < T < 37 models (upper panels) and 2 forthe colder model with σ R, = 0 . 27 (bottom panels). For each m , the P profile scales with σ ( R ) ρ ( R ) − / ,where σ ( R ) is the disk radial velocity dispersion pro-file and ρ ( R ), the disk volume density profile, c.f.eq.(23). The plotted curves are the azimuthally av-eraged time means of P as a function of radius forthe total time interval over which the models wereevolved, from T = 0 to T = 200 = 2 Gyr . Wehave draw two horizontal limits for P , P = 1 (lowerstraight line), and P = 1 . 15 (upper straight line).These two limits are used to estimate the radial limitoutwards of which a given m component is expectedto show a clear increase in its amplitude, as predictedby the developments of the previous sections.Now, to calculate the amplitude (or power) ofeach m component as a function of radius, we di-vide the disk in 50 rings (from R = 0 to R = 5,∆ R = 0 . F ( m ) = 1 N N X i =1 e − imθ i (29)where N is the number of disk particles in thering, m is the m -armed component and θ the polarcoordinate of each particle in the ring. The real andimaginary parts for each m component are I c ( m ) = N X i =1 cos ( mθ i ) I s ( m ) = N X i =1 sin ( mθ i )(30)and the power is then P ow m = I c ( m ) + I s ( m ) (31)The phase Φ m of each component is then:Φ m = atan ( I s ( m ) /I c ( m )) /m. (32)We have calculated the power and phase of thefirst ten modes, as a function of radius and time. Infigure (5), we show the power of the m = 2 to 7modes, as a function of radius. We present the meanvalue of the power for 4 time intervals. Superim-posed on the curves, we have draw shaded areas rep-resenting the radial interval resulting from the ana-lytic developments of section (2), as inferred througheq.(23) applied to figure (4). There is a clear agree-ment between the radial positions from figure (4) andthe radius at which we see an increase of power forthe different m ’s, lending empirical support to thetheoretical model presented in this study.Next, we present some plots of phase as a func-tion of radius and time in figure (6). We present0 X. HERNANDEZ AND I. PUERARI Fig. 6. The phase for the m = 3 (upper panels) and m = 6 (bottom panels) modes for two models appearing in Fig.7, as a function of time and radius (mkd s27 z10, left; mkd s37 z05, right). The phases range from 0 (black) to 2 π/m (white). As all features are almost vertical, we see that the perturbations appearing are quite radial, far from tightlywound. The horizontal line for each m of each model is the mean of the radial range show in figure (7). Note the clearchange in the behavior of the phase on crossing the radial position marked, where we expect an increase of power interms of the theory presented. the phase of m = 3 (upper panels) and m = 6 (bot-tom panels) for the models mkd s27 z10 (left) andmkd s37 z05 (right). In these graphs, the phasesrange from 0 to 2 π/m . The amplitudes of the per-turbations in our isolated models are quite low com-pared to those associated with bars and strong m = 2spiral arms in unstable models (Puerari, in prepara-tion). Even with these low amplitudes, the phaseshows a clear behaviour, quite ordered, represent-ing quasi radial perturbations on the disk, as evi-dent from the almost vertical form of the featuresappearing. We also mark in figure (6) the mean ra-dial position of the corresponding P < p , R space where power appears forthe different modes in one of the simulations run.We see again a clear agreement with the theoreti-cal expectations developed, with the amplitude ofthe modes constrained to the region where Ω p < Ω ± κ/m . Comparing with the analytical predic-UASI-RADIAL DENSITY WAVES IN GALAXIES 11 Fig. 7. The contour plots give the zones where power ispresent for the first 6 modes for one of the models, as afunction of frequency and radius. The three curves giveΩ and Ω ± κ/m tions of eq. (21), and given the values of Q of thesimulated galaxies of ∼ 2, we see again good a agree-ment. If one wanted to calibrate the details of thetheory through a comparison of eq.(21) with the re-sults of figure(7), we see than a small numerical ad-justment, e.g. in the first order scale height estimateof eq.(14) of a factor of 2, would furnish an accurateaccordance of the limit Ω p predictions of eq.(21) andfigure (7).In summary, we see that even though in the the- oretical developments of sections (2) and (3) we ig-nored the vertical structure of the galactic disk, theradial motion of the stars, and the radical variationsin surface density, epicycle frequency and velocitydispersion, the predictions for the critical radii atwhich quasi-radial perturbations appear, the way inwhich these increase progressively for higher modes,and the limit frequencies which these exhibit, areall in good accordance with the results of fully self-consistent 3D numerical models.5. CONCLUSIONSWe have developed a simple description of aquasi-radial density wave in a galactic disk, assum-ing only tangential forces due to the perturbation,we reach a description of the problem which is com-plementary to the classical tight-winding approxima-tion. The main features are the substitution of thedifferential rotation for the velocity dispersion of thedisk in the resulting dispersion relation. The result-ing dispersion relation allows a clear understandingof why lower m modes are more unstable than higher m modes, and of why the higher m values appearprogressively at larger radii.Weak bars in harmonic potential central regionsappear in the oscillatory regime, an expression for Ω p is found which shows a divergence towards R = 0,perhaps helping to understand the phenomenon ofincreasingly rapidly rotating bars within bars seenin numerical simulations. The extent of this weakbars is seen to be naturally limited by the corotationcondition.The expectations for quasi-radial spiral arms aretested in detail through extensive numerical simu-lations, which confirm to good accuracy the mainresults of the theory developed here.acknowledgementsThis work was supported in part through CONA-CyT (45845E), and DGAPA-UNAM (PAPIIT IN-114107 and IN103011) grants.REFERENCES Athanassoula E., Puerari I., Bosma A., 1997, MNRAS286, 284Barnes J.E., 1988, ApJ 331, 699Binney J., Tremaine, S., 1987, Galactic Dynamics(Princeton University Press, Princeton, NJ)Binney J., Tremaine, S., 2008, Galactic Dynamics(Princeton University Press, Princeton, NJ)Dehnen W., 2000, ApJ 536, L39Dehnen, W., 2002, JCP 179, 27Dopita M. A., Ryder S. D., 1994, ApJ, 430, 163Evans N.W., 1993, MNRAS 260, 1912 X. HERNANDEZ AND I. PUERARI