Mon. Not. R. Astron. Soc. 000, 000–000 (0000) Printed 18 December 2019 (MN LATEX style file v2.

2)

Violent Buckling Benefits Galactic Bars

Angela
1
Collier1?
JILA and Department of Astrophysical and Planetary Sciences, CU Boulder, Boulder, CO 80309, USA
arXiv:1912.08190v1 [[Link]] 17 Dec 2019

Accepted ?; Received ??; in original form ???

ABSTRACT

Galactic bars are unstable to a vertical buckling instability which heats the disk
and in some cases forms a boxy/peanut shaped bulge. We analyze the buckling insta-
bility as an application of classical Euler buckling followed by nonlinear gravitational
Landau damping in the collisionless system. We find that the buckling instability is
dictated by the kinematic properties and geometry of the bar. The analytical result is
compared to simulations of isolated galaxies containing the disk and dark matter com-
ponents. Our results demonstrate that violent buckling does not destroy bars while a
less energetic buckling can dissolve the bar. The disks that undergo gentle buckling
remain stable to bar formation which may explain the observed bar fraction in the
local universe. Our results align with the results from recent surveys.
Key words: methods: numerical — galaxies: evolution, galaxies: interactions —
galaxies: kinematic & dynamics

1 INTRODUCTION shaped bulge (B/P bulge) often attributed to the buck-

ling instability (e.g., Pfenniger & Norman 1990; Raha et al.
Over 70% of disk galaxies host galactic bars. For a review
1991; Athanassoula & Misiriotis 2002; Martinez-Valpuesta
of observations and bar fraction measurements, see Erwin
& Shlosman 2004; Debattista et al. 2006; Saha et al. 2013).
(2018) and the references therein. The galactic bar is im-
Athanassoula & Misiriotis (2002) found that the B/P ap-
portant for galaxy morphology and evolution. It not only
pearance continues to grow in size along with the bar, long
shapes the disk; the bar is responsible for angular momen-
after buckling has subsided. The Milky Way bar is believed
tum transfer throughout the galaxy. The earliest simulations
to have buckled in the past and formed the observed ’X-
of isolated systems show bars forming spontaneously in ax-
shaped’ bulge (Shen et al. 2010; Martinez-Valpuesta & Ger-
isymmetric disks (e.g., Ostriker & Peebles 1973; Hohl 1975).
hard 2011). Surveys have found the buckling could plausibly
Simulated bars can also form from tidal interactions (e.g.,
account for the large fraction of B/P bulges in the local uni-
Noguchi 1987; Gerin et al. 1990). The ease with which simu-
verse (Erwin & Debattista 2016).
lated disks form bars and their pervasiveness in observations
raises the question, how can some disks avoid bar formation? Therefore, buckling is of great importance to galactic
Regardless of the origin of the bar, the increase in radial evolution and deserves further study. Recent works have be-
dispersion velocities along the growing bar makes the disk gun to highlight how halo angular momentum can contribute
increasingly unstable to buckling modes. The bar quickly to stellar disk evolution, particularly its effect on the long
undergoes a second instability which thickens and heats the term evolution of bars. While the angular momentum found
disk vertically. The buckling instability first appeared in sim- in dark matter halos is inconsequential dynamically (Bul-
ulations when three dimensional models were evolved (e.g., lock et al. 2001), there are many effects seen in bar evolu-
Combes et al. 1990; Pfenniger & Friedli 1991; Raha et al. tion (e.g., Saha & Naab 2013; Collier et al. 2018, 2019 a,b).
1991). Following the discovery of buckling in thin disks, the First, the more rotation in the halo the shorter the timescale
structure of the orbit families of the three dimensional bar of the initial bar instability. Secondly and more related to
were also studied (e.g., Pfenniger & Friedli 1991; Skokos et this work, the buckling instability leads to a larger reduc-
al. 2002; Martinez-Valpuesta & Shlosman 2004). Bar orbits tion in bar strength when halos have angular momentum.
develop vertical structure (such as banana and anti-banana Finally, in halos of high spin the bar is dissolved during the
orbits) to support the thickened bar. buckling and does not reform within the disk for the rest of
Edge-on the bars sometimes have a ’boxy/peanut’ its evolution. Collier et al. (2019 a) studied the long term
evolution of the resulting bar dissolution in great detail and
found it to be due to the formation of a relatively strong dark
? E-mail: [Link]@[Link] matter or ’ghost’ bar found in spinning halos. The trapping

c 0000 RAS
2 Angela Collier
persion velocities (σz ) to radial dispersion velocities (σr )
is >
∼ 0.3. Outside this range the buckling modes appear in
the thin sheets and are compared to the fire-hose instabil-
ity of plasma physics. This idealized analysis does not apply
strictly to bars but it is appropriate to follow the radial and
vertical dispersion velocities in the bar due to the drastic in-
crease in vertical dispersion velocity during buckling. In iso-
lated systems, the bar forms in the initially thin disk which
increases the radial dispersion velocities while the vertical
dispersion velocities remain unchanged until the eventual
buckling. This has led to other works discussing the impor-
tance of σr alone in kinematic fractiation of the edge on disk
(Debattista et al. 2017).
Figure 1. a. Indicates the direction of compression and tension An alternate instability was suggested by Combes et
on the bar just before buckling. These forces cause the bar to al. (1990) and Pfenniger & Friedli (1991) after studying the
bend and the magnitude of these forces determine the buckling three-dimensional orbits of stars in barred potentials. They
moment of force (M ). 2. Cartoon of buckled bar. Direction of found that orbits close to a 2:1 vertical resonance are un-
buckling moment of force is shown. Note at the instantaneous stable causing the bar to bend out of the plane. Merritt &
time of buckling the bars edges remain pinned to the galactic Sellwood (1994) however, stressed that this type of instabil-
plane conserving Rb . The maximum deflection is marked as δo .
ity cannot fully explain buckling. While not all stars lie close
to the resonance, all stars trapped in the bar are perturbed
of dark matter material by the stellar bar changes the mass, by the initial buckling mode.
shape and orbital structure of the galactic bar which in turn Once the buckling perturbation appears the bar oscil-
changes the initial conditions of the buckling instability. lates above and below the plane until the buckling is washed
The discussion of buckling is usually limited to the re- out. Both phases of the buckling instability are described be-
sulting bar structure as there is no complete analytical de- low. The second phase of buckling will be analyzed as a form
scription of buckling instability that predicts the time or of Landau damping (Landau 1946). Landau damping as a
mechanics of its onset within a barred galaxy. Additionally, method of collisionless decay of oscillations was initially pre-
the timescale of the instability is quite short compared to scripted for plasmas but can be converted to describe grav-
bar lifetimes. Here we discuss the detailed physical processes itational systems as well (e.g., Lynden-Bell (1962) Maoz
involved in the buckling instability by studying what bar (1991) and Kandrup (1998)). Gravitational Landau damp-
properties effect the initial buckle as well as how the pertur- ing will arise when individual particles interact with a wave
bation moves through the system. We provide a description potential (such as a buckling wave). With Landau damping,
of the mixing that washes out the buckling perturbation. particles with a velocity larger than the wave will lose energy
We ask; what bar properties predict the initial buckling de- to the wave while particles with a lower velocity will gain
flection from the plane? How does the shape of the buck- energy which damps the wave. Binney & Tremaine (2008)
ling perturbation effect the time scale and energetics of the make the analogy of a surfer riding an ocean wave. When
event? And finally, how do different buckling profiles effect the surfer passes behind the crest of the wave they do work
the long term evolution of the disk? on the wave but when the surfer points the board in the
It must be mentioned that Martinez-Valpuesta & Shlos- downhill direction of the wave the wave does work on the
man (2006) have shown that a bar that grows substantially surfer. To extend this analogy, the closer the surfer’s veloc-
in length may undergo a second buckling instability which ity is to the velocity of the wave the longer the surfer can
has structurally different evolution than the first buckling. hold the energy gained by the wave. This is true for Lan-
The second buckling instability will be discussed in a future dau damped particles as well. The Maxwellian distribution
work. of particle velocities implies that some particles will have
This paper is structured as follows. The two phases of a velocity near the velocity of the perturbation and some
the buckling instability are described in Section 2. Next, will be very far away from it. Those far away will gain and
simulations are outlined in Section 3 and the results are lose energy quickly while those particles with energy near
presented in Section 4. We discuss the simulations in the the perturbation will ’hold’ the energy given to them by the
context of the analytical results in Section 5 followed by wave for a longer time. Therefore, a perturbation with a ve-
a brief summary of conclusions and comparison to recent locity near the mean of the system velocity will take longer
observations in Section 6. to damp than a perturbation with a much larger velocity.
We separate the buckling instability into two phases,
an instantaneous breaking of symmetry in the x-z plane fol-
lowed by the oscillation and subsequent washing out of the
2 THE BUCKLING INSTABILITY density perturbation.
After the bar instability breaks the symmetry of the disk
the bar grows in length and strength by trapping additional
2.1 Symmetry breaking by buckling
stellar orbits. The initial edge-on buckling perturbation in
the disk develops spontaneously. Toomre (1966) analyzed The galactic bar feels increasing internal pressures until it
buckling modes in idealized, thin sheets with uniform den- self-buckles under self-gravity. Just before the buckling in-
sity and found the sheets stable if the ratio of vertical dis- stability the bar is aligned with the plane of the thin disk

c 0000 RAS, MNRAS 000, 000–000
 Violent Buckling Benefits Galactic Bars 3

so the initial buckling perturbation is that of a bar pinned Ω. Both of these properties will vary with initial conditions
at both ends with internal pressure forces acting axially. of the buckling wave described in Section 2.1. This analysis
Figure 1 depicts the coordinate system and position of the is done in the rotating frame and we ignore motion of the
compression and tension forces on the bar at the moment of stars in the radial direction. The equation of motion of a
buckling. The bottom of Figure 1 exaggerates the bending single oscillator in the presence of the buckling oscillation is
perturbation of the bar as it buckles. In order to continue
to to grow in strength and length the bar must thicken via
z̈ + ω 2 z = δo cos(Ωt) (4)
the buckling instability to support itself.
The axial stress force on the bar that eventually leads where ω is the vertical oscillation frequency of the star
to buckling can be approximated from dispersion velocities. above and below the plane. Every harmonic oscillator has
The stress tensor of a stellar fluid (τ ) is related to the dis- a unique ω due to the continuous distribution of velocities
persion velocities by Binney & Tremaine (2008) where in the z direction. The general solution to the equation of
2 motion is
τ = −ρ(r)σi,j . (1)
The volume density is ρ(r) and σi,j, is defined as sin(ωt)
2 z(t) = zo cos(ωt) + z˙o +
σi,j = h(vi − hvi i)(vj − hvj i)i (2) ω (5)
δo
where i, j = r, φ, z are in cylindrical coordinates. We [cos(Ωt) − cos(ωt)].
ω 2 − Ω2
separate τ into the normal and shear stresses
To study the response of the stars to the buckling we
2 2
τ = τn + τs = −ρ(r)(σii + σij ). (3) ignore the first two terms without Ω. These terms relate
It follows that for the galactic bar, the ratio of stress to to decoherence among the stellar ensemble while we focus
strain is τs /τn . We compare the buckling of the galactic bar solution to the response ro the perturbation (zres ):
to the classical Euler buckling of columns where the critical
stress that buckles the column is proportional to Young’s δo
zres (t) = [cos(Ωt) − cos(ωt)]. (6)
modulus, the ratio of stress to strain in the material. We ω 2 − Ω2
suggest the energy of the buckling instability should depend We calculate the alignment along the plane for the en-
on the ratio of internal stresses felt by the bar at the time tire bar by averaging over all stars
of buckling.
The geometric properties of the bar should play a role Z ∞
ρ(ω)
in the buckling as well. In the classical Euler buckling in- hzres (t)i = δ(t) = δo dω [cos(Ωt) − cos(ωt)]. (7)
−∞ ω 2 − Ω2
stability of columns the force of buckling is inversely related
to the square of the slenderness of the column. The slen- As the perturbation gets damped δ(t) goes to zero. The
derness ratio of the bar (Rb /zb ) is also important in galactic distribution of frequencies (ρ(ω)) is centered around some
buckling but we predict that a more slender galactic bar will ω̄ and the frequency of the perturbation must be on the
have a more dramatic/stronger buckling when compared to order of this average oscillation. Consequently, we can ex-
a thicker bar. This prediction is based on a similar orbital pand Equation 7 with ω = Ω + (ω − Ω) and approximate the
inclination instability found in thin disks. Madigan et al. average displacement from the plane as
(2018) found that the disk inclination instability is driven
by torques between orbits which are stronger when the disks δo
 Z ∞
1 − cos(ω − Ω)t

are thinner. Comparing this to a stellar bar, we expect the δ(t) ≈ cos(Ωt) dωρ(ω)
2ω̄ −∞ ω−Ω
orbits within a thicker bar to experience less torques result-  Z ∞  (8)
ing in an overall weaker buckling. The energy of the buckling δo sin(ω − Ω)t
+ sin(Ωt) dωρ(ω) .
is reflected in the deflection from the plane, the more slen- 2ω̄ −∞ ω−Ω
der bar will have a larger deflection from the plane (due to We assume that the distribution around ω̄ is small
the increased torques increasing the energy of the buckling) which reduces Equation 8 to
when compared to a thicker bar. The maximum deflection
of the bar, δo , is labeled in Figure 1. This parameter should ∞
sin 12 (ω − Ω)t
Z 
increase with buckling energy as energy increases with the δo sin(ω̄t)
δ(t) ≈ dωρ(ω) . (9)
square of the amplitude of the sinusoidal buckling force. We ω̄ −∞ ω−Ω
predict the kinematic and geometric properties of the galac- A single oscillator with frequency ω will find itself vi-
tic bar have an effect on the initial buckling wave with the brating due the buckling perturbation according to
energy of buckling increasing with internal stresses and bar
slenderness.
δo sin 12 (ω − Ω)t
z(t) ≈ . (10)
ω̄ (ω − Ω)
2.2 Landau damping of the buckling wave
Hence, a star with frequency ω will be excited after the
The buckling causes a sinusoidal force to be exerted onto the initial buckling and will move to a position z = δ0 /[ω̄(ω−Ω)]
stars and then the buckling perturbation evolves as a wave on a timescale of t = π/[ω − Ω] and then will align with
oscillating above and below the plane. To analyze the damp- the plane after t ≈ 2π/[ω − Ω]. This star will experience
ing of the system we describe each star as a harmonic oscil- repeat oscillations of lower amplitude as the perturbation
lator perturbed by this force of amplitude δ and frequency gets damped. For stars with ω very near or equal to Ω the

c 0000 RAS, MNRAS 000, 000–000
4 Angela Collier

energy given to the star by the wave is held for a longer

time compared to stars with ω very far from Ω which is
a fundamental property of gravitational Landau damping.
(Lynden-Bell 1962) Initially, nearly all stars in the bar con-
tribute to the oscillation but this fraction gets lower as t goes
to infinity. Eventually, only stars with ω = Ω contribute to
the oscillation—if those stars exist.
At the time of the initial perturbation all oscillators are
driven in phase with the initial displacement. As the insta-
bility progresses more oscillators become out phase and no
longer interact with the perturbation which causes destruc-
tive interference. This leads to the perturbation decaying
rapidly. The energy of the perturbation is converted into in-
coherent motion associated with the oscillations about the
plane. The timescale of the Landau damping will be effected
by the bar parameters at the time of buckling which deter-
mine δo and Ω. A higher energy wave will have a larger
amplitude and frequency and with each oscillation a larger
fraction of energy will be moved from the wave motion to
the individual particle motion. The opposite is true for lower
energy waves. If the initial Ω of the perturbation is lower and
closer to the ω̄ of the system the change in energy per oscil-
lation will be smaller. In this case, the individual oscillators
will spend more time in phase with the perturbation which
lengthens the damping time scale. The time scale is also
lengthened by the lower frequency as the time to complete
one oscillation increases.

3 NUMERICS
The issue of bar buckling in the context of N-body simu-
lations of disk galaxies is now addressed. We model stellar
disks inside a spherical Navarro-Frenk-White halos (Navarro Figure 2. Top: Evolution of the bar strength parameter (A2 /A0 )
et al. 1996, hereafter NFW) using the N -body part of in the x-y plane which compares the growth of the Fourier m = 2
the tree-particle-mesh Smoothed Particle Hydrodynamics mode to m = 0 mode the for three analyzed models. Bottom:
(SPH/N -body) code GIZMO (Hopkins 2015). Our code Evolution of the ratio of the Fourier m = 1 mode to the m = 0
units for mass, distance, and time are 1010 M , 1 kpc, and mode (A1z /A0 ) in the x-z plane for three galaxy models. This
1 Gyr. The time step between snapshots is shortened to quantity represents the growth in asymmetry of the edge-on bar.
0.001 Gyr during the buckling in order to better study the
instability which happens on quite a short timescale.
The simulations presented in this work are adapted
from the models described in Collier et al. (2018, 2019 a,b). iterative method from Rodionov & Sotnikova (2006), see also
We refer the reader to these works for more analysis on the Rodionov et al. (2009). The iterations create a halo with an
long term evolution of these models. isotropic velocity distribution while cosmological halos are
found to have a range of spin which can be fit by a lognormal
distribution,
3.1 Initial Conditions of Galaxy Models
We have created a series of two-component isolated galaxies ln2 (λ/λ0 )
 
1
with identical disks inside halos of different rotation. The P (λ) = exp − , (12)
λ(2πσλ )1/2 2σλ2
halo density follows the NFW profile,
where λ0 = 0.035 ± 0.005 and σλ = 0.5 ± 0.03 are the fitting
−(r/rt )2 parameters (Bullock et al. 2001).
ρs e
ρh (r) = , (11) To add angular momentum to the halo a fraction of
[(r + rc )/rs ](1 + r/rs )2
randomly chosen halo particles have their tangential veloc-
where ρ(r) is the DM density in spherical coordinates, ρs is ities reversed. The new velocity distribution maintains the
the fitting density parameter, and rs = 9 kpc is the charac- solution to the Boltzmann equation and does not alter the
teristic radius, where the power law slope is −2, and rc is a velocity profile (Lynden-Bell 1960; Weinberg 1985; Collier
central density core where rc = 1.4 kpc. The Gaussian cutoff et al. 2018, 2019 a,b), so the equilibrium state is preserved.
is applied at rt = 86 kpc for the halo. The DM halo contains While our halos vary in spin, each simulation begins
7.2 × 106 particles and the halo mass is Mh = 6.3 × 1011 M . with an identical disk. The volume density of the exponential
The halo velocities are found by using a version of the stellar disk is;

c 0000 RAS, MNRAS 000, 000–000
 Violent Buckling Benefits Galactic Bars 5

Nd
A2 1 X
= mi e2iφi , (14)
A0 A0 i=1

which is achieved by summing over all bar particles with

R ≤ 14 kpc, and mass m = mi at azimuthal angle φi .
We plot the bar strength for all three models in the top
of Figure 2 for the entirety of the simulation. While this
paper strictly focuses on buckling event itself, note that the
long term evolution of these models varies after buckling.
This effect is discussed in Collier et al. (2018, 2019 a) in great
detail. Briefly, these models vary only in λ and the disks are
identical at t = 0. The resulting changes in bar strength
evolution can be seen in the top of Figure 2. The disks start
the simulations axisymmetric. When the disk undergoes the
bar instability there is an exponential growth in bar strength
that includes the trapping of stellar and halo particles in the
bar. The bar reaches some peak strength and then undergoes
the buckling instability which reduces the strength of the bar
and in some models dissolves the bar, i.e. P90. In halos of
lower spin the bar recovers after buckling and continues to
grow in strength for the remainder of the simulation.
Halo spin effects the evolution of the bar in the following
ways. First, the timescale of the bar instability is shortened
in spinning halos. For example, the bar in the P90 model
Figure 3. Top: Evolution of the radial dispersion velocities nor-
malized by σr (t = 0), Middle: Evolution of vertical dispersion appears nearly two Gyrs before the bar in the P00 model.
velocities normalized by σz (t = 0), Bottom: Evolution of the ra- Secondly, change in strength during the buckling in stabil-
tio of vertical to radial dispersion velocities (σz /σr ). For all plots ity (∆A2 /A0 ) increases with increasing halo spin. The P90
the data was measured for particles at the radii between 3 − 6 bar loses ∼ 80% of its strength after buckling while the P00
kpc. bar loses only ∼ 40%. Third, we see the time scale of the
buckling instability increases with increasing halo spin. By
buckling time scale (tb ) we are referring to the time between
the maximum value of A2 /A0 just before the minimum of
A2 /A0 . For the P90 model, tb = 0.72 Gyrs and for the P00
Md
z
model, tb = 0.2 Gyrs. Finally, the post buckling evolution
ρd (R, z) = exp(−R/h) sech2 , (13)
4πh2 z0 z0 depends greatly on the spin of the halo, with the bar dis-
solving in P90. Obviously, this effect must be intensely in-
where Md is the disk mass, h = 2.85 kpc is its radial scale- tertwined with the dynamics of the buckling instability.
length, and z0 = 0.6 kpc is the scaleheight. The stellar The growth of the Fourier m = 1 mode compared to
disk has 0.8 × 106 particles and the disk mass is Md = the m = 0 mode represents the growth in asymmetry in the
6.3 × 1010 M . The initial disk velocities depend on the po- x-z plane. As the bar buckles out of the plane we expect to
tential of both components of the galaxy and are assigned see this parameter increase in strength. We plot this for each
last. The radial and vertical dispersion velocities are as- model in the bottom of Figure 2. The peak in this parameter
signed as exponentials. generally matches the time of the peak in A2 /A0 just before
The galaxy simulations presented here vary only in halo the buckling instability.
spin. Following the notation of Collier et al. (2018, 2019 a,b) As discussed in the Section 1, the dispersion velocities
the models are labeled as P for prograde rotation and then in the radial direction increase as more and more stars be-
multiplied by λ of the halo. For example, the model P00, is come trapped in the bar. In the z-direction the orbits have
the nonspinning halo with λ = 0. random motions and oscillate above and below the plane.
While the bar is growing in the disk and increasing the ra-
dial dispersion velocities; the vertical dispersion velocities
do not increase until the onset of the buckling instability
in these isolated models. Figure 3 shows the time evolution
4 RESULTS FROM NUMERICAL MODELS
of the radial and vertical dispersion velocities normalized
The following analysis is done considering the trapped dark by their t = 0 values, as well as the ratio between these two
matter component as part of the total galactic bar. The measures. This plot was taken for all time steps and includes
dark matter component contributes to the mass, strength, all stellar and DM particles in the bar at the radii between
and height of the bar. (e.g. Collier et al. (2019 a,b)) |3 − 6| kpc to avoid anything odd inside the very central bar.
All three models produce strong bars before buckling. In the top of Figure 3 we plot the evolution of σr nor-
The strength of the bar is defined by the ratio of the Fourier malized by σr (t = 0) for the three models. The growth of
m = 2 mode to the m = 0 mode, σr is closely tied to the bar evolution as trapping stars in

c 0000 RAS, MNRAS 000, 000–000
6 Angela Collier

the bar increases the fraction of radial orbits. The evolution

of σr appears to mirror the evolution of the bar strength
parameter (see top of Figure 2). Radial dispersion velocities
increase after the onset of the bar instability and there is a
dip during the buckling instability. For non rotating and low
spin halos the bars inside recover from buckling and continue
to increase radial dispersion velocities, while the P90 disk is
stagnant after buckling. The radial dispersion velocity does
not reduce to its initial value in the P90 model after the
buckling instability which explains why this disk does not
under go a second bar instability—the disk is too hot.
The middle of Figure 3 shows the time evolution of the
vertical dispersion velocities normalized by σz (t = 0). For
the P00 and P45 models, σz is relatively constant until the
onset of the buckling instability where there is a swift, large
increase. σz for these models appears to saturate and reach
some equilibrium after the buckling instability. The steep-
ness of the growth from σz at t = 0 to σz post buckling
decreases with λ. This change in slope completely washes
out the buckling instability in this plot for the P90 model.
Between t ∼ 1.8 Gyr and t ∼ 4.4 Gyrs the vertical dis-
persion velocities linearly increase until they reach the final
equilibrium value for the P90 model. This is a dramatic dif-
ference when comparing to the P00 and P45 models which
look more like step functions rather than a slow and steady
increases over large times. Figure 4. Time series of mean vz for all models. Outside the
The bottom of Figure 3 displays the time evolution of buckling instability the vertical velocities of the disk and bar
the ratio of σz /σr . When the bar develops and σr increase should average to zero. The spike seen each model coincides the
buckling instability.
we see a dip in the ratio of dispersion velocities. Associated
with the buckling instability is a swift increase in this ratio
as the vertical dispersion velocities increase. For the P00 and
P45 models the ratio reaches a maximum and then begins Model P00 P45 P90
to decrease again, hinting at a future second buckling. For ∆A2 /A0 0.19 0.29 0.37
the P90 model we see the dip in the ratio followed by a small Rb (kpc) 9.4 6.2 6.1
increase before the final jump with the buckling instability. zb (kpc) 1.2 1.7 2.2
Due to σz slowly and linearly increasing before the buckling Mb (1010 M ) 2.79 3.10 3.54
instability in the P90 model, at the time of buckling σz /σr tb (Gyr) 0.2 0.2 0.7
is larger when compared to the ratio for P00 and P45 just
before buckling. Over longer scales for the P90 model we Table 1. Bar parameters for three models. ∆A2 /A0 is the drop in
see no evolution in the dispersion velocity ratio — no second bar strength during buckling. Rb , zb , and Mb are the bar length,
buckling or bar instability will happen for this stagnant disk. height, and mass at the time of buckling. Finally, tb is the time
Having looked at the variance in the velocities in the required for the bar strength to go from the last pre-buckling
models we now turn our attention toward the evolution of maximum to the minimum bar strength.
the mean velocities plotted in Figure 4. While the buckling
wave is coherent we see a brief increase in mean value of vz . 2018). The bar mass (stellar plus dark matter contribution)
Outside the buckling instability the mean velocity about the is measured in Collier et al. (2019 a).
is zero. The spike of the maximum average vz coincides with
the time of maximum buckling. We suggest that this mea-
sure acts as a better indicator of time of the actual buckling
when compared to the bottom of Figure 2 which is a measure 5 DISCUSSION
of disk asymmetry. Note the trend seen in the height of the We now apply the buckling instability analysis from Section
peaks for each model. Increasing halo spin limits the velocity 2 to our models. For the development of the initial buckling
reached by the buckled stars. This implies that the energy perturbation we found the geometric and kinematic proper-
of the buckling is decreasing with increasing halo spin. ties of the bar can effect the energy of the buckling.
To end this section, we compile a table of the resulting From the bottom Figure 3 we note that the ratio of
bar parameters from our simulations. ∆A2/A0 and tb are dispersion velocities for each model at the time of buckling
calculated from the data in Figure 2. To measure Rb and zb is different, with the P90 model having the largest ratio of
we calculate isodensity contours in the xy-plane and rz-plane σz /σr . We can relate the value of the buckling stress (τs /τn )
respectively. The bar extends to the point were the elliptic- to the ratio of vertical and radial dispersion velocities. The
ity of the contour as decreased from its maximum by 15%. ratio of dispersion velocities measured at the time of buck-
(e.g., Martinez-Valpuesta & Shlosman 2006; Collier et al. ling is ∼ 0.4 for the P90 model, ∼ 0.32 for the P45 model,

c 0000 RAS, MNRAS 000, 000–000
 Violent Buckling Benefits Galactic Bars 7

Figure 5. Fourier transform of mean Vz distribution along radii ±20 kpc in the disk with data binned into 0.5 kpc bins. This measurement
was taken at all times relevant to buckling for each individual model. Color represents the buckling power at that position and frequency.
More powerful buckling at higher frequencies is taken to be more energetic.

and the ratio is ∼ 0.29 for the P00 model. Looking only at
the kinematic component we would expect the energy of the
P00, and P45 models to be about equivalent with the en-
ergy being lower in the P90 model. However, the geometric
properties should also effect the buckling energy. From Ta-
ble 1, we calculate the slenderness of the simulated bars at
the time of buckling and find that the most slender bar re-
sides in the P00 model, followed by the P45 model and then
the P90 model. From the dynamical argument in Section
2.1, the more slender bars should buckle with more energy.
By analyzing the energy of the buckling in our models we
can determine the relative importance of the geometric and
kinematic properties of the bar to the buckling energy.

To quantify the energy of the buckling wave in each

simulation, we start by finding the average vz in 0.5 kpc
bins for every relevant time step during buckling for each
model. Next, we apply a Fourier transform to this data and
Figure 6. Calculated change in kinetic energy (KE) in the bar
plot the color map of the power spectrum which is shown in
of each model normalized by the pre-buckling KE. The buckling
Figure 5. The frequency of the buckling wave is plotted for
instability happens at different times for each model. t = 0 in
all relevant radii. Color denotes the power of the buckling this plot is the initial buckling time for easier comparison between
wave on scale from blue (least powerful) to red (most power- models.
ful). Energetic bucklings will be more powerful at more radii
than a less energetic buckling. If slenderness were the only
predictor of buckling energy we would expect the bucklings the ratio of σz /σr at the moment of buckling the stronger
of the P45 model, with a slenderness ratio of 3.6, to be on the pressure of the buckling force.
a similar scale to the P90 model, with a slenderness ratio We now turn to the second phase of buckling and study
of 2.8. However, we see that the bucklng is much less pow- the Landau damping in the simulated bars. We note that we
erful in the P90 model compared to the P45 model. These do not expect these simulations (or the mixing within them)
models with similar geometric properties but very different to be dependent on the particle number, N . Although too
kinematic profiles allow us to see that the kinematic proper- few particles in the disk will allow the bar to decay after
ties are much more important to the initial buckling profile. buckling due to phase mixing effects we do not expect that
We find the most powerful buckling in the P00 model, show- to be relevant here. Particle number was explored in Dubin-
ing deeper reds and larger extent than the other two models. ski, Berentzen & Shlosman (2009). The particle minimum
Note, the bar in the P00 model is longer than the other two needed to avoid those saturation effects is most readily seen
models. As noted above the entire length of the bar is in- in Figure 17 of their paper. We greatly exceed this minimum
volved in the buckling perturbation which is why this model in our simulations.
is spread along more radii than other models. From this fig- In Figure 6 we compare the change in kinetic energy in
ure, we find the kinematic properties of the bar are the most each simulated bar from the time step just before buckling
important for determining the buckling energy. The smaller to long after the buckling has subsided. The buckling time of

c 0000 RAS, MNRAS 000, 000–000
8 Angela Collier

each model is set to t = 0 in this plot for better comparison required for Landau damping to remove the perturbation.
between the different models. The frequency and amplitude This allows the chaos that develops at the end of the bar to
of the wave is higher in the P00 model (blue) and P45 model move to smaller radii, further dissociating bar orbits. In our
(green) when compared to the P90 model (red). The larger models, the bar parameters were changed by the initial con-
the initial energy of the perturbation the more quickly it is ditions of the halo. There is evidently some turnover point
damped out. The P00 and P45 models see more oscillations for these models between λ = 0.045 − 0.09 where the buck-
while the P90 model takes much more time to complete a ling timescale is lengthened and a large enough fraction of
single oscillation. bar orbits become dissociated preventing the bar from recov-
Why should the timescale of buckling effect bar evo- ering from buckling. The work of Collier et al. (2018) allows
lution? A second process develops at the edge of the bar us to limit the value of λ for these models to λ = 0.045−0.06
near corotation while the buckling instability is in progress. However, changing halo spin is not the only way to effect the
It has been known for some time that strong bars gener- initial conditions of the buckling instability.
ate chaos (Contopoulos 1981) and the development of chaos We have shown that a nonviolent buckling can dissolve
during the buckling instability has been well documented in a strong bar. This result provides a possible answer to the
Martinez-Valpuesta & Shlosman (2004). The authors found question in the introduction. How can some disk galaxies
the weakening of the bar happens at large radii and then avoid the bar instability for Hubble time? Perhaps, they
propagates inward due to the increase in chaos during buck- formed a bar and following a weak buckling the disk was
ling. The strength of the bar sharply decreases due to the too hot to undergo a second bar instability. This idea is
process which is happening simultaneously to the Landau supported by observations of unbarred disks which are found
damping. If buckling is lengthened more chaos will be intro- to be kinematically hotter than their barred counterparts
duced to the bar which will weaken the bar further when (Sheth et al. 2012).
compared to bars that had a shorter buckling phase. More While a weak buckling results in bar dissolution a vio-
chaos induced in the outer region of the bar makes the radial lent buckling produces an edge-on B/P bulge. The first ob-
bar orbits less coherent further decreasing the bar strength servational measurements of B/P bulges with redshift have
as seen in Figure 2. The important result here is that a more been published by Kruk et al. (2019). The authors found
volatile buckling instability found in the P00 model, and to the fraction of bars with B/P bulges increased from 10% at
a lesser extent in the P45 model, is not given enough time to z = 1 to 70% at z = 0. They find that buckling is sup-
dissolve the bar. In fact, the bar quickly recovers and con- pressed in thick disks, not because of disk thickness but
tinues to grow. The slow and gentle buckling found in the because of higher vertical dispersion velocities. This find-
P90 model does, however, result in a dissolved galactic bar. ing agrees with our results; in models with larger vertical
dispersion velocities at the time of buckling we find a less
obvious B/P bulge. We end by showing the observational
results of differing buckling instabilities. In Figure 7 we plot
6 CONCLUSIONS
isodensity contours of the stellar disk surface densities for
We have analyzed the buckling instability of a self- the face-on and edge-on disks. Bars that undergo a violent
gravitating bar followed by the Landau damping that washes buckling host a B/P shaped bulge, while those that have an
out the perturbation. We compare these results to three sim- elongated weak buckling host a dissolved bar and no B/P
ulations with different kinematic and geometric properties bulge edge-on.
in an attempt to discover what effect these properties have This work demonstrates the diverse stellar bar evolu-
on the buckling instability and resulting disk evolution. Our tion that can follow buckling of varying initial conditions.
results are summarized in the following paragraphs. However, we have not developed a comprehensive tool to
We found that the slenderness ratio of the galactic bar predict the buckling instability parameters or the time of
effects the energy of the buckling instability, with more slen- the initial symmetry breaking. Additional theoretical work
der bars having a more energetic buckling. For these models must be done in order to gain a complete understanding of
we considered the trapped dark matter component to con- the buckling instability.
tribute to the thickness of the galactic bar but other factors
could play a roll in changing the bar parameters before buck-
ling. For example, the thickness of the stellar disk can delay
the onset of the bar instability and change the results of the
buckling instability (e.g. Aumer & Binney (2017)).
ACKNOWLEDGEMENTS
Of more importance to the energy of the buckling per-
turbation, we found that there is is no simple correlation We gratefully acknowledge support from NASA grant
between a certain value of σz /σr and the buckling insta- 80NSSC17K0720. This manuscript was greatly improved
bility. As seen in the bottom of Figure 3 the time a stable from suggestions made by Isaac Shlosman, Ann-Marie Madi-
barred disk can live below the ratio of σz /σr < 0.4 depends gan, and the reviewer Daniel Pfenniger. Thanks to Jeremiah
on many bar properties. We do find that a larger ratio of Hudson for helpful mechanics discussions. This work utilized
dispersion velocities will inhibit the energy of the buckling the RMACC Summit supercomputer, which is supported
instability and plays a much larger roll when compared to by the National Science Foundation (awards ACI1532235
the geometric properties. and ACI-1532236), the University of Colorado Boulder, and
Finally, an increase in the buckling time scale leads to Colorado State University. The Summit supercomputer is a
a larger bar strength loss (∆A2 /A0 ) during the buckling in- joint effort of the University of Colorado Boulder and Col-
stability. A weaker buckling energy will increase the time orado State University.

c 0000 RAS, MNRAS 000, 000–000
 Violent Buckling Benefits Galactic Bars 9

Figure 7. Linear isodensity contours of a 1 kpc slice of the stellar component for each simulation. Top row is face-on, middle is a view
of the major axis and the bottom row is the bar viewed from the minor axis. λ increases to the right. We see a reduction in bar strength,
disk radii, and spiral arm activity with the increase of λ in the face on view. Edge on the B/P bulge shape decreases with the P90 models
major and minor axis being nearly indistinguishable.

REFERENCES AJ, 156, 141

Maoz, E. 1991, ApJ, 375, 687
Athanassoula E., Misiriotis A., 2002, MNRAS, 330, 35
Martinez-Valpuesta I., Shlosman I., 2004, ApJL, 613, L29
Aumer M., Binney J., 2017, MNRAS, 470, 2113
Martinez-Valpuesta, I., Shlosman, I., & Heller, C. 2006,
Binney, J., Tremaine, S. 2008, Galactic Dynamics, Prince-
ApJ, 637, 214
ton Univ. Press
Martinez-Valpuesta I., Gerhard O., 2011, ApJL, 734, L20
Bullock, J.S., Dekel, A., Kolatt, T.S., Kravtsov, A.V.,
Merritt D., Sellwood J. A., 1994, ApJ, 425, 551
Klypin, A.A., Porciani, C., & Primack, J.R. 2001, ApJ,
Navarro, J.F., Frenk, C.S., White, S.D.M. 1996, ApJ, 462,
555, 240
563 (NFW)
Collier A., Shlosman I., Heller C., 2018, MNRAS, 476, 1331
Noguchi, M. 1987, MNRAS, 228, 635
Collier A., Shlosman I., Heller C., 2019, MNRAS, 488, 5788
Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467
Collier A., Shlosman I., Heller C., 2019, MNRAS, 489, 3102
Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75
Contopoulos, G. 1981, A&A, 102, 265
Pfenniger D., Norman C., 1990, ApJ, 363, 391
Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D.
Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D.
1990, A&A, 233, 82
1991, Nature, 352, 411
Debattista V. P., Mayer L., Carollo C. M., Moore B., Wad- Rodionov, S.A., Athanassoula, E., Sotnikova, [Link]. 2009,
sley J., Quinn T., 2006, ApJ, 645, 209 MNRAS, 392, 904
Debattista V. P., Ness M., Gonzalez O. A., Freeman K., Rodionov, S.A., Sotnikova, [Link]. 2006, [Link]., 50,
Zoccali M., Minniti D., 2017, MNRAS, 469, 1587 983
Dubinski J., Berentzen I., Shlosman I., 2009, ApJ, 697, 293 Saha K., Naab T., 2013, MNRAS, 434, 1287
Erwin P., Debattista V. P., 2016, ApJL, 825, L30 Saha K., Pfenniger D., Taam R. E., 2013, ApJ, 764, 123
Erwin P., 2018, MNRAS, 474, 5372 Shen, J., Rich, R. M., Kormendy, J., et al. 2010, ApJL,
Gerin, M., Combes, F., & Athanassoula, E. 1990, A&A, 720, L72
230, 37 Sheth K., Melbourne J., Elmegreen D. M., Elmegreen
Hohl F., 1975, IAUS, 349, IAUS...69 B. G., Athanassoula E., Abraham R. G., Weiner B. J.,
Hopkins, P.E., 2015, MNRAS, 450, 53 2012, ApJ, 758, 136
Kandrup H. E., 1998, ApJ, 500, 120 Skokos C., Patsis P. A., Athanassoula E., 2002, MNRAS,
Kruk S. J., Erwin P., Debattista V. P., Lintott C., 2019, 333, 847
arXiv, arXiv:1910.04768 Toomre A., 1966, in Geophys. Fluid Dyn., 66, 111
Landau, L. D. 1946,J. Phys.(USSR), 10, 25-34 Weinberg, M.D. 1985, MNRAS, 213, 451
Lynden-Bell, D. 1960, MNRAS, 120, 204
Lynden-Bell D., 1962, MNRAS, 124, 279
Madigan A.-M., Zderic A., McCourt M., Fleisig J., 2018,

c 0000 RAS, MNRAS 000, 000–000