Orbital Structure Evolution in Self-Consistent N-body Simulations

We use the N-body simulations described in Valencia-Enríquez et al. (2019) (hereafter VE19) to analyze the orbital frequencies. VE19 presented two sets of simulations: four isolated models and thirteen perturbed ones. In this study, we focus on conducting a frequency analysis for three isolated models. The models begin as an axisymmetric disk, which, due to internal dynamical instabilities (Athanassoula, 2002, 2003, and references therein), develop strong bars that persist until the end of the simulations. The isolated models employed an exponential stellar disk $\Sigma(r)=\Sigma_{0}e^{r/r_{d}}$ where $\Sigma_{0}=M_{d}/(2\pi r_{d})$, and an NFW dark matter halo (Navarro et al., 1996, 1997). The vertical mass disk distribution is composed of an isothermal sheet with a radially constant vertical scale length $z_{0}$. Therefore, the three-dimensional stellar density in the disk is $\rho_{d}(r,z)=\Sigma(r)/(2z_{0})sech^{2}(z/(2z_{0}))$.

The models were built as an equilibrium N-body realization (Springel & White, 1999) of $7\times 10^{6}$ particles ($2\times 10^{6}$ for the disk and $5\times 10^{6}$ for the halo). The halo has a mass of $5.11\times 10^{11}$ $M_{\odot}$ and a concentration of 8.0. The disk has a mass of $2.55\times 10^{10}$ $M_{\odot}$ and radial scale-length changes for the different models. For the $A\lambda 03$, $A\lambda 04$, and $A\lambda 05$ models discussed here, the radial scales are $r_{d}=1.99$, $r_{d}=2.80$, and $r_{d}=3.58$ kpc, respectively, and the models have $z_{0}=0.2r_{d}$. We performed collisionless N-body simulations with the Gadget-2 code (Springel et al., 2001; Springel, 2005).

Models described in VE19 were chosen to represent galaxies where the disk dominates (the Critical Spin parameter $\lambda_{c}$ is larger than the Spin parameter of the disk $\lambda_{d}$) and models where the halo dominates ($\lambda_{c}<\lambda_{d}$). VE19 showed that the isolated model keeps that configuration during the whole evolution, and the disturbed models change such parameters at the interaction. Therefore, models form a bar if $\lambda_{c}>\lambda_{d}$ for both isolated and interaction models. Besides, disk-dominate models develop a bar relatively quickly, and halo-dominate ones do not form that structure. The disk-dominate models evolve strong bars due to the exchange of angular momentum between the components of the galaxy (e.i. halo and disk) and between the inner and outer region of the disk (Athanassoula et al., 2013). This process is more efficient if the Critical Spin parameter is relatively much larger than the Spin parameter of the disk.

In this article, we recalculate the simulations of the isolated models $A\lambda 03$, $A\lambda 04$, and $A\lambda 05$ from VE19, now with a higher time resolution. This is very important to have a better determination of the orbital frequencies (see Section 3).

The orbits of stars in a galaxy describe oscillations in three dimensions. In Cartesian coordinates ($x$, $y$ and $z$), the $\omega_{x}$, $\omega_{y}$, and $\omega_{z}$ frequencies represent those oscillations. Likewise, in cylindrical coordinates, the $\Omega$ and $\omega_{R}$ frequencies represent the azimuthal and radial oscillations; $\omega_{z}$ represents the vertical oscillation for both coordinate systems.

In order to analyze these frequencies in our models, we take the bar frame of reference where the bar is in a horizontal position for all snapshots (see Figure 1). We calculate the frequencies of the orbits from each disk particle using the phase space coordinates as a function of discrete time ($x(t)$, $y(t)$, $z(t)$, $v_{x}(t)$, $v_{y}(t)$ and $v_{z}(t)$ for Cartesian coordinates, and the azimuthal position $\theta(t)$ and the radial velocity $v_{r}(t)$, for cylindrical coordinates).

We use an algorithm similar to the numerical analysis of fundamental frequencies, developed by Laskar (1990). To apply the frequency analysis, we find the peaks of amplitude from the Fourier Transform of a time series of a given coordinate, e.g., $x(t)$. The continuous Fourier Transform coefficients, $H_{x}(w)$, have the following form

Because $x(t)$ results from a simulation, we have a finite number of sampled points where we have $N$ consecutive sampled values $x_{k}=x(t_{k})$, equally spaced in the interval time $\Delta t$, so $t_{k}=k\Delta t$, and $k=0,1,2,...,N-1$. The integral approximates as the sum:

To perform the analysis, we first subtract the average value of $x(t)$ from the signal to eliminate the spurious amplitudes that can appear at the zero frequency. Next, we filter the signal by multiplying it with the Hanning window to reduce side lobes. Then, we apply the Fast Fourier Transform (FFT) for a real-valued sample using the realft subroutine of Press et al. (1992). We search for a maximal value of the $H_{x}(w_{k})$ then we save the points $(w_{k},H_{x}(w_{k}))$, the previous value $(w_{k-1},H_{x}(w_{k-1}))$, and the next value $(w_{k+1},H_{x}(w_{k+1}))$. With these three points, we fit a parabolic curve to get the vertex position. That value corresponds to the frequency $\omega_{x}$. We check our code reproducing the orbits showed in Carpintero & Aguilar (1998), and we obtain the same frequencies; likewise, we test our code with a periodic function $x^{\prime}(t_{k})=H_{x}\sin(w_{k}t_{k}+\phi)$, and the calculated $w_{k}$ agrees very well with the input $w_{k}$.

We plot in Figure 3, the frequency ratio maps $\Omega/\omega_{z}$ versus $\omega_{R}/\omega_{z}$ for all disk particles of our models, and for all intervals of time. The frequency $\Omega$ here is just $\omega_{\theta}$, the frequency related to the cylindrical coordinate $\theta$. As discussed in Athanassoula (2002),Ceverino & Klypin (2007), and Gajda et al. (2016), the frequency $\omega_{\theta}$ is much more complex to calculate and in some cases, the direct Fourier Transform may result in erroneous values. In our case, we only plot this frequency in Figures 3 and 4, and we continue our main analysis using the Cartesian coordinates. We observe that the bar structure initiates at different times for the various models. The 2:1, 5:2, and 3:1 resonances are distinguished in their population as the models evolve. In the case of the disk-dominated model $A\lambda 03$, these resonances are swiftly populated. Conversely, for the other models, the population of these resonances occurs at later times.

For $A\lambda 03$ model (panels in the left column in Figure 3), in the first Gyr of the evolution, the resonance 2:1 begins to be populated with particles with low $\Omega/\omega_{z}$ and $\omega_{R}/\omega_{z}$: these particles have high $\omega_{z}$. Out of this resonance, the particles are distributed main at low $\Omega/\omega_{z}$ and at a large range of $\omega_{R}/\omega_{z}$. As the model evolves, particles begin to populate other resonances. For the time interval T=[1-2], both 2:1 and 3:1 are populated. From T=2 to the end of the evolution, a large number of particles are trapped in the 2:1 resonance and around the 5:2 one. Furthermore, note how these particles move on these resonances from low to high $\Omega/\omega_{z}$ and $\omega_{R}/\omega_{z}$. This movement is due to the decreasing of the vertical frequency $\omega_{z}$ of these particles as the buckling of the bar structure appears. Orbits not belonging to the main resonances do not suffer this effect.

The panels in the middle column in Figure 3 show the evolution for model $A\lambda 04$. Here, at the early stages of the model evolution, particles are caught at some resonances (note the radial structures in T=[0-1] and T=[1-2]). At the time interval T=[2-3], the resonances 2:1 and 3:1 begin to be populated. As in the model $A\lambda 03$, particles here also move from the 3:1 to the 5:2 resonances and suffers the same effect we have noticed: the particles move on these resonances from low to high $\Omega/\omega_{z}$ and $\omega_{R}/\omega_{z}$ as the buckling appears.

Finally, panels in the right column of Figure 3 represent the evolution of the halo-dominated $A\lambda 05$ model. In this case, numerous particles stay around the 5:2 and the 3:1 resonances for the first 3 Gyrs of the simulation. Only then the main 2:1 resonance begins to be populated. Moreover, the particles on these resonances also undergo the effect of moving to large values of $\Omega/\omega_{z}$ and $\omega_{R}/\omega_{z}$, but the effect is much less important. This is attributed to the very late appearance of the buckling phase in this model.

To understand the number of azimuthal periods used in the calculation of particles frequencies, we constructed an accumulative figure representing the number of orbits normalized to one as a function of the azimuthal rotations of the particle orbits (see Figure 4). It is important to note that, in order to obtain an accurate frequency calculation using the Fourier Transform, the signal must contain at least one complete cycle within the analyzed time interval. This implies that the minimum number of periods required in the signal depends on the highest frequency present in the signal and the duration of the analyzed time interval. The highest frequency in the signal is related to the sampling frequency used in the data acquisition process. According to the Nyquist-Shannon sampling theorem, to avoid the phenomenon known as "aliasing" and accurately recover the frequencies present in the signal, the sampling frequency must be at least twice the highest frequency present in the signal (Oppenheim et al., 1997).

In the case of models $A\lambda 03$, $A\lambda 04$, and $A\lambda 05$, the frequency calculations are accurate because over 97% of the particles show two or more rotations. Specifically, for model $A\lambda 03$, more than 50% of the particles exhibit over twenty cycles; for model $A\lambda 04$, more than 50% of the particles have approximately eight cycles or more; and for model $A\lambda 05$, 50% of the particles have five rotations (see Figure 4).

To better quantify the number of particles that are grabbed by the bar potential as the models evolve, we describe the evolution of the frequencies ratio $\omega_{R}/\omega_{x}$ as the bar growth. Figure 5 shows the histograms for the $\omega_{R}/\omega_{x}$ quotient at different time intervals for our models and the Frozen potential experiments. The left column of this figure shows the results of the N-body models, while the center and right columns show the results for the frozen potential. The frozen potentials were calculated using a snapshot from the middle of each time interval. The orbits are then integrated during 1 Gyr (center column) and 6 Gyr (right column). The details of the Frozen potential calculations are described in Appendix A.

For the disk-dominated model $A\lambda 03$, during the T=[1-2] interval (marked by the red line), the model exhibits a significant peak at $\omega_{R}/\omega_{x}\approx 2.0$ and a secondary peak at $\omega_{R}/\omega_{x}\approx 2.8$. The latter peak corresponds to particles trapped around the 3:1 resonance. Notably, particles trapped around the 2:1 and 3:1 resonances already display a prominent bar structure in the model (refer to Figure 1 and Figure 10). As the model evolves, the peak around $\omega_{R}/\omega_{x}\approx 2.0$ grows, and the 3:1 resonance particles transition into a 5:2 configuration (indicated by the peaks at $\omega_{R}/\omega_{x}\approx 2.4$). Additionally, we observed a decrease in the significance of the peak around the 5:2 resonance, while the 2:1 resonance increasingly traps more particles.

The evolution of the $\omega_{R}/\omega_{x}$ ratio for the $A\lambda 04$ model is shown in the middle row of Figure 5. It is only at the interval T=[2-3] when the resonance 2:1 begins to be populated (green line). For latter time intervals, the evolution is similar to the $A\lambda 03$ model the peak at $\omega_{R}/\omega_{x}\approx 2.0$ increases, while particles from the 3:1 resonance move to the 5:2 one.

The evolution for the $\omega_{R}/\omega_{x}$ ratio for the halo-dominated model $A\lambda 05$ (bottom row of Figure 5) is quite different. For instance, the histogram shows a considerable amount of particles, which are around $\omega_{R}/\omega_{x}\approx 3.0$, but these particles do not exhibit any barred structure (see Figure 1). Then, the amount of these orbits decrease (see the cyan line) as the peak moves for lower $\omega_{R}/\omega_{x}$, up to the appearance of the oval distortion about the fourth Gyr (cyan line). In the next time interval (T=[5-6]), the bar seems to lost particles in the $\omega_{R}/\omega_{x}\approx 2.0$ quotient; however, at the next time interval, this peak begins to be more significant again. On the other hand, the behavior depicted in the results of the frozen potentials (see second and third columns of Figure 5) is similar to the N-body models. The importance of the peaks appears in the same resonances $\omega_{R}/\omega_{x}$ that the N-body models for both time integrations (1 and 6 Gyr), showing that 1 Gyga-year is sufficient to get the frequencies in an N-body model. The slight change in the pattern speed rotation for an interval of 1 Gyr does not considerably affect the results when we measure the frequencies from an N-body model.

Portail et al. (2015) showed that orbits with an elongated shape along the bar major axis are identified by the frequency ratio $\omega_{R}/\omega_{x}\approx 2$, which are related to orbits on the 2:1 resonance. Therefore, we choose the orbits in the interval $1.9<\omega_{R}/\omega_{x}<2.1$ and study their evolution. Figure 6 shows the disk mass fraction of orbits between $1.9<\omega_{R}/\omega_{x}<2.1$ ratio, hereafter B-type orbits (top panel). Additionally, we show the disk mass fraction to orbtis between $2.4<\omega_{R}/\omega_{x}<2.6$ (middle panel), and $2.9<\omega_{R}/\omega_{x}<3.1$ (button panel). The filled triangle, star, and circle symbols represent the $A\lambda 03$, $A\lambda 04$, and $A\lambda 05$ models, respectively. The color of the symbols depicts the time interval as in Figure 5.

The growth of B-type orbits in our models, depicted in the top panel of Figure 6, is different; it depends on the initial spin parameter of the disk $\lambda_{d}$. For the $A\lambda 03$ model, these orbits reach a value around 30% of the total disk mass fraction at the bar saturation phase¹¹1Bar saturation phase is the time when the amplitude of 1-D Fourier Transform coefficients for mode $m=2$ for disk particles is the maximum (see Figure 8 from VE19). (T=[1-2]), decreasing to around 24% (T=[2-3]) and increasing again to achieve 40% at the end of the simulation (T=[5-6]). The $A\lambda 04$ model has a larger number of B-type orbits at the bar saturation phase showing a maximum of 32% of the disk mass fraction. Similar to the $A\lambda 03$ model, the number of B-type orbits decreases after the bar saturation phase and later increases approximately up to 35% of the disk mass fraction at the end of the simulation. The $A\lambda 05$ model gets less than 30% of B-type orbits at the bar saturation phase (T=[4-5]) and later decreases drastically to around 12% of the disk mass fraction (T=[5-6]). To gain further insights, we continue monitoring the simulation of this model for an additional Gyr, aiming to determine whether the bar structure is ultimately destroyed or if it regrows (Bournaud & Combes, 2002). In Figure 6, we observe that the fill circles, which represent model $A\lambda 05$, extends up to 7 Gyr. We note that the disk mass fraction of B-type orbits increase in this last Gyr up to around 16%.

Furthermore, the middle panel of Figure 6 illustrates the evolution of the ratio $2.4<\omega_{R}/\omega_{x}<2.6$, revealing a significant increase as the saturation phase of the bar concludes. We observe an approximately 10% increase for the $A\lambda 03$ model, a 17% increase for the $A\lambda 04$ model, and a similar growth pattern in this type of orbit for the $A\lambda 05$ model. On the other hand, in the lower panel, corresponding to orbits with a frequency ratio of $2.9<\omega_{R}/\omega_{x}<3.1$, we observe a slight increase during the saturation phase of the bar for the $A\lambda 03$ and $A\lambda 04$ models. However, this increase is more pronounced in the case of the $A\lambda 05$ model, indicating a slower formation of the bar in this configuration.

Now, we focus on the frequency ratio $\omega_{z}/\omega_{x}$ to identify the orbital evolution which leads to the formation of the peanut structure. Figures 7, 8, and 9 show the number of B-type disk particles ($1.9<\omega_{R}/\omega_{x}<2.1$) as a function of $\omega_{z}/\omega_{x}$ ratio measured in different intervals of time. The left column shows the frequencies measured from the N-body model. The middle and right columns show the frequency ratio $\omega_{z}/\omega_{x}$ for the frozen potentials over 1 and 6 Gyr of integration time. In general, we notice that the shape of the curves of the frozen potentials results are very similar to the results of the N-body models. For instance, the highest peaks contribution appears in the same ratio frequency as the measures calculated in the N-body model, which means that the ratio frequency obtained from the N-body simulation for an integration time of 1 Gyr is enough to get a good approximation of the frequencies of the orbits.

The $A\lambda 03$ model, depicted in Figure 7, shows that the bar structure is formed by particles spread in the frequency ratio $2.0<\omega_{z}/\omega_{x}<3.5$ during the earliest times (top left panel). Thus, it corresponds to the initial stages of bar formation before the buckling phase occurs (see Figure 1). Similarly, orbit scattering occurs for calculations of frozen potential by 1 Gyr (top middle panel); however, orbits with higher frequency ratios are lost when calculating the frozen potential by 6 Gyr. The aforementioned effect could happen because it is challenging to maintain the frequency of orbits for more than 1 Gyr in the potential where the bar is still weak. Later, when the bar is fully established, and the orbits start to swing from the disk plane (the buckling phase begins), all the stellar orbits fall within the frequency ratio interval $1.4<\omega_{z}/\omega_{x}<2.2$ from the second to the sixth Gyr.

We can relate Figure 7 with the classification developed by Portail et al. (2015). They classified the orbits according to the frequencies ratios and the shapes using the letters A, B, C, D, E, and F for the intervals $1.45<\omega_{z}/\omega_{x}<1.55$, $1.55<\omega_{z}/\omega_{x}<1.65$, $1.65<\omega_{z}/\omega_{x}<1.75$, $1.75<\omega_{z}/\omega_{x}<1.85$, $1.85<\omega_{z}/\omega_{x}<1.95$, and $1.95<\omega_{z}/\omega_{x}<2.05$, respectively. Therefore, it is possible to identify these frequency intervals in the histograms of Figure 7 for the N-body model (left column) and its frozen potential (middle and right columns). In this manner, the first peak coincides with a frequency ratio from 1.4 to 1.5, which aligns with the A classification of Portail et al. (2015). This growth is throughout the evolution of the bar structure. A second peak, which is close the frequency ratio $\omega_{z}/\omega_{x}=1.8$ from the fourth to sixth Gyr, corresponds to the D classification. And the last peak corresponds to the F classification ($\omega_{z}/\omega_{x}=2.0$); this peak keeps almost equal during the evolution of the bar structure. We can observe that during the buckling phase, around the third Gyr, a considerable number of F and A orbits appears. Then, the peak for the F orbits is more or less constant, and the first and second peaks increase (A and D orbits) as the bar evolves. The $1.4<\omega_{z}/\omega_{x}<1.5$ interval (D orbits) has the highest peak. Furthermore, all other types of orbits increase the contribution to the bar structure in some extent. It can be noted unlike the N-body model, in the calculations made for the frozen potential integrated by 6 Gyr, the peak located in the ratio 1.4 has a greater contribution than the other frequency intervals, which corresponds to orbits populating in the center of the galaxy, (see Figure 12).

Figure 8 shows the frequency ratio $\omega_{z}/\omega_{x}$ evolution for the $A\lambda 04$ model. As we can observe, in this case, the outcomes of frequency measures from the frozen potential are very similar to the N-body model, where the form and the contribution of the ratio of frequencies appear in the same intervals (see middle and right columns of this figure). At earlier stages of the bar formation, we observe that our B-type orbits have a histogram similar to the one from model $A\lambda 03$, but delayed in time and have larger $\omega_{z}/\omega_{x}$ (compare T=[1-2] for model $A\lambda 03$ and T=[2-3] for model $A\lambda 04$). In the next Gyr interval, the histogram for model $A\lambda 04$ moves to lower $\omega_{z}/\omega_{x}$. Then, when the bar begins to develop the buckling phase, the orbits are established in the frequency ratio interval $1.4<\omega_{z}/\omega_{x}<2.2$, similar to the $A\lambda 03$ model. At the end of the simulation, a high peak in the $1.7<\omega_{z}/\omega_{x}<1.9$ ratio appears, which aligns with the D classification, according to Portail et al. (2015).

The halo-dominated model ($A\lambda 05$) shows quite different behavior. In the first instance, a much lower number of disk particles stand in the $1.9<\omega_{R}/\omega_{x}<2.1$ range (B-type orbits). That is why the maximum scale in the left column of Figure 9 reaches only half of the value plotted for the other two models (see Figures 7 and 8). At the early stages of the bar formation, the orbits locate around the frequency ratio $\omega_{z}/\omega_{x}\approx 3$. These high $\omega_{z}/\omega_{x}$ values represent a morphology of a boxy structure. Similar to the $A\lambda 03$ and $A\lambda 04$ model, the histogram moves toward the left as the bar evolves, but we notice that the $A\lambda 05$ model does not enter the buckling phase. For this model, it is only at the end of the simulation when particles participating in the bar structure present a distribution with lower $\omega_{z}/\omega_{x}$; the orbits are beginning to oscillate outside the disk plane with lower $\omega_{z}/\omega_{x}$, and only then they will be able to evolve from a boxy to a peanut shape. We should note that the results from the frozen potential of this model are also quite similar to the N-body model: the histograms appear in the same quotient of frequencies, and they have very similar behavior.

We have performed a similar analysis in the GravPot16 potential that resembles the main complexities of the Milky-Way galaxy (Chaves-Velasquez et al., in preparation). This potential has been widely used in works related to the dynamics of the Galaxy (see for instance, Fernández-Trincado et al., 2020, 021a, 021b). Our histograms show distributions of frequencies at high values for the ratio $\omega_{z}/\omega_{x}$ coefficient, similar to models $A\lambda 03$ and $A\lambda 04$ at the early stages of the simulation before the galactic buckling. We have studied the orbital families that shape the B/P morphology, and we found that orbits at higher vertical bifurcations of the planar x1 family (x1v3, x1v5, x1v7 according to the notation of Skokos et al., 2002b) contribute to shaping the edge-on view of the bar.

We now discuss the morphologies that arise from assembling 2:1 resonance orbits considering their face-on, edge-on, and end-on views and their characteristic frequency ratios. To distinguish the contribution to the shape of the bar structure, Figure 10 shows orbits assembly density profile (OADP) in their different views for B-type orbits of different $\omega_{z}/\omega_{x}$; the range of $\omega_{z}/\omega_{x}$, which is used to construct these maps, are displayed in the histograms of Figures 7, 8, and 9 (thick vertical lines). In Figure 10, the left, middle, and right columns represent the $A\lambda 03$, $A\lambda 04$, and $A\lambda 05$ models, respectively. The range of $\omega_{z}/\omega_{x}$ used for each model at each time interval is written in the bottom left of the XY views.

We observe that at the early stages of the bar formation, the orbits that participate in the oval/boxy distortion have high values of $\omega_{z}/\omega_{x}$. The disk-dominated model develops mainly an oval shape, and the other models display a more boxy morphology. The models arrive at this stage at different times; T=[1-2] Gyr for the disk-dominated model, T=[2-3] Gyr for the intermediate one, and T=[3-4] Gyr for the halo-dominated model.

The next stage in the evolution of the bar structure is the buckling phase. This phase can also be recognized in Figures 7, 8, and 9. The histograms for the frequency ratio $\omega_{z}/\omega_{x}$ move towards lower values of this ratio. That means the particles will oscillate less in the vertical direction while undergoing the oscillation following the bar potential. The $A\lambda 03$ and $A\lambda 04$ models evolve to that phase at T=[2-3] Gyr and T=[4-5] Gyr, respectively. For the halo-dominated $A\lambda 05$ model, this phase is not yet fully achieved at the end of the time simulation (T=7 Gyr).

In Appendix B, we present the OADP maps for our models in bins of 0.1*$\omega_{z}/\omega_{x}$ for the reader better understand the contribution of different kinds of orbits to the morphology of the bar structure.

We presented in Figures 7, 8, and 9 the evolution of the $\omega_{z}/\omega_{x}$ ratio distribution, which means that the orbits populated different values of this ratio during the growth of the bar structure. To illustrate this effect, in Figure 11, we show the evolution of five orbits for the N-body model $A\lambda 03$ which are B-type ($1.9<\omega_{R}/\omega_{x}<2.1$) at least in the last time interval of analysis (T=[5-6] Gyr). In that figure, we show xy, xz, and yz views of each orbit for five different time intervals. We can recognize several orbits morphology, which appears in models using a Ferrers analytical bar (e.g., Skokos et al., 2002a, b; Patsis et al., 2002, among others).

The ‘a’ orbit of Figure 11 is swinging on the bar region between $2.3<\omega_{z}/\omega_{x}<4.0$ to finally be part of the bar 2:1 resonance at the last time interval. We can perceive that this type of orbit contributes to the bar structure, but not to the peanut shape. It arrives later to the 2:1 resonance and stays on the disk plane. Likewise, we notice that these types of orbits do not add so much to the bar potential because the contribution of these orbits is low (see Figure 7).

The ‘b’ orbit is part of the B-type in the whole of its evolution. This orbit first appears on $\omega_{z}/\omega_{x}=2.66$ frequency ratio. Next, it moves to $\omega_{z}/\omega_{x}=1.93$ at T=[2-3]. Finally, for the following time intervals, this orbit contributes to the peanut shape of the bar in the edge-on view. Moreover, at the last time interval, this orbit shows a pretzel shape in the face-on view, a peanut shape in the edge-on view, and an X shape in the end-on view. Besides, this orbit contributes to the bump of the distribution located in $\omega_{z}/\omega_{x}\approx 2.0$ in Figure 7. This orbit contributes to the bin located at about $\omega_{z}/\omega_{x}\approx 1.8$ in the rest of the simulation.

The $\omega_{z}/\omega_{x}$ ratio in the ‘c’ orbit decreases from 3.32 to 2.00, contributing to the buckling of the bar (first two panels, T=[1-2] and T=[2-3]). This orbit shows a banana shape in the edge-on view. Then, the $\omega_{z}/\omega_{x}$ ratio decreases to 1.85 and remains almost constant during the evolution. At the last time interval, this orbit shows a peanut shape in the edge-on view and an elongated rectangular shape in the face-on and end-on views.

For the ‘d’ trajectory, the $\omega_{z}/\omega_{x}$ frequency ratio experiences a slight reduction from 1.93 to 1.83. Initially, it closely aligns with the bin distribution around $\omega_{z}/\omega_{x}\approx 2.0$. This orbit exhibits a clear banana shape in the xz view and a pretzel shape in the xy view (T=[2-3]). Over time, the pretzel shape evolves into a fish shape in the face-on view. Eventually, the orbit displays a boxy contour in the face-on view, a peanut shape in the edge-on view, and an oval shape in the end-on view.

The fifth orbit contributes to the central part of the bar; the $\omega_{z}/\omega_{x}$ ratio moves from 2.17 to 1.51 during the evolution. This orbit has an elongated shape in the face-on and edge-on view. At the last time-interval, the orbit depicts a boxy/peanut shape in the end-on view. These kinds of orbits contribute to the central boxy structure of the model (see Figure 10).

Considering the previous description of the orbits in Figure 11, it can be said that the particles change their frequencies during the evolution of the model. In general, the $\omega_{z}/\omega_{x}$ ratio decreases from high values before the buckling phase to lower values after the buckling (see Figure 7). This effect happens earlier for disk-dominated models. The orbits change their morphologies as they evolve with the bar potential; therefore, different orbital families in the models, which form the backbone of the bar structure in a disk galaxy, are assembled by distinct particles during the evolution. We must note that this effect is not exhibited when we study orbits using Ferrers analytical potentials or when we ‘freeze’ a simulation, calculate the potential produced by the particle distribution, and calculate orbits in this numerical potential. These techniques do not permit the self-consistent evolution of the bar potential.