# Orbital support and evolution of flat profiles of bars (shoulders)

Leandro Beraldo e Silva Department of Astronomy, University of Michigan, 1085
S. University Ave., Ann Arbor, MI, 48109 USA  Jeremiah Horrocks Institute,
University of Central Lancashire, Preston PR1 2HE, UK  Victor P. Debattista
Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1
2HE, UK  Stuart Robert Anderson Jeremiah Horrocks Institute, University of
Central Lancashire, Preston PR1 2HE, UK  Monica Valluri Department of
Astronomy and Astrophysics, University of Michigan, Ann Arbor, MI, USA  Peter
Erwin Max Planck Institut für Extraterrestrische Physik, Giessenbachstrasse,
D-85748 Garching, Germany  Kathryne J. Daniel Department of Astronomy &
Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  Center for
Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA
Nathan Deg Department of Physics, Engineering Physics, and Astronomy, Queen’s
University, Kingston, ON, K7L 3N6, Canada Leandro Beraldo e Silva
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Many barred galaxies exhibit upturns (shoulders) in their bar major-axis
density profile. Simulation studies have suggested that shoulders are
supported by looped $x_{1}$ orbits, occur in growing bars, and can appear
after bar-buckling. We investigate the orbital support and evolution of
shoulders via frequency analyses of orbits in simulations. We confirm that
looped orbits are shoulder-supporting, and can remain so, to a lesser extent,
after being vertically thickened. We show that looped orbits appear at the
resonance
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$
(analogous to the classical Inner Lindblad Resonance, and here called ILR)
with vertical-to-radial frequency ratios
$1\lesssim\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 3/2$ (vertically warm
orbits). Cool orbits at the ILR (those with
$\mbox{$\>\Omega_{z}/\Omega_{R}$}>3/2$) are vertically thin and have no loops,
contributing negligibly to shoulders. As bars slow and thicken, either
secularly or by buckling, they populate warm orbits at the ILR. Further
thickening carries these orbits towards crossing the vertical ILR [vILR,
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{z}$}=1/2$],
where they convert in-plane to vertical motion, become chaotic, kinematically
hotter and less shoulder-supporting. Hence, persistent shoulders require bars
to trap new stars, consistent with the need for a growing bar. Since buckling
speeds up trapping on warm orbits at the ILR, it can be followed by shoulder
formation, as seen in simulations. This sequence supports the recent
observational finding that shoulders likely precede the emergence of BP-
bulges. The python module for the frequency analysis, naif, is made available.

Galactic dynamics — Galactic bars — Frequency analysis

††software: numpy (Harris et al., 2020), scipy (Virtanen et al., 2020),
pynbody (Pontzen et al., 2013), Agama (Vasiliev, 2019), naif (this work).

##  1 Introduction

Stellar bars are present in $\approx 60-70\%$ of nearby spiral galaxies
(Eskridge et al., 2000; Menéndez-Delmestre et al., 2007; Sheth et al., 2008;
Nair & Abraham, 2010; Erwin, 2018). Barred galaxies have surface-brightness
profiles along the bar major axis which are traditionally classified as either
exponential or as having, on top of an exponential profile, a nearly-flat part
in the bar outskirts (on both sides) – hereafter the “shoulders” (e.g.
Elmegreen & Elmegreen, 1985; de Carvalho & da Costa, 1987; Elmegreen et al.,
1996; Prieto et al., 2001; Gadotti et al., 2007, but see Erwin et al. (2023)
for an updated classification scheme).

Elmegreen et al. (1996), Regan & Elmegreen (1997) and Elmegreen et al. (2011)
found that shoulders are generally present in profiles of early-type barred
spiral galaxies, while late-type galaxies normally have bars with exponential
profiles. In a sample of 144 barred galaxies from the Spitzer Survey of
Stellar Structure in Galaxies ($S^{4}G$, Sheth et al., 2010), Kim et al.
(2015) found that massive and bulge-dominated galaxies tend to have flat bar
profiles, while bulgeless galaxies tend to have exponential profiles, in
accordance with early predictions based on $N$-body simulations (Combes &
Elmegreen, 1993). The same trends were found by Kruk et al. (2018) in a sample
of 3461 barred galaxies from SDSS, while Erwin et al. (2023) find that the
presence of shoulders is more strongly correlated with stellar mass. Shoulders
in the bar major axis profile have also been observed in M31 (Athanassoula &
Beaton, 2006) and in several simulation studies (e.g. Athanassoula &
Misiriotis, 2002; Bureau & Athanassoula, 2005; Anderson et al., 2022,
hereafter Paper I).

Fig. 1 shows examples of bar major-axis profiles for the galaxies NGC 3681 and
NGC 4340 from $S^{4}G$ (Sheth et al., 2010) and for the simulation SD1
(described in Sec. 2) at 7 Gyr (dotted black/red) and at 10Gyr (dash-dotted
blue/orange). The bars are aligned with the $x$-axis and for the simulation we
select all star-particles within $|y|/\mbox{$a_{\rm bar}$}<0.145$, where
$a_{\rm bar}$ is the bar length – see Sec. 2 for details. For the observed
galaxies, we show one-pixel-wide profiles along the bar major axis from
Spitzer 3.6µm images. NGC 3681 has a pure-exponential profile, while NGC 4340
and the model SD1 have profiles with shoulders (red), as detected by Paper I.

The galaxy NGC 4340 was chosen as an example of a bar profile with
particularly pronounced shoulders, and these look much more pronounced than
the ones in the model SD1 at 7Gyr. The shoulder strength, estimated as the
excess mass of the shoulder over a purely exponential profile, normalized by
the total mass within the shoulder region, typically grows over time, as shown
in Paper I. In fact, we see that the shoulders of SD1 get stronger at 10Gyr
(blue/orange dash-dotted lines), although its length (in units of $a_{\rm
bar}$) is approximately constant and smaller than in the observed galaxy NGC
4340.

Figure 1: Normalised observed bar-major-axis brightness profiles for the
Spitzer $3.6\micron$ band for NGC 3681 and NGC 4340, in arbitrary units (S4G,
Sheth et al., 2010). Also shown the bar-major-axis surface density profiles
for simulation SD1 at $7\mbox{$\>{\rm Gyr}$}$ (dotted black) and at
$10\mbox{$\>{\rm Gyr}$}$ (dash-dotted blue) – see Sec. 2 and Fig. 3. NGC 3681
has a pure-exponential profile, while NGC 4340 and model SD1 have shoulders
(red and orange), as detected in Paper I.

In principle, profiles with shoulders might simply signal the presence of a
strong and long bar. However, exploring a diverse set of simulations, Paper I
found only a weak dependence on bar strength or bar length, while Erwin et al.
(2023) found no evidence of such correlation in their observational sample.
Importantly, Paper I demonstrated that shoulders tend to be present in
galaxies with growing bars, although their results suggest that this is a
necessary but not a sufficient condition. Furthermore, Paper I found that in
secularly evolving (non-buckling) bars, once formed, shoulders tend to
persist. On the other hand, in bars that buckle, shoulders tend to appear
right after the buckling, and to be erased (at least temporarily) if a second
buckling happens. Finally, Paper I suggested that shoulders are mostly
supported by looped $x_{1}$ orbits.

In this paper, we further investigate the shoulders in bar major-axis profiles
of galaxies, focusing on their orbital support and the time evolution of these
orbits. The main questions we aim to answer are: how and why do some bars
develop shoulders? Which orbits support the shoulders? Why do bars need to
grow in order to retain the shoulders once formed? What is behind the role of
buckling events in creating or erasing shoulders?

We answer these questions by selecting particles from $N$-body simulations,
integrating their orbits in the “frozen” potentials of different snapshots,
and performing frequency analyses of these orbits. The $N$-body simulations
are described in Sec. 2. In Sec. 3, we explain the procedures for orbit
integration and frequency analyses. The resulting frequency maps and stacked
density profiles for different orbital groups, as well as their time
evolution, are presented in Sec. 4. We demonstrate the somewhat surprising
result that, once shoulders have formed, vertical resonances act as the main
culprit for erasing them, shedding light on the suggestion of Paper I that, in
order for shoulders to be present for long times, the bar needs to be growing
and trapping new stars. In Sec. 5 we discuss our results, and the conclusions
are summarized in Sec. 6.

##  2 Simulations

Here we describe the $N$-body simulations explored in this paper. These
simulations are chosen to represent a diverse subset of the models analysed in
Paper I. Briefly, one model (SD1, our fiducial model) has clear and persistent
shoulders appearing around 4Gyr; one model (SD1S) does not develop shoulders;
model HG1 is a fully self-consistent star-forming $N$-body+SPH simulation with
weak shoulders; and Model 4 has a bar that buckles around $3.8\mbox{$\>{\rm
Gyr}$}$, developing persistent shoulders after that. These simulations are
described in detail in Paper I, but we summarise them here – see Table 1.

Model SD1 is a pure $N-$body model, built using a modified version of
GalactICS (Kuijken & Dubinski, 1995; Widrow & Dubinski, 2005) that generates
discs using a Sérsic surface density profile:

$\Sigma_{d}=\frac{M_{d}}{2\pi
nR_{d}^{2}\Gamma(2n)}e^{-(R/R_{d})^{1/n}}\mathrm{sech}^{2}(z/z_{d})~{},$ (1)

where $M_{d}$ is the disc mass, $R_{d}$ is the scale length, $n$ is the Sérsic
index, and $\Gamma$ is the gamma function. The Sérsic disc has
$M_{d}=5.74\times 10^{10}\>{\rm M_{\odot}}$, $R_{d}=0.265\mbox{$\>{\rm
kpc}$}$, $z_{d}=0.25\mbox{$\>{\rm kpc}$}$ and $n=2.05$. The model has a
Hernquist dark matter (DM) halo, defined with GalactICS parameters
$\sigma_{h}=550$ kms-1, $R_{h}=30$ kpc, $\alpha=1$ and $\beta=4$. It has
$5\times 10^{6}$ DM particles and $4.4\times 10^{6}$ disc particles. In this,
and in the other pure $N$-body simulations described below, the softening
length is $\epsilon=50$ pc and $\epsilon=100$ pc for star- and DM-particles,
respectively.

Model SD1S suppresses most secular bar growth by setting the halo of model SD1
in full prograde rotation (Debattista & Sellwood, 2000; Long et al., 2014;
Collier et al., 2018). This results in a large halo spin ($\lambda\simeq
0.091$), rare in cosmological simulations (Bullock et al., 2001). Our goal is
to contrast this model with its counterpart SD1 whose bar is growing, so we
are not concerned by its fully rotating halo being unrealistic.

HG1 is a star-forming simulation described in Cole et al. (2014), Gardner et
al. (2014) and Debattista et al. (2017). It starts out with a gas corona but
no stars, hence the evolution is totally self-consistent. It is evolved with
the smooth particle hydrodynamics code GASOLINE (Wadsley et al., 2004). The
simulation uses high resolution ($\epsilon=50$ pc for the gas, and gas
particles with initial mass $2.7\times 10^{4}\>{\rm M_{\odot}}$) and after
being evolved for 10 Gyr has $\sim 1.1\times 10^{7}$ stellar particles, and a
total stellar mass of $\sim 6.5\times 10^{10}\mathrm{M}\odot$. It uses the gas
cooling, star formation and stellar feedback prescriptions of Stinson et al.
(2006). The DM halo consists of $5\times 10^{6}$ particles with $\epsilon=103$
pc, where $4.5\times 10^{6}$ of them have mass $8.5\times 10^{4}\>{\rm
M_{\odot}}$, and $1.7\times 10^{6}\>{\rm M_{\odot}}$ for the remainder.

Finally, Model 4 is a disc $N$-body model (no gas or star formation) – see
also Debattista et al. (2020). It is a baryon-dominated disc Milky Way-like
model, again set up using GalactICS. It has a Navarro-Frenk-White DM halo
(Navarro et al., 1996), and an exponential disc with an isothermal vertical
profile. The disc has a central radial velocity dispersion of 128 kms-1, a
disc scale radius of 2.4 kpc and a scale height of 150 pc. It has $6\times
10^{6}$ particles in the disc and $4\times 10^{6}$ in the halo.

All $N$-body models were evolved using PKDGRAV (Stadel, 2001) for 10 Gyr. In
what follows, we take model SD1 as our fiducial model, because of its clear
trends illustrating the main results of this paper, although model HG1 may be
more realistic.

###  2.1 Simulation global properties

We use pynbody (Pontzen et al., 2013) to read, center and align the
simulations such that, at each snapshot analysed, the bar is initially along
the $x$-axis. We estimate the gravitational potentials $\Phi$ and integrate
orbits using agama (Vasiliev, 2019). Smooth potentials due to star- and gas-
particles (when present) are separately modelled as triaxial distributions and
estimated with the “cylspline” potential type, while the halo contribution is
assumed to be axisymmetric and is estimated with the “multipole” type. For
each snapshot analysed, we then estimate the total potential.

We compute the bar length (radius) $a_{\rm bar}$ as the mean of the
cylindrical radius at which the amplitude of the $m=2$ Fourier moment of the
$(x,y)$-plane surface density distribution reaches half its maximum value
after its peak, and the radius at which the phase of the $m=2$ component
deviates from a constant by more than $10^{\circ}$ – see Paper I for further
details. We fit a straight line to the bar length as a function of time and
use the fitted function to represent $a_{\rm bar}$ at each snapshot. Both
estimated and fitted bar lengths are shown in the upper panel of Fig. 2 as
thin and thick lines, respectively, and we see that a straight line is a
reasonable fit.

Figure 2: Evolution of bar lengths (top) and pattern speeds (bottom) for our
simulations. In the upper panel, thick lines show the linear fit to the bar
length evolution (thin), which provides $a_{\rm bar}$ used in our analysis. In
all models the bars slow down (but less significantly for SD1S and Model 4),
and grow (but less significantly for models SD1S and HG1). For HG1, the values
are shifted down by $100\mbox{$\>{\rm km/s/kpc}$}$.

The bar pattern speeds $\>\Omega_{\mathrm{P}}$ are computed as the time
derivative of the phase angle of the $m=2$ Fourier mode of star-particles, as
done by Debattista & Sellwood (2000). We use 21 snapshots spaced by 5 Myr
around the snapshots of interest. We then fit a straight line to the bar phase
angle as a function of time in each of these time intervals, where the phase
is assumed to increase monotonically with time. The pattern speeds are shown
in Fig. 2 (bottom panel) as functions of time. As often observed in
simulations (e.g. Hernquist & Weinberg, 1992; Debattista & Sellwood, 1998;
O’Neill & Dubinski, 2003; Athanassoula, 2003), the bars slow down
significantly over the time span of the simulation – see also Fragkoudi et al.
(2021) for a recent study of cosmological simulations and Hamilton et al.
(2022) for a recent theoretical treatment.

The bar in the Model SD1 starts with a slightly higher $\>\Omega_{\mathrm{P}}$
in comparison to Model 4, but the pattern speeds in both models evolve to very
similar, almost constant, values towards the end of the simulations. The model
HG1 has a very centrally concentrated mass distribution, with a steep rotation
curve in the central parts. Thus it has a very fast bar, which decelerates
quickly without increasing in length significantly. However, these differences
do not seem to alter the general trends of the shoulders we investigate.

Figure 3: Edge-on (top row) and face-on (middle row) normalised surface
density maps for the four simulations at the selection times. Shaded
rectangles show the shoulder regions $|y|/\mbox{$a_{\rm bar}$}<0.145$ and
$x$-edges identified in Paper I (except for SD1S which does not have detected
shoulders), where particles are selected. The bottom row shows the surface
density along the stripe $|y|/\mbox{$a_{\rm bar}$}<0.145$, with the red
snippets showing the shoulders (where detected).

Fig. 3 shows edge-on (top row) and face-on (middle row) views of the four
models at the times when particles are selected for orbit integration. The
coordinates are normalized by the respective bar lengths. Colors represent
surface densities in log-scale. The $x-y$ ($x-z$) contours show 9 (7) equally-
spaced percentiles, from 85 to 99 (84 to 99). The edge-on view shows clear
signs of a BP-bulge in all models, although the “X” shape is weaker in models
SD1S and HG1.

In the face-on view, we see the bars along the $x$-axis, and hints of density
excess at their outskirts (except for SD1S, and barely for HG1) shown by the
contours. The shaded rectangles represent the shoulder region, defined as
$|y|/\mbox{$a_{\rm bar}$}<0.145$ and $x$-edges detected with the algorithm of
Paper I. Model SD1S has no shoulders, and we define its shoulder region using
the same fractions of the bar length as those of SD1. In all models with
shoulders, these are detected in the transition region between the thick and
thin parts of the bar. The bottom row shows the surface density along the
stripe defined by $|y|/\mbox{$a_{\rm bar}$}<0.145$ in each model, with red
snippets showing where shoulders are detected.

Table 1 shows the number of star-particles selected in the shoulder region,
the bar length $a_{\rm bar}$ and bar speed parameter
$\mathcal{R}=R_{\mathrm{CR}}/\mbox{$a_{\rm bar}$}$, where $R_{\mathrm{CR}}$ is
the co-rotation radius (Debattista & Sellwood, 2000) at the selection times,
and whether the model has gas or a strongly buckling bar (and when it
buckles). Models SD1S and HG1 have fast bars ($\mathcal{R}\lesssim 1.4$) at
these snapshots, while SD1 and Model 4 have slow bars ($\mathcal{R}\gtrsim
1.4$).

Table 1: Selection times ($t_{s}$), bar speed parameter $\mathcal{R}=R_{\mathrm{CR}}/\mbox{$a_{\rm bar}$}$, bar length $a_{\rm bar}$, number of selected star-particles and whether the model has gas or a strongly buckling bar. Model SD1 and model 4 have reasonably strong shoulders. Model HG1 has weak shoulders and SD1S has a barely evolving bar and no shoulders. The simulation naming convention is the same as in Paper I. Model | $t_{s}$ | $\mathcal{R}$ | $\mbox{$a_{\rm bar}$}/\mbox{$\>{\rm kpc}$}$ | $N$ | gas? | Buckling?
---|---|---|---|---|---|---
SD1 | $7\mbox{$\>{\rm Gyr}$}$ | 1.99 | 6.58 | 106032 | N | N
SD1S | $7\mbox{$\>{\rm Gyr}$}$ | 1.31 | 3.78 | 46973 | N | N
HG1 | $8\mbox{$\>{\rm Gyr}$}$ | 1.15 | 2.46 | 177426 | Y | N
4 | $5\mbox{$\>{\rm Gyr}$}$ | 1.60 | 6.85 | 136488 | N | Y, 3.8 Gyr

###  2.2 Epicyclic frequencies, resonance names and selection region

Fig. 4 shows the rotation curve $\Omega(R)$ (solid black), and the curves
${\Omega(R)-\kappa(R)/2}$ (dashed) and ${\Omega(R)-\kappa(R)/4}$ (dotted),
where

$\Omega(R)\equiv\sqrt{\frac{1}{R}\frac{\partial\Phi}{\partial R}},$ (2)

and

$\kappa(R)\equiv\sqrt{\frac{\partial^{2}\Phi}{\partial
R^{2}}+\frac{3}{R}\frac{\partial\Phi}{\partial R}}$ (3)

is the radial epicyclic frequency (Binney & Tremaine, 2008), evaluated along
the bar major-axis, at the times when star-particles are selected. The times
are chosen when clear and stable shoulders are present (except for SD1S which
does not have shoulders), and which are far from buckling events (in Model 4).

Figure 4: Rotation curves $\Omega(R)$ (solid black), $\Omega(R)-\kappa(R)/2$
(dashed) and $\Omega(R)-\kappa(R)/4$ (dotted) for the different simulations at
the time when star-particles are selected. Gray areas show the average radial
interval where shoulders are detected (SD1S has no shoulders, and we define
this region using the same fractions of the bar length as those of SD1),
within which we select particles in the simulations. Horizontal red lines show
the bar pattern speeds, with the thick part representing the bar lengths, at
the respective times. Model SD1 and Model 4 have slow bars
($\mathcal{R}\gtrsim 1.4$), while models SD1S and HG1 have fast bars. The
shoulder region is around the Ultra Harmonic radius where
$\Omega-\kappa/4=\mbox{$\>\Omega_{\mathrm{P}}$}$.

The bar pattern speeds at these times are shown as thin red lines, with the
thick part representing the bar length. Gray rectangles show the area within
the mean shoulder $x$-edges, as detected in Paper I, within which particles
are selected. For simulation SD1S, we adopt the same fraction of the bar
length as that of the shoulder extent of SD1 at the same time.

It is interesting to note that for all simulations with shoulders, these start
outside the outer Inner Lindblad radius (where $\Omega-\kappa/2=\Omega_{P}$).
This is the region where the (main bar-supporting) $x_{1}$ orbits can start
developing loops at their ends (Contopoulos & Papayannopoulos, 1980; Binney &
Tremaine, 2008), and this is the first hint of the importance of these orbits
to shoulders, as we confirm below. On the other hand, this might suggest that
the frequency ratio
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$,
where $\>\Omega_{\varphi}$ and $\Omega_{R}$ are the azimuthal and radial
frequencies for each orbit, is not relevant for the shoulders. Conversely,
Fig. 4 shows that the shoulders typically appear near the Ultra Harmonic
radius (where $\Omega-\kappa/4=\Omega_{P}$), which might suggest an important
role of the resonance
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/4$.
However, our frequency analysis (Sec. 4) reveals a fundamental role for the
resonance
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$
in supporting the shoulders, and no important role for the resonance
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/4$.

The resonance
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$
is the analog of the classical Inner Lindblad Resonance, but applies to orbits
of arbitrary eccentricity. Keeping in mind that resonances are defined in the
space of frequencies (rather than physical space), in this paper we will use
actual frequency ratios to define the main resonances. For instance, the
resonance
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$
will be called the Inner Lindblad Resonance (ILR), and the resonance
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/4$
will be called Ultra Harmonic Resonance (UHR). When identifying the
approximate locations of these resonances, based on the epicyclic frequencies,
we will refer to their radii, e.g. the Inner Lindblad radius and the Ultra
Harmonic radius – see Athanassoula (2003) for a discussion on the pros and
cons of different ILR definitions.

Figure 5: Effective potential along the bar major-axis (thick black solid) and
Jacobi integral (points) of the particles selected in the shoulder region from
the different simulations, at the selection time. Also shown the effective
potential values at L1 (co-rotation, thin solid), at the Ultra Harmonic radius
(dotted) and outer Inner Lindblad radius (dashed).

With the bar aligned with the $x$-axis, we select star-particles within
$|y|/\mbox{$a_{\rm bar}$}<0.145$ and in the $x$-intervals where shoulders are
identified (with no cut in $z$) – see Figs. 3 and 4. Since the detected
shoulders are not exactly symmetric around $x=0$, we define the “shoulder
region” averaging their edges on the $x<0$ and $x>0$ sides. This selection is
convenient for identifying orbits supporting the shoulder, but also because,
being in the bar outskirts, this region is in the transition to the disc, and
particles supporting both structures are selected. Table 1 shows the total
number of star-particles selected in each model.

Fig. 5 illustrates the selection (shoulder) region in terms of the effective
potential in the rotating frame,

$\mbox{$\>\Phi_{\rm
eff}$}=\Phi-\frac{1}{2}\mbox{$\>\Omega_{\mathrm{P}}$}^{2}x^{2},$ (4)

evaluated along the bar major-axis $x$. The points show the Jacobi integral

$\mbox{$\>E_{\rm J}$}=E-\mbox{$\>\Omega_{\mathrm{P}}$}L_{z},$ (5)

which is conserved in the frame rotating with angular speed
$\>\Omega_{\mathrm{P}}$, where $E$ is the energy in the inertial frame and
$L_{z}$ is the angular momentum. These points are color-coded by the
eccentricity $\mbox{$\>{\rm
ecc}$}\equiv(r_{\mathrm{apo}}-r_{\mathrm{per}})/(r_{\mathrm{apo}}+r_{\mathrm{per}})$,
where $r_{\mathrm{apo}}$ and $r_{\mathrm{per}}$ are the apo- and pericenter
radii. These quantities are calculated after the orbit integration explained
in Sec. 3.1. The horizontal lines show the effective potential evaluated at
the co-rotation radius, $R_{\mathrm{RC}}$ (solid thin), at the Ultra Harmonic
radius (dotted), and at the outer Inner Lindblad radius, $R_{\mathrm{ILR}}$
(dashed). In agreement with Fig. 4, the selected particles have
$E_{J}>\mbox{$\>\Phi_{\rm eff}$}(R_{\mathrm{ILR}})$. We note that for the
simulations with shoulders the most bound orbits have high eccentricities
(thus small $r_{\mathrm{per}}$) and these orbits extend right up to the Ultra
Harmonic radius. Besides that, there are orbits with $E_{J}>\mbox{$\>\Phi_{\rm
eff}$}(R_{\mathrm{CR}})$ and reasonably large eccentricities, which can visit
both the disc and the bar regions – see Sec. 4.1.3.

##  3 Orbits

The coordinates of particles selected at each snapshot set the initial
conditions for orbit integration in the frozen potential rotating with the
corresponding pattern speed (see Athanassoula, 2013). This means that in the
rotating frame the potential is static for the whole integration time. We then
estimate the frequencies of motion in the three directions (as detailed
below), producing a frequency map for each snapshot analysed.

###  3.1 Orbit integration

For particles selected at a time $t_{s}$, we integrate orbits in the frozen
potentials of five snapshots, $[t_{s}-2,t_{s}-1,t_{s},t_{s}+1,t_{s}+2]$ Gyr,
with a relative error tolerance of $10^{-15}$ for each coordinate, which is
the best we could achieve for all particles at all snapshots. When plotting
orbits, we use coordinates in the bar rotating frame, but the frequency
analysis uses coordinates in the inertial frame, as detailed below.

We integrate orbits for a (pseudo-)time of $5\mbox{$\>{\rm Gyr}$}$ (storing
$10^{5}$ points per orbit), which is a good compromise between two limiting
time-scales: a lower limit of, at least, tens of orbital periods for precise
frequency estimates (Valluri & Merritt, 1998), and an upper limit given by the
predictability horizon for chaotic orbits with finite precision integration
(e.g. Roy, 1991). Note that this is not the real time of the simulation, but a
parameter in the determination of instantaneous orbital properties in the
frozen potential. This integration time represents at least $20\times
T_{\mathrm{circ}}(E)$, where $T_{\mathrm{circ}}(E)$ is the period of a
circular orbit with energy $E$, for almost all particles, and this is
typically the case for all selections. Furthermore, among the cylindrical
coordinates $(R,\varphi,z)$, the azimuthal period is typically the longest, so
the integration time is also large enough for the other coordinates.

###  3.2 Frequency analysis

We employ the Numerical Analysis of Fundamental Frequencies (NAFF) algorithm
originally proposed by Laskar et al. (1992) and further developed by Valluri &
Merritt (1998); Valluri et al. (2010) to extract the leading frequencies of
orbits. Regular orbits are quasi-periodic, so the time-series of each
coordinate $q$ can be approximated by

$q(t)=\sum_{k=1}^{k_{\mathrm{max}}}a_{k}e^{i\omega_{k}t},$ (6)

where $a_{k}$ is a complex amplitude, and $\omega_{k}$ is an integer
combination of the fundamental frequencies $\Omega_{1}$, $\Omega_{2}$ and
$\Omega_{3}$ (time-derivatives of the canonical angle variables). The power
spectrum of a Discrete Fourier Transform (DFT) of a regular orbit has discrete
peaks at frequencies $\omega_{k}$ (spectral lines, see Binney & Spergel,
1982). Typically, although not always, the frequency with the largest
amplitude (hereafter the leading frequency) corresponds to the fundamental
frequency.

In non-integrable potentials, quasi-periodicity cannot be assumed a priori,
which might hinder the extraction of meaningful fundamental frequencies.
However, most realistic galactic potentials, although non-integrable, still
host a large number of regular or weakly chaotic orbits, a picture (formally)
supported by the KAM theorems Kolmogorov (1954); Arnol’d (1963); Moser (1962).
Additionally, in non-integrable potentials resonance trapping significantly
affects the system’s dynamical evolution. Once the fundamental frequencies are
identified, resonances stand out in a plot showing ratios of the frequencies,
i.e. a frequency map (Laskar, 1990, 1999), such as the one shown in Fig. 6.

We have written a new pure-python implementation of the NAFF algorithm, which
we dub naif111The documentation and installation instructions can be accessed
at http://naif.readthedocs.io/en/latest/index.html., based on the Fortran
implementation of Valluri & Merritt (1998). A summary of the algorithm is
presented in Appendix A. Variants of this method have been developed and
applied in different contexts (e.g. Binney & Spergel, 1982; Carpintero &
Aguilar, 1998; Valluri et al., 2016; Price-Whelan et al., 2016; Beraldo e
Silva et al., 2019; Dodd et al., 2022; Lucey et al., 2022). The analysis
involves important decisions on implementation details, such as the coordinate
system (e.g. Cartesian or cylindrical), the use of real or complex time-
series, the reference frame, the window function etc. In the study of bars
(e.g. Portail et al., 2015; Valluri et al., 2016; Gajda et al., 2016; Parul et
al., 2020; Smirnov et al., 2021; Sellwood & Gerhard, 2020), it is common to
adopt Cartesian coordinates (or mixed Cartesian and cylindrical) in the bar
rotating frame, where the frequency of oscillation along the major axis,
$\Omega_{x}$, is used as a proxy for the azimuthal frequency
$\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$}$. Furthermore,
spatial coordinates as (real) time-series are often used. See Sec. 5.1 for a
discussion on different choices and techniques.

To get some guidance in these decisions, and to test our algorithm, in
Appendix B we compare our frequency estimates with analytical expressions for
orbits in the isochrone potential. We compare the use of real and complex
time-series (Appendix A), and show that the radial frequency $\Omega_{r}$ is
estimated with relative errors $|\delta\Omega_{r}/\Omega_{r}|\approx 10^{-6}$
in both cases. For the azimuthal frequency $\>\Omega_{\varphi}$, we show that
the complex time-series
$f_{\varphi}=\sqrt{2|L_{z}|}(\cos\varphi+i\sin\varphi)$ (Papaphilippou &
Laskar, 1996, 1998), where $L_{z}$ is the $z$-component of the angular
momentum, does offer advantages, compared with $f_{\varphi}=\varphi$: it
improves the accuracy and it does not require setting
$\operatorname{sign}(\mbox{$\>\Omega_{\varphi}$})=\operatorname{sign}(L_{z})$,
which is the case for the real time-series. After excluding the precessing
frequency $\Omega_{r}-\mbox{$\>\Omega_{\varphi}$}$, which is the leading one
for orbits with large apocenter, we estimate $\>\Omega_{\varphi}$ with
relative errors
$|\delta\mbox{$\>\Omega_{\varphi}$}/\mbox{$\>\Omega_{\varphi}$}|\lesssim
10^{-3}$ for almost all orbits.

All frequency estimates in this paper use the Hanning window, i.e. $p=1$ in
Eq. A2. In the following analysis, for an easier identification of important
known resonances, we use cylindrical coordinates ($R,\varphi,z$) stored in the
inertial frame of the galaxy, rather than the bar rotating frame. This allows
identification of the co-rotation resonance, which would require an infinite
integration time in the rotating frame – see Appendix C. Based on the
discussions above, we use as input in the frequency analysis the complex time-
series:

$\begin{array}[]{l}f_{R}=R+iv_{R}\\\
f_{\varphi}=\sqrt{2|L_{z}|}(\cos\varphi+i\sin\varphi)\\\
f_{z}=z+iv_{z},\end{array}$ (7)

where $v_{R}$ and $v_{z}$ are the radial and vertical velocities,
respectively. For the radial and vertical components we take the absolute
value of the leading frequency, while we retain the original sign for
$\>\Omega_{\varphi}$. Since our cuts do not select particles with
significantly large apocenter distances, we expect the number of orbits with
misidentified fundamental frequencies to be negligible.

##  4 Results

###  4.1 Frequency maps: overview

Having demonstrated the accuracy of our numerical procedure, we move on to the
analysis of orbits in the frozen potentials of our $N$-body simulations. Fig.
6 shows the frequency map obtained for star-particles selected in the shoulder
region of model SD1 at $7\mbox{$\>{\rm Gyr}$}$ (see Figs. 3 and 4). The upper
panel is color-coded by the normalized Jacobi integral (see Fig. 5)

$\frac{\mbox{$\>E_{\rm J}$}-\mbox{$\>\Phi_{\rm eff}$}(0)}{\mbox{$\>\Phi_{\rm
eff}$}(R_{\mathrm{CR}})-\mbox{$\>\Phi_{\rm eff}$}(0)},$ (8)

where $\Phi_{\mathrm{eff}}$ and $\>E_{\rm J}$ are defined in Eqs. (4) and (5),
respectively, and $R_{\mathrm{CR}}$ is the co-rotation radius. The bottom
panel is color-coded by the eccentricity, and numbered stars indicate the
example orbits presented in Figs. 7, 9 and 10.

Figure 6: Frequency map for the shoulder region of model SD1 at
$7\mbox{$\>{\rm Gyr}$}$, with dashed lines showing prominent resonances. We
indicate the Inner Lindblad Resonance (ILR), the vertical ILR (vILR) and the
3D resonance
$\Omega_{R}+4(\Omega_{\varphi}-\mbox{$\>\Omega_{\mathrm{P}}$})-2\Omega_{z}=0$.
The upper panel is color-coded by the normalized Jacobi integral, Eq. 8. The
upper-left corner is dominated by typical orbits in the central region, and
moving towards the lower-right corner, orbits get more disc-like. The vILR
cloud is the high-eccentricity region around the vILR. Numbered stars are the
sample orbits shown in Figs. 7, 9 and 10.

Given the typical shape of rotation curves (cf. Fig. 4), we can anticipate
that orbits with high
$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$ in
Fig. 6 are typical of the galaxy central regions. Moving down in Fig. 6 we go
towards the end of the bar. For fast bars, this happens near the co-rotation,
at $\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$}=0$, and
transitioning to the disc we have
${\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$}\lesssim 0}$,
corresponding to retrograde stars in the bar frame. Evaluating the functions
$\Omega(R)$ and $\kappa(R)$, Eqs. (2)-(3), for this model at the average outer
shoulder edge (the outermost region of our selection), we obtain
$(\Omega-\mbox{$\>\Omega_{\mathrm{P}}$})/\kappa\approx 0.25$, in broad
agreement with the lower values ($\approx 0.15$) on the $y$-axis of Fig. 6,
and with the fact that the shoulder outer edge in this model matches the Ultra
Harmonic radius – see Fig. 5.

The ratio $\>\Omega_{z}/\Omega_{R}$ is a proxy for the vertical thinness of
the orbit: it is significantly larger than unity for typical orbits in the
disc (in the Solar neighborhood, $\mbox{$\>\Omega_{z}/\Omega_{R}$}\approx 2$,
Binney & Tremaine, 2008), and can be $\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim
1$ in thicker or more spheroidal regions. Thus, the upper-left corner in Fig.
6 is dominated by typical orbits that can visit the very central regions of
the galaxy, and orbits get more disc-like when we move towards the lower-right
corner, in agreement with the eccentricity trends seen in the bottom panel.
With the radial $\kappa(R)$, Eq. 3, and vertical
$\nu(R)\equiv\sqrt{\partial^{2}\Phi/\partial z^{2}}$ epicyclic frequencies
evaluated at the outer shoulder edge, we get $\nu/\kappa\approx 1.70$, also in
broad agreement with the rightmost limit of Fig. 6 (neglecting orbits at the
ILR, which are well inside the bar region and thus less well approximated by
the linear expressions).

We identify the main regions in the frequency map that we explore in detail
below. The horizontal line at
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$
characterizes the $x_{1}$ orbits, which are known to be the backbone orbits of
bars (Contopoulos & Grosbol, 1989; Skokos et al., 2002; Athanassoula, 2003).
In principle, $x_{2}$ orbits, which are aligned with the bar minor-axis, can
also be present along this resonance. However, these orbits are only expected
within the Inner Lindblad radius (or between two Inner Lindblad radii, see van
Albada & Sanders, 1982) and, as Fig. 4 shows, we select particles outside this
region.

As already mentioned in Sec. 2.1, we refer to the frequency ratio
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$
as the ILR (see e.g. Athanassoula, 2003; Weinberg & Katz, 2007). We emphasize
that this differs from the classical definition of the ILR based on the
epicyclic frequencies, Eqs. (2)-(3), or on the existence of $x_{2}$ orbits
(van Albada & Sanders, 1982), but can be seen as their generalizations, in the
frequency space, to orbits of arbitrary eccentricities. The UHR,
${\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/4}$,
is also well populated in Fig. 6. Also highlighted are the vertical ILR
(vILR), where
$(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{z}=1/2$;
the resonance $\mbox{$\>\Omega_{z}/\Omega_{R}$}=1$ and the 3D resonance
${\Omega_{R}+4(\Omega_{\varphi}-\mbox{$\>\Omega_{\mathrm{P}}$})-2\Omega_{z}=0}$,
denoted by (1,4,-2).

We call the cloud of highly eccentric orbits that are to the left and above
the vILR the vILR cloud. We also note two other diffuse populations more
clearly identified in the bottom panel: a large blue/green region with low
eccentricities and permeated by several resonances (the region of orbits 6 and
9), and an orange stripe below the ILR, with $0.65\lesssim\mbox{$\>{\rm
ecc}$}\lesssim 0.85$ and
$0.4\lesssim\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}\lesssim
0.45$ (the region of orbit 8).

We now illustrate the orbital morphologies in these different regions, noting
that Fig. 6 does not represent the whole bar region, but only orbits selected
at its outermost parts at a particular time. The whole menu of orbits in bars
has been extensively studied in other works and is beyond the scope of this
paper (see e.g. Athanassoula et al., 1983; Sparke & Sellwood, 1987; Combes et
al., 1990; Skokos et al., 2002; Patsis et al., 2002; Valluri et al., 2016;
Patsis & Athanassoula, 2019). For a review, see Sellwood & Wilkinson (1993).

####  4.1.1 Orbital shapes at the ILR

We first select orbits at the ILR by requiring
$|\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}-1/2|<0.0025$.
By inspection this cut is a good choice to ensure minimal contamination from
neighboring orbits. Unless otherwise stated, we use the same tolerance in
selections of other resonances. Fig. 6 suggests grouping the orbits at the ILR
into three groups, which we will refer to as vertically cool, warm and hot
orbits at the ILR (in some instances we will simply refer to cool, warm and
hot orbits for conciseness). Vertically cool orbits are those at the ILR and
with $\mbox{$\>\Omega_{z}/\Omega_{R}$}>3/2$. They are the least bound and are
generally less eccentric. A typical orbit from this group is shown in the
right column (first and third rows) of Fig. 7. It has an elliptical shape that
is elongated in the direction of the major axis, and it is vertically very
thin (see third row).

Figure 7: First and third rows: three typical orbits at the ILR selected at
different $\Omega_{z}/\Omega_{R}$-ranges from simulation SD1 at
$7\mbox{$\>{\rm Gyr}$}$. Orbits are color-coded by time, with dots showing the
initial and final locations. Numbered stars identify these orbits in Fig. 6.
Second and fourth rows: Stacked surface density maps of all orbits at the ILR
in the respective $\Omega_{z}/\Omega_{R}$ ranges. The orbits and stacked
density maps are shown in the bar rotating frame, and the thick gray lines
show the bar size. Thin gray lines are the same contours shown in Fig. 3 for
the total mass. The vertically warm orbits
${(1<\mbox{$\>\Omega_{z}/\Omega_{R}$}<3/2)}$ produce a marked density peak in
the outskirts of the bar. A similar, but diluted, shape is also produced by
the hot orbits $(\mbox{$\>\Omega_{z}/\Omega_{R}$}<1)$.

Transitioning to the middle section of the ILR, the vertically warm orbits
have ${1<\mbox{$\>\Omega_{z}/\Omega_{R}$}<3/2}$. They are more bound (lower
$\>E_{\rm J}$) and eccentric ($\mbox{$\>{\rm ecc}$}\approx 0.85$). The upper-
middle panel in Fig. 7 shows a typical orbit in this region (orbit 2). These
orbits still avoid the center and have characteristic loops at their ends,
which have previously been identified in both analytical (e.g. Contopoulos,
1978; Papayannopoulos & Petrou, 1983) and numerical studies (e.g. Sparke &
Sellwood, 1987; Petersen et al., 2021). In particular, Contopoulos (1988)
found that these orbits typically appear near the Ultra Harmonic radius, which
corresponds to our selection region – see Fig. 4. We refer to these as
“looped” orbits.

Furthermore, for all simulations analyzed in this paper, we find that this
transition from elliptical-like to this looped shape occurs at
$\mbox{$\>\Omega_{z}/\Omega_{R}$}\approx 3/2$, suggesting some kind of
influence of vertical resonances. In fact, we note that the crossing of the
ILR with $\mbox{$\>\Omega_{z}/\Omega_{R}$}=3/2$ is the convergence point of
several 3D resonances in Fig. 6. In the next sections we demonstrate the
fundamental role of this looped morphology in building the shoulders in the
bar major-axis density profiles.

Finally, the vertically hot orbits are those at the ILR and with
${\mbox{$\>\Omega_{z}/\Omega_{R}$}<1}$ (see Fig. 6). They are slightly deeper
in the potential, and thus have slightly smaller $r_{\mathrm{apo}}$. Despite
this, they are significantly more eccentric ($\mbox{$\>{\rm ecc}$}\approx
0.98$, see also Fig. 5), because they can get very close to the center, as
illustrated in the upper-left panel in Fig. 7 (orbit 1).

Fig. 8 shows a normalized histogram of
$\log_{10}(r_{\mathrm{per}}/\mathrm{kpc})$ for cool, warm and hot orbits at
the ILR of the model SD1 at $7\mbox{$\>{\rm Gyr}$}$. It is clear that a
significant fraction of the hot orbits can get as close as 10-100 pc to the
center. Moreover, Fig. 7 shows that typical orbits in this region do not have
the coherent shapes observed before, scattering more in the $x-y$ plane (a
signature of their chaotic nature, as demonstrated in Sec. 4.2). Finally,
their vertical thickness is significantly larger than in orbits 2 and 3.

Figure 8: Normalized histogram of $\log_{10}(r_{\mathrm{per}}/\mathrm{kpc})$
for orbits along the ILR, in the different $\>\Omega_{z}/\Omega_{R}$
intervals, and for orbits at the vILR cloud – see main text. A large fraction
of the vertically hot orbits at the ILR (and the vILR cloud in general) can
get as close as $\sim 10-100$ pc from the center. The small secondary peaks
for cool and warm orbits at the ILR are due to cross-contamination of these
orbits near the edge $\mbox{$\>\Omega_{z}/\Omega_{R}$}=3/2$

.

The second and fourth rows in Fig. 7 show the orbit-averaged surface densities
of orbits at the ILR in each of the three groups in $\>\Omega_{z}/\Omega_{R}$.
These plots are 2D histograms built by stacking all points of a single orbit
(a time average) and all orbits (an ensemble average), normalized by the total
number of orbits in each region to emphasize characteristic shapes rather than
relative contributions of each group. We see that the sample orbits are
representative of the populations.

The thin gray lines in the second and fourth rows are the same face-on and
edge-on contours shown in Fig. 3 for the whole star-particle distribution. We
anticipate the role of the looped morphology of vertically warm orbits at the
ILR to produce the density excess in the bar outskirts. A similar morphology
is also present for the hot orbits at the ILR, but the density excess appears
more diluted – see Secs. 4.4 and 4.6. We also note the X-shape in the $x-z$
plane of the orbit-averaged density of the hot orbits, suggesting their
contribution to the BP-bulge visible in the contours.

####  4.1.2 Orbital shapes at vertical resonances

We now select orbits at the three vertical resonances highlighted in Fig. 6
and mentioned before: the vILR, the resonance
$\mbox{$\>\Omega_{z}/\Omega_{R}$}=1$, and the resonance (1,4,-2). In this
selection we exclude the intersections with the ILR and orbits above it by
requiring
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}-1/2<-0.01$.
Typical orbits (orbits 4, 5 and 6 in Fig. 6) are shown in the first and third
rows of Fig. 9, while the panels below show the stacked density maps of all
the orbits selected along these resonances. From right to left, we observe a
similar trend of orbits becoming more elongated, less regular, and vertically
thicker, with orbit 4 reaching all the way to the center. This orbit, at the
vILR, has the typical “banana” shape in the $x-z$ plane (Pfenniger & Friedli,
1991). Overall, orbits along these resonances imprint noticeable features in
the $x-z$ stacked density maps, but not as distinctly concentrated in the
$x-y$ plane as those shown in Fig. 7.

Figure 9: Similar to Fig. 7, but for typical orbits selected at three vertical
resonances (titles), and requiring
$(\Omega_{\varphi}-\Omega_{\mathrm{P}})/\Omega_{R}<1/2$. The numbered stars
indicate these orbits in Fig. 6.

####  4.1.3 Orbital shapes at the vILR cloud and other regions

Finally, we select orbits in three additional regions mentioned before: the
vILR cloud (with sample orbit 7); the stripe with intermediate eccentricities
(with sample orbit 8), requiring $0.65<\mbox{$\>{\rm ecc}$}<0.85$ and not
belonging to any of the previously selected regions; and the low-eccentricity
region (with sample orbit 9), requiring $\mbox{$\>{\rm ecc}$}<0.65$ and not
belonging to any of the previously selected regions.

Figure 10: Similar to Figs. 7-9, but for orbits in the vILR cloud and other
locations in the frequency map (see main text). The numbered stars indicate
these orbits in Fig. 6.

Fig. 10 shows the sample orbits and the orbit-averaged surface density maps
for these regions. From right to left, i.e. in order of decreasing
$\>\Omega_{z}/\Omega_{R}$, we note again the tendency of orbits to become more
elongated and less regular, finally approaching the center and becoming
significantly thicker. Orbit 9 circulates a little bit further out in the bar
(in comparison to orbits 6 and 3). Its orbital group has
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}\lesssim
0.4$ and these resemble disc orbits in the vicinity of the bar. Orbit 8 has a
peculiar shape, with a strong resemblance to the looped orbits at the ILR
(e.g. orbit 2), but oscillating between this and a close-to resonant
triangular shape. Its orbital group has
$0.4\lesssim\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}\lesssim
0.45$ and it produces an excess in the outskirts of the bar similar to that
produced by the looped orbits, but less prominent. These orbits are located
near the convergence of several 3D resonances and near the ILR and, as will be
demonstrated in Sec. 4.2, are chaotic. These orbits have
$0.4\lesssim\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}\lesssim
0.45$, similar to the location of a peak already observed by Martinez-
Valpuesta et al. (2006); Smirnov et al. (2021). Some of these orbits (but not
all) have $\mbox{$\>E_{\rm J}$}>\mbox{$\>\Phi_{\rm eff}$}(R_{\mathrm{RC}})$
and reasonably large eccentricities (see Fig. 5), allowing them to visit both
the bar and the disc region, as those identified in $N$-body simulations by
Sparke & Sellwood (1987) and Pfenniger & Friedli (1991). Fig. 11 is similar to
Fig. 5, but only shows points for the orbits in this group. We see that in
models SD1S and HG1 a large fraction of these orbits have $\mbox{$\>E_{\rm
J}$}>\mbox{$\>\Phi_{\rm eff}$}(R_{\mathrm{RC}})$.

Figure 11: Similar to Fig. 5, but only showing $\>E_{\rm J}$ (points) for the
particles selected in the group of orbit 8 in Fig. 6. These orbits typically
have
$0.4\lesssim\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}\lesssim
0.45$ and in some simulations a significant fraction of these orbits have
$\mbox{$\>E_{\rm J}$}>\mbox{$\>\Phi_{\rm eff}$}(R_{\mathrm{RC}})$, i.e. can
visit both the bar and the disc region, as those identified by Sparke &
Sellwood (1987) and Pfenniger & Friedli (1991)

.

An interesting aspect of Figs. 7, 9 and 10 is the similarity of the orbit-
averaged density maps in the left column of each figure, especially regarding
the X-shape in the $x-z$ projection supporting the BP-bulge (contours). All
these orbits are within the vILR cloud, now specifically defined as
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{z}$}-1/2>-0.01$,
i.e. including hot orbits at the ILR. Fig. 8 shows that approximately half of
these orbits can get as close as $10-100$ pc to the galactic center.

Having shown the typical spatial distribution of orbits selected in the
shoulder region, we now compare the frequency maps of the different
simulations.

###  4.2 Frequency maps: comparing simulations

Here we compare the frequency maps for orbits selected in the shoulder regions
for the four simulations (see Sec. 2), at the snapshots where particles are
selected (cf. Fig. 4). Fig. 12 shows the results for the model SD1, with the
same data as in Fig. 6, but with the upper panel color-coded by $z_{\rm max}$
and the bottom panel by the frequency-drift $\log\Delta\Omega$, where

$\Delta\Omega=\mathrm{max}_{i=1,2,3}\left|\frac{\Omega_{i}(T_{2})-\Omega_{i}(T_{1})}{\Omega_{i}(T_{1})}\right|,$
(9)

and $\Omega_{i}(T_{2})$ and $\Omega_{i}(T_{1})$ are the frequencies estimated
in two separate halves of the time-series (Valluri et al., 2010). Since
regular orbits conserve frequencies, a large $\log\Delta\Omega$ indicates a
chaotic orbit. The side-panels show histograms of the frequency ratios color-
coded by the mean $z_{\rm max}$ (or $\log\Delta\Omega$) in bins of width
0.0025.

Figure 12: Frequency map for orbits in the shoulder region of simulation SD1
at $7\mbox{$\>{\rm Gyr}$}$. Smaller panels show histograms (in log-scale) with
bin-widths of 0.0025. The upper panel is color-coded by the vertical
excursion, and the bottom panel by the frequency drift (high values indicating
chaotic orbits). The ILR,
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$,
is the most prominent resonance. Orbits in the vILR cloud, to the left and
above the vILR
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{z}$}=1/2$,
have significantly larger $z_{\rm max}$ and are strongly chaotic.

Figs. 13-15 show the equivalent frequency maps for the models SD1S, HG1 and
Model 4. The most noticeable feature in the four models is the well-populated
ILR, producing a very prominent peak in the histogram to the right (note the
logarithmic scale). This strong prominence of the ILR could be a selection
effect since we chose particles which are predominantly part of the bar.
However, in a random selection of $10^{6}$ star-particles in HG1 at
$10\mbox{$\>{\rm Gyr}$}$, we identify more than half of the particles at the
ILR; such a strong ILR has also been observed in previous works (e.g.
Athanassoula, 2003). In our selection, within $\pm 0.0025$ ($\pm 0.01$) from
the ILR we identify $19.8\%$ ($22.2\%$) of the orbits for model SD1, $4.5\%$
($4.5\%$) for SD1S, $13.5\%$ ($18.3\%$) for HG1, and $15.1\%$ ($18.1\%$) for
Model 4, meaning the ILR is well populated in all the models, except for model
SD1S.

Figure 13: Similar to Fig. 12, but for simulation SD1S (fast bar and no
detected shoulders). The most prominent resonances are the ILR, the co-
rotation and $\mbox{$\>\Omega_{z}/\Omega_{R}$}=1$. The ILR is not as strongly
populated as in Fig. 12, and there is no vILR cloud. Orbits with
$\mbox{$\>\Omega_{z}/\Omega_{R}$}<1$ still have larger $z_{\rm max}$ than the
rest, but not anywhere as large as in Fig. 12.

In the models SD1S and HG1, we identify the co-rotation resonance at
$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=0$.
This is in agreement with SD1S and HG1 having fast bars, and the selection
region, always in the outermost parts of the bar, being close to co-rotation
in these models– see Fig. 4. We emphasize that identifying the co-rotation is
only possible when using coordinates in the inertial frame – see Fig. 31 for a
comparison with the bar frame. We also note, in all simulations, big gaps
around the co-rotation and ILR, testifying to their influence: orbits that get
too close may be trapped or scattered by them.

Figure 14: Similar to Fig. 12, but for model HG1, which has gas and a fast
bar. We see several well populated resonances, particularly the ILR, the co-
rotation and the vertical resonances $\mbox{$\>\Omega_{z}/\Omega_{R}$}=1$ and
$\mbox{$\>\Omega_{z}/\Omega_{R}$}=3/2$. The vILR cloud is strongly populated
with chaotic orbits and relatively high $z_{\rm max}$, but not as high as in
the model SD1, Fig. 12.

Another remarkable feature in all models, except SD1S, is the vILR cloud, with
significantly large $z_{\rm max}$ and $\log\Delta\Omega$. In HG1, which is the
only model with gas, the $z_{\rm max}$ at the vILR cloud is smaller than in
the other models, although $\log\Delta\Omega$ is still comparably large. In a
certain sense, the vILR clearly separates the orbits that support the thin and
thick parts of the bar. The looped orbits at the (vertically warm) ILR mostly
support the thin part, and the orbits at the vILR cloud mostly support the
thick part – see also Fig. 16.

Finally, we see a plethora of additional resonances, particularly in HG1 (Fig.
14). Orbits at the main resonances (and away from any crossing resonances)
tend to have smaller $\log\Delta\Omega$, i.e. to be regular, as expected since
these resonances are stable. On the other hand, orbits in areas with multiple
strong resonances, such as the crossing of the ILR and the vILR, or on the
floor below the ILR (the group of orbit 8 in Fig. 6), tend to be chaotic, as
is generally the case at resonance crossings (Chirikov, 1979). These results
illustrate the richness of the dynamics in the region between the bar and the
disc (see e.g. Contopoulos, 1981), and also the power and precision of the
frequency map technique we employ.

Figure 15: Similar to Fig. 12, but for Model 4, whose bar buckles at
$4\mbox{$\>{\rm Gyr}$}$ and is slow at $5\mbox{$\>{\rm Gyr}$}$. This figure is
similar to that of model SD1, with $z_{\rm max}$ values in the vILR cloud
intermediate between those in SD1 and HG1, Figs. 12 and 14.

###  4.3 More on morphologies along the ILR

Figure 16: Top: vertical excursion versus $\>\Omega_{z}/\Omega_{R}$ for
particles selected in the shoulder region and at the ILR for all models.
Bottom: the fraction of time spent in retrograde motion (in the rotating
frame). The horizontal dashed line is at 0.5. The different morphologies
identified in Sec. 4.1 are clearly seen here (except for SD1S): cool orbits at
the ILR have low $z_{\rm max}$, low eccentricity and are almost exclusively
prograde; warm orbits at the ILR have slightly larger $z_{\rm max}$,
relatively high eccentricity and expend a significant time in retrograde
motion (a symptom of the strong loops at their ends); hot orbits at the ILR
have significantly larger $z_{\rm max}$ and eccentricity, and moderate
fractional time in retrograde motion. The looped morphology is not detected in
the shoulder-less model SD1S. We also note the overall lower $z_{\rm max}$ for
the SPH model HG1, in comparison to models SD1 and Model 4.

We explore in greater depth the different morphologies for orbits along the
ILR. Fig. 16 shows, for the different models, the vertical excursion $z_{\rm
max}$ (top panels) versus $\>\Omega_{z}/\Omega_{R}$ for orbits at the ILR
($|\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}-1/2|<0.01$)
at the selection times, color-coded by the eccentricity. All models (except
SD1S) show three distinct groups: the cool orbits
($\mbox{$\>\Omega_{z}/\Omega_{R}$}\gtrsim 3/2$) have low $z_{\rm max}$ and
low-eccentricity (blue cloud) – the elliptical orbits in Sec. 4.1. The warm
orbits ($1\lesssim\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 3/2$) have slightly
larger $z_{\rm max}$ and relatively high eccentricity (orange cloud) – the
looped orbits. And the hot orbits ($\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim
1$) have significantly larger $z_{\rm max}$ and eccentricity (red cloud) –
these are chaotic orbits within the vILR cloud. We also see some hot orbits at
the ILR (orbits with high $z_{\rm max}$) contaminating the
$\>\Omega_{z}/\Omega_{R}$-interval of the warm orbits, which is not surprising
given the crowd of orbits around the crossing of the ILR and the vILR – see
e.g. Fig. 6.

As suggested by Fig. 7 and demonstrated below, the loops at the ends of warm
orbits at the ILR play an important role in producing the shoulders along the
bar major-axis density profiles. As a consequence of these loops, these orbits
spend a considerable time in retrograde motion in the rotating frame of the
bar. The bottom panels of Fig. 16 show the fraction of time the orbits at the
ILR are in retrograde motion, i.e. with
$L_{z}<\mbox{$\>\Omega_{\mathrm{P}}$}R^{2}$. We see that this fraction is zero
for all cool orbits, except for some contamination of the warm orbits in model
HG1. In the range $1<\mbox{$\>\Omega_{z}/\Omega_{R}$}<3/2$ of the warm orbits,
the fraction of time in retrograde motion is typically high ($\approx
0.2-0.55$), a signature of the loops developed by these orbits. Finally, the
vertically hot orbits also have large fractions of time in retrograde motion,
but slightly smaller than the warm orbits.

Fig. 16 demonstrates the common $\>\Omega_{z}/\Omega_{R}$-dependence of the
different morphologies seen in Sec. 4.1 for the models with shoulders:
vertically warm orbits at the ILR have loops and spend a considerable time in
retrograde motion; and hot orbits at the ILR jump to large $z_{\rm max}$ and
also spend a considerable time in retrograde motion, but less than warm orbits
on average. The consequence of this for the shoulders is discussed in Sec.
4.6. Moreover, in agreement with the importance of the looped morphology to
the shoulders, as suggested in Paper I and confirmed in this work, we note
that this morphology is absent in the shoulder-less model SD1S.

###  4.4 The orbital support of shoulders

Here we quantify the contribution of the main orbital groups to the density
profiles in the shoulder region. Since the ILR is the main bar-supporting
resonance, and inspired by the morphological study in Secs. 4.1 and 4.3, we
investigate three groups: the vertically cool and warm orbits at the ILR, and
the vILR cloud (which includes the hot orbits at the ILR).

Figure 17: Surface density of particles in the shoulder region (gray shaded
areas) of model SD1 at $7\mbox{$\>{\rm Gyr}$}$, with coordinates normalized by
the bar length at $7\mbox{$\>{\rm Gyr}$}$. Panel a: instantaneous density
profile along the major axis ($|y|/\mbox{$a_{\rm bar}$}<0.145$) in the
simulation. The solid black curve shows the total profile, and colored curves
show contributions from the groups in panels b-d (see main text). The dashed
curve is the total profile after subtracting the warm orbits at the ILR
($1<\mbox{$\>\Omega_{z}/\Omega_{R}$}<3/2$), and the dotted one results from
also subtracting orbits at the vILR cloud. Panels b, c, d show the 2-D orbit-
averaged maps for the respective orbital groups, with colors normalized by the
total number of particles in the shoulders. Panel e shows the 1D density
profile, similar to panel a, but for the orbits, and restricting to
$|y|/\mbox{$a_{\rm bar}$}<0.145$. In this panel, the large contribution of
warm orbits to the profile in the shoulder region becomes more evident, with a
smaller but considerable contribution from the vILR cloud. The contribution
from the cool orbits (blue) is negligible.

Fig. 17 (panel a) shows the $N$-body snapshot surface density profile for SD1
at 7 Gyr, for particles in the shoulder region. The density profile is
calculated as the number of particles in $(x/\mbox{$a_{\rm bar}$})-$bins of
width 0.015 and within $|y|/\mbox{$a_{\rm bar}$}<0.145$, divided by $(\delta
x/\mbox{$a_{\rm bar}$})\,(\Delta y/\mbox{$a_{\rm bar}$})\,N$, where $\Delta
y=2\cdot 0.145$ and $N$ is the total number of selected particles. Since these
profiles are almost perfectly symmetric around $x=0$, for better
visualization, we present them averaged between the two sides $x<0$ and $x>0$.
The gray areas show the average shoulder region, where particles are selected.
The black solid curve shows the total density profile, and the colored solid
curves show the relative contribution of particles in the different orbital
groups, in the same colors as in the titles of panels b-d.

In this model the contributions from particles at the vILR cloud (red) and the
warm ones at the ILR (green) are similarly large. The black dashed curve shows
the total density profile excluding the contribution of the latter. Removing
the warm orbits significantly reduces the density profile in the shoulder
region, demonstrating their important role in the shoulders. The dotted curve
shows the profile after also subtracting particles at the vILR cloud, and we
see that they mostly contribute to the innermost parts of the shoulders.
Together, these two groups contribute so much to the density profile that
their subtraction produces significant deficits in it.

The picture provided by panel a in Fig. 17 is limited because it shows
particles currently in the shoulder region. To evaluate the global profile
produced by these particles, one can calculate the profile produced by them at
a later time, allowing the particles to phase-mix, as done in Paper I.
However, in the mean time the bar pattern speed and the orbits themselves may
have evolved into different ones. For a global and instantaneous view of the
density profiles, we stack the orbits integrated for these particles in the
frozen potential and plot the surface density in the $x-y$ plane in Fig. 17,
panels b, c, d. To emphasize each group’s relative contribution rather than
the orbital shapes, we now normalize colors by the total number of particles
in the shoulder region.

The orbit-averaged density of the orbits at the vILR cloud (panel b) resembles
that of looped orbits, although it is more diffuse in the $x-y$ plane, which
arises because of the chaotic nature of orbits at the vILR cloud – see bottom
panels of Figs. 12, 14, 15. Furthermore, its orbit-averaged density is less
extended along the major axis, peaking outside of the shoulder region. Panel c
shows the orbit-averaged density of the warm orbits at the ILR. We see that
this profile is more elongated and has a very pronounced peak in the shoulder
region. On the other hand, the cool orbits (panel d) do not represent any
visible excess in the shoulder region.

Panel e shows the density profiles of these stacked orbits (for points within
$|y|/\mbox{$a_{\rm bar}$}<0.145$) along the $x$-axis, analogous to the left
panel, and using the same color/style convention. For a comparison with panel
a, these profiles are further normalized by the number of points per orbit.
The smaller values, in comparison to panel a, are due to the “dilution” of
each particle along its orbit out of the selection region.

In these orbit-averaged profiles, the warm orbits at the ILR contribute much
more significantly to the shoulder because: i) they are more elongated than
those at the vILR cloud; ii) particles stay longer near apocenters; and iii)
the sharply defined loops at their ends make them spend even more time near
apocenters. This confirms the important contribution of the looped morphology
of the warm orbits at the ILR to the shoulders, but also demonstrates the
importance of the vILR cloud (compare the green and red curves in panels a and
e). Finally, the contribution of the cool orbits at the ILR is negligible also
in the orbit-averaged profiles.

Figure 18: As in Fig. 17, but for model SD1S, with no detected shoulders. In
this case, warm orbits at the ILR do not have loops and their contribution to
the density profile (solid green) is much smaller than in the model SD1, Fig.
17.

Figs. 18-20 show the corresponding profiles for models SD1S, HG1 and Model 4,
with the same color-normalization as in Fig. 17. Recall that model SD1S has no
shoulders, and that its “shoulder region” is set as the same fraction of the
bar length as in SD1 at the equivalent snapshot. The warm orbits at the ILR
(green) in this model do not develop loops, as already shown in Fig. 16, and
these orbits’ contribution to the density profile in the shoulder region is
negligible, as is the contribution from the vILR cloud. In fact, just a
handful of orbits populate the vILR cloud in this model (Fig. 13).

Fig. 19 shows analogous profiles for model HG1, which has weak shoulders – see
Paper I, Fig. 11. In this case, the contribution of the warm orbits is
relatively small, and its subtraction produces a profile (dashed) which is
close to the total one. On the other hand, in HG1 the vILR cloud orbits are
very numerous and produce an overall shape in the $x-y$ plane that resembles
the looped morphology, contributing significantly to the density in the
shoulder region. Fig. 20 shows the equivalent profiles for Model 4. This model
is intermediate between SD1 and HG1 in terms of bar length and pattern speed
evolution (see Fig. 2), and so are the relative contributions of warm orbits
at the ILR and of the vILR cloud. We also observe significant contributions
from the warm orbits and the vILR cloud, with the latter producing a more
diluted looped-like shape in the $x-y$ plane, in comparison with the warm
orbits.

Figure 19: As in Fig. 17, but for model HG1. The contribution from warm orbits
at the ILR (solid green) to the density profile is much smaller than in SD1.
However, in this case the number orbits at the vILR cloud is much larger, and
they produce an orbit-averaged density in the $x-y$ which resembles the one
produced by looped orbits. Figure 20: As in Fig. 17, but for Model 4. The
contributions from warm orbits at the ILR (green) and the vILR cloud (red) are
comparable, and a looped-like morphology is produced by both groups, as in the
model HG1, Fig. 19.

In summary, we have confirmed the important role of looped orbits to the
shoulders, as previously suggested in Paper I, and have shown that this
morphology is found in vertically warm orbits at the ILR. We also demonstrate
that, although individual orbits in the vILR cloud are chaotic and do not have
a well defined looped shape, they still produce an overall shape in the $x-y$
plane that can be seen as a diluted version of the shape produced by the
looped orbits. We note that the morphology of vILR cloud orbits in SD1 is the
most diffuse in the $x-y$ plane, and also has the largest $z_{\rm max}$, among
models with shoulders. Conversely, the vILR cloud of HG1 produces the least
diffuse shape in the $x-y$ plane, and the smallest $z_{\rm max}$, while Model
4 is intermediate between the two. This connection between $z_{\rm max}$ and
the looped shape is explored via the evolution of the frequency maps presented
in the next subsections.

###  4.5 Evolution of the frequency maps

Figure 21: Time evolution of the frequency map for the same set of star-
particles selected at $7\mbox{$\>{\rm Gyr}$}$ in the shoulder region of model
SD1. Points are color-coded by the vertical excursion $z_{\rm max}$. We see an
overall shift upwards in this plot mainly due to the bar slowing down, i.e.
decreasing $\>\Omega_{\mathrm{P}}$, and to the left, with orbits becoming
vertically thicker over time. Orbits crossing the vILR from right to left have
vertical motion dramatically amplified.

Having discussed the orbital support for shoulders, and having shown the
frequency maps at the time of selection in Sec. 4.2, we now show their time
evolution, shedding light on an important mechanism for the shoulders’
evolution. We select the same particles from the model snapshots discussed
above and integrate their orbits at four additional snapshots, with the
corresponding frozen potentials and bar pattern speeds, and using the
particles’ coordinates at each snapshot as initial conditions (for HG1, which
has star formation, we only use star-particles born before the selection
time). The results are shown in Figs. 21-24, with the selection times shown in
the central panels.

Fig. 21 shows the evolution of the frequency maps for the fiducial model SD1.
The main trends, which also occur in the other models with shoulders, are a
global shift upwards, i.e. to higher
$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$,
and to the left, i.e. to lower $\>\Omega_{z}/\Omega_{R}$. One obvious reason
for the shift upwards is the decline in $\>\Omega_{\mathrm{P}}$ as the bar
slows down (see Fig. 2), with its main resonances sweeping the disc outwards
and trapping new particles. Once trapped by the strong ILR, the orbits tend
not to shift any longer upwards in the frequency maps.

The ratio $\>\Omega_{z}/\Omega_{R}$ is a proxy for the vertical thinness of
the orbit (Sec. 4.1). Loosely speaking, the net shift to the left in the
frequency maps indicates evolution to thicker orbits, which is confirmed by
their large $z_{\rm max}$ (colors). While the thickness grows gradually for
decreasing $\>\Omega_{z}/\Omega_{R}$, the vertical motion is amplified when
orbits reach the vertical resonance $\mbox{$\>\Omega_{z}/\Omega_{R}$}=1$, and
a more dramatic amplification happens at the vILR. This is particularly
important near the crossing point of these two resonances and the ILR. While
large vertical excitation by the vILR for orbits close to the ILR
[$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$]
has been predicted by Binney (1981) and confirmed by e.g. Pfenniger & Friedli
(1991), Fig. 21 demonstrates this in a particularly transparent way. A
comparison with Fig. 12 also suggests that, once orbits cross the vILR, their
quasi-periodicity is destroyed, in agreement with Valluri et al. (2016) who
found that “banana” orbits are commonly chaotic – see also Manos et al.
(2022).

Note that Fig. 21 shows the evolution of a fixed set of particles selected in
the shoulders at $7\mbox{$\>{\rm Gyr}$}$, instead of the evolution of the
shoulders themselves. In particular, starting from $\approx 4\mbox{$\>{\rm
Gyr}$}$ this simulation persistently manifests shoulders which move outwards
following the bar growth – see Fig. 8 of Paper I.

Figure 22: Similar to Fig. 21, but for model SD1S. In this case, we do not
observe a significant time evolution of the frequency map, in agreement with
the very mild evolution of the bar parameters – see Fig. 2.

Fig. 22 shows the frequency map evolution for the model SD1S. In this case, a
small overall shift upwards and to the left is barely noticeable, in agreement
with the very slow bar evolution in this simulation (Fig. 2). We also note a
significant increase in $z_{\rm max}$ for orbits crossing to
$\mbox{$\>\Omega_{z}/\Omega_{R}$}<1$, but the vILR seems relevant only very
close to the crossing point with the ILR in this model, and there is no vILR
cloud.

Figure 23: Similar to Fig. 21, but for model HG1. We observe the same trends
as in model SD1, i.e. a net shift upwards and to the left. Orbits crossing the
vertical resonances $\mbox{$\>\Omega_{z}/\Omega_{R}$}=1$ and the vILR have
their vertical motion significantly amplified.

For HG1, the evolution of the frequency map is shown in Fig. 23. Once more, we
note a shift upwards and to the left. In this case, the resonance
$\mbox{$\>\Omega_{z}/\Omega_{R}$}=1$ also seems more effective in producing
large vertical excursions. In the vILR cloud, $z_{\rm max}$ is not as
significantly amplified as in the model SD1 (Fig. 21). However, the vILR cloud
is more populated in this model.

Figure 24: Similar to Fig. 21, but for Model 4. Besides the same trends
observed for SD1, we note the abrupt shift from $3\mbox{$\>{\rm Gyr}$}$ to
$4\mbox{$\>{\rm Gyr}$}$ towards smaller $\>\Omega_{z}/\Omega_{R}$. This is
produced by the buckling bar, slightly before $4\mbox{$\>{\rm Gyr}$}$, and
after which shoulders are detected in the model.

Finally, Fig. 24 shows the frequency map evolution for Model 4. The bar in
this model buckles slightly before $4\mbox{$\>{\rm Gyr}$}$, and shoulders are
observed right after that – see Sec. 2. The abrupt shift to lower
$\>\Omega_{z}/\Omega_{R}$ from $3\mbox{$\>{\rm Gyr}$}$ to $4\mbox{$\>{\rm
Gyr}$}$ is a consequence of this buckling, which vertically thickens the bar
region significantly. This interpretation and the frequency map at
$4\mbox{$\>{\rm Gyr}$}$ itself must be considered with caution, since this is
too close to the buckling and the gravitational potential may lack the
triaxial symmetry we assume in order to estimate it. On the other hand, the
main features in the frequency map at $4\mbox{$\>{\rm Gyr}$}$ are still
observed at $5\mbox{$\>{\rm Gyr}$}$, by which point the potential is more
symmetric, and therefore our fit is more reliable, with the same trend of a
shift upwards and to the left. This suggests that the frequencies at
$4\mbox{$\>{\rm Gyr}$}$ are reasonably accurate, despite the disturbance from
the buckling. In this model, the strong role of the vILR in exciting vertical
motion is again observed, producing a well populated vILR cloud with high
$z_{\rm max}$.

###  4.6 Evolution of shoulder-supporting orbits

In Sec. 4.5 we showed that orbits selected in the shoulder region were either
already trapped by the ILR at previous snapshots or previously had smaller
$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$,
with the latter evolving towards larger
$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$,
until possibly being trapped by the ILR. The orbits also evolve towards
smaller $\>\Omega_{z}/\Omega_{R}$, until possibly crossing the vILR, thereby
being excited to high $z_{\rm max}$ and becoming chaotic (e.g. Fig. 12). In
Sec. 4.4 we demonstrated the important role of vertically warm orbits at the
ILR, and of the vILR cloud, in building the shoulders. We now analyse the
morphological evolution of these warm orbits.

Figure 25: Evolution of the orbit-averaged density maps in the $x-y$ (top) and
$x-z$ planes (middle) for star-particles selected to be on warm orbits at the
ILR at $7\mbox{$\>{\rm Gyr}$}$ in model SD1. The coordinates are normalized by
the bar length at the selection time. Vertical shaded gray areas show the
shoulder region at the selection time. The bottom row shows the same orbits in
the frequency map. Those that at $7\mbox{$\>{\rm Gyr}$}$ have the looped-
shape, and support the shoulder morphology, evolve from less elongated orbits
and typically away from smaller $\>\Omega_{z}/\Omega_{R}$. Once they cross the
vILR, they spread in the vILR cloud, with much larger $z_{\rm max}$, diluting
the over-densities in the $x-y$ plane.

Fig. 25 shows, for the model SD1, the orbit-averaged surface density maps at
different times for the same star-particles selected to be warm orbits at the
ILR at $7\mbox{$\>{\rm Gyr}$}$, with the central column representing the
selection time. The upper (middle) row shows the maps in the $x-y$ ($x-z$)
plane, with coordinates normalized by the bar length at the selection time.
Here again, colors are normalized by the total number of star-particles in the
shoulder region. The shaded areas show the normalized shoulder region at the
selection time. The bottom panels show the same orbits in frequency maps.

The upper row shows that the warm orbits at the ILR at $7\mbox{$\>{\rm Gyr}$}$
mostly come from less elongated orbits which gradually become more elongated,
until being trapped by the ILR. A fraction of these orbits had
$\mbox{$\>\Omega_{z}/\Omega_{R}$}>3/2$ before the selection time, and either
before or after being trapped by the ILR. Before they cross the vILR
($t<7\mbox{$\>{\rm Gyr}$}$), these trends are accompanied by negligible
evolution in the $x-z$ plane. Once they cross the vILR ($7-8\mbox{$\>{\rm
Gyr}$}$), these orbits occupy the whole vILR cloud region, with large $z_{\rm
max}$. This is accompanied by a dissolution of the density excess in the $x-y$
plane. Additionally, the orbits start to shrink along the major axis, and
become confined to the inner parts at $9\mbox{$\>{\rm Gyr}$}$, no longer
falling in the original ($7\mbox{$\>{\rm Gyr}$}$) shoulder region. This is
closely connected to the abrupt change in the $x-z$ plane from $7\mbox{$\>{\rm
Gyr}$}$ to $8\mbox{$\>{\rm Gyr}$}$. As shown in Fig. 12, this is associated
with the chaotic nature of the orbits in the vILR cloud, which dilutes the
density distribution in the $x-y$ plane and weakens their contribution to the
shoulders. But perhaps most importantly, we see that, to the extent that the
vILR increases the vertical excursion, it decreases the in-plane amplitude,
converting in-plane motion to vertical motion, like a pole vaulter.

Figure 26: Similar to Fig. 25, but for model SD1S. As the bar barely evolves
in this model, the morphology of the warm orbits selected at 7 Gyr at ILR does
not change significantly.

For completeness, Fig. 26 shows the corresponding plot for model SD1S.
However, since the bar barely evolves in this model (see Fig. 2), we see no
significant evolution in this plot.

Fig. 27 shows a corresponding plot for the model HG1, where we observe the
same trends. However, in HG1 the vILR is less vigorous in exciting vertical
motion and, accordingly, the orbit-averaged density in the $x-y$ plane is less
disturbed for particles crossing the vILR to the vILR cloud. Fig. 28 shows the
equivalent plots for Model 4, and we see the same connection between the
dilution of the density excess in the $x-y$ plane and the increase in the
vertical motion after the orbits cross the vILR, populating the vILR cloud.

These results demonstrate the role of the vILR in perturbing the shoulder-
supporting orbits: once these orbits cross this resonance, their in-plane
motion is converted to vertical motion and they become chaotic, contributing
to diluting the shoulders.

Figure 27: Similar to Fig. 25, but for model HG1. For orbits crossing the
vILR, the amplification of $z_{\rm max}$ and the dilution of the over-
densities in the $x-y$ plane are smaller than in SD1, but we still observe the
vertical heating by the vILR. Figure 28: Similarly to Fig. 25, but for Model
4. We see the same trends observed in the models SD1 and HG1, but with the
abrupt change from $3\mbox{$\>{\rm Gyr}$}$ to $4\mbox{$\>{\rm Gyr}$}$ due to
the buckling bar.

##  5 Discussion

###  5.1 Comparison with different techniques

We have performed frequency analyses for orbits integrated in frozen
potentials of barred galaxy simulations. Frequency analysis of the coordinates
in ongoing $N$-body simulations (without introducing frozen potentials) has
also been proven possible (e.g. Ceverino & Klypin, 2007; Gajda et al., 2016;
Parul et al., 2020; Smirnov et al., 2021). While this has the apparent
advantage of tracking the orbital evolution in a more self-consistent manner,
the short time-series (typically $\approx 1\mbox{$\>{\rm Gyr}$}$), the
unavoidable changes in the potential in the considered time interval, and all
sources of noise in $N$-body potentials introduce uncertainties in the
estimated frequencies. On the other hand, orbit integration in frozen
potentials provides instantaneous orbital properties, and it allows the use of
densely sampled time-series (allowing investigation of higher frequencies),
with the narrowness of the main resonances in the frequency maps testifying to
the high accuracy of our procedure.

The frequency analysis involves the choice of the coordinate system, the
reference frame and the use of real or complex time-series. Regarding the last
point, Hunter (2002), using Cartesian coordinates for a few orbits, reports no
improvements in the accuracy of estimated frequencies when using complex time
series. Here, we chose cylindrical coordinates for their easier
interpretation, given the approximate symmetries of the potential, and complex
time-series. In particular, we use the complex Poincaré polar coordinates for
the azimuthal angle.

Difficulties in estimating the azimuthal frequency $\Omega_{\varphi}$ have
been reported in several works (e.g. Athanassoula, 2002, 2003; Ceverino &
Klypin, 2007), but often not in explicit terms. In Appendix B, exploring the
Isochrone Model, we show how to overcome these difficulties, demonstrating
that the use of the complex Poincaré polar coordinates improves the accuracy
of the azimuthal frequency and facilitates its extraction. The accuracy of
estimated frequencies in more generic potentials with (semi-)analytic
expressions for the fundamental frequencies deserves further investigation.

Unlike the common practice of performing frequency analysis in the rotating
frame of the bar, we opted to store the coordinates in the inertial frame of
the galaxy. As demonstrated in Appendix C, this choice allows a precise
identification of the co-rotation resonance, which is not possible, at least
with the techniques employed in this work, if the rotating frame of the bar is
adopted.

###  5.2 Implications for the evolution of shoulder profiles

Paper I demonstrated that shoulders are present in the density profiles of
growing bars. It also suggested the important contribution of looped orbits to
the shoulders, which is confirmed with the detailed orbital study performed
here (Sec. 4.4). We now return to the questions that motivated this work and
that we posed in Sec. 1.

The evolution of the density profile is complex and must depend to some extent
on all orbital groups, but we draw a simple picture that is based on the
evolution and prevalence of looped orbits along the ILR. In Sec. 4.4 we
demonstrated that the looped morphology in the $x-y$ plane is very well
defined for vertically warm orbits at the ILR (but not in the model SD1S,
where orbits are not looped), while cool orbits at the ILR are more
elliptical-like and contribute negligibly to the shoulders. Hot orbits at the
ILR (and at the vILR cloud) are chaotic and do not share the same looped
morphology of warm orbits, but overall still spend a considerable amount of
time in retrograde motion and produce a similar, more diluted, looped-like
shape in the $x-y$ plane. Their contribution to the shoulders depends on the
strength of the vILR in exciting vertical motion.

Thus, in order for shoulders to emerge in the density profiles, the ILR needs
to trap a significant number of orbits as the bar slows and its resonances
sweep through the disc. Additionally, a significant fraction of these orbits
need to develop loops at their ends, which typically happens in the warm
orbits at the ILR (${1\lesssim\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 3/2}$).
Since orbits in the disc and in the bar’s thinnest parts typically have
$\mbox{$\>\Omega_{z}/\Omega_{R}$}>3/2$, to support shoulders these orbits need
to decrease $\>\Omega_{z}/\Omega_{R}$, either before or after being trapped by
the ILR, which occurs when the orbits thicken vertically. While the bar
thickens, either via secular processes or via buckling, the shoulder-
supporting warm orbits at the ILR can sooner or later be populated. This is
compatible with the findings of Paper I that shoulders can appear in either
buckling or non-buckling bars and that, in buckling bars, shoulders commonly
appear just after a buckling event – see their Fig. 7.

Once shoulders have emerged in the density profile, their supporting orbits
can depart from the looped morphology in the $x-y$ plane if they move out of
the range $1<\mbox{$\>\Omega_{z}/\Omega_{R}$}<3/2$. This typically happens
with further bar thickening and orbits at the ILR crossing the vILR towards
$\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 1$ – see Figs. 25-28. Once more,
this thickening can happen either secularly or via a secondary buckling. Thus,
if the bar is not replenished with looped orbits, the shoulders can eventually
be erased. This explains why Paper I noted several instances of sudden erasure
of shoulders after a secondary buckling. It is worth noting that the newly
vertically excited orbits will contribute to a thicker potential, thus
decreasing the ratio $\>\Omega_{z}/\Omega_{R}$ of nearby orbits, bringing them
towards the vILR, in a positive feedback process.

Thus, the common thickening of orbits in the bar region would make shoulders
transient features, unless new orbits are trapped at the ILR with the
appropriate $\>\Omega_{z}/\Omega_{R}$ range and develop the looped morphology,
a process that seems to happen often for bars evolving secularly – see Fig. 8
of Paper I.

The fact that the looped morphology tends to appear in orbits at the ILR and
with $1\lesssim\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 3/2$ suggests some
connection with the vertical motion. This is intriguing since those orbits
were first found in 2D bars (Contopoulos, 1988). An interesting possibility is
that this looped morphology might be somehow triggered by vertical resonances
– for instance, comparing Figs. 12-15, we see that the model SD1S, the only
one without looped orbits or shoulders, is particularly unpopulated in the
region near the crossing of the ILR and
${\mbox{$\>\Omega_{z}/\Omega_{R}$}=3/2}$.

###  5.3 Connection with bulges

Investigating mechanisms for bar thickening, Sellwood & Gerhard (2020)
contrast three scenarios: a bar that buckles, one in which orbits are
vertically excited by staying trapped for significant times at the vILR
(Quillen, 2002), and another one in which orbits cross the vILR, get
vertically heated and remain that way after leaving the resonance (Quillen et
al., 2014). Having found evidence for the long-lasting trapping mechanism only
in an artificially vertically symmetrized model, Sellwood & Gerhard (2020)
conclude that the third mechanism is probably more relevant in real galaxies
without bar buckling.

Our results, in particular e.g. Figs. 21 and 25 (bottom row), confirm the
rapid crossing of the vILR, as opposed to a long trapping time, and the rapid
vertical excitation with a long-lasting X-morphology. It is worth mentioning
that the orbits at the vILR cloud, after crossing and being perturbed by the
vILR, can get very close to the galactic center (Fig. 8), possibly within the
sphere of influence of a supermassive black hole (SMBH). This can potentially
represent an important mechanism to persistently replenish the sphere of
influence of a SMBH with new stars. This possibility is even more interesting
if galactic bars can co-evolve in harmony with black holes, as recently
demonstrated by Wheeler et al. (2023).

In Sec. 4, we have shown that orbits crossing the vILR populate what we call
the vILR cloud, characterized by chaotic orbits with large $z_{\rm max}$ and
X-shapes in the $x-z$ plane. This morphology being produced by the vILR agrees
with several previous works (e.g. Pfenniger & Friedli, 1991), where the
contribution of the “banana” orbits to BP-bulges has been emphasized. On the
other hand, Portail et al. (2015) find that these orbits only contribute in
the bulge’s outermost regions, while orbits with $\Omega_{x}/\Omega_{z}\approx
3/5$ (“brezel” orbits) are more relevant for the bulge as a whole (in their
study, $\Omega_{x}$ is the frequency of oscillation about the bar major-axis
and can be seen as a proxy for our
$\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$}$) – see also
Abbott et al. (2017).

Secs. 4.3 and 4.6 demonstrate that after crossing the vILR the banana orbits
can also retain, to a certain extent, the loops near apocenters inherited from
warm orbits at the ILR. In fact, these loops have already been found in orbits
supporting either the thin (Patsis & Katsanikas, 2014) or the thick (Combes et
al., 1990) parts of bars.

Altogether, our results suggest that: i) the banana orbits can simultaneously
contribute to the outskirts of the bulge peanut shape (in the $x-z$ plane)
and, depending on the strength of the vertical perturbation, to the shoulders
(in the $x-y$ plane) – see e.g. Fig. 27 ($9\mbox{$\>{\rm Gyr}$}$); ii) in the
study of the orbital support of box-peanut (BP) bulges, the importance of the
vILR may be underestimated if we restrict to orbits trapped at it and exclude
the vILR cloud, since the latter is populated by orbits crossing the vILR and
acquiring a long-term BP-morphology.

In a detailed observational study with the $S^{4}G$ survey (Sheth et al.,
2010), Erwin et al. (2023) find some cases of galaxies with shoulders but no
BP-bulge, and no clear case of bars with BP-bulges and no shoulders,
suggesting that shoulders appear before BP-bulges. Given the role of the vILR
in producing BP-bulge-supporting orbits and the evolution of shoulder-
supporting orbits along the ILR, from the vertically warm ones to the vILR
cloud after crossing the vILR, the suggestion of Erwin et al. (2023) emerges
as a natural outcome.

The presence of shoulders has been found to correlate with massive (classical)
bulges, or centrally concentrated halos (Kim et al., 2015; Kruk et al., 2018;
Athanassoula & Misiriotis, 2002). Furthermore, Athanassoula (1992) found that
looped orbits tend to appear in galaxies with massive bulges (or high central
concentration, in general), and with more elongated and slowly rotating bars
(looped orbits in a slow bar have also been detected by Chaves-Velasquez et
al., 2017). This is in agreement with the important contribution of looped
orbits to shoulders, suggested by Paper I and confirmed here. Additionally,
among our three models with shoulders, two (SD1 and Model 4) have slow bars
and one (HG1, with weak shoulders) has a fast bar, suggesting that shoulders
are not restricted to slow bars, but can also be present, with varying
strength, in fast bars. All models have clear signs of BP-bulges, although
they are weaker for models SD1S and HG1, which suggests that the correlation
shoulder-bulge is not restricted to either BP- or classical bulges.

A possible interpretation of the shoulders-massive bulges correlation might be
that massive bulges induce thicker orbits in their outskirts, facilitating
occupation of the shoulder-supporting warm orbits at the ILR, while bars in
bulgeless galaxies would tend to trap stars in thinner orbits, i.e. with
larger $\>\Omega_{z}/\Omega_{R}$ such as the cool orbits at the ILR, with a
morphology that does not support shoulders – see Fig. 7. This line of thought
would shift a massive bulge towards being the cause for shoulders, which would
be in contradiction with the apparent precedence of shoulders with respect to
BP-bulges suggested by Erwin et al. (2023).

On the other hand, warm orbits at the ILR crossing the vILR are vertically
excited and change into chaotic orbits. As we have shown, this is accompanied
by these orbits shrinking in the $x-y$ plane and passing close to the center
(see e.g. Fig. 17, panels $b$ and $c$), overall increasing the density of the
central regions. This may explain the sudden increase in the central surface
density after buckling events observed in simulated galaxies (Raha et al.,
1991; Debattista et al., 2006), since this brings orbits towards the vILR –
see Fig. 24. Therefore, this mass transfer to the center may contribute to the
growth of massive bulges which, instead of causing the appearance of
shoulders, would be the consequence of the evolution of shoulder supporting
orbits towards the vILR.

Kim et al. (2015) suggest that shoulders are indicative of dynamically evolved
bars, which should be formed with exponential profiles inherited from the
disc. They also suggest that this might explain the correlation with the
galaxy’s stellar mass and Hubble type, with more massive early-type spirals
forming bars earlier. Although the eventual erasure of shoulders found by
Paper I complicates this narrative, this would agree with the strong
dependence of shoulders on stellar mass found by Erwin et al. (2023). Our
finding that the most shoulder-supporting morphology appears in orbits at the
ILR with intermediate thickness (thus requiring a minimum thickness) is also
in line with shoulders being associated with dynamically evolved bars.
Additionally, the evolution of these orbits towards vertically thicker orbits
supporting bulges also agrees with the suggestion of Erwin et al. (2023) that
shoulders seem to precede BP-bulges.

###  5.4 Shoulder definition and possible observational implications

Figure 29: Similar to Fig. 5, but only showing $\>E_{\rm J}$ (points) for warm
orbits at the ILR. These orbits are confined to within the UHr.

Our results show that the most relevant orbits to shoulders are the looped
ones at the vertically warm ILR. In this way, one can expect shoulders to be
restricted to the region where these orbits are present. Similarly to Fig. 5,
in Fig. 29 we show the effective potential in the bar-major axis, as well as
the approximate locations of the co-rotation, UHR and the ILR. We also show
the Jacobi integral, but now restricting to warm orbits at the ILR.

It is clear that these orbits are restricted to the region inside the Ultra
Harmonic radius. In Fig. 5, we see that for the model SD1 the outer edge of
the shoulder approximately coincides with the Ultra Harmonic radius. However,
for the other simulations, the shoulder outer edge (as detected in Paper I)
can be significantly outside this radius.

This observation is relevant for two reasons: on one hand, this suggests that
identifying the outer edge of the shoulders based on a minimum of the
curvature radius of the density profile (as implemented in Paper I) might be
overestimating the shoulder extent, if this is to be identified with the
region of looped orbits at the ILR. If that is the case, our current selection
may be contaminated by a significant number of orbits that are not associated
with the shoulders. Thus, the contribution of the looped orbits to shoulders
we inferred in Sec. 4.4 represents a lower bound to their true importance to
shoulders. On the other hand, if we can infer that the shoulders should end at
the Ultra Harmonic radius, this can potentially be translated into constraints
on the mass distribution and/or bar pattern speed of galaxies where shoulders
are detected, or even in the Milky Way.

##  6 Summary

In this paper, we explored three pure $N$-body models and one $N$-body+SPH
model which form bars self-consistently. Three of these models, including the
one with gas, develop shoulders in the bar-major axis density profile. Among
these, one has a bar that buckles. We select star-particles in snapshots and
regions where shoulders are identified, and integrate their orbits at
different times in frozen potentials rotating with the corresponding bar
pattern speeds. We perform frequency analysis of these orbits, using
cylindrical coordinates in the inertial frame and selecting the leading
frequency, presenting results in the frequency space
$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$ vs
$\>\Omega_{z}/\Omega_{R}$. We evaluate the contribution of different orbital
types to the density in the shoulders region, and follow their evolution. Our
main conclusions are:

  * •

We confirm the suggestion of Paper I of the importance of orbits with loops at
their ends (in the $x-y$ plane) to the shoulders in bar density profiles.

  * •

We show that this orbital shape typically happens for orbits at the ILR
[$\mbox{$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$}=1/2$]
and with $1\lesssim\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 3/2$ (vertically
warm). Among the orbits at the ILR, those with
$\mbox{$\>\Omega_{z}/\Omega_{R}$}\gtrsim 3/2$ (cool orbits) have a mildly
elliptical shape and contribute negligibly to the shoulders.

  * •

As the bar slows down, the main resonances sweep the disc and, for a fixed set
of particles, the ratio
$\>(\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$})/\Omega_{R}$
increases, until they are trapped by the ILR. While the bar thickens, the
ratio $\>\Omega_{z}/\Omega_{R}$ tends to decrease, causing a significant
number of orbits to cross the vILR (Figs. 21-24).

  * •

In agreement with Quillen et al. (2014) and Sellwood & Gerhard (2020), we find
that the vILR does not trap stars for long times. However, orbits crossing it
are permanently vertically excited, developing an $X$ shape in the $x-z$ plane
– see also Wheeler et al. (2023).

  * •

Warm orbits at the ILR crossing the vILR towards
$\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 1$ are vertically excited at the
expense of in-plane motion, become chaotic and populate a region around the
vILR in the frequency map that we dub the vILR cloud. In comparison to warm
orbits at the ILR, those at the vILR cloud are less extended along the
$x-$axis and produce a more diffuse net figure in the $x-y$ plane (Figs.
25-28).

  * •

Orbits at the vILR cloud also contribute significantly to shoulders, but the
larger their vertical excursions, the less shoulder-supporting they are.

  * •

As the bar thickens, the orbital support for shoulders can increase with
orbits at the ILR entering the range of warm orbits, and potentially decrease
once they move to $\mbox{$\>\Omega_{z}/\Omega_{R}$}\lesssim 1$ (where they
cross the vILR). This explains why shoulders can develop right after a
buckling event but can also be erased right after a second buckling (as noted
in Paper I), since these can suddenly either populate or de-populate those
orbits. This simple evolutionary sequence supports the suggestion of Erwin et
al. (2023) that shoulders seem to appear before BP-bulges (with a possible co-
formation of shoulders and BP-bulges in buckling bars).

  * •

An analysis with the isochrone model (Appendix B) shows that for orbits with
large apocenters, the leading frequency in the azimuthal spectrum is the
precession one, instead of $\>\Omega_{\varphi}$. This has possible
implications for the analysis of Milky Way halo stars and deserves further
investigation. Moreover, frequency estimates with the time-series
$f_{\varphi}=\sqrt{2|L_{z}|}(\cos\varphi+i\sin\varphi)$ provide higher
precision (compared to $f_{\varphi}=\varphi$) and facilitates its extraction –
see Fig. 30.

In conclusion, bar thickening and vertical resonances can dilute the shoulders
in the bar major-axis density profiles, unless the bar keeps trapping new
orbits in the appropriate frequency ratios, with the looped morphology. This
sheds light on the conclusion of Paper I that long-term shoulders require a
growing bar. Our results help to build a coherent picture tying together known
analytical results with several recent observational and simulation-based
results on the properties and evolution of shoulders in the density profiles
of bars.

## Acknowledgements

We thank Eugene Vasiliev for his help with agama and for always enriching
discussions. We also thank Rubens Machado and João Amarante for a careful
reading and comments. We are also truly grateful to the anonymous referee,
whose comments contributed to improve this paper. V.P.D. and L.B.S. were
supported by STFC Consolidated grant ST/R000786/1. L.B.S and M.V. acknowledge
the support of NASA-ATP award 80NSSC20K0509 and U.S. National Science
Foundation AAG grant AST-2009122. KJD is grateful for support from the Simons
Foundation for supporting her sabbatical stay at the CCA, Flatiron Institute.
The orbital analysis was carried out on Stardynamics, which was funded from
Newton Advanced Fellowship # NA150272 awarded by the Royal Society and the
Newton Fund. The simulations were run at the High Performance Computing
Facility of the University of Central Lancashire and at the DiRAC Shared
Memory Processing system at the University of Cambridge, operated by the
COSMOS Project at the Department of Applied Mathematics and Theoretical
Physics on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This
equipment was funded by BIS National E-infrastructure capital grant
ST/J005673/1, STFC capital grant ST/H008586/1 and STFC DiRAC Operations grant
ST/K00333X/1. DiRAC is part of the National E-Infrastructure

##  Appendix A Frequency analysis algorithm

Here we present the algorithm used to extract frequencies of orbits, which is
essentially the same as NAFF, proposed by Laskar et al. (1992) and further
developed by Valluri & Merritt (1998) – see also Price-Whelan et al. (2016).
We refer the reader to these papers for a complete description. We also
comment on some details and improvements implemented in the python package
used in this work, which we dub naif (Beraldo e Silva & Valluri, 2023) – the
documentation is available at http://naif.readthedocs.io/en/latest/index.html.
This preliminary version extracts peak frequencies very accurately, as
demonstrated in Appendix B, and will be improved in the future with several
tools for frequency analysis.

Assume an orbit is integrated for a total time $T$ and let, for simplicity of
notation, the time variable be defined in the symmetric interval $-T/2\leq
t\leq T/2$. In practice, $t$ is an array with $N$ elements $t_{n}$, and we
have a (real or complex) time-series $f(t_{n})$ associated with each
coordinate. We start by performing a windowed Discrete Fourier Transform

$F_{j}=\frac{1}{N}\sum_{n=0}^{N-1}f(t_{n})\chi_{p}(t_{n})e^{-2\pi inj/N},$
(A1)

where $j=-N/2,\dots,0\dots N/2$, and

$\chi_{p}(t)=\frac{2^{p}(p!)^{2}}{(2p)!}\left[1+\cos\left(\frac{2\pi
t}{T}\right)\right]^{p}$ (A2)

is the window function, with $p\in\mathbb{N}$. In this work we use $p=1$, but
other values might improve the estimate (Laskar, 1999; Hunter, 2002).

The spectrum is calculated at frequencies $\omega_{j}=2\pi j/T$, and the
minimum frequency that can be probed is $\omega_{\mathrm{min}}=2\pi/T$. On the
other hand, the maximum (Nyquist) frequency is
$\Omega_{\mathrm{Ny}}=\Omega_{j=N/2}=N\pi/T$. If we conservatively require
$\Omega_{\mathrm{Ny}}\approx 10^{3}\mbox{$\>{\rm Gyr}$}^{-1}$, this requires
$N\gtrsim(10^{3}/\pi)(T/\mbox{$\>{\rm Gyr}$})$. For an orbit with period of
circular motion $T_{\mathrm{circ}}\approx 0.25$ Gyr, integrated for $T\approx
20T_{\mathrm{circ}}\approx 5$ Gyr ($T\approx 40T_{\mathrm{circ}}\approx 10$
Gyr), we have $\omega_{\mathrm{min}}\approx 1.25\mbox{$\>{\rm Gyr}$}^{-1}$
($\omega_{\mathrm{min}}\approx 0.63\mbox{$\>{\rm Gyr}$}^{-1}$), and we need
$N\gtrsim 1500$ ($N\gtrsim 3000$) points. Using $N\sim 10^{4}$ points is thus
safe for most applications.

A rough estimate of the leading frequency is the spectral line with maximum
amplitude $|F_{j}|$, $\>\omega_{\mathrm{max}}$. Since the spectrum is defined
in the discrete set $\omega_{j}=2\pi j/T$, this estimate’s precision is
limited by the resolution $\Delta\omega=2\pi/T$. This estimate is then refined
by calculating the projection $\phi(\omega)=\langle f(t),e^{i\omega t}\rangle$
of the time-series $f(t)$ onto the “frequency vector” $e^{i\omega t}$, where

$\langle
f(t),g(t)\rangle\equiv\frac{1}{T}\int_{-T/2}^{T/2}f(t)g^{*}(t)\chi_{p}(t)\,\mathrm{d}t.$
(A3)

The refined frequency is found by maximizing $|\phi(\omega)|$ around
$\>\omega_{\mathrm{max}}$. The maximum is determined with the optimized
Brent’s algorithm, as implemented in scipy. If the maximum cannot be
identified in this first shot, a local maximum is identified via a brute force
search. The interval around $\>\omega_{\mathrm{max}}$ used in the original
NAFF implementation is $\mbox{$\>\omega_{\mathrm{max}}$}\pm 2\pi/T$, but we
have found that the interval $\mbox{$\>\omega_{\mathrm{max}}$}\pm 4\pi/T$ is
more appropriate to guarantee the Brent’s algorithm will find the maximum, and
this is the interval used in naif. Having thus identified this first leading
frequency $\omega_{1}$, we calculate the complex amplitude $a_{1}$ associated
with it, where

$a_{k}=\langle f(t),u_{k}\rangle,$ (A4)

and $u_{k}=e^{i\omega_{k}t}$. If we want to extract additional frequencies, we
“subtract” this spectral line from the time-series, defining in general

$f_{k}(t)=f_{k-1}(t)-a_{k}e_{k},$ (A5)

where $e_{1}=u_{1}$, $f_{k}(t)$ is the time-series after extracting the $k$-th
frequency $\omega_{k}$, and $f_{0}(t)$ is the original one. Here, $e_{k}$
denotes an orthonormal basis, built from the $u_{k}$’s via a Gram–Schmidt
process:

$e_{k}=u_{k}-\sum_{j=1}^{k-1}\langle u_{k},e_{j}\rangle\,e_{j}.$ (A6)

We then iterate through Eqs. A1-A6 to extract as many frequencies as desired.

Compared to the original NAFF algorithm, a small (but relevant) improvement of
naif is the following: NAFF first applies a DFT without windowing to look for
the rough estimate $\>\omega_{\mathrm{max}}$, and then uses the window
function $\chi_{p}(t)$ in the refinement step when maximizing
$|\phi(\omega)|$. This different procedure to look for the same frequency ends
up changing the order of the highest amplitudes in some cases. This is later
corrected by sorting the frequencies by amplitude and taking the leading one,
but requires the extraction of some frequencies even if we are only interested
in the leading one. On the other hand, naif performs a windowed DFT, Eq. A1,
already in the first search for $\>\omega_{\mathrm{max}}$, which ends up being
always the leading frequency. This removes the need to extract additional
frequencies if we are only interested in the leading one.

Finally, naif provides easy access to the spectrum before the extraction of
each frequency and also to $|\phi(\omega)|$ around each extracted frequency,
which can be useful for detailed investigations of individual orbits.

##  Appendix B Numerical frequencies for the isochrone model

We use agama to generate a self-consistent sample of the isochrone model and
integrate each orbit for $\approx 100$ azimuthal periods. A self-consistent
ensemble is not strictly necessary, since we simply compare frequencies of
orbits individually and do not focus on their relative numbers. However, this
guarantees a broad phase-space region for investigation. Since this is a
spherical model, we only estimate the radial and azimuthal frequencies,
$\Omega_{r}$ and $\>\Omega_{\varphi}$, respectively (the transversal frequency
$\Omega_{\vartheta}=\pm\mbox{$\>\Omega_{\varphi}$}$). For this model, we set
$GM=1$.

Fig. 30 shows the relative errors of the estimated frequencies for $10^{3}$
orbits, in comparison to the analytic expressions for the isochrone model
(Binney & Tremaine, 2008). Points are color-coded by the apocenter radii (in
units of the scale radius). The left column shows results obtained using as
input the real time-series

$\begin{array}[]{l}f_{r}=r\\\ f_{\varphi}=\varphi,\end{array}$ (B1)

and the right column shows results for the complex time-series (Papaphilippou
& Laskar, 1996, 1998)

$\begin{array}[]{l}f_{r}=r+iv_{r}\\\
f_{\varphi}=\sqrt{2|L_{z}|}(\cos\varphi+i\sin\varphi),\end{array}$ (B2)

where $v_{r}=dr/dt$ is the radial velocity and $L_{z}$ is the $z$-component of
the angular momentum. For the radial coordinate (first row), the absolute
value of the leading frequency has a relative error
$|\delta\Omega_{r}/\Omega_{r}|\approx 10^{-6}$, using either the real or the
complex time-series.

Figure 30: Relative errors in the fundamental frequencies for $10^{3}$ orbits
in the isochrone potential. Left column: real time-series, Eq. B1. Right
column: complex time-series, Eq. B2. The radial frequency (top row) is
estimated as the absolute value of the leading frequency. Real and complex
time-series give similar errors. For $\>\Omega_{\varphi}$ (middle row), the
real time series requires setting
$\operatorname{sign}(\mbox{$\>\Omega_{\varphi}$})=\operatorname{sign}(L_{z})$,
while this is automatically set for the complex time-series. Requiring
$\mbox{$\>\Omega_{\varphi}$}>\Omega_{r}/2$ (see main text) and taking the
leading frequency, we have
$|\delta\mbox{$\>\Omega_{\varphi}$}/\mbox{$\>\Omega_{\varphi}$}|\lesssim
10^{-3}$ (bottom row). The complex time-series provides $\approx 10\times$
smaller errors for $\>\Omega_{\varphi}$.

For the azimuthal angle (second row), using the real time-series
${f_{\varphi}=\varphi}$ the leading frequencies provide errors
$|\delta\mbox{$\>\Omega_{\varphi}$}/\mbox{$\>\Omega_{\varphi}$}|\lesssim
10^{-3}$ for a fraction of the orbits, but only after imposing
${\operatorname{sign}(\mbox{$\>\Omega_{\varphi}$})=\operatorname{sign}(L_{z})}$,
while the complex time-series (Eq. B2) does not require this intervention.
Additionally, there are a few orbits with small analytic $\>\Omega_{\varphi}$
and $|\delta\mbox{$\>\Omega_{\varphi}$}/\mbox{$\>\Omega_{\varphi}$}|\approx 1$
for both real and complex time-series. These orbits have large apocenters
(${r_{\mathrm{apo}}/a\gtrsim 15}$), and expend most of the time far from the
center, effectively behaving as precessing ellipses in a Kepler-like
potential. As a consequence, the leading frequency in the spectrum is
$\Omega_{r}-\mbox{$\>\Omega_{\varphi}$}$ (the precessing frequency), although
$\>\Omega_{\varphi}$ is also present, with a smaller amplitude. In a Kepler-
like potential, $\Omega_{r}\approx\Omega_{\varphi}$, thus the leading
frequency is $\approx 0$, which implies
$|\delta\mbox{$\>\Omega_{\varphi}$}/\mbox{$\>\Omega_{\varphi}$}|\approx 1$. In
realistic galactic potentials, typical stars have
$\Omega\lesssim\kappa\lesssim 2\Omega$. Thus, if we heuristically require
$|\mbox{$\>\Omega_{\varphi}$}|>|\Omega_{r}|/2$, and take the leading frequency
respecting this restriction, the fundamental $\>\Omega_{\varphi}$ is recovered
with good accuracy
($|\delta\mbox{$\>\Omega_{\varphi}$}/\mbox{$\>\Omega_{\varphi}$}|\lesssim
10^{-3}$) for almost all orbits (third row). Finally, the complex time-series
provides a $\approx 10\times$ better precision (although the estimates are
almost identical to those of the real series if orbits are integrated for only
$\approx 50$ circular periods).

##  Appendix C Inertial vs rotating frame

Here we compare the frequencies obtained using coordinates in the inertial
frame and in the bar rotating frame. We use the model SD1S, which has a fast
bar at $7\mbox{$\>{\rm Gyr}$}$ and the region where star-particles are
selected is near the co-rotation radius – see Fig. 4. Thus, we can expect the
co-rotation resonance to be populated in the frequency map, as is the case
when we use the inertial frame (Fig. 13). In this case, the frequency analysis
outputs $\>\Omega_{\varphi}$ and we literally subtract the bar pattern speed
$\>\Omega_{\mathrm{P}}$ before plotting the frequency map in Fig. 13.

Fig. 31 shows the frequency map for the same snapshot and the same orbits, but
now using coordinates in the bar rotating frame. With $\>\Omega_{\varphi}$
still representing the azimuthal frequency in the inertial frame, the output
of the frequency analysis is now
$\mbox{$\>\Omega_{\varphi}$}-\mbox{$\>\Omega_{\mathrm{P}}$}$.

Figure 31: Frequency maps for the same orbits as in Fig. 13, but here with
frequencies estimated using coordinates in the rotating frame. The two
frequency maps are almost identical, except for the co-rotation resonance,
which is not correcly identified in this case.

Comparing Figs. 13 and 31, we see that the two frequency maps are almost
identical, except for the co-rotation
($\mbox{$\>\Omega_{\varphi}$}=\mbox{$\>\Omega_{\mathrm{P}}$}$), which is
sharply defined in the inertial frame, but is completely lost if we use
coordinates in the rotating frame. This is probably due to the leading
frequency in this case being the libration around the co-rotation, which might
be relevant for detailed investigations of orbits trapped at co-rotation, but
it is not of interest for this work.

## References

  * Abbott et al. (2017) Abbott, C. G., Valluri, M., Shen, J., & Debattista, V. P. 2017, MNRAS, 470, 1526, doi: 10.1093/mnras/stx1262
  * Anderson et al. (2022) Anderson, S. R., Debattista, V. P., Erwin, P., et al. 2022, MNRAS, 513, 1642, doi: 10.1093/mnras/stac913
  * Arnol’d (1963) Arnol’d, V. I. 1963, Russ. Math. Surv., 18, 9, doi: 10.1070/rm1963v018n05ABEH004130
  * Athanassoula (1992) Athanassoula, E. 1992, MNRAS, 259, 328, doi: 10.1093/mnras/259.2.328
  * Athanassoula (2002) —. 2002, ApJL, 569, L83, doi: 10.1086/340784
  * Athanassoula (2003) —. 2003, MNRAS, 341, 1179, doi: 10.1046/j.1365-8711.2003.06473.x
  * Athanassoula (2013) —. 2013, in Secular Evolution of Galaxies, ed. J. Falcón-Barroso & J. H. Knapen, 305
  * Athanassoula & Beaton (2006) Athanassoula, E., & Beaton, R. L. 2006, MNRAS, 370, 1499, doi: 10.1111/j.1365-2966.2006.10567.x
  * Athanassoula et al. (1983) Athanassoula, E., Bienayme, O., Martinet, L., & Pfenniger, D. 1983, A&A, 127, 349
  * Athanassoula & Misiriotis (2002) Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35, doi: 10.1046/j.1365-8711.2002.05028.x
  * Beraldo e Silva et al. (2019) Beraldo e Silva, L., de Siqueira Pedra, W., Valluri, M., Sodré, L., & Bru, J.-B. 2019, ApJ, 870, 128, doi: 10.3847/1538-4357/aaf397
  * Beraldo e Silva & Valluri (2023) Beraldo e Silva, L., & Valluri, M. 2023, naif: Frequency analysis package, Astrophysics Source Code Library, record ascl:2303.004.  http://ascl.net/2303.004
  * Binney (1981) Binney, J. 1981, MNRAS, 196, 455, doi: 10.1093/mnras/196.3.455
  * Binney & Spergel (1982) Binney, J., & Spergel, D. 1982, ApJ, 252, 308, doi: 10.1086/159559
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bullock et al. (2001) Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240, doi: 10.1086/321477
  * Bureau & Athanassoula (2005) Bureau, M., & Athanassoula, E. 2005, ApJ, 626, 159, doi: 10.1086/430056
  * Carpintero & Aguilar (1998) Carpintero, D. D., & Aguilar, L. A. 1998, MNRAS, 298, 1, doi: 10.1046/j.1365-8711.1998.01320.x
  * Ceverino & Klypin (2007) Ceverino, D., & Klypin, A. 2007, MNRAS, 379, 1155, doi: 10.1111/j.1365-2966.2007.12001.x
  * Chaves-Velasquez et al. (2017) Chaves-Velasquez, L., Patsis, P. A., Puerari, I., Skokos, C., & Manos, T. 2017, ApJ, 850, 145, doi: 10.3847/1538-4357/aa961a
  * Chirikov (1979) Chirikov, B. V. 1979, PhR, 52, 263, doi: 10.1016/0370-1573(79)90023-1
  * Cole et al. (2014) Cole, D. R., Debattista, V. P., Erwin, P., Earp, S. W. F., & Roškar, R. 2014, MNRAS, 445, 3352, doi: 10.1093/mnras/stu1985
  * Collier et al. (2018) Collier, A., S., I., & Heller, C. 2018, MNRAS, 476, 1331, doi: 10.1093/mnras/sty270
  * Combes et al. (1990) Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A&A, 233, 82
  * Combes & Elmegreen (1993) Combes, F., & Elmegreen, B. G. 1993, A&A, 271, 391
  * Contopoulos (1978) Contopoulos, G. 1978, A&A, 64, 323
  * Contopoulos (1981) —. 1981, A&A, 102, 265
  * Contopoulos (1988) —. 1988, A&A, 201, 44
  * Contopoulos & Grosbol (1989) Contopoulos, G., & Grosbol, P. 1989, A&ARv, 1, 261, doi: 10.1007/BF00873080
  * Contopoulos & Papayannopoulos (1980) Contopoulos, G., & Papayannopoulos, T. 1980, A&A, 92, 33
  * de Carvalho & da Costa (1987) de Carvalho, R. R., & da Costa, L. N. 1987, A&A, 171, 66
  * Debattista et al. (2020) Debattista, V. P., Liddicott, D. J., Khachaturyants, T., & Beraldo e Silva, L. 2020, MNRAS, 498, 3334
  * Debattista et al. (2006) Debattista, V. P., Mayer, L., Carollo, C. M., et al. 2006, ApJ, 645, 209, doi: 10.1086/504147
  * Debattista et al. (2017) Debattista, V. P., Ness, M., Gonzalez, O. A., et al. 2017, MNRAS, 469, 1587, doi: 10.1093/mnras/stx947
  * Debattista & Sellwood (1998) Debattista, V. P., & Sellwood, J. A. 1998, ApJL, 493, L5, doi: 10.1086/311118
  * Debattista & Sellwood (2000) —. 2000, ApJ, 543, 704, doi: 10.1086/317148
  * Dodd et al. (2022) Dodd, E., Helmi, A., & Koppelman, H. H. 2022, A&A, 659, A61, doi: 10.1051/0004-6361/202141354
  * Elmegreen & Elmegreen (1985) Elmegreen, B. G., & Elmegreen, D. M. 1985, ApJ, 288, 438, doi: 10.1086/162810
  * Elmegreen et al. (1996) Elmegreen, B. G., Elmegreen, D. M., Chromey, F. R., Hasselbacher, D. A., & Bissell, B. A. 1996, AJ, 111, 2233, doi: 10.1086/117957
  * Elmegreen et al. (2011) Elmegreen, D. M., Elmegreen, B. G., Yau, A., et al. 2011, ApJ, 737, 32, doi: 10.1088/0004-637X/737/1/32
  * Erwin (2018) Erwin, P. 2018, MNRAS, 474, 5372, doi: 10.1093/mnras/stx3117
  * Erwin et al. (2023) Erwin, P., Debattista, V. P., & Anderson, S. R. 2023, MNRAS, stad1972, doi: 10.1093/mnras/stad1972
  * Eskridge et al. (2000) Eskridge, P. B., Frogel, J. A., Pogge, R. W., et al. 2000, AJ, 119, 536, doi: 10.1086/301203
  * Fragkoudi et al. (2021) Fragkoudi, F., Grand, R. J. J., Pakmor, R., et al. 2021, A&A, 650, L16, doi: 10.1051/0004-6361/202140320
  * Gadotti et al. (2007) Gadotti, D. A., Athanassoula, E., Carrasco, L., et al. 2007, MNRAS, 381, 943, doi: 10.1111/j.1365-2966.2007.12295.x
  * Gajda et al. (2016) Gajda, G., Łokas, E. L., & Athanassoula, E. 2016, ApJ, 830, 108, doi: 10.3847/0004-637X/830/2/108
  * Gardner et al. (2014) Gardner, E., Debattista, V. P., Robin, A. C., Vásquez, S., & Zoccali, M. 2014, MNRAS, 438, 3275, doi: 10.1093/mnras/stt2430
  * Hamilton et al. (2022) Hamilton, C., Tolman, E. A., Arzamasskiy, L., & Duarte, V. N. 2022, arXiv e-prints, arXiv:2208.03855, doi: 10.48550/arXiv.2208.03855
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hernquist & Weinberg (1992) Hernquist, L., & Weinberg, M. D. 1992, ApJ, 400, 80, doi: 10.1086/171975
  * Hunter (2002) Hunter, C. 2002, SSRv, 102, 83, doi: 10.1023/A:1021360731798
  * Kim et al. (2015) Kim, T., Sheth, K., Gadotti, D. A., et al. 2015, ApJ, 799, 99, doi: 10.1088/0004-637X/799/1/99
  * Kolmogorov (1954) Kolmogorov, A. 1954, Dokl. Akad. Nauk. SSR, 98, 2
  * Kruk et al. (2018) Kruk, S. J., Lintott, C. J., Bamford, S. P., et al. 2018, MNRAS, 473, 4731, doi: 10.1093/mnras/stx2605
  * Kuijken & Dubinski (1995) Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341
  * Laskar (1990) Laskar, J. 1990, Icar, 88, 266, doi: 10.1016/0019-1035(90)90084-M
  * Laskar (1999) Laskar, J. 1999, Introduction to Frequency Map Analysis, ed. C. Simó (Dordrecht: Springer Netherlands), 134–150, doi: 10.1007/978-94-011-4673-9_13
  * Laskar et al. (1992) Laskar, J., Froeschlé, C., & Celletti, A. 1992, PhyD, 56, 253, doi: 10.1016/0167-2789(92)90028-L
  * Long et al. (2014) Long, S., Shlosman, I., & Heller, C. 2014, ApJ, 783, 18
  * Lucey et al. (2022) Lucey, M., Pearson, S., Hunt, J. A. S., et al. 2022, arXiv, arXiv:2206.01798.  https://arxiv.org/abs/2206.01798
  * Manos et al. (2022) Manos, T., Skokos, C., & Patsis, P. A. 2022, MNRAS, 509, 1995, doi: 10.1093/mnras/stab2987
  * Martinez-Valpuesta et al. (2006) Martinez-Valpuesta, I., Shlosman, I., & Heller, C. 2006, ApJ, 637, 214, doi: 10.1086/498338
  * 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
  * Moser (1962) Moser, J. 1962, Nachr. Akad. Wiss. Gött., II. Math.-Phys. Kl., 1962, 1
  * Nair & Abraham (2010) Nair, P. B., & Abraham, R. G. 2010, ApJ, 714, L260, doi: 10.1088/2041-8205/714/2/L260
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * O’Neill & Dubinski (2003) O’Neill, J. K., & Dubinski, J. 2003, MNRAS, 346, 251, doi: 10.1046/j.1365-2966.2003.07085.x
  * Papaphilippou & Laskar (1996) Papaphilippou, Y., & Laskar, J. 1996, A&A, 307, 427
  * Papaphilippou & Laskar (1998) —. 1998, A&A, 329, 451
  * Papayannopoulos & Petrou (1983) Papayannopoulos, T., & Petrou, M. 1983, A&A, 119, 21
  * Parul et al. (2020) Parul, H. D., Smirnov, A. A., & Sotnikova, N. Y. 2020, ApJ, 895, 12, doi: 10.3847/1538-4357/ab76ce
  * Patsis & Athanassoula (2019) Patsis, P. A., & Athanassoula, E. 2019, MNRAS, 490, 2740, doi: 10.1093/mnras/stz2588
  * Patsis & Katsanikas (2014) Patsis, P. A., & Katsanikas, M. 2014, MNRAS, 445, 3546, doi: 10.1093/mnras/stu1970
  * Patsis et al. (2002) Patsis, P. A., Skokos, C., & Athanassoula, E. 2002, MNRAS, 337, 578, doi: 10.1046/j.1365-8711.2002.05943.x
  * Petersen et al. (2021) Petersen, M. S., Weinberg, M. D., & Katz, N. 2021, MNRAS, 500, 838, doi: 10.1093/mnras/staa3202
  * Pfenniger & Friedli (1991) Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75
  * Pontzen et al. (2013) Pontzen, A., Roškar, R., Stinson, G., & Woods, R. 2013, pynbody: N-Body/SPH analysis for python, Astrophysics Source Code Library, record ascl:1305.002.  http://ascl.net/1305.002
  * Portail et al. (2015) Portail, M., Wegg, C., & Gerhard, O. 2015, MNRAS, 450, L66, doi: 10.1093/mnrasl/slv048
  * Price-Whelan et al. (2016) Price-Whelan, A. M., Johnston, K. V., Valluri, M., et al. 2016, MNRAS, 455, 1079, doi: 10.1093/mnras/stv2383
  * Prieto et al. (2001) Prieto, M., Aguerri, J. A. L., Varela, A. M., & Muñoz-Tuñón, C. 2001, A&A, 367, 405, doi: 10.1051/0004-6361:20000426
  * Quillen (2002) Quillen, A. C. 2002, AJ, 124, 722, doi: 10.1086/341753
  * Quillen et al. (2014) Quillen, A. C., Minchev, I., Sharma, S., Qin, Y.-J., & Di Matteo, P. 2014, MNRAS, 437, 1284, doi: 10.1093/mnras/stt1972
  * Raha et al. (1991) Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Nature, 352, 411, doi: 10.1038/352411a0
  * Regan & Elmegreen (1997) Regan, M. W., & Elmegreen, D. M. 1997, AJ, 114, 965, doi: 10.1086/118527
  * Roy (1991) Roy, A. 1991, NATO Advanced Study Institute (ASI) Series B, Vol. 272, Predictability, stability, and chaos in N-body dynamical systems. Proceedings.
  * Sellwood & Gerhard (2020) Sellwood, J. A., & Gerhard, O. 2020, MNRAS, 495, 3175, doi: 10.1093/mnras/staa1336
  * Sellwood & Wilkinson (1993) Sellwood, J. A., & Wilkinson, A. 1993, RPPh, 56, 173, doi: 10.1088/0034-4885/56/2/001
  * Sheth et al. (2008) Sheth, K., Elmegreen, D. M., Elmegreen, B. G., et al. 2008, ApJ, 675, 1141, doi: 10.1086/524980
  * Sheth et al. (2010) Sheth, K., Regan, M., Hinz, J. L., et al. 2010, PASP, 122, 1397, doi: 10.1086/657638
  * Skokos et al. (2002) Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847, doi: 10.1046/j.1365-8711.2002.05468.x
  * Smirnov et al. (2021) Smirnov, A. A., Tikhonenko, I. S., & Sotnikova, N. Y. 2021, MNRAS, 502, 4689, doi: 10.1093/mnras/stab327
  * Sparke & Sellwood (1987) Sparke, L. S., & Sellwood, J. A. 1987, MNRAS, 225, 653, doi: 10.1093/mnras/225.3.653
  * Stadel (2001) Stadel, J. G. 2001, Ph.D. Thesis, University of Washington
  * Stinson et al. (2006) Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, 373, 1074, doi: 10.1111/j.1365-2966.2006.11097.x
  * Valluri et al. (2010) Valluri, M., Debattista, V. P., Quinn, T., & Moore, B. 2010, MNRAS, 403, 525, doi: 10.1111/j.1365-2966.2009.16192.x
  * Valluri & Merritt (1998) Valluri, M., & Merritt, D. 1998, ApJ, 506, 686, doi: 10.1086/306269
  * Valluri et al. (2016) Valluri, M., Shen, J., Abbott, C., & Debattista, V. P. 2016, ApJ, 818, 141, doi: 10.3847/0004-637X/818/2/141
  * van Albada & Sanders (1982) van Albada, T. S., & Sanders, R. H. 1982, MNRAS, 201, 303, doi: 10.1093/mnras/201.2.303
  * Vasiliev (2019) Vasiliev, E. 2019, MNRAS, 482, 1525, doi: 10.1093/mnras/sty2672
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wadsley et al. (2004) Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New Astronomy, 9, 137, doi: 10.1016/j.newast.2003.08.004
  * Weinberg & Katz (2007) Weinberg, M. D., & Katz, N. 2007, MNRAS, 375, 425, doi: 10.1111/j.1365-2966.2006.11306.x
  * Wheeler et al. (2023) Wheeler, V., Valluri, M., Beraldo e Silva, L., Dattathri, S., & Debattista, V. P. 2023, arXiv e-prints, arXiv:2306.07424, doi: 10.48550/arXiv.2306.07424
  * Widrow & Dubinski (2005) Widrow, L. M., & Dubinski, J. 2005, ApJ, 631, 838, doi: 10.1086/432710