Unveiling the Influence of Dark Matter in Axially Symmetric Galaxies

Abstract

We investigate the regular or chaotic nature of orbits of stars moving in the meridional plane (R, z) of an axially symmetric galactic model with a dense, massive spherical nucleus and a dark matter halo component. In particular, we study the influence of the fractional portion of the dark matter, by computing in each case the percentage of chaotic orbits, as well as the percentages of orbits of the main regular resonant families. In an attempt to distinguish between regular and chaotic motion, we use the fast Lyapunov indicator method to extensive samples of orbits obtained by integrating numerically the equations of motion as well as the variational equations. Furthermore, a technique which is based mainly on the field of spectral dynamics that utilises the Fourier transform of the time series of each coordinate is used for identifying the various families of regular orbits and also to recognise the secondary resonances that bifurcate from them. Two cases are studied in our work: (i) the case where we have a disk galaxy model and (ii) the case where our model represents an elliptical galaxy. A comparison with early related work is also made.

1 INTRODUCTION

From time immemorial, the movements of stars have attracted the interest of the first astronomers. Several decades have been past, since astronomers have turned their interest to the kinematics of stars in galaxies and in galaxy clusters in general. Undoubtedly, galaxies are the basic units of the Universe. An interesting task is to determine the amount of matter in galaxies and, consequently, the amount of matter of the Universe. In particular, there are two key events that are closely associated with this interesting endeavour. The first is the pioneer research of J. Oort who studied the motion of stars in the neighbourhood of the Sun and realised that in fact the stars were pulled by forces stronger than those that would be caused only by the visible matter (Oort 1924). On the other hand, the astronomer F. Zwicky discovered that the fast motion of galaxies in clusters cannot be justified only by the visible mass (Zwicky 1933). In other words, clusters of galaxies should be scrapped unless the their mass was much greater. The existence of dark matter is the most acceptable scenario to resolve these concerns. Things were clarified when Rubin & Ford (1970) presented observational evidence supporting the existence of dark matter in galaxies.

Knowing the rate ratio between luminous and dark matter in galaxies (disk or ellipticals), as well as the relations connecting the fundamental quantities which characterise both components is of great importance, since it allow us to understand how galaxies born and evolve and also how dark matter influences these procedures. Today, it is widely accepted that dark matter is the dominant element in galaxies, taking into account that the vast majority of the total matter of the Universe, 80% according to today measurements, is indeed dark. Dark matter seems to interact through gravity. Moreover, apart from gravitational, no other type of dark matter interaction has been observed. This strongly indicates that dark matter interactions should be very weak, probably much more weaker than the particle physics weak interactions. The basic proposed candidates corresponding to dark matter are: (i) neutrinos (hot dark matter), (ii) warm dark matter, (iii) cold dark matter (CDM) and finally (iv) weak interactive massive particles.

A strong indication for the presence of dark matter in galaxies is derived from their flattened rotation curves at large radii. Using Kepler's third law, we have $\Theta (R) = \sqrt{GM(R)/R}$ , where Θ is the rotational velocity at a radius R, G is the gravitational constant, while M(R) is the total mass within radius R. Conducting observations at large galactocentric distances, where no luminous galactic component was present, astronomers found that, instead of declining at the expected rate Θ(R)∝R ^−1/2, which is true if only $M = \text{const.}$ , the velocity curves Θ(R) appear to be flattened. However, $\Theta (R) = \text{const.}$ implies that M(R)∝R. This strongly suggests that the mass of galaxies continues to grow, even when there is no luminous component to account for this increase. Moreover, the profiles of the velocity curves indicate that the distribution of light does not match the distribution of mass. In other words, in each galaxy the mass-to-light ratio (M/L) increases with the radius R. The circular velocity needed for the construction of the rotation curve is obtained, as usually, by measuring the 21-cm emission line from neutral hydrogen H i (Rubin, Ford, & Thonnard 1980; Bosma 1981; Clemens 1985; Persic & Salucci 1995; Honma & Sofue 1997).

The existence of dark matter in elliptical galaxies has been confirmed using observational data derived either from hot gas or their X-ray emission (e.g. Loewenstein & White 1999; Humphrey et al. 2006; Johnson et al. 2009; Das et al. 2010) or even using strong lensing methods (e.g. Rusin & Kochanek 2005; Gavazzi et al. 2007; Koopmans et al. 2009; Faure et al. 2011). On the other hand, for the disk galaxies there is a plethora of scenarios describing the formation and the evolution of disk galaxies in correlation with dark matter (e.g. Dalcanton, Spergel, & Summers 1997; Firmani, Avila-Reese, & Hérnandez 1997; Avila-Reese, Firmani, & Hernández 1998; van de Bosch 1998; Avila-Reese & Firmani 2000). A large number of measurements of luminous and dark matter have become available over the last years (e.g. Cappellari et al. 2006, 2012; Wegner et al. 2012). Using these data in combination with dynamically modelling, we can reconstruct and therefore study the orbital structure of galaxies.

Today, we know that galaxies contain large amounts of dark matter and, therefore, a further study could provide important information about this invisible matter. A simple way to do this is to construct models of galaxies containing dark matter. The study of the dynamical properties of these models will provide interesting and useful information, which combined with data from observation will aid significantly in finding a solution for the problem of dark matter. The reader can find interesting information in the field of fitting mass models to the kinematics of disk and elliptical galaxies in a series of papers (see e.g. Barnes, Sellwood, & Kosowsky 2004; Chaktabarty 2007; de Blok et al. 2008; Gebhardt & Thomas 2009). Another important point that needs to be emphasised is the difference between the distribution of dark matter in galaxies and clusters. Observational data (e.g. Sahni 2004) suggest that dark matter increases as we move away from the centre to the outer parts of the galaxies. On the contrary, in galaxy clusters the dark matter distribution decreases with increasing distance from the galactic centre.

Even today, dark matter is still an open and controversial issue in astronomy. This is true because the standard model of cosmology seems to be incompatible with a large amount of data derived from extragalactic observations and modified gravity theories. The reader can find more interesting and detailed information regarding this issue in Kroupa (2012).

Taking into account all the above, there is no doubt that dark matter plays an important role in the dynamical behaviour of galaxies. On this basis, it seems of particular interest to build an analytical dynamical model describing the motion of stars both in disk and elliptical galaxies containing dark matter and then study how the presence and the amount of dark matter affects the regular or chaotic nature of orbits as well as the behaviour of the different families of orbits.

The present article is organised as follows. In Section 2, we present in detail the structure and the properties of our gravitational galactic model. In Section 3, we describe the computational methods we used in order to determine the character of orbits. In the following section, we investigate how the parameter corresponding to the fractional portion of the dark matter in galaxies influences the character of the orbits, in both disk and elliptical galaxy models. Our paper ends with Section 5, where the discussion and the conclusions of this research are presented.

2 DESCRIPTION OF THE GALACTIC MODEL

In order to study the dynamical properties of galaxies, astronomers often construct galactic models. A galactic model is usually a mathematical expression giving the potential or the mass density of the galaxy, as a function of the distance. The reader can find interesting models, describing motion in galaxies, in Binney & Tremaine (2008). Potential density pairs for galaxies were also presented by Vogt & Letelier (1996). Over the years, many galactic models have been proposed in order to model the orbital properties in axially symmetric systems. A simple yet realistic axisymmetric logarithmic potential was introduced in Binney (1981) for the description of galactic haloes at which the mass density drops like R ⁻² (see also Evans 1993). However, the most well-known model for CDM haloes is the flattened cuspy Navarro Frenk White (NFW) model (Navarro, Frenk, & White 1996, 1997), where the density at large radii falls like R ⁻¹. This model being self-consistent has a major advantage and that is why it is mainly used for conducting N-body simulations. Our model, on the other hand, is not self-consistent but simple and contrived, in order to give us the ability to study in detail the orbital behaviour of the galactic system. Nevertheless, contrived models can provide an insight into more realistic stellar systems, which unfortunately are very difficult to be studied if we take into account all the astrophysical aspects. Due to the fact that our gravitational model consists of Plummer-type potentials, at large galactocentric distances the mass density decreases following the R ⁻⁵ law.

In the present work, we shall investigate how the presence of the dark matter influences the character of the orbits in the meridional plane of an axially symmetric galaxy model with a spherical nucleus and a dark matter halo component. We shall use the usual cylindrical coordinates (R, φ, z), where z is the axis of symmetry.

The total potential V(R, z) in our model consists of three components: the main galaxy potential V _g, the central spherical component V _n, and the dark matter halo component V _h. The first one is represented by the new mass model introduced in Caranicolas (2012) and is given by

(1) \begin{equation} V_{\rm g}(R,z) = \frac{-G M_{\rm g}(1 + \delta )}{R_{\rm g}}, \end{equation}

where

(2) \begin{equation} R_{\rm g} = \sqrt{\left(1+\delta \right)\beta ^2 + R^2 + \left(\alpha + \sqrt{h^2 + \left(1+\delta ^2\right) z^2}\right)^2}. \end{equation}

Here, G is the gravitational constant, M _g is the mass of the galaxy, δ is the fractional portion of the dark matter in the galaxy, while β represents the core radius of the halo of the disk. The shape of the galaxy is controlled by the parameters α and h which correspond to the horizontal and vertical scale length of the galaxy, respectively. Therefore, this potential allow us to describe a variety of galaxy types from a flat disk galaxy when α, β≫h to an oblate elliptical galaxy when h≫α, β.

For the description of the spherically symmetric nucleus, we use a Plummer potential (e.g. Binney & Tremaine 2008):

(3) \begin{equation} V_{\rm n}(R,z) = \frac{-G M_{\rm n}}{\sqrt{R^2 + z^2 + c_{\rm n}^2}}, \end{equation}

where M _n and c _n are the mass and the scale length of the nucleus, respectively. This potential has been used successfully in the past in order to model and therefore interpret the effects of the central mass component in a galaxy (see e.g. Hasan & Norman 1990; Hasan, Pfenniger, & Norman 1993; Zotos 2012a). At this point, we must make it clear that Equation (3) is not intended to represent the potential of a black hole nor that of any other compact object, but just the potential of a dense and massive nucleus; thus, we do not include relativistic effects. The dark matter halo is modelled by a similar spherically symmetric potential

(4) \begin{equation} V_{\rm h}(R,z) = \frac{-G M_{\rm h}}{\sqrt{R^2 + z^2 + c_{\rm h}^2}}, \end{equation}

where, in this case, M _h and c _h are the mass and the scale length of the halo, respectively. The spherical shape of the dark halo is simply an assumption, due to the fact that galactic haloes may have a variety of shapes.

In this work, we use the well-known system of galactic units, where the unit of length is 1 kpc, the unit of mass is 2.325×10⁷ M_⊙ and the unit of time is 0.9778×10⁸ yr. The velocity unit is 10 km s⁻¹, the unit of angular momentum (per unit mass) is 10 km kpc s⁻¹, while G is equal to unity. Finally, the energy unit (per unit mass) is 100 km² s⁻². In these units, the values of the involved parameters are M _g = 8200, M_n = 400, c _n = 0.25, M _h = 7560, and c _h = 15. For the disk model, we choose β = 6, α = 3, and h = 0.3, while for the elliptical model, we have set β = 2, α = 0.5, and h = 11.5. The fractional portion of dark matter δ, on the other hand, is treated as a parameter and its value varies in the interval 0 ≤ δ ≤ 0.5.

It is well known that in disk galaxies, the circular velocity in the galactic plane z = 0,

(5) \begin{equation} \theta (R) = \sqrt{R\left|\frac{\partial V}{\partial R}\right|_{z=0}}, \end{equation}

is a very important physical quantity. A plot of θ(R) for both disk and elliptical galactic models when δ = 0.1 is presented in Figures 1(a, b), as a black curve. Moreover, in the same plot, the green curve is the contribution from the luminous matter, while the red line corresponds to the contribution from the dark matter. It is seen that in both cases the component of the rotational curve generated by the dark matter remains flat for large galactocentric distances, while, on the other hand, the luminous matter component continues to decline with increasing distance from the galactic centre. We also observe the characteristic local minimum of the rotation curve due to the massive nucleus, which appears at small values of R, when fitting observed data to a galactic model (e.g. Irrgang et al. 2013; Gómez et al. 2010).

Figure 1. A plot of the rotation curve in our (a) disk and (b) elliptical galactic models. We can distinguish the total circular velocity (black) and the contributions from luminous matter (green) and that of the dark matter (red).

At this point, we must clarify that the mass density in our new galaxy model obtains negative values when the distance from the centre of the galaxy described by the model exceeds a minimum distance d _min, which strongly depends on the parameter δ. Figures 2(a, b) show a plot of d _min versus δ for the both disk and elliptical galaxy models. We see that even when δ = 0.5, the first indication of negative density occurs only when d _min>40 kpc, that is almost at the theoretical boundaries of a real galaxy. Here, we must point out that our gravitational potential is truncated at R _max = 30 kpc for both reasons: (i) otherwise, the total mass of the galaxy modelled by this potential would be infinite, which is obviously not physical, and (ii) to avoid the existence of negative values of density.

Figure 2. The evolution of the minimum distance d _min where negative density appears as a function of the parameter δ for (a) the disk galaxy model and (b) the elliptical galaxy model.

Taking into account that the total potential V(R, z) is axisymmetric, the z-component of the angular momentum L_z is conserved. With this restriction, orbits can be described by means of the effective potential

(6) \begin{equation} V_{\rm eff}(R,z) = V(R,z) + \frac{L_z^2}{2R^2}. \end{equation}

The L ² _z /(2R ²) term represents a centrifugal barrier; only orbits with small L_z are allowed to pass near the axis of symmetry. The three-dimensional (3D) motion is thus effectively reduced to a 2D motion in the meridional plane (R, z), which rotates non-uniformly around the axis of symmetry according to

(7) \begin{equation} \dot{\phi } = \frac{L_z}{R^2}, \end{equation}

where of course the dot indicates derivative with respect to time.

The equations of motion on the meridional plane are

(8) \begin{eqnarray} \ddot{R} &=& - \frac{\partial V_{\rm eff}}{\partial R}, \nonumber \\ \ddot{z} &=& - \frac{\partial V_{\rm eff}}{\partial z}, \end{eqnarray}

while the equations describing the evolution of a deviation vector $\delta {{\bm w}} = (\delta R, \delta z, \delta \dot{R}, \delta \dot{z})$ which joins the corresponding phase space points of two initially nearby orbits, needed for the calculation of the standard chaos indicators (the FLI in our case), are given by the variational equations:

(9) \begin{eqnarray} \dot{(\delta R)} &=& \delta \dot{R}, \nonumber \\ \dot{(\delta z)} &=& \delta \dot{z}, \nonumber \\ (\dot{\delta \dot{R}}) &=& - \frac{\partial ^2 V_{\rm eff}}{\partial R^2} \delta R - \frac{\partial ^2 V_{\rm eff}}{\partial R \partial z}\delta z, \nonumber \\ (\dot{\delta \dot{z}}) &=& - \frac{\partial ^2 V_{\rm eff}}{\partial z \partial R} \delta R - \frac{\partial ^2 V_{\rm eff}}{\partial z^2}\delta z. \end{eqnarray}

Consequently, the corresponding Hamiltonian to the effective potential given in Equation (6) can be written as

(10) \begin{equation} H = \frac{1}{2} \left(\dot{R}^2 + \dot{z}^2 \right) + V_{\rm eff}(R,z) = E, \end{equation}

where $\dot{R}$ and $\dot{z}$ are momenta per unit mass, conjugate to R and z respectively, while E is the numerical value of the Hamiltonian, which is conserved. Therefore, an orbit is restricted to the area in the meridional plane satisfying E ≥ V _eff.

3 COMPUTATIONAL METHODS

In order to obtain the mass profiles of galaxies, we have to construct dynamical models describing the main properties of the galaxies. These models can be generated by deploying two main techniques: (i) using superposition of libraries of orbits (e.g. Gebhardt et al. 2003; Thomas et al. 2004; Cappellari et al. 2006) or (ii) using distributions functions (e.g. Dejonghe et al. 1996; Gerhard et al. 1998; Kronawitter et al. 2000). In the literature, there are also other more specialised dynamical models combining kinematic and photometric data. For instance, axially symmetric Schwarzschild models were used by Bridges et al. (2006), while Hwang et al. (2008) used Jeans models in order to fit observational data in the X-ray potential introduced by Humphrey et al. (2006). Furthermore, axisymmetric Schwarzschild models were also used by Shen & Gebhardt (2010) to fit data derived from the Hubble Space Telescope.

In the recent years, the Schwarzschild superposition method (Schwarzschild 1979) has been heavily utilised by several researchers (e.g. Rix et al. 1997; Gebhardt et al. 2003; Thomas et al. 2004, 2005; Valluri, Merritt, & Emsellem 2004; Krajnović et al. 2005) in order to study dark matter distributions in elliptical galaxies; therefore, we deem it necessary to recall and describe briefly the basic points of this interesting method. Initially, N closed cells define the configuration space, while K orbits extracted from a given mass distribution construct the phase space. Then, integrating numerically the equations of motion, we calculate the amount of time spent by each particular orbit in every cell. Thus, the mass of each cell is directly proportional to the total sum of the stay times of orbits in every cell. Using this technique, we manage to compute the unknown weights of the orbits, assuming they are not negative.

In our study, we want to know whether an orbit is regular or chaotic. Several chaos indicators are available in the literature; we chose the FLI method. The FLI (Froeschlé, Gonczi, & Lega 1997; Lega & Froeschlé 2001) has been proved a very fast, reliable and effective tool, which is defined as

(11) \begin{equation} {\rm FLI}(t) = \log \Vert {{\bm w}(t)} \Vert , \ \ \ \ t \le t_{\rm max}, \end{equation}

where ${{\bm w}}(t)$ is a deviation vector. For distinguishing between regular and chaotic motion, we need to compute the FLI for a relatively short time interval of numerical integration t _max. In particular, we track simultaneously the time evolution of the orbit itself as well as the deviation vector ${{\bm w}}(t)$ in order to compute the FLI. The variational equations (9), as usual, are used for the evolution of the deviation vector. The main advantage of the FLI method is that only one deviation vector is required to be computed, while in the case of other dynamical indicators such as SALI (Skokos et al. 2004) or GALIs (Skokos, Bountis, & Antonopoulos 2007) more than one deviation vectors are needed. Therefore, by using the FLI method we need considerably less computation time for integrating and classifying massive sets of initial conditions of orbits.

The particular time evolution of the FLI allows us to distinguish fast and safely between regular and chaotic motion as follows: when an orbit is regular the FLI exhibits a linear increase, while, on the other hand, in the case of chaotic orbits the FLI evolution is super-exponential. The time evolution of a regular (R) and a chaotic (C) orbit for a time period of 10⁴ time units is presented in Figure 3. We observe that both regular and chaotic orbits exhibit an increase in the value of FLI but with a complete different rate. Unfortunately, this qualitative criterion is applicable only when someone wants to check the character of individual orbits by plotting and then inspecting by eye the evolution of the FLI. Nevertheless, we can easily overcome this drawback by establishing a numerical threshold value, in order to quantify the FLI. After conducting extensive numerical experiments in different types of dynamical systems, we conclude that a safe threshold value for the FLI, taking into account the total integration time of 10⁴ time units, is 10. The horizontal, blue, dashed line in Figure 3 corresponds to that threshold value which separates regular from chaotic motion. In order to decide whether an orbit is regular or chaotic, we may use the usual method, according to which we check after a certain and predefined time interval of numerical integration whether the value of the FLI has become greater than the established threshold value. Therefore, if FLI ≥10 the orbit is chaotic, while if FLI <10 the orbit is regular.

Figure 3. Evolution of the FLI of a regular orbit (green colour, R), a sticky orbit (orange colour, S), and a chaotic orbit (red colour, C) in our model for a time period of 10⁴ time units. The horizontal, blue, dashed line corresponds to the threshold value of 10 which separates regular from chaotic motion. The chaotic orbit needs only about 130 time units in order to cross the threshold value, while, on the other hand, the sticky orbit requires a vast integration time of about 4400 time units so as to reveal its chaotic nature.

In order to investigate the orbital properties (chaoticity or regularity) of the dynamical system, we need to establish some samples of initial conditions of orbits. The best approach, undoubtedly, would have been to extract these samples of orbits from the distribution function of the model. Unfortunately, this is not available so, we followed another course of action. For determining the chaoticity of our models, we chose, for each set of values of the parameters of the potential, a dense grid of initial conditions in the $(R,\dot{R})$ phase plane, regularly distributed in the area allowed by the value of the energy E. The points of the grid were separated 0.1 units in R and 0.5 units in the $\dot{R}$ direction. For each initial condition, we integrated the equations of motion (8) as well as the variational equations (9) with a double-precision Bulirsch–Stoer algorithm (e.g. Press et al. 1992). In all cases, the energy integral (Equation (10)) was conserved better than one part in 10⁻¹⁰, although for most orbits, it was better than one part in 10⁻¹¹.

Each orbit was integrated numerically for a time interval of 10⁴ time units (10 billion years), which corresponds to a time span of the order of hundreds of orbital periods but of the order of one Hubble time. The particular choice of the total integration time is an element of paramount importance, especially in the case of the so-called ‘sticky orbits’ (i.e. chaotic orbits that behave as regular ones during long periods of time). A sticky orbit could be easily misclassified as regular by any chaos indicator ¹ , if the total integration interval is too small, so that the orbit does not have enough time in order to reveal its true chaotic character. A characteristic example of a sticky orbit (S) in our galactic system can be seen in Figure 3, where we observe that the chaotic character of the particular sticky orbit is revealed after a considerably long integration time of about 4400 time units. Thus, all the sets of orbits of the grids of the initial conditions were integrated, as already stated, for 10⁴ time units, thus avoiding sticky orbits with a stickiness at least of the order of a Hubble time. All the sticky orbits which do not show any signs of chaoticity for 10⁴ time units are counted as regular ones, because vast sticky periods are completely out of scope of our research.

A first step towards the understanding of the overall behaviour of our galactic system is knowing the regular or chaotic nature of orbits. Of particular interest, however, is also the distribution of regular orbits into different families. Therefore, once the orbits have been characterised as regular or chaotic, we then further classified the regular orbits into different families, by using the frequency analysis of Carpintero & Aguilar (1998) and Muzzio, Carpintero, & Wachlin (2005). Initially, Binney & Spergel (1982, 1984) proposed a technique, dubbed spectral dynamics, for this particular purpose. Later on, this method has been extended and improved by Carpintero & Aguilar (1998) and Šidlichovský & Nesvorný (1996). In a recent work (Zotos & Carpintero 2013), the algorithm was refined even further, so it can be used to classify orbits in the meridional plane. In general terms, this method calculates the Fourier transform of the coordinates of an orbit, identifies its peaks, extracts the corresponding frequencies and searches for the fundamental frequencies and their possible resonances. Thus, we can easily identify the various families of regular orbits and also recognise the secondary resonances that bifurcate from them.

Before closing this section, we would like to make a short note about the nomenclature of orbits. All the orbits of an axisymmetric potential are in fact 3D loop orbits, i.e. orbits that rotate around the axis of symmetry always in the same direction. However, in dealing with the meridional plane, the rotational motion is lost, so the path that the orbit follows onto this plane can take any shape, depending on the nature of the orbit. We will call an orbit according to its behaviour in the meridional plane. Thus, if, for example, an orbit is a rosette lying in the equatorial plane of the axisymmetric potential, it will be a linear orbit in the meridional plane, etc.

4 ORBIT CLASSIFICATION

In this section, we will numerically integrate several sets of orbits in an attempt to distinguish the regular or chaotic nature of motion. We use the initial conditions mentioned in Section 3 in order to construct the respective grids of initial conditions, taking always values inside the zero velocity curve (ZVC) defined by

(12) \begin{equation} \frac{1}{2} \dot{R}^2 + V_{\rm eff}(R,0) = E. \end{equation}

In all cases, the value of the angular momentum of the orbits is L_z = 15. We chose for both disk and elliptical galaxy models such energy levels which correspond to R _max≃10 kpc, where R _max is the maximum possible value of R on the $(R,\dot{R})$ phase plane. Once the values of the parameters were chosen, we computed a set of initial conditions as described in Section 3 and integrated the corresponding orbits calculating the value of the FLI and then classifying the regular orbits into different families.

4.1 Disk galaxy model

For the disk galaxy models, we choose the energy level E = −1200, which is kept constant. Our investigation reveals that in our disk galaxy model there are six main types of orbits: (a) box orbits, (b) 1:1 linear orbits, (c) 2:1 banana-type orbits, (d) 3:2 resonant orbits, (e) 4:3 resonant orbits, and (f) chaotic orbits. Note that every resonance n:m is expressed in such a way that m is equal to the total number of islands of invariant curves produced in the $(R,\dot{R})$ phase plane by the corresponding orbit. In Figures 4(a–f), we present an example of each of the five basic types of regular orbits, plus an example of a chaotic one. In all cases, we set δ = 0.5. The orbits shown in Figures 4(a, f) were computed until t = 100 time units, while all the parent periodic orbits were computed until one period has completed. The black thick curve circumscribing each orbit is the limiting curve in the meridional plane (R, z) defined as V _eff(R, z) = E. Table 1 shows the types and initial conditions for each of the depicted orbits; for the resonant cases, the initial conditions and the period T _per correspond to the parent periodic orbits.

Table 1. Types and initial conditions of the disk galaxy model orbits shown in Figures 4(a–f). In all cases, z ₀ = 0 and $\dot{z_0}$ is found from the energy integral, Equation (10), while T _per is the period of the resonant parent periodic orbits.

Figure 4. Orbit collection of the six basic types of orbits in the disk galaxy model: (a) box orbit, (b) 2:1 banana-type orbit, (c) 1:1 linear orbit, (d) 3:2 boxlet orbit, (e) 4:3 boxlet orbit, and (f) chaotic orbit.

To study how the fractional portion of the dark matter δ influences the level of chaos, we let it vary while fixing all the other parameters of our disk galaxy model. As already said, we fixed the values of all the other parameters and integrate orbits in the meridional plane for the set δ = {0, 0.1, 0.2,. . ., 0.5}. In all cases, the energy was set to −1200 and the angular momentum of the orbits was L_z = 15. Once the values of the parameters were chosen, we computed a set of initial conditions as described in Section 3 and integrated the corresponding orbits computing the FLI of the orbits and then classifying regular orbits into different families.

In Figures 5(a–f) we present six grids of orbits that we have classified for different values of the fractional portion of the dark matter δ. Here, we can identify all the different regular families by the corresponding sets of islands which are formed in the phase plane. In particular, we see the five main families already mentioned: (i) 2:1 banana-type orbits surrounding the central periodic point; (ii) box orbits are situated mainly outside of the 2:1 resonant orbits; (iii) 1:1 open linear orbits form the double set of elongated islands in the outer parts of the phase plane; (iv) 3:2 resonant orbits form the double set of islands above the box orbits; and (v) 4:3 resonant orbits correspond to the outer triple set of islands shown in the phase plane. Apart from the regions of regular motion, we observe the presence of a unified chaotic sea which embraces all the islands of stability. The outermost black thick curve is the ZVC defined by Equation (12).

Figure 5. Orbital structure of the $(R,\dot{R})$ phase plane of the disk galaxy model for different values of the fractional portion of the dark matter δ.

Figure 6 shows the resulting percentages of the chaotic orbits and of the main families of regular orbits as δ varies. It can be seen that there is a strong correlation between the percentage of most orbits and the value of δ. As the portion of the dark matter increases, there is a gradual decrease in the percentage of chaotic orbit, although this tendency is reversed in models with significant amount of dark matter (δ>0.3). In particular, we observe that chaotic motion is always the dominant type of motion and when δ>0.3 the amount of chaotic orbits grows at the expense of box orbits and 1:1 linear orbits. The meridional 2:1 banana-type orbits, on the other hand, are almost unperturbed by the shifting of the portion of dark matter. Moreover, the 4:3 resonant orbits exhibit a constant increase, while the percentage of the 3:2 resonant orbits remains at very low values. From the diagram shown in Figure 6, one may conclude that the fractional portion of the dark matter affects mostly the 1:1, 4:3 resonant orbits and chaotic orbits in disk galaxy models.

Figure 6. Evolution of the percentages of the different kinds of orbits in our disk galaxy model, when varying the fractional portion of the dark matter δ.

Of particular interest is to investigate how the variation in the fractional portion of the dark matter influences the position of the periodic points of the different families of periodic orbits shown in the grids of Figure 5. For this purpose, we use the theory of periodic orbits (Meyer & Hall 1992) and the algorithm developed and applied in Zotos (2013b). In Figures 7(a–d) we present the evolution of the starting position of the parent periodic orbits of the four basic families of resonant orbits. The evolution of the 2:1 and 4:3 families shown in Figures 7(a) and (b) respectively is two dimensional since the starting position (R ₀, 0) of both families lies on the R axis. On the contrary, studying the evolution of the 1:1 and 3:2 families of periodic orbits is indeed a real challenge due to the peculiar nature of their starting position $(R_0,\dot{R_0})$ . In order to visualise the evolution of these families, we need 3D plots such as those presented in Figures 7(c) and (d), taking into account the simultaneous relocation of R ₀ and $\dot{R_0}$ . The stability of the periodic orbits can be obtained from the elements of the monodromy matrix X(t) as follows:

(13) \begin{equation} K = {\rm Tr} \left[X(t)\right] - 2, \end{equation}

where Tr stands for the trace of the matrix and K is called the stability index. For each set of values of δ, we first located, by means of an iterative process, the position of the parent periodic orbits. Then, using these initial conditions, we integrated the variational equations in order to obtain the matrix X, with which we computed the index K. Our numerical calculations indicate that in the disk galaxy models, all the different families of periodic orbits remain stable throughout the entire range of the values of δ.

Figure 7. Evolution of the starting position $(R_0,\dot{R_0})$ of the periodic orbits as a function of the fractional portion of the dark matter δ. (a) 2:1 resonant family, (b) 4:3 resonant family, (c) 1:1 resonant family, and (d) 3:2 resonant family.

4.2 Elliptical galaxy model

In the case of the elliptical galaxy model, we choose the energy level E = −1100, which is kept constant. Our numerical investigation shows that in our elliptical galaxy model, there are seven main types of orbits: (a) box orbits, (b) 1:1 linear orbits, (c) 2:1 banana-type orbits, (d) 3:2 resonant orbits, (e) 4:3 resonant orbits, (f) 5:3 resonant orbits, (g) 8:5 resonant orbits, and (h) chaotic orbits. It is worth noting that the basic resonant families, that is the 2:1, 1:1, 3:2, and 4:3, are common in both disk and elliptical galaxy models. However, in the case of the elliptical galaxy, additional secondary resonances (i.e. 5:3 and 8:5) appear. Again, every resonance n:m is expressed in such a way that m is equal to the total number of islands of invariant curves produced in the $(R,\dot{R})$ phase plane by the corresponding orbit. In Figures 8(a–h), we present an example of each of the seven basic types of regular orbits, plus an example of a chaotic one. In all cases, we set δ = 0.5. The orbits shown in Figures 8(a) and (f) were computed until t = 100 time units, while all the parent periodic orbits were computed until one period has completed. Table 2 shows the types and the initial conditions for each of the depicted orbits; for the resonant cases, the initial conditions and the period T _per correspond to the parent periodic orbits.

Table 2. Types and initial conditions of the elliptical galaxy model orbits shown in Figures 8(a–h). In all cases, z ₀ = 0 and $\dot{z_0}$ is found from the energy integral, Equation (10), while T _per is the period of the resonant parent periodic orbits.

Figure 8. Orbit collection of the eight basic types of orbits in the elliptical galaxy model: (a) box orbit, (b) 2:1 banana-type orbit, (c) 1:1 linear orbit, (d) 3:2 boxlet orbit, (e) 4:3 boxlet orbit, (f) 5:3 boxlet orbit, (g) 8:5 boxlet orbit, and (h) chaotic orbit.

In order to study how the fractional portion of the dark matter δ influences the level of chaos, we let it vary while fixing all the other parameters in our elliptical galaxy model. As already said, we fixed the values of all the other parameters and integrate orbits in the meridional plane for the set δ = {0, 0.1, 0.2, . . ., 0.5}. In all cases, the energy was set to −1100 and the angular momentum of the orbits was L_z = 15. Once the values of the parameters were chosen, we computed a set of initial conditions as described in Section 3 and integrated the corresponding orbits computing the FLI of the orbits and then classifying regular orbits into different families.

Six grids of initial conditions $(R_0,\dot{R_0})$ that we have classified for different values of the fractional portion of the dark matter δ are shown in Figures 9(a–f). By inspecting these grids, we can identify all the different regular families by the corresponding sets of islands which are produced in the phase plane. In particular, we see the seven main families of regular orbits already mentioned: (i) 2:1 banana-type orbits correspond to the central periodic point; (ii) box orbits situated mainly outside of the 2:1 resonant orbits; (iii) 1:1 open linear orbits form the double set of elongated islands in the outer parts of the phase plane; (iv) 3:2 resonant orbits form the double set of islands; (v) 4:3 resonant orbits correspond to the outer triple set of islands shown in the phase plane; (vi) 5:3 resonant orbits form the set of the three small islands inside the region of box orbits; and (vii) 8:5 resonant orbits produce the set of five islands. It is evident that the structure of the phase plane in the elliptical galaxy models differs greatly from that of the disk models. We observe that when the amount of dark matter in the elliptical galaxy is low, almost the entire phase plane is covered by different types of regular orbits. On the other hand, chaotic motion appears, mainly at the outer parts of the phase plane, only when the galaxy possesses a significant amount of dark matter (δ ≥ 0.4).

Figure 9. Orbital structure of the $(R,\dot{R})$ phase plane of the elliptical galaxy model for different values of the fractional portion of the dark matter δ.

In Figure 10, we present the resulting percentages of the chaotic orbits and of also the main families of regular orbits as δ varies. It can be seen that the motion of stars in elliptical galaxies is almost entirely regular, the box orbits being the all-dominant type. The percentage of box orbits is however reduced as the portion of the dark matter is increased, although they always remain the most populated family. It is also seen that the percentages of the 2:1 banana-type orbits exhibit a minor decrease with increasing δ. On the other hand, the chaotic orbits start to grow rapidly as soon as the galaxy contains significant amount of dark matter (δ>0.3). Moreover, all the other resonant families of orbits are immune to changes of the portion of the dark matter, since their percentages remain almost unperturbed and at very low values (less than 5%). From Figure 10 one may conclude that dark matter in elliptical galaxies mostly affects box and chaotic orbits. In fact, a portion of box orbits turns into chaotic as the galaxy gains more dark matter.

Figure 10. Evolution of the percentages of the different kinds of orbits in our elliptical galaxy model, when varying the fractional portion of the dark matter δ.

We close this section by presenting how the variation in the fractional portion of the dark matter influences the position of the periodic points of the different families of periodic orbits shown in the grids of Figures 9(a–f). We follow the same method as used previously in the case of the disk galaxy. In Figures 11(a–f), the evolution of the starting position of the parent periodic orbits of the six basic families of resonant orbits is given. Once again, the evolution of the 2:1, 4:3, and 8:5 families shown in Figures 11(a–c) is two dimensional since the starting position (R ₀, 0) of these families lies on the R axis. On the contrary, the evolution of the 1:1, 3:2, and 5:3 families is shown in the 3D plots in Figures 11(d–f), thus following simultaneously the relocation of R ₀ and $\dot{R_0}$ as the fractional portion of the dark matter δ varies. Our numerical calculations suggest that all the computed resonant periodic orbits were found to be stable. Furthermore, we should point out that several families of periodic orbits in the case of the elliptical galaxy do not cover the entire range of the values of δ

Figure 11. Evolution of the starting position $(R_0,\dot{R_0})$ of the periodic orbits as a function of the fractional portion of the dark matter δ. (a) 2:1 resonant family, (b) 4:3 resonant family, (c) 8:5 resonant family, (d) 1:1 resonant family, (e) 3:2 resonant family, and (f) 5:3 resonant family.

5 DISCUSSION AND CONCLUSIONS

There is no doubt that the determination of the nature of dark matter is one of the most interesting and challenging open problems that scientists try to solve. In the present paper, we used the analytic, axisymmetric, mass model which was introduced in Caranicolas (2012) and embraces the general features of a dense, massive nucleus and a spherical dark matter halo. We made this choice because observations show that the assumption of a spherical halo seems to be reasonable. On the other hand, non-spherical and triaxial haloes are also possible in some galaxies (see e.g. Gentile et al. 2004; Trachternach et al. 2008; Vera-Ciro et al. 2011). A galaxy with a dark matter halo is undoubtedly a very complex entity and, therefore, we need to assume some necessary simplifications and assumptions in order to be able to study mathematically the orbital behaviour of such a complicated stellar system. For this purpose, our model is simple and contrived, in order to give us the ability to study different aspects of the dynamical model. Nevertheless, contrived models can provide an insight into more realistic stellar systems, which unfortunately are very difficult to be studied, if we take into account all the astrophysical aspects. On the other hand, self-consistent models are mainly used when conducting N-body simulations. However, this is entirely out of the scope of the present paper. Once again, we have to point out that the simplicity of our model is necessary; otherwise, it would be extremely difficult, or even impossible, to apply the extensive and detailed dynamical study presented in this study. Similar gravitational models with the same limitations and assumptions were used successfully several times in the past in order to investigate the orbital structure in much more complicated galactic systems (see e.g. Zotos 2012b, 2013a).

In this work, we investigated how influential is the parameter corresponding to the portion of the dark matter δ on the level of chaos and also on the distribution of regular families among its orbits in both disk and elliptical galaxy models. The main results of our research can be summarised as follows.

1. 1. In disk galaxy models, the fractional portion of the dark matter affects mostly the 1:1 and 4:3 resonant orbits and the chaotic orbits, while the effect on all the other resonant families is very weak compared with them. In particular, chaotic motion is always the prevailing type of motion but when the amount of dark matter is high enough, the amount of chaotic orbits grows at the expense of box orbits and 1:1 linear orbits.

2. 2. That portion of dark matter in elliptical galaxy models influences mainly the box and the chaotic orbits. In fact, box orbits are the dominant family when the amount of dark matter is low, but the percentage of chaotic orbits quickly grows as the dark matter is being accumulated, by collapsing the percentage of box orbits, although they always remain the most populated family.

3. 3. The percentage of the observed chaos in disk galaxy models with dark matter is significantly larger than that in elliptical galaxy models. This result coincides with similar conclusions obtained using different types of dynamical models in order to model disk and elliptical galaxies (see e.g. Papadopoulos & Caranicolas 2005; Zotos 2011).

The outcomes of the present research are considered as a promising first step in the task of exploring the orbital structure in both disk and elliptical galaxy models containing dark matter. Taking into account that our results are encouraging, it is in our future plans to modify properly our dynamical model in order to expand our study in three dimensions.