# Effects of Halo Spin on the Formation and Evolution of Bars in Disk Galaxies

Dajeong Jang Department of Physics $\&$ Astronomy, Seoul National University,
Seoul 08826, Republic of Korea  Woong-Tae Kim Department of Physics $\&$
Astronomy, Seoul National University, Seoul 08826, Republic of Korea  SNU
Astronomy Research Center, Seoul National University, Seoul 08826, Republic of
Korea

###### Abstract

The spin of dark halos has been shown to significantly affect bar formation
and evolution in disk galaxies. To understand the physical role of the halo
spin on bar formation, we run $N$-body simulations of isolated, Milky Way-
sized galaxies by varying the halo spin parameter in the range
$-0.16\leq\lambda\leq 0.16$ and the bulge mass. We find that our adopted halo
_alone_ is subject to swing amplification of an $m=2$ non-axisymmetric mode
rotating in the same sense as the halo, which assists or inhibits the bar
formation in a disk depending on its sense of rotation. The $m=2$ mode in the
disk, growing via swing amplification, interacts constructively
(destructively) with the $m=2$ mode in the prograde (retrograde) halo,
promoting (delaying) bar formation. A bar grows by losing its angular momentum
primarily to a halo. Since the halo particles inside (outside) the corotation
resonance with the bar can emit (absorb) angular momentum to (from) the bar,
the bar pattern speed decays slower for larger $\lambda>0$, while it decreases
relatively fast almost independent of $\lambda\leq 0$. Models with a strong
bar develop a boxy peanut-shaped bulge. In models without a bulge, this occurs
rapidly via buckling instability, while the bars with a bulge thicken
gradually without undergoing buckling instability. Among the models considered
in the present work, the bar in the $\lambda=0.06$ model with a bulge of 10%
of the disk mass best describes the Milky Way in terms of the bar length and
pattern speed.

Disk Galaxies (391), Milky Way Galaxy (1054), Galaxy dark matter halos (1880),
Galaxy Bulges (578), Galaxy Disks (589), Barred Spiral Galaxies (136), Galaxy
Bars (2364)

BPS

    boxy peanut-shaped

##  1 Introduction

Nearly $2/3$ of disk galaxies in the local universe possess a weak or strong
bar in optical and near-infrared images (de Vaucouleurs, 1963; Whyte et al.,
2002; Laurikainen et al., 2004; Marinova & Jogee, 2007; Menéndez-Delmestre et
al., 2007; Aguerri et al., 2009; Méndez-Abreu et al., 2012; Buta et al., 2015;
Díaz-García et al., 2016, 2019). Bars are one of the main drivers of the
dynamical and secular evolution of disk galaxies, responsible for stellar
migration (Di Matteo et al., 2013; Kawata et al., 2017; Halle et al., 2018;
Iles et al., 2024), angular momentum transport from disk to halo (Saha & Naab,
2013; Petersen et al., 2016; Collier et al., 2018; Collier & Madigan, 2021),
and formation of pseudo bulges (Combes & Sanders, 1981; Combes et al., 1990;
Athanassoula, 2005; Debattista et al., 2006; Gadotti, 2011). Despite the
importance of bars, however, there has been no clear answer as to why some
galaxies are barred or non-barred and what controls the physical properties of
bars.

Many numerical studies have attempted to find one-parameter conditions for bar
formation using the initial galaxy properties such as the ratio of the kinetic
to gravitational potential energy and the ratio of the bulge to disk mean
density (Ostriker & Peebles, 1973; Efstathiou et al., 1982; Kataria & Das,
2018; Saha & Elmegreen, 2018). While the proposed conditions explain the bar
formation in their own adopted models, they fail to be applicable to diverse
galaxy models. The most likely reasons for this may be that the proposed
conditions are based on a limited set of galaxy models by fixing either bulge
or halo parameters and that the one-parameter conditions cannot capture the
complicated processes of bar formation. Recently, Jang & Kim (2023, hereafter
JK23) proposed the two-parameter condition

$\left(\frac{{Q_{T,\text{min}}}}{1.2}\right)^{2}+\left(\frac{\text{CMC}}{0.05}\right)^{2}\lesssim
1$ (1)

for bar formation. In Equation 1, ${Q_{T,\text{min}}}$ is the minimum value of
the Toomre stability parameter of the disk and CMC denotes the central mass
concentration. This condition is physically motivated: bar formation requires
several cycles of swing amplification and feedback loops (Sellwood, 1980;
Toomre, 1981; Binney & Tremaine, 2008), and perturbations grow more easily in
a disk with smaller ${Q_{T,\text{min}}}$ and CMC. Equation 1 implies that the
presence of a (massive) halo tends to suppress bar formation by increasing
both ${Q_{T,\text{min}}}$ and CMC, which is qualitatively consistent with the
established result that the density of a halo has a negative influence on the
bar formation (Ostriker & Peebles, 1973; Hohl, 1976; Efstathiou et al., 1982;
Debattista & Sellwood, 2000).

While Equation 1 accounts for bar formation in diverse galaxy models with
various degrees of the central concentration of halo and bulge, it is limited
to models with a non-rotating halo. Cosmological simulations of galaxy
formation found that dark halos acquire angular momentum from tidal torque
(e.g., Hoyle, 1951; Fall & Efstathiou, 1980; White, 1984) and satellite
accretion (e.g., Vitvitska et al., 2002; Dillamore et al., 2023). As a
dimensionless measure of halo spin, Bullock et al. (2001) introduced

$\lambda\equiv\frac{L_{\mathrm{vir}}}{\sqrt{2}M_{\mathrm{vir}}r_{\mathrm{vir}}V_{c}},$
(2)

where $L_{\rm vir}$ and $M_{\mathrm{vir}}$ are the angular momentum and mass
within the virial radius $r_{\mathrm{vir}}$, respectively, and
$V_{c}=(GM_{\mathrm{vir}}/r_{\mathrm{vir}})^{1/2}$ is the circular velocity at
$r_{\mathrm{vir}}$.111Equation 2 was derived from
$\lambda_{P}=L|E|^{1/2}/(GM^{5/2})$, where $L$, $E$, and $M$ refer to the
total angular momentum, energy, and mass of a spherical system within radius
$r$, respectively (Peebles 1969, 1971; see also Ansar & Das 2023). At the
virial radius where $E=-GM_{\mathrm{vir}}/(2r_{\mathrm{vir}})$, $\lambda$ in
Equation 2 becomes identical to $\lambda_{P}$. Note that positive (negative)
$\lambda$ corresponds to a spinning halo in the same (opposite) sense as the
disk rotation. Cosmological simulations showed that $\lambda$ follows a
lognormal distribution, with the mean value $\lambda\simeq 0.03$–$0.04$
(Bullock et al., 2001; Hetznecker & Burkert, 2006; Bett et al., 2007; Jiang et
al., 2019). The Milky Way halo is also thought to spin with $\lambda=0.061$
(Obreja et al., 2022).

$N$-body simulations for bar formation with a spinning halo commonly found
that a stellar disk embedded in a prograde halo with larger $\lambda\;(>0)$
develops a bar faster (Saha & Naab, 2013; Long et al., 2014; Collier et al.,
2018; Collier & Madigan, 2021; Kataria & Shen, 2022; Li et al., 2024), while
bar formation in a disk under a retrograde halo with smaller $\lambda\;(<0)$
is progressively delayed (Collier et al., 2019a). To explain this, Saha & Naab
(2013) invoked the condition $t_{\rm OP}\equiv T_{\mathrm{rot}}/|W|>0.14$ of
Ostriker & Peebles (1973) for bar formation, where $T_{\mathrm{rot}}$ and $W$
refer to the rotational and gravitational potential energies of a galaxy,
respectively. They argued that the halo spin increases $T_{\mathrm{rot}}$ to
form a bar faster. Since Ostriker & Peebles (1973) employed models with a non-
rotating, fixed halo, however, it is questionable whether $T_{\mathrm{rot}}$
should include the disk rotation only or the halo rotation as well. Even if
the halo rotation can be included, $t_{\rm OP}$ is independent of the sense of
halo rotation, as Saha & Naab (2013) noticed, which is not compatible with the
numerical results: a bar forms earlier in a model with $\lambda\;(>0)$ than in
the $-\lambda$ counterpart. This makes $t_{\rm OP}$ inadequate for explaining
the effect of the halo spin on bar formation.

Once a bar forms, the angular momentum transfer between it and its surrounding
halo affects the bar growth significantly. When $\lambda=0$, a bar grows by
losing its angular momentum to a halo (Athanassoula 2002; Martinez-Valpuesta
et al. 2006; Kataria & Das 2018; Saha & Elmegreen 2018; JK23). The angular
momentum transfer appears more complicated when the halo spin is considered.
Long et al. (2014) found that the inner part of a prograde halo loses its
angular momentum to a bar, while the outer part acquires it. Collier et al.
(2019b) suggested that the angular momentum transfer by the resonances is more
active under a halo with higher $\lambda>0$. Collier et al. (2019a) and
Collier & Madigan (2021) proposed that orbital reversals of halo particles are
a crucial factor of angular momentum transfer between disk and halo with
$\lambda<0$. Kataria & Shen (2022) found that it is the spin of the inner halo
(with $R<30\;{\rm kpc}$) rather than the total halo spin that is important for
the bar formation and its ensuing evolution.

In this paper, we use $N$-body simulations to investigate the effects of halo
spin and the presence of a classical bulge on the bar formation and evolution.
Our primary purpose is to understand how the halo spin physically affects the
bar formation. For this, we consider Milky Way-sized galaxy models with
differing bulge mass and halo spin. To delineate the effect of halo spin, we
consider a wide range of the spin parameter $-0.16\leq\lambda\leq 0.16$, and
compare the results from halo-only models with $\lambda=\pm 0.16$ or 0. We
monitor the temporal changes of the bar strength and pattern speed, study the
vertical buckling of bars in our simulations, and analyze their dependencies
on the spin parameter.

This paper is structured as follows. In Section 2, we describe our galaxy
models and numerical methods. In Section 3, we run simulations without a disk
and bulge to study the gravitational susceptibility of our adopted halos. In
Section 4, we present the results of full galaxy models, focusing on the
dependence on $\lambda$ of the bar formation, evolution, pattern speed, and
angular momentum transfer. We also present the results for vertical buckling
instability of our models and its outcomes. In Section 5, we discuss our
results in comparison with the previous results and constrain the model
parameters for the Milky Way. Finally, we conclude our work in Section 6.

##  2 Galaxy Model and Method

Figure 1: Radial distributions of the total circular velocity $v_{c}$ (black
solid) as well as the contributions of a halo (blue dotted), a disk (red
dashed), and a bulge (green dot-dashed). The C00 and C10 series in the upper
and lower panels refer to the models without and with a classical bulge,
respectively.

To study the effects of halo spin on bar formation and evolution, we consider
isolated galaxy models whose mass and size are similar to those of the Milky
Way. Our galaxy models consist of a dark matter halo, a stellar disk, and a
central supermassive black hole. Some models do not have a classical bulge
(C00 series), but we also consider models with a bulge (C10 series).

The stellar disk has the conventional exponential-secant hyperbolic density
distribution

$\rho_{d}(R,z)=\frac{M_{d}}{4\pi
z_{d}R_{d}^{2}}\exp\left(-\frac{R}{R_{d}}\right){\rm
sech}^{2}\left(\frac{z}{z_{d}}\right),$ (3)

where $R$ is the cylindrical radius, $z$ is the vertical height, $R_{d}$ is
the disk scale radius, $z_{d}$ is the disk scale height, and $M_{d}$ is the
total disk mass. We take $R_{d}=3\;{\rm kpc}$, $z_{d}=0.3\;{\rm kpc}$, and
$M_{d}=5\times 10^{10}\;{\rm M}_{\odot}$, similar to the Milky Way (Bland-
Hawthorn & Gerhard 2016, Helmi 2020).

Table 1: Model parameters  Name |  $M_{b}/M_{d}$ |  $\lambda$ |  $f_{p}$
---|---|---|---
C00P16 |  $\cdots$ | $0.16$ | 1.00
C00P10 |  $\cdots$ | $0.10$ | 0.80
C00P06 |  $\cdots$ | $0.06$ | 0.67
C00P00 |  $\cdots$ | $0.00$ | 0.50
C00R06 |  $\cdots$ | $-0.06$ | 0.32
C00R10 |  $\cdots$ | $-0.10$ | 0.20
C00R16 |  $\cdots$ | $-0.16$ | 0.00
C10P16 | 0.1 | $0.16$ | 1.00
C10P10 | 0.1 | $0.10$ | 0.80
C10P06 | 0.1 | $0.06$ | 0.67
C10P00 | 0.1 | $0.00$ | 0.50
C10R06 | 0.1 | $-0.06$ | 0.32
C10R10 | 0.1 | $-0.10$ | 0.20
C10R16 | 0.1 | $-0.16$ | 0.00

Both the halo and classical bulge have a density distribution that follows
Hernquist (1990) profile

$\rho(r)=\frac{M}{2\pi}\frac{a}{r(r+a)^{3}},$ (4)

where $r=(R^{2}+z^{2})^{1/2}$ is the spherical radius, and $M$ and $a$
indicate the mass and the scale radius of each component, respectively. For
the halo, we fix its mass and scale radius to $M_{h}=1.3\times 10^{12}\;{\rm
M}_{\odot}=26M_{d}$ and $a_{h}=30\;{\rm kpc}$, similar to the Milky Way (e.g.,
Bland-Hawthorn & Gerhard 2016, Helmi 2020). For the bulge, we set its scale
radius to $a_{b}=0.4\;{\rm kpc}$ and mass to $M_{b}=0.1M_{d}$ (e.g., Shen et
al. 2010, Bland-Hawthorn & Gerhard 2016, Helmi 2020). We place a supermassive
black hole with mass $M_{\mathrm{BH}}=4\times 10^{6}\;{\rm M}_{\odot}$ at the
galaxy center (Ghez et al. 2008; Gillessen et al. 2009a, b). Figure 1 plots
the total circular velocity $v_{c}$ calculated from the gravitational
potential and the contribution of each component for the models in the C00 and
C10 series. Due to the bulge, $v_{c}$ in the C10 series increases rapidly at
small $R$ compared to that in the C00 series.

To assign the spin parameter $\lambda$ to the dark halo, we follow the method
of Long et al. (2014) and Collier et al. (2018, 2019a). That is, we randomly
select a fraction of the halo particles and change the directions of their
tangential velocities, while keeping their amplitudes and radial velocities
intact. While Kataria & Shen (2022) explored the models in which the fraction
$f_{p}$ of the halo particles on prograde orbits with respect to the disk
rotation varies with radius, we keep $f_{p}$ spatially constant for
simplicity. It turns out that our halo has $\lambda=0.16$ when $f_{p}=1$. The
prograde fraction $f_{p}$ decreases as $\lambda$ decreases, leading to
$\lambda=0$ for $f_{p}=0.5$, and $\lambda=-0.16$ for $f_{p}=0$. This simple
method of the velocity reversals does not change the Boltzmann equation and
thus makes the density profile of the halo unchanged by varying $\lambda$
(Lynden-Bell, 1960; Weinberg, 1985; Long et al., 2014; Collier et al., 2018,
2019b). In our halo model, the maximum and minimum value of $\lambda$ is $\pm
0.16$, corresponding to either all prograde or all retrograde particles.

Table 1 lists the names, bulge-to-disk mass ratio, spin parameter $\lambda$,
and prograde fraction $f_{p}$ of our models. The models in the C00 series do
not possess a bulge, while those in the C10 series have a bulge whose mass is
10% of the disk mass. The infixes P and R stand for a prograde and retrograde
halo, respectively, and the numbers after the infixes denote $\lambda$
multiplied by 100. Note that the models C00P00 and C10P00 have a non-spinning
halo, and are identical to models C00 and C10 presented in JK23, respectively.
In the present work, we do not consider the rotation of a bulge: the effect of
bulge rotation was recently studied by Li et al. (2024).

Figure 2: Radial distributions of the halo and disk rotational velocity
$v_{\mathrm{rot}}$ for C00P16, C00P00, and C00R16.

Figure 2 compares the rotational velocities $v_{\mathrm{rot}}=\langle
v_{\phi}\rangle$ for the disk and halo in models C00P16, C00P00, and C00R16.
Here, $v_{\phi}$ is the azimuthal velocities of the particles and the angle
brackets denote the azimuthal and vertical averages. Note that
$v_{\mathrm{rot}}$ is slightly smaller than $v_{c}$ shown in Figure 1,
corresponding to asymmetric drift (e.g., Oh et al., 2008). In model C00P00,
the rotational velocity of the halo is nearly zero. However, when all halo
particles move in the same azimuthal direction, the rotational velocity of the
halo particles is almost flat, amounting to $v_{\mathrm{rot}}\sim\pm 120\;{\rm
km}\;{\rm s}^{-1}$ at $R\gtrsim 1\;{\rm kpc}$. We will show in Section 3 that
the halo rotation and velocity shear results in the creation of spirals,
evolving into a weak bar, via swing amplification of density perturbations
even in a halo-only system with $\lambda=\pm 0.16$.

To construct our galaxy models, we utilize the GALIC code (Yurin & Springel,
2014) which finds a desired equilibrium by adjusting the velocities of
individual particles to the level that satisfies the collisionless Boltzmann
equations. The number of particles used is $N_{d}=1.0\times 10^{6}$,
$N_{b}=1.0\times 10^{5}$, and $N_{h}=2.6\times 10^{7}$ for the disk, bulge,
and halo, respectively. The mass of a single particle is set to $\mu=5\times
10^{4}\;{\rm M}_{\odot}$, which is equal for all three components.

We evolve our galaxy models until $t=10\;{\rm Gyr}$ using a public version of
the Gadget-4 code (Springel et al., 2021). We take the multipole expansion of
order $p=4$ for fast force evaluation and a hierarchical time-integration
scheme to reduce the computation time. We fix the force accuracy parameter to
$\alpha=3\times 10^{-4}$ and the softening parameter to $0.05\;{\rm kpc}$,
$0.01\;{\rm kpc}$, and $0.01\;{\rm kpc}$ for the halo, disk, and bulge
particles, respectively, which we check conserves the total angular momentum
within $\sim 0.1\%$.

Figure 3: Evolution of the Fourier coefficients with dimensionless wavenumber
$p$, defined in Equation 5, of the $m=2$ logarithmic arms in the regions with
$R\leq 5\;{\rm kpc}$ and $|z|\leq 5\;{\rm kpc}$ of the halo-only model with
$\lambda=0.16$. The dominant mode at $t=0.9,1.0,1.1\;{\rm Gyr}$ has
$p=-0.6,0.6,0.8$, respectively.  Figure 4: Projected density distributions of
the $m=2$ mode in the regions with $|z|\leq 5\;{\rm kpc}$ and $R\leq 5\;{\rm
kpc}$ for the halo-only model with $\lambda=0.16$ at (_a_) $t=0.9\;{\rm Gyr}$,
(_b_) $t=1.0\;{\rm Gyr}$, (_c_) $t=1.1\;{\rm Gyr}$, and (_d_) $t=3.0\;{\rm
Gyr}$. In each panel, the concentric rings indicate the radii from 1 to
$5\;{\rm kpc}$. The dashed lines draw the crest of the dominant modes
determined by Equation 5. Since the halo is rotating in the counterclockwise
direction, the dominant waves in (_a_) are leading, while they are trailing in
(_b_) and (_c_). The dominant mode is a bar in (_d_).

##  3 Halo-only Models

A rotating halo may be prone to forming weak, non-asymmetric structures via
swing amplification, just as in a rotating disk. To explore this possibility,
we run $N$-body simulations of halo-only models with $\lambda=0.0$ (no spin),
$0.16$ (prograde), and $-0.16$ (retrograde), constructed in Section 2. By
evolving these isolated halo models up to $t=2\;{\rm Gyr}$, we confirm that
the azimuthally-averaged, radial distributions of halo density, velocities,
and velocity dispersions remain unchanged, indicating that the system is
overall in dynamical equilibrium. At the same time, we also find that the
rotating halo allows small-amplitude, non-axisymmetric perturbations to grow
in regions with $R\leq 5\;{\rm kpc}$ and $|z|\leq 5\;{\rm kpc}$.

To quantify the strength of $m=2$ perturbations that grow in the halo, we
define the Fourier coefficients in $\ln R$ and $\phi$ as

$A(m,p)=\frac{1}{N}\sum_{j=1}^{N}\exp[i(m\phi_{j}+p\ln R_{j})],$ (5)

where $N$ is the total number of particles and $(R_{j},\phi_{j})$ are the
radial and azimuthal coordinates of the $j$-th particle in the regions with
$R\leq 5\;{\rm kpc}$ and $|z|\leq 5\;{\rm kpc}$. In Equation 5, the
dimensionless radial wavenumber $p$ is related to the pitch angle
$i=\tan^{-1}(m/p)$ of $m$-armed logarithmic spirals, and $p=0$ corresponds to
a bar mode (e.g., Sellwood & Carlberg 1984; Sellwood & Athanassoula 1986; Oh
et al. 2008).

Figure 3 plots the temporal evolution of $|A(2,p)|$ of the $m=2$ logarithmic
spiral waves in the halo-only model with $\lambda=0.16$ at $t=0.9$–$1.1\;{\rm
Gyr}$. The perturbations are dominated by the mode with $p=-0.6$ at
$t=0.9\;{\rm Gyr}$ and $p=0.8$ at $t=1.1\;{\rm Gyr}$. The corresponding
projected distributions of the $m=2$ spiral waves, and the loci of the
dominant modes are presented in Figure 4(_a_)–(_c_). The waves grow as they
swing from leading (negative $p$) to trailing (positive $p$) configurations.

Figure 5 plots the temporal changes in $|A(2,0)|$, the strength of a bar mode
with $p=0$, in the halo-only models with $\lambda=0.16$ (red), $\lambda=0$
(black), and $\lambda=-0.16$ (blue), as well as its position angle

$\psi(R)\equiv\frac{1}{2}\tan^{-1}\left[\frac{\sum_{j}\sin(2\phi_{j})}{\sum_{j}\cos(2\phi_{j})}\right],$
(6)

measured at $R=2\;{\rm kpc}$ for the $\lambda=0.16$ model. When the halo does
not spin, the perturbations remain very weak with $|A(2,0)|\lesssim 5\times
10^{-3}$, primarily due to Poisson noise in the particle distribution.222By
running models with differing $N_{h}$, we have confirmed that $|A(2,0)|\propto
N_{h}^{-1/2}$ in the models with $\lambda=0$. However, density perturbations
in the rotating halos are swing-amplified to achieve $|A(2,0)|\sim 0.04$ at
$t\sim 2.4$–$2.8\;{\rm Gyr}$ in models with $\lambda=\pm 0.16$. The small
difference in $|A(2,0)|$ between the prograde and retrograde halos is
presumably caused by the amplitudes of the initial leading perturbations most
susceptible to swing amplification.333Since the initial particle distribution
is random, the power spectrum of the corresponding density perturbations is
not symmetric in the radial wavenumber $p$. The leading waves with negative
$p$ in the prograde halo become trailing waves with positive $p$ in the
retrograde halo.

Figure 5:  (_a_) Temporal changes in the strength of the bar mode in the
regions with $R\leq 5\;{\rm kpc}$ and $|z|\leq 3\;{\rm kpc}$ of the halo-only
models with $\lambda=0.16$ (red), $\lambda=0.0$ (black), and $\lambda=-0.16$
(blue). (_b_) Temporal evolution of the bar position angle $\psi$ measured at
$R=2\;{\rm kpc}$ in the $\lambda=0.16$ model.

Figure 4(_d_) shows that the spirals eventually grow into a weak bar. By using
the cross-correlation of the halo surface densities at two different epochs in
the annular regions with width $\Delta R=0.1\;{\rm kpc}$ at $R=1\;{\rm kpc}$
for $t=2$–$3\;{\rm Gyr}$ (see Seo et al. 2019, JK23), we find that the bar in
the halo-only model with $\lambda=0.16$ has a pattern speed of $\Omega_{h}\sim
6\;{\rm km}\;{\rm s}^{-1}\;{\rm kpc}^{-1}$, rotating in the counterclockwise
direction. Note that $\Omega_{h}\approx d\psi/dt$ for $t=2$–$3\;{\rm Gyr}$
shown in Figure 5(_b_). Similarly, the perturbed density in the retrograde
halo with $\lambda=-0.16$ grows first into trailing spirals and then to a bar
with a pattern speed of $\Omega_{h}\sim-6\;{\rm km}\;{\rm s}^{-1}\;{\rm
kpc}^{-1}$. These indicate that a spinning halo _alone_ can form a weak bar
that rotates in the same direction as the halo.

Figure 6: Snapshots of the disk surface density at every $0.1\;{\rm Gyr}$
interval from $t=0.1$ to $1.8\;{\rm Gyr}$ in models C00P00 with non-spinning
halo (top), C00P16 with a prograde halo (middle), and C00R16 with a retrograde
halo (bottom).

Since swing amplification requires self-gravity, rotational shear, and
epicyclic shaking (Toomre, 1981), it has been studied extensively in the
context of galactic disks. However, it can also operate in a spherical system
as long as the three agents mentioned above are present. As Figure 2 shows,
our spinning halo with $\lambda=\pm 0.16$ follows an almost flat rotation
curve, indicative of shear and epicyclic motions. By analyzing the orbits of
all the halo particles, we find that the faction of the particles whose orbits
are relatively circular with $\max(v_{\phi})/\max(v_{R})\geq 2$ and remain in
the inner regions with $R\leq 5\;{\rm kpc}$ close to the midplane with
$|z|\leq 5\;{\rm kpc}$ is about 0.15% of the halo mass, corresponding to $\sim
10\%$ of the disk mass at $R\leq 5\;{\rm kpc}$. This suggests that swing
amplification in the spinning halo can still allow perturbations to grow to
produce $m=2$ spirals and a weak bar, albeit much weaker than a strong bar
that forms in the disk.

To summarize, our rotating halo in itself is vulnerable to the formation of
non-axisymmetric structures. Small-amplitude (internal) perturbations in a
spinning halo grow due to swing amplification into spirals and then to a bar
which rotates in the same sense as the halo. As we will show in Section 4,
this tendency of a spinning halo forming non-axisymmetric structures even
without external perturbations promotes or delays the bar formation in a
stellar disk depending on its direction of rotation relative to the halo.

Figure 7: Temporal evolution of $A_{2}/A_{0}$ (solid) and $A_{3}/A_{0}$
(dashed) for models C00P00 (black), C00P16 (red) and C00R16 (blue) at early
time ($t\leq 2\;{\rm Gyr}$).

##  4 Full Galaxy Models

We now present the simulation results of full galaxy models with halo and disk
listed in Table 1. We first focus on the bar formation and evolution in our
models, and then explore the bar pattern speed, angular momentum transport,
and vertical buckling instability.

###  4.1 Bar Formation and Evolution

All of our models produce bars in the disks. The bar formation is due to
successive swing amplification aided by feedback loops (e.g., Binney &
Tremaine, 2008): small-amplitude perturbations in a stellar disk grow as they
transition from a leading to a trailing configuration, and trailing waves
appear as leading waves after they pass through the galaxy center. In the
growing stage, they interact gravitationally with the particles in the halo
and bulge, and the halo spin certainly affects the wave growth. When the waves
grow sufficiently, they shape into a bar. Figure 6 plots snapshots of the disk
surface density for models at every $0.1\;{\rm Gyr}$ interval from $t=0.1$ to
$1.8\;{\rm Gyr}$ in models C00P00 (top), C00P16 (middle), and C00R16 (bottom).

We follow the method of Kataria & Das (2018) to measure the strength of a bar
that forms in the disk (see also JK23). We consider an annulus centered at
radius $R$ with width $\Delta R=1\;{\rm kpc}$ in the regions with $R\leq
10\;{\rm kpc}$ of the disk, and calculate the amplitudes of the Fourier modes
as

$\displaystyle a_{m}(R)$ $\displaystyle=\sum_{j}\mu_{j}\cos{(m\phi_{j})},$
(7a) $\displaystyle b_{m}(R)$
$\displaystyle=\sum_{j}\mu_{j}\sin{(m\phi_{j})},$ (7b)

where $\phi_{j}$ and $\mu_{j}$ are the azimuthal angle and mass of the $j$-th
disk particle in the annulus, respectively. Then we define the strength of the
density perturbations with mode $m$ as the maximum wave amplitude across all
annuli:

$\frac{A_{m}}{A_{0}}=\max_{R}\left\\{\frac{[a_{m}(R)^{2}+b_{m}(R)^{2}]^{1/2}}{\sum_{j}\mu_{j}}\right\\}.$
(8)

Note that the term inside the curly brackets in Equation 8 is equal to
$|A(m,0)|$ defined in Equation 5 if the particles have the same mass.

Figure 7 plots the temporal evolution of $A_{2}/A_{0}$ and $A_{3}/A_{0}$ at
early time ($t\leq 2\;{\rm Gyr}$) for models C00P00 (black), C00P16 (red), and
C00R16 (blue). In model C00P00 with no halo spin, the disk favors the growth
of $m=2$ and 3 spirals, with the $m=3$ mode dominating at $t\lesssim 0.8\;{\rm
Gyr}$ (see the top row of Figure 6). These modes with different $m$ interact
nonlinearly with each other, and one arm of the $m$ = 3 spirals merges with
the other two arms, producing a bar at $t\sim 0.9\;{\rm Gyr}$ (see also, Seo
et al., 2019). The $m=3$ mode in these models decays secularly to become
$A_{3}/A_{0}\sim 3\times 10^{-2}$ at the end of the runs ($t=10\;{\rm Gyr}$).

In model C00P16 with $\lambda=0.16$, the growth of the $m=2$ mode in the disk
is assisted by the prograde halo which has an intrinsic tendency of amplifying
the $m=2$ perturbations rotating in the same sense as the halo. Consequently,
the $m=2$ spirals in this model become stronger and thus transform to a bar
earlier ($t\sim 0.5\;{\rm Gyr}$) than in model C00P00 (see the middle row of
Figure 6). We note that the regions outside the bar in model C00P16 are more
strongly perturbed by the $m=2$ spirals than in model C00P00, limiting further
growth of the bar. In model C00R16 with $\lambda=-0.16$, in contrast, the
growth of the $m=2$ mode in the disk is opposed to some extent by the counter-
rotating halo, which in turn allows for the $m=3$ mode to grow more strongly
than in model C00P00. Therefore, it takes longer for the $m=3$ spirals to
transform eventually into a bar than under the non-spinning halo.

Figure 8: Temporal evolution of the bar strength $A_{2}/A_{0}$ for the C00
series (top) and C10 series (bottom) for $t\leq 10\;{\rm Gyr}$. A bar tends to
form early under a prograde halo and without bulge.  Figure 9: Snapshots of
the disk surface density (top) and the projected disk density along the bar
semi-minor axis (bottom) at the end of the runs ($t=10\;{\rm Gyr}$) for all
models. The upper and lower rows correspond to the models in the C00 and C10
series, respectively. The $x$\- and $y$-axis correspond to the semi-major and
semi-minor axis, respectively. In the edge-on views, the dotted contours
denote $\int\rho\,dy=10^{3.5},10^{3.0},10^{2.5},10^{2.0}\;{\rm
M}_{\odot}\;{\rm pc}^{-2}$ from inside to outside.

Figure 8 plots the temporal evolution of $A_{2}/A_{0}$ for C00 series (upper
panel) and C10 series (lower panel) over entire evolution. The bar in model
C00P16 stops growing at $t\sim 1.0\;{\rm Gyr}$. At this time, the inner parts
of the spirals dominating the outer regions become in phase with the bar,
making $A_{2}/A_{0}$ increase temporarily to $\sim 0.4$. As the spirals rotate
relative to the bar and become out of phase, $A_{2}/A_{0}$ decreases to $\sim
0.32$ at $t\sim 2\;{\rm Gyr}$. Since the regions outside the bar are strongly
perturbed by the $m=2$ spirals, they do not participate in the bar growth,
making $A_{2}$ stay almost constant for a long period. The bar undergoes a
weak buckling instability at $t\sim 9.5\;{\rm Gyr}$ (see Section 4.3).

In model C00P00 with no halo spin, $A_{2}/A_{0}$ peaks near $t\sim{1.1}\;{\rm
Gyr}$ when the spirals are in phase with the bar, after which $A_{2}/A_{0}$
decreases as the bar becomes out of phase with the spirals. Since the regions
outside the bar possess weaker spirals and are thus less disturbed than in
model C00P16, the bar in this model can grow in size and strength by gathering
outer particles until it undergoes a buckling instability at $t\sim 5.2\;{\rm
Gyr}$. In model C00R16 with a counter-rotating halo, the bar grows slower but
more strongly than the other models since the spiral perturbations in the
outer regions are weakest. The bar becomes so strong that the buckling
instability occurs at $t\sim 2.5\;{\rm Gyr}$, earlier than the other models.

The models in the C10 series with a classical bulge also form a bar, but later
than the no-bulge counterparts. This is because the gravity of a bulge reduces
the feedback loop and swing amplification. The bars in the C10 series are
weaker and shorter at early time ($t\lesssim 2\;{\rm Gyr}$) than those in the
C00 series, indicating that the presence of a non-rotating bulge, in general,
resists the bar formation (Sellwood 1980; Binney & Tremaine 2008; Kataria &
Das 2018; Saha & Elmegreen 2018; JK23).

The bars in models C10P16, C10P10 and C10P06 are always weaker than those in
the no-bulge counterparts. However, the bars in the other models with a bulge
do not undergo buckling instability and can thus be stronger than those
without a bulge. For example, the bar in model C10P00 keeps growing to reach
$A_{2}/A_{0}\sim 0.64$, while the bar in model C00P00 suddenly becomes weaker
after experiencing buckling instability at $t\sim 5.2\;{\rm Gyr}$.

Figure 9 plots the snapshots of the disks projected along the vertical
direction and the bar semi-minor axis at $t=10\;{\rm Gyr}$ for all models. The
bar rotating in the counterclockwise direction is oriented so that the $x$\-
and $y$-axes correspond to the semi-major and semi-minor axes, respectively.
The bar in model C00P16 shows a twist of isodensity contours in the outer
parts. The bar in model C10P16 is so weak that it can be regarded as an oval.
The surface density near the center has a dumbbell-like distribution only in
model C00P00. The models with relatively strong bars (other than models
C00P16, C00P10, C10P16, C10P10 and C10P06) possess a boxy peanut-shaped (BPS)
bulge as well as an inner ring just outside the bar.

Figure 10: Temporal changes in the angular momentum $\Delta L\equiv L(t)-L(0)$
of the entire disk (top), halo (middle), and bulge (bottom) for the models in
the C10 series. The ordinate is normalized by the initial angular momentum of
the entire disk.

###  4.2 Angular Momentum and Bar Pattern Speed

After formation, a bar embedded in a disk transfers its angular momentum to a
halo (and a bulge). To study the radial and temporal dependence of the angular
momentum transfer, we bin each component radially into concentric cylindrical
shells with width $\Delta R=1\;{\rm kpc}$ and calculate the vertical component
of the angular momentum

$L=\sum_{i}\mu_{i}(x_{i}v_{y,i}-y_{i}v_{x,i})$ (9)

in each shell. Figure 10 plots $\Delta L\equiv L(t)-L(0)$ summed over the
entire shells as a function of time for the models in the C10 series. Overall,
the amount of angular momentum absorbed by the bulge is much less than that by
the halo. The disk under a counter-rotating halo loses more angular momentum
than under the prograde halo, making the bar in the former grow stronger and
longer at the end of the runs (see Figure 9). Of course, the actual temporal
change of $\Delta L$ depends on not only $\lambda$ but also the bar properties
such as strength, length, pattern speed, mass distribution, etc.

Figure 11 displays the rate of angular momentum transfer, $dL/dt$, as a
function of $R$ and $t$ in the halo and disk for models C00P16, C00P00, and
C00R16. When $\lambda>0$, the halo particles inside the corotation resonance
$R_{\textrm{CR}}$ with the bar rotate faster than, and thus lose their angular
momenta to the bar, while the particles in the outer halo gain angular
momentum (Long et al., 2014; Collier et al., 2018). Since $R_{\textrm{CR}}$ is
larger for larger $\lambda$, the angular momentum transfer from disk to halo
tends to be larger for smaller $\lambda>0$. In model C00P16, most angular
momentum transfer occurs near $t\sim 1\;{\rm Gyr}$ when the bar is strongest.
The bar angular momentum is also transferred to the outer disk at
$R>R_{\textrm{CR}}$.

When $\lambda\leq 0$, in contrast, all the halo particles lag behind the bar
and thus allow the efficient angular momentum transfer. Figure 11 shows that
the halo in model C00P00 absorbs angular momentum mainly through three
resonances which extend radially outward as the bar grows in size over time
(e.g., Villa-Vargas et al. 2009; Collier et al. 2018, 2019a, 2019b). Again,
the largest $dL/dt$ occurs near $t\sim 3$–$5\;{\rm Gyr}$ when the bar is
strongest. The outer disk absorbs the angular momentum emitted by the bar. The
halo in model C00R16 absorbs angular momentum strongly at $t\sim 1$–$3\;{\rm
Gyr}$ and $\sim 5$–$7\;{\rm Gyr}$, corresponding to a strong bar before the
buckling instability and in the regrowth stage (see Figure 8).

Figure 11:  Rates of the angular momentum change, $dL/dt$, as a function of
cylindrical radius and time for models C00P16 (left), C00P00 (middle), and
C00R16 (right). The cylindrical shells have a width $\Delta R=1\;{\rm kpc}$
and height $|z|=\infty$ for the halos (upper panels) and the stellar disks
(lower panels).

We calculate the bar pattern speed, $\Omega_{b}$, using the cross-correlation
of the disk surface density in the annular regions with width $\Delta
R=0.1\;{\rm kpc}$ at $R=2\;{\rm kpc}$: we check that the resulting
$\Omega_{b}$ agrees within $\sim 1\%$ with that from the method that utilizes
the temporal changes in the bar position angle
$\psi=0.5\tan^{-1}(b_{2}/a_{2})$ (JK23). Figure 12 plots the temporal changes
in the bar pattern speed for all models. As the bar loses its angular momentum
to the halo and bulge, its pattern speed decreases over time. The initial bar
pattern speed tends to be higher in models with a prograde halo and with a
bulge. This is consistent with the numerical results that a short bar tends to
rotate fast (Aguerri et al. 2015; Roshan et al. 2021; Lee et al. 2022; Frankel
et al. 2022; JK23). Note that the bar slows down very slowly in models with
$\lambda=0.16$ since the angular momentum absorption by the halo outside
$R_{\mathrm{CR}}$ is almost balanced by the angular momentum emission by the
halo inside $R_{\mathrm{CR}}$. Models with smaller $\lambda\;(>0)$ have
smaller $R_{\mathrm{CR}}$ and thus causes $\Omega_{b}$ to decay faster over
time. This is qualitatively consistent with the analytic result of Chiba &
Kataria (2024) who showed that the dynamical friction of a bar due to halo
particles weakens with increasing halo spin. Models with $\lambda<0$ have
$\Omega_{b}$ decreasing more or less similarly to that in the $\lambda=0$
models, largely consistent with the result of Collier et al. (2019a) in that
the amount of angular momentum transport is insensitive to $\lambda$ as long
as the halo has a retrograde spin.

###  4.3 Buckling Instability

All the bars in our models thicken vertically over time, evolving into a BPS
bulge if strong. The BPS strength is usually defined as

$P_{s}=\max_{R}\left(\frac{\widetilde{|z|}}{\widetilde{|z_{0}|}}\right),$ (10)

where the tilde denotes the median and $z_{0}$ is the initial height for the
disk particles at $R\leq 10\;{\rm kpc}$ (Iannuzzi & Athanassoula, 2015;
Fragkoudi et al., 2017; Seo et al., 2019). Figure 13 plots temporal evolution
of $P_{s}$ for all models. Note that the models in the C00 series with no
bulge (except for model C00P16) go through a short period where $P_{s}$
increases rapidly due to buckling instability. In contrast, $P_{s}$ in the
models in the C10 series with a bulge does not experience such a rapid
increase. Instead, it increases relatively slowly due to the vertical heating
of disk particles via gravitational interactions with the bar without
undergoing buckling instability.

Figure 14 compares the projected disk densities along the bar semi-minor axis
in models C00P00 (left) and C10P00 (right) for $t=4.5$–$7.0\;{\rm Gyr}$.
Clearly, the bars in both models produce a BPS bulge. Note that the projected
disk densities in model C00P00 become asymmetric with respect to the $z=0$
plane at $t\sim 5.5\;{\rm Gyr}$, indicative of the buckling instability.
However, the bar in model C10P00 almost always maintains a mirror symmetry
about the midplane, while it gradually thickens vertically. This relatively
slow thickening of the bars can be caused by 2:1 and 4:1 vertical resonances
(Combes et al., 1990; Pfenniger & Friedli, 1991; Patsis et al., 2002; Sellwood
& Gerhard, 2020). This indicates that the buckling instability is not a
necessary but sufficient condition for forming a BPS bulge (Li et al., 2023).

Figure 12: Temporal changes of the bar pattern speed $\Omega_{b}$ for the C00
series (top) and C10 series (bottom).

Figure 13: Time evolution of the BPS strength, $P_{s}$, defined in Equation 10
for the models in the C00 series with no bulge (top) and C10 series with a
bulge (bottom). Figure 14: Contours of the logarithm of the projected disk
densities at $t=4.5$–$7.0\;{\rm Gyr}$ for models C00P00 (left) and C10P00
(right). The buckling instability occurs at $t\sim 5.5\;{\rm Gyr}$ in model
C00P00. The $x$\- and $z$-axis correspond to the bar semi-major and vertical
axes, respectively. The dotted contours denote
$\int\rho\,dy=10^{3.5},10^{3.0},10^{2.5},10^{2.0}\;{\rm M}_{\odot}\;{\rm
pc}^{-2}$ from inside to outside.  Figure 15: Temporal evolution of the mean
vertical velocity $|\langle v_{z}\rangle|$ at $R=2.0\;{\rm kpc}$ of the disk
for the models in the C00 series (top) and the C10 series (bottom).  Figure
16: Temporal evolution of the vertical asymmetry parameter
$A_{\mathrm{asym}}$, defined in Equation 11, from all the disk particles
(left) and the particles at $R<5\;{\rm kpc}$ (right) for the models in the C00
series (top) and the C10 series (bottom).

Since the buckling instability deforms the disk in the vertical direction, it
naturally involves the asymmetry in the vertical velocities of disk particles.
One way to quantify the onset and strength of the buckling instability is to
monitor the evolution of $|\langle v_{z}\rangle|$, where the angle brackets
denote the spatial average at $R=2\;{\rm kpc}$ (Kwak et al., 2017, 2019).
Figure 15 plots $|\langle v_{z}\rangle|$ against time for all models. A sharp
peak in $|\langle v_{z}\rangle|$ corresponds to the buckling instability,
which starts to occur at $t=2.5,2.9,3.7,5.2,6.2,6.8\;{\rm Gyr}$ for models
C00R16, C00R10, C00R06, C00P00, C00P06, C00P10, respectively. Model C00P16
also undergoes a buckling instability at $t\sim 9.5\;{\rm Gyr}$, although it
is very mild. The bottom panel of Figure 15 shows that $|\langle
v_{z}\rangle|$ of the models in the C10 series exhibit only noises, suggesting
that they do not experience buckling instability at all.

Another way to measure the strength of the buckling instability is to use the
vertical asymmetry parameter of the disk introduced by Smirnov & Sotnikova
(2018) as

$A_{\textrm{asym}}=\left|\frac{A_{2}(z>0)-A_{2}(z<0)}{A_{0}}\right|,$ (11)

where $A_{2}(z>0)$ and $A_{2}(z<0)$ denote the Fourier amplitudes of the $m=2$
mode applied to the disk particles above and below the $z=0$ plane,
respectively (see also Li et al., 2024). Figure 16 plots $A_{\textrm{asym}}$
from all the disk particles (left) as well as the particles located at $R\leq
5\;{\rm kpc}$ (right). The models with no bulge have $A_{\textrm{asym}}$
varying a lot over time, while $A_{\textrm{asym}}$ for the models with a bulge
is close to zero. For the models without a bulge, $A_{\textrm{asym}}$ from the
entire disk exhibits longer-term variations than that from the regions with
$R\leq 5\;{\rm kpc}$, suggesting that the buckling instability initiates from
smaller radii and propagates radially outward. Compared to $|\langle
v_{z}\rangle|$, $A_{\textrm{asym}}$ for the models in the C00 series are more
widely distributed over time. This is because $A_{\textrm{asym}}$ captures the
asymmetry in a wider radial range and the buckling instability occurs in the
inner disk first and then propagates outward (Martinez-Valpuesta et al.,
2006). Note that $|\langle v_{z}\rangle|$ can be small if the nodes of the
vertical bending in the course of buckling are located close to $R=2\;{\rm
kpc}$, as in models C00P00, C00P06, and C00P10. This suggests that the peak
values of $A_{\textrm{asym}}$ represent the overall strength of the buckling
instability better than $|\langle v_{z}\rangle|$.

As Figures 15 and 16 show, the buckling instability occurs earlier in our
models with smaller $\lambda$, that is, as the fraction of retrograde halo
particles increases. This trend is not observed in previous investigations
(Long et al., 2014; Collier et al., 2018, 2019a, 2019b; Kataria & Shen, 2022).
In Collier et al. (2018, 2019a), the maximal pre-buckling amplitude of
$A_{2}/A_{0}$ is almost insensitive to $\lambda>0$, while it becomes smaller
as $\lambda$ decreases in the models with a retrograde halo. In our models,
the pre-buckling strength is most significant in the models with $\lambda=0$
and decreases as $|\lambda|$ increases in models with either a prograde or
retrograde halo (see Figure 8).

It is well known that the buckling instability occurs when the radial velocity
dispersion $\sigma_{R}$ of disk particles exceeds a critical value relative to
the vertical velocity dispersion $\sigma_{z}$ (e.g., Binney & Tremaine, 2008).
Toomre (1966) and Araki (1987) found that the critical value is at
$\sigma_{z}/\sigma_{R}\sim 0.3$ for non-rotating, razor-thin disks. For
realistic disks with spatially-varying $\sigma_{z}/\sigma_{R}$, Raha et al.
(1991) found that the disks undergo buckling instability when
$\sigma_{z}/\sigma_{R}\lesssim 0.25$–$0.55$. For barred galaxies, $N$-body
simulations of Martinez-Valpuesta et al. (2006) and Kwak et al. (2017) found
that buckling occurs when $\sigma_{z}/\sigma_{R}\lesssim 0.6$.

Figure 17 plots the temporal changes in $\sigma_{z}/\sigma_{R}$ of the disk
particles at $R=2\;{\rm kpc}$ for the models in the C00 series (top) and the
C10 series (bottom). For most of the models without a bulge, the onset of the
buckling instability corresponds to $\sigma_{z}/\sigma_{R}\simeq 0.47$–$0.60$.
Note that model C00P16 has $\sigma_{z}/\sigma_{R}\sim 0.75$ at $t=9.5\;{\rm
Myr}$, but it still exhibits mild buckling instability, as evidenced by the
asymmetry parameter shown in Figure 16. All the models with a bulge have
$\sigma_{z}/\sigma_{R}\gtrsim 0.55$, and remain stable against buckling
instability. Instead, the bars in these models thicken gradually to produce
BPS bulges (see Figure 14). This indicates that the critical value of
$\sigma_{z}/\sigma_{R}$ for buckling instability is applicable only to the
models without a bulge.

Figure 17: Temporal variations of the radio $\sigma_{z}/\sigma_{R}$ of the
vertical to radial velocity dispersion of the disk particles at $R=2\;{\rm
kpc}$ for the models in the C00 series (top) and the C10 series (bottom).

##  5 Discussion

We have investigated the formation and evolution of bars in disk galaxies
similar to the Milky Way, with a spinning halo. In this section, we discuss
our results compared to the previous investigations. We also use the bar
properties to constrain the models for the Milky Way.

###  5.1 Bar Formation under a Spinning Halo

In our models, a disk embedded in a faster corotating halo develops a bar
earlier, while the bar formation is delayed under a counter-rotating halo,
consistent with the results of the previous studies (Saha & Naab, 2013; Long
et al., 2014; Collier et al., 2018, 2019a). As mentioned in Section 1, Saha &
Naab (2013) tried to explain this in terms of the Ostriker-Peebles parameter
$t_{\rm OP}$, arguing that the halo spin increases $t_{\rm OP}$ and thus makes
the disk more susceptible to bar formation. Since $t_{\rm OP}$ does not
consider the sense of halo rotation, however, it cannot explain why a bar
forms faster in a disk under a prograde halo than under the retrograde
counterpart.

To account for the disk response to a spinning halo, we examine the growth of
various modes in an isolated halo and a disk-halo system. In Section 3, we
present that our spinning halo alone is vulnerable to the growth of $m=2$
spirals due to swing amplification, transforming to a weak bar rotating in the
same sense as the halo. When a disk is inserted into the halo, $m=2$ spirals
in the disk growing also via swing amplification supply gravitational
perturbations for the halo. If the halo spin is prograde, the perturbations
are well maintained and amplified in the halo, giving positive feedback to the
spirals in the disk. Therefore, the disk and the prograde halo work
constructively to make the spirals grow faster than in the case with no halo
spin, promoting the bar formation (see Figure 7). If the halo spin is instead
retrograde, $m=2$ perturbations in the disk and halo become out of phase and
interact deconstructively. In this case, the $m=2$ spirals in the disk grow
more slowly, delaying the bar formation.

###  5.2 Buckling Instability

Numerical simulations commonly found that bars undergo vertical buckling
instability when sufficiently strong, resulting in weakened and shorted bars
and eventually leading to BPS bulges. Our simulations also show that strong
bars are subject to buckling instability. There are, however, some
discrepancies between the results of our models and the previous simulations
(Saha & Naab, 2013; Long et al., 2014; Collier et al., 2018, 2019b). First, in
models of Collier et al. (2018, 2019a, 2019b), there is no apparent
correlation between $\lambda$ and the epoch of the buckling instability, while
the buckling in our models is progressively delayed with increasing $\lambda$.
Second, the bar strength before the buckling instability is almost independent
of $\lambda$ in models with a prograde halo (Collier et al., 2018, 2019b) and
decreases with $\lambda$ in models with a retrograde halo (Collier et al.,
2019a), while it is a decreasing function of $|\lambda|$ for the whole range
of $\lambda$ in our models (see Figure 8). Third, the buckling instability
almost destroys the bars when $\lambda\gtrsim 0.06$ in Collier et al. (2018,
2019b), while the bars, albeit becoming weaker, still remain strong after the
buckling in our models. Lastly, our models with a bulge do not undergo
buckling instability at all, while the models even with a bulge in Li et al.
(2024) experience buckling instability, although the onset and amplitude of
the buckling decrease with increasing bulge mass.

It is uncertain what causes these differences in the numerical results, but
the most likely reason may be differences in the galaxy models adopted. Our
models consider a Milky Way-sized galaxy, while the previous authors (Long et
al., 2014; Collier et al., 2018, 2019a, 2019b; Collier & Madigan, 2021; Li et
al., 2024) employed the models with a twice less massive halo and a twice
thicker disk than our models. Also, their halo modeled by a truncated NFW
profile (Navarro et al., 1996) has a different degree of central concentration
and a quite different rotation curve from ours (e.g., Figure 1 in the present
work can be compared with Figure 1 of Lieb et al. 2022). These suggest that
one should be cautious about drawing a general conclusion by running a limited
set of numerical simulations.

The firehose instability has been invoked as the physical mechanism behind the
buckling instability (Toomre, 1966; Binney & Tremaine, 2008). An alternative
mechanism may be the trapping of overlapped planar and vertical 2:1
resonances, proposed recently by Li et al. (2023). In the firehose
instability, bending perturbations become unstable if the centrifugal force on
stars traveling over corrugations exceeds the restoring gravitational force
from the other stars, which is achieved when $\sigma_{z}/\sigma_{R}$ becomes
less than a critical value (see e.g., Combes & Sanders, 1981; Combes et al.,
1990; Raha et al., 1991; Merritt & Sellwood, 1994; Martinez-Valpuesta et al.,
2006; Kwak et al., 2017, 2019; Li et al., 2024, and references therein). Our
simulations show that the bucking instability in strongly barred galaxies
occurs when $\sigma_{z}/\sigma_{R}\lesssim 0.47$–$0.60$, roughly consistent
with the results of Martinez-Valpuesta et al. (2006) and Kwak et al. (2017).

The buckling instability in most of our models is consistent with the
conventional picture that it needs $\sigma_{z}/\sigma_{R}<0.6$ for operation
and that it increases $\sigma_{z}/\sigma_{R}$ rapidly, while decreasing
$A_{2}/A_{0}$, as Figures 8 and 17 show. However, the buckling instability in
model C00P16 is exceptional in that it works even with
$\sigma_{z}/\sigma_{R}\sim 0.75$ and does not cause a sharp increase in
$\sigma_{z}/\sigma_{R}$ and bar weakening. Figure 18 compares the face-on
snapshots of the mean vertical velocities, $\bar{v}_{z}$, and the mean
vertical positions, $\bar{z}$, of the disk particles for models C00R16 (left)
and C00P16 (right) each at the epoch of the maximum buckling. Model C00R16
which undergoes strong buckling at $t\sim 2.7\;{\rm Gyr}$ exhibits
characteristic quadrupole patterns in the face-on maps, with
$|\bar{v}_{z}|\lesssim 30\;{\rm km}\;{\rm s}^{-1}$ and $|\bar{z}|\lesssim
0.5\;{\rm kpc}$ (see, e.g., Łokas, 2019; Xiang et al., 2021; Li et al., 2023).
Note that model C00P16 also shows similar quadrupole patterns, although they
are weaker with $|\bar{v}_{z}|\lesssim 10\;{\rm km}\;{\rm s}^{-1}$ and
$|\bar{z}|\lesssim 0.2\;{\rm kpc}$, demonstrating that it indeed suffers the
buckling instability near $t\sim 9.7\;{\rm Gyr}$. One cannot rule out the
possibility of overlapping planar and vertical 2:1 resonances, as envisaged by
Li et al. (2023), for the nature of the bucking instability in model C00P16.

Figure 18: Face-on maps of the mean vertical velocities (top) and mean
vertical positions (bottom) of the stars at $t=2.7\;{\rm Gyr}$ for model
C00R16 (left) and at $t=9.7\;{\rm Gyr}$ for model C00P16 (right). The contour
in each panel outlines the bar boundaries defined as the regions with
$\Sigma=2\times 10^{2}\;{\rm M}_{\odot}\;{\rm pc}^{-2}$.

###  5.3 Milky Way Bar

We now discuss which of our models best describes the properties of the Milky
Way bar. Observations show that the Milky Way bar is long with size $R_{b}\sim
4.5$–$5\;{\rm kpc}$ and rotates slowly with pattern speed
$33<\Omega_{b}<45\;{\rm km}\;{\rm s}^{-1}\;{\rm kpc}^{-1}$ (Wegg et al., 2015;
Sormani et al., 2015; Portail et al., 2017; Bland-Hawthorn & Gerhard, 2016;
Clarke & Gerhard, 2022). Measuring the spin parameter $\lambda$ of a galaxy is
a difficult task observationally. Recently, Obreja et al. (2022) found a
correlation between the angular momenta of dark halos and disks from a sample
of galaxies identified in cosmological simulations, and used the relation to
estimate $\lambda_{\rm MW}$ of the Milky Way after measuring its disk angular
momentum. They obtained $\lambda_{\rm MW}\sim 0.061$ when the dark halo of the
Milky Way follows a contracted NFW profile, and $\lambda_{\rm MW}\sim 0.088$
for an uncontracted profile.

Figure 19: Temporal evolution of the bar length $R_{b}$ in model C10P06. The
blue line corresponds to $R_{b}$ calculated by imposing the changes of the
position angle $\Delta\psi\leq 0.1$ of the $m=2$ mode, while the orange line
is obtained by defining the bar as the regions with $\Sigma\geq 2\times
10^{2}\;{\rm M}_{\odot}\;{\rm pc}^{-2}$.

JK23 found that models C10 and C20 with $\lambda=0$, containing a classical
bulge with mass $10\%$ and $20\%$ of the disk mass, respectively, match the
properties of the Milky Way bar quite well. This suggests that model C10P06 in
the present paper best represents the Milky Way when the halo spin is
considered. Indeed, Figure 12 shows that the bar in model C10P06 has
$\Omega_{b}\sim 30$–$35\;{\rm km}\;{\rm s}^{-1}\;{\rm kpc}^{-1}$ for a long
period of time. Figure 19 plots the temporal changes of the bar length
calculated from two methods: (1) the maximum radius where the change in the
bar position angle $\Delta\psi<0.1$ and (2) the semi-major axis with
$\Sigma\geq 2\times 10^{2}\;{\rm M}_{\odot}\;{\rm pc}^{-2}$. The bar length of
$R_{b}\sim 4$–$5\;{\rm kpc}$ at $t\gtrsim 3\;{\rm Gyr}$ is consistent with the
observed size of the Milky Way bar. That $R_{b}$ from $\Delta\psi$ is longer
than that from $\Sigma$ at $t\lesssim 4\;{\rm Gyr}$ is due to the temporal
alignment of the inner parts of the spirals with the bar ends.

We remark on some caveats in concluding that the Milky Way bar is similar to
that in model C10P06. First, as Obreja et al. (2022) suggested, a contracted
NFW profile most likely describes the halo of the Milky Way, while our galaxy
models adopt the Hernquist profile. Second, our models consider only isolated
galaxies, while non-axisymmetric perturbations by external triggers such as
accretion and minor mergers with satellite galaxies can impact the bar
evolution in real situations (Zana et al., 2018; Grand et al., 2019; Peschken
& Łokas, 2019; Ghosh et al., 2021; Cavanagh et al., 2022). The current halo
spin of the Milky Way has likely been affected by the most recent major
merger. Finally, while we adopt gas-free galaxy models, the presence of gas in
the disk tends to reduce the bar slow-down rate (Seo et al., 2019; Beane et
al., 2023). To understand the bar evolution of the Milky Way more accurately,
it is desirable to run more realistic models by allowing for the halo density
distribution, accretion, minor mergers, and the presence of gas and associated
star formation.

##  6 Conclusions

We have presented the results of $N$-body simulations to investigate the
effects of halo spin and the presence of a classical bulge on the bar
formation and evolution. For this, we have borrowed two galaxy models from
JK23 that have similar properties to the Milky Way: model C00 without a bulge
and model C10 with a bulge whose mass is 10% of the disk mass. We have varied
the spin parameter $\lambda$ of the halo from $0.16$ (prograde) to $-0.16$
(retrograde) in both series and run the models up to $t=10\;{\rm Gyr}$. Our
main conclusions are summarized below.

  1. 1.

Simulations of the halo-only models show that our adopted halo is
gravitationally stable if it does not rotate ($\lambda=0$). When it rotates
with $\lambda=\pm\,0.16$, however, the halo itself is susceptible to forming
$m=2$ spirals via swing amplification, eventually growing into a weak bar. The
spirals and bar rotate in the same sense as the halo spin, with a pattern
speed of $\Omega_{h}\sim\pm\,7\;{\rm km}\;{\rm s}^{-1}\;{\rm kpc}^{-1}$. This
indicates that our rotating halo _alone_ is vulnerable to the formation of
non-axisymmetric structures rotating in the same direction as the halo.

  2. 2.

The tendency of our spinning halo to form $m=2$ non-axisymmetric structures
even without external perturbations affects the bar growth time in an embedded
disk. When $\lambda>0$, $m=2$ spirals growing via swing amplification in the
disk interact constructively with the $m=2$ mode in the prograde halo,
promoting the bar formation in the inner disk. When $\lambda$ is large, the
spirals in the outer regions are strong enough to limit the further bar
growth. When $\lambda<0$, in contrast, the $m=2$ mode in the disk interacts
destructively with the $m=2$ mode in the retrograde halo, delaying the bar
formation. Since the outer disk is relatively less perturbed by the spirals,
the bar in the disk under a counter-rotating halo can grow stronger and longer
than that in the models with large $\lambda$.

  3. 3.

A bar grows by losing its angular momentum to both a halo and bulge, although
the amount of angular momentum absorbed by the bulge is much less than that by
the halo. The halo particles inside (outside) the corotation resonance
$R_{\mathrm{CR}}$ with the bar can emit (absorb) angular momentum to (from)
the bar. Since $R_{\mathrm{CR}}$ is larger for larger $\lambda$, the angular
momentum transfer from disk to halo tends to be larger, causing the bar
pattern speed $\Omega_{b}$ to decay faster, for smaller $\lambda>0$. In models
with $\lambda=0.16$, the angular momentum absorption and emission are almost
balanced, keeping $\Omega_{b}$ and bar strength nearly constant for a long
period. Under a retrograde or non-spinning halo, in contrast, all the halo
particles lag behind the bar and thus can absorb angular momentum from the bar
efficiently. This causes $\Omega_{b}$ to decay relatively fast, almost
independent of $\lambda\leq 0$.

  4. 4.

All our models form a bar, and models with a strong bar develop a BPS bulge.
In our models, the presence of a classical bulge reduces the bar growth and
suppresses the buckling instability completely. In all models without a
classical bulge (C00 series), a BPS bulge is produced rapidly through the
buckling instability. In models with a classical bulge (C10 series), however,
bars slowly thicken vertically via gravitational interactions with the disk
particles.

  5. 5.

In models with no classical bulge, the buckling instability tends to occur
earlier along the sequence with smaller $\lambda$: a bar under a retrograde
halo deforms earlier in the vertical direction than the prograde counterpart.
For most of our models, the onset of the buckling instability corresponds to
the ratio of the velocity dispersion $\sigma_{z}/\sigma_{R}\lesssim
0.47$–$0.60$ at $R=2\;{\rm kpc}$, roughly consistent with the previous results
(e.g., Martinez-Valpuesta et al., 2006; Kwak et al., 2017), except model
C00P16 which undergoes mild buckling instability even at
$\sigma_{z}/\sigma_{R}\sim 0.75$. The models with a bulge have
$\sigma_{z}/\sigma_{R}\gtrsim 0.55$ and remain stable to buckling instability.

  6. 6.

Among our models, model C10P06 is similar to the Milky Way in terms of the
halo spin and bar properties. It has the halo spin parameter of
$\lambda=0.06$, identical to the observationally inferred value (Obreja et
al., 2022), and a classical bulge with mass 10% of the disk mass. The bar in
this model has a semi-major axis of $R_{b}\sim 4$–$5\;{\rm kpc}$ and a pattern
speed of $\Omega_{b}\sim 30$–$35\;{\rm km}\;{\rm s}^{-1}\;{\rm kpc}^{-1}$ for
an extended period of time, which are consistent with the observed properties
of the Milky Way bar.

## acknowledgments

The work of D.J. was supported by Basic Science Research Program through the
National Research Foundation of Korea (NRF) funded by the Ministry of
Education (RS-2023-00273275). The work of W.-T. K. was supported by a grant of
the National Research Foundation of Korea (2022R1A2C1004810). Computational
resources for this project were provided by the Supercomputing Center/Korea
Institute of Science and Technology Information with supercomputing resources
including technical support (KSC-2023-CRE-0175).

## References

  * Aguerri et al. (2009) Aguerri, J. A. L., Méndez-Abreu, J., & Corsini, E. M. 2009, A&A, 495, 491, doi: 10.1051/0004-6361:200810931
  * Aguerri et al. (2015) Aguerri, J. A. L., Méndez-Abreu, J., Falcón-Barroso, J., et al. 2015, A&A, 576, A102, doi: 10.1051/0004-6361/201423383
  * Ansar & Das (2023) Ansar, S., & Das, M. 2023, arXiv e-prints, arXiv:2311.11998, doi: 10.48550/arXiv.2311.11998
  * Araki (1987) Araki, S. 1987, AJ, 94, 99, doi: 10.1086/114451
  * Athanassoula (2002) Athanassoula, E. 2002, ApJ, 569, L83, doi: 10.1086/340784
  * Athanassoula (2005) —. 2005, MNRAS, 358, 1477, doi: 10.1111/j.1365-2966.2005.08872.x
  * Beane et al. (2023) Beane, A., Hernquist, L., D’Onghia, E., et al. 2023, ApJ, 953, 173, doi: 10.3847/1538-4357/ace2b9
  * Bett et al. (2007) Bett, P., Eke, V., Frenk, C. S., et al. 2007, MNRAS, 376, 215, doi: 10.1111/j.1365-2966.2007.11432.x
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529, doi: 10.1146/annurev-astro-081915-023441
  * Bullock et al. (2001) Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240, doi: 10.1086/321477
  * Buta et al. (2015) Buta, R. J., Sheth, K., Athanassoula, E., et al. 2015, ApJS, 217, 32, doi: 10.1088/0067-0049/217/2/32
  * Cavanagh et al. (2022) Cavanagh, M. K., Bekki, K., Groves, B. A., & Pfeffer, J. 2022, MNRAS, 510, 5164, doi: 10.1093/mnras/stab3786
  * Chiba & Kataria (2024) Chiba, R., & Kataria, S. K. 2024, MNRAS, 528, 4115, doi: 10.1093/mnras/stae288
  * Clarke & Gerhard (2022) Clarke, J. P., & Gerhard, O. 2022, MNRAS, 512, 2171, doi: 10.1093/mnras/stac603
  * Collier & Madigan (2021) Collier, A., & Madigan, A.-M. 2021, ApJ, 915, 23, doi: 10.3847/1538-4357/ac004d
  * Collier et al. (2018) Collier, A., Shlosman, I., & Heller, C. 2018, MNRAS, 476, 1331, doi: 10.1093/mnras/sty270
  * Collier et al. (2019a) —. 2019a, MNRAS, 489, 3102, doi: 10.1093/mnras/stz2327
  * Collier et al. (2019b) —. 2019b, MNRAS, 488, 5788, doi: 10.1093/mnras/stz2144
  * Combes et al. (1990) Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A&A, 233, 82
  * Combes & Sanders (1981) Combes, F., & Sanders, R. H. 1981, A&A, 96, 164
  * de Vaucouleurs (1963) de Vaucouleurs, G. 1963, ApJ, 138, 934, doi: 10.1086/147696
  * Debattista et al. (2006) Debattista, V. P., Mayer, L., Carollo, C. M., et al. 2006, ApJ, 645, 209, doi: 10.1086/504147
  * Debattista & Sellwood (2000) Debattista, V. P., & Sellwood, J. A. 2000, ApJ, 543, 704, doi: 10.1086/317148
  * Di Matteo et al. (2013) Di Matteo, P., Haywood, M., Combes, F., Semelin, B., & Snaith, O. N. 2013, A&A, 553, A102, doi: 10.1051/0004-6361/201220539
  * Díaz-García et al. (2019) Díaz-García, S., Salo, H., Knapen, J. H., & Herrera-Endoqui, M. 2019, A&A, 631, A94, doi: 10.1051/0004-6361/201936000
  * Díaz-García et al. (2016) Díaz-García, S., Salo, H., Laurikainen, E., & Herrera-Endoqui, M. 2016, A&A, 587, A160, doi: 10.1051/0004-6361/201526161
  * Dillamore et al. (2023) Dillamore, A. M., Belokurov, V., Evans, N. W., & Font, A. S. 2023, MNRAS, 519, L87, doi: 10.1093/mnrasl/slac158
  * Efstathiou et al. (1982) Efstathiou, G., Lake, G., & Negroponte, J. 1982, MNRAS, 199, 1069, doi: 10.1093/mnras/199.4.1069
  * Fall & Efstathiou (1980) Fall, S. M., & Efstathiou, G. 1980, MNRAS, 193, 189, doi: 10.1093/mnras/193.2.189
  * Fragkoudi et al. (2017) Fragkoudi, F., Di Matteo, P., Haywood, M., et al. 2017, A&A, 606, A47, doi: 10.1051/0004-6361/201630244
  * Frankel et al. (2022) Frankel, N., Pillepich, A., Rix, H.-W., et al. 2022, ApJ, 940, 61, doi: 10.3847/1538-4357/ac9972
  * Gadotti (2011) Gadotti, D. A. 2011, MNRAS, 415, 3308, doi: 10.1111/j.1365-2966.2011.18945.x
  * Ghez et al. (2008) Ghez, A. M., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 1044, doi: 10.1086/592738
  * Ghosh et al. (2021) Ghosh, S., Saha, K., Di Matteo, P., & Combes, F. 2021, MNRAS, 502, 3085, doi: 10.1093/mnras/stab238
  * Gillessen et al. (2009a) Gillessen, S., Eisenhauer, F., Fritz, T. K., et al. 2009a, ApJ, 707, L114, doi: 10.1088/0004-637X/707/2/L114
  * Gillessen et al. (2009b) Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009b, ApJ, 692, 1075, doi: 10.1088/0004-637X/692/2/1075
  * Grand et al. (2019) Grand, R. J. J., van de Voort, F., Zjupa, J., et al. 2019, MNRAS, 490, 4786, doi: 10.1093/mnras/stz2928
  * Halle et al. (2018) Halle, A., Di Matteo, P., Haywood, M., & Combes, F. 2018, A&A, 616, A86, doi: 10.1051/0004-6361/201832603
  * Helmi (2020) Helmi, A. 2020, ARA&A, 58, 205, doi: 10.1146/annurev-astro-032620-021917
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hetznecker & Burkert (2006) Hetznecker, H., & Burkert, A. 2006, MNRAS, 370, 1905, doi: 10.1111/j.1365-2966.2006.10616.x
  * Hohl (1976) Hohl, F. 1976, AJ, 81, 30, doi: 10.1086/111849
  * Hoyle (1951) Hoyle, F. 1951, in Problems of Cosmical Aerodynamics, 195
  * Iannuzzi & Athanassoula (2015) Iannuzzi, F., & Athanassoula, E. 2015, MNRAS, 450, 2514, doi: 10.1093/mnras/stv764
  * Iles et al. (2024) Iles, E. J., Pettitt, A. R., Okamoto, T., & Kawata, D. 2024, MNRAS, 527, 2799, doi: 10.1093/mnras/stad3377
  * Jang & Kim (2023) Jang, D., & Kim, W.-T. 2023, ApJ, 942, 106, doi: 10.3847/1538-4357/aca7bc
  * Jiang et al. (2019) Jiang, F., Dekel, A., Kneller, O., et al. 2019, MNRAS, 488, 4801, doi: 10.1093/mnras/stz1952
  * Kataria & Das (2018) Kataria, S. K., & Das, M. 2018, MNRAS, 475, 1653, doi: 10.1093/mnras/stx3279
  * Kataria & Shen (2022) Kataria, S. K., & Shen, J. 2022, ApJ, 940, 175, doi: 10.3847/1538-4357/ac9df1
  * Kawata et al. (2017) Kawata, D., Grand, R. J. J., Gibson, B. K., et al. 2017, MNRAS, 464, 702, doi: 10.1093/mnras/stw2363
  * Kwak et al. (2017) Kwak, S., Kim, W.-T., Rey, S.-C., & Kim, S. 2017, ApJ, 839, 24, doi: 10.3847/1538-4357/aa674c
  * Kwak et al. (2019) Kwak, S., Kim, W.-T., Rey, S.-C., & Quinn, T. R. 2019, ApJ, 887, 139, doi: 10.3847/1538-4357/ab5716
  * Laurikainen et al. (2004) Laurikainen, E., Salo, H., & Buta, R. 2004, ApJ, 607, 103, doi: 10.1086/383462
  * Lee et al. (2022) Lee, Y. H., Park, M.-G., Hwang, H. S., et al. 2022, ApJ, 926, 58, doi: 10.3847/1538-4357/ac3bc1
  * Li et al. (2023) Li, X., Shlosman, I., Pfenniger, D., & Heller, C. 2023, MNRAS, 520, 1243, doi: 10.1093/mnras/stad076
  * Li et al. (2024) —. 2024, MNRAS, 527, 11026, doi: 10.1093/mnras/stad3907
  * Lieb et al. (2022) Lieb, E., Collier, A., & Madigan, A.-M. 2022, MNRAS, 509, 685, doi: 10.1093/mnras/stab2904
  * Łokas (2019) Łokas, E. L. 2019, A&A, 629, A52, doi: 10.1051/0004-6361/201936056
  * Long et al. (2014) Long, S., Shlosman, I., & Heller, C. 2014, ApJ, 783, L18, doi: 10.1088/2041-8205/783/1/L18
  * Lynden-Bell (1960) Lynden-Bell, D. 1960, MNRAS, 120, 204, doi: 10.1093/mnras/120.3.204
  * Marinova & Jogee (2007) Marinova, I., & Jogee, S. 2007, ApJ, 659, 1176, doi: 10.1086/512355
  * Martinez-Valpuesta et al. (2006) Martinez-Valpuesta, I., Shlosman, I., & Heller, C. 2006, ApJ, 637, 214, doi: 10.1086/498338
  * Méndez-Abreu et al. (2012) Méndez-Abreu, J., Sánchez-Janssen, R., Aguerri, J. A. L., Corsini, E. M., & Zarattini, S. 2012, ApJ, 761, L6, doi: 10.1088/2041-8205/761/1/L6
  * Menéndez-Delmestre et al. (2007) Menéndez-Delmestre, K., Sheth, K., Schinnerer, E., Jarrett, T. H., & Scoville, N. Z. 2007, ApJ, 657, 790, doi: 10.1086/511025
  * Merritt & Sellwood (1994) Merritt, D., & Sellwood, J. A. 1994, ApJ, 425, 551, doi: 10.1086/174005
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * Obreja et al. (2022) Obreja, A., Buck, T., & Macciò, A. V. 2022, A&A, 657, A15, doi: 10.1051/0004-6361/202140983
  * Oh et al. (2008) Oh, S. H., Kim, W.-T., Lee, H. M., & Kim, J. 2008, ApJ, 683, 94, doi: 10.1086/588184
  * Ostriker & Peebles (1973) Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467, doi: 10.1086/152513
  * Patsis et al. (2002) Patsis, P. A., Skokos, C., & Athanassoula, E. 2002, MNRAS, 337, 578, doi: 10.1046/j.1365-8711.2002.05943.x
  * Peebles (1969) Peebles, P. J. E. 1969, ApJ, 155, 393, doi: 10.1086/149876
  * Peebles (1971) —. 1971, A&A, 11, 377
  * Peschken & Łokas (2019) Peschken, N., & Łokas, E. L. 2019, MNRAS, 483, 2721, doi: 10.1093/mnras/sty3277
  * Petersen et al. (2016) Petersen, M. S., Weinberg, M. D., & Katz, N. 2016, MNRAS, 463, 1952, doi: 10.1093/mnras/stw2141
  * Pfenniger & Friedli (1991) Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75
  * Portail et al. (2017) Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017, MNRAS, 470, 1233, doi: 10.1093/mnras/stx1293
  * Raha et al. (1991) Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Nature, 352, 411, doi: 10.1038/352411a0
  * Roshan et al. (2021) Roshan, M., Ghafourian, N., Kashfi, T., et al. 2021, MNRAS, 508, 926, doi: 10.1093/mnras/stab2553
  * Saha & Elmegreen (2018) Saha, K., & Elmegreen, B. 2018, ApJ, 858, 24, doi: 10.3847/1538-4357/aabacd
  * Saha & Naab (2013) Saha, K., & Naab, T. 2013, MNRAS, 434, 1287, doi: 10.1093/mnras/stt1088
  * Sellwood (1980) Sellwood, J. A. 1980, A&A, 89, 296
  * Sellwood & Athanassoula (1986) Sellwood, J. A., & Athanassoula, E. 1986, MNRAS, 221, 195, doi: 10.1093/mnras/221.2.195
  * Sellwood & Carlberg (1984) Sellwood, J. A., & Carlberg, R. G. 1984, ApJ, 282, 61, doi: 10.1086/162176
  * Sellwood & Gerhard (2020) Sellwood, J. A., & Gerhard, O. 2020, MNRAS, 495, 3175, doi: 10.1093/mnras/staa1336
  * Seo et al. (2019) Seo, W.-Y., Kim, W.-T., Kwak, S., et al. 2019, ApJ, 872, 5, doi: 10.3847/1538-4357/aafc5f
  * Shen et al. (2010) Shen, J., Rich, R. M., Kormendy, J., et al. 2010, ApJ, 720, L72, doi: 10.1088/2041-8205/720/1/L72
  * Smirnov & Sotnikova (2018) Smirnov, A. A., & Sotnikova, N. Y. 2018, MNRAS, 481, 4058, doi: 10.1093/mnras/sty2423
  * Sormani et al. (2015) Sormani, M. C., Binney, J., & Magorrian, J. 2015, MNRAS, 454, 1818, doi: 10.1093/mnras/stv2067
  * Springel et al. (2021) Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871, doi: 10.1093/mnras/stab1855
  * Toomre (1966) Toomre, A. 1966, inn Geophysical Fluid Dynamics Ref. No. 66-46, ed. W. V. R. Malkus, (Woods Hole, MA: Woods Hole Oceanographic Institute), 111
  * Toomre (1981) Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, ed. S. M. Fall & D. Lynden-Bell, 111–136
  * Villa-Vargas et al. (2009) Villa-Vargas, J., Shlosman, I., & Heller, C. 2009, ApJ, 707, 218, doi: 10.1088/0004-637X/707/1/218
  * Vitvitska et al. (2002) Vitvitska, M., Klypin, A. A., Kravtsov, A. V., et al. 2002, ApJ, 581, 799, doi: 10.1086/344361
  * Wegg et al. (2015) Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050, doi: 10.1093/mnras/stv745
  * Weinberg (1985) Weinberg, M. D. 1985, MNRAS, 213, 451, doi: 10.1093/mnras/213.3.451
  * White (1984) White, S. D. M. 1984, ApJ, 286, 38, doi: 10.1086/162573
  * Whyte et al. (2002) Whyte, L. F., Abraham, R. G., Merrifield, M. R., et al. 2002, MNRAS, 336, 1281, doi: 10.1046/j.1365-8711.2002.05879.x
  * Xiang et al. (2021) Xiang, K. M., Nataf, D. M., Athanassoula, E., et al. 2021, ApJ, 909, 125, doi: 10.3847/1538-4357/abdab5
  * Yurin & Springel (2014) Yurin, D., & Springel, V. 2014, MNRAS, 444, 62, doi: 10.1093/mnras/stu1421
  * Zana et al. (2018) Zana, T., Dotti, M., Capelo, P. R., et al. 2018, MNRAS, 473, 2608, doi: 10.1093/mnras/stx2503

  *[BPS]: boxy peanut-shaped