© The Authors 2024

1. Introduction

When they are viewed edge-on, the central regions of many disk galaxies exhibit a prominent structure that bulges out from the midplane of the disk and resembles either a box or, in some cases, a peanut (e.g., NGC 128, NGC 4013). These structures, usually called boxy/peanut-shaped (B/PS) bulges, are directly observed in 30%−40% of the edge-on disk galaxies in the local Universe (Lütticke et al. 2000; Erwin & Debattista 2013, 2017; Yoshino & Yamauchi 2015; Kruk et al. 2019; Marchuk et al. 2022), and some candidates have been identified up to redshift z ∼ 1 (Kruk et al. 2019). Early numerical (Combes & Sanders 1981; Combes et al. 1990; Pfenniger & Friedli 1991; Raha et al. 1991) and kinematic (Kuijken & Merrifield 1995; Bureau & Freeman 1999; Merrifield & Kuijken 1999; Veilleux et al. 1999; Chung & Bureau 2004) studies showed that B/PS bulges originate from the bars and are vertically extended parts of the bars. Combes & Sanders (1981) were the first to show that a flat bar, shortly after its appearance, begins to grow in the vertical direction and acquires the B/P shape. Extensive orbital studies of numerical models carried out later (Skokos et al. 2002; Patsis & Katsanikas 2014; Portail et al. 2015; Abbott et al. 2017; Parul et al. 2020) proved that B/PS bulges are built by the same orbits that populate the bars, that is, B/PS bulges share their building blocks with bars¹.

Although the connection between bars and B/PS bulges is well known, the debate is ongoing as to why initially flat orbits rise above the disk plane and what causes this (Sellwood & Gerhard 2020; Li et al. 2023). There are two different approaches. The first approach is based on the commensurability of three orbital frequencies (azimuthal Ω, radial κ, and vertical ω_z), which leads to resonant interaction. As a result, the orbits experience an average nonzero lifting force (Combes et al. 1990; Quillen 2002; Quillen et al. 2014). The second concept considers the emergence of instability due to the collective response of the orbits to a vertical disturbance (Raha et al. 1991). This instability is often called buckling instability in the case of a growing bar and is similar to the firehorse instability. In the case of resonant interactions, the important resonance is the vertical inner Lindblad resonance (vILR). In this region of the disk, ω_z ≈ 2(Ω − Ω_p), and Ω_p is the bar pattern speed. Combes et al. (1990) noted that the location of this particular resonance almost coincides with the location of another resonance, namely the inner Lindblad resonance (ILR, κ ≈ 2(Ω − Ω_p)), which is known to play an important role in the formation and growth of the bar. Because the locations of resonances coincide, the elongated orbits of a bar, which are initially almost completely confined to the disk plane, sooner or later rise above the disk plane, and the orbits themselves become part of the B/PS bulge.

A simple Hamiltonian model describing orbital trapping by vILR was first presented by Quillen (2002) and was based on an analogous model developed for in-plane resonance by Contopoulos (1975). The ideas of Quillen (2002) were later expanded based on the fact that the secular evolution of a bar, driven by interaction with the disk and halo particles at the corotation resonance (CR; Ω = Ω_p) leads to an increase in the bar size (Athanassoula 2003). This suggests that the resonance moves outward and elevates particles from the midplane along the path (Quillen et al. 2014). As the resonance moves farther from the centre, increasingly more orbits are lifted from the disk plane and eventually form an impressive B/PS bulge. This process can briefly be described as resonant heating or a resonant passage, when we consider the orbits in the bar reference frame.

Quillen et al. (2014) showed how resonant heating works using limited model examples. Sellwood & Gerhard (2020) studied several self-consistent N-body models and one model with an artificially imposed vertical symmetry. The authors followed the secular evolution of the orbital frequencies of a random sample of bar particles in each model and concluded that although peanut shapes formed in all models, each model individually exhibited the effect of buckling, resonant heating, and resonant trapping.

Although three distinct mechanisms have clearly been demonstrated that each contribute to the growth of the B/PS bulge, several points regarding them are still unclear. It is unclear whether each of these mechanisms operates only when certain conditions are met, and if not, whether they can operate simultaneously or if each is associated with a specific stage in the bar life. It also remains to be determined which of the two is more important for the B/PS bulge growth. It is critical to answer these questions if we wish to understand how B/PS bulges grow in real galaxies.

In this work, we demonstrate that at least the long-term resonant trapping and resonant passage (heating) that were studied in particular by Sellwood & Gerhard (2020) occur simultaneously in a typical N-body model. In addition, resonant trapping can lock the orbits for very long periods of time, up to several billion years, while the bar itself continues to evolve slowly. To demonstrate this, we analyzed the evolution of a sample of orbits that were trapped into the bar in such a model in terms of both actions and resonant angles. We used the techniques developed in our previous work (Zozulia et al. 2024). Although the resonant angles arise naturally in theoretical studies, it is an entirely new approach to employ them to describe the evolution of orbits in N-body simulations. This allows us to directly identify orbits that are captured or escape the resonance. As we show below, the commensurability of the frequencies at some point in time is clearly insufficient to state that the orbit is in resonance, and it is therefore important to study the corresponding angles. We would also like to note that by their theoretical definition, the actions (Lichtenberg & Lieberman 1992) should be constants of motion. However, when a system is perturbed by a strong bar, the actions defined for an axisymmetric potential are no longer conserved and are not really proper actions in the original sense. This arises because in this study and in some recent studies that considered the influence of the bar (Binney 2018, 2020; Trick et al. 2021; Trick 2022), the actions were calculated using the polar coordinates of particular orbits. These actions are sometimes called unperturbed actions (Binney 2018, 2020) because they correspond to actions of the angularly averaged system (i.e., an axisymmetric system). To investigate the behavior of unperturbed actions on different timescales, we time-averaged these actions over different time intervals (Sect. 2.2).

Our model contains an early buckling episode. We excluded this episode and its role in the uplift of orbits in the central regions in the early stages of bar evolution here and focused on the later stages, when buckling no longer affects the orbits that join the bar.

2. Simulations and action variables

2.1. Numerical model

To study how the orbits acquire their vertical extensions, we used the same numerical model as in Zozulia et al. (2024), where we studied the process of orbital trapping in terms of the actions and frequencies of the orbits. Here, we briefly describe the properties of the model and refer to the cited work for more details.

The model initially consisted of two physical components: a pure stellar disk (exponential) and a spherical dark halo with an Navarro–Frenk–White density profile (NFW, Navarro et al. 1996). A numerical implementation of the model was prepared using the subroutine mkgalaxy (McMillan & Dehnen 2007). The components were then evolved in a self-consistent manner for about 8 Gyr using the gyrfalcON code (Dehnen 2002), and a strong and fast bar emerged at the center of the disk at about 1 Gyr. Almost immediately after the formation of a bar, its inner part become vertically extended, that is, a B/PS bulge emerged. After the bar² formation, both the flat bar and the B/PS bulge continued to increase in size, capturing new disk orbits and subsequently lifting them out of the disk midplane (Zozulia et al. 2024).

Throughout this work, we use the natural system of units of the N-body model, in which the gravitational constant G = 1, the disk mass M_d = 1, and the initial exponential scale of the disk R_d = 1. One time unit corresponds to 13.8 Myr when the disk total mass is M_d = 5 × 10¹⁰ M_⊙ and the scale length R_d = 3.5 kpc.

2.2. Frequencies and action-angle variables

Following our previous work (Zozulia et al. 2024), we estimated three unperturbed (see Binney 2018, 2020) actions J = (J_R,  J_z,  L_z) and frequencies Ω = (κ,  ω_z,  Ω) on different timescales, namely medium-term and secular values. Medium-term values were calculated either on the scale of one radial oscillation (J_R,  κ, and L_z,  Ω) or of one oscillation along the z-axis (J_z,  ω_z). On the same timescale, we found the angles (θ_R,  θ_z,  θ_ϕ) by integrating the frequencies over time. The secular values (J and Ω) were obtained on the timescale of one oscillation of the medium-term value. Almost all quantities in this work are secular unless otherwise noted.

We explored the relation between the vertical bar structure and the vILR. To do this, we defined four different behaviors of the resonant angle θ_z − θ_R:

• 1, 2.

  ${\dot{\theta }}_z>{\dot{\theta }}_{\mathrm{R}}$ or ${\dot{\theta }}_z<{\dot{\theta }}_{\mathrm{R}}$. Circulation of the resonance angle with a positive or negative angular speed.

• 3.

  “In resonance”. This designation is used for orbits that librate near the resonance angles 0 or π (for banana orbits) and π/2 or 3π/2 (for anti-banana orbits).

• 4.

  “Passage”. This notation is used for orbits that are in the process of changing their behavior (changing the direction of circulation) and are also near the resonance (but do not become stuck in it).

In Sect. 3.3 we give specific examples of the orbits for each of the cases described above. Using our classification of orbital behavior, we not only determine whether the orbit is in resonance or in the passage, but can also estimate the time it spends in one or the other state.

A detailed description of the calculations of the action-angle variables and the determination of the different behavior of the resonant angle is presented in Appendix A.

3. Evolution of the vertical bar structure in action-frequency space

3.1. General structure of a bar in action-frequency space

In this section, we explore the structure of our model in action–frequency space at the stage of a mature bar (t = 400) that gradually grows in both directions, in-plane and vertical. We specifically focus on the vertical action J_z, which is tied to the vertical excursions of the orbits, and which changes when the orbits become trapped or heated by the resonance. We also note that by the time moment t = 400, the episode of early buckling (t ≈ 180) has long passed and does not affect the evolution of the vertical bar structure in the period under study.

The first three top plots of Fig. 1 demonstrate color-coded values of the vertical action J_z in the maps (L_z,  J′_R), (L_z,  ω_z/κ), and (2(Ω − Ω_p)/κ,  ω_z/κ), from left to right, respectively. In the three bottom plots, on similar maps as in the top plots, color gradations show the changes ΔJ_z over the time period from t = 300 to t = 400. Areas with an average increase in J_z are marked in red, while areas in which J_z decreases are shown in blue. The displacements (flows) of the orbits are presented in the form of a vector field in the same maps. These displacements were calculated based on the changes in the density of the orbits in a pixel over a period of time from t = 400 to t = 405, and they indicate where the ensemble of orbits tends to move in phase space over short periods of time.

Fig. 1. Vertical action and its evolution. Left and middle columns: 2D maps of the vertical action Jz plotted in different coordinates, (Lz, JR′=JR − θ(−Lz)Lz) (first column), (Lz, ωz/κ) (second column), and (2(Ω − Ωp)/κ, ωz/κ) (third column). Top row: Average value of the action in the pixel. Bottom row: Mean changes in Jz from t = 300 to t = 400 (actions and frequencies refer to the moment t = 400). The vector fields in the bottom rows demonstrate the average flow of orbits in given planes over five time units (from 400 to 405). The thin gray lines in each plot represent isodensity contours. The thick gray lines in the top left plot correspond to the lines with the same value of the adiabatic invariant Jv = JR′+Jz + Lz/2. The vertical green lines in the middle column indicate the location of the CR and ILR. Right column: Snapshots (xy and xz) of the bar orbits (abnormal orbits; see Zozulia et al. 2024) at t = 400 with different frequency ratios ωz/κ > 1 (top two plots) and ωz/κ < 1 (bottom two plots).

In the top left plot in Fig. 1, most bar orbits are located below an outline resembling an inverted parabola and have 0 ≲ L_z ≲ 1.5. On the right, the parabola is limited by the envelope that corresponds to the constant value of the adiabatic invariant J_v = 0.8 (one of the gray lines). Flat orbits with low values of J_z are located here (dark blue and black in the color bar). These orbits have the highest values of the adiabatic invariant J_v of the orbits that only recently became trapped into the bar³. Orbits that joined the bar earlier (with lower values of J_v) have already increased J_z. In the bottom left plot, we show in the same coordinates that the orbits move toward lower values of L_z and higher values of J′_R in general, and their J_z increases.

The encounter of the orbit with the vILR leads to an increase in J_z. In terms of frequencies, the vertical resonance condition⁴ for the bar orbits can be written as ω_z/κ = 1. The second upper plot in Fig. 1 shows that the orbits in the region 0 ≲ L_z ≲ 1.5 with high values of J_z⁵ (purple and orange) lie in the map below the horizontal line ω_z/κ = 1, that is, they have already passed through the vILR. The orbital flows are generally directed below the horizontal line ω_z/κ = 1 (bottom plot in the same phase space).

In the third top map in Fig. 1, most of the orbits in the bar are located along the vertical stripe near the line 2(Ω − Ω_p)/κ = 1. Orbits with the highest J_z have ω_z/κ < 1. For these, the encounter with the vILR has already taken place. Flat orbits (low values of J_z) trapped into a bar have ω_z/κ > 1. Some of the orbits have ω_z/κ = 1 and appear to be in both in-plane and vertical resonances simultaneously. However, based on this fact alone, we cannot conclude whether they are stuck in the resonance for a long time or are in the process of passage. The most interesting changes are visible in the vector map (2(Ω − Ω_p)/κ,  ω_z/κ). Here, the flows are generally directed from left to right toward the line 2(Ω − Ω_p)/κ = 1, and then the arrows indicate a descent to the vILR (ω_z/κ = 1) and immediately below it, where the region with the greatest change in J_z is located.

The bar orbits that are in the process of passage through the vILR (ω_z/κ > 1), and those that have already passed through it (ω_z/κ < 1), are assembled into morphologically different structures (Fig. 1, fourth column of the four snapshots). The latter orbits constitute a more compact structure in the in-plane projection and have a pronounced B/P shape in the vertical direction. Orbits with ω_z/κ > 1 join the bar later. They are more extended and have not yet had time to elevate significantly. Their B/P shape is weakly expressed.

It should be noted that the maps of the J_z distributions in a particular time moment in Fig. 1 are the result of the action of all possible mechanisms of orbital elevation over a long time interval. These mechanisms can operate on different timescales. In Sect. 3.2 we analyze the fate of initially flat orbits (those with J_z(t = 300) < 0.05) that were only recently trapped into the bar to distinguish between the mechanisms of vertical elevation of orbits in momento.

3.2. Trapped orbits in action-frequency space

Figure 2 (the first three maps in each row) shows the density of all orbits by yellow-green gradations in color in a pixel for the same phase coordinates as in Fig. 1. We overlaid the location of the trapped orbits in each map, and the black-orange-yellow color gradations (as in Fig. 1) correspond to J_z for these orbits. The maps are given for three time moments. t = 300: The orbits in question have not yet fallen into the trap of the bar. On the right, they are limited by the isoline J_z ≈ 0.8, and in the bar itself, there are no orbits yet with these values of J_z (first plot on the left). They are also enclosed left of the vertical line 2(Ω − Ω_p)/κ = 1 in a black region (flat orbits), outlined by light blue isodenses (third plot from the left). t = 400: The orbits have already been captured by a bar. t = 500: The orbits have been in the bar for a long time and continue to evolve.

Fig. 2. Vertical action and its evolution for captured orbits. Left and middle columns: 2D density maps (bar and disk) in the background and 2D maps (superimposed in the foreground) of Jz for orbits that were flat (Jz < 0.05) at t = 300 and trapped into the bar over the time interval of t = 300 − 400. Different time moments are shown from top to bottom, t = 300 (top row), 400 (middle row), and 500 (bottom row). The thin gray lines in each plot represent the isodensity contours of trapped orbits. In the left column, the thick gray lines correspond to isolines of the same value of the adiabatic invariant Jv = JR′+Jz + Lz/2, and the thick cyan line corresponds to the frequency ratio isoline ωz/κ = 1. The vertical blue lines in the middle column indicate the location of the CR and ILR. Right column: Snapshots (xy and xz) of the trapped orbits at the corresponding time moments.

The changes in the phase maps (Fig. 2) correspond to the general trends noted in Sect. 3.1. The orbits of the entire bar and the orbits trapped into the bar over the time interval t = 300 − 400 both increase in J′_R (first column). The area overlaid on the maps, which only consists of trapped orbits, moves toward lower L_z values (first and second columns). With increasing J_z, the trapped orbits move below the line ω_z/κ = 1 (in the first column of Fig. 2, this line is plotted in blue). The plots in the third column show that the orbits gradually propagate toward the in-plane resonance (from left to right) and then descend along a vertical line 2(Ω − Ω_p)/κ = 1, increase in J_z, with the highest values of J_z achieved for ω/κ < 1 (t = 500).

The change in the morphology of the entire ensemble of trapped orbits is visible in the plots in the rightmost column of Fig. 2. The edge-on snapshot for t = 400 closely resembles the edge-on snapshot in Fig. 1 for ω_z/κ > 1, which corresponds to orbits that mostly have not yet entered into vILR. Although these orbits have increased J_z, they do not protrude too much from the midplane.

Since ω_z/κ = 1 is a condition for vILR, the orbits with ω_z/κ < 1 are those that have already passed through the vILR. At the same time, the blue isodenses at t = 500 (Fig. 2, third column) show that pixels with the highest density for the entire ensemble of captured orbits are not located in the area with the highest J_z values, but rather fall in the ω_z/κ = 1 region. This is the first indication that even at the stage of a mature bar, there may be many orbits among the orbits that newly joined the bar that are directly trapped in vILR. In the next section, we study orbital trapping in vILR in more detail based on the long-term behavior of the resonant angle rather than based on the frequencies. We classify the four main types of orbits with respect to the vILR, and we provide statistics on orbital types.

3.3. Orbital classification with respect to the vILR

In Sect. 2.2 we described how four different types of orbital behavior with respect to the vILR can be described. Considering the vILR and the resonant angle θ_res = θ_z − θ_R, these four types correspond to (1) circulations with a monotonic increase in θ_res; (2) the same, but with monotonically decreasing θ_res; (3) libration around the resonant angles θ_res = 0 or π (banana orbits) and π/2 or 3π/2, (anti-banana orbits); and (4) the passage through the resonance (θ_res changes sign in a particular manner). By studying the behavior of the resonant angle, it is possible (a) to determine the status of the orbit at any point in time by tracking the history of its resonant angle back and forth in time; (b) to analyze how many orbits are in one state or another in relation to vILR at a given time; and (c) to determine the time the orbit stays in resonance or in passage throughout the entire interval under consideration, and to give the statistics of this parameter.

To illustrate how the process of trapping in vILR and passage through it occurs in a realistic N-body galactic model, we determined the mode of the resonant angle behavior in the time interval t = 250 − 550, that is, after buckling. We did this for bar orbits and for recently trapped orbits.

Figure 3 gives specific examples illustrating the orbital behavior described above. For each orbit, we plot the projections on the xy and xz planes, the time evolution of the resonant angle θ_res modulo 2π, and the evolution of the orbit in the plane (θ_res,  J_z)⁶. The orbit, evolving, can change the modes of behavior of the resonant angle. Two upper and lower left sets of plots, each consisting of three panels, show three orbits that are trapped into vILR. As the galactic potential evolves, the orbits in the upper subplots escape from resonance and go into circulation mode (green), while their vertical action increases significantly. The orbit in the bottom left subplot shows a different behavior. It continues to librate around the resonance until the end of the simulation time, increasing and decreasing J_z while maintaining the average value of J_z. In bottom right plot, we plot the orbit that rapidly passes through vILR without being locked into libration mode and leaves the disk midplane, increasing its vertical action J_z.

Fig. 3. Four examples of different orbital behaviors in our N-body simulations. For each orbit, from left to right, we show snapshots (xy and xz), the change in the resonant angle θz − θR with time, the orbit in the plane (θres,  Jz), where Jz is the midterm vertical action. The top left and right orbits remain a long time in vILR and librate around θz − θR ∼ π (banana orbit) and π/2 (anti-banana orbit), respectively, and then leave it. The bottom right orbit is trapped into vILR and also librates around θz − θR ∼ π until the end of the integration time. The bottom left orbit passes through the resonance around θz − θR = π. Orange and green in each plot correspond to the time when θz − θR increases and decreases, respectively. Blue and purple correspond to the stage when the orbit is in resonance or passes through it, respectively.

The orbits shown in Fig. 3 are taken as an example. They are trapped into the bar at t ≈ 200 − 250, that is, after the buckling episode. In the time interval 250 − 550, the average value of the adiabatic invariant for them is close to 0.65−0.71. While they are flat, they are practically indistinguishable in morphology and extent along the x-axis, but their evolution in the vertical direction is different.

Thus, in the N-body model that experienced buckling in an early stage of evolution, at a later stage we observe (at the level of individual orbits) the evolution of the vertical structure due to other mechanisms: trapping into resonance, long-term trapping and subsequent exit from resonance, and rapid passage through it.

Further, we considered the general structure assembled from the bar orbits and orbits trapped into the bar in t = 300 − 400 (the same as in Sect. 3) on the xy and xz planes. We distinguished all orbits according to their behavior with respect to the vILR at time moment t = 500. To determine the type of the corresponding orbital behavior at t = 500, we do not need to know the entire history of the orbital evolution in the phase plane in the interval t = 250 − 550. We track the orbit behavior over a certain period of time backward in time and, if necessary, forward, up to t = 550. Figure 4 shows snapshots for each orbit type identified at time t = 500. The top two rows refer to all orbits in the bar, and the two lower rows only show recently trapped orbits.

1. Orbits that only increase their angle θ_res are collected into a flat bar (Fig. 4, second column) consisting of classical x1 and elliptical orbits. This is especially noticeable for orbits that are recently trapped into the bar. There are quite a few orbits with increasing θ_res in a mature bar (9%), but they are most numerous among the recently trapped orbits (35%) and among all orbits with increasing θ_res (trapped orbits of this type make up almost a quarter). These orbits continue to actively replenish the previously formed flat part of the bar. They have not yet entered the vILR, but they may lift out of the disk plane and enter the vILR later.

2. Orbits in the vILR that are trapped into the bar at the time interval t = 300 − 400 account for 15% of all orbits of this type in a bar that continue to be in the vILR for more than one libration period (Fig. 4, third column) at time t = 500. Trapped orbits that are stuck in the vILR continue to participate in building this part of the B/PS bulge, which has an X-shaped structure that extends along the radius. Following the shape of the banana orbits, the distribution itself resembles two intersecting bananas. Among the orbits trapped into a bar in the interval t = 300 − 400, this is the second most numerous type of orbit (30%), although among all orbits in the bar, this type accounts for 12% of the cases.

3. Orbits in passage (Fig. 4, fourth column) generally repeat the distribution of the resonant orbits. Near the moment t = 500, they experience an episode of rapid passage through the vILR without turning around on the plane (θ_res,  J_z), but with a sharp increase in J_z. The bar and the ensemble of trapped orbits contain about 5% and 21% of these orbits, respectively.

4. Orbits with a decreasing resonant angle (Fig. 4, fifth column), that is, those that are already out of the vILR, are the most numerous type of orbits included in the B/PS bulge (74%), and in its thickest part. In the past, they may have been in resonance for some time or may have quickyl passed through the vILR. We do not exclude the possibility that for some orbits of this type, the increase in J_z was associated with an episode of buckling at t ≈ 180. However, this does not apply to orbits trapped into the bar in t = 300 − 400. The number of recently trapped orbits among orbits of this type is relatively small (14%), but they contribute precisely to the thickest and most compact part of the B/PS bulge. The structure into which they assemble is morphologically similar to the structure assembled from all orbits of this type.

Fig. 4. Snapshots (xy and xz) of the bar (top two rows) and the structure assembled from the orbits trapped into the bar in the time interval t = 300 − 400 (bottom two rows) at t = 500. The first column shows all orbits under consideration. Structures assembled from orbits that increased and decreased their resonant angle θz − θR are displayed in the second and fifth columns, respectively. The third and forth columns show orbits in the vILR and passing through it, respectively.

While the orbit is librating, J_z oscillates around the mean value. Leaving the vILR, the orbit enters the circulation mode with the maximum value of J_z (Fig. 3). Thus, the thickest part of the B/PS bulge grows due to orbits that have passed through the vILR and begin to circulate, decreasing the resonant angle. However, orbits that have recently been trapped into the mature bar and entered circulation mode with a decreasing resonant angle are few, and they do not effectively replenish the thick parts of the bar. Among all the orbits of a bar of this type, they make up 1%. In 100 − 200 time units, about 2 ⋅ 10⁴ of 14 ⋅ 10⁴ orbits were added to the B/PS bulge in this manner. About half of the newly trapped orbits are still in or pass trough resonance (about 7 ⋅ 10⁴ of 14 ⋅ 10⁴). From the given statistics, it follows that the inner thickest part of the bar was formed long ago, in contrast to the thinner extended parts, which are effectively replenished by resonant orbits.

3.4. Evolution of the orbital types

In Sect. 3.3 all orbits trapped into the bar from t = 300 to t = 400 were classified into four types with respect to the vILR at t = 500. We present the change in ω_z/κ and J_z for these types from t = 300 to t = 500. Figure 5 shows the 1D distribution over these quantities⁷. All types of orbits on average reduce the ratio ω_z/κ and increase J_z.

Fig. 5. 1D distributions of the orbits trapped into the bar over the time interval t = 300 − 400, which were flat (Jz < 0.05) at t = 300 (the same orbits as in Fig. 2): The three left plots display the secular frequency ratio ωz/κ at t = 300, 400, and 500, respectively. The two right plots correspond to the secular vertical action Jz at t = 400 and t = 500. The gray lines show the distributions of all trapped orbits, and the other colors indicate orbits with different behaviors (the same as in Fig. 3) at t = 500.

Orbits that increase their resonant angle indeed have ω_z/κ > 1 and are flat in all time moments (orange lines in Fig. 5). For these, J_z < 0.1 even at t = 500 with an explicit peak at J_z = 0. In turn, orbits that decrease their angle at t = 500 have a frequency ratio ω_z/κ < 1 and the widest distribution over J_z, spanning the range up to J_z ≈ 0.3 and beyond (green lines in Fig. 5).

Orbits in resonance and passage through it (in Fig. 5, blue and violet lines, respectively) are distributed around ω_z/κ ≈ 1 at t = 500, but because these orbits librate, the distributions are smeared, and both distributions are poorly distinguishable in ω_z/κ. For this reason, almost resonant orbits are more accurately identified by the evolution of their resonant angle and not by their frequency ratio⁸. Both orbit types have a wide range of secular J_z from 0 to 0.25, with an average value of approximately 0.15.

We considered the entire history of the orbital evolution and determined how long orbits stay near resonance in the time interval from t = 250 to t = 550 (the stage of a mature bar; after buckling). The distributions over this time (Δt) are shown in Fig. 6 separately for all orbits in the bar and for orbits trapped in the time interval between t = 300 and t = 400. The classification of orbits by type (passage through resonance and libration near resonance) is made at t = 500. We also distinguish between banana orbits (θ_z − θ_R ∼ π in Fig. 6) and anti-banana orbits (θ_z − θ_R ∼ π/2 in Fig. 6). Anti-banana orbits are unstable (Pfenniger & Friedli 1991), and we expect Δt for them to be shorter than for banana orbits. They also leave resonance and enter the circulation mode earlier.

Fig. 6. Distributions of the time that the bar (four left plots) and trapped (four right plots) orbits spend in resonance or pass through it, around the angles 0,  π, and π/2,  3π/2. Orbits with resonant angles of 0,  π are combined and are designated in the plots as θz − θR ∼ π. The same holds for the angles π/2,  3π/2. They are denoted as θz − θR ∼ π/2. The number of each type of orbit is indicated in the top right corner of each plot. The resonant and passing orbits are identified at t = 500, while the time interval t = 250 − 550 is being considered. The orange distributions correspond to orbits that enter and leave resonance or pass through it during this time period. The blue part of distributions corresponds to orbits that at the time t = 250 were already in resonance and/or were still in resonance at moment t = 550, and thus, they show a lower estimate for the time spent near the vILR.

The orange part of distributions corresponds to the orbits that enter the resonant area (vILR) after t = 250 and leave it after t = 500, but before t = 550, that is, we exactly know the time of entry and exit from resonance. These orbits cannot spend more than 300 time units of time near resonance.

The blue part represents the remaining orbits. These orbits could be already in resonance at t = 250 and/or were still in resonance at t = 550. For orbits trapped into the bar in the time interval 300 − 400, we can determine the moment of the vertical resonance encounter, but we cannot determine the moment of exit from vILR after t = 550. It follows from the above that for all bar orbits and for orbits that were recently captured into the bar, the blue distributions show a lower estimate for the time spent near resonance because we do not analyze the model before t = 250 for all bar orbits and after t = 550 and cannot fix the moment of entry and/or exit from resonance. The sum of the orange and blue distributions gives all orbits in the bar or all trapped orbits (t = 300 − 400).

The plots show that near the resonant angle ∼π/2 (3π/2), the orbits spend little time (about 50 time units, or about 0.7 Gyr). In contrast, orbits that are in resonance near ∼0 (π) are stable banana orbits and can be in a state of long-term libration around the resonance. Half of all bar orbits with a resonant angle ∼0 (π) spend more than 180 time units in vertical resonance over a time interval from t = 250 to t = 550. Around 3 ⋅ 10⁴ orbits in the bar spend more than 300 time units in resonance, or more than 4 Gyr, and some of them could have entered this mode from a buckling episode. However, in our simulations, ‘recently’ trapped orbits cannot spend more than 200 − 250 time units (2.8−3.5 Gyr) in resonance, because before time t = 300, they were not in the bar and had a low J_z value. Nevertheless, half of trapped orbits stay in the vILR for more than 135 time units (1.9 Gyr), that is, we observe them in the vILR most of the time.

4. Discussion

An important preliminary step in describing resonance phenomena in a stellar disk is to use the frequency characteristics of the disk itself, which are determined by the model potential. However, the dynamic properties of individual orbits and their evolution in a changing potential can hardly be described within the framework of this approach. This is the general idea that we continue to develop in this work.

Ceverino & Klypin (2007) were perhaps the first to clearly demonstrate with examples of completely self-consistent models without imposed restrictions in the form of a frozen potential that resonances in the stellar disk (ILR, CR, and OLR) are not special regions in the disk, but collections of orbits with certain frequency characteristics. Parul et al. (2020) used a similar approach to show that the brightest part of BP/S bulges, the so-called X-structures, is built from orbits with different spectral characteristics that are determined directly in a self-consistent N-body potential.

Having supplemented frequencies with action variables, Zozulia et al. (2024) introduced the concept of orbital mode, normal and abnormal, by analogy with Lynden-Bell (1979) normal and abnormal regions. Zozulia et al. (2024) concluded that an orbit is trapped into a bar when the orbital mode changes from normal to abnormal, that is, when the angular momentum of the orbit and the rate of orbital precession begin to change synchronously, and not simply when the particle crosses the line of zero Lynden-Bell (1979) derivative.

Finally, we included into the analysis of self-consistent N-body simulations the concept of a resonant angle in addition to actions. The behavior of the resonant angle enables us to describe the vILR more accurately than via the commensurability of frequencies. We showed that the orbit that becomes stuck in or passes through the vILR is also a characteristic of the orbit itself, and not a consequence of its falling into a special region in the disk. Even at late stages of the bar evolution, the orbit can enter a long-term libration regime in the (θ_res,  J_z) plane. A significant fraction of orbits retain this behavior from the early stages of bar formation. Our findings are consistent with the results by Parul et al. (2020). In the models considered there, banana-shaped orbits (resonant orbits) occurred in a wide range of the Jacobi integral values, that is, both in the inner areas of the bar (formed earlier) and in the outer areas (joined to the bar later).

When we qualitatively compare our results with the results of Hamiltonian analysis in Quillen (2002), Quillen et al. (2014), individual particles indeed show the behavior described in both these works, that is, 2:1 orbits that librate around the resonance and orbits that pass through it. However, judging by the results of our study, the vILR moves and traps new orbits not along the disk, but in the space (J_R,  L_z,  J_z), as in Fig. 2 (left column, blue line). Sellwood & Gerhard (2020) gave an example of an orbit in resonance only for the case with artificially imposed vertical symmetry (their Fig. 13, upper plots). In its behavior, this orbit is analogous to our two upper orbits in Fig. 3. At the same time, we did not impose any restrictions on the model, but our orbits show a similar behavior at later stages, and most surprisingly, despite a long stay in resonance, the orbit can leave it (Fig. 3, upper left plots). The orbits in Fig. 3, including the orbits in passage, have similar values of the adiabatic invariant J_v (0.6−0.7) and similar initial values of J_z, close to zero, but their fate is different.

5. Conclusions

Considering the ensemble of orbits in the bar at the late stage of B/PS bulge building, we came to the following conclusions:

• We identified the orbits that are being trapped into the bar at a particular time interval and tracked their evolution both in the plane and in the vertical direction. We studied the properties of corresponding orbits, including the action variables, the frequency ratio ω_z/κ, and the resonant angle θ_res = θ_z − θ_R. An analysis of the behavior of the resonant angle for all orbits in the N-body model on a long timescale was made for the first time. This allowed us to determine the various mechanisms whose action leads to the vertical bar shape in the stage of a mature bar, when the bar experiences steady growth.

• We found that some orbits quickly pass through resonance (resonance heating) while others remain in the vILR for a long time. Surprisingly, the latter mechanism appears to be dominant in lifting the orbit from the disk midplane. Both types of interactions are observed for a significant number of orbits. Thus, we conclude that both mechanisms work simultaneously at least in our model. Since at the stage of a mature bar the episode of early buckling is long past, we conclude that orbit lifting arises only from resonant interaction. Because the properties of the orbits were calculated in a self-consistent model without imposed symmetry or in a frozen potential (Sellwood & Gerhard 2020), we assume that both mechanisms work in real galaxies.

• The orbits spend different times in resonance. The time of the resonance passage is very short for all orbit types (0.7 Gyr on average). Banana-like orbits can become stuck in the vILR for several billion years, however. A small group of orbits resides more than 4 Gyr in the vILR, but they do exist even though the properties of the bar change, that is, the bar slows down.

• Although our model includes a significant fraction of orbits in vILR in the late stages of bar evolution (Fig. 4), the thickest part of the bar is dominated by orbits in circulation mode with decreasing θ_res, that is, orbits that have left the vILR with a significant increase in J_z (Fig. 5).

• In terms of their dynamical properties (the value of the adiabatic invariant, the initial vertical action, and the frequencies), the orbits trapped in the vILR and passing through it are not practically different. However, they choose their own path, that is, they evolve differently.

Acknowledgments

We acknowledge financial support from the Russian Science Foundation, grant no. 24-22-00376. We thank the anonymous referee for his/her thorough review and highly appreciate the very useful comments and suggestions that contributed to improving the quality of the article. We also acknowledge the use of the AGAMA (Vasiliev 2019) and mpsplines (Ruiz-Arias 2022) python packages, without which the present work would not be possible.

References

Appendix A: Action-angle variables, the resonances

We estimate three unperturbed actions J = (J_R,  J_z,  L_z) and frequencies Ω = (κ,  ω_z,  Ω) on different time scales using a subroutine similar to the one we used in our previous work (Zozulia et al. 2024). First of all, the instantaneous actions and frequencies (calculated via AGAMA’s package subroutines, Vasiliev 2019) experience periodic short-term fluctuations due to the influence of the bar on time scales of the order of one radial or vertical oscillation. We calculate the actions on these timescales by averaging them along the orbit from apocenter to apocenter (J_R and L_z) or from one z-maximum/minimum to another to obtain J_z^max and J_z^min, respectively. The frequencies κ, ω_z^max, and ω_z^min are defined as the reciprocal of the time between two adjacent apocentres (T_R) or two z-maxima/minima (T_z^max and T_z^min), multiplied by 2π⁹. Mean angular velocity Ω can be obtained through the increase in the polar angle calculated between two adjacent apocentres Δϕ/T_R. After this, we have six piecewise functions (four for actions and four for frequencies) that can be smoothed using mean-preserving spline (Ruiz-Arias 2022). This procedure allows one to obtain middle-term actions J(t) and frequencies Ω(t) (see Zozulia et al. 2024 for details). Knowing Ω(t), we also calculate the dependence of angles θ = (θ_R,  θ_z^max,  θ_z^min,  θ_ϕ) on time, integrating expression for frequencies θ = ∫_t0^tΩ(t)dt + θ₀, where θ₀ is the angle phase at the time t₀. We use the fact that at the time of passage of the first apocenter θ_R(t₀^apo) = 0 and θ_ϕ(t₀^apo) = ϕ₀ (ϕ₀ is the polar angle of the first apocentre). Also at the time of passage of the first z-maxima/minima θ_z^max(t₀^z_max) = 0 and θ_z^min(t₀^z_min) =  ± π. It should be noted that the use of mean-preserving spline (Ruiz-Arias 2022) allows one to preserve the angle between two apocentres equal to 2π for θ_R, Δϕ for θ_ϕ, and the angle between the two z-maxima/minima equal to θ_z. We find the vertical variables as follows: J_z = (J_z^max + J_z^min)/2, ω_z = (ω_z^max + ω_z^min)/2 and θ_z = (θ_z^max + θ_z^min)/2.

In addition to the middle-term oscillations, there are also secular orbital changes. To study them, we average middle-term actions and frequencies over the corresponding oscillation period (Zozulia et al. 2024). In this case, we repeated the averaging procedure, but it is carried out between the maximum and minimum values of all quantities (Ω and J).

In this paper, we consider the adiabatic invariant of the following type J_v = J_R′+J_z + L_z/2, where J_R′=J_R − L_zθ(−L_z) (θ(x) is the Heaviside step function). The introduction of J_R′ variable is due to the fact, that the precession rate for flat ILR orbits with negative Ω (Δϕ < 0) is equal to Ω_p = Ω + κ/2, and is equal to Ω − κ/2 for orbits with positive one. With such a frequency modification, the adiabatic invariant for the flat orbits takes the form J_f = J_R − L_z/2.

When an orbit is trapped into resonance, it starts to librate around some average angle (Lichtenberg & Lieberman 1992; Binney 2020), or near some stable point in the phase space. Let θ_res be the resonant angle. Since we explore the relationship between vertical bar structure and vILR in this work, we define the resonant angle θ_res as θ_z − θ_R. Without loss of generality, we can select this angle in such a way that the stable point corresponds to θ_res = 0 (or π), stable banana orbits (the orbit’s tips are pointing downwards and upwards, respectively), or π/2 (or 3π/2), unstable anti-banana orbits, (Pfenniger & Friedli 1991).

Orbits can exhibit a rich variety of behavior. We were able to clearly identify four types of characteristic behavior. In practice, to identify resonant orbits, we assume that such an orbit can not leave the area between two adjacent stable angles, 0 ± π/2 (π ± π/2) or π/2 ± π/2 (3π/2 ± π/2), respectively.

• 1, 2.

  The first two types of behavior are associated with a resonant angle circulation. If the resonant angle θ_res consistently takes values that are multiples of π/2 three times in ascending or descending order, it corresponds to the circulation of the resonant angle with positive ${\dot{\theta }}_{\mathrm{res}}>0$ or negative ${\dot{\theta }}_{\mathrm{res}}<0$ angular velocities, respectively.

• The orbit librates near the stable point in the phase space. In this case, the orbit’s resonant angle θ_res crosses some angle multiple of π/2 at least 3 time without going beyond ±π/2, we mark such an orbit as “in resonance” during the period of time between the first and last intersections of 0 (π, π/2 or 3π/2). In this case, the orbit makes at least one full revolution (libration) around a stable point on the phase plane.

• Passage through the resonance. In this case, the value of ${\dot{\theta }}_{\mathrm{res}}$ changes its sign. When the resonant angle θ_res intersects 0 (π, π/2 or 3π/2) two times and then crosses −π/2 (π/2, 0 or π) without completing a revolution on the phase plane, this means that the orbit changes its behavior and passes through resonance. We mark such an event as “passage”.

Using our criteria for classifying orbital behavior, we not only determine whether the orbit is in resonance or in the passage, but we can also estimate the time it spends in one state or another (in terms of resonant angle changes).

It should be noted that in this paper the notation ${\dot{\theta }}_z\lessgtr {\dot{\theta }}_{\mathrm{R}}$ and ω_z/κ ≶ 1 are not the same thing. The first inequality characterizes the behavior of only circulating orbits. In the second case, ω_z and κ are calculated numerically using the procedure described at the beginning of this section for all orbits, not just circulating ones.

All Figures

Fig. 1. Vertical action and its evolution. Left and middle columns: 2D maps of the vertical action Jz plotted in different coordinates, (Lz, JR′=JR − θ(−Lz)Lz) (first column), (Lz, ωz/κ) (second column), and (2(Ω − Ωp)/κ, ωz/κ) (third column). Top row: Average value of the action in the pixel. Bottom row: Mean changes in Jz from t = 300 to t = 400 (actions and frequencies refer to the moment t = 400). The vector fields in the bottom rows demonstrate the average flow of orbits in given planes over five time units (from 400 to 405). The thin gray lines in each plot represent isodensity contours. The thick gray lines in the top left plot correspond to the lines with the same value of the adiabatic invariant Jv = JR′+Jz + Lz/2. The vertical green lines in the middle column indicate the location of the CR and ILR. Right column: Snapshots (xy and xz) of the bar orbits (abnormal orbits; see Zozulia et al. 2024) at t = 400 with different frequency ratios ωz/κ > 1 (top two plots) and ωz/κ < 1 (bottom two plots).

In the text

Fig. 2. Vertical action and its evolution for captured orbits. Left and middle columns: 2D density maps (bar and disk) in the background and 2D maps (superimposed in the foreground) of Jz for orbits that were flat (Jz < 0.05) at t = 300 and trapped into the bar over the time interval of t = 300 − 400. Different time moments are shown from top to bottom, t = 300 (top row), 400 (middle row), and 500 (bottom row). The thin gray lines in each plot represent the isodensity contours of trapped orbits. In the left column, the thick gray lines correspond to isolines of the same value of the adiabatic invariant Jv = JR′+Jz + Lz/2, and the thick cyan line corresponds to the frequency ratio isoline ωz/κ = 1. The vertical blue lines in the middle column indicate the location of the CR and ILR. Right column: Snapshots (xy and xz) of the trapped orbits at the corresponding time moments.

In the text

Fig. 3. Four examples of different orbital behaviors in our N-body simulations. For each orbit, from left to right, we show snapshots (xy and xz), the change in the resonant angle θz − θR with time, the orbit in the plane (θres,  Jz), where Jz is the midterm vertical action. The top left and right orbits remain a long time in vILR and librate around θz − θR ∼ π (banana orbit) and π/2 (anti-banana orbit), respectively, and then leave it. The bottom right orbit is trapped into vILR and also librates around θz − θR ∼ π until the end of the integration time. The bottom left orbit passes through the resonance around θz − θR = π. Orange and green in each plot correspond to the time when θz − θR increases and decreases, respectively. Blue and purple correspond to the stage when the orbit is in resonance or passes through it, respectively.

In the text

Fig. 4. Snapshots (xy and xz) of the bar (top two rows) and the structure assembled from the orbits trapped into the bar in the time interval t = 300 − 400 (bottom two rows) at t = 500. The first column shows all orbits under consideration. Structures assembled from orbits that increased and decreased their resonant angle θz − θR are displayed in the second and fifth columns, respectively. The third and forth columns show orbits in the vILR and passing through it, respectively.

In the text

Fig. 5. 1D distributions of the orbits trapped into the bar over the time interval t = 300 − 400, which were flat (Jz < 0.05) at t = 300 (the same orbits as in Fig. 2): The three left plots display the secular frequency ratio ωz/κ at t = 300, 400, and 500, respectively. The two right plots correspond to the secular vertical action Jz at t = 400 and t = 500. The gray lines show the distributions of all trapped orbits, and the other colors indicate orbits with different behaviors (the same as in Fig. 3) at t = 500.

In the text

Fig. 6. Distributions of the time that the bar (four left plots) and trapped (four right plots) orbits spend in resonance or pass through it, around the angles 0,  π, and π/2,  3π/2. Orbits with resonant angles of 0,  π are combined and are designated in the plots as θz − θR ∼ π. The same holds for the angles π/2,  3π/2. They are denoted as θz − θR ∼ π/2. The number of each type of orbit is indicated in the top right corner of each plot. The resonant and passing orbits are identified at t = 500, while the time interval t = 250 − 550 is being considered. The orange distributions correspond to orbits that enter and leave resonance or pass through it during this time period. The blue part of distributions corresponds to orbits that at the time t = 250 were already in resonance and/or were still in resonance at moment t = 550, and thus, they show a lower estimate for the time spent near the vILR.

In the text