Orbital Support and Evolution of CX/OX Structures in Boxy/Peanut Bars
Behzad Tahmasebzadeh Shashank Dattathri Monica Valluri Juntai Shen Ling Zhu Vance Wheeler Ortwin Gerhard Sandeep Kumar Kataria Leandro Beraldo e Silva Kathryne J. Daniel
























































behzad@umich.edu,
Department of Astronomy, University of Michigan, Ann Arbor, MI, 48109, USA



Department of Astronomy, Yale University, Kline Tower 266 Whitney Avenue, New Haven, CT 06511, USA


Department of Astronomy, University of Michigan, Ann Arbor, MI, 48109, USA


Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China
Key Laboratory for Particle Astrophysics and Cosmology (MOE) / Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai 200240, China


 
Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China


Department of Physics, University of Chicago, Chicago, Illinois 60637, USA



Max-Planck-Institut für Extraterrestrische Physik, Gießenbachstraße 1, 85748 Garching, Germany





Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China


Department of Astronomy & Steward Observatory, University of Arizona, Tucson, AZ 85721, USA



Department of Astronomy & Steward Observatory, University of Arizona, Tucson, AZ 85721, USA
Center for Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA













§ ABSTRACT 



Barred galaxies exhibit boxy/peanut or X-shapes (BP/X) protruding from their disks in edge-on views. Two types of BP/X morphologies exist depending on whether the X-wings meet at the center (CX) or are off-centered (OX). Orbital studies indicate that various orbital types can generate X-shaped structures. Here, we provide a classification approach that identifies the specific orbit families responsible for generating OX and CX-shaped structures. Applying this approach to three different N-body bar models, we show that both OX and CX structures are associated with the x1 orbit family, but OX-supporting orbits possess higher angular momentum (closer to x1 orbits) than orbits in CX structures. Consequently, as the bar slows down the contribution of higher angular momentum OX-supporting orbits decreases and that of lower angular momentum orbits increases resulting in an evolution of the morphology from OX to CX. If the bar does not slow down, the shape of the BP/X structure and the fractions of OX/CX supporting orbits remain substantially unchanged. Bars that do not undergo buckling but that do slow down initially show the OX structure and are dominated by high angular momentum orbits, transitioning to a CX morphology. Bars that buckle exhibit a combination of both OX and CX supporting orbits immediately after the buckling, but become more CX dominated as their pattern speed decreases. This study demonstrates that the evolution of BP/X morphology and orbit populations strongly depends on the evolution of the bar angular momentum.

















§ INTRODUCTION


In numerous observations of external disk galaxies viewed edge-on, boxy/peanut or X-shaped bulges (hereafter referred to as BP/X bulges) have been identified <cit.>. Within the Milky Way, a distinct BP/X-shaped structure was first discerned in the multi-parameter model of COBE/DIRBE images of the Galactic Bulge <cit.>. 



The BP/X bulges have been consistently observed in a variety of N-body simulations. 
Early simulations showed that BP/X bulges can form following a short-lived buckling event in a bar which results from an asymmetric bending of the bar out of the disk midplane  <cit.>.  
However, recent evidence suggests that orbital resonances may play a crucial role in forming and enhancing BP/X structures. One scenario for the formation of BP/X structures without a buckling event is 'resonant trapping,' where orbits are vertically excited by being trapped at the vertical Inner Lindblad Resonance (vILR) for significant periods <cit.>. Another mechanism is 'resonant sweeping,' in which orbits cross the vILR, become vertically heated, and remain that way after leaving the resonance <cit.>.  <cit.> reviewed all three mechanisms for the formation of BP/X structures and found evidence for the resonant trapping mechanism only in an artificially vertically symmetrized model. They concluded that the resonant sweeping mechanism is likely more applicable in real galaxies without bar buckling. This conclusion was later confirmed by <cit.> and <cit.>.

BP/X bulges are classified into two categories based on their morphology, as described by <cit.>. The two categories are the off-centered X-shape (OX), where the X-wings do not intersect at the center (which look like >-<), and the centered X-shape (CX), where the X-wings cross at the center and in the disk plane (><). Figure <ref> presents examples of an OX bulge (NGC 1381) and a CX bulge (NGC 4710), displaying their S^4G images with the corresponding Gaussian-filtered unsharp masked versions.
Orbital analysis is a crucial tool for deciphering the building blocks of barred galaxies. While the orbital composition of BP/X bulges has been extensively studied, several aspects remain contentious. Key among these is the identification of specific orbit families that independently support CX and OX structures. Another is understanding how the evolution of the BP/X structure correlates with changes in the bar's characteristics, such as its pattern speed.




In many studies of three-dimensional analytical bar models <cit.>, the periodic orbits bifurcating from the x_1 family (referred to as x_1-tree orbits) are considered the backbone of X-structures. This includes orbit families such as x_1v_1 (⌣ or ⌢) and x_1v_2 (∞), which are associated with the (Ω_z:Ω_x = 2:1) resonance. 



However, some other studies questioned these conclusions. <cit.> analyzed the orbital structure of Made-to-Measure (M2M) models for the Milky Way bar. They found 3D resonant orbits with (Ω_z:Ω_x = 2:1) or higher vertical resonances of the x_1 orbit family cannot fully explain the X-shape of the bar. In their models, the fraction of (Ω_z:Ω_x = 2:1) resonant orbits is relatively small, and they are predominantly located in the outer regions of the bar. Consequently, they proposed an alternative family of resonant boxlet orbits, termed `brezel orbits' associated with the (Ω_z:Ω_x = 5:3) resonance. These brezel orbits are posited to generate the X-shape in the inner regions of the bar. 



<cit.> studied an N-body bar model to explore the orbits responsible for building the BP/X bulge. They also found a few x_1v_1 orbits (∼ 3 %) which are in the outer half of the bar (similar to <cit.> and <cit.>). The most populated resonant boxlet family in their N-body bar is fish/pretzel orbits (∼ 6 %) that are associated with the (Ω_x:Ω_y:Ω_z = 3:-2:0) resonance. They argued that no individual orbit family is the backbone of the BP/X bulge. They found that non-resonant box orbits  (∼ 63 %), banana orbits  (∼ 3 %), fish/pretzel orbits  (∼ 6 %), and brezel orbits (∼ 1.5 %), each contributes to the formation of X-shaped structures.  This finding underscores the relatively minor role of resonant orbits in making up the BP/X structure.
<cit.> confirmed that various types of orbits can support the BP/X structures. They examined the orbital structures of two N-body bar models. They classified the bar orbits only based on the ratio of the vertical oscillation frequency to the in-plane frequency (Ω_z / Ω_x). Three groups of orbits are considered: 
1.55<Ω_z / Ω_x≤ 1.75, 1.75<Ω_z / Ω_x≤ 1.95, and  1.95<Ω_z / Ω_x≤ 2.05. Their findings demonstrated that each group is capable of forming an X-shaped structure, as illustrated in Fig. 9 in <cit.>. 
They concluded that X-shaped structures are not formed by specific orbits serving as the backbone. Instead, these structures arise from the assembling of high-density regions of different types of orbits at their highest points. 
The orbital composition of the CX and OX structures is discussed only in a few studies. <cit.> showed that periodic orbits around x_1v_2 support a CX profile (see Fig. 6 in <cit.>), while the sticky-chaotic orbits with initial conditions close to x_1v_1 and 3D quasi-periodic orbits around x_1 support an OX structure (see Fig. 7 in <cit.>). They also analyzed an N-body bar model presented in <cit.>, and found that sticky chaotic orbits support an OX profile, and periodic orbits build a CX structure. <cit.> showed the orbits with higher Ω_z / Ω_x are generating the OX shape, and those with lower Ω_z / Ω_x are contributing to the CX structure (see Fig. 9 in <cit.>). <cit.> conducted an orbital frequency analysis of live N-body models and confirmed that orbits with varying ranges of Ω_z / Ω_x contribute to the formation of an X-shaped structure. They found that during the initial stages of bar formation, prior to the buckling phase, the ratios Ω_z / Ω_x are higher. As the models evolve, the distribution of Ω_z / Ω_x shifts towards lower values. See Also <cit.>.







In this work, we revisit the orbital composition in BP/X bulges. While it is now understood that various orbit types can give rise to X-shaped structures, yet a universal classification applicable across models with diverse characteristics has not been established. Moreover, there is no clear explanation of what determines the differences between OX and CX structures, nor an understanding of how their orbital compositions differ. 


Here, we employ an automated orbit classification method based on frequency analysis <cit.> to systematically investigate the orbit families that give rise to the OX and CX-shaped structures. Additionally, we explore the influence of the bar pattern speed on the development of OX/CX structures, as well as the types of orbits that constitute bars with different pattern speeds.

The paper is organized as follows. In Section <ref>, we provide
the details of the simulations used in this study. Section <ref> presents the description of the orbital analysis and classification method. The orbital structure of the BP/X bulges is studied in Section <ref>, and the discussion follows in Section <ref>.




§ N-BODY BAR MODELS
 
We study three N-body Milky Way-like models evolved using the GALAXY code <cit.>. Model A, was kindly provided by J. Sellwood and is referenced as model C in <cit.>. This model consists of 10^6 equal-mass particles representing an exponential disk with a total mass of 4.21 × 10^10 M_⊙, and 10^6 particles as a live isotropic spherical Hernquist halo with a total mass of 4.21 × 10^11 M_⊙. This model is evolved for 6 Gyr, and the buckling instability is inhibited by imposing reflection symmetry about the mid-plane at each step of the simulation <cit.>.

We also study model run6000 and model C from <cit.>, which we refer to as model B and model C, respectively. Both models evolved under the same simulation setup. The first 9 Gyr of the evolution of Model B is used, after which it soon reaches a steady state, and model C evolved for 7.5 Gyr. See <cit.> for details. 

Both model B and model C comprised 6 × 10^6 particles representing an exponential disk with a total mass of 5.37 × 10^10 M_⊙, and 4 × 10^6 particles in a live Navarro–Frenk–White (NFW) halo <cit.> with a total mass of 6.77 × 10^11 M_⊙. Model B is characterized by a higher initial radial velocity dispersion compared to model C.

To quantify and compare the properties of the bars, we utilize standard measures, bar amplitude (A_m = 2) and buckling amplitude (A_ buck), following e.g., <cit.>. These parameters are defined by employing the m = 2 symmetry mode in an azimuthal Fourier expansion of the disk viewed face-on. These quantities are  normalized by the m = 0 mode and defined as follows:


    A_m=2=|∑_k m_k e^2 i ϕ_k/∑_k m_k|


and 


    A_buck =|∑_k z_k m_k e^2 i ϕ_k/∑_k m_k|.

where m_k, ϕ_k, and z_k are the mass, azimuth, and vertical position of the k th particle, respectively. We measure the pattern speed of the models over time using the approach developed by <cit.>, which makes it feasible to determine the pattern speed using a single snapshot.  Fig. <ref> presents the evolution of several bar properties: buckling amplitude (A_buck), bar amplitude (A_m=2),  bar length, and bar pattern speed (Ω_p), for model A (orange), model B (red), and model C (blue).

Model A does not undergo a buckling event, and the BP/X structure gradually forms, becoming more pronounced by the end of the simulation. The bar strengthens and grows until ∼ 5 Gyr, and continuously slows down by the end of the simulation. Model B experiences two buckling events, approximately at 3.8 Gyr, and 7.7 Gyr, followed by a gradual increase in bar amplitude and associated decrease in pattern speed throughout the evolution. In contrast, model C undergoes only one strong buckling event early in its evolution at around 2.8 Gyr after which the bar ceases to strengthen and grow and maintains a steady pattern speed until the end of the simulation.

Figure <ref> illustrates the face-on and edge-on projected surface densities alongside the unsharp masked images for model A (left), model B (middle), and model C (right). The top row displays snapshots taken soon after the formation of the BP/X structure, while the bottom row presents snapshots from around the end of the simulation. At first glance, it is evident that the BP/X shape has undergone significant evolution in model A and model B, but shows minimal alteration in model C.  






§ ORBIT ANALYSIS METHODS




 §.§ Orbit integration

In this section, we examine the orbital structure of all the snapshots shown in Fig. <ref>. Detailed analysis plots will be provided exclusively for model B at t=4.5 Gyr. For the remaining models, only the final results will be presented to avoid an excess of figures.

We employ AGAMA [<https://github.com/GalacticDynamics-Oxford/Agama>] <cit.> to compute the potential and orbits in our N-body models. The N-body system is frozen at the specified snapshots, and then potentials are calculated from the particle distribution using  expansion <cit.>. 


We randomly selected 15,000 initial conditions corresponding to the positions and velocities of particles from the snapshots. These orbits are then integrated over a duration of 20 Gyr (∼ 200 orbital period at the end of the bar region) in the presence of the rotating bar, with the given Ω_p for each snapshot. We store 10,000 points per orbit with equal time intervals. 
We choose a long integration time to better visualize the different bar structures discussed in the following section and to ensure that orbits in the outer regions are integrated for a sufficient period to compute frequencies accurately. However, we also test our results by performing the integration over a much shorter time of 2 Gyr (∼  20 orbital periods at the end of the bar region) and demonstrate that the integration time does not significantly affect our analysis (See Appendix <ref>).





 §.§ Orbits classification

Frequency analysis of orbits plays a crucial role in deciphering the characteristics of orbital structures within a large sample. To compute the fundamental frequencies in Cartesian and cylindrical coordinates, we utilize the NAFF (Numerical Analysis of Fundamental Frequencies) software [<https://bitbucket.org/cjantonelli/naffrepo/src/master/>]  <cit.>.
Automated classification of orbits, based on fundamental frequencies, is a well-developed technique that has been widely utilized in the literature <cit.>. As discussed in <cit.>, orbit classification based solely on frequency ratios is insufficient for a comprehensive orbital structure framework. This is because orbit families with identical frequency ratios can exhibit divergent properties, including variations in energy ranges, stability, extent, and shapes. In contrast, the auto-classification method in NAFF incorporates multiple quantities in addition to the orbital frequencies, such as the apocenter radius of orbits, the maximum values of the x and y coordinates, and the net rotation parameter. See appendix B in <cit.> for details. A new implementation of NAFF code named naif [<https://naif.readthedocs.io/en/latest/index.html>] has been publicly released in a python package with some new features by <cit.>. The orbit auto-classifier and resonance auto-finder for frequency maps are currently being developed for inclusion in the naif package (R. Ranjan et al., 2024, in preparation).
To determine whether an orbit is prograde or retrograde in the bar's reference frame, we use the following expression, which we call the net rotation parameter:


    L_zn =  Σ_i=1^i=Nsign(L_z)_i/N,
	
where N is the number of time steps. For a positive and negative value of L_z, sign(L_z) is +1 and -1, respectively. A short-axis orbit with L_zn∼ 1.0 is prograde with the highest angular momentum, while L_zn∼ -1.0 indicates a retrograde orbit around the z-axis in the bar rotating frame. An orbit with no net rotation around the z-axis has L_zn∼ 0.0. A long-axis tube orbit has L_xn∼∓ 1.0 since it has net rotation around the x-axis.
In this approach, disk orbits are determined by those that have apocenter radii greater than the bar radius.

Then the orbits within the bar are classified into:   boxo (ordinary box), boxp (periodic or resonant box, i.e. boxlet), ztub (z-axis tube), ztup (periodic or resonant z-tube), xtub (x-axis tube), xtup (periodic or resonant x-tube), x2++ (similar to z-axis tube but elongated along the y-axis of the bar and prograde), x4++ (similar to z-axis tube but elongated along the y-axis and retrograde), bo32 (box with 3:2 fish resonance), bobr (box brezel with 5:3 resonance), x1++ (resonant or near resonant x_1 boxlet), x1cl (x_1 with 3:1 resonance), x122 (x_1 with 2:2 resonance), x132 (x_1 with 3:2 resonance), x1bn (x_1 banana). We lumped all of the periodic x_1 resonant orbits together since there are very few.
<cit.> demonstrated that this automated orbit classification is generally consistent with visual classification by examining 20,000 individual orbits. The differences in classification accuracy range from approximately ∼ 1% to ∼ 4 % depending on the types of orbits being classified.
NAFF employs the frequency drift parameter (diffusion rate) to measure the fraction of chaotic orbits. The frequency drift parameter log_10(Δ f) is computed by the change of the frequency in two equal orbital time segments <cit.>.
We use log_10 (Δ f) >-1.2, which is a good empirical criterion for determining chaotic orbits as tested in <cit.> and <cit.>. All orbits classified into the aforementioned families are additionally tagged with regl (regular) or chao (chaotic).
We found that some orbits initially classified as periodic boxes exhibit chaotic behavior after a short time. Consequently, those orbits demonstrating chaotic characteristics were reclassified from the resonant orbit classes  boxp, bo32, and bobr to non-periodic box class boxo. Furthermore, we noted that some orbits initially classified as periodic z-axis tubes (ztup), characterized by L_zn < 0.0 and -0.85 < L_xn < 0.85, should be reclassified into the x4++ group. Similarly, those with L_zn < 0.0 and L_xn∼± 1 are more accurately categorized as xtub orbits. Although these adjustments result in minor changes to orbital classes, they are crucial for our subsequent analysis of the morphology of different orbit classes.




 §.§ Orbital decomposition of CX and OX structure

We further categorize the orbit types mentioned above by grouping together orbital families that exhibit similar morphologies.

(1) OX: These are prograde short-axis tube orbits that are elongated along the bar, including 
 x1++, x1cl, x122, x132, x1bn, ztub, and ztup. They make up an OX structure.
(2) CX: These are box orbits that are elongated along the bar, contributing to the formation of the CX structure. To select these groups of orbits, we begin by aggregating both periodic and non-periodic box orbits, including  boxo,  boxp, bo32, and bobr. Box orbits exhibit a diverse range of morphologies. We have determined that the parameter L_zn serves as an effective criterion for distinguishing between box orbits that support a CX shape and those that exhibit a more rounded boxy morphology. Fig. <ref> presents the surface densities of box orbits across various L_zn bins. The columns from left to right display the surface density in the x-z planes for different L_zn bins, which is plotted for model B at t=4.5 Gyr.

Box orbits possessing higher values of L_zn prominently exhibit an OX structure, which transitions to a CX shape as the value decreases. With further reduction in L_zn, the X wings gradually disappear. By examining the structures within various bins, we have identified a specific value of L_zn to serve as a threshold for visually distinguishing between OX and CX structures across all box orbits. For Model B, as depicted in Fig. <ref>, box orbits with L_zn > 0.3 are categorized as the OX group as they clearly exhibit a bridge connecting the two X wings. Those with 0.3 > L_zn > -0.2 are classified as the CX group, while the remaining orbits, which form a round shape, are included in the general box orbit group. These boundaries are determined visually and may vary from model to model. It should be noted that although altering the threshold value of L_zn can result in minor variations in the contribution of OX/CX orbits, these changes are typically limited to a few percent and are not significant.

(3) box: These orbits are dominated by retrograde motion and are not included in CX or OX groups.
(4) x_4: These orbits are retrograde short-axis tube orbits elongated perpendicular to the bar, referred to as x4++ in the NAFF classification. It is worth noting that x2++ orbits are their prograde counterparts, also short-axis tubes with a similar perpendicular orientation. However, our models do not contain any x2 orbits.
(5) LAT:  These are long-axis tube orbits - orbits that have net angular momentum about the long x-axis of the bar model.
(6) disk: Orbits with apocenter radii exceeding half the length of the bar are classified as disk orbits.







§ RESULTS


Fig. <ref> presents the extracted surface densities in the x - y and x - z planes for six orbit categories. The columns, from left to right, display surface densities constructed from orbits in the OX, CX, box, x_4, LAT, and disk categories, respectively. In all models, the OX and CX orbits contribute to structures elongated along the bar, forming OX- and CX-shaped structures respectively (shown in the 2nd and 3rd columns). The box, x_4, and LAT orbits do not contribute to the BP/X shape, instead supporting rounder or slightly boxy structures. In the following, we explore the orbital origins and evolution of the OX and CX structures in detail.

Furthermore, we demonstrate that orbits supporting the OX/CX vertical structures contribute to outer/central bar-elongated structures when the disk is viewed face-on. Our results confirm the lack of a relationship between X-shaped structures and the inner rounded parts of the face-on bar (e.g., barlenses). <cit.> illustrated the importance of orbits with loops at their ends (in the x-y plane) to the bar shoulders in face-on density profiles. Additionally, they showed that bar thickening and vertical resonances can dilute the shoulders in the bar major-axis density profiles. We confirm that it is the OX rather than the CX orbits that contribute to the bar shoulders. This conclusion aligns with the findings from <cit.>, since the OX orbits have higher angular momentum and are therefore closer to their parent x_1 orbits with loops at their ends, while the CX orbits have lower angular momentum and are vertically thicker.







 §.§ Orbital origin of CX and OX structure


The fraction of periodic orbits, including closed x_1 orbits, is less than <3% in all our models, and the majority of orbits that make up the BP/X-shaped structure are non-periodic. This is consistent with previous studies, which show that periodic orbits (such as banana, pretzel, etc.) constitute only a small fraction of orbits in N-body bars. Fig. <ref> depicts typical x_1 orbits (left panel), non-periodic orbits generating the OX (middle panel), and CX structures (right panel) in the x-y and x-z planes.



<cit.> argued that assembling non-periodic orbits into a pattern resembling a `bowtie' can create X-structures. These orbits occupy large areas on the x-z plane and can be linked to quasi-periodic orbits surrounding the plane of the x_1 orbital family. For examples of such orbits, see Fig. 15 in <cit.>, Fig. 2 in <cit.>, and Fig. 9 in <cit.>. 

Based on the computed frequencies from NAFF, we confirm that the orbits contributing to the OX and CX structures can exhibit a wide range of Ω_x/Ω_z values, extending beyond those typical of resonant orbits, such as banana or pretzel orbits. This finding is in agreement with previous studies.

As we demonstrated in Fig. <ref>, the typical non-periodic orbits within the OX and CX groups exhibit similar bowtie-like shapes. However, OX orbits tend to be vertically thinner bowties, whereas CX orbits are characterized by centrally thicker bowties. <cit.> suggested that further investigations using the Poincaré surface of sections are required to elucidate the origins of such orbits.


To elucidate the origin and characteristics of orbits within each OX and CX structure, in the following, we examine their phase-space and investigate their origin using Poincaré surface of sections plots.



  §.§.§ Poincaré surface of section


A useful technique for exploring the distribution of orbits in phase space is to visualize their Poincaré surfaces of section. Typically, the surfaces-of-section (SoS) is plotted for a set of orbits that share the same Jacobi integral value E_J. For orbits in two-dimensional bars, the SoS is plotted by mapping the velocity component V_y against y at each time step when an orbit intersects the x-axis with a negative V_x value <cit.>.
 


Plotting the SoS of orbits in an N-body simulation is complicated and presents two challenges: firstly, since all the orbits possess three-dimensional structures, they cannot be adequately represented on a two-dimensional SoS. Secondly, the distribution of E_J values is broad and continuous (rather than discrete), leading to fuzzy curves on the SoS. However, <cit.> demonstrated that a two-dimensional SoS can still offer valuable insights into the underlying three-dimensional orbital structures. 


We explore the SoS of 100 orbits in our models. Fig. <ref> illustrates V_y versus y when orbits cross the y-z plane with a positive V_x. This is plotted for orbits within the energy range of -1.0 < E_J < -0.98 in Model B at t=4.5 Gyr. The spread in E_J values results in fuzzy curves in the SoSs, particularly when plotted for a large number of orbits, as highlighted in Fig. 10 of <cit.>. In our analysis of frozen potentials, most orbits traverse a broader area rather than being confined to narrow curves in the SoSs.

 
The left panel of Fig. <ref> displays the SoS for x_1 orbits (in red), OX orbits (in cyan), and CX orbits (in blue). To avoid overcrowding, the SoS for box orbits (in green) is presented separately in the right panel. The x_1 orbits are the ovals at the center of the bull's eye of the SoS. It is similar to what is shown in  <cit.> (see Fig. 9) and <cit.> (see Fig. 10). OX orbits are primarily clustered around V_y = 0 at positive y values, confirming they are parented by x_1 orbits and exhibit prograde motion around the z axis. Orbits in the red and cyan regions contribute to the formation of the OX structure. CX orbits, while also lying on oval curves surrounding the x_1 orbits, cover a wider region, indicating their inclusion in the same sequence. CX orbits are predominantly found at positive y values. However, as CX orbits identified with 0.3>L_zn>-0.2, SoS of those orbits with negative L_z encompass negative y values. Box orbits lack a dominant direction of rotation, implying that as orbits evolve into a box-like shape, their angular momentum approaches zero on average.



In summary, we confirm that orbits closer to the x_1 family exhibit an OX structure. As they move away from the x_1 orbits and become thicker, they transition to a CX shape. Eventually, they evolve into completely boxy or rounded shapes as they lose angular momentum.













 §.§ Evolution of CX and OX orbital structures

Fig. <ref> presents contour densities of OX and CX orbital structures in model A (left), model B (middle), and model C (right). The first row displays these orbital structures shortly after the formation of the BP/X structure, while the second row illustrates these structures at the end of the simulations. 
The proportions of OX and CX orbits in the models at various times are indicated in Fig. <ref>. In Model A, the BP/X bulge does not undergo buckling, and its pattern speed decreases significantly over time. Initially, the dominance of the OX orbits is evident at 25% of all orbits, while the CX orbits represent only 8%. By the end of the simulation, the OX orbits have almost completely disappeared, dropping to 3%. In contrast, the BP/X bulge becomes predominantly characterized by the CX orbits, which increase to 36%.

Model B and Model C both undergo buckling, leading to the formation of the BP/X bulge where both OX and CX orbits coexist. In Model B, as the bar experiences a slowdown by the end of the simulation, the proportion of OX orbits decreases while that of CX orbits increases. The first row in the middle panel shows the OX/XC structure of model B after the first buckling event, and the bottom row displays the OX/XC structure after the second buckling. In contrast, in Model C, where the bar pattern speed remains constant, the proportions of CX and OX orbits stay the same by the end of the simulation. This suggests that as long as the bar's pattern speed does not decrease, the proportions of CX and OX orbits, and therefore the morphology of the BP/X structure remain unchanged. 

We defer to a future study the exploration of whether there is a relationship between the intensity of buckling and the proportions of CX and OX orbits after the formation of the BP/X structure.



 §.§ Photometric Parametrization of BP/X Bulges


We use the method presented by <cit.> and IMFIT software <cit.> to parameterize the shape of the BP/X bulges. In this parametrization, the bar is modeled using a sech^2 profile, applied to a dimensionless, scaled radius given by

    ρ_pea  = ρ_0 sech^2(-R_s),

where

    R_s=([(x/X_bar)^c_⊥+(y/Y_bar)^c_⊥]^c_ / c_⊥+(z/Z_bar)^c_)^1 / c_ .

The coordinate system is centered at the galaxy's nucleus, with X_bar,
Z_bar and Y_bar denoting the semi-major, intermediate, and semi-minor axis lengths of the bar, respectively. c_ and c_⊥ control the diskiness/boxiness of the bar introduced by <cit.>. Here, we adopt 
c_,c_⊥∈ [1.5,5], based on the results of <cit.>. The BP/X feature's morphology is defined by a scale height perpendicular to the disk plane, which varies across the x - y  plane and is modeled by a double Gaussian distribution centered on the galactic center:

    Z_bar (x, y)=     A_pea exp(-(x-R_pea )^2/2 σ_pea ^2-y^2/2 σ_pea ^2) + 
        A_pea exp(-(x+R_pea )^2/2 σ_pea ^2-y^2/2 σ_pea ^2)+z_0,

where R_pea and σ_pea denote the distance of the peanut's center from the galactic center and the width of each peanut, respectively. A_pea quantifies the vertical extent of the peanut feature above the ellipsoidal bar's scale height z_0. 

As demonstrated by <cit.>, these three parameters provide substantial flexibility in capturing the diverse morphologies exhibited by the BP/X feature. Although these values do not directly quantify the OX/CX structures, they are useful for demonstrating the relationship between the evolution of BP/X morphology and the bar pattern speed independent of orbital analysis. The quantification of OX/CX structures through photometric parameterization will be addressed in future studies. 

Figure <ref> displays the variation of BP/X shape parameters R_pea, ρ_pea, and σ_pea versus the bar pattern speed (left column) and versus the dimensionless bar rotation parameter ℛ≡ R_ cor/R_ bar (right column) for model A (orange), model B (red), and model C (blue). The x-axis of the left column is the bar pattern speed in the reverse direction. Models A and B exhibit significant changes in BP/X morphology over time as the bar pattern speed decreases. In contrast, the BP/X parameters in model C are relatively constant or change very little, which is attributed to its nearly constant bar pattern speed.
In models A and B, the relation between the BP/X parameters and the bar pattern speed is monotonic (although not linear). This holds true even after the second buckling in model B (represented by the red point number 7).
As the bar pattern speed decreases, all three parameters R_pea, ρ_pea, and σ_pea increase. However, we do not find a monotonic relationship between the BP/X parameters and  ℛ. 
This indicates that the evolution of BP/X morphology is correlated with the evolution of the bar pattern speed rather than ℛ.










 §.§ OX/CX bulges and the bar pattern speed in Auriga simulations

The Auriga simulations <cit.> are a suite of 30 magneto-hydrodynamical cosmological zoom-in simulations. The Auriga barred galaxies exhibit fast bars primarily because they are baryon-dominated disks, experiencing less dynamical friction, which prevents the bars from slowing down <cit.>.


<cit.> analyzed the structural and photometric properties of 21 barred galaxies from the Auriga simulations. Using edge-on unsharp-masked images at z=0, they identified BP/X structures in the inner parts of 6 out of 21 galaxies. They visually determined whether these BP/X structures are OX or CX, and found that only one of these bulges has CX morphology, while the other five have OX structures. See Table 3 in <cit.>.


<cit.> presented the evolution of the pattern speed for five Auriga galaxies with high cadence outputs as a function of lookback time. See Fig B.1 in <cit.>. Three of these galaxies, A17, A18, and A26, are identified as having BP/X bulges by <cit.>. A18 galaxy, the only one with a CX structure, has a relatively lower pattern speed of around  27 km s^-1 kpc^-1. A17 and A26, which have OX structures, exhibit significantly higher pattern speeds around 45 and 37 km s^-1 kpc^-1, respectively. Our results explain the findings in the Auriga simulations: bars with higher pattern speeds tend to exhibit BP/X bulges with more dominant OX structures. 



§ SUMMARY



We present a methodology that enables us to decompose the orbital structures supporting OX and CX shapes in BP/X bulges. Our method relies on auto-classification with NAFF software, focusing on the morphological and kinematic characteristics of orbits.

To understand the origin and evolution of OX and CX structures in BP/X bulges, we applied our method to classify orbits using a frozen potential from three N-body bar models, each with distinct features, at two different times, beginning shortly after the formation of the BP/X bulge. The main results are as follows:
(1) We demonstrate that both OX and CX structures are composed of non-periodic orbits associated with the x_1 orbit family. OX orbits, being closer to x_1 orbits, possess higher angular momentum compared to those in CX structures. 
(2) We found that orbits supporting OX/CX structures in the edge-on view also form a similar OX/CX-like shape in the face-on view. 
(3) We showed that as the bar pattern speed decreases, the contribution from OX orbits decreases while that from CX orbits increases. If the bar remains at a constant speed, the shape of the BP/X structure stays unchanged until the end of the simulation. This indicates a strong dependence of the BP/X bulge shape on the bar pattern speed.
(4) We show that Model A (the bar that does not buckle) is initially dominated by OX orbits. These orbits surrounding the x_1 orbits begin to thicken. As the bar loses angular momentum, the OX structure transitions to a CX structure by the end of the simulation. In contrast, our bar models that experience buckling events exhibit a combination of both OX and CX structures immediately after the buckling and evolve toward CX dominance if the pattern speed decreases.
(5) Using photometric parameterization of the BP/X structure, we find that the evolution of BP/X morphology is correlated with the bar pattern speed rather than the bar rotation parameter.








§ ACKNOWLEDGEMENTS

We thank J. A. Sellwood for providing the data for Model A. BT, SD, MV, VW, LBeS gratefully acknowledge funding from the National Science Foundation (grant NSF-AST-2009122 to MV). LBeS is supported by the Heising Simons Foundation through the Barbara Pichardo Future Faculty Fellowship from grant # 2022-3927.




Software: AGAMA <cit.>, NAFF <cit.>, Jupyter Notebook <cit.>,  matplotlib <cit.>, numpy <cit.>,  scipy <cit.>.  




§ DATA AVAILABILITY

The simulation data for model B and model C are available at doi:[10.5281/zenodo.11237062]https://doi.org/10.5281/zenodo.11237062 Readers wishing to obtain the simulation data for Model A should contact J. A. Sellwood. Any additional data will be shared on reasonable request to the corresponding author.




aasjournal











 
 



§ THE ORBIT INTEGRATION TIME
 

To investigate the effect of integration time on our orbit classification, we repeated our analysis using an integration time of 2 Gyr, which corresponds to approximately 20 orbital periods at the end of the bar region. We employed an orbit integration with the required accuracy of 10^-8. One primary concern is that a very long integration time with low integration accuracy may result in an overproduction of chaotic orbits. We examined the fraction of chaotic orbits under integration times of 20 Gyr and 2 Gyr, using the frequency drift parameter log _10(Δ f) as described by <cit.>. We adopted log _10(Δ f) > -1.2 as the criterion for defining chaotic orbits, consistent with other studies <cit.>.

Figure <ref> presents the distributions of the frequency drift parameter for 15,000 orbits integrated over durations of 20 Gyr (red) and 2 Gyr (blue). We found that the fraction of chaotic orbits is similar in both cases, with the orbits integrated for 2 Gyr showing approximately 2% more chaotic orbits than those integrated for 20 Gyr. This difference could also be due to the randomness in selecting initial conditions, which vary for each orbit sample. This result suggests that the accuracy we are using for orbit integration is sufficient. The overall fraction of chaotic orbits in our orbit sample, as used in this paper, is 10-13% across different models.

Another concern when adopting a longer integration time is achieving more accurate orbital frequencies. However, our orbit classification is not solely based on orbit frequencies; we also consider parameters such as angular momentum, maximum z, and the ratio of maximum y to maximum x. Figure <ref> shows a Cartesian frequency map, colored by number density, for orbits integrated for 20 Gyr (left panel) and 2 Gyr (right panel). The frequency maps are very similar, and we did not find any systematic changes between them. Figure <ref> illustrates the effect of integration time on our classification, indicating only minor changes in each orbital group (1-3%), which could also be due to randomness in selecting initial conditions. Overall, these tests indicate that our analysis outcomes throughout this paper are not significantly dependent on the integration time.












