2 THE POTENTIAL FORMULA

The gravitational potential in stellar systems (subsystems), such as, for instance bulges and discs of spiral galaxies, is often represented by a variant of the Green formula (e.g. Cuddeford Reference Cuddeford1993). This would mean that in general the gravitational potential Φ should have the following form:

(1) $$\begin{equation} {\Phi }={G {\mathcal M} \over {\mathcal R}}, \end{equation}$$

where G and ${\mathcal M}$ are constants, whereas ${\mathcal R}$ is a function. The first constant (G) is the universal gravitation constant, the second one (${\mathcal M}$) is the total mass of the system (subsystem) appearing as the source of the gravitational field. ${\mathcal R}$ is a function of the generalised coordinates and time. It must satisfy the following conditions: to have dimension of length and to be approximately equal to the distance from the centre of the source system r at the points lying very far from the system centre.

The potential of a flattened stellar system, as here of interest, is assumed to depend on two arguments: R, distance to the axis of symmetry z, and |z|, distance to the midplane. The function proposed by MN is a well-known example

(2) $$\begin{equation} {\mathcal R}_{MN}=\bigg [R^2 + \Big (a+ \sqrt{z^2+b^2}\Big )^2\bigg ]^{1/2} . \end{equation}$$

In this function, a and b are constants and for significantly flattened systems satisfy b/a < 1. When the function ${\mathcal R}_{MN}$ is substituted in Equation (1), one obtains the Miyamoto–Nagai potential formula. The present author proposes here a potential formula which contains a new term ${\mathcal R}_N$, as follows:

(3) $$\begin{equation} \Phi ={G{\mathcal M} \over {{\mathcal R}_{MN}-{\mathcal R}_N}}. \end{equation}$$

In a general case, the new term ${\mathcal R}_N$ is also a function of the same two variables, R and |z| and is everywhere positive. Since ${\mathcal R}_{MN}$ satisfies the condition of being approximately equal to the distance r at the points lying very far from the system centre, the new term at such points must be negligibly small compared to ${\mathcal R}_{MN}$. In this way, the difference ${\mathcal R}_{MN}-{{\mathcal R}_N}$ will satisfy the condition of being approximately equal to the distance from the centre r at very distant points.

3 THE ROTATION CURVE

The surface density of the disc of a spiral galaxy, bearing in mind the observational constraints, is usually assumed to be exponential (often referred to as Freeman law), for instance its luminosity or mass density ρ obeys the following formula [e.g. Deg & Widrow Reference Deg and Widrow2013, their Equation (13)]:

(4) $$\begin{equation} \rho (R,z)=\rho (0) exp\bigg (-{R \over R_d}\bigg ) sech^2\bigg ({z \over z_d}\bigg ), \end{equation}$$

where R_d and z_d are constants. As easily seen, the surface density following from this expression will depend on R as a simple exponential function with R_d as the scale length.

Unfortunately, Equation (4) yields no analytical solutions for the potential via the Poisson equation. This is important because with the known potential, it is possible to obtain the circular speed u_c by using the well-known relation

(5) $$\begin{equation} u_c=\sqrt{-R{\partial \Phi \over {\partial R}}} , \ \ z=0 . \end{equation}$$

The plot of circular speed versus distance R according to Freeman’s solution (Figure 1) is characterised by existence of radius at which the circular speed is the same for the Keplerian and exponential disc cases of the same mass. As a consequence, at infinity the circular speed for the exponential disc approaches the Keplerian curve from ‘above’, i.e. from higher values of u_c. This is not the case with the potential formula of Miyamoto and Nagai. Substituting Equation (2) in Equation (1) and applying Equation (5), one obtains the dependence of the circular speed on R corresponding to the Miyamoto–Nagai potential. Since z = 0, the parameters a and b enter the circular-speed dependence always through their sum a + b. In other words, the flattening expressed by means of b/a, does not affect the circular-speed dependence. The maximum circular speed occurs at about 1.4(a + b) and the circular-speed value resulting from the Miyamoto–Nagai formula is everywhere smaller than the value yielded by the Keplerian dependence for the same total mass. Due to this, a good fit between the circular speed following from the potential of Miyamoto and Nagai and that corresponding to the exponential disc is achieved at the cost of different total masses. Usually the total mass in the Miyamoto–Nagai formula is about 1.5 times as large as that of the exponential disc (e.g. Ninković Reference Ninković1992).

Figure 1.

Rotation curves from Equation (7) as described in the text (blue curve γ₁ = −0.1, red curve γ₁ = −0.3) together with Freeman’s (Reference Freeman1970) curve; ξ = R/R_d, $\chi =u_c /\sqrt{{G{\mathcal M}\over R_d}}.$

The purpose of introducing the new term ${{\mathcal R}_N}$ [equation (3)] is to reproduce the property of intersection with the Keplerian curve and to have the same total mass as for the exponential disc. As shown earlier (Ninković Reference Ninković, Boily, Patsis, Portegies Zwart, Spurzem and Theis2003a), this is possible even with a constant substituted for ${{\mathcal R}_N}$. The constant introduced here will be R_d, the scale length. Then in fitting the rotation curve, one should determine the ratio (a + b)/R_d. According to the definition assumed in the present paper [Equation (1)], the potential cannot have negative values, so this ratio must be greater than 1. We find a best fit value of 2.1. The fit can be further improved. This is done by generalising the term ${{\mathcal R}_N}$ in the following way:

(6) $$\begin{equation} {\mathcal R}_N={1 \over 2}R_d\Bigg [\bigg (1+{R^2 \over c_1^2}\bigg )^{\gamma _1}+\bigg (1+{z^2 \over c_2^2}\bigg )^{\gamma _2}\Bigg ] , \end{equation}$$

where γ₁ < 0.5, γ₂ < 0.5. Though it may seem that the number of parameters now tends to be too large, this is not the case in practice. For instance, R_d is quite acceptable to be substituted for c ₁. The second term in Equation (6) (the function of |z|) is foreseen because of the density distribution the potential returns (see next section). With regard to all equations written above except Equation (4), one obtains the following expression which yields the circular speed:

(7) $$\begin{eqnarray} u_c&=&\sqrt{{G{\mathcal M}R \over \Bigg [\sqrt{R^2 + (a+b)^2}-{\frac{1}{2}R_d\bigg (1+\Big (1+{R^2 \over R_d^2}\Big )^{\gamma _1}\bigg )}\Bigg ]^2}}\cdot \nonumber \\ &&\cdot \sqrt{{{R \over \sqrt{R^2 + (a+b)^2}}}-{\gamma _1 {R \over R_d} \Bigg (1+{R^2 \over R_d^2}\Bigg )^{\gamma _1-1}}}. \end{eqnarray}$$

The curves presented in Figure 1 follow from Equation (7) with a + b = 2.1R_d, γ₁ = −0.1, and γ₁ = −0.3. They are presented together with the corresponding curve from Freeman’s (Reference Freeman1970) paper. The dimensionless variables, ξ and χ, are defined as: ξ = R/R_d, $\chi =\sqrt{G{\mathcal M}/R_d}$. The comparison requires the scale lengths R_d and the total masses ${\mathcal M}$ to be equal in both formalisms.

4 THE DENSITY

The density which generates a gravitational potential and the potential are related through the Poisson equation

(8) $$\begin{equation} \nabla ^2 \Phi =-4\pi G \rho . \end{equation}$$

∇² is the Laplacian, ρ is the density. The full formula is given in the Appendix.

In the density calculation, the constants c ₁ and c ₂ [Equation (6)], as well as the ratio b/a, must be specified. The fit for the rotation curve has yielded a + b = 2.1R_d, under the condition b ≪ a (flattened mass distribution), so it is easy to conclude that a is approximately 2R_d. In such conditions, the first term within the brackets in Equation (6) should not differ substantially from 1. As for the second one, a too small c ₂ can contribute to significant values of ${\mathcal R}_N$, for instance c ₂ = b, leads to such a high ${\mathcal R}_N$ that the denominator in Equation (3) becomes negative and, as a consequence, the potential, contrary to the convention assumed here, becomes also negative. Thus, it is reasonable to assume c ₁ = c ₂ = R_d.

In spite of introducing the new term ${\mathcal R}_N$ [Equations (3) and (6)], negative density values cannot be avoided. In the midplane, the density is nowhere negative. Its dependence on R in z = 0 is a monotonously decreasing function. As for the dependence ρ(z), R = const, there are three types of profiles: (i) with a minimum ρ_min < 0 and then approaching zero from the negative side (Figure 2); (ii) a wavy profile, with one minimum and one maximum and then approaching zero, along R = 0 both extrema can be positive, in general the minimum is negative, the maximum positive (Figure 3); (iii) a sound profile, without negative values and monotonously decreasing, but obtained along R = 0 only (Figure 4). In order to be more clearly visible, the most essential part of the curve is magnified and given in the same figure (Figures 2 and 3).

Figure 2.

Density dependence on height |z|, for R = 3R_d, γ₁ = γ₂ = 0 [see Equation (6)]; distance unit R_d, density unit $\mathcal {M}$R ⁻³_d.

Figure 3.

Density dependence on height |z|, for R = R_d, γ₁ = −0.3, γ₂ = 0.3 [see Equation (6)]; distance unit R_d, density unit $\mathcal {M}$R ⁻³_d.

Figure 4.

Density dependence on height |z|, for R = 0, γ₁ = −0.3, γ₂ = 0.3 [see Equation (6)]; distance unit R_d, density unit $\mathcal {M}$R ⁻³_d.

In the simplest case—γ₁ = γ₂ = 0 [equation (6)]—the density dependence ρ(z), R = const, is always of type (i) (Figure 2). The z value at which the density minimum occurs depends on R, the smaller R is, the closer is the minimum to the midplane. Finally, when R tends to infinity, the z value also tends to infinity, because the density tends to zero.

When the term ${\mathcal R}_N$ [Equations (3) and (6)] takes a functional form, i.e. both γ₁ and γ₂ are non-zero, the dependence ρ(z), R = const, changes the form, towards the types (ii) and (iii). For γ₁, some kind of critical value appears at which the dependence ρ(z), R = 0, obtains the form (iii) (Figure 4). It should be γ₁ < 0, the modulus is affected by the ratio b/a. For instance, for the pair of values b = 0.18, a = 1.92, the critical value is γ₁ = −0.3, for b = 0.08, a = 2.02, the corresponding is γ₁ = −0.5, for b = 0.02, a = 2.08, γ₁ = −0.6, etc. The influence of γ₂ is not so strong; it is better to be γ₂ > 0, but increasing γ₂ leads to weaker decreasing of ρ with increasing z. The most important is that increasing the modulus of γ₁ further does not improve the situation, as the negative density values remain and the wavy profiles survive (Figure 3). However, the moduli of both exponents γ₁ and γ₂ are limited. As for the former, the main limitation comes from the fit of the rotation curve, in the case of the latter, since it is supposed to be positive, the physical limitation (γ₂ < 0.5) becomes essential. The final conclusion is that the best achievable result is to obtain a monotonously decreasing density along R = 0 and wavy density profiles along cylinders R > 0. For the purpose of giving a more clear explanation the following example is chosen. The values of the parameters in Equations (3) and (6) are: a + b = 2.1R_d (b = 0.18R_d), c ₁ = c ₂ = R_d, γ₁ = −0.3, γ₂ = 0.3. The rotation curve (Figure 1—red curve) and the density profiles (Figures 3 and 4) correspond to this example.

In Figure 2, we show another density comparison, typical for γ₁ = γ₂ = 0, corresponding to the same a + b and the ratio b/a. The wave (Figure 3) can be explained by the rather rapid density decrease with increasing distance to the midplane. Beyond the centre, the density values in the midplane are generally low, significantly lower than at the centre, so that at reasonable distances to the midplane they are practically zero. Here, one obtains the density values by using the Poisson Equation (8) which in the case of axial symmetry means as the algebraic sum of three terms (see Appendix). Clearly, in a numerical procedure instead of density values of exactly zero we usually have values of low moduli with both signs. In this way, the wavy density profiles obtained here can be easily understood as due to numerical precision issues.