2.1. The Smaller Alignment Index (SALI)

In order to define SALI, let us consider a Hamiltonian flow of N degrees of freedom, an orbit in the 2N-dimensional phase space with initial condition x0=x10,⋯,x2N0 and two normalized deviation vectors w1^0, w2^0 from the initial condition x0.

We define

where the quantities w^1t-w^2t and w^1t+w^2t are called Parallel Alignment Index and Antiparallel Alignment Index, respectively.

It is evident that SALIt∈0,2 and when SALI=0 the two normalized vectors have the same direction, being equal or opposite.

The SALI value is a very useful tool for detecting chaos in Hamiltonian systems. Chaotic or regular motions are easily distinguishable applying the SALI method. In the case of chaotic orbits, the deviation vectors w^1t and w^2t align in the direction defined by the Maximum Lyapunov Exponent (MLE) and SALIt falls exponentially to zero:

with L ₁ and L ₂ the two largest Lyapunov exponents.

Furthermore, for regular motions the orbits develop on a phase space torus and eventually the vectors w^1t and w^2t fall in the torus tangent space, following a t-1 time dependence. In this case, SALI oscillates at nonzero values:

We have a clear distinction between ordered and chaotic behaviors using the SALI method in Hamiltonian systems. For mathematical SALI details, we recommend reading the papers ^{Skokos (2001)}; ^{Skokos et al. (2002a}, ²⁰⁰³, ²⁰⁰⁴⁾.

2.2. Mathematical Modeling of the Gravitational Potential

In this investigation, we used the gravitational potential divided into three basic components: bulge, disk and bar, according to the following equation:

Each component of equation (4) was mathematically modeled according to a classical potential: the ^{Plummer potential was used for the bulge (Plummer 1911}), Miyamoto-Nagai’s for the disk (^{Miyamoto & Nagai 1975}) and ^{Ferrers’ for the bar (Ferrers 1877}). This way of representing the total gravitational potential has been extensively used in many articles, such as ^{Patsis (2002)}; ^{Manos & Athanassoula (2011)}; ^{Skokos et al. (2002c},^d); ^{Patsis et al. (2002}, ²⁰⁰³⁾ and ^{Caritá et al. (2017)}.

The Plummer potential is written as:

where ϵ is the scale-length of the bulge, M _S is its total mass, and G is the gravitational constant.

The Miyamoto-Nagai’s potential is written as:

where M _D is the total disk mass, A and B are its horizontal and vertical scale-lengths, and G is the gravitational constant.

The Ferrers’ potential is written as:

where m2u=x2a2+u+y2b2+u+z2c2+u, Δ2u=a2+ub2+uc2+u, (λ) = 1 for the region outside the bar (m ≥ 1) and λ = 0 for the region inside the bar (m < 1).

In this last potential, the density is given by

where the central density is ρc=10532πGMBabc, MB is the bar mass and m2=x2a2+y2b2+z2c2, where a > b > c > 0 are the semi-axes of the ellipsoid which represents the bar.

In order to implement this bar model computationally, we used the analytical version given by Dr. Pfenniger, who kindly provided us with his Fortran 77 routine of the Ferrers potential. In this routine, the polynomial form of the Ferrers potential (^{Pfenniger 1984}; ^{Caritá et al. 2017}) was used.

In the course of this work, the SALI method was applied to study stellar orbits in a gravitational potential of barred galaxies, since the motion of a test particle in a rotating 3-dimensional model of a barred galaxy is given by the Hamiltonian:

where the bar rotates around the z-axis; x and y respectively are the major and minor galactic bar axes, Φ _Total is the total gravitational potential given by equation (4) and Ω _b is the bar pattern speed.

We emphasize that in order to follow the evolution of the orbits and that of their deviation vectors (for SALI computation), it is necessary to know the equations of motion and the variational equations linked to the Hamiltonian (9). The corresponding motion and variational equations can be checked in ^{Manos & Machado (2014)}.

To study orbit stability in models with 2 degrees of freedom, in our calculations, z = 0 and p _z = 0 were adopted in the Hamiltonian shown in equation (9).

2.3. Implementation and Computation Using the LP-VIcode

To perform the orbital integrations and SALI calculation, the LP-VIcode (^{Carpintero et al. 2014}) was employed, which is freely available at http://lp-vicode.fcaglp.unlp.edu.ar/.

LP-VIcode is an operational code in Fortran 77 that efficiently calculates 10 chaos indicators for dynamic systems, including SALI. The program reads the initial conditions for one or more orbits, integrates them (using a Bulirsch-Stoer integrator), and calculates the equations of the chosen chaos indicators. More details about the structure and operation of the LP-VIcode can be found in ^{Carpintero et al. (2013)} and ^{Carpintero et al. (2014)}.

In order to integrate orbits using the program and to study their stability, the user must provide the potential expressions and the motion and variational equations. That is, there is an external routine where these equations must be written in Fortran 77 by the user.

Two actions performed in the LP-VIcode implementation and adaptation stage should be highlighted in this section: the first is the adjustments that were made in order to implement a rotating coordinate system in the code (since the general equations (of motion and variational) present in the original main program consider only a static reference frame; it is known that in order to model a rotating galactic bar it is necessary to consider a coordinate system that rotates along with the bar. To do this we inserted the pattern speed Ω _b in the motion and variational equations in the LP-VIcode main program transforming the original equations into the form which ^{Manos & Athanassoula (2011)} displayed in their text. The second important action was bar growth. As previously explained, in this work we modeled the emergence of a galactic bar. The idea was to create a system that started as a barless galaxy and that would become a barred galaxy later on, where the bar grew over time. In order to model this evolution of the bar potential, a linear function of time was assumed for the mass of the bar in the Ferrers potential (this was done in an external routine of the LP-VIcode, the same routine where the user provides the potential and the motion and variational equations).

With these two actions we were able to use the LP- VIcode to study the dynamics of a barred galaxy, with the mass of the bar growing over time, and remaining constantly rotating around the z-axis.

3.1. Models and Parameters

We studied three models X, Y and Z whose parameters are shown in Table 1. Although it is not necessary to understand the origin of these parameters for our study, the reader can find a brief description of the procedures adopted in Appendix A, where we explain that the inspiration for these parameters came from the galaxy NGC 936, and where we describe how we computed the parameters based on the works of ^{Kent & Glaudell (1989)} and ^{Merrifield & Kuijken (1995)}. Each of the models was divided into two more specific models, where the rate of bar growth was varied (in one of them, the bar totally evolves with 5 turns around itself and in the other with 10 turns) generating in total six models: X₅,X₁₀,Y₅,Y₁₀,Z₅ and Z₁₀.

In Table 1, and all along this paper, the model system of units was defined considering the gravitational constant G = 1. We adopted 1kpc for length, 10³ kms⁻¹ for velocity, 10³ kms⁻¹ kpc⁻¹ for pattern speed, 1 Myr for time, and 2 × 10¹¹ M _⊙ for mass.

The total mass G(M _S +M _D +M _B ) was always equal to 1. For the energy, the unit is 10⁶ km²s⁻². The integration time was 10000 Myr.

For all models (X, Y and Z) the masses of the bulge, disk and bar components, as well as their other parameters do not change. Therefore, the difference between these three models is basically the galactic bar pattern speed; Model X has a slower Ω _b , model Y is intermediate and model Z has a faster Ω _b .

The formation of a bar is a secular process that can have several histories. Bars can form, change (increase, decrease, rotate etc.) and extinguish themselves with time, in processes that depend on the parameters of the galaxies that host such bars ^{Bournaud & Combes (2002)}. No analytical studies on gravitational bar potentials that evolve over time are known. In some recent works, such as like ^{Manos & Machado (2014)}, ^{Machado & Manos (2016)} and ^{Chaves-Velasquez et al. (2017)}, the authors wrote about the barred galaxy stability using time-dependent potentials. However, these studies were based on N-body simulations of barred galaxies by extracting parameters of the simulation for certain times in the system evolution, and treating each snapshot as a time-independent model.

Therefore, our intention is to carry out a study where the gravitational potential that represents the bar evolves over time in an analytical way. The idea is that our system starts as a barless galaxy and becomes a barred galaxy, whose bar grows over time.

We began the integrations with totally axisymmetric potentials, without bar, and over time we transformed these potentials into non-axisymmetric ones, with a bar.

Let us recall that the effective potential is given by Φeff(x)=Φ(x)-12|Ω×x|2 and the Lagrange points are five points where ▽Φeff=0. Writing the potential like this, we have a rotating system representation. The quantity EJ=12|v|2+Φeff is called the Jacobi energy and is conserved in the rotating system . Figure 1 shows the initial effective potentials and the final effective potentials for Models X, Y, and Z, where the emergence of the bar can be clearly seen.

Fig. 1 Effective potential contours of Models X, Y and Z. The top three images illustrate the initial effective potential contours, when the models did not yet have a bar, so the potential is axisymmetric. The three bottom images illustrate the final effective potential, when the bar is fully grown. Although there is no bar formed yet in the models of the first row, the effective potential, in a coordinate system that rotates with the pattern speed of the forthcoming bar, defines a radius of corrotation which is displayed in red. All images in the bottom row display the L1 − L5 Lagrange points. Twenty contours between energies −0.25 and −0.18 are displayed for each model. The color figure can be viewed online. 

In order to perform this bar evolution (shown in Figure 1), it was decided to implement a linear time dependent function of the mass of the bar in the Ferrers potential. In this process, two specific cases were created for each model: the evolution is completed in a time corresponding to 5 or 10 complete turns of the bar around itself. With this approach, we created the X ₅ ,X ₁₀ ,Y ₅ ,Y ₁₀ ,Z ₅ and Z ₁₀ notations. The bar evolution time for each model obviously depends on the pattern speed of the bar Ω _b and on the length of the bar. These times were calculated to serve as parameters of the growth function of the bar and are displayed in Table 2.

The Jacobi energy EJ=12|v|2+Φeff(x) is conserved in a rotating potential system representative of a fixed bar. However, this energy is not conserved during the evolution. While the mass of the bar is growing linearly, the E _J value of any particle also decreases linearly, and it is conserved again as soon as the bar growth finishes and the system bears a fixed bar (after evolution). To exemplify this statement, a random orbit was used, with an initially circular motion, integrated in Model X ₅. The integration begins without the bar, the bar structure starts to emerge and its mass increases linearly until the time 614.35 Myr (as shown in Table 2). After this evolution, the system has a fixed bar until the end of the integration at 10,000 Myr. This whole process is displayed in Figure 2.

Fig. 2 This image displays the E_J behavior for a random orbit as the bar grows in Model X5. It can be observed that E_J is not conserved during bar growth, but it is conserved after the bar has evolved. Notice that the evolution time is completely in agreement with that shown in Table 2 for Model X₅ (614.35 Myr). The integrations were made up to 10,000 Myr, and for this illustration we plot the time until 3,000 Myr. The E_J behavior for this orbit is not unique: for all integrated orbits this energy decrease occurs during the time of bar evolution. 

3.2. Initial Conditions

Galactic bars behave like rigid bodies, that is, Ω _b is always constant. However, galactic disks do not behave this way, their pattern speed Ω is a function of the radial coordinate. Thus, it is natural to imagine that resonances will appear between the bar and disk. An example is the corrotation resonance (CR), which occurs where Ω = Ω _b .

There will also be resonances when the following condition is satisfied:

where m is an integer related to the symmetry of the structure in which we are interested (m-armed spiral structures, bars etc.), and κ2=d2ΦdR2+3RdΦdR is the epicyclic frequency. In this case, there will be two resonances, the Lindblad resonances. In Equation (10), for the negative sign, there is the Inner Lindblad Resonance (ILR); for the positive sign, the Outer Lindblad Resonance (OLR). For a galactic bar potential we have m = 2 because of the bisymmetric structure. Figure 3 displays the curves Ω, Ω + 12κ and Ω - 12κ for Models X, Y and Z.

Fig. 3 Ω, Ω + 12κ and Ω - 12κ f curves for Models X, Y e Z and the corresponding CR, ILR and OLR resonances. The color figure can be viewed online. 

As the galactic bar is expected to always be contained in the CR radius, its influence does not exceed the OLR radius. Therefore, we only consider orbits with initial positions inside the OLR radius. For this study, we launched particles in initially circular orbits, distributed randomly from the center up to the OLR resonance, with 10,000 prograde orbits and 10,000 retrograde orbits for each model, starting from the positions shown in Figure 4. By prograde and retrograde orbits we mean orbits launched in the direct and opposite directions, respectively, in the bar corrotating non inertial reference frame. It is important to stress that we are dealing with motions of individual particles in a gravitational potential, i.e., this is not a self-consistent N-body simulation. The same number of prograde and retrograde orbits does not mean that we are weighting them equally, since it is known that prograde orbits play a much more important role in a barred galaxy potential. It just means we are exploring possible prograde and retrograde orbits with different initial conditions.

Fig. 4 Initial particle distributions for Models X, Y and Z. All orbits were launched with initial circular velocity and were distributed randomly inside the OLR resonance for each model. It is clear in the images that Model X has more scattered orbits than Model Y and, in turn, Model Y has more scattered orbits than Model Z; this happens because the greater the bar pattern speed Ω_b, the smaller becomes the OLR radius. This image shows dots representing 10,000 initial positions. Indeed, 10,000 prograde orbits and 10,000 retrograde orbits were computed for each model, starting from these same positions. 

Notice that the greater the bar pattern speed Ω _b , the smaller will be the OLR radius. According to our criterion for the choice of initial conditions, Model X has more scattered orbits than Model Y and, in turn, Model Y has more scattered orbits than Model Z. This phenomenon is clearly shown in Figure 4.

3.3. Results and Discussion

For efficiency, we inserted a condition in the SALI calculation on the LP- VIcode program to show us the moment when SALI < 10⁻⁸, which we consider close enough to zero to classify the orbit as chaotic. With this, we were able to create a classification for the chaos level of an orbit. The categories are as follows:

where t is measured in Myr ⁵.

Figures 5 to 7 show the amount of chaos that arose at each integration time for each model, for both prograde and retrograde orbits. The levels of chaoticity are represented by different colors. Apparently, orbits launched as retrogrades are more conducive to chaos. This observation is consistent with ^{Caritá et al. (2017)}; however it contradicts results of some classic investigations, (^{Athanassoula et al. 1983}; and ^{Pfenniger 1984}), where there is more order in the retrograde parts of the surfaces of section of these fixed potentials. Certainly, this will lead us to further investigations in a future article.

Fig. 5 Integration time × chaotic orbit number for Models X₅ and X₁₀. Model X₅, in which the bar grows faster, has slightly more chaos (mainly Level 1) than Model X₁₀. The colors of the line segments are defined in the legend inside the upper plot. The color figure can be viewed online. 

Fig. 6 Integration time × chaotic orbit number for Models Y₅ and Y₁₀. Model Y₅, in which the bar grows faster, has slightly more chaos (mainly Level 1) than Model Y₁₀. The colors of the line segments are defined as in Figure 5. The color figure can be viewed online. 

Fig. 7 Integration time × chaotic orbit number for Models Z₅ and Z₁₀. Model Z₅, in which the bar grows faster, has slightly more chaos (mainly Level 1) than Model Z₁₀. The colors of the line segments are defined as in Figure 5. The color figure can be viewed online. 

All models studied presented strong dominant Level 1 chaos for the retrograde orbits and an apparent domain of Level 2 chaos for the prograde orbits. Some Level 1 chaos is generated in the prograde orbits, and specifically in Models X ₁₀ and Z ₅ this type of chaos is null or practically negligible.

Figure 8 shows the cumulative number of chaotic orbits for prograde and retrograde orbits, and for the total number of orbits of the two types in all models. The distributions between order and chaos, considering the initial positions for each model, are shown in Figures 9 to 11.

Fig. 8 Time × cumulative number of chaotic orbits for Models X₅, X₁₀, Y₅, Y₁₀, Z₅ and Z₁₀. Looking at the sequence of images, one realizes that the number of orbits that become chaotic at some point is closely related to the bar pattern speed. Clearly in Model Z, in which the bar rotates faster, there is a greater amount of chaotic motions when compared to Model X, in which the bar is slower. While Figures 5 to 7 show a slight difference indicating that models in which the bar grows faster have slightly more Level 1 chaos, here by analyzing the general context, for the cumulative number of chaotic orbits nothing can be said. 

Fig. 9 Each plot shows the initial positions of all Model X particles, color-coded according to their chaos level, as in the legend above. This image again reinforces the result that orbits launched as retrogrades are more conductive to chaos. The outer limit for the inicial positions is the OLR, as stated in § 3.2. The initial positions of the orbits that have become chaotic are located in a well-defined ring-like region, confined between the ILR and CR resonances. ILR, CR and OLR circles are shown in black, and identified in the upper left plot. Almost no orbit presents chaos with initial conditions outside these regions. The color figure can be viewed online. 

Fig. 10 Each plot shows the initial positions of all Model Y particles. All plots are organized and color-coded as in Figure 9. The color figure can be viewed online. 

Fig. 11 Each plot shows the initial positions of all Model Z particles. All plots are organized and color-coded as in Figure 9. The color figure can be viewed online. 

In general, for the prograde orbits, the number of chaotic orbits generated was quite low, ranging from 5% to 10% of the total number of prograde orbits launched (depending on the model analyzed). On the other hand, for the retrograde orbits, this percentage increased considerably, to between 15% and 25% of the total number of retrograde orbits launched (depending on the model). In total numbers, considering the prograde and retrograde orbits, the percentage of chaos in the integrated orbits was always between 10% and 18% (depending on the model).

Figure 8 does not show appreciable changes with the rate of bar growth. However, in Figures 5 to 7 dissimilarities appear. The models where the bar grows over 5 turns seem to provide slightly more chaos than the models where the bar grows over 10 turns. This is an indication that an abrupt appearance of the bar causes more disturbance in the system.

Figures 9 through 11 show the stability of the orbits according to the initial positions. They all present a common feature: very well-defined ringlike regions of chaos. For the prograde orbits, there is only one ring of chaos for each model. For the retrograde orbits, there are two rings of chaos, a large and thick outer ring, surrounding a subtle inner ring, with one exception for Model Z ₅. The CR resonance limits these rings; in this context, interestingly, the most prominent rings are confined between the ILR and CR resonances, with a single exception for Model Z ₅. Very few orbits presented chaos with initial conditions outside these ring-like regions. In Figures 9 through 11 the difference in the amount of chaos for prograde and retrograde orbits is also clearly shown.

The greater the bar pattern speed Ω _b , the smaller the OLR radius. Figures 5 to 11 also show that the greater Ω _b , the more chaos the orbits that are within the OLR radius will present. In fact, Figures 5 through 7 make it clear that the main difference in this respect is especially in the retrograde orbits in the Level 1 of chaos. This indicates that the bar pattern speed also influences the system sensitivity to the bar appearance, since the orbits with Level 1 of chaos presented chaos in a time close to the appearance of the bar. On the other hand, for the prograde orbits, no significant differences are displayed at this point.

As expected, the E _J of all particles is not conserved during bar evolution. As already mentioned, while the mass of the bar is growing linearly, the value E _J also decreases linearly, and it is conserved again from the moment the bar growth finishes and the system becomes fixed (after evolution). To visualize this phenomenon, Appendix B presents some images where the number of orbits for certain times during and after the growth of the bar for each model is arranged by E _J . In these images, changes in the E _J values can be seen until the times listed in Table 2 are reached. Afterwards, E _J is conserved.