# Orbital Structure Evolution in Self-Consistent N-body Simulations

Diego Valencia-Enríquez1, Ivânio Puerari2, and Leonardo Chaves-Velasquez3,4,5
1Universidad Mariana, Calle 18 No. 34 - 104, 52001 Pasto, Colombia
2Instituto Nacional de Astrofísica, Óptica y Electrónica, Calle Luis Enrique
Erro 1, Santa María Tonantzintla, 72840 Puebla, Mexico
3Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de
México, Apdo. Postal 3-72, Morelia, 58089 Michoacán, Mexico
4Astronomical Observatory, Universidad de Nariño, Sede VIIS, Avenida
Panamericana, Pasto, Nariño, Colombia
5Departamento de Física de la Universidad de Nariño, Torobajo Calle 18 Carrera
50, Pasto, Nariño, Colombia E-mail<EMAIL_ADDRESS>(COL)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The bar structure in disk galaxies models is formed by different families of
orbits; however, it is not clear how these families of orbits support the bar
throughout its secular evolution. Here, we analyze the orbital structure on
three stellar disk N-body models embedded in a live dark matter halo. During
the evolution of the models, disks naturally form a bar that buckles out of
the galactic plane at different ages of the galaxy evolution generating boxy,
X, peanut, and/or elongated shapes. To understand how the orbit families hold
the bar structure, we evaluate the orbital evolution using the frequency
analysis on phase space coordinates for all disk particles at different time
intervals. We analyze the density maps morphology of the 2:1 family as the bar
potential evolves. We showed that the families of orbits providing bar support
exhibit variations during different stages of its evolutionary process,
specifically prior to and subsequent to the buckling phase, likewise in the
secular evolution of the bar. The disk-dominated model develops an internal
boxy structure after the first Gyr. Afterwards, the outer part of the disk
evolves into a peanut-shape, which lasts till the end of the simulation. The
intermediary model develops the boxy structure only after 2 Gyr of evolution.
The peanut shape appears 2 Gyr later and evolves slowly. The halo-dominated
model develops the boxy structure much later, around 3 Gyr, and the peanut
morphology is just incipient at the end of the simulation.

###### keywords:

galaxies: bar – galaxies: structure – galaxies: evolution – methods: numerical

††pubyear: 2015††pagerange: Orbital Structure Evolution in Self-Consistent
N-body Simulations–B

##  1 Introduction

In the context of the cosmological evolution of galaxies, baryonic matter
disks assemble at the center of dark matter halos. In the absence of large
perturbations, these disks undergo secular evolution, giving rise to several
exciting structures, such as grand design spiral arms, flocculent spiral
patterns, bars, lenses, and rings. Bars significantly impact the dynamics of
the disks as they strongly alter the axisymmetric potential of the original
disks (Athanassoula, 2013).

In the local Universe, observers have noted that two-thirds of disk galaxies
display bars (Buta et al., 2015). Optical surveys reveal that approximately 30
percent of these galactic bars are strong. If we also consider weak bars
visible in the Fourier decomposition of the light distribution, the fraction
increases to 50 percent or more (Marinova & Jogee, 2007). This fraction even
rises to $60-70$ percent when observing in the near-infrared (Eskridge et al.,
2000; Aguerri et al., 2009). Furthermore, numerical simulations have
demonstrated that bars grow in the vertical direction by a buckling
instability, which then settles down into a boxy or peanut-like (B/P) shape
when viewed edge-on (e.g. Raha et al., 1991; Merritt & Sellwood, 1994;
Debattista et al., 2004; Martinez-Valpuesta & Shlosman, 2004; Debattista et
al., 2006; Martinez-Valpuesta et al., 2006; Saha et al., 2013). These
structures are certainly found in observations (Lütticke et al., 2000; Erwin &
Debattista, 2017). Yet barred galaxies show a pronounced B/P-shape (e.g., SDSS
image of NGC 128), other galaxies do not display any boxy or peanut-like shape
or X-shaped structure in their isophotes. These face-on bars only demonstrate
elongated elliptical isophotes, e.g., IC 5240 (Laurikainen & Salo, 2017). In
this regard, Li et al. (2017) conducted a statistical study of bar galaxies
based on the Carnegie-Irvine Galaxy Survey (CGS), classifying such bars as
unbuckled bars, B/P-shaped bars, and bars with barlenses, which are buckled
bars. The study found that unbuckled bars do not exhibit any distinct
features, except for elongated elliptical isophotes and are predominantly
found in later-type galaxies, while buckled bars develop more frequently in
early-type galaxies with a massive disk (Li et al., 2017; Erwin & Debattista,
2017).

Researchers have studied bar formation in disk galaxies for several decades.
Essentially, a bar-like structure is triggered by two processes: the first one
stems from internal causes, such as dynamical instabilities within individual
galaxies, and the second one is influenced by external (e.g., tidal)
influences. Lynden-Bell (1996) and Polyachenko & Polyachenko (2003) reviewed
different mechanisms for bar formation that are all connected by the exchange
of angular momentum. In numerical simulations, the evolution of bars in N-body
models that are unstable during their formation depends on the distribution of
stellar orbits. The different components of a galaxy (halo, bulge, and disk)
interact with each other and cause the angular momentum exchange between disk
resonances and their other components (Athanassoula (2002) and references
therein).

N-body numerical simulations show that the vertical buckling instability
influences the secular growth of stellar bars (Combes & Sanders, 1981; Raha et
al., 1991; Martinez-Valpuesta et al., 2006). Martinez-Valpuesta et al. (2006)
showed the vertical structure of the bar buckles twice during the simulation.
The first buckling appears in the central area of the bar, and the second
buckling comes into view on the peripheral part of the bar. In that respect,
Martinez-Valpuesta et al. (2006) analyzed the self-consistent evolution of the
stellar bar to understand the three-dimensional structure by studying periodic
orbits in arbitrary gravitational potentials. They chose periodic orbits,
which characterize the overall orbital structure of the bar phase space, and
showed that the general shape evolves from a peanut/boxy to a vertically
asymmetric shape and then to a peanut/X-shape.

Those studies show that the general morphology of a barred galaxy changes
during its evolution; therefore, the bar potential changes as well. So,
different orbital families, which are the backbone of a barred structure,
change their importance during the evolution of the bar. Smirnov & Sotnikova
(2018, 2019) demonstrate a relationship between the face-on and edge-on
morphology of bars and the parameters of the underlying galaxy. Additionally,
in their study, Smirnov et al. (2021) showed boxy/peanut bars are supported by
different families of orbits. In the early work by Contopoulos &
Papayannopoulos (1980), they showcased the basic families of periodic orbits
($x_{1}$, $x_{2}$, and $x_{3}$) in different isochrone bar model potentials.
They showed that in weak bars, inside the inner ILR and outside the outer ILR,
the orbits follow the periodic $x_{1}$ family, and between the ILRs, most
orbits follow the $x_{2}$ family. However, if the bar is sufficiently strong,
the $x_{2}$ family of orbits does not exist. Another study by Athanassoula et
al. (1983) investigated the orbits in prolate heterogeneous bar models and
classified them into A (confined to small radii, anti-aligned), B (elongated
orbits aligned with the bar), and R (retrograde orbits). They found that in
the case of a massive bar, most quasi-periodic orbits belong around either the
B or R family. Subsequent studies by Skokos et al. (2002a, b), Patsis et al.
(2002), and Patsis & Athanassoula (2019) revealed various bifurcations of the
$x_{1}$ orbit families in 2D and 3D (known as the $x_{1}$-tree) that support
the bar. They utilized Ferrers bar potentials with different pattern speeds to
investigate these bifurcations. Similar studies have used frozen potentials of
N-body simulations, as demonstrated by Athanassoula (2013) and Patsis &
Harsoula (2018). Other studies, such as those by Wozniak (1994), Kaufmann &
Contopoulos (1996), Patsis et al. (1997), Wozniak & Pfenniger (1999), Muzzio
et al. (2005), Harsoula & Kalapotharakos (2009), Patsis et al. (2010),
Contopoulos & Harsoula (2013), Patsis & Katsanikas (2014), Tsigaridi & Patsis
(2015), and Chaves-Velasquez et al. (2017), have demonstrated that sticky-
chaotic orbits are capable of supporting galactic bars. Additionally, Patsis &
Athanassoula (2019) considered distinct types of periodic orbits, which do not
belong to $x_{1}$ family, but these could probably be suitable foundation
building blocks for bars.

Figure 1: Surface logarithm density maps for our models and their evolution.
The left, middle, and right panels show the models $A\lambda 03$, $A\lambda
04$, and $A\lambda 05$, respectively. Time increases downward. The color scale
at the bottom right of the figure corresponds to the logarithm of the surface
density ($\Sigma$ in $1.0\times 10^{10}$ M⊙/kpc2). Each panel shows the face-
on, edge-on, and end-on views, with the horizontal lines below the face-on
view indicating a distance of 10 kpc. From the figure, we can observe the bar
growth and the formation of buckling. The buckling appears at different times
in the models; around 2 Gyr for $A\lambda 03$ and 4 Gyr for $A\lambda 04$. For
$A\lambda 05$, the buckling is incipient at the end of the simulation (7 Gyr).

The computed orbital frequency of the disk particles determines the orbital
structure in barred galaxies. According to Voglis et al. (2007), not only 2:1
resonance orbits but also orbits resonances close to 3:1 and 5:2 support the
bar with regular and chaotic orbits. Furthermore, Harsoula & Kalapotharakos
(2009) presented the distribution of frequencies for the most significant
populations in various N-body simulations with frozen potentials. They
demonstrated that B/P (boxy-peanut) bulges are composed of a superposition of
several peanut shapes produced by different orbits, and peanut-shaped bulge is
primarily assembled of families of pretzel-like orbits. Additionally, Valluri
et al. (2016), and Abbott et al. (2017) showed that the fraction of orbits
associated with different families may somewhat depend on the details of the
bar potential or the parameters of the underlying galaxy (Parul et al., 2020).
Conducting a frequency analysis directly from the N-body model, Gajda et al.
(2016) found that stellar orbits have boxy shapes with significant elongation
degree, exhibiting banana or infinity symbol shapes when viewed edge-on.
Finally, Łokas (2019) suggests that the initial distortion of the bar is
likely caused by the resonant trapping of $x_{1}$ orbits, which subsequently
transition into banana-like orbits. In general, all these works highlight the
existence of different orbit families that support the bar, but remains
unclear which specific families of orbits support the bar as it evolves.

Figure 2: Slowdown of the bar pattern speed for $A\lambda 03$ (full line),
$A\lambda 04$ (dashed line), and $A\lambda 05$ (dot-dash line) models.

In this study, we show the evolution of the orbital structure in three fully
self-consistent N-body models taken from Valencia-Enríquez et al. (2019). To
assess this evolution, we measure a range of frequencies for all disk
particles at different time intervals. We examine the distribution of orbital
frequencies across the face-on, edge-on, and end-on perspectives from the
barred disk galaxy models as the bar potential undergoes changes. Section 2
provides an overview of the N-body models and the frequency analysis
methodology employed in this research. We present our findings on the
evolution of the orbital structure in Section 3. Finally, in Sections 4 and 5,
we discuss our results and draw conclusions based on our analysis.

##  2 N-Body Models and Methodology

This section includes a description of galaxy models and the methodology to
study the evolution of orbital structure through frequency analysis.

###  2.1 Disk Galaxy Models

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).

Figure 3: The plots show the ratio frequency maps for all disk particles where
the time intervals in which the frequencies were calculated increase from top
to bottom. The left, middle, and right columns represent $A\lambda 03$,
$A\lambda 04$, and $A\lambda 05$ models, respectively. The straight lines in
the plots show the location of the resonances 2:1, 5:2, and 3:1. We note for
$A\lambda 03$ that the 2:1 resonance is quickly populated, and the particles
move to larger $\omega_{R}/\omega_{z}$ on this resonance as they evolve. The
same resonance is populated later on in the $A\lambda 04$ model and at much
later times for $A\lambda 05$ model. As noted in Section 3.1, the value of
$\Omega$ may have larger errors than the other frequencies.

###  2.2 Frequency analysis method

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

$H_{x}(\omega_{x})=\left|\int_{-\infty}^{\infty}x(t)e^{i\omega_{x}t}dt\right|.$
(1)

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:

$H_{x}(\omega_{x})=\left|\dfrac{1}{N}\sum_{k=0}^{N-1}x_{k}e^{i\omega_{x}k\Delta
t}\right|$ (2)

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}$.

##  3 Orbital Structure Evolution

For a more precise frequency determination, we recalculated the simulations
$A\lambda 03$, $A\lambda 04$, and $A\lambda 05$ presented in VE19 with higher
time resolution; instead of having 612 snapshots as in VE19, the simulations
have 6120 steps to represent 6.0 Gyr in the simulations. Furthermore, we
evolve the $A\lambda 05$ model for one more Gyr because in this model, the bar
formation is triggered much later compared to the other two models (Figure 1).
By having more data points for the same time interval, we ensure that the
Nyquist critical frequency improves, and then the estimation of the
frequencies of the orbits of the disk particles are better. In this way, for
our simulations, we successfully calculated the frequencies using time
intervals of 1 Gyr. In Gajda et al. (2016), the bar is tidally induced and
presents a constant pattern speed during the time interval. The bars in our
models are evolving, and their pattern speed changes slowly (see Figure 2).
This figure shows that the slowdown for the $A\lambda 03$ model is
approximately three km/s/kpc/Gyr, while for $A\lambda 04$ and $A\lambda 05$ is
one km/s/kpc/Gyr approximately. So, the frequencies we calculated must be
understood as the most predominant frequency of a given particle within a
given time interval. In the following sections, we present our results for the
three models, covering 6 (or 7 in the case of $A\lambda 05$) time intervals.

###  3.1 Frequency evolution of all disk orbits

Figure 4: Normalized cumulative number of orbits as a function of the number
of periods of the particles for orbits within the frequency range
$1.9<w_{R}/w_{x}<2.1$.

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.

Figure 5: Histograms for the ratio $\omega_{R}/\omega_{x}$ for all disk
particles for different time intervals. Top, middle, and bottom rows show the
results for $A\lambda 03$, $A\lambda 04$, and $A\lambda 05$ models,
respectively, together with the results from the frozen potential
calculations. We notice how the resonances are populated. The $A\lambda 03$
shows a large number of particles on the 2:1 resonance and some around 5:2,
and $A\lambda 04$ follows a similar behaviour. Instead, $A\lambda 05$ displays
particles in a large range of $\omega_{R}/\omega_{x}$.

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.

###  3.2 Frequency evolution of Bar type orbits

Figure 6: Disk mass fraction of particles with frequencies showing the ratio
of $1.9<\omega_{R}/\omega_{x}<2.1$ (upper panel),
$2.4<\omega_{R}/\omega_{x}<2.6$ (middle panel), and
$2.9<\omega_{R}/\omega_{x}<3.1$ (bottom panel). The filled triangle, star, and
circle symbols represent the $A\lambda 03$, $A\lambda 04$, and $A\lambda 05$
models, respectively. The colors represent the time interval as in Figure 5.
We can notice that the $A\lambda 04$ model has the highest amount of B-type
particles between $2-4$ Gyrs, which corresponds to the bar saturation phase in
that model. At the same time, the bar formation is just starting for the
$A\lambda 05$ 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.

Figure 7: Histograms were generated to analyze the ratio
$\omega_{z}/\omega_{x}$ for orbits within the range of
$1.9<\omega_{R}/\omega_{x}<2.1$ (B-type orbits). The left column illustrates
the evolutionary behavior of orbital frequencies in the $A\lambda 03$ N-body
model. The middle column depicts frequency analysis performed on the frozen
potentials integrated throughout 1 Gyr. The right column showcases the
frequency measures integrated over a time span of 6 Gyr. In this figure, time
progresses downwards. The two prominent vertical lines in the histogram
represent the range of $\omega_{z}/\omega_{x}$ used to construct the density
maps (Figures 10, 12, 13, and 14). Figure 8: As in Figure 7, but for model
$A\lambda 04$. Figure 9: As in Figure 7, but for model $A\lambda 05$.

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 phase111Bar
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.

###  3.3 Evolution of the orbital assembly

Figure 10: Orbits assembly density profile (OADP) maps in face-on, edge-on,
and end-on views for particles with ratio $1.9<\omega_{R}/\omega_{x}<2.1$ and
$\omega_{z}/\omega_{x}$ as indicated in the bottom left of the face-on views.
The left, middle, and right panels show the models $A\lambda 03$, $A\lambda
04$, and $A\lambda 05$, respectively. The time increases downwards. Each panel
represents a different model at different time intervals and a different range
of the $\omega_{z}/\omega_{x}$ ratio (see left column of Figures 7, 8, and 9).
The x and y axes are from $-5.0$ to $5.0$ kpc, and the vertical direction z,
from $-2.5$ to $2.5$ kpc. Each map is normalized to its maximum, and six
contours are plotted (0.025, 0.05, 0.1, 0.2, 0.3 and 0.6). We can easily
follow the buckling formation in the $A\lambda 03$ model. In the beginning, it
is quite asymmetric concerning the disk plane, and later on, this structure
becomes symmetric. Instead, at the end of our simulation, $A\lambda 04$ does
not fully achieve the buckling process. In the case of $A\lambda 05$, the full
bar formation process is delayed, and the bisymmetric structure is just
incipient at the end of the simulation.

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.

###  3.4 Orbit shapes

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).

Figure 11: Evolution of five orbits in their different views from the
$A\lambda 03$ model is presented. For all the particles, the time increases to
the right in units of one Gyr. The horizontal line in the bottom left location
in the left panels represents 5.0 kpc for particles a) to d). For particle e),
the horizontal line represents 1.0 kpc. The numbers at the bottom right of
each panel are the values for $\omega_{R}/\omega_{x}$ (upper value) and the
$\omega_{z}/\omega_{x}$ (bottom value).

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.

##  4 Discussion

The present work extends the studies by Portail et al. (2015) and Gajda et al.
(2016). In Portail et al. (2015), they performed frequency analysis on orbits
calculated in frozen potentials. On the other hand, Gajda et al. (2016) used
orbits taken from a fully self-consistent N-body simulation, and the bar, in
their model, was tidally induced. They present the frequency analysis in a
single interval of around 1.6 Gyrs. During this interval, the bar in their
model exhibits a constant pattern speed.

In this work, we studied the orbital structure evolution for three N-body
models from VE19 using the dominant frequencies of the disk particles. The
models differ on the values of the initial Spin parameter of the disk
$\lambda_{d}$. We used simulations with high time resolution, and then we
could analyze the orbital frequency distribution of the bar structure in a
disk galaxy as the bar potential evolves. In the process of bar formation, the
main resonances 2:1, 5:2, and 3:1 are populated (Voglis et al., 2007). Our
disk-dominated model replenishes these resonances at earlier stages of the
model evolution, while it gets much longer for a halo-dominated model. In the
present study, we focused on B-type orbits, which fill the 2:1 resonance
($1.9<\omega_{R}/\omega_{x}<2.1$).

On the frequency ratio plots $\Omega/\omega_{z}$ vs. $\omega_{R}/\omega_{z}$
(Figure 3), we can perceive how the main resonances replenish as the models
evolve. Once the bar forms, a large number of particles occupy the 2:1
resonance, with a lower contribution of orbits with frequencies between the
5:2 and 3:1 resonances (see Figure 6). For this work, we focused our analysis
on the 2:1 resonance (B-type orbits).

We show that the distribution of B-type orbits in the $\omega_{z}/\omega_{x}$
ratio differs among all models at different time intervals. During the bar
formation in the disk-dominated model (T=[1-2]), the B-type orbits exhibit a
ratio interval of $2<\omega_{z}/\omega_{x}<3$. In the subsequent Gyrs, as the
bar reaches its mature state (boxy/peanut shape and buckling), the B-type
orbits show a ratio interval of $1.5<\omega_{z}/\omega_{x}<2.2$ (see Figures
7, 8, and 9). In the intermediate $A\lambda 04$ model, high values for
$\omega_{z}/\omega_{x}$ are achieved only in the interval T=[2-3]. The orbits
transition to lower $\omega_{z}/\omega_{x}$ much slower than the disk-
dominated model and reach the ratio interval of
$1.5<\omega_{z}/\omega_{x}<2.2$ only after the fourth Gyr. At that moment, the
orbits begin to leave the disk plane; in other words, the buckling phase
initiates. Finally, the halo-dominant $A\lambda 05$ model evolves in a
significantly different manner. The B-type orbits exhibit high values for the
ratio $\omega_{z}/\omega_{x}$ throughout the entire simulation (Figure 9). In
this case, the distribution gradually shifts to the left as the model evolves
from the fifth to the seventh Gyr. In this halo-dominated model, the buckling
of the bar is just incipient at the end of our simulation (T=7 Gyr).

As we can observe, the frequency ratios drift as the bar evolves. These shifts
are more noticeable during the buckling phase of the bar. After this event,
some orbits remain within a specific range of these frequency ratios (see the
bins in Figures 7, 8 or 9, and the orbits in Figure 11). These particular
orbits can be potentially categorized as regular or quasiperiodic orbits due
to their minimal deviations, as demonstrated by Wang et al. (2016). In that
study, the orbits were classified into regular and chaotic orbits using the
frequency drift technique, revealing that this classification depends on the
orbit projection, where the orbit appears to be more regular in the xy view
than in the other perspectives. Hence, it can infer that orbits persisting
within the same bin of the distributions over time can be classified as
regular orbits.

Furthermore, consistent findings are obtained when employing the frozen
potential. The frequency distribution exhibits a similar pattern to that
observed in the N-body simulations. However, for a more comprehensive
elucidation of the dynamics of models, an alternative approach is available
whereby orbits can be computed using a frozen potential derived from Basis
Function Expansion (e.g. Hernquist & Ostriker, 1992; Kuijken & Dubinski, 1995;
Vasiliev & Athanassoula, 2015; Lilley et al., 2018; Vasiliev, 2019; Wang et
al., 2020; Wang et al., 2022).

In a recent study by Wang et al. (2022), they demonstrated that there exist
periodic orbits exhibiting diverse pairs of morphological families with
multiplicities222Multiplicity represents the number of revolutions needed for
an orbit to close in phase space. of at least two ($\mathscr{M}=2$) or more.
Notably, certain orbits, such as the pretzel and fish orbits, exhibit a
multiplicity of two. Therefore, in N-body simulations, the orbits supporting
the bar structure can manifest distinct multiplicities of quasi-periodic
types. As a consequence, it changes as the bar potential evolves, like we
illustrate in Figure 11. For instance, orbit ‘d’ demonstrates two different
$\mathscr{M}=2$ shapes during two investigated time intervals. Furthermore, it
is plausible to consider a correlation between the multiplicity, and the
frequencies of the orbits. Thus, an extension of this study is warranted to
explore the relationship between multiplicity and orbit frequencies in greater
detail.

We also plot the orbits assembly density profile (OADP) maps to study the
morphology of the B-type orbits at different $\omega_{z}/\omega_{x}$ intervals
(Figures 10, 12, 13 and 14). In the $A\lambda 03$ and $A\lambda 04$ models,
the orbits located at the interval $1.5<\omega_{z}/\omega_{x}<2.2$ are pulled
out from the disk plane. The orbits showing lower values of
$\omega_{z}/\omega_{x}$ are located in the inner part of the disk and are
responsible for the internal boxy/X shape. The orbits having larger
$\omega_{z}/\omega_{x}$ participated in the evolution of the buckling
structure. These orbits buckle first showing a ‘smile down’ figure (Gajda et
al., 2016) and then evolve to a peanut shape (see Figure 11, orbit ‘d’). The
halo-dominated model $A\lambda 05$ develops a bar structure much later
compared to the other models, and the buckling phase is just starting at the
end of the simulation.

In Erwin & Debattista (2013), the authors present several moderately inclined
galaxies as visual examples of boxy bars. The study highlights that galaxies
such as NGC 5377 and IC5240 exhibit strong bars with broad and distinctly boxy
shapes in their isophotes, resembling the $A\lambda 03$ and $A\lambda 04$
models. Additionally, they identified other galaxies like NGC 3049 and IC 676,
which possess bars that have not vertically thickened, similar to the
$A\lambda 05$ model. Roughly speaking, based on our findings, we can conclude
that the boxy/peanut-shaped bar structures for these observed galaxies are
supported by a frequency range of $1.45<\omega_{z}/\omega_{x}<2.55$, with
smaller ratios predominantly located in the central region of the bar, while
larger quotients find in the peripheral region. Conversely, bars in galaxies
exhibiting any boxy shape sustain likely frequency ranges
$\omega_{z}/\omega_{x}>2.55$.

The higher frequencies maintained in the disk plane are probably related to
the velocity dispersions. A lower velocity dispersion indicates stronger
interactions among particles compared to those with a higher velocity
dispersion. The stronger interactions experienced by particles with low
velocity dispersion cause them to deviate from the disk plane. This
instability in the disk-spin parameter results in particles orbiting out of
the disk plane with lower $\omega_{z}$ frequency (VE19). The calculation of
the evolution of a set of models with different initial vertical velocity
dispersions should bring light to this subject.

##  5 Summary

In summary, the frequency analysis of the N-body models reveals a distinction
in the vertical behavior of the orbits. The models, which are dominated by the
disk component, catch around 30% of disk particles in the 2:1 resonance.
During the initial stages of bar formation, prior to the buckling phase, large
ratios of $\omega_{z}/\omega_{x}$ are observed. Later on in the evolution of
the models, the distribution of $\omega_{z}/\omega_{x}$ moves to lower values.
The buckling process is fully evolved in the $A\lambda 03$ model, partially
developed in the $A\lambda 04$ one, and just beginning in the $A\lambda 05$
case within our models. The main results of this study are as follows:

  * •

The orbits that form the bar structure prior to the buckling phase are within
the range of $2.25<\omega_{z}/\omega_{x}<3.05$ for the $A\lambda 03$ model,
which has a smaller scale radius, and $2.75<\omega_{z}/\omega_{x}<3.55$ for
the $A\lambda 04$ and $A\lambda 05$ models, which have a larger scale radius.

  * •

The orbits that form the bar structure after the buckling phase transition are
within the range of $1.45<\omega_{z}/\omega_{x}<2.25$ for the $A\lambda 03$
model and $1.55<\omega_{z}/\omega_{x}<2.35$ for the $A\lambda 04$ model. The
$A\lambda 05$ model shows a weak buckling in the range of
$1.75<\omega_{z}/\omega_{x}<2.55$ during the final time interval.

  * •

Lastly, we observe that after the buckling phase, the orbits with smaller
ratios of $\omega_{z}/\omega_{x}$ are located in the central region of the bar
structure, and these orbits contribute the most to the potential of the bar.
On the other hand, orbits with larger ratios of $\omega_{z}/\omega_{x}$ are
found in the peripheral regions of the bar (see Figures in Appendix B).

## DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the
corresponding author.

## Acknowledgements

We thank the referees and the scientific editor for comments and suggestions
that greatly improved the manuscript. D.V.E thanks the Facultad de Ingeniería
and Dirección de Investigaciones of Universidad Mariana, Colombia (Project
239). L.C.V. acknowledges the support of the postdoctoral Fellowship of DGAPA-
UNAM, Mexico. I.P. thanks the Mexican Foundation Conahcyt.

## References

  * Abbott et al. (2017) Abbott C. G., Valluri M., Shen J., Debattista V. P., 2017, MNRAS, 470, 1526
  * Aguerri et al. (2009) Aguerri J. A. L., Méndez-Abreu J., Corsini E. M., 2009, A&A, 495, 491
  * Athanassoula (2002) Athanassoula E., 2002, ApJ, 569, L83
  * Athanassoula (2003) Athanassoula E., 2003, MNRAS, 341, 1179
  * Athanassoula (2013) Athanassoula E., 2013, in Falcón-Barroso J., Knapen J. H., eds, , Secular Evolution of Galaxies. Cambridge University Press, Cambridge, UK, p. 305
  * Athanassoula et al. (1983) Athanassoula E., Bienayme O., Martinet L., Pfenniger D., 1983, A&A, 127, 349
  * Athanassoula et al. (2013) Athanassoula E., Machado R. E. G., Rodionov S. A., 2013, MNRAS, 429, 1949
  * Bournaud & Combes (2002) Bournaud F., Combes F., 2002, A&A, 392, 83
  * Buta et al. (2015) Buta R. J., et al., 2015, ApJS, 217, 32
  * Carpintero & Aguilar (1998) Carpintero D. D., Aguilar L. A., 1998, MNRAS, 298, 1
  * Ceverino & Klypin (2007) Ceverino D., Klypin A., 2007, MNRAS, 379, 1155
  * Chaves-Velasquez et al. (2017) Chaves-Velasquez L., Patsis P. A., Puerari I., Skokos C., Manos T., 2017, ApJ, 850, 145
  * Combes & Sanders (1981) Combes F., Sanders R. H., 1981, A&A, 96, 164
  * Contopoulos & Harsoula (2013) Contopoulos G., Harsoula M., 2013, MNRAS, 436, 1201
  * Contopoulos & Papayannopoulos (1980) Contopoulos G., Papayannopoulos T., 1980, A&A, 92, 33
  * Debattista et al. (2004) Debattista V. P., Carollo C. M., Mayer L., Moore B., 2004, ApJ, 604, L93
  * Debattista et al. (2006) Debattista V. P., Mayer L., Carollo C. M., Moore B., Wadsley J., Quinn T., 2006, ApJ, 645, 209
  * Erwin & Debattista (2013) Erwin P., Debattista V. P., 2013, Monthly Notices of the Royal Astronomical Society, 431, 3060
  * Erwin & Debattista (2017) Erwin P., Debattista V. P., 2017, MNRAS, 468, 2058
  * Eskridge et al. (2000) Eskridge P. B., et al., 2000, AJ, 119, 536
  * Fernández-Trincado et al. (2020) Fernández-Trincado J. G., Chaves-Velasquez L., Pérez-Villegas A., Vieira K., Moreno E., Ortigoza-Urdaneta M., Vega-Neme L., 2020, MNRAS, 495, 4113
  * Fernández-Trincado et al. (021a) Fernández-Trincado J. G., et al., 2021a, A&A, 647, A64
  * Fernández-Trincado et al. (021b) Fernández-Trincado J. G., et al., 2021b, ApJ, 918, L37
  * Gajda et al. (2016) Gajda G., Łokas E. L., Athanassoula E., 2016, ApJ, 830, 108
  * Harsoula & Kalapotharakos (2009) Harsoula M., Kalapotharakos C., 2009, MNRAS, 394, 1605
  * Hernquist & Ostriker (1992) Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375
  * Kaufmann & Contopoulos (1996) Kaufmann D. E., Contopoulos G., 1996, A&A, 309, 381
  * Kuijken & Dubinski (1995) Kuijken K., Dubinski J., 1995, MNRAS, 277, 1341
  * Laskar (1990) Laskar J., 1990, Icarus, 88, 266
  * Laurikainen & Salo (2017) Laurikainen E., Salo H., 2017, A&A, 598, A10
  * Li et al. (2017) Li Z.-Y., Ho L. C., Barth A. J., 2017, ApJ, 845, 87
  * Lilley et al. (2018) Lilley E. J., Sanders J. L., Evans N. W., Erkal D., 2018, MNRAS, 476, 2092
  * Łokas (2019) Łokas E. L., 2019, A&A, 629, A52
  * Lütticke et al. (2000) Lütticke R., Dettmar R. J., Pohlen M., 2000, A&AS, 145, 405
  * Lynden-Bell (1996) Lynden-Bell D., 1996, in , Barred galaxies and circumnuclear activity. Springer, pp 7–18
  * Marinova & Jogee (2007) Marinova I., Jogee S., 2007, ApJ, 659, 1176
  * Martinez-Valpuesta & Shlosman (2004) Martinez-Valpuesta I., Shlosman I., 2004, ApJ, 613, L29
  * Martinez-Valpuesta et al. (2006) Martinez-Valpuesta I., Shlosman I., Heller C., 2006, ApJ, 637, 214
  * Merritt & Sellwood (1994) Merritt D., Sellwood J. A., 1994, ApJ, 425, 551
  * Muzzio et al. (2005) Muzzio J. C., Carpintero D. D., Wachlin F. C., 2005, Celestial Mechanics and Dynamical Astronomy, 91, 173
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Oppenheim et al. (1997) Oppenheim A. V., Willsky A. S., Nawab S. H., Ding J.-J., 1997, Signals and systems. Second Edition Vol. 2, Prentice hall Upper Saddle River, NJ
  * Parul et al. (2020) Parul H. D., Smirnov A. A., Sotnikova N. Y., 2020, ApJ, 895, 12
  * Patsis & Athanassoula (2019) Patsis P. A., Athanassoula E., 2019, MNRAS, 490, 2740
  * Patsis & Harsoula (2018) Patsis P. A., Harsoula M., 2018, A&A, 612, A114
  * Patsis & Katsanikas (2014) Patsis P. A., Katsanikas M., 2014, MNRAS, 445, 3546
  * Patsis et al. (1997) Patsis P. A., Athanassoula E., Quillen A. C., 1997, ApJ, 483, 731
  * Patsis et al. (2002) Patsis P. A., Skokos C., Athanassoula E., 2002, MNRAS, 337, 578
  * Patsis et al. (2010) Patsis P. A., Kalapotharakos C., Grosbøl P., 2010, MNRAS, 408, 22
  * Polyachenko & Polyachenko (2003) Polyachenko V. L., Polyachenko E. V., 2003, Astronomy Letters, 29, 447
  * Portail et al. (2015) Portail M., Wegg C., Gerhard O., 2015, MNRAS, 450, L66
  * Press et al. (1992) Press W. H., Teukolsky S. A., Flannery B. P., Vetterling W. T., 1992, Numerical recipes in Fortran 77: volume 1, volume 1 of Fortran numerical recipes: the art of scientific computing. Cambridge university press
  * Raha et al. (1991) Raha N., Sellwood J. A., James R. A., Kahn F. D., 1991, Nature, 352, 411
  * Richmond-Navarro et al. (2017) Richmond-Navarro G., Barquero-Mena T. G., Solís-Villalta O. M., Palma-Quirós D. M., 2017, Revista Tecnología en Marcha, 30, 14
  * Saha et al. (2013) Saha K., Pfenniger D., Taam R. E., 2013, ApJ, 764, 123
  * Skokos et al. (2002a) Skokos C., Patsis P. A., Athanassoula E., 2002a, MNRAS, 333, 847
  * Skokos et al. (2002b) Skokos C., Patsis P. A., Athanassoula E., 2002b, MNRAS, 333, 861
  * Smirnov & Sotnikova (2018) Smirnov A. A., Sotnikova N. Y., 2018, MNRAS, 481, 4058
  * Smirnov & Sotnikova (2019) Smirnov A. A., Sotnikova N. Y., 2019, MNRAS, 485, 1900
  * Smirnov et al. (2021) Smirnov A. A., Tikhonenko I. S., Sotnikova N. Y., 2021, MNRAS, 502, 4689
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel & White (1999) Springel V., White S. D. M., 1999, MNRAS, 307, 162
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79
  * Tsigaridi & Patsis (2015) Tsigaridi L., Patsis P. A., 2015, MNRAS, 448, 3081
  * Valencia-Enríquez et al. (2019) Valencia-Enríquez D., Puerari I., Rodrigues I., 2019, AJ, 157, 175
  * Valluri et al. (2016) Valluri M., Shen J., Abbott C., Debattista V. P., 2016, ApJ, 818, 141
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 482, 1525
  * Vasiliev & Athanassoula (2015) Vasiliev E., Athanassoula E., 2015, MNRAS, 450, 2842
  * Voglis et al. (2007) Voglis N., Harsoula M., Contopoulos G., 2007, MNRAS, 381, 757
  * Wang et al. (2016) Wang Y., Athanassoula E., Mao S., 2016, MNRAS, 463, 3499
  * Wang et al. (2020) Wang Y., Athanassoula E., Mao S., 2020, A&A, 639, A38
  * Wang et al. (2022) Wang Y., Athanassoula E., Patsis P., Mao S., 2022, A&A, 668, A55
  * Wozniak (1994) Wozniak H., 1994, Regular Orbits and Cantori in the Potential of the Barred Galaxy NGC 936. Springer Berlin Heidelberg, p. 264, doi:10.1007/BFb0058117
  * Wozniak & Pfenniger (1999) Wozniak H., Pfenniger D., 1999, Celestial Mechanics and Dynamical Astronomy, 73, 149

##  Appendix A Frozen potential and integration

The methodology employed to determine the 3D ’frozen’ potential $V(x,y,z)$ was
based on creating a three-dimensional cubic grid that encompasses the model.
At each grid point, we calculated the potential, which results from all the
particles present. We obtained the potential at the specific position of the
particle being integrated to determine its orbital trajectory. This potential
was obtained through an interpolation process, while the force components were
derived from this interpolation. In our analysis, we operate within the co-
rotating frame of reference aligned with the bar. Utilizing a frozen
potential, our model becomes an independent Hamiltonian system, and we can
formulate its Hamiltonian as follows:

$H=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2}+\dot{z}^{2})+V(x,y,z)-\frac{1}{2}\Omega_{p}^{2}(x^{2}+y^{2})$
(3)

where $(x,y,z)$ denotes the coordinates in a Cartesian frame of reference that
rotates clockwise around the z-axis with an angular velocity of $\Omega_{p}$.
Finally, the temporal evolution was determined using the fourth order Runge-
Kutta method.

The proposed approach for three-dimensional interpolation involves the
utilization of a cubic three-dimensional grid (Richmond-Navarro et al., 2017),
see equation 4. This entails that every point in the grid is surrounded by
eight neighboring points positioned at the vertices. Each vertex represents a
known value, thereby enabling the interpolation process. The method considers
the proximity to the known points, implying that a point located closer to the
point of interest will exert a more significant influence on the resulting
interpolation outcome.

$\begin{split}A(x,y,z)=&\left(\frac{x_{1}-x}{x_{1}-x_{0}}\right)\left(\frac{y_{1}-y}{y_{1}-y_{0}}\right)\left(\frac{z_{1}-z}{z_{1}-z_{0}}\right)A(x_{0},y_{0},z_{0})\\\
&+\left(\frac{x_{1}-x}{x_{1}-x_{0}}\right)\left(\frac{y_{1}-y}{y_{1}-y_{0}}\right)\left(\frac{z-z_{0}}{z_{1}-z_{0}}\right)A(x_{0},y_{0},z_{1})\\\
&+\left(\frac{x_{1}-x}{x_{1}-x_{0}}\right)\left(\frac{y-y_{0}}{y_{1}-y_{0}}\right)\left(\frac{z_{1}-z}{z_{1}-z_{0}}\right)A(x_{0},y_{1},z_{0})\\\
&+\left(\frac{x_{1}-x}{x_{1}-x_{0}}\right)\left(\frac{y-y_{0}}{y_{1}-y_{0}}\right)\left(\frac{z-z_{0}}{z_{1}-z_{0}}\right)A(x_{0},y_{1},z_{1})\\\
&+\left(\frac{x-x_{0}}{x_{1}-x_{0}}\right)\left(\frac{y_{1}-y}{y_{1}-y_{0}}\right)\left(\frac{z_{1}-z}{z_{1}-z_{0}}\right)A(x_{1},y_{0},z_{0})\\\
&+\left(\frac{x-x_{0}}{x_{1}-x_{0}}\right)\left(\frac{y_{1}-y}{y_{1}-y_{0}}\right)\left(\frac{z-z_{0}}{z_{1}-z_{0}}\right)A(x_{1},y_{0},z_{1})\\\
&+\left(\frac{x-x_{0}}{x_{1}-x_{0}}\right)\left(\frac{y-y_{0}}{y_{1}-y_{0}}\right)\left(\frac{z_{1}-z}{z_{1}-z_{0}}\right)A(x_{1},y_{1},z_{0})\\\
&+\left(\frac{x-x_{0}}{x_{1}-x_{0}}\right)\left(\frac{y-y_{0}}{y_{1}-y_{0}}\right)\left(\frac{z-z_{0}}{z_{1}-z_{0}}\right)A(x_{1},y_{1},z_{1})\end{split}$
(4)

The initial conditions for the calculation of the orbits were taken from the
positions and velocities of the N-body disk particles at the same snapshot of
the calculated potential.

##  Appendix B Individual density maps

In this appendix, we present the density maps (see Figure 10), but now for
each bin on $\omega_{z}/\omega_{x}$ as indicated in Figure 7, 8, and 9.

Figure 12 shows orbits assembly density profile (OADP) maps in their different
views for B-type orbits from the $A\lambda 03$ model. The frequency ratio
interval ($\omega_{z}/\omega_{x}$) shows on the bottom left of each panel. We
note that before the buckling phase, all the B-type orbits in
$\omega_{z}/\omega_{x}>2.25$ contribute to the shape of the bar structure in
its initial state, which is the boxy shape (the first column, T=[1-2]). When
the bar buckles, these types of orbits no longer contribute to the form of the
bar. The orbits evolve and accommodate in the frequency interval
$1.45\leq\omega_{z}/\omega_{x}\leq 2.25$ (T=[2-3] to T=[5-6]), but the largest
contributions to the bar structure in this model come from orbits in the lower
part of this interval (see Figures 7). These orbits with low
$\omega_{z}/\omega_{x}$ are located in the internal part of the model, where
the density of the model is high. They evolve to a small and symmetric peanut
shape in the face-on and edge-on views. These orbit families all evolve in a
similar way. The evolution of orbits with larger $\omega_{z}/\omega_{x}$ are a
bit different. The size of these stellar orbits increases, and they are
responsible for the formation of the buckle. But in this case, firstly, the
orbits bend up showing a “down smile” figure. Later on, these orbits evolve to
a more symmetric figure, and the buckle becomes symmetric with respect to the
disk plane.

We present the OADP for the $A\lambda 04$ model in Figure 13. This figure
shows the OADP from 2 to 6 Gyrs. For this model, the B-type orbits that begin
to contribute to the bar structure appear in the interval
$\omega_{z}/\omega_{x}>2.75$ (the first column, T=[2-3]). Later on, the bar
evolves and the buckle begins to appear. Similar to the $A\lambda 03$ model,
the orbits with a frequency ratio larger than 1.85 participate in the buckle
structure. But here, as the bar formation was delayed compared to the
$A\lambda 03$ model, the buckle has not achieved its final configuration at
the time of the simulation.

Finally, we present the OADP for the $A\lambda 05$ model, depicted in Figure
14. In this model, significantly fewer particles exhibit
$1.9<\omega_{R}/\omega_{x}<2.1$ (refer to Figure 9). This halo-dominated model
lacks peanut shapes in the OADP maps; instead, it displays boxy and X-shaped
patterns in OADP edge-on views. Moreover, due to the dominant contribution of
the halo, the model does not reach the buckling phase during our simulation.
As we note in the former models, stellar orbits need to have a frequency ratio
$\omega_{z}/\omega_{x}$ around 2.0 for the buckling appearance. Although some
orbits show this frequency ratio, they take longer to buckle out of the disk
plane.

Figure 12: Orbits assembly density profile (OADP) in face-on, edge-on, and
end-on views for disk particles of the $A\lambda 03$ model having
$1.9<\omega_{R}/\omega_{x}<2.1$ and $\omega_{z}/\omega_{x}$ as indicated in
each panel. The time increases from left to right in units of one Gyr; the x
and y axes are from -5.0 to 5.0 kpc, while z is from -2.5 to 2.5 kpc. The
eight bins of the ratio $\omega_{z}/\omega_{x}$ shown here are those selected
in Figure 7. For this model, the bins from time interval T=[2-3] to T=[5-6]
are at the same frequencies, so we can see the time evolution of the structure
created by that particular $\omega_{z}/\omega_{x}$ bin (see text). Figure 13:
As in Figure 12, but for the $A\lambda 04$ model (see text). Figure 14: As in
Figure 12, but for the $A\lambda 05$ model (see text).