Stability of general relativistic Miyamoto-Nagai galaxies

The stability of a recently proposed general relativistic model of galaxies is studied in some detail. This model is a general relativistic version of the well known Miyamoto-Nagai model that represents well a thick galactic disk. The stability of the disk is investigated under a general first order perturbation keeping the spacetime metric frozen (no gravitational radiation is taken into account). We find that the stability is associated with the thickness of the disk. We have that flat galaxies have more not-stable modes than the thick ones i.e., flat galaxies have a tendency to form more complex structures like rings, bars and spiral arms.


I N T RO D U C T I O N
The natural shape of an isolated self-gravitating fluid is axially symmetric. For this reason, exact axially symmetric solutions of Einstein field equations are good candidates to model astrophysical bodies in general relativity. In the last decades, several exact solutions were studied as possible galactic models. Static thin disc solutions were first studied by Bonnor & Sackfield (1968) and Morgan & Morgan (1969), where they considered discs without radial pressure. Discs with radial pressure and with radial tension were considered by Morgan & Morgan (1970) and González & Letelier (1999), respectively. Self-similar static discs were studied by Lynden-Bell & Pineault (1978) and Lemos (1989). Moreover, solutions that involve superpositions of black holes with static discs were analysed by Lemos & Letelier (1993, 1994, 1996 and Klein (1997). Also, relativistic counter-rotating thin discs as sources of the Kerr type metrics were found by Bičák & Ledvinka (1993). Counter-rotating models with radial pressure and dust discs without radial pressure were studied by González & Espitia (2003) and García & González (2004), respectively; while rotating discs with heat flow were studied by González & Letelier (2000). Furthermore, static thin discs as sources of known vacuum space-times from the Chazy-Curzon metric (Chazy 1924;Curzon 1924) and Zipoy-Voorhees metric (Zipoy 1966;Voorhees 1970) were obtained by Bičák, Lynden-Bell & Katz (1993). Also, Bičák, Lynden-Bell & Pichon (1993) found an infinite number of new relativistic static solutions that correspond to the classical galactic disc potentials of Kuzmin & Toomre (Kuzmin 1956;Toomre 1963) and Mestel & Kalnajs (Mestel 1963;Kalnajs 1972). Stationary disc models including electric fields (Ledvinka, Zofka & Bičák 1999), magnetic fields (Letelier 1999) and both electric and magnetic fields (Katz, Bičák & Lynden-E-mail: mujevic@ufabc.edu.br (MU); letelier@ime.unicamp.br (PSL) Bell 1999) were studied. In recent years, exact solutions for thin discs made with single and composite haloes of matter (Vogt & Letelier 2003), charged dust (Vogt & Letelier 2004a) and charged perfect fluid (Vogt & Letelier 2004b) were obtained. For a survey on relativistic gravitating discs, see Semerák (2002) and Karas, Huré & Semerák (2004). Most of the models constructed above were found using the metric to calculate its energy momentum tensor, that is, an inverse problem. Several exact disc solutions were found using the direct method that consists of computing the metric for a given energy momentum tensor representing the disc (Neugebauer & Meinel 1995;Klein & Richter 1999;Frauendiener & Klein 2001;Klein 2001Klein , 2002Klein , 2003a. In a first approximation, the galaxies can be thought to be thin, which usually simplifies the analysis and provides very useful information. But, in order to model real physical galaxies, the thickness of the discs must be considered. Exact axially symmetric relativistic thick discs in different coordinate systems were studied by González & Letelier (2004). Also, different thick discs were obtained from the Schwarzschild metric in different coordinate systems with the 'displace, cut, fill and reflect' method (Vogt & Letelier 2005a).
The applicability of these disc models to any structure found in nature lies in its stability. The study of the stability, analytically or numerically, is vital to the acceptance of a particular model. Also, the study of different types of perturbations, when applied to these models, might give an insight on the formation of bars, rings or different stellar patterns. Moreover, a perturbation can cause the collapse of a stable object with the later appearance of a different kind of structure. An analytical treatment of the stability of discs in Newtonian theory can be found in Binney & Tremaine (1987), Fridman & Polyachenko (1984), and references therein. In general, the stability of discs in general relativity is evaluated in two ways. One way is to study the stability of the particle orbits along geodesics. This kind of study was made by Letelier (2003) Lifshitz 1987) into a general relativistic formulation. Using this criterion, the stability of orbits around black holes surrounded by discs, rings and multipolar fields was analysed (Letelier 2003). Also, this criterion was employed by Vogt & Letelier (2003) to study the stability of the isotropic Schwarzschild thin disc, and thin discs of single and composite haloes. The stability of circular orbits in stationary axisymmetric space-times was studied by Bardeen (1970) and Abramowicz & Prasanna (1990). Moreover, the stability of circular orbits of the Lemos-Letelier solution (Lemos & Letelier 1994) for the superposition of a black hole and a flat ring was considered by Semerák &Žáček (2000a,b) and Semerák (2003). Also, Bičák et al. (1993) analysed the stability of several thin discs without radial pressure or tension, studying their velocity curves and specific angular momentum. Another way of studying the stability of discs is perturbing the energy momentum tensor. This way is more complete than the analysis of particle motions along geodesics, because we are taking into account the collective behaviour of the particles. However, there are few studies in the literature performing this kind of perturbation. A general stability study of a relativistic fluid, with both bulk and dynamical viscosity, was done by Seguin (1975). He considered the coefficients of the perturbed variables as constants, that is, local perturbations. Usually, this condition is too restrictive. Stability analyses of thin discs from the Schwarzschild metric, the Chazy-Curzon metric and the Zipoy-Voorhees metric, perturbing their energy momentum tensor with a general first-order perturbation, were made by Ujevic & Letelier (2004), finding that the thin discs without radial pressure are not stable. Moreover, stability analyses of the static isotropic Schwarzschild thick disc as well as the general perturbation equations for thick discs were studied by Ujevic & Letelier (2007).
In Newtonian gravity, models for globular clusters and spherical galaxies were developed by Plummer (1911) and King (1966). In the case of disc galaxies, important thick disc models were obtained by Miyamoto & Nagai (1975) and Nagai & Miyamoto (1976) from the prior work of Kuzmin (1956) and Toomre (1963) on thin disc galaxies. Miyamoto and Nagai 'thickened-up' Toomre's series of disc models, and obtained pairs of three-dimensional potential and density functions. Also, Satoh (1980) obtained a family of threedimensional axisymmetric mass distributions from the higher order Plummer models. The Miyamoto-Nagai potential shares many of the important properties of actual galaxies, especially the contour plots of the mass distribution which are qualitatively similar to the light distribution of disc galaxies (Binney & Tremaine 1987). Recently, two different extensions of the Miyamoto-Nagai potential appeared in the literature: a triaxial generalization (Vogt & Letelier 2007) which has as a particular case the original axially symmetric model, and a relativistic version (Vogt & Letelier 2005b) which has as a Newtonian limit the same original model.
In order to have a general relativistic physical model for galaxies, we must consider, first of all, the thickness of the disc and its stability under perturbations of the fluid quantities. The purpose of this work is to study numerically the stability of the general relativistic Miyamoto-Nagai disc under a general first-order perturbation. The perturbation is done in the temporal, radial, axial and azimuthal components of the quantities involved in the energy momentum tensor of the fluid. In the general thick disc case (Ujevic & Letelier 2007), the number of unknowns is larger than the number of equations. This opens the possibility of performing several types of combinations of the perturbed quantities. In this paper, we search for perturbations in which a perturbation in a given direction of the pressure creates a perturbation in the same direction of the four velocity. The energy momentum perturbation consid-ered in this paper is treated as 'test matter', so it does not modify the background metric obtained from the solution of Einstein equations.
The article is organized as follows. In Section 2, we present the general perturbed conservation equations for the thick disc case. The energy momentum tensor is considered diagonal with all its elements different from zero. Also, in particular, we discuss the perturbations that will be considered in some detail in the next sections of this work. In Section 3, we present the thick disc model whose stability is analysed, that is, the general relativistic Miyamoto-Nagai disc. The form of its energy density and pressures, as well as the restrictions that the thermodynamic quantities must obey to satisfy the strong, weak and dominant energy conditions are shown. Later, in Section 4, we perform the perturbations to the general relativistic Miyamoto-Nagai disc; in particular, we study its stability. Finally, in Section 5, we summarize our results.

P E RT U R B E D E Q UAT I O N S
The thick disc considered is a particular case of the general static axially symmetric metric ds 2 = −e 2 1 dt 2 + e 2 2 R 2 dθ 2 + e 2 3 (dR 2 + dz 2 ), where 1 , 2 and 3 are functions of the variables (R, z). (Our conventions are: G = c = 1, metric signature +2, partial and covariant derivatives with respect to the coordinate x μ denoted by , μ and ;μ, respectively.) In its rest frame, the energy momentum tensor of the fluid T μν is diagonal with components (−ρ, p R , p θ , p z ), where ρ is the total energy density and (p R , p θ , p z ) are the radial, azimuthal and axial pressures or tensions, respectively. So, in this frame of reference, the energy momentum tensor can be written as where U μ , X μ , Y μ and Z μ are the four vectors of the orthonormal tetrad: which satisfy the orthonormal relations. Note that with the above definitions, the time-like four velocity of the fluid is U μ and the quantities X μ , Y μ and Z μ are the space-like principal directions of the fluid. Furthermore, the energy momentum tensor satisfies Einstein field equations, G μν = 8πT μν . Moreover, the quantities involved in the energy momentum tensor and the coefficients of the perturbed conservation equations are functions of the coordinates (R, z) only. Let us consider a general perturbation A μ P of a quantity A μ of the form: where A μ (R, z) is the unperturbed quantity and δ A μ (R, z)e i(k θ θ−wt) is the perturbation. Replacing equation (4)  energy momentum equations, δT μν ;ν = 0, we obtain: where are the Christoffel symbols. In finding equations (5)-(8), we assumed that the perturbed energy momentum tensor does not modify the background metric. Also, we disregard terms of order greater than or equal to δ 2 (for details see Ujevic & Letelier 2004.
Besides the four equations furnished by the energy momentum conservation equations, T μν ;ν = 0, there is another important conservation equation, the equation of continuity, where n is the proper number density of particles. The proper number density of particles n and the total energy density ρ are related through the relation where m b is the constant mean baryon mass and ε is the internal energy density. Multiplying equation (12) by U μ , performing the covariant derivative (;μ) and using equation (11), we obtain that Now, from the relation U ν T μν ;μ = 0 and the energy momentum tensor (equation 2), we obtain an expression for (ρU μ ) ;μ . Substituting this last expression into equation (13), we finally arrive at which is a first-order differential equation for ε. Therefore, with ε given by equation (14) the equation of continuity (11) is satisfied. For this reason, the continuity equation can be omitted in our stability analysis because, in principle, we can always find a solution for ε.
Hereafter, the contribution of nm b and ε to the total energy density is taken into account in ρ. In the case in which the internal energy density of the fluid is given, the equation of continuity must be considered. The thermodynamic properties of the system can be obtained from observations or theoretically, for example, from the Fokker-Planck equation, where we obtain the particle distribution function of the disc. Solving the three-dimensional Fokker-Planck equation is not an easy task, but some progress in Newtonian gravity had been made (Ujevic & Letelier 2005, 2006. The four equations, (5)- (8), contain seven independent unknowns, say δU R , δU θ , δU z , δρ, δp R , δp θ and δp z . So, at this point, the number of unknowns is greater than the number of equations. This opens the possibility to perform different kind of perturbations. In this paper, we are interested in perturbations in which the velocity perturbation in a certain direction leads to a pressure perturbation in the same direction. For example, if we perturbed the axial component of the velocity, δU z , then we must perturb δp z . With the above criterion, and without imposing any extra conditions on the individual perturbations, only four perturbation combinations are allowed and will be considered in our thick disc model. Furthermore, we perform the perturbation δU R , δp R , δU z and δp z with the extra imposed condition δp R ≡ δp z . In this particular case, the system of equations reduces to a second-order partial differential equation.

G E N E R A L R E L AT I V I S T I C M I YA M OTO -NAG A I G A L A X I E S
A static general relativistic version of the Miyamoto-Nagai disc was constructed by Vogt & Letelier (2005b) by making a correspondence between the general isotropic line element in cylindrical coordinates and the Miyamoto-Nagai model (Miyamoto & Nagai 1975;Binney & Tremaine 1987). These general relativistic discs are obtained with equation (1) and the specializations , m is the mass of the disc and (a, b) are constants that control the shape of the density curves. With this metric, the energy density and pressures for the general relativistic Miyamoto-Nagai disc are where ξ = √ z 2 + b 2 and χ = R 2 + (a + ξ ) 2 . Without losing generality, we set m = 1 in equations (17)-(19). To satisfy the strong energy condition (gravitational attractive matter), we must have that the 'effective Newtonian density' = ρ + p R + p θ + p z 0. The weak energy condition requires ρ 0 and the dominant energy condition requires |p R /ρ| 1, |p θ /ρ| 1 and |p z /ρ| 1. The parameters used in this paper satisfy all energy conditions. Furthermore, the level curves show that it is physically acceptable. We remark that these are not the only parameters in which the level curves are physically acceptable. In the next section, we apply the selected perturbations of Section 2 to the general relativistic Miyamoto-Nagai disc mentioned above and study its stability.

P E RT U R BAT I O N S
Before applying the different kinds of perturbations to the general relativistic Miyamoto-Nagai disc, we must make some considerations. Note that the general relativistic Miyamoto-Nagai disc is infinite in the radial and the axial directions. We want to study the stability of a finite disc. So, in order to achieve this requirement, we need a cut-off in the radial coordinate. In equations (17)-(19), we see that the thermodynamic quantities decrease rapidly enough to define a cut-off in both coordinates. The radial cut-off R cut and the axial cut-off Z cut are set by the following criterion: the energy density within the disc formed by the cut-off parameters has to be more than 90 per cent of the infinite thick disc energy density. The above criterion, and the parameters used in the article, leads to a radial cut-off of R cut = 10 units and an axial cut-off of |Z cut | = 5 units. The other 10 per cent of the energy density that is distributed from outside the cut-off parameters to infinity can be treated, if necessary, as a perturbation in the outermost boundary condition.

Perturbation with δU θ , δp θ , δU R and δp R
We start perturbing the four velocity in its components θ and R. From the physical considerations mentioned in Section 2, we also expect variations in the thermodynamic quantities p θ and p R . The set of equations (5)-(8) reduces to a second-order ordinary differential equation for the perturbation δp R , say where (F A , F B , F C ) are functions of (R, z, w, k θ ) (see Appendix A). For this particular case, we have δp θ = −δp R . Note that in equation (20), the coordinate z only enters as a parameter. Moreover, the equation for δp R is independent of the parameter w, but w needs to be different from zero to reach that form. The second-order equation (20) is solved numerically with two boundary conditions, one at R ≈ 0 and the other at the radial cut-off. At R ≈ 0, we set the perturbation δp R to be ≈10 per cent of the unperturbed pressure p R (equation 18). In the outer radius of the disc, we set δ p R | R=Rcut = 0 because we want our perturbation to vanish when approaching the edge of the disc and, in that way, to be in accordance with the applied linear perturbation. We say that our perturbations are valid if their values are lower than, or of the same order of magnitude as, the 10 per cent values of its unperturbed quantities.
In Fig. 1, we present the amplitude profile of the radial pressure perturbation in the plane z = 0 for different values of the parameters a and b. As in the Newtonian case, the less the ratio b/a, the flatter is the mass distribution. We see that the perturbation δp R for (a = 1, b = 1) decreases rapidly with R and has oscillatory behaviour. At first sight, the perturbation δp R appears to be stable for all R, but in order to make a complete analysis, we have to compare at each radius the values of the perturbations with the values of the radial pressure. For this purpose, we included in the same graph a profile of the 10 per cent value of p R . We see that the perturbations of δp R for different values of k θ are always lower or, at least, of the same order of magnitude when compared to these 10 per cent values. In the flatter case (a = 1, b = 0.5), the perturbation δp R shows the same qualitative behaviour, but the amplitudes of the oscillations are slightly higher. In both cases, the amplitudes are well below the 10 per cent values of p R . If we consider a very flat galaxy (a = 1, b = 0.1) with w = 1, we find that some modes are not stable in a small region near the centre of the disc, from R ≈ 0 to R ≈ 0.3, because the perturbation amplitude is bigger than the 10 per cent value of p R and our general linear perturbation is no longer valid. We also performed stability analyses for the physical radial velocity perturbationδ U R = √ g R R δU R and the physical azimuthal velocity perturbationδ U θ = √ g θθ δU θ . Note that our four velocity U μ (equation 3) has only components in the temporal part, so we do not have values of U R and U θ to make comparisons with the perturbed valuesδ U R andδ U θ . For that reason we compared, in first approximation, the amplitude profiles of these perturbations with the value of the escape velocity in the Newtonian limit. In the Newtonian limit of general relativity, f 1, we have the wellknown relation g 00 = η 00 + 2 . So, the Newtonian escape velocity V esc = √ 2| | can be written as V esc = 2 √ f (see Vogt & Letelier 2005b). With this criterion, the perturbationsδ U r andδ U θ are stable because their values are always well below the escape velocity value. Recall that the perturbation δp R does not depend on the parameter w, but the perturbationsδ U R andδ U θ do. We performed numerical solutions for the perturbationsδ U R andδ U θ with different values of the frequency w, and we find that when we increase the value of w the perturbations become more stable.
In this section, we set the value of the parameter z = 0. We performed the same analysis for different values of the parameter −5 < z < 5, and we found that the perturbations show the same qualitative behaviour. Therefore, we can say that the general relativistic Miyamoto-Nagai disc shows some non-stable modes for very flat galaxies, for example, (a = 1, b = 0.1). Otherwise the disc is stable under perturbations of the form presented in this section.
Nevertheless, if we treat the 10 per cent of the energy density as a perturbation in the outermost radius of the disc by setting δ p R | R=Rcut = , where < 10 per cent of p R | R=Rcut , the qualitative behaviour of the mode profiles is the same. In the case of flat galaxies, when they present non-stable modes, more complex structures like rings, bars or haloes can be formed. Moreover, if we set the frequency w → iw, we obtain the same equation for the perturbation δp R , say, equation (20). In this case, the real part of the general perturbation diverges with time and the perturbation is not stable. These last considerations can be applied to every perturbation in the following sections.

Perturbation with δU θ , δp θ , δU z and δp z
In this section, we perturb the four velocity in its components θ and z, and we expect variations in the thermodynamic quantities p θ and p z . The set of equations (5)-(8) reduces to a second-order ordinary differential equation for the perturbation δp z given by where (F A , F B , F C ) are functions of (R, z, w, k θ ) (see Appendix B). Note that in equation (21), the coordinate R only enters as a parameter. Like the previous case, equation (21) is independent of Figure 2. Profiles of the amplitude perturbation for the axial pressure, axial physical velocity and azimuthal physical velocity of the fluid for the cases when R = 0.1 and k θ = 0, 5, 10, 15, 20. The graphs at the left-and the right-hand sides correspond to the cases when (a = 1, b = 1) and (a = 1, b = 0.5), respectively. The 10 per cent p z profile and the escape velocity are depicted for better stability comparisons. Note, in the axial pressure perturbation graph, that some nodes are not stable because they do not satisfy our stability criterion. In order to be stable, a mode must have the correct behaviour in all the perturbed quantities. For example, in the case (a = 1, b = 1), the mode with w = 1 and k θ = 5 seems to be stable in δp z , but looking into theδ U θ perturbation graph, we note that this statement is not true. In the flatter galaxies, the modes are not stable on larger regions of the domain.
the parameter w, but in order to reach that form we must have w different from zero. The second-order equation (21) is solved numerically with two boundary conditions, one in z = 0 and the other in z = Z cut . At z = 0, we set the perturbation δp z to be ≈10 per cent of the unperturbed pressure p z (equation 19). In the outer plane of the disc, we set δ p z | z=Zcut = 0 because we want our perturbation to vanish when approaching the edge of the disc, and in that way, to be in accordance with the linear perturbation applied. In Fig. 2, we present the amplitude profiles of the axial pressure perturbation, the physical axial velocity perturbationδ U z = √ g zz δU z and the physical azimuthal velocity perturbation for R = 0.1 and different values of the parameters a and b. For comparison reasons, we included in the graphs the amplitude profile that corresponds to 10 per cent of the value of p z and the escape velocity profile. Note that for (a = 1, b = 1) some modes of the axial pressure perturbation are above the 10 per cent profile of p z , for example, the modes with k θ = 0 and 5. In these cases, we can say that the mode with k θ = 0 is not stable and that the mode with k θ = 5 is near the validity criterion used for the perturbations. These modes are also present in the flatter galaxy (a = 1, b = 0.5) and have the same behaviours. The mode k θ = 5 is actually not stable. This can be seen in the azimuthal velocity perturbation profiles, where its amplitude is greater than the escape velocity. Note that in the velocity perturbation graphs, the mode k θ = 0 is also not stable. The azimuthal pressure perturbation, not depicted in Fig. 2, has all the modes well below the 10 per cent profile of p θ , and therefore is stable.
The perturbations δp z and δp θ do not depend on the parameter w, but the perturbationsδ U z andδ U θ do. We performed numerical solutions for the perturbationsδ U z andδ U θ with different values of the frequency w, and we find that when we increase the value of w the perturbations become more stable.
We have performed the above analysis for different values of the parameter 0 < R < 10, and we found that the qualitative behaviour is the same. We see from Fig. 2 that the non-stable modes are more pronounced for the flatter galaxy. Furthermore, for very flat galaxies some modes like k θ = 10 become non-stable. In general, for nonstable modes, more complex structures like rings, bars or spiral arms may be formed.

Perturbation with δU R , δp R and δρ
In this section, we perturb the radial component of the four velocity, the radial pressure and the energy density of the fluid. The set of equations (5)-(8) reduces to a second-order ordinary differential equation for the perturbation δp R of the form equation (20). The forms of the functions (F A , F B , F C ) are given in Appendix C. In this case, the coordinate z only enters as a parameter. Due to the fact that we are not considering perturbations in the azimuthal axis, the coefficients of the second-order ordinary differential equation do not depend on the wavenumber k θ . This second-order equation is solved numerically with the same boundary conditions described in Section 4.1.
In Fig. 3, we present the amplitude profiles for different perturbation modes of the radial pressure in the plane z = 0 for different values of the parameters a and b. We see in the graph that the perturbation profiles decrease rapidly in a few units of R. Also, the values of the radial velocity perturbation and energy density perturbation, not depicted, are well below the escape velocity and the 10 per cent energy profile, respectively.
We performed the above analysis for different values of −5 < z < 5, and we found that the quantities involved have the same qualitative behaviour. From these results, we can say that the general linear perturbation applied is highly stable and, for that reason, the perturbations do not form more complex structures. δp R w=1 w=5 w=10 w=15 w=20 10% of p R Figure 3. Profiles of the amplitude perturbation for the radial pressure of the fluid for the cases when z = 0 and w = 1, 5, 10, 15, 20. The graphs at the left-and the right-hand sides correspond to the cases when (a = 1, b = 1) and (a = 1, b = 0.5), respectively. The 10 per cent p R profile is depicted for better stability comparisons. These modes are highly stable.

Perturbation with δU z , δp z and δρ
In this section, we perturb the axial component of the four velocity, the axial component of the pressure and the energy density of the fluid. The set of equations (5)-(8) reduces to a second-order ordinary differential equation for the perturbation δp z of the form of equation (21). The functions (F A , F B , F C ) are given in Appendix D. Note that, like in Section 4.2, the coordinate R only enters as a parameter.
In this case, we are not considering azimuthal perturbations, and therefore the quantities involved do not depend on the parameter k θ . The second-order equation is solved following the procedure of Section 4.2.
In Fig. 4, we present the amplitude profiles of the axial pressure perturbation and the physical axial velocity perturbation, for R = 0.1 and for different values of the parameters a and b. We see that the axial pressure perturbation modes for (a = 1, b = 1) are always of the same order of magnitude or lower when compared to the 10 per cent profile. In the flatter case (a = 1, b = 0.5), note that the amplitude of the mode w = 1 is greater in some region of the domain. This fact is reflected in the axial velocity perturbation profile where the mode w = 1 has a strange behaviour. All of the modes, including the mode with w = 1, are stable because they are well below the escape velocity, which is not depicted. The modes that correspond to the energy density perturbation are all stable. For highly flat galaxies, the mode w = 1 is not stable and may form more complex structures.
For higher values of the parameter w the modes are more stable. We performed the above analysis for different values of the parameter 0 < R < 10, and we found that the quantities involved have the same qualitative behaviour.

Perturbation with δU R , δp R , δU z , δp z and δp R ≡ δp z
In this section, we perturb the radial component of the four velocity, the axial component of the four velocity, the radial pressure and the axial pressure. As we said in Section 2, we need an extra condition to set the number of unknowns equal to the number of equations. In this case, we set δp R ≡ δp z ≡ δp. Therefore, the set of equations (5)-(8) reduces to a second-order partial differential equation for the pressure perturbation δp, say where (F A , F B , F C , F D , F E ) are functions of (R, z, w) (see Appendix E). The partial differential equation (22) 1,5,10,15,20. The graphs at the left-and the right-hand sides correspond to the cases when (a = 1, b = 1) and (a = 1, b = 0.5), respectively. The 10 per cent p z profile is depicted in the graph of δp z for better stability comparisons. In the graph ofδ U z , the escape velocity is not depicted because it is several orders of magnitude greater.
For these examples, all the modes are stable, but for a highly flat galaxy some modes, like the mode with w = 1, are not stable.  Figure 5. Profiles of the amplitude perturbation for the pressure, the radial physical velocity and the axial physical velocity of the fluid for the case when w = 1 and for a rod perturbation in R ≈ 0. With this kind of perturbation, the disc tends to form some kind of ring around the centre of the disc. This phenomenon is greater for highly flat galaxies and lower for more spherical systems. R ≈ 0 and R = R cut . They are different ways in which we can set the boundary conditions in order to simulate various kinds of pressure perturbations. Here, we treat only the case when we have a pressure perturbation at R ≈ 0 and along the z axis, that is, some kind of a rod perturbation. We set the value of the rod pressure perturbation to be 10 per cent of the axial pressure. We set the values of the other boundary conditions equal to zero, because we want the perturbation to vanish when approaching the edge of the disc. We choose the 10 per cent of the value of the axial pressure instead of the radial pressure, because it has the lowest value near R ≈ 0. In that way, the perturbation values are also below the 10 per cent values of the radial pressure and the general linear perturbation is valid.
In Fig. 5, we present the perturbation amplitudes for the pressure, the physical radial velocity and the physical axial velocity, for w = 1 and for different values of the parameters a and b. We see in the pressure perturbation graph that the perturbation rapidly decays to values near zero when we move out from the centre of the disc. This behaviour is the same for every galaxy considered. In the velocity perturbation profiles, we can see a phenomenon that is more clear in the flatter galaxy. Note that in the lower domain of the disc [−5,0) the axial velocity perturbation is positive and in the upper domain (0,5] the axial velocity perturbation is negative. This means that due to the linear perturbation, the disc tries to collapse to the plane z = 0. Now, if we look to the radial velocity perturbation graph, we note that the upper and the lower parts depart from the centre of the disc due to the positive radial perturbation. So, with these considerations, we may say that the disc tends to form some kind of ring around the centre of the disc. This phenomenon is greater for highly flat galaxies and lower for more spherical systems.

C O N C L U S I O N S
In this paper, we studied the stability of the recently proposed general relativistic Miyamoto-Nagai model (Vogt & Letelier 2005b) by applying a general first-order perturbation. We can say that the stability analysis performed is more complete than the stability analysis of particle motion along geodesics because we have taken into account the collective behaviour of the particles. However, this analysis can be said to be incomplete because the energy momentum perturbation tensor of the fluid is treated as a test fluid and does not alter the background metric. This is a second degree of approximation to the stability problem in which the emission of gravitational radiation is considered.
The different stability analyses made of the general relativistic Miyamoto-Nagai disc show that this disc is stable for higher values of the wavenumber k θ and the frequency w. For lower values of k θ and w, the disc presents non-stable modes that may form more complex structures like rings, bars or haloes, but in order to study them we need a higher order perturbation formalism. In general, nonstable modes appear more for flatter galaxies and less for spherical systems.

AC K N OW L E D G M E N T S
MU and PSL thank FAPESP for financial support; PSL also thanks CNPq.

A P P E N D I X A : F U N C T I O N S F A , F B A N D F C O F S E C T I O N 4 . 1
The general form of the functions (F A , F B , F C ) appearing in the second-order ordinary differential equation (20) is given by F A = A 1 α 1 , F B = A 1 α 1,R + A 1 α 2 + A 3 α 1 , where α 1 , α 2 and α 3 are In equations (A1) and (A2), we denote the coefficients of equation (5) by A i , the coefficient of equation (6) by B i , the coefficient of equation (7) by C i , the coefficient of equation (8) by D i , for example, the first term in equation (5) has the coefficient A 1 multiplied by the factor of δU R ,R , the second term has the coefficient A 2 multiplied by the factor of δU R , etc. The explicit form of the above equations is obtained replacing the fluid variables (ρ, p R , p θ , p z ) of the isotropic Schwarzschild thick disc.

A P P E N D I X B : F U N C T I O N S F A , F B A N D F C O F S E C T I O N 4 . 2
The general form of the functions (F A , F B , F C ) appearing in the second-order ordinary differential equation (21) is given by where α 1 , α 2 and α 3 are and the meaning of the coefficients (A i , B i , C i , D i ) is explained in Appendix A.

A P P E N D I X C : F U N C T I O N S F A , F B A N D F C O F S E C T I O N 4 . 3
The general form of the functions (F A , F B , F C ) is given by where α 1 , α 2 and α 3 are and the meaning of the coefficients (A i , B i , D i ) is explained in Appendix A.

A P P E N D I X D : F U N C T I O N S F A , F B A N D F C O F S E C T I O N 4 . 4
The general form of the functions (F A , F B , F C ) is given by where α 1 , α 2 and α 3 are and the meaning of the coefficients (A i , B i , D i ) is explained in Appendix A.

A P P E N D I X E : F U N C T I O N S F A , F B , F C , F D A N D F E O F S E C T I O N 4 . 5
The general form of the functions (F A , F B , F C , F D , F E ) appearing in the partial second-order differential equation (22) is given by F A = A 1 α 1 , F B = A 1 α 1,R + A 1 α 2 + A 3 α 1 , F C = A 2 α 3 , F D = A 2 α 3,z + A 2 α 4 + A 5 α 3 , where α 1 , α 2 , α 3 and α 4 are and the meaning of the coefficients (A i , B i , D i ) is explained in Appendix A.
This paper has been typeset from a T E X/L A T E X file prepared by the author.