Evolution of Stellar Bars in Spinning Dark Matter Halos and Stellar Bulges

We use high-resolution numerical simulations to follow the barred disk evolution in a suite of models with progressively more massive stellar bulges, with bulge-to-total (disk$+$bulge) mass ratios of $B/T\sim 0-0.25$, embedded in dark matter (DM) halos with the spin $\lambda \sim 0 - 0.09$. We focus on models with a sequence of initial rotational support for bulges, and analyze their spinup and spindown. We find that (1) the presence of a bulge affects the evolution of stellar bars, i.e., the timescale of bar instability, bar pattern speed and its decay, and the vertical buckling instability. The bar strength is nearly independent of $B/T$ in halos with spin $\lambda = 0$, and is suppressed by a factor $\sim 2$ for halos with $\lambda = 0.09$; (2) The main effect of the bulge is the destruction of the harmonic core which affects the buckling; (3) The bulge plays a minor role in the exchange of angular momentum between the barred disk and the DM halo, during its spinup and spindown; (4) Most interestingly, the buckling process triggers different response above/below the disk midplane, which anti-correlates with the bulge mass; (5) In spinning halos, the buckling process has a prolonged amplitude tail, extending by few Gyr, as verified by measuring distortions in the Laplace plane; (6) Furthermore, as verified by orbital spectral analysis, the bulge gains its spin from the bar mainly via the inner Lindblad resonance, while losing it via a number of resonances lying between the outer and inner Lindblad resonance.


INTRODUCTION
Galaxies are commonly thought to reside in dark matter (DM) halos whose properties follow from the larger scale filamentary structure in the universe -the cosmic web.These halos can be characterized by a number of parameters, such as mass, spin, baryon fraction, environment, and more.While halos are far from being supported by rotation, this is a net characterization, and various components within halos can still harbor a substantial amount of angular momentum which can potentially be transferred most easily towards the halo.This angular momentum appears to have significant dynamical and secular effects on the evolution of embedded stellar disks, especially when the latter ones are barred (Collier et al. 2018(Collier et al. , 2019a,b;,b;Li et al. 2023b).
Galaxies possess an additional spheroidal component to the DM component -the classical stellar bulges.While formation and evolution of these bulges have been extensively investigated over the last 50 years, their interaction with barred disks still has a number of unanswered questions.While work on the spinup of these bulges has been performed (e.g., Sellwood 1980;Saha et al. 2012Saha et al. , 2016)), the parameter space involving the range of the bulge mass fraction of the total parent galactic disk+bulge, /, deserves more careful study.
★ E-mail: xingchen.li@uky.edu† E-mail: isaac.shlosman@uky.eduMoreover, most the classical bulges studied so far have been nonrotating initially, and it is not clear how the angular momentum evolution of spinning bulges behaves.The focus of our analysis, therefore, is (1) to investigate the angular momentum flow in the barred galaxies with a range of / and the initial spin of these bulges, and (2) to investigate evolution of stellar bars in spinning bulges and DM halos under various initial conditions.
A large sample of local galaxies using the optical -band imaging, including the early and late-type disks have shown that the bar fraction is anticorrelating with the bulge-to-total mass ratio, and rising from 40% to 70% (Barazza et al. 2008).On the other hand, Hoyle et al. (2011) found that bulge-dominated galaxies host longer stellar bars compared to bulgeless galaxies.This point has been also addressed by Sellwood (1980) and Athanassoula & Martinet (1980).
Numerical simulations have demonstrated that bulges can exert a substantial influence both on the bar instability and on the vertical buckling instability in stellar bars (e.g., Sellwood & Gerhard 2020).Kataria & Das (2018) found that massive bulges prohibit the bar formation when the ratio of the radial force due to bulge to that of the galaxy at the disk scalelength radius is larger than 0.35.Large velocity dispersion weakens the bar instability in the disk with such massive bulges (e.g., Athanassoula & Sellwood 1983).Furthermore, Kataria & Das (2019) found that the slowdown rate of the stellar bar is stronger in the presence of a more massive bulge, possibly due to the massive bulge gaining more angular momentum from the barred disk.
It was also demonstrated that the boxy/peanut bulges can be formed by a variety of processes and may even bypass the break of the vertical symmetry in the presence of the vertical inner Lindblad resonance (vILR) (Friedli & Pfenniger 1990).Sellwood & Gerhard (2020) analyzed three such mechanisms, namely, the usual vertical buckling of a bar, the heating of the bar by the vILR in the presence of a nuclear cluster and vertical heating when the vertical symmetry is enforced.Note that the nuclear cluster had a radius of ∼ 200 pc, which was also the gravitational softening in the simulations.Sellwood (1980) was the first to simulate the angular momentum transfer between a stellar bar and a live DM halo in numerical simulations, while Weinberg (1985) studied the interaction between a rigid disk bar with a spheroidal component.Saha et al. (2012) analyzed the interaction between a low bulge-to-disk mass ratio, / ∼ 0.067, with a bar, with initially nonrotating bulge, and the angular momentum channeled to the bulge via resonances, especially via the 1:1 resonance (Lagrange points) and the inner Linblad resonance (ILR) 2:1.The bulge was found to spin up and to become an anisotropic triaxial object embedded in a boxy bulge, developing a near cylindrical rotation.Collier et al. (2019a,b) has analyzed the retrograde orbits in spheroidal components and found that such orbits can absorb the retrograde angular momentum.
The bulge has been modeled with / ∼ 6.7% for 3 Gyr, with the initial spin measured in rotational-to-dispersion velocities, rather than in the angular momentum units (Saha & Gerhard 2013).It has affirmed that initially spinless bulge gains more angular momentum.The goal of this work was to study the boxy/peanut-shaped bulge evolution when its classical progenitor is spinning.The bar amplitude and size have been found to be reduced by ∼ 10 − 40%.Finally, Saha et al. (2016) studied a spinup of initially spinless bulges, using bulges with / ∼ 11 − 43%.
Recent modeling of isolated disk galaxies have shown that parent DM halos have a two-fold effect on the stellar bar development, which depends on the halo spin λ and its DM density (Li et al. 2023b).The halo spin has been confirmed to speed up the bar instability as noted before (Saha & Gerhard 2013;Long et al. 2014;Collier et al. 2018).The decrease in the DM density has led to the appearance of a plateau in the bar strength during which the bar pattern speed remains constant.The buckling instability has been found to be separated from the maximum of the bar strength by this plateau whose length can be few Gyrs.
In this work, we study the effect of initial conditions in the bulgedisk-halo system on the stellar bar evolution.We use a sequence of bulge masses and spins in tandem with varying the parent DM halo spin.We aim at obtaining the general picture for bar evolution in the presence of a range of bulge and DM halo properties.
This paper is structured as following.Section 2 describes the numerical methods used, our results are presented in section 3, which is followed by discussion and conclusions.

NUMERICAL METHODS
We use the -body part of the mesh-free hydrodynamics code gizmo (Hopkins 2015), an extension of the gadget-2 code (Springel 2005).The models of stellar disk galaxies with central spherical bulges embedded in spherical DM halos have been constructed with different bulge-to-disk mass ratios and different spin in two spherical components, the stellar bulges and DM halos.
The units of mass, length, and velocity have been chosen as 10 10 M ⊙ , 1 kpc, and 1 km s −1 , respectively.As a result, the unit of time is approximately 1 Gyr.The number of DM halo particles has been taken as  DM = 7.2 × 10 6 , and the total stellar particles (disk plus bulge) has been fixed to  S = 8 × 10 5 .The mass-per-particle of disc and bulge are the same, and similar to that of the halo.The gravitational softening length for both DM and stellar particles is  DM =  S = 25 pc.The opening angle  of the tree code has been reduced from 0.7 used in cosmological simulations to 0.4 for a better quality of force calculation.All models have been run for 10 Gyr, with the angular momentum conservation within 0.5% and energy conservation within 0.2%.

Initial Conditions
The initial DM component consists of a spherical halo with the modified Navarro, Frenk & White (1996, hereafter NFW) density profile as a function of a spherical radius , where  s is a normalization parameter,  c = 1.4 kpc is the size of the flat density core, and  s = 10 kpc is a characteristic radius.We use a Gaussian cut-off radius  t = 180 kpc to obtain a finite mass.
The initial conditions of an exponential stellar disk with a density profile are given by where  is the cylindrical radius,  d is the mass of the disk ,  0 = 2.85 kpc is the disk radial scalelength, and  0 = 0.6 kpc is the disk scaleheight.The disk is truncated at 6 0 ∼ 17 kpc, i.e., at 98% of its mass.Velocities of disc particles have been assigned using the epicycle approximation and asymmetric drift correction.The disc dispersion velocities are where  ,0 = 120 km/s, and  ,0 is determined by setting the minimal Toomre parameter  (e.g., Binney & Tremaine 2008) , to  = 1.5, at ∼ 2 0 in all models.
To obtain the initial mass and velocity distributions in spheroidal components, we adopt the iteration method introduced by Rodionov & Sotnikova (2006), see also Rodionov et al. (2009); Collier et al.
Table 1.From left to right: DM halo masses, total disk + bulge mass (), bulge mass (), bulge-to-disk mass ratio (/), and bulge-to-total mass ratio (/).Each of these four bulge models has been run with three different spins, namely, S50, S75, and S95, as explained in the text.For each of the bulge models, two halo spins have been run, namely, λ = 0 (P00 model) and λ = 0.09 (P90 model).Li et al. (2023a).At each iteration, we freeze the disk and the bulge, and release the DM particles from the initial density distribution for 0.3 Gyr.Next, we assign to each of the unevolved particles, the nearest evolved particle velocity, and iterate.A total of 70 iterations have been executed until the virial ratio and the velocity distribution of DM particles converge.

Models
For the stellar bulge, we use the Plummer model with a scale length 1 kpc.We utilize the aforementioned iteration method to relax the bulge while the halo and the disc are frozen.About 50 iterations, each lasting 0.03 Gyr, are implemented until the virial ratio and the velocity distribution of the bulge particles converge.
We use four bulge mass models which differ only by the bulge mass fraction of the total mass, i.e., disk + bulge mass being kept fixed, /.The bulge models have / = 0, 0.1, 0.17 and 0.25 bulge mass fractions and are given in Table 1, where we follow the following notation: the model name, e.g., B10, represents the mass fraction of the bulge, i.e., bulge-to-total mass ratio.
For the spherical components, the DM halo and the stellar bulge, the angular momenta are zero after the iterations.To spin up the halo and the bulge, we randomly sample the retrograde particles (with respect to the disc rotation axis) in each component, with a constant fraction at all radii, and reverse the azimuthal velocities,   , of these particles to achieve the desired spin.In comparison with the bulgeless models in Li et al. (2023a,b), we obtain spinning halos with λ = 0.09 (P90 model) and additional three sets of initially nonspinning bulge and two initially spinning bulges.For the spinning bulges we have reversed 50% and 90% of originally retrograde particles.
The initial conditions for the density profiles in DM halo, stellar disk and stellar bulges are given in Figure 1, and the circular velocities in each component are shown in Figure 2. From left to right the Figures correspond to B10, B17 and B25 bulge models.

Maximal angular momentum in bulge-disk-halo systems
The angular momentum redistribution in the bulge-disk-halo systems is one of the dominant factors in their evolution.Following the work of Sundman (1913) (see Appendix), the total angular momentum of any steady collection of particles is bound by the product of the particles moment of inertia and kinetic energy.Initially the mass distribution of each of the bulge, disk and halo components are the same for all models.Since the mass distribution fixes the potential energy, the virial theorem fixes the total kinetic energy, but not the components respective kinetic energies, so Sundman's inequality may provide different bounds for each models.We elaborate on this issue in the Appendix, where we use tighter inequalities derived by Pfenniger (2019) when the spin vector orientation is known.These inequalities are tighter than Sundman's because they do no use the  and   quantities.
In Table 2, we calculate the tightest bounds on the maximal angular momentum of the above three components, i.e., the DM halo, disk and bulge, for our models using Eq.(A7).Table 3 presents the initial angular momentum in our models.
Since the average velocities do not change much in galaxy components, we see that the maximum angular momentum in the bulge component, which is small both in mass and size, must be small in regard of the disk, and the bulge cannot absorb much angular momentum.For the disk and halo there is a competition between the masses, radii, and velocities, but in view of the larger mass and size of the halo, it is the largest potential angular momentum reservoir by far.
Plugging in numbers, the bounds on the maximum angular momentum of each component of our models are, in 10 13 M ⊙ kpc km s −1 units: The disk is initially near its maximum angular momentum content, so we can expect that most of the angular momentum will migrate from the disk to the halo, because the bulge is unable to absorb much of it.This does not mean that the bulge is irrelevant.Since its dynamical time-scale is the shortest of the three components, any angular momentum exchange, even if small in magnitude, can trigger first instabilities like a stellar bar, or a bending, which may propagate into the larger components, determining the later evolution of the system. .Initial circular velocities for the DM halo (dashed), stellar disk (dotted), stellar bulge (dashed-dot), and the total (solid) in models with three different bulge masses.From left to right: B10, B17 and B25.Note that the disk + bulge masses have been kept the same in all models.

The stellar bar characteristics
In all models, stellar bars develop from the initially axisymmetric mass distribution.We quantify the bar strength by the Fourier components of the surface density.For  mode, the Fourier amplitude is where and Σ(, ) is the surface stellar density.To quantify the bar strength, we use the normalized  2 amplitude which is defined as We choose the upper limit of integration,  max , as the radius which contains 98% of the disc mass at a given time.
The bar length has been measured using orbital analysis, which is the most reliable method (Heller & Shlosman 1996;Berentzen et al. 1998;Martinez-Valpuesta et al. 2006;Collier et al. 2018).We start by computing the fundamental orbit family, the x 1 family in the notation of Contopoulos & Papayannopoulos (1980), which constitutes the backbone of a stellar bar.Next, we determine the extent of this family to the highest value of Jacobi energy.This family is terminated inside the CR.
We measure the vertical buckling strength,  1 , i.e. the vertical asymmetry, by calculating the the  = 1 Fourier amplitude in the -plane, where the major axis of the bar is aligned with the -axis, and the rotation axis is along -axis, The integral is over the region || < 12 kpc, || < 3 kpc, || < 5 kpc.
The phase of the bar,  bar , is obtained from Generally,  bar displays small variations, because of the noise in  2 and  2 .We take an average  bar in a range of  which defines the bar size.

Measuring buckling amplitude using an alternative method
The above method of measuring the vertical buckling of stellar bars has provided inconclusive results for some models, due to noise.To remove this ambiguity, we have used a modified vertical asymmetry parameter  asym introduced by Smirnov & Sotnikova (2018) to measure the buckling strength.The  asym is defined as where  2 ( > 0) and  2 ( < 0) are the disk Fourier amplitudes,  2 , from the numerator in Equation 7 calculated above and below the disk mid-plane, and  0 is calculated from the whole disk.The origin ( = 0,  = 0,  = 0) is set to the center-of-mass of the disk and the  = 0 plane is perpendicular to the vector of total angular momentum of the disk component within its 90%-mass-radius.
Our definition of  asym differs from Smirnov & Sotnikova (2018) which normalized it by  2 .This normalization measures  asym in terms of the bar strength.We find that it is advantageous to measure it in units of  0 , which does not vary in the simulations.

Measuring the bulge shape
To measure the bulge shape evolution we apply the method of Heller et al. (2007) by computing the eigenvalues of the moment of inertia tensor for the mass within a specified radius.This determines the axes  >  >  of a uniform spheroid with the same eigenvalues.The ratios of these axes are used to characterize the shape of the bulge, as follows, where Δ =  2 −  1 +  3 for the eigenvalues  3 >  2 >  1 .

RESULTS
We start by describing the evolution of stellar bars in all the models with the bulge and DM halos with λ = 0 and 0.09.As a next step we analyze the bulge evolution for various bulge masses and spin.Finally, we address the interactions between the stellar bars and bulges.

Evolution of disk bars in nonspinning halos
The evolution of stellar disks with nonspinning bulge models immersed in the nonspinning halos is shown in Figure 3, at  ∼ 6.25 Gyr, the time when the developed stellar bars have similar strength.The disk bars appear shorter in the more massive bulge models.Also, in all models with bulges in nonspinning DM halos, the surface density in the central few kpc exhibits a dumbbell distribution.The dumbbell does not develop in models with spinning halos of λ = 0.09.The evolution of stellar bar strength in all our models is displayed in Figure 4. We start by analyzing the bars in nonrotating halos (top row in this Figure corresponding to λ = 0 halo).The maximal  2 amplitude before buckling is ∼ 0.3 and similar in all the models, showing only a slight decrease, less than 10%, with increasing bulge mass.We also notice that the bar growth is slower for more massive and faster spinning bulges.This is to be compared with  2 ∼ 0.4 for the bulgeless model.The end product after 10 Gyr of evolution resulted in strong bars with  2 ∼ 0.35 − 0.45.
The growth rate of the stellar bar displays a two-fold effect: it is slower for more massive bulges, and slower for faster spinning bulges as well.The second effect contradicts the growth of stellar bars in DM halos with larger spin, where bars formed earlier with increasing λ (Long et al. 2014;Collier et al. 2018).
The bar pattern speeds in nonspinning halos are shown in the top row of Figure 5 and display a monotonic decay with time.Due to the increased mass concentration from B0 to B25, the initial pattern speeds decrease along this sequence.To avoid the noise, we only calculate the pattern speeds for  2 > ∼ 0.1.Hence, one can observe that the bar instability has been delayed along the B0 → B25 sequence as well.
We have calculated the evolution of bar sizes and present the example in Figure 6 for B25 models.These sizes are based on the extent of stable x 1 orbits on the characteristic diagrams -a family of orbits which form the backbones of stellar bars (e.g., Heller & Shlosman 1996;Collier et al. 2018;Li et al. 2023b).In this Figure we compare stellar bars in λ = 0 and λ = 0.09 DM halos.Note that while bars in the former halos grow monotonically, bars in fast spinning halos essentially grow only till the buckling, and their sizes remain constant after this, in agreement with Collier et al. (2018).Of course, bars in fast spinning bulges form later, and hence have a limited time to grow during the runs.
All stellar bar models in nonrotating halos experience a vertical buckling.This instability is delayed for massive bulges, in tandem  with the bar instability itself, and in models with a larger bulge spin (Figure 7).The buckling amplitude,  1 weakens with increasing bulge mass.
Some models display the secondary buckling of stellar bars, which has been first observed in simulations by Martinez-Valpuesta et al. (2006).The secondary buckling is well separated in time from the first buckling.While the first buckling occurs in the central region of a stellar bar, where the harmonic core lies, the secondary buckling happens at larger radii and proceeds much slower than the first one (see Figures 2,.The buckling timescales of the first and secondary bucklings are ∼ 0.5 Gyr and 2 − 3 Gyr, respectively, so they differ substantially.While secondary buckling is not the subject of this work, we note nevertheless that it  can be seen in B10-S75, B10-S95 and B17-S95 bulge models with a nonrotating halo.We use the vertical asymmetry parameter,  asym , introduced in section 2.3.1 to clarify some prime and secondary bucklings, especially in S75 and S95 bulge models.Figure 8 can be directly compared to Figure 7.Both types of bucklings can be now clearly verified as the noise has been reduced dramatically in  asym .We still observe that the buckling amplitude has decreased with increasing bulge mass.The B0 models exhibit the maximal buckling.What is interesting and was immersed in the noise for  1 , is the prolonged timescale of the secondary buckling with  asym .In all three models, the length of this buckling is almost 3 Gyr.In other words, this is the typical extent of the vertical asymmetry in secondary buckling.In fact, this timescale agrees perfectly with Martinez-Valpuesta et al. ( 2006).

Evolution of disk bar in spinning halos
As a next step, we discuss the evolution of stellar bars with bulges immersed in spinning halos.Figure 4 exhibits the evolution of  2 Fourier parameters of the disk.A number of conclusion can be made here.First, the maxima of  2 just before the buckling appear lower than in nonspinning halos, even for the B0 model.This is in ac-cordance with previous simulations (Long et al. 2014;Collier et al. 2018).With the exception of the B0 model, all models form rather an oval distortion than a bar with  2 ∼ 0.2, which represents a reasonable boundary between bars and oval distortions (e.g., Bi et al. 2022).
Second, the subsequent evolution of  2 also differs substantially from those in nonspinning halos.The B0 model displays a prolonged plateau after buckling, for about 6.5 Gyr, as observed already by Li et al. (2023b).Other models which experience buckling show this plateau as well.The length of this plateau decreases with increasing bulge mass.
Next, the models B0 and B10-S50 are the only ones that experience vertical buckling, at ∼ 8 Gyr and 6 Gyr, respectively.B17 and B25 and the rest of B10 models do not buckle.The B10-S50 buckling amplitude is very low,  1 ∼ 0.035.The bar in B25 model shows a very slow growth rate and only reaches  2 ∼ 0.2 late in the run.B25-S75 and B25-S95 reach this value only after 10 Gyr.
To verify the authenticity of the buckling, we turn to the alternative method described in section 2.3.1.This method appeared to be more sensitive to the vertical asymmetry as shown by the models in the nonrotating halos (section 3.1).The buckling in B10-S50 has been confirmed (see the lower row in Figure 8).The  asym parameter also underlines the prolonged width of the buckling in this mode, ∼ 4 − 5 Gyr!
The alternative method finds a clear buckling asymmetry in B10-S75, B17-S75, B10-S95 and B17-S95, where the standard Fourier method (Eq.( 8)) found no buckling at all.Moreover, the width of the buckling is characterized by a long timescale, ∼ 4 − 5 Gyr as well.The amplitude of the buckling is similar in all these models.For the bulge models that experience buckling, this instability develops fast but decays much more slowly than that in the models of nonspinning halos.The B25 models show no buckling in this method as well.Using the alternative method, displays a lower buckling amplitude in spinning halos compared to the nonspinning ones.
We emphasize that the height of the peaks in Figure 7 and Figure 8 have a different meaning.The  1 Fourier parameter detects the asymmetry in the edge-on disk density distribution.While,  asym parameter detects the difference in the face-on density, between  > 0 and  < 0. Additional question is why we do not observe the prolonged peaks in  1 while they appear in the  asym .This happens because  1 by construction is more sensitive to the central region of the disk, while  asym is equally sensitive at all radii.
In summary, the stellar bar evolution in spinning halos with the mass sequence of bulges differs dramatically from those embedded in the nonspinning halos.The stellar bar amplitude  2 which is substantially lower than for nonspinning halo and lingers around  2 ∼ 0.2 which represents rather an oval distortion.

Vertical asymmetry of stellar bars using the alternative method
To understand and visualize the alternative method to determine the vertical asymmetry more fully, we present the face-on snapshots of the disk for all the models with and without the bulge.Figure 3 displays the face-on stellar disks and bulges for B10-S50 to B25-S50 models immersed in nonrotating DM halos.All the snapshots are given at  ∼ 6.25 Gyr, after the stellar bars have reached the amplitude close to the maximal one.
We test the method of Smirnov & Sotnikova (2018) for sensitivity to the vertical buckling.The traditional method of measuring the buckling amplitude  1 has difficulty to detect the instability in some of the models with nonspinning and spinning halos.Using this alternative method, the noise is dramatically suppressed.This allows to recognize the buckling process more clearly.For example, for the bulgeless model B0 with a spinning halo,  1 ∼ 0.075, while the alternative method results in ∼ 1.25 in Figure 8.Moreover, the B10 nonspinning bulge within the halo of λ = 0.09, has  1 ∼ 0.035, while the new method results in 0.65 buckling amplitude.Factor of 20 higher leads to a clear detection of the buckling in B10 and B17 intermediate and high rotating bulge.Not only the buckling amplitude is increased in the new method, but also the duration of the buckling points a slowly restoring vertical symmetry over a few Gyr due to a higher sensitivity of the method.
Most interesting is observing the actual response of the disk to the buckling in Figure 9, for a bulgeless B0 model in a nonrotating halo.While the left frame shows the usual surface density distribution with a fully developed stellar bar, the middle and right frames display the separated disks for  > 0 and  < 0 surface densities.The  > 0 frame differs from the full surface density in the left frame by exhibiting a dumbbell shape in addition to having a very elongated distribution.At the same time, the  < 0 frame displays no dumbbell and being inclined by about 30 • to the bar which is horizontal.Its ellipticity is much smaller than that of the dumbbell frame.In this Figure, the bar is horizontal and the disk rotates counter-clockwise.
One can understand the appearance of the dumbbell shape during downward buckling -such a buckling is always associated with a double maxima response in the opposite direction, i.e., formation of the bell-shape and a dumbbell.As an example, we point to the image of an edge-on disk shown at the time of buckling downwards at  = 2.4 Myr in Figure 3 of Martinez-Valpuesta et al. (2006).What is new and unexpected is that the response above and below the midplane of the disk different orientation in the  −  plane.Based on the bar position and the disk rotation, we conclude that the response of the  < 0 part of the disk is trailing the bar and the response of the upper hemisphere.
Such a response will induce the shear across the disk midplane.This phenomenon has been indeed observed by Li et al. (2023a) in their Figures 6-9.They have observed formation of circulation cells and the associated vorticity in response to the shear.The associated stellar orbits giving rise to this circulation have presented in the Appendix of Li et al. work, and also discussed by Heller & Shlosman (1996).
Figure 10 describes the same phenomenon for the B17-S75 model embedded in the λ = 0.09 halo.Here we do not observe the misaligned response between the upper and lower disk hemispheres.The difference with the Figure 9 is that now a B17-S75 bulge is present and the DM halo has a spin.The amplitude of the asymmetric response is much lower,  asym ∼ 0.6 compared to ∼ 1.35 in the B0 model.Hence, the bulge acts to damps this asymmetry and this damping ability is expected to be stronger for more massive bulges.
We have tested the appearance of the shear and the associated phenomena for all models, but show only two models described above.To summarize these results, the shear and misaligned response has been detected in all models with nonspinning DM halos.Only B10-S50 models with λ = 0.09 halo exhibits this behavior.The rest of the models do not display it.The presence or absence of this behavior corresponds to appearance of  asym peak in Figure 8.The  1 Fourier coefficient which is typically used to measure the buckling strength of stellar bars is not so sensitive to buckling.
Using the  asym amplitude can be supplemented by mapping of the Laplace plane for the barred disk (e.g., Dekel & Shlosman 1983;Binney & Tremaine 2008).In Figure 11, we present an example of B10-S50 model during the maximal buckling, when the twin negative peak formed along the stellar bar (horizontal), and a positive peak  along the minor axis.This configuration is supported by the formation of a circulation cell of compression/stretching during the buckling analyzed in Li et al. (2023a).
We do observe that the  > 0 forming a bar is shorter than the  < 0 response, about 5 kpc versus 6 kpc radius bar, respectively.

Evolution of the bulge in the presence of stellar bar
We turn to the bulge models, both nonspinning and spinning initially.We start by discussing the evolution of their  = 2 modes in -plane.For the nonspinning bulges, this mode grows in response to the stellar bar development.But for bulges with an initial spin, the  = 2 mode can be triggered spontaneously, at least in principle, and evolve not in isolation but in interaction with the large-scale stellar bar.In all cases, the bulge and the bar possess gravitational quadrupole moments, and therefore can exchange their angular momenta both between them and with the host DM halo, again in principle.
Figure 12 displays the evolution of the bulge  2 Fourier coefficient for bulges embedded in λ = 0 nonspinning halos (top row) and λ = 0.09 spinning halos.For nonspinning bulges, the amplitude of  = 2 remains low,  2 < ∼ 0.1, which in our notation would correspond to the oval distortion.
The bulge  2 grows in response and completely mimics the development of the bar instability in the stellar disk (left column of Figure 12, S50 models).It responds by weakening during the buckling of the parent stellar bar (if this occurs).In nonrotating halo, the maximal bulge  2 barely reaches ∼ 0.1, while in the spinning halo, it is even weaker, again mimicking the evolution of the parent stellar bar.
For the middle column, the bulge S75 models, develop a much stronger  = 2 mode.Yet, it can be barely considered a "bulge bar" because its amplitude is ∼ 0.2.Again it seems that the bulge  = 2 mode has been induced by the parent bar instability in the stellar disk.Finally, the right column displays development of a strong bulge bar in S95 bulge models, reaching the amplitude of ∼ 0.3.13, of each component, namely, the DM halo (dashed lines), disk (solid lines), and bulge (dotted lines), normalized by the disk angular momentum of the bulgeless model B0 at  = 0 Gyr, in B0 (red), B10 (blue), B17 (black), and B25 (orange) models.The upper row displays the models with λ = 0 halos and the bottom row shows models immersed in λ = 0.09 halos.Note that the -axis is in linear scale in the range [ −10 −2 , 10 −2 ], otherwise in logarithmic scale.
In the -plane the bulge responds to the stellar bar, but in the vertical plane, it flattens in response to rotation and the disk potential (Figure 13).Thus the resulting Figure of a spinning bulge is that of a triaxial ellipsoid.This flattening is usually measured by ellipticity, i.e.,  = 1 − /, where  is directed along the -axis, and  and  are the semimajor and semiminor axes in the -plane (see section 2.4 for our method to calculate these axes using the moment of inertia of the bulge).

Angular momentum redistribution in the disk-bulge-halo system
As a next step we analyze the redistribution of angular momentum between the three components, namely, the bulge, disk and parent DM halo.We start by following the evolution of angular momentum in each of these components, then following the resonant interactions between them.Figure 14 exhibits the angular momentum gain and loss by each of the three morphological components, namely, We have normalized the angular momenta of each component by the initial angular momentum of stellar disk,  B0,disk ( = 0), of the bulgeless model B0.
In all cases, the lion share of angular momentum redistribution happened between the disk and its parent halo.The former lost and the latter gained the angular momenta, i.e., the stellar disk is the source of the angular momentum in the system and the halo is its sink.This exchange is much larger in the models with λ = 0 halos (top row in this Figure ), by a factor of ∼ 10 compared to models with spinning halos, λ = 0.09.As formation of stellar bars requires a substantial loss of angular momentum in the bar region, its transfers mostly to the parent halo.
The decreased angular momentum transfer between the disk and the halo explains in part the weaker bars forming inside such halos (e.g., Figure 4).In addition, this Figure also points to a difference between the bar strength in the bulgeless B0 model and B10-B25 models, and therefore to an additional angular momentum transfer, this time between the stellar disks and their bulges.
Our next question is to determine the role of the bulge in the angular momentum redistribution.All nonspinning bulges in λ = 0 halos gain angular momentum, but this amount is small and on the linear scale compared to the angular momentum gained by the halos on the logarithmic scale.All the bulges inside λ = 0 halos gain angular momentum until the associated stellar bars have reached the maximal amplitude at the buckling.After this, the transfer of angular momentum has almost ceased until the end of the run at 10 Gyr.For initially nonspinning bulges immersed inside λ = 0.09 halos, the bulges continue to gain angular momentum till the end of the simulations, and the saturation effect is not visible almost at all.
Overall, all initially nonspinning bulges gained angular momentum by a factor of ∼ 50 less than their parent halos.Interesting that more massive bulges in this left column of Figure 14 have gained more angular momentum.We have checked and found that this result holds for the specific angular momentum of bulge particles as well.
For initially intermediate spin bulge, i.e., S75 models, the change of angular momentum is basically non-existent, independent of the bulge mass.For the initially fast spinning bulges, i.e., S95 models, in λ = 0 halos, the bulges appear to lose their angular momentum, although compared to the parent stellar disks, this amount is not significant.The same bulges in λ = 0.09 halos, in contrast, lose very little angular momentum.
Based on our understanding of the source and sink of angular momentum in the bulge-disk-halo system, we turn to the rate of the angular momentum transfer to and away from the bulge, attempting to answer the question, what is the role of the bulge in this evolution.We use the method developed by Villa-Vargas et al. (2009), see also Collier et al. (2018) and Li et al. (2023b).The cylindrical radius has been binned into cylindrical shells with the height of || = 6 kpc.We calculate the angular momentum change in these shells,  (), and obtain the rate of change,  = [ ( + Δ) −  ()]/Δ.
The left panel in Figure 15 exhibits the rate of change of  in the all the bulge models immersed in the λ = 0 halos.It must be compared with Figures 12 and 14.The B0 model can be found in Li et al. (2023b).We observe that initially nonrotating bulges (S50 models) absorb the angular momentum almost all the time.After the associated stellar bar buckling, the angular momentum absorption rate decreases sharply and some emission of  appears.This explains that the termination of the bulge spinup saturates in Figure 14.
For S75 models, the  redistribution differs.The inner bulge emits its angular momentum at small radii and absorbs it at larger radiias a result almost no net -transfer occurs.This absorption increases even more in S95 models, leading to the loss of  by these bulges.
The right panel in Figure 15 displays the same phenomenon in λ = 0.09 halos.Pure spin-up of bulges in S50 models, appearance of absorption in the inner bulges in S75 models, and dominant absorption in S95 models.

Orbital spectral analysis
Finally, we performed the orbital spectral analysis of the bulge-diskhalo systems, with a particular emphasis on the bulge, in order to examine the role of resonances in angular momentum transfer among the morphological sub-components.We followed the method developed and used by us (Martinez-Valpuesta et al. 2006;Dubinski et al. 2009;Collier et al. 2019a).The dimensionless frequency ratio is defined as  = (Ω−Ω bar )/, where Ω and  are the angular and epicyclic frequencies characterizing the individual orbits.The frequency ratio  is binned in Δ = 0.005, and we choose a random sample of 50% of bulge particles from each model for this analysis.In this notation, the main resonances are the inner Lindblad resonance (ILR), the corotation resonance (CR), and the outer Lindblad resonance (OLR).They correspond to  = 0.5, 0, and −0.5, respectively.
Our goal is to find the fraction of bulge orbits trapped at particular resonances and whether these resonances gain or loose angular momentum.The orbital spectral analysis has been performed at different times before and after buckling of stellar bars in models with a nonrotating bulge and nonrotating halo.This allowed the measure of the angular momentum gain or loss at these resonances.The resulting distribution of bulge particles with  has shown that the main resonances which trap these particles are  = 0.5,  = 0, and  = −0.5.
We have calculated the trapping efficiency and angular momentum gain/loss for initially nonrotating bulges of B10, B17 and B25 embedded in λ = 0 halos.In these models, the bulge ILR, has always shown absorption of angular momentum, with an efficiency of trapping increasing with the bulge mass fraction from ∼ 14% to ∼ 24%.This is opposite to what we observe for the barred disks which always emit angular momentum at this resonance.The second resonance corresponds to the OLR with a trapping efficiency ∼ 10% to ∼ 15%, but it involves the retrograde orbits with Ω < 0. These orbits have shown an absorption of negative angular momentum.To understand this behavior, we looked at these orbits and found that they lie deep within the bulge and circulate in the direction opposite to that of the stellar bar and disk.They appear elongated orthogonal to the bar.
We have also detected small trapping of bulge orbits at  = −1, which corresponds to the Lagrange points.But the efficiency of trapping by this resonance has been low, ∼ 1%, and we did not follow it.In comparison, Saha et al. (2012) recorded an absorption of angular momentum at this resonance, although with a varying efficiency.
In our analysis, the CR resonance for the bulge and additional higher resonances have been found to trap a negligible fraction of bulge orbits, < ∼ 1%.This happens typically at very late times of the runs.The reason is that these bulge orbits lie deep inside the total gravitational potential well.

DISCUSSION
We run a suite of numerical simulations of stellar disks with bulges immersed in DM halos.A range of bulge masses of 10% -25% of the total disk + bulge mass and various degrees of initial rotation of the bulge and DM halos have been used in our models.
Our results are as follows: • In nonspinning DM halos, i.e., λ = 0, we find a very small decrease of the maximal amplitude of stellar bars before the vertical buckling.The decrease in bar strength during the buckling, Δ 2 , however, becomes smaller with increasing bulge mass, and so is the buckling amplitude,  1 , independently of the initial bulge spin.All modeled bars have restarted their growth after buckling and reached about the same amplitude of  2 after 10 Gyr.
• In contrast, in spinning DM halos with λ = 0.09, all models with the bulge mass fractions of 10 − 25% have reduced the prebuckling amplitude  2 from ∼ 0.3 to ∼ 0.2.The buckling amplitude,  1 , has been reduced dramatically in all these models by a factor of ∼ 2 compared to the bulgeless model.Moreover, some models do not show any buckling of stellar bars, in agreement with Sellwood & Gerhard (2020), although our bulge parameters differ.Furthermore, stellar bar amplitude does not recover after buckling in spinning halos, in agreement with previous results (Long et al. 2014;Collier et al. 2018;Li et al. 2023b) • Stellar bar instability is delayed both in more massive bulges and even more in faster spinning bulges, i.e., along the bulge mass sequence B0 → B25, and along the bulge spin sequence S50 → S95.The stellar bars form with higher pattern speeds in higher mass bulges and even faster in bulges with a higher spin.We reiterate, as noted above, the formation of stellar bars is delayed for more massive and faster spinning bulges.This behavior is irrespective of the parent DM halo spin.
• In order to verify that indeed the buckling of stellar bars has been suppressed in some models with spinning halos, we have used a modified method introduced by Smirnov & Sotnikova (2018), and discussed here in section 2.3.1.This method has been found to be much more sensitive to the buckling process.Moreover, it clearly demonstrated that the buckling amplitude in spinning halos,  asym , although reduced by a factor of ∼ 2 in bulges within the analyzed mass fraction, has a prolonged duration, as confirmed by a 2-D Laplace plane distortion map.In fact, stellar bars in these models have never fully restored their vertical symmetry within the run time.We show that this behavior emphasizes the different morphology above and below the bar midplane.Note also that the most massive bulge model B25 in λ = 0.09 halo has never buckled, because the timescale of the bar instability is ∼ 10 Gyr, irrespective of the bulge spin.
• The angular momentum flow in the bulge-disk-halo systems is heavily dominated by the barred disk -DM halo interactions.The lion share of angular momentum flow is from the barred disk to DM halo, as is well known from previous investigations.In this triple system, the bulges play a secondary role in -transfer.Yet, the initially nonrotating bulges, S50, absorb the angular momentum from the barred disk, the intermediate rotating bulges, S75, do not change their spin, while the fast spinning bulges, S95, slowdown their initial spin.This conclusion applies both for λ = 0 and λ = 0.09 halos equally.
• By applying the orbital spectral analysis for the bulge orbits, we find two competing resonances, at  = 0.5 and  = −0.5 trapping the bulge orbits, where  = (Ω bar − Ω)/ is the dimensionless frequency ratio.The former resonance is the inner Lindblad resonance (ILR) and is the dominant one.The latter resonance involves retrograde bulge orbits only (with respect to the disk and bulge rotations), elongated perpendicular to the stellar bar.In this respect they mimic the outer Lindblad resonance (OLR) orbits, but they circulate in the opposing direction.The efficiency of this latter resonance in trapping the bulge orbits and, therefore, in transferring the angular momentum, varies with time significantly.Furthermore, both trapped prograde (at ILR) and retrograde (at OLR) bulge orbits have been found to gain angular momentum during the evolution.In S50 models, most of the angular momentum transfer to/from the bulge proceeds via these trapped prograde and retrograde bulge orbits.Their angular momentum becomes more positive or negative, respectively.
• In S95 models, orbits trapped at the ILR gain angular momentum.Orbits trapped at the OLR, are retrograde and gain negative angular momentum.Moreover, orbits trapped in the range of −0.5 <  < 0, gain negative angular momentum.While orbits in the range of 0 <  < 0.5 lose positive angular momentum.In balance, for S95 models, the angular momentum is lost.
We now return to the focal aspect of this work.Namely, how does the bulge affect the stellar bar evolution and how does the bar reciprocate.For this purpose, we constructed sequences of the / mass ratios and the initial bulge spin.Additional important parameters is the spin of the parent DM halo.Previous works have studied many aspects of this bulge-disk-halo system.Our main contribution in this work is to quantify the effect of the spin in both spheroidal components, understanding the flow of angular momentum between these components, and especially the role of the orbital resonances in the interaction between them.
An interesting outcome of our analysis of these simulations is the application of the modified method of Smirnov & Sotnikova (2018) to quantify the vertical asymmetry of stellar bars.Not only it provides a more sensitive way to detect the vertical buckling of bars, but it reveals the long-term effect of the buckling, which has been found to extend for additional few Gyr in spinning halos.This result is supported by calculation of the distorted Laplace plane.Such a response of a barred disk was not observed before in numerical simulations and brings forward the question of what is the underlying mechanism of this prolonged violation of the vertical symmetry.
To demonstrate the prolonged buckling process in spinning halos, we refer to the face-on Laplace plane in B10-S50 model immersed in λ = 0.09 Dm halo (Fig. 16).Note that 1 Gyr after the buckling, the double peak is still visible at slightly reduced amplitude.
The first application of the orbital dynamics to the formation of peanut/boxy shapes of galactic bulges has determined the population of the so-called BAN/ABAN orbits (Pfenniger & Friedli 1991).These are 3-D orbits and their origin lies in the action of the combined horizontal ILR and vertical ILR (vILR).These orbits are 2:1 orbits (two vertical oscillations per one turn in the plane) and bifurcate vertically from the x 1 family, and in projection on the  − -plane they look like the x 1 orbits.Particles populating these orbits are involved in the vertical oscillations which take them away from the disk midplane.Due to the vILR occurring at some point of the x 1 orbit family, a vertical orbit family fork-bifurcation forms, and the possibility of a transverse break of symmetry occurs.Triggered by noise, symmetry may be broken, one side of the BAN/ABAN family pair is then more populated than the other, and a -asymmetry grows.This vertical instability grows until saturation in the non-linear regime.Under these conditions, the Laplace plane, which is determined by the potential minimum along the -axis or, alternatively, by the zeroacceleration along this axis, is bent away from the disk midplane.
The behavior of this population of orbits is closely related to the appearance of KAM surfaces 1 which partly confine the chaotic motion.In 2-D systems (i.e., systems with two degrees of freedom) these KAM surfaces isolate the phase space volumes.But they cannot do this perfectly in non-integrable systems with more than two degrees of freedom, due to the Arnold's diffusion (e.g., Lichtenberg & Lieberman 1992).When close to the local integrability, the Arnold's diffusion can be slow.
Both BAN and ABAN orbits populate different regions in phase space.When the vertical asymmetry of the potential exists, the populations of these orbits initially differ, but slowly diffuse through chaotic phase space by Arnold's diffusion until they become equal.This process can support the vertical asymmetry for prolonged time periods as discussed above (see also Figures 8 and 9).
Next, we focus on the corollaries of a bulge presence in disk galaxies.We find that whether the bulge is a recipient or donor of angular momentum in the system depends on its initial spin.Previous studies have shown that a nonrotating bulge absorbs the angular momentum from the barred disk via gravitational torques (e.g., Saha et al. 2012;Saha & Gerhard 2013;Saha et al. 2016).Such bulge correspond to our S50 models.Figure 14  Moreover, it shows that more massive bulges absorb more angular momentum.
However, our S75 models, with an intermediate rotating bulges, show that they both emit and absorb the angular momentum (Figures 14 and 15), either in the DM halos with λ = 0 or 0.09.S95 models, with fast spinning bulges, actually slow down by losing their angular momentum.
For the range of / in our models and irrespective of the initial spin of the bulge, the amount of angular momentum absorbed or emitted by the bulge is small compared to the angular momentum exchange between the barred disk and the DM halo.This conclusion is expected as the maximal amount of angular momentum that can be stored in the bulge of fixed shape is negligible compared to the disk or DM halo (see Tables 2 and 3, section 2.2 and the Appendix).Hence, the angular momentum stored in the bulge is not expected to directly affect the stellar bar evolution.However, the bulge can affect the disk evolution by other means, e.g., by destroying the harmonic motion in the core which plays an important role in the buckling process.Furthermore, the bulge can provide an initial perturbation since its dynamical timescale is the shortest in the system.Under these conditions, the development of a  = 2 mode there can affect the development of this mode in the larger disk.
The presence of a bulge affects the pattern speed of a bar, as was investigated by Kataria & Das (2019), concluding that the slowdown is faster for more massive bulges.We confirm this result, but have noticed that the bars form also with higher Ω bar in the presence of more massive bulges (e.g., Figure 5).The net slowdown of bars is such that the Ω bar for a spectrum of bulge masses never intersect.So, bars in models with higher bulge masses will end up with a higher Ω bar even after a 10 Gyr run.Note also that the bar pattern speeds in spinning halos become flat quickly, stopping their decay, as expected (Collier et al. 2018).
Figure 12 displays the evolution of  2 for S50, S75 and S95 bulge models.When compared with the stellar bar amplitudes in Figure 4, we observe a one-to-one correspondence.Yet, for S50 models, the bulge  = 2 amplitude is smaller by a factor of 4 than the bar  = 2 amplitude.For S75 models, this difference is only a factor of 2. And for S95 models, there is almost no difference between the  = 2 mode of the bulge and the parent stellar bar for parent DM halos with λ = 0.
For spinning halos with λ = 0.09, already S75 bulge models have similar amplitudes of  = 2 modes with the bar.And for S95 models, the bulge  2 is stronger compared to the associated bar by about 50%.
The spin and mass sequences of bulges have clear effects on the bulge shapes, both in the  −  and vertical planes.In the  −  plane, the bulge shape measures its strength via the amplitude of the  = 2 mode.This shape can be alternatively measured using ellipticity (see section 2.4 for definitions).
Figure 13 displays the bulge ellipticity evolution in the vertical plane.In both nonspinning and spinning halos, the bulge becomes more eccentric with increasing initial spin reaching  ∼ 0.3 for S95 models.This appears to be the maximal ellipticity reached by a bulge nearly maximally spinning, corresponding to the axial ratio of / ∼ 0.7.Such axial ratios are observed in maximally flattened bulges and elliptical galaxies (e.g., Binney & Merrifield 1998).

CONCLUSIONS
We have studied the evolution of barred disks immersed in spheroidal components which include DM halos and stellar bulges.The initial conditions include a sequence of bulge masses, as well as sequences of bulge and halo spins.The analyzed / ∼ 0 − 25% mass ratio range, corresponds to the one typically observed in disk galaxies, while the DM halo spin range, λ ∼ 0 − 0.09 to the most probable one based on the numerical simulations.
We find that the presence of a bulge affects evolution of stellar bars, including the characteristic timescale of the bar instability, the bar pattern speed and its decay, bar length, and its vertical buckling instability.The main effect of the bulge is the destruction of the harmonic core in the disk which affects the vertical buckling instability in the bar, as shown by Li et al. (2023a).Moreover, the bulge does not preclude the secondary buckling of stellar bars, as expected, as it happens away from the central regions dominated by the bulge (Martinez-Valpuesta et al. 2006).However, as a reservoir or sink of angular momentum, the bulge plays a minor role, and the angular momentum exchange proceeds mainly between the DM halo and the barred disk.
In the modeled range of / ∼ 0 − 25% in nonspinning halos, we do not find any damping of the amplitude of the bar instability, but find a progressively increasing delay in the characteristic timescale of instability.In contrast, within spinning halos with λ = 0.09, the bar amplitude has decreased by ∼ 1/3 compared to halos with λ = 0.
We confirm the previous result by Kataria & Das (2019) that pattern speeds of newly formed bars appear higher for more massive bulges.Moreover, we find that, although such bars formed with higher Ω bar and experiencing a slowdown, will always end up with higher pattern speeds even after a 10 Gyr.But we do not find such a sensitivity of the stellar bar amplitude to the presence of the bulge as claimed by Saha & Gerhard (2013) and Sellwood & Gerhard (2020).
We have applied a modified version of the parameter proposed by Smirnov & Sotnikova (2018) measuring the vertical buckling amplitude,  asym , and found it to be substantially more sensitive to buckling compared to the Fourier  1 parameter used in the literature.Using  asym , we find that the buckling amplitude has been decreased by a factor of ∼ 2 in halos with λ = 0.09 compared to halos with λ = 0.
Most interestingly, we find that the buckling process triggers different responses above and below the disk midplane.This difference can be diluted by progressively more massive bulges.In spinning halos, the buckling distortion has a prolonged amplitude tail which extends by a few Gyr.We have tested this phenomenon and verified it by using the distortion of the Laplace plane.We relate this pro-longed asymmetry with the evolution of the population of the 3-D BAN/ABAN orbits and slow Arnold's diffusion in our systems.This prolonged asymmetry can potentially have observational corollaries.
The angular momentum flow in the barred disk-bulge-halo systems proceeds mainly between the barred disk and the DM halo, with bulge in the mass range of / ∼ 0−25% playing a minor role.As verified using orbital spectral analysis, non-spinning bulges gain their spin from the stellar bar, and do it mainly via the ILR, while fast-spinning bulges lose their spin via prograde orbits which lie between the CR and the ILR, and retrograde orbits trapped by the OLR.

Figure 1 .
Figure1.Initial density profiles for the DM halo (the NFW profile, blue dashed line) and the bulge (the Plummer sphere, dashed-dot line) models.also shown is the total mass profile (continuous black line).The disk + bulge masses have been kept the same in all models.

Figure 3 .
Figure3.Disk surface density contours for all nonrotating bulge models embedded in the λ = 0 halo, at  = 6.25 Gyr when the disk  2 are similar.The color corresponds to the surface density and is given by the logarithmic palette.The isodensity contours are separated by a factor of 1.5 in surface density, and the unit is in 10 10 M ⊙ kpc −2 .The stellar bar is horizontal and the disk rotates counter-clockwise.The vertical dashed lines have been introduced to compare the bar sizes.

Figure 5 .
Figure 5. Evolution of the stellar bar pattern speeds, Ω bar .

Figure 6 .
Figure 6.Evolution of bar size of B25 models.

Figure 7 .Figure 8 .
Figure 7. Evolution of the normalized Fourier amplitude  1 of the disk, using the traditional method.

Figure 9 .Figure 10 .
Figure 9. Left frame:The disk surface density contours of the B0 model embedded in the nonrotating halo with λ = 0 at the time of maximal buckling,  = 3.35 Gyr.Middle and right frames: same as left frame but separating  > 0 and  < 0 surface densities.The  = 0 plane is defined in Eq. (10) and the color palette is in the same level as Figure3.

Figure 11 .
Figure11.The face-on Laplace plane for the model B10-S50 with the λ = 0 halo, at the time of the maximum buckling,  = 3.30 Gyr.The major axis stellar bar is horizontal.During this buckling, the stellar bar center bends upward.

Figure 14 .
Figure 14.Evolution of the angular momenta exchange, Δ  in Equation13, of each component, namely, the DM halo (dashed lines), disk (solid lines), and bulge (dotted lines), normalized by the disk angular momentum of the bulgeless model B0 at  = 0 Gyr, in B0 (red), B10 (blue), B17 (black), and B25 (orange) models.The upper row displays the models with λ = 0 halos and the bottom row shows models immersed in λ = 0.09 halos.Note that the -axis is in linear scale in the range [ −10 −2 , 10 −2 ], otherwise in logarithmic scale.

Figure 15 .
Figure 15.The rate of the angular momentum change in the bulge,   , in the models with λ = 0 halo (left) and λ = 0.09 halo (right).  is calculated by every 0.05 Gyr within cylindrical shells Δ = 1 kpc and |  | < 6 kpc.The red and blue color palette corresponds to the gain and loss of angular momentum, respectively.  is given in units of 10 10 M ⊙ kpc 2 Gyr −2 .

Figure 16 .
Figure16.Same plot as Figure11for the model B10-S50 in λ = 0.09 halo, at the time of the maximal buckling (left) and 1 Gyr later (right).Note that the buckling footprint is still visible as a double peak along the stellar bar (horizontal).

Table 3 .
Initial angular momentum of halo, disk, and bulge in all models.The units are in 10 13 M ⊙ kpc km s −1 .