1 INTRODUCTION
The gravitational potential proposed by Miyamoto and Nagai (Reference Miyamoto and Nagai1975) (MN) is well known. Its applications to models of the Milky Way and other disc galaxies are numerous (e.g. Allen & Santillan Reference Allen and Santillan1991, Ninković Reference Ninković1992, Dinescu, Girard, & van Altena Reference Dinescu, Girard and van Altena1999, Binney & MacMillan Reference Binney and McMillan2011). The observational evidence is in favour of exponential discs characterised by the dependence of the circular speed on the distance to the centre (rotation curve) as given in Freeman’s paper (Reference Freeman1970). As to the fitting of the circular-speed dependence resulting from the Miyamoto–Nagai formula to Freeman’s curve, there are some difficulties indicated in an earlier paper of the present author (Ninković Reference Ninković, Boily, Patsis, Portegies Zwart, Spurzem and Theis2003a). In the present paper, a new formula for the gravitational potential is proposed. It appears as a generalisation of MN and it contains the formula analysed in the earlier paper (Ninković Reference Ninković, Boily, Patsis, Portegies Zwart, Spurzem and Theis2003a) as a special case. The density generating this potential as function of the coordinates is also studied.
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:
$$\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
$$\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:
$$\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)]:
$$\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 Rd and zd are constants. As easily seen, the surface density following from this expression will depend on R as a simple exponential function with Rd 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 uc by using the well-known relation
$$\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 uc. 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 γ1 = −0.1, red curve γ1 = −0.3) together with Freeman’s (Reference Freeman1970) curve; ξ = R/Rd, 
$\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 Rd, the scale length. Then in fitting the rotation curve, one should determine the ratio (a + b)/Rd. 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:
$$\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 γ1 < 0.5, γ2 < 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, Rd is quite acceptable to be substituted for c 1. 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:
$$\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.1Rd, γ1 = −0.1, and γ1 = −0.3. They are presented together with the corresponding curve from Freeman’s (Reference Freeman1970) paper. The dimensionless variables, ξ and χ, are defined as: ξ = R/Rd, 
$\chi =\sqrt{G{\mathcal M}/R_d}$. The comparison requires the scale lengths Rd 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
$$\begin{equation}
\nabla ^2 \Phi =-4\pi G \rho .
\end{equation}$$
∇2 is the Laplacian, ρ is the density. The full formula is given in the Appendix.
In the density calculation, the constants c 1 and c 2 [Equation (6)], as well as the ratio b/a, must be specified. The fit for the rotation curve has yielded a + b = 2.1Rd, under the condition b ≪ a (flattened mass distribution), so it is easy to conclude that a is approximately 2Rd. 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 2 can contribute to significant values of ${\mathcal R}_N$, for instance c 2 = 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 1 = c 2 = Rd.
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 = 3Rd, γ1 = γ2 = 0 [see Equation (6)]; distance unit Rd, density unit $\mathcal {M}$R −3d.
Figure 3. Density dependence on height |z|, for R = Rd, γ1 = −0.3, γ2 = 0.3 [see Equation (6)]; distance unit Rd, density unit $\mathcal {M}$R −3d.
Figure 4. Density dependence on height |z|, for R = 0, γ1 = −0.3, γ2 = 0.3 [see Equation (6)]; distance unit Rd, density unit $\mathcal {M}$R −3d.
In the simplest case—γ1 = γ2 = 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 γ1 and γ2 are non-zero, the dependence ρ(z), R = const, changes the form, towards the types (ii) and (iii). For γ1, some kind of critical value appears at which the dependence ρ(z), R = 0, obtains the form (iii) (Figure 4). It should be γ1 < 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 γ1 = −0.3, for b = 0.08, a = 2.02, the corresponding is γ1 = −0.5, for b = 0.02, a = 2.08, γ1 = −0.6, etc. The influence of γ2 is not so strong; it is better to be γ2 > 0, but increasing γ2 leads to weaker decreasing of ρ with increasing z. The most important is that increasing the modulus of γ1 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 γ1 and γ2 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 (γ2 < 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.1Rd (b = 0.18Rd), c 1 = c 2 = Rd, γ1 = −0.3, γ2 = 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 γ1 = γ2 = 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.
5 DISCUSSION AND CONCLUSIONS
The present author proposes an analytical form for the gravitational potential which would correspond to the exponential disc, the most luminous subsystem of a spiral galaxy. This is done by modifying the well-known formula of Miyamoto and Nagai [Equations (1)–(3) and (6)]. In this modification, there are three essential parameters: the total mass $\mathcal M$, the exponential scale length Rd, and the other scale length b. The other three parameters—a, γ1, and γ2 are auxiliary. Their role is to improve the fit of the rotation curve particularly, their values are limited. Since the system (subsystem) under study is flattened, it must satisfy a ≫ b, and bearing in mind the constraint based on the rotation curve, one obtains a ≈ 2Rd. The exponents γ1 and γ2 cannot exceed 0.5 and γ1 is constrained to be small and negative, while γ2 lies in 0 < γ2 < 0.5. Thus, an exponential disc modelled in the way as proposed here should be generally characterised by its total mass and two scale lengths.
Though negative density values cannot be avoided completely, the functions in the new term [Equation (6)] contribute to a significant reducing of this effect. A discussion concerning the analogous spherical-symmetry model can be found in Ninković (Reference Ninković2003b).
ACKNOWLEDGEMENTS
This research has been supported by the Serbian Ministry of Education, Science and Technological Development (Project No 176011 ‘Dynamics and Kinematics of Celestial Bodies and Systems’). The author wants to thank an anonymous referee for valuable suggestions which contributed to a better presentation of the results.
A APPENDIX
Equations (1) and (3) from the main text are given again
$$\begin{eqnarray*}
\Phi & =& {G{\mathcal M} \over {{\mathcal R}}} ,\\
{\mathcal R}& =& {\mathcal R}_{MN}-{{\mathcal R}_N} .
\end{eqnarray*}$$
The designations used are: Φ, the gravitational potential, G the universal gravitation constant, whereas ${\mathcal R}_{MN}$ and ${\mathcal R}_N$ are two functions depending on the arguments R and z (axial symmetry).
The Poisson Equation [(8) from the main text] relating the density and potential is also rewritten
$$\begin{equation*}\nabla ^2 \Phi =-4\pi G \rho .
\end{equation*}$$
With regard to the axial symmetry, the Laplacian has the form
$$\begin{equation*}
\nabla ^2\Phi ={\partial ^2\Phi \over {\partial R^2}}+{1\over R}{\partial \Phi \over {\partial R}}+{\partial ^2\Phi \over {\partial z^2}} .
\end{equation*}$$
The derivatives will be
$$\begin{equation*}{\partial \Phi \over {\partial R}} = -{G{\mathcal M} \over {\mathcal R}^2} \Bigg [ {\partial {\mathcal R}_{MN} \over {\partial R}} - {\partial {\mathcal R}_N \over {\partial R}}\Bigg ] ;
\end{equation*}$$
$$\begin{eqnarray*}
{\partial ^2 \Phi \over {\partial R^2}}&=& {2G{\mathcal M} \over {\mathcal R}^3} \Bigg [ {\partial {\mathcal R}_{MN} \over \partial R}- {\partial {\mathcal R}_N \over \partial R}\Bigg ]^2\\
&&-{G{\mathcal M} \over {{\mathcal R}^2}} \Bigg [{\partial ^2{\mathcal R}_{MN} \over {\partial R^2}}- {\partial ^2{\mathcal R}_N \over {\partial R^2}}\Bigg ] ;
\end{eqnarray*}$$
$$\begin{eqnarray*}
{\partial ^2 \Phi \over {\partial z^2}}&=& {2G{\mathcal M} \over {\mathcal R}^3} \Bigg [ {\partial {\mathcal R}_{MN} \over \partial z}- {\partial {\mathcal R}_N \over \partial z}\Bigg ]^2\\
&&-{G{\mathcal M} \over {{\mathcal R}^2}} \Bigg [{\partial ^2{\mathcal R}_{MN} \over {\partial z^2}}- {\partial ^2{\mathcal R}_N \over {\partial z^2}}\Bigg ] .
\end{eqnarray*}$$
The first function in the denominator ${\mathcal R}_{MN}$ and its necessary derivatives are
$$\begin{equation*}{\mathcal R}_{MN}=\bigg [R^2 + \Big (a+ \sqrt{z^2+b^2}\Big )^2\bigg ]^{1/2} ;
\end{equation*}$$
$$\begin{equation*}{\partial {\mathcal R}_{MN} \over {\partial R}}={R \over {{\mathcal R}_{MN}}} ;
\end{equation*}$$
$$\begin{equation*}{\partial ^2{\mathcal R}_{MN} \over {\partial R}^2}={{{\mathcal R}_{MN}-{R^2 \over {\mathcal R}_{MN}}} \over {{\mathcal R}_{MN}^2}} ;
\end{equation*}$$
$$\begin{equation*}q=a+\sqrt{z^2+b^2} ;
\end{equation*}$$
$$\begin{equation*}{\partial {\mathcal R}_{MN} \over \partial z}= {q \over {{\mathcal R}_{MN}}} {\partial q \over \partial z} ;
\end{equation*}$$
$$\begin{equation*}{\partial ^2{\mathcal R}_{MN} \over {\partial z^2}}= {q_z{{\mathcal R}_{MN}-q_z ({q^2 / {\mathcal R}_{MN}})} \over {{\mathcal R}_{MN}^2}}q_z + {q \over {\mathcal R}_{MN}}q_{zz} ;
\end{equation*}$$
$$\begin{equation*}q_z \equiv {\partial q \over {\partial z}}={z \over {\sqrt{b^2+z^2}}} , \ \ q_{zz}\equiv {\partial ^2q \over {\partial z^2}}={b^2 \over {(b^2+z^2)^{3/2}}} .
\end{equation*}$$
The second function in the denominator ${\mathcal R}_N$ and its necessary derivatives are
$$\begin{equation*}{\mathcal R}_N={1 \over 2}R_d\Bigg [\bigg (1+{R^2 \over R_d^2}\bigg )^{\gamma _1}+\bigg (1+{z^2 \over R_d^2}\bigg )^{\gamma _2}\Bigg ] ;
\end{equation*}$$
$$\begin{equation*}{\partial {\mathcal R}_N \over {\partial R}}= \gamma _1{R \over {R_d}}\bigg (1+{R^2 \over {R_d^2}}\bigg )^{\gamma _1-1} ;
\end{equation*}$$
$$\begin{eqnarray*}
{\partial ^2{\mathcal R}_N \over {\partial R^2}}&=& 2\gamma _1(\gamma _1-1){R^2 \over {R_d}^3}\bigg (1+{R^2 \over {R_d^2}}\bigg )^{\gamma _1-2}\\
&&+{\gamma _1 \over {R_d}}\bigg (1+{R^2 \over {R_d^2}}\bigg )^{\gamma _1-1} ;
\end{eqnarray*}$$
$$\begin{equation*}{\partial {\mathcal R}_N \over {\partial z}}=\gamma _2{z \over {R_d}}\bigg (1+{z^2 \over {R_d^2}}\bigg )^{\gamma _2-1} ;
\end{equation*}$$
$$\begin{eqnarray*}
{\partial ^2{\mathcal R}_N \over {\partial z^2}}&=& 2\gamma _2(\gamma _2-1){z^2 \over {R_d}^3}\bigg (1+{z^2 \over {R_d^2}}\bigg )^{\gamma _2-2}\\
&&+{\gamma _2 \over {R_d}}\bigg (1+{z^2 \over {R_d^2}}\bigg )^{\gamma _2-1} .
\end{eqnarray*}$$
