Introduction
Since the seminal paper by Miyamoto & Nagai (1975), the literature on three-dimensional analytical models for the gravitational field of different types of galaxies has grown considerably. In this respect, particular attention deserve the models proposed by Jaffe (1983) and Hernquist (1990), who derived analytical models which closely approximate the light distribution for spherical and elliptical galaxies, respectively. A few years later, Long & Murali (1992) presented an analytical potential for barred galaxies that reduces to the Miyamoto-Nagai disk by an appropriate setting of the free parameters, while Dehnen (1993) generalized the Jaffe and Hernquist models by means of a family of density profiles with different central slopes. More recently, Vogt & Letelier (2005) derived an analytical expression for the gravitational field of galaxies, based on the multipolar expansion up to the quadrupole term. Using a different approach, González et al. (2010) obtained a family of finite thindiscs models for four galaxies in the Ursa major cluster in which the circular velocities were adjusted to fit the observed rotation curves.
One advantage of an analytical galaxy model is the possibility to study the dynamics (regular or chaotic) of orbits. This can be considered one of the standing problems in galactic dynamics because it could allow us to understand the formation and evolution of galaxies (Contopoulos, 1979), as shown by the pioneer simulations of Lindblad (1960). Despite the fact that early papers on this topic studied only regular orbits in the meridional plane (Martinet & Mayer, 1975, Manabe, 1979, Greiner, 1987, Lees & Schwarzschild, 1992), soon after, the existence of chaos on the orbital motion started to beconsidered by Caranicolas (1996) and Caranicolas & Papadopoulos (2003). In the majority of cases all these studies focused on the distinction between regular and chaotic orbits (Manos & Athanassoula, 2011, Bountis ET. AL., 2012, Manos et al., 2013) or the influence of the galaxy components (nucleus, bulge, disk, halo) on the character of orbits, see e.g. (Zotos, 2012, Zotos & Caranicolas, 2013, Zotos, 2014). Notwithstanding the evidence that both chaotic and regular motions are possible in many axisymmetric potentials, recent studies on generalized axisymmetric potentials suggest that a third integral of motion seems to exist for energy values closer to the escape energy (Dubeibe et al., 2018, Zotos et al., 2018). Hence, such apparent ambiguity might only be solved by performing systematic studies of each particular model.
In this paper, we are interested in meridional motions of free test particles (stars) in presence of analytical realistic galaxy models. Our models possess axial symmetry, which is a good approximation given the morphology of galaxies that are mainly approximate figures of revolution. Additionally, the galaxy components were not added one by one, instead of this, we derived a generalized Miyamoto-Nagai model that can be adjusted very accurately to fit the observed rotation curve and hence it is assumed that all (or most of) the components are taken into account. The determination of the specific values of the coefficients of the series expansion let us calculate the corresponding surface densities and all the kinematic quantities characterizing the particular galaxy models. Unlike the models derived by González et al. (2010), which exhibit instabilities to small vertical perturbations (see e.g. the cases of NGC 3877 and NGC 4010), our models satisfy the stability conditions for radial and vertical perturbations. On the other hand, the dynamics of the orbits is studied through the Poincaré surfaces of section, showing that the orbital motion exhibits a strong dependence on the angular momentum and energy of the test particles (stars).
The paper is organized as follows: in the first section, we derive the generalized Miyamoto-Nagai model. Next, from the new potential the explicit expressions for the physical quantities of interest are determined. In the second section we adjust the observed rotation curves of three specific spiral galaxies (NGC 3726, NGC 3877 and NGC 4010) to the analytical circular velocities derived with our model. Then, the massdensity profiles are calculated, along with the vertical and epicyclic frequencies, showing that our model not only is well-behaved but also satisfy the stability conditions. A dynamical analysis in terms of the Poincaré surfaces of section is performed in the third section. Finally, in the fourth section, we summarize our main conclusions.
Generalized Miyamoto-Nagai model
Let us start considering the axially symmetric Laplace’s equation in spherical coordinates
whose general solution reads as
where A l and B l are constants to be determined, Pl are the Legendre polynomials, and the notation (r, θ, φ) means (radial, polar, azimuthal) coordinates, respectively.
Since Φ(r, θ) denotes the gravitational potential of an axisymmetric finite distribution of mass, the boundary condition limr→∞ Φ(r, θ) = 0 must be satisfied, thus the solution (2) takes the form
Following Vogt & Letelier (2005), in order to obtain a generalized Miyamoto-Nagai model and for the sake of simplicity, we shall consider terms up to l = 3 in (3), therefore, transforming to cylindrical coordinates (R, z) by means of the relations
and applying the additional transformation (Satoh, 1980),
with α and b two arbitrary parameters, the generalized potential takes the form∗
Once the potential has been specified, the mass-density distribution Σ can be calculated directly from Poisson equation,
while the circular velocity v of particles in the galactic plane, the epicyclic frequency k, and the vertical frequency ν of small oscillations about the equilibrium circular orbit, can be obobtained from the following expressions evaluated at z = 0 (Binney & Tremaine, 2011)
From (8-10), it is important to emphasize that a feasible model must satisfy the constraints set by the conditions v2 ≥ 0, k2 ≥ 0, and ν2 ≥ 0, where the last two inequalities are understood as stability conditions (Vogt & Letelier, 2005).
As is evident from the preceding paragraphs, the galactic models and its associated physical quantities are uniquely determined by the set of constants a, b, B0, B1, B2, and B3, which (taking a pragmatic approach) can be estimated from the observational data of the corresponding rotation curves, as we will discuss in detail in the next section.
Rotation curves fitting
The observational data were taken from Verheijen & Sancisi (2001) for three specific galaxies in the Ursa Major cluster: NGC 3726, NGC 3877 and NGC 4010. Following the procedure outlined in González et al. (2010), we take the galaxy radius Rd as the given by the largest tabulated value of the data.
Thus, introducing dimensionless variables and setting the nonlinear least square curve fitting method allows us to calculate the numerical values of the parameters for each particular galaxy. The resulting values of for the three galaxies under consideration, are give in Table 1.
In panels (a) of Figures 1, 2, and 3, we show the observational data (points) of the rotation curve with the corresponding velocity dispersions (error bars) as reported by Verheijen & Sancisi (2001) for NGC 3726, NGC 3877 and NGC 4010. The solid lines correspond to the analytical expressions (8) fitted to the rotation curves. As can be seen, in each case the model fits the observed data with good accuracy. Additionally, in panels (b) of Figures 1, 2, and 3, we plot the normalized mass-density distribution (7) at z = 0 for the three galaxies, as a function of the dimensionless radial coordinate . Here, we obtain a well-behaved mass-density function, showing a maximum value at the center that decreases to zero at the edge of the disk. On the other hand, in panels (c) of Figures 1, 2, and 3, we present four isodensity curves of the mass-density distribution (7) in the meridional plane , showing that each model corresponds to a very different mass distribution. Finally, from panels (d) and (e) of the same figures, it is noteworthy that in the three cases the stability conditions are fully satisfied.
Stellar Dynamics
It is a well-known fact that using rough estimates of the dimensions of typical stars and galaxies, the collision interval between stars is about 108 times longer than the average age for most galaxies (Binney & Tremaine, 2011). This implies that the star’s motion can be determined solely by the gravitational attraction of the galaxy and that collisions between stars are so rare that are irrelevant (Maoz, 2016). Therefore, as a first approximation, the orbital dynamics of a star in a given galaxy can be studied following the usual Lagrangian and Hamiltonian approaches for the motion of a test particle in the presence of an estimated gravitational potential.
The orbital motion of a test particle in an axisymmetric potential is governed by the Lagrangian
with (R, ø, z) the usual cylindrical coordinates. The generalized canonical momenta read as
and the Hamiltonian takes the form
With
Here, Lz = pφ =constant, denotes the conserved component of angular momentum about the z-axis.
From (13), the resulting Hamilton’s equations of motion can be expressed as
where Φ(R, z) is given by Eq. (6) and its respective parameters should be taken from Table 1.
Since the Hamiltonian is autonomous, H is an integral of motion
with h the energy of an orbit.
The existence of an analytic integral of motion reduces the phase space dimensionality, and hence the Poincaré surface of section is an appropriate and well-established method to analyze the dynamics of the system. Taking into account the axial symmetry associated to the system, it is customary to choose the equatorial plane z = 0 as the Poincaré plane in order to represent the surface of sections in the ( )-plane. The orbits were numerically integrated forward in time for 1000 units of time by using a Runge-Kutta-Fehlberg Method (RKF45), with this setting the numerical error related to the conservation of the energy is at most 10−14. In all cases we set z0 = pR0 = 0 and we scan the phase space with a large number of initial conditions for the radii R0, these three values allow us to determine the values of pz0 through the relation (19).
The transition from regularity to chaos (or viceversa) that takes place for the three considered galaxy models was inspected through the Poincaré sections in Figs. 4-9, by using different values of Lz (Figs. 4, 6 and 8) and h (Figs. 5, 7 and 9). It can be observed that the orbital motion exhibits a strong dependence on the angular momentum Lz and energy h of the test particle. In particular, from the surfaces of section presented in Figs. 4, 6, and 8, we may infer that there exists an increase in the regularity of the system for larger values of the angular momentum Lz , i.e. if there exists a chaotic sea the increase of Lz will fill the phase with KAM islands, while the opposite effect is observed for larger values of energy h (see Figs. 5, 7, and 9), where the KAM islands deform and shrink giving place to larger regions of chaos.
Concluding remarks
In the present paper, using the general solution to the Laplace equation, we have derived a generalized Miyamoto-Nagai potential. By means of the nonlinear least square fitting, the analytical velocity curves were adjusted to the observed ones of three specific spiral galaxies: NGC 3726, NGC 3877 and NGC 4010. The resulting analytical models were used to determine the mass-density distributions and the vertical and epicyclic frequencies, showing that unlike the results presented in González et al. (2010) for NGC 3877 and NGC 4010, our models satisfy the stability conditions for radial and vertical perturbations. Even though the set of models presented here should be considered as a rough approximation, the circular velocities were shown to fit very accurately to the observed rotation curves and in the three cases the stability conditions were fully satisfied. Here, it is important to note that contrary to the observed tendency in the Miyamoto-Nagay model, where the limit α → 0 reduces to the Plummer sphere, our models exhibit a tendency to an spherical mass distribution with increasing of the parameter a.
On the other hand, by using the Poincaré section method we have also studied the dynamics of the meridional orbits of stars in presence of the gravitational field of the galaxy models. From our results it may be inferred that there exists an increase in the regularity of the orbits for larger values of the angular momentum, while for larger values of energy the orbits tend to be more chaotic. Our toy models suggest that in the three galaxy models chaotic orbits are possible, however the chaotic behavior is very weak for the NGC 3877 model in comparison to NGC 3726 and NGC 4010. It should be noted that none of the studied models showed a fully chaotic phase space. Our results could have significant implications for the study of the dynamics and kinematics of these three specific galaxies, since the regular or chaotic behaviors could shed lights into the evolution and structure of these galaxies, i.e., in phase space, regular orbits are trapped in the vicinity of neighbor orbits, while chaotic orbits, by its own nature, will diverge exponentially in time from its neighbors by filling the phase space in an erratic manner.