2.1. Galactic Models

The selected gravitational potentials are composed by a sum of the axisymmetric and bar components. The axisymmetric component itself is a superposition of several elements: a core and a bulge, both modeled by a Plummer potential (^{Plummer 1911}), a disk represented by a Miyamoto-Nagai model (^{Miyamoto & Nagai 1975}), and a halo described by a logarithmic potential (^{Richstone 1980}). The full axisymmetric potential is then

where r ² = x ² + y ² + z ², while R ² = x ² + y ². M _c1 and M _c2 are the masses of the bulge and the core, respectively, and r _c1 and r _c2 are their radial scale lengths. M _D is the disk mass and a and b its structural parameters. For the halo, v _H is the asymptotic velocity and r _H its radial scale length. By comparing the circular velocity profiles obtained from different parameter sets with the observed velocities in the Milky Way within the first 40 kpc, as reported by ^{Sofue (2020)}, and by employing the gradient descent method, we identified three parameter sets with nearly identical χ ² values. This indicates that the potentials corresponding to these parameters are degenerate. However, these potentials exhibit distinct characteristics from one another. Figure 1 displays the circular velocities derived from the three models, which we will henceforth refer to as Model 1 (left panel), Model 2 (middle panel), and Model 3 (right panel), based on the influence of their respective disks. The parameter values corresponding to each model are provided in Table 1.

Fig. 1 The circular velocity curves corresponding to the axisymmetric potential for three different sets of parameter values. CC means Central Component (bulge+core). These curves are compared with the circular velocity of the Milky Way (^{Sofue 2020}) within the first 40 kpc, along with their associated χ ². The color figure can be viewed online. 

2.2. Bar Models

For the bar component, we employed the ellipsoidal Ferrers potential. The mass density associated with this potential is defined as:

where ρBc=10532πMBaBbBcB is the central density and m=(x/aB)2+(y/bB)2+(z/cB)2, while a _B , b _B and c _B are the semi-axes of the ellipsoid. The index n is chosen to be n = 2. The forces generated by this potential are described in ^{Pfenniger (1984)}.

To simplify the models, we fixed the ratios between the semi-axes, setting b _B = a _B /3 and c _B = a _B /6. However, the major semi-axis a _B was varied, considering three distinct values: 6, 4.5, and 3 kpc. We introduced the bar component as a smooth time-dependent function by gradually transferring mass from the disk to the bar in the following way:

Here, MBf represents the final mass of the bar after its growth is completed at T _max which is set to 1 Gyr out of a total simulation time of 11.25 Gyr. Then, our models have a transient phase of 1 Gyr and are time independent after that. The evolution of the disk mass is then expressed as M _D (t) = M _Di M _B (t). Rather than directly referring to MBf, we will use μB=MBfMDi, which denotes the final fraction of mass transferred from the disk to the bar.

The parameter R=RCRaB (R _CR is the corotation radius, i.e., the radius at which the stars have the same angular speed as the pattern speed of the bar) is used to characterize the bar rotation rate. Bars are kinematically classified as ‘slow’ if ℛ > 1.4, ‘fast’ if 1.0 < ℛ < 1.4 or ‘ultrafast’ if ℛ < 1.0 (see ^{Debattista & Sellwood (2000)}; ^{Rautiainen et al. (2008)}; ^{Aguerri et al. (2015)}; ^{Lee et al. (2022)}). Concerning the ℛ found for our barred models, we must note that each of the models has different axisymmetric backgrounds (Model 1, 2 or 3), different bar semi major axis a _b , different μB=MBfMDi and different Ω _B . However, after the bar growth phase, with all parameters set, ℛ depends solely on R _CR (determined by Ω _B ) and a _B .

We are now able to analyze the bar rotation rate in a two-dimensional parameter space for ℛ corresponding to each axisymmetric model, along-side the three selected values for the bar length, as illustrated in Figure 2 for Model 1. This methodology allows for the identification of models featuring either fast or slow bars, with varying values of μ _B and Ω _B . Therefore, we focused on models with μ _B = 0.1, 0.2, 0.3, 0.4, 0.5 and Ω _B = 20, 25, 30, 35, 40, 50, 60 km/s: Furthermore, only models with 1 ≤ ℛ ≤ 3 were analyzed (indicated by red dots in Figure 2), corresponding to bars classified as fast and slow by ^{Aguerri et al. (2015)}.

Fig. 2 Parameter space for ℛ with respect to the parameters μ _B and Ω _B for model 1 and three values of bar length. The red dots are the models we studied in this article. Here we show the values only for Model 1, because there is no significant change of ℛ when considering the three axisymmetric models (Model 1, Model 2 and Model 3), once they are time independent (see text). From bottom to top, the straight lines in each plot denote equal ℛ: 3 and 2 (white solid lines) and 1.4 and 1 (black solid lines). The color figure can be viewed online. 

2.3. Initial Conditions and Orbit Integrations

Given our focus on stars within the disk, we selected the initial positions of the test particles following the axisymmetric distribution of the Kuzmin model whose surface density is (^{Binney & Tremaine 2008}):

We have chosen to use KT models to create the initial particle distribution since our Miyamoto-Nagai models have a scale length ratio (a/b) of 20, i.e. the disks are very thin. Furthermore, all the particles are initially in the galactic plane (z(t = 0) = 0). Considering the relationship between surface density and the number of stars, we observe that at a given radius, the number of stars can be expressed as N _K (R) = 2πΣ _K (R). By employing the Monte Carlo method, we can achieve a well-distributed arrangement of stars within the galactic plane. The maximum radius selected was given by R _max = 1.5R _CR . Although all particles are initially positioned within the galactic plane, the three-dimensional aspect of the simulations is established through the introduction of an initial velocity component along the z-axis.

For the velocity initial conditions, we decided to make v-θ=vc=RdΦaxdR, while v-R=v-z=0, since all the stars in an axisymmetric galaxy (as in our models for t = 0) have near-circular orbits. On the other hand, the velocity dispersion for each coordinate decreases exponentially with radius as noted by ^{Lewis & Freeman (1989)}:

where σ _R , σ _T and σ _z are the radial, tangential and vertical velocity dispersions, σ _R (0), σ _T (0) σ _z (0) are their central values, and h _R , h _T and h _z are the scale lengths. In their work, ^{Lewis & Freeman (1989)} computed the specific values of h _R and h _T for the Milky Way, which were found to be 4.37 kpc and 3.36 kpc, respectively. In the same work, Lewis and Freeman also assumed that the ratio of radial velocity dispersion σ _R to vertical velocity dispersion σ _z remains constant. Consequently, we have chosen h _z = h _R .

For our research, we have chosen to equate the three central velocity dispersions (σ _R (0) = σ _T (0) = σ _z (0) = σ _D ). Additionally, we explore three distinct values for this new parameter: 100, 80 and 50 km/s.

Hence, in total, we have 5 different galactic parameters: three different axisymmetric models, five distinct bar/disk mass ratios, seven angular velocities, three sets of bar lengths and three disk central velocity dispersion. This resulted in 945 different models. However, as mentioned before, only models with 1 ≤ ℛ ≤3 were studied, and then our - final set of simulations discussed here included 708 simulations. Each simulation comprised a total of 30,000 test particles, yielding a total of 21,240,000 calculated orbits. A concise summary of the model parameter space is presented in Table 2.

The integration was carried out using a fourth-order Runge-Kutta integrator, employing Fortran subroutines. To assess stability, we monitored the Jacobi energy of test particles following bar growth. Typically, ΔE _J is better than 10^-10 for t > T _max , since only then the models become time independent. The total simulation time was 11.25 Gyr, during which 1225 snapshots were captured at equidistant intervals. The bar mass growth ceases at T _max = 1 Gyr to ensure 1024 snapshots after the bar mass evolution. This number of points (2¹⁰) was chosen to simplify the Fast Fourier Transform analysis.

3.1. Detecting the Periodic Orbits x ₁ and x ₂

As highlighted by several authors (e.g., ^{Contopoulos (1970)}, ^{Athanassoula et al. (1983)}, ^{Sellwood & Wilkinson (1993)}, ^{Skokos et al. (2002a)} and ^{Patsis & Athanassoula (2019)}, periodic orbits play a crucial role in shaping the structure of barred galaxies. These periodic orbits are categorized into four primary families, namely the x ₁, x ₂, x ₃ and x ₄ family orbits, following the classification by ^{Contopoulos & Papayannopoulos (1980)}, with the most significant being the x ₁ and x ₂ families. Identifying these periodic orbits is essential for understanding the dynamics of barred galaxies.

To detect these orbits in our simulations, we applied a Fourier transform on the particle coordinates projected onto the equatorial plane: x(t), y(t), and R’(t) = R(t) -R-, where R- is the mean radius of the particle orbit. This allowed us to extract the dominant frequencies, i.e., those with the highest amplitudes, in each coordinate. These frequencies are denoted as ω _x , ω _y , and ω _R , respectively. Additionally, the corresponding amplitudes were determined and labeled as A _x , A _y , and A _R .

It is crucial to emphasize that this analysis was conducted after the bar has reached its final mass, i.e., for t > T _max , as in our approach the stellar orbits stabilize and maintain a consistent pattern after the bar growth. In addition, prior to performing the Fourier transform, we applied a Hanning window function to the orbit positions. The purpose of this window function was to mitigate signal ‘leakage’ in the Fourier spectra.

We can now identify sticky orbits around the x ₁ and x ₂ families of periodic orbits (^{Contopoulos & Harsoula 2008}; ^{Katsanikas et al. 2013}), but for simplicity, we will refer to them as members of the x ₁ or x ₂ families. An orbit is classified as part of the x ₁ family if it satisfies the condition 1.9 ≤ ω _R /ω _x ≤ 2.1 and A _x /A _y ≥ 2. On the other hand, orbits that meet the criteria 1.9 ≤ ω _R /ω _x ≤ 2.1 and A _x /A _y ≤ 0.5 are associated with either the x ₂ or x ₃ family. However, it is important to note that the x ₃ family is significantly less stable than the x ₂ family (^{Skokos et al. 2002b)}, and as such, the presence of sticky-chaotic orbits around x ₃ is expected to be minimal, making them insignificant for classification purposes.

Having identified the orbits belonging to the x ₁ and x ₂ families, we can now quantify the number of orbits in each family, denoted as Nx1 and Nx2, respectively. Alternatively, we can calculate the proportion of x ₁ and x ₂ orbits relative to the total number of elliptical orbits, i.e., orbits that satisfy 1.9 ≤ ω _R /ω _x ≤ 2.1. These proportions are denoted as Px1 and Px2, respectively.

3.2. Bar Strength

We also performed a Fourier analysis on the stars positions to calculate the bar strength. For this, we computed the m = 2 mode Fourier coefficients (a ₂ and b ₂) based on the particle positions located within an annulus of width ΔR at a radius R. Hence, as highlighted by ^{Chantavat et al. (2024)}, the amplitude of the bar is:

where the bar strength corresponds to the maximum value à ₂ within R _CR , expressed as:

We applied this method to each snapshot, making A ₂ a time-dependent parameter A ₂(t). This time evolution of A ₂(t) was then used to analyze the properties and dynamics of the orbits in our galaxy models. In Figure 3, we observe the temporal evolution of A ₂(t) for two given simulations. In the same figure, a smoothed version of this quantity is shown (the smoothing process employed a Savitzky-Golay filter with a sixth-order polynomial).

Fig. 3 Evolution of the bar strength A ₂(t) over time (black line) for two different models. The cyan line corresponds to the smoothed versions. Additionally, the white/black dot indicates the point of saturation in A₂(t), while the horizontal yellow line represents the secular value associated. Notice the different behavior of A ₂(t); in one simulation, A ₂(t) decreases after the time of saturation, while in the other model, A ₂(t) remains more or less constant. The color figure can be viewed online. 

Even with the smoothed curves of A ₂(t), analyzing each curve individually is impractical since we have a very large number of galaxy models; therefore, we seek for characteristic values that help us to evaluate the model.

One such value is the saturation point in A ₂(t). As observed in Figure 3, the A ₂(t) curves initially experience a rapid growth, followed by a decline. Notably, several curves exhibit this decline after a specific time. We designate the values at this inflection point as A _sat and t _sat .

Another characteristic value becomes evident toward the end of the simulations for A ₂(t). At this point, the curves exhibit minimal changes over time. This behavior is expected since the structure of rotating Ferrers bars is primarily supported by the stable portion of the x ₁ family. In a response model, any changes in A ₂ after a couple of bar revolutions beyond T _max can therefore be attributed to the influence of chaotic or escaping orbits. To quantify this stability in the A ₂ values, we calculated the average values over the last 1 Gyr. specifically, we denote these values as < A _sec >.

An additional observation is the difference between A _sat and < A _sec >. Consequently, we designate this difference as another characteristic value, denoted by ΔA = A _sat - < A _sec >.

5.1. Evolution of x ₁ and x ₂ Orbits and Double Barred Galaxies

It is essential to note that the methods introduced in § 4 are directly applied to the response model, allowing for the tracking of changes introduced whenever the parameter combinations depicted in Figure 4 occur. This approach provides significant practical value in the analysis of fully self-consistent N-body models, enabling a comprehensive assessment of the evolution of key orbital structures, particularly the x ₁ and x ₂ families.

In Figure 5, it is evident that the ℛ parameter is the second most influential feature for determining Nx1. Notably, the correlation coefficient between these two parameters in Figure 4 is strongly negative. Conversely, the most significant feature for Nx2 is Ω _B , which also exhibits a negative correlation.

Previous studies using N-body simulations (e.g., ^{Athanassoula 2003}; ^{Manos & Machado 2014}) have demonstrated that in barred galaxies, Ω _B tends to decrease over time, consequently ℛ increases. When combined with our findings, this suggests that barred galaxies may initially feature a fast bar, characterized by a high number of x ₁ orbits, and gradually lose some of these orbits as the system evolves. Simultaneously, as Ω _B decreases, the bar could be acquiring a larger number of x ₂ orbits.

Moreover, the presence of a significant concentration of x ₂ orbits could offer an explanation for the existence of a secondary bar in certain galaxies (^{Friedli et al. 1996}; ^{Maciejewski et al. 2002}; ^{Erwin 2004}; ^{Wozniak 2015}; ^{Erwin 2024}). This indicates that double-barred galaxies may be dynamically more evolved systems.

5.2. Characterizing Bar Strength Evolution Through ΔA

Figure 3 presents two distinct bar strength curves over time, A ₂(t). The first curve exhibits rapid growth until it reaches a saturation point. After saturation, the bar rapidly weakens until it reaches an equilibrium (as expected for response models) with < A _sec >. In contrast, in the second curve the bar also weakens beyond the saturation point, although at a much slower rate. Looking at the other A(t) curves we can note that there are intermediate cases.

The parameter ΔA plays a crucial role in understanding this behavior. A large positive value of ΔA indicates that A _sat is larger than < A _sec >, aligning with the behavior observed in the first case in Figure 3. Conversely, a small ΔA corresponds to a behavior more related to the second case shown in the same figure.

As previously discussed, thex ₁​ family of orbits forms the backbone of rotating Ferrers bars, remaining stable in the region that supports the bar. In a response model, the decreasing of ΔAcould be attributed to the presence of chaotic or escape orbits within the system. To verify this last statement, a study using an index to quantify the chaotic behavior of the orbits as GALI2 (^{Skokos et al. 2007}; ^{Chaves-Velasquez et al. 2017}; ^{Caritá et al. 2019}) could be performed. However, this is beyond the scope of the present paper.

The permutation feature importance analysis for ΔA(Figure 6) shows that the primary feature parameter, with a significantly larger importance compared to other feature parameters, is ℛ. Furthermore, the Spearman correlation between ΔAand ℛ (as shown in Figure 4) is negative. From this, we can infer the following:

5.3. Impact of Disk-to-Halo Ratio on Bar Formation

In § 2.1, we constructed three distinct axisymmetric models. These models exhibit the circular velocity closest to that of the Milky Way within the first 40 kpc, as shown in Figure 1. Despite exhibiting degeneracy in terms of circular velocity, discernible differences exist among them, primarily in the disk-to-halo ratio. Model 1 shows a strong disk influence compared to the halo for R < 8 kpc, whereas Model 3 shows a stronger halo influence compared to the disk from R > 3 kpc. Model 2 represents an intermediate case.

Since earlier works by ^{Athanassoula & Sellwood (1986)}, it is known that galaxies with a stronger disk influence, as seen in Model 1, tend to form bars more rapidly than those in which the halo is predominant, like in Model 3 (see also ^{Valencia-Enríquez et al. 2023}). However, the permutation importance for all eight result values shown in Figures 5 and 6 indicate that the variations among the axisymmetric models have a negligible impact on the result values. This finding is corroborated by the Spearman correlation coefficients presented in Figure 4, where it is evident that most of the coefficients related to the axisymmetric model are nearly zero. The most significant correlation is for ΔA with a coefficient of -0.20, which is insufficient to establish a strong anticorrelation.

The apparent contradiction presented here arises from the context of this work, which belongs to rigid potential models. It is noteworthy that when one performs a N-body fully self-consistent models, the structural parameters of the galactic components exhibit temporal evolution, and there is a transfer of angular momentum among the components. In our research, we are imposing the same bar model characterized by parameters μ _B , a _B , and Ω _B across all three axisymmetric models.

Consequently, we can conclude that for degenerate models, using rigid potentials, the variation in disk-to-halo does not significantly affect the the formation of x ₁ and x ₂ orbits.

Finally, it is crucial to highlight that combining the Spearman correlation coefficients with the feature importance derived from a Random Forest Regressor significantly enhances the analysis of the effects of different input parameters on output results. The Spearman correlation provides insight into the monotonic relationships between values, while the feature importance in a Random Forest Regressor evaluates the overall significance of each parameter. Using both methods allows for a comprehensive understanding of the importance and behavior of various parameters.

In conclusion, our work provides a detailed investigation into the dynamics of barred galaxies, offering insights into the interplay between various galactic parameters and the formation and evolution of galactic structures. The methodologies employed and the findings derived from this study contribute to the broader understanding of galactic dynamics and serve as a foundation for future research in this field.