# Oscillating dynamical friction on galactic bars by trapped dark matter

Rimpei Chiba1 and Ralph Schönrich2
1Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, Parks
Road, Oxford OX1 3PU, UK
2Mullard Space Science Laboratory, University College London, Holmbury St.
Mary, Dorking, Surrey, RH5 6NT, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The dynamic evolution of galactic bars in standard $\Lambda$CDM models is
dominated by angular momentum loss to the dark matter haloes via dynamical
friction. Traditional approximations to dynamical friction are formulated
using the linearized collisionless Boltzmann equation and have been shown to
be valid in the fast regime, i.e. for rapidly slowing bars. However, the
linear assumption breaks down within a few dynamical periods for typical
slowly evolving bars, which trap a significant amount of disc stars and dark
matter in resonances. Recent observations of the Galactic bar imply this slow
regime. We explore the fundamental mechanism of dynamical friction in the slow
regime with analytical and test-particle methods. Here, angular momentum
exchange is dominated by resonantly trapped orbits which slowly librate around
the resonances. In typical equilibrium haloes, the initial phase-space density
within the trapped zone is higher at lower angular momentum. Since the
libration frequency falls towards the separatrix, this density contrast winds
up into a phase-space spiral, resulting in a dynamical friction that
oscillates with $\sim$Gyr periods and damps over secular timescales. We
quantify the long-term behaviour of this torque with secular perturbation
theory, and predict two observable consequences: i) The phase-space spirals
may be detectable in the stellar disc where the number of windings encodes the
age of the bar. ii) The torque causes oscillations in the bar’s pattern speed
overlaying the overall slowdown – while not discussed, this feature is visible
in previous simulations.

###### keywords:

Galaxy: kinematics and dynamics – Galaxy: evolution – methods: numerical

††pubyear: 2021††pagerange: Oscillating dynamical friction on galactic bars by
trapped dark matter–G

##  1 Introduction

Dark matter haloes strongly affect the evolution of galactic bars: dynamical
friction transfers angular momentum and energy from the galactic bar to the
dark matter halo, allowing bars to slow and grow. This angular momentum
transfer depends both on the density and on the kinematic state of the dark
halo. Consequently, we can constrain the otherwise poorly known phase-space
distribution and nature of dark haloes from the measured slowing rate of bars.

While a wide range of numerical simulations have addressed the angular
momentum transfer from bars to haloes (e.g. Hernquist & Weinberg, 1992;
Debattista & Sellwood, 2000; Athanassoula, 2003), the analytical treatment is
still incomplete. Tremaine & Weinberg (1984, hereafter TW84) derived a formula
of dynamical friction in an inhomogeneous spherical system using linear
perturbation theory. Their formula generalises the earlier LBK torque formula
(Lynden-Bell & Kalnajs, 1972), which describes the angular momentum transport
in two-dimensional disc galaxies. Unlike the dynamical friction in a
homogeneous system (Chandrasekhar, 1943), the LBK formula shows that in an
inhomogeneous system, where orbits are quasi-periodic, the dominant angular
momentum transfer arises from discrete resonances. TW84 also investigated the
non-linear motion of orbits near resonances and showed that the LBK formula is
only valid in the so-called fast regime. Here, fast means that the resonances
sweep across the phase space of the halo rapidly enough to preclude the growth
of non-linear responses, in particular trapping at resonances. This allows the
system to remain in the linear regime. Weinberg (1985) showed that the LBK
torque formula predicts a strong friction on the bar, extracting most of its
angular momentum in a few rotation periods. This problem was revisited by
Weinberg (2004, hereafter W04), where he pointed out that the assumption of
the perturbation being turned on adiabatically in the distant past (the time-
asymptotic limit) is problematic since the evolution timescale of the bar is
only several times the characteristic dynamical period. Relaxing this
assumption enabled W04 to model the transient effect and reduce the overall
torque. Yet, the W04 treatment remains restricted to the fast regime and
cannot be applied to the slow regime where resonant trapping comes into play.

This slow regime (resonant trapping by the Galactic bar) is, however, implied
by the stellar kinematics of the Solar neighbourhood observed by the Gaia
satellite (e.g. Katz et al., 2019; Hunt et al., 2019; Monari et al., 2019;
Binney, 2020; Trick et al., 2021). The resonantly trapped stars manifest as
stellar streams which occupy a sizeable fraction of the local phase-space
volume. Modelling particularly the Hercules stream, Chiba et al. (2021)
inferred that the bar is decelerating at a moderate rate keeping a good
portion of local stars trapped. Furthermore, signatures of trapping is found
in the metallicity distribution of local stars: the corotation resonance of
the slow/long bar exhibit a monotonic increase of metallicity towards the core
of the resonance, indicating a tree-ring like growth as predicted by a
naturally slowing bar that sequentially captured stars from small
Galactocentric radii (Chiba & Schönrich, 2021).

Resonant trapping by the bar is also commonly seen in self-consistent N-body
simulations: spectral analysis of orbits in barred galaxies shows significant
clustering near resonances (e.g. Ceverino & Klypin, 2007). Simulations by
Halle et al. (2018) show that these trapped stars adhere to the resonance as
the bar slows down. Some works have also discussed direct trapping of dark
matter into a ‘shadow bar’ (Athanassoula, 2007; Petersen et al., 2016; Collier
& Madigan, 2021).

Thus both observations and simulations indicate that bars generally evolve
slowly enough to exert resonant trapping and make the formalism of TW84 and
W04 inapplicable. The dynamics of this resonant trapping is best viewed in the
slow angle-action plane which exhibits a phase flow akin to that of a simple
pendulum: trapped orbits librate around the resonance centre (and their slow
angle oscillates around the resonance midpoint) while untrapped orbits
circulate above and below the separatrix (i.e. their slow angle moves freely
across the whole $2\pi$ range). Before bar formation, the phase-space density
typically declines towards large slow action (i.e. the $z$-angular momentum).
When the bar forms, the trapped region grows and captures orbits from both
above and below the resonance, reconnecting phase space across regions of
different density. The newly trapped orbits are thus inevitably non-uniform in
their libration angle distribution. As these trapped orbits librate around the
resonance along contours of constant libration actions, their innate angle
imbalance results in a periodic dynamical friction on the bar. However, due to
the monotonic decrease of libration frequency towards the separatrix, the
density inside the trapped zone gradually winds up into a phase-space spiral,
and the net torque from trapped orbits slowly attenuates. The phase mixing
timescale is set by the libration period (typically $>1\,{\rm Gyr}$), which
prevents complete mixing in a Hubble time. This leads us to predict that bars
are subject to pattern speed oscillations on the timescale of the libration,
which is an order of magnitude longer than the short-period fluctuations due
to bar-spiral interactions (Wu et al., 2016; Hilmi et al., 2020).

In this paper, we study this behaviour both analytically and in test particle
simulations. The analytical approach uses the standard secular perturbation
theory founded on the method of averaging. Our formulation is initially done
for a constantly rotating bar. The resulting torque terms provide a necessary,
but still incomplete step towards a fully self-consistent model where the bar
slows in response to the torque by the halo. For this, one has to further
address the processes of dragging by the resonance, jumping over the
resonance, as well as capture into and escape from it (e.g. Henrard, 1982;
Sridhar & Touma, 1996).

This paper is organized as follows. In section 2, we introduce our model of
the halo and the bar. Section 3 reviews the torque formula in the fast regime
and discuss its problems when applied to systems in the slow regime. In
section 4, we lay out the standard method to treat the non-linear dynamics
near resonance and formulate the dynamical friction in the slow regime which
we verify against test-particle simulations. We will discuss the consequence
of an evolving bar in section 5 and sum up in section 6.

##  2 Model

###  2.1 Model of dark halo

We model the dark halo with an isotropic Hernquist (1990) model:

$\displaystyle\rho_{0}(r)$ $\displaystyle=\frac{M}{2\pi}\frac{r_{\rm
s}}{r(r_{\rm s}+r)^{3}},$ (1) $\displaystyle\Phi_{0}(r)$
$\displaystyle=-\frac{GM}{r_{\rm s}+r},$ (2)

where $G$ is the gravitational constant, $M$ is the total halo mass, and
$r_{\rm s}$ the halo scale radius. By default, we set $M=1.5\times
10^{12}\,{\rm M}_{\odot}$ and $r_{\rm s}=20\,{\rm kpc}$. The corresponding
distribution function and setup of the test-particle simulations are described
in Appendix A.

###  2.2 Model of Galactic bar

We model the bar potential with a quadrupole rotating at pattern speed
$\Omega_{\rm p}$:

$\displaystyle\Phi_{\rm b}(r,\vartheta,\varphi,t)=\Phi_{\rm
b}(r)\sin^{2}\vartheta\cos 2\left(\varphi-\Omega_{\rm p}t\right),$ (3)

where $(r,\vartheta,\varphi)$ are the usual spherical coordinates. The radial
dependence of the bar potential is

$\displaystyle\Phi_{\rm b}(r)=-\frac{Av_{\rm c}^{2}}{2}\left(\frac{r}{r_{\rm
CR}}\right)^{2}\left(\frac{b+1}{b+r/r_{\rm CR}}\right)^{5},$ (4)

where $A$ is the dimensionless strength of the bar and $b$ describes the ratio
between the bar scale length and the corotation radius $r_{\rm CR}$. This
potential is designed to lengthen as the bar slows. We fit our model to that
of Sormani et al. (2015) at the Galactic plane ($\vartheta=\pi/2$) and obtain
$A=0.02,b=0.28$. By default, we use the current bar pattern speed $\Omega_{\rm
p}=35\,{\rm km}\,{\rm s}^{-1}\,{\rm kpc}^{-1}\simeq 35.8\,{\rm Gyr}^{-1}$
(Binney, 2020; Chiba & Schönrich, 2021; Clarke & Gerhard, 2021).

###  2.3 Angle-action variables in spherical potential

The description of orbits in an integrable system is greatly simplified with
the use of angle-action variables $({\bm{\theta}},{\bm{J}})$ which are
canonical coordinates on phase space (e.g. Binney & Tremaine, 2008). By
construction, the Hamiltonian $H_{0}({\bm{J}})$ does not depend on the angles,
so the equations of motion are trivial in these coordinates:

$\displaystyle\dot{{\bm{J}}}=-\frac{{\partial}H_{0}}{{\partial}{\bm{\theta}}}=0,~{}~{}~{}~{}~{}~{}\dot{{\bm{\theta}}}=\frac{{\partial}H_{0}}{{\partial}{\bm{J}}}\equiv\bm{\Omega}({\bm{J}}),$
(5)

i.e., the actions are conserved and the conjugate angles increase linearly
with time at constant rate $\bm{\Omega}({\bm{J}})$.

In a spherically symmetric potential, a convenient set of actions are
${\bm{J}}=(J_{r},L,L_{z})$ where $J_{r}$ is the radial action, $L$ is the
total angular momentum, and $L_{z}$ is the $z$-component of $L$. The ratio
between $L$ and $L_{z}$ defines the orbital inclination
$\beta\equiv\cos^{-1}\left(L_{z}/L\right)$. The conjugate angle variables
${\bm{\theta}}=(\theta_{r},\theta_{\psi},\theta_{\varphi})$ describe,
respectively, the radial phase, the azimuthal phase in the orbital plane, and
the (fixed) longitude of the ascending node, i.e. the intersection of the
orbital plane and the galactic equatorial plane.

##  3 Torque in the fast regime

To compare our approach (section 4) with the conventional approximations, we
here briefly review the torque formula in the fast regime developed by TW84
and W04. Following W04, we solve for the second-order torque using the linear
solution to the collisionless Boltzmann equation. We will demonstrate the
breakdown of the formula in the slow regime and discuss the cause of problem.
Note that throughout the paper we ignore the halo’s self-gravitational
perturbation.

We assume that the halo is in a steady state prior to bar formation. The
unperturbed distribution is then a function of the actions alone
$f_{0}({\bm{J}})$ (Jeans’ theorem). When the bar forms, the perturbing
potential changes the halo’s distribution by
$f_{1}({\bm{\theta}},{\bm{J}},t)$. We can then write the Hamiltonian $H$ and
the distribution function $f$ as

$\displaystyle
H({\bm{\theta}},{\bm{J}},t)=H_{0}({\bm{J}})+\Phi_{1}({\bm{\theta}},{\bm{J}},t),~{}f({\bm{\theta}},{\bm{J}},t)=f_{0}({\bm{J}})+f_{1}({\bm{\theta}},{\bm{J}},t),$
(6)

where $\Phi_{1}$ consists of the bar potential (3) and the self-gravitational
perturbation induced by the response of the halo $f_{1}$, although here we
ignore the latter. The angle variables are by construction $2\pi$ periodic, so
one can expand the perturbations into a Fourier series:

$\displaystyle\Phi_{1}({\bm{\theta}},{\bm{J}},t)=\sum_{{\mathbf{n}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)\mathrm{e}^{i{\mathbf{n}}\cdot{\bm{\theta}}},~{}~{}~{}f_{1}({\bm{\theta}},{\bm{J}},t)=\sum_{{\mathbf{n}}}\hat{f}_{{\mathbf{n}}}({\bm{J}},t)\mathrm{e}^{i{\mathbf{n}}\cdot{\bm{\theta}}},$
(7)

where ${\mathbf{n}}=({\rm n}_{r},{\rm n}_{\psi},{\rm n}_{\varphi})$ is a
vector of integers. $\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)$ is given in
Appendix B. The torque exerted on a single orbit is then

$\displaystyle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}$
$\displaystyle=-\frac{{\partial}H}{{\partial}\theta_{\varphi}}=-i\sum_{{\mathbf{n}}}{\rm
n}_{\varphi}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)\mathrm{e}^{i{\mathbf{n}}\cdot{\bm{\theta}}},$
(8)

and the total torque integrated over the full distribution is

$\displaystyle\left\langle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\right\rangle=\int\mathrm{d}^{3}{\bm{J}}\mathrm{d}^{3}{\bm{\theta}}f({\bm{\theta}},{\bm{J}},t)\frac{\mathrm{d}L_{z}}{\mathrm{d}t}$
$\displaystyle=\\!\int\\!\mathrm{d}^{3}{\bm{J}}\mathrm{d}^{3}{\bm{\theta}}\left[f_{0}({\bm{J}})+\sum_{\mathbf{n}}\hat{f}_{{\mathbf{n}}}({\bm{J}},t)\mathrm{e}^{i{\mathbf{n}}\cdot{\bm{\theta}}}\right]\left[-i\sum_{{\mathbf{n}}^{\prime}}{\rm
n}_{\varphi}^{\prime}\hat{\Phi}_{{\mathbf{n}}^{\prime}}({\bm{J}},t)\mathrm{e}^{i{\mathbf{n}}^{\prime}\cdot{\bm{\theta}}}\right].$
(9)

The first term vanishes upon integration over ${\bm{\theta}}$ and the second
term is non-zero only when ${\mathbf{n}}^{\prime}=-{\mathbf{n}}$, so

$\displaystyle\left\langle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\right\rangle=i\left(2\pi\right)^{3}\sum_{\mathbf{n}}{\rm
n}_{\varphi}\int\mathrm{d}^{3}{\bm{J}}\hat{f}_{{\mathbf{n}}}({\bm{J}},t)\hat{\Phi}_{{\mathbf{n}}}^{*}({\bm{J}},t),$
(10)

where we used $\hat{\Phi}_{-{\mathbf{n}}}=\hat{\Phi}_{{\mathbf{n}}}^{*}$ as
required by ${\rm Im}[\Phi_{1}]=0$. We see from the above equation that the
torque is second order in the perturbation. To proceed, we need to develop the
response $\hat{f}_{{\mathbf{n}}}({\bm{J}},t)$ as a function of the potential
perturbation $\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)$. This is described by the
collisionless Boltzmann equation (CBE):

$\displaystyle\frac{\mathrm{d}f}{\mathrm{d}t}=\frac{{\partial}f}{{\partial}t}+\frac{{\partial}f}{{\partial}{\bm{\theta}}}\cdot\frac{{\partial}H}{{\partial}{\bm{J}}}-\frac{{\partial}f}{{\partial}{\bm{J}}}\cdot\frac{{\partial}H}{{\partial}{\bm{\theta}}}=0.$
(11)

Using the Fourier decomposition for both $f$ and $H$, we have, for each
${\mathbf{n}}$,

$\displaystyle\frac{{\partial}\hat{f}_{{\mathbf{n}}}}{{\partial}t}$
$\displaystyle+i{\mathbf{n}}\cdot\bm{\Omega}\hat{f}_{{\mathbf{n}}}-i{\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\hat{\Phi}_{{\mathbf{n}}}$
$\displaystyle+\sum_{{\mathbf{n}}^{\prime}}\left[i({\mathbf{n}}-{\mathbf{n}}^{\prime})\cdot\frac{{\partial}\hat{\Phi}_{{\mathbf{n}}^{\prime}}}{{\partial}{\bm{J}}}\hat{f}_{{\mathbf{n}}-{\mathbf{n}}^{\prime}}-i{\mathbf{n}}^{\prime}\cdot\frac{{\partial}\hat{f}_{{\mathbf{n}}-{\mathbf{n}}^{\prime}}}{{\partial}{\bm{J}}}\hat{\Phi}_{{\mathbf{n}}^{\prime}}\right]=0.$
(12)

The standard prescription to solve this for $\hat{f}_{{\mathbf{n}}}$ is to
ignore the non-linear terms inside the square bracket assuming that both
$|f_{1}/f_{0}|$ and $|\Phi_{1}/H_{0}|$ are sufficiently small. This yields the
linearized CBE, which has the following solution for an external perturbation
imposed after $t=0$:

$\displaystyle\hat{f}_{{\mathbf{n}}}({\bm{J}},t)=i{\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-i{\mathbf{n}}\cdot\bm{\Omega}(t-t^{\prime})}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t^{\prime}).$
(13)

Substituting equation (13) to (10) yields the general torque formula

$\displaystyle\left\langle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\right\rangle=$
$\displaystyle-\left(2\pi\right)^{3}\sum_{\mathbf{n}}{\rm
n}_{\varphi}\int\mathrm{d}^{3}{\bm{J}}{\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}$
$\displaystyle\hskip
28.45274pt\times\left[\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-i{\mathbf{n}}\cdot\bm{\Omega}(t-t^{\prime})}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t^{\prime})\right]\hat{\Phi}_{{\mathbf{n}}}^{*}({\bm{J}},t).$
(14)

For a perturbation with constant amplitude and pattern speed $\Omega_{\rm p}$,
the time dependence separates as
$\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)=\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\mathrm{e}^{-i{\rm
n}_{\varphi}\Omega_{\rm p}t}$, and we have (W04)

$\displaystyle\left\langle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\right\rangle\\!=$
$\displaystyle\\!-\\!\left(2\pi\right)^{3}\\!\sum_{\mathbf{n}}{\rm
n}_{\varphi}\\!\int\\!\mathrm{d}^{3}\\!{\bm{J}}{\mathbf{n}}\\!\cdot\\!\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}|\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})|^{2}\frac{\sin\left[\left({\mathbf{n}}\cdot\bm{\Omega}-{\rm
n}_{\varphi}\Omega_{\rm p}\right)t\right]}{{\mathbf{n}}\cdot\bm{\Omega}-{\rm
n}_{\varphi}\Omega_{\rm p}},$ (15)

which, in the time-asymptotic limit $t\rightarrow\infty$, reduces to the LBK
formula (Lynden-Bell & Kalnajs, 1972, TW84)

$\displaystyle\left\langle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\right\rangle=$
$\displaystyle-\left(2\pi\right)^{3}\sum_{\mathbf{n}}{\rm
n}_{\varphi}\\!\int\mathrm{d}^{3}{\bm{J}}{\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}|\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})|^{2}\pi\delta\left({\mathbf{n}}\cdot\bm{\Omega}-{\rm
n}_{\varphi}\Omega_{\rm p}\right).$ (16)

W04 and more recently Banik & van den Bosch (2021) emphasized the problem of
taking the time-asymptotic limit as the age of the galaxy (and yet less that
of the perturbation) may not be longer enough than the relevant dynamical
time.

It is obvious from equations (15) and (16) that the torque is dominated by
contribution from phase space near resonances, i.e.
${\mathbf{n}}\cdot\bm{\Omega}-{\rm n}_{\varphi}\Omega_{\rm p}\approx 0$.
However, the vicinity of resonances is precisely where the linearized CBE
becomes invalid: the linear response (13) grows indefinitely at resonances
(Appendix D), and this will soon violate the assumption that the non-linear
terms of the CBE (12) are negligible, which is the basis of this formalism. A
correct physical description must respect that the orbits near resonances
become trapped, a non-linear secular behaviour that cannot be described by
linear theory.

Figure 1: Total torque on the halo (divided by the total halo mass) perturbed
by a bar with constant pattern speed and amplitude switched on at $t=0$.
Linear theory (blue, equation 15) agrees with the torque of the test-particle
simulation (black) for a few dynamical period (top panel) but fails to predict
the long-term evolution (bottom panel). The LBK formula (green, equation 16)
gives the time-asymptotic value of linear theory.

To illustrate the problem, we show in Fig. 1 the torque exerted on the halo by
a constantly rotating bar. We compare the torque measured from a test particle
simulation (black) with linear theory (equation 15, blue) and the LBK
approximation (equation 16, dashed green). For the two approximations, we
restrict the summation to the indices with ${\rm n}_{r}\in[-10,10]$, ${\rm
n}_{\psi}\in-2,0,2$, and ${\rm n}_{\varphi}=2$. Since both the halo density
and the perturbing potential decrease towards large radii, the contributions
are largely limited to the main resonances with small $r$, the strongest being
the corotation resonance ${\mathbf{n}}=(0,2,2)$. The top panel shows the
evolution for the first $1\,{\rm Gyr}$ which we see a fair agreement between
simulation (black) and time-dependent linear theory (blue) as long as $t$ is
sufficiently small compared to the libration periods (of order $1-2\,{\rm
Gyr}$). However, the linear theory qualitatively fails to predict the
subsequent long-term evolution shown in the bottom panel: while the torque of
linear theory rapidly converges to a non-zero constant value predicted by the
LBK formula, the ‘true’ torque fluctuates over a much longer timescale
gradually settling down to zero111We note that this limit depends on the time-
dependence of the perturbation. Here the approach to a steady state (zero
torque) is the result of the assumed constant bar pattern speed and
amplitude.. This relaxation process is driven by phase mixing (a collisionless
relaxation), which we describe in the next section.

In this experiment, we have kept the bar’s pattern speed constant. If the bar
decelerates sufficiently fast that the resonances pass over the orbits before
they can respond non-linearly, then resonant trapping will not occur and the
formula will predict the right torque for an extended period of time as
demonstrated in W04. For this reason, the formula is said to be valid in the
fast regime ( TW84, ). The validity of the formula is characterized by the
speed parameter222In the limit of epicycle approximation, the speed parameter
reduces to $s\sim\eta/A$ where $\eta\equiv-\dot{\Omega}_{\rm p}/\Omega_{\rm
p}^{2}$ is the bar’s dimensionless slowing rate and $A$ the dimensionless bar
strength (Chiba et al., 2021, equation 27).( TW84, equation 89):

$\displaystyle s\equiv\frac{\left|{\rm n}_{\varphi}\dot{\Omega}_{\rm
p}\right|}{\omega_{0}^{2}},$ (17)

where $\dot{\Omega}_{\rm p}$ measures the secular evolution of the pattern
speed (i.e. not the short-term fluctuations as discussed in Wu et al., 2016)
and $\omega_{0}$ is the libration frequency of trapped orbits at the core of
the resonance (equation 25). In the fast regime $s>1$, orbits cannot stay
trapped in resonance because the resonant potential no longer features a local
minimum (TW84, Chiba et al., 2021, see also Appendix G).

If the perturbation vanishes before any non-linear effect develops, the
formula will be valid regardless of $s$. This could be the case for e.g.
transient spiral arms. Hence, the traditional torque formulas (14)-(16) are
appropriate to study the evolution of perturbations that either decelerate
rapidly or decay rapidly. Since neither applies to long-lived galactic bars,
we must fully model the non-linear response of the halo, i.e. trapping at
resonances.

##  4 Torque in the slow regime

###  4.1 Angle-action variables near resonances

The alternative to solving the non-linear terms in the CBE (12) is to change
the coordinates. The majority of orbits near each resonance
${\bm{N}}=(N_{r},N_{\psi},N_{\varphi})$ behave quasi-periodically regardless
of trapping and thus there exist a new set of angle-action coordinates for
them (e.g. TW84; Kaasalainen, 1994; Sridhar & Touma, 1996; Binney, 2016, 2018;
Monari et al., 2017). In these ‘true’ angle-action coordinates, the time
evolution of the distribution function can be modelled straightforwardly from
the non-linearized CBE (section 4.2). In this section, we first lay out the
known coordinate transformations, and from there develop our equations for the
torque on both trapped and untrapped orbits.

We start with a canonical transformation to the so-called slow-fast angle-
action variables $({\bm{\theta}}^{\prime},{\bm{J}}^{\prime})=(\theta_{{\rm
f}_{1}},\theta_{{\rm f}_{2}},\theta_{\rm s},J_{{\rm f}_{1}},J_{{\rm
f}_{2}},J_{\rm s})$ (e.g. TW84, ):

$\displaystyle\theta_{{\rm f}_{1}}=\theta_{r},\hskip 41.25648pt\theta_{{\rm
f}_{2}}=\theta_{\psi},\hskip 36.98857pt\theta_{\rm
s}={\bm{N}}\cdot{\bm{\theta}}-N_{\varphi}\Omega_{\rm p}t,$ (18) $\displaystyle
J_{{\rm f}_{1}}=J_{r}-\frac{N_{r}}{N_{\varphi}}L_{z},\hskip 8.53581ptJ_{{\rm
f}_{2}}=L-\frac{N_{\psi}}{N_{\varphi}}L_{z},\hskip 8.53581ptJ_{\rm
s}=\frac{L_{z}}{N_{\varphi}},$ (19)

using the generating function

$\displaystyle W({\bm{\theta}},{\bm{J}}^{\prime},t)=\theta_{r}J_{{\rm
f}_{1}}+\theta_{\psi}J_{{\rm
f}_{2}}+\left({\bm{N}}\cdot{\bm{\theta}}-N_{\varphi}\Omega_{\rm
p}t\right)J_{\rm s}.$ (20)

The Hamiltonian (6) transforms to

$\displaystyle H({\bm{\theta}}^{\prime},{\bm{J}}^{\prime})$
$\displaystyle=H_{0}({\bm{J}}^{\prime})-N_{\varphi}\Omega_{\rm p}J_{\rm
s}+\sum_{{\bm{k}}}\hat{\Psi}_{{\bm{k}}}({\bm{J}}^{\prime})~{}\mathrm{e}^{i{\bm{k}}\cdot{\bm{\theta}}^{\prime}},$
(21)

where we have expanded the perturbation $H_{1}$ into a Fourier series
(Appendix B) with indices ${\bm{k}}=(k_{{\rm f}_{1}},k_{{\rm f}_{2}},k_{\rm
s})$. The purpose of transforming to the slow-fast variables is to separate
the motion into slow and fast components: near the resonance, the slow angle
$\theta_{\rm s}$ evolves much slower than the fast angles ${\bm{\theta}}_{\rm
f}=(\theta_{{\rm f}_{1}},\theta_{{\rm f}_{2}})$. This allows one to average
the Hamiltonian over ${\bm{\theta}}_{\rm f}$ while holding $\theta_{\rm s}$
fixed which removes all perturbation terms in (21) with ${\bm{k}}_{\rm
f}=(k_{{\rm f}_{1}},k_{{\rm f}_{2}})\neq{\bm{0}}$. Since only terms with
$k_{\rm s}=\pm 1$ are non-zero for resonances with $N_{\varphi}=2$ (Appendix
B), we may write

$\displaystyle\bar{H}(\theta_{\rm s},{\bm{J}}^{\prime})$
$\displaystyle=H_{0}({\bm{J}}^{\prime})-N_{\varphi}\Omega_{\rm p}J_{\rm
s}+\Psi({\bm{J}}^{\prime})\cos\left(\theta_{\rm s}-\theta_{\rm s,res}\right),$
(22)

where $\Psi\equiv 2|\hat{\Psi}_{(0,0,1)}|$. In this averaged system, the two
fast actions ${\bm{J}}_{\rm f}=(J_{{\rm f}_{1}},J_{{\rm f}_{2}})$ are
conserved, i.e. $\dot{\bm{J}}_{\rm
f}=-\frac{\partial\bar{H}}{\partial{\bm{\theta}}_{\rm f}}={\bm{0}}$. Thus this
Hamiltonian is integrable: it has three isolating integrals, the two fast
actions and the Hamiltonian itself, meaning that angle-action coordinates for
$\bar{H}$ exist.

A convenient analytical transformation to the angle-action coordinates of
$\bar{H}$ is available if we Taylor expand $\bar{H}$ around the resonance
$J_{\rm s,res}$ up to second order and ignore terms smaller than
$|\Psi/H_{0}|$:

$\displaystyle\bar{H}(\theta_{\rm
s},{\bm{J}}^{\prime})\\!\simeq\\!\frac{1}{2}G\\!\left({\bm{J}}_{\rm f},J_{\rm
s,res}\right)\\!\left(J_{\rm s}\\!-\\!J_{\rm
s,res}\right)^{2}+\Psi\\!\left({\bm{J}}_{\rm f},J_{\rm
s,res}\right)\cos\left(\theta_{\rm s}\\!-\\!\theta_{\rm s,res}\right),$ (23)

where

$\displaystyle
G\left({\bm{J}}^{\prime}\right)\equiv\frac{\partial^{2}H_{0}}{\partial J_{\rm
s}^{2}}=\frac{\partial\Omega_{s}}{\partial J_{\rm s}}.$ (24)

The first order derivative of $H_{0}({\bm{J}}^{\prime})-N_{\varphi}\Omega_{\rm
p}J_{\rm s}$ vanishes at the resonance and we have dropped the constant terms.
The approximated Hamiltonian takes the form of a classical pendulum
Hamiltonian apart from a different sign. In galactic dynamics, $G$ (equation
24) is in many cases negative ( TW84, Figure 5) so in what follows we will
assume $G<0$ although cases with $G>0$ can be treated similarly by adjusting
$\theta_{\rm s,res}$ as appropriate. For subsequent use, we define

$\displaystyle\omega_{0}^{2}\equiv-G\Psi,~{}~{}~{}J_{0}^{2}\equiv-\Psi/G,~{}~{}~{}\varepsilon\equiv\sqrt{\frac{1}{2}\left(1-\frac{\bar{H}}{\Psi}\right)}.$
(25)

Depending on the value of $\bar{H}$, parametrized by the dimensionless energy
of the pendulum $\varepsilon$, there are two types of qualitatively different
motions as depicted in Fig. 2: orbits librate when $0\leq\varepsilon<1$
(trapped) and circulate above or below the resonance when $\varepsilon>1$
(untrapped). $\varepsilon=1$ defines the separatrix.

Figure 2: Phase flow for the approximate $\bar{H}$ (23). Level curves of
$\bar{H}$ parametrized by $\varepsilon$ are plotted. For $0\leq\varepsilon<1$,
the motion of orbits in $\theta_{\rm s}$ is bounded (libration) whereas for
$\varepsilon>1$, $\theta_{\rm s}$ is unbounded (circulation).

In the following, we summarize the canonical transformation from the slow
angle-action variables $(\theta_{\rm s},J_{\rm s})$ to the angle-action
coordinates of libration/circulation (Brizard, 2013), and formulate the
secular torque on trapped/untrapped orbits.

####  4.1.1 Trapped orbits ($0\leq\varepsilon<1$)

The angle-action coordinates of libration are

$\displaystyle\theta_{\ell}=\theta_{\ell
0}+\Omega_{\ell}t,~{}~{}\Omega_{\ell}=\frac{\pi}{2K(\varepsilon)}\omega_{0},$
(26) $\displaystyle J_{\ell}=\frac{1}{2\pi}\oint\mathrm{d}\theta_{\rm s}J_{\rm
s}=\frac{8J_{0}}{\pi}\left[E(\varepsilon)-\left(1-\varepsilon^{2}\right)K(\varepsilon)\right],$
(27)

where $K(\varepsilon)$ and $E(\varepsilon)$ are the complete elliptic integral
of the first and second kinds (Appendix E). $\omega_{0}$ and $J_{0}$ are
defined in equation (25). Note that $J_{0}>0$ but $\omega_{0}<0$. The inverse
transformation is

$\displaystyle\theta_{\rm s}-\theta_{\rm s,res}=2\arcsin\left[\varepsilon{\rm
sn}\left(\left.\frac{2K(\varepsilon)}{\pi}\theta_{\ell}\right|\varepsilon\right)\right],$
(28) $\displaystyle J_{\rm s}-J_{\rm s,res}=2J_{0}\varepsilon{\rm
cn}\left(\left.\frac{2K(\varepsilon)}{\pi}\theta_{\ell}\right|\varepsilon\right),$
(29)

where sn, cn and dn are the Jacobi elliptic functions (Appendix E). The torque
exerted on trapped orbits is then

$\displaystyle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}$
$\displaystyle=N_{\varphi}\Omega_{\ell}\frac{\mathrm{d}J_{\rm
s}}{\mathrm{d}\theta_{\ell}}=2N_{\varphi}\varepsilon\Psi{\rm
sn}\\!\left(\\!\left.\frac{2K(\varepsilon)}{\pi}\theta_{\ell}\right|\varepsilon\right){\rm
dn}\\!\left(\\!\left.\frac{2K(\varepsilon)}{\pi}\theta_{\ell}\right|\varepsilon\right),$
(30)

which becomes harmonic near the centre of the resonance:

$\displaystyle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\simeq
2N_{\varphi}\varepsilon\Psi\sin\theta_{\ell}~{}~{}~{}~{}{\rm
for}~{}~{}~{}~{}\varepsilon\ll 1.$ (31)

####  4.1.2 Non-trapped orbits ($\varepsilon>1$)

The angle-action coordinates of circulation are

$\displaystyle\theta_{\rm c}=\theta_{\rm c0}+\Omega_{\rm c}t,~{}~{}\Omega_{\rm
c}=\frac{\pi\varepsilon}{K(\varepsilon^{-1})}\omega_{0},$ (32) $\displaystyle
J_{\rm c}=\frac{1}{2\pi}\oint\mathrm{d}\theta_{\rm s}J_{\rm s}=J_{\rm
s,res}\pm\frac{4J_{0}\varepsilon}{\pi}E(\varepsilon^{-1}),$ (33)

where the $\pm$ signs correspond to the upper and lower circulating regimes,
respectively. Note that the subscript ‘c’ denotes ‘circulating’ (untrapped)
orbits and not ‘circular’ orbits. The inverse transformation is

$\displaystyle\theta_{\rm s}-\theta_{\rm s,res}=\pm 2\arcsin\left[{\rm
sn}\left(\left.\frac{K(\varepsilon^{-1})}{\pi}\theta_{\rm
c}\right|\varepsilon^{-1}\right)\right],$ (34) $\displaystyle J_{\rm s}-J_{\rm
s,res}=\pm 2J_{0}\varepsilon{\rm
dn}\left(\left.\frac{K(\varepsilon^{-1})}{\pi}\theta_{\rm
c}\right|\varepsilon^{-1}\right).$ (35)

In the limit $\varepsilon\rightarrow\infty$, the slow action is conserved and
the slow angle $\theta_{\rm s}$ asymptotically approaches the angle of
circulation $\theta_{\rm c}$. The torque exerted on untrapped orbits is

$\displaystyle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}$ $\displaystyle=\pm
2N_{\varphi}\Psi{\rm sn}\left(\left.\frac{K(\varepsilon^{-1})}{\pi}\theta_{\rm
c}\right|\varepsilon^{-1}\right){\rm
cn}\left(\left.\frac{K(\varepsilon^{-1})}{\pi}\theta_{\rm
c}\right|\varepsilon^{-1}\right).$ (36)

Analogous to (31), far from the resonance, the torque becomes harmonic:

$\displaystyle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\simeq\pm
N_{\varphi}\Psi\sin\theta_{\rm c}~{}~{}~{}~{}{\rm
for}~{}~{}~{}~{}\varepsilon\gg 1.$ (37)

The behaviour of the torque is plotted in Fig. 3. As seen in (31) and (37),
the torque fluctuates sinusoidally near the core $(\varepsilon\ll 1)$ of the
resonance, as well as far from it $(\varepsilon\gg 1)$. As predicted from
(31), the torque on trapped orbits diminishes towards $\varepsilon=0$ where it
vanishes. Near the separatrix $(\varepsilon\approx 1)$, the torque forms a dip
at $\theta_{\ell}=\pi/2,3\pi/2$ since orbits stagnate near the unstable
singular points $(\theta_{\rm s}-\theta_{\rm s,res},J_{\rm s}-J_{\rm
s,res})=(\pm\pi,0)$: in the case of the CR at $J_{{\rm f}_{1}}=0$ (circular,
$J_{r}=0$) and $J_{{\rm f}_{2}}=0$ (in-plane, $L=L_{z}$), these points
correspond to the unstable Lagrange points on the bar’s major axis. As the
orbit approaches the separatrix, the dip deepens and eventually the torque
becomes zero over the entire phase except at $\theta_{\ell}=0$ and $\pi$. This
situation is analogous to a real pendulum in transition from oscillation to
circulation that spends most of it time/phase stagnating at the top of the
pivot with hardly any change in momentum.

Figure 3: Torque applied to trapped $(\varepsilon<1)$ and untrapped
$(\varepsilon>1)$ orbits as a function of the libration/circulation angle
$\theta_{\ell/{\rm c}}$. The behaviour deviates from a sinusoidal wave near
the separatrix $\varepsilon=1$.

###  4.2 Evolution of the distribution function

The averaged Hamiltonian $\bar{H}$ (23) now only depends on the new actions
${\bm{J}}^{\prime\prime}=(J_{{\rm f}_{1}},J_{{\rm f}_{2}},J_{\ell/{\rm c}})$
that are conjugate to the new angles
${\bm{\theta}}^{\prime\prime}=(\theta_{{\rm f}_{1}},\theta_{{\rm
f}_{2}},\theta_{\ell/{\rm c}})$ where the subscript ‘$\ell/{\rm c}$’ indicates
libration or circulation depending on $\varepsilon$. The collisionless
Boltzmann equation (11) for the averaged system is then

$\displaystyle\frac{{\partial}f}{{\partial}t}+\bm{\Omega}^{\prime\prime}\cdot\frac{{\partial}f}{{\partial}{\bm{\theta}}^{\prime\prime}}=0,~{}~{}{\rm
where}~{}~{}~{}\bm{\Omega}^{\prime\prime}({\bm{J}}^{\prime\prime})\equiv\frac{{\partial}\bar{H}}{{\partial}{\bm{J}}^{\prime\prime}}.$
(38)

This equation has the trivial solution

$\displaystyle
f({\bm{\theta}}^{\prime\prime},{\bm{J}}^{\prime\prime},t)=f({\bm{\theta}}^{\prime\prime}-\bm{\Omega}^{\prime\prime}t,{\bm{J}}^{\prime\prime},0).$
(39)

One can thus readily obtain the phase space density at any time $t=\tau$ from
the original distribution at $t=0$ by simply winding the angles backward by
$\tau$. Since the initial distribution of the halo (43) is uniform in the fast
angles $(\theta_{{\rm f}_{1}},\theta_{{\rm f}_{2}})$ but not in the
libration/circulation angle $\theta_{\ell/{\rm c}}$, the system is out of
equilibrium after the emergence of the bar and starts phase mixing in
$\theta_{\ell/{\rm c}}$.

Figure 4: Time evolution of the phase-space density $f(\theta_{\rm s},J_{\rm
s})$ near the corotation resonance ${\bm{N}}=(0,2,2)$ at $(J_{{\rm
f}_{1}},J_{{\rm f}_{2}})=(10,0)\,{\rm kpc}^{2}\,{\rm Gyr}^{-1}$ plotted every
$1\,{\rm Gyr}$. Left column: test particle simulation. Right column:
analytical (equation 39). Black curves mark the separatrix.

Figure 4 shows the evolution of the phase-space density near the corotation
resonance (CR) at $(J_{{\rm f}_{1}},J_{{\rm f}_{2}})=(10,0)\,{\rm
kpc}^{2}\,{\rm Gyr}^{-1}$ plotted every $1\,{\rm Gyr}$ from top to bottom. The
bar is instantaneously switched on at $t=0$. The left-hand column shows the
distribution of the test-particle simulation which is well reproduced by the
analytical distribution shown on the right. Inside the separatrix (black
curves), the density, initially uniform in the slow angle, phase mixes in the
libration angle. Since the libration period increases towards the separatrix,
the phase space differentially rotates resulting in a phase-space spiral.
Similarly, the phase space outside the separatrix winds into stripes since the
circulation period is also maximal (infinite) at the separatrix and drops away
from there.

Figure 5: Phase-space density of the halo perturbed by a bar growing slowly in
$4\,{\rm Gyr}$. The distribution outside the resonance is smooth due to the
adiabatic emergence of the bar, while that inside the resonance displays a
winding spiral since the expanding separatrix simultaneously captures orbits
from above (low density) and below (high density).

The development of phase-space spiral inside the resonance does not hinge on
the instantaneous emergence of the bar. It also happens in a slowly growing
bar, as shown in Fig. 5, which displays the same phenomenon for a bar that
slowly grows over $4\,{\rm Gyr}$ following a polynomial prescription (Dehnen,
2000, equation 4). Far outside the separatrix, the change is generally
adiabatic since the circulation periods are small compared to the timescale on
which the potential changes. Here, the circulation action $J_{\rm c}$ is
conserved (see Appendix C) and so the distribution in the circulation angle is
kept uniform as it was at $t=0$, thus showing no stripes. However, the
situation is qualitatively different in the trapped region. The separatrix
simultaneously engulfs the high-density region from below and the low-density
region from above, thus rendering the distribution in the libration angle
inevitably nonuniform (almost like a rectangular wave distribution). Due to
this abrupt transition between maximum and minimum density, the phase-space
spiral in a growing bar has a hard boundary compared to the non-growing case
(Fig. 4) which shows a smoother transition between peaks and troughs.

The conservation of $J_{\rm c}$ also implies that the untrapped orbits bend
around the advancing separatrix, i.e. near the midpoint $\theta_{\rm
s}-\theta_{\rm s,res}=0$ they bend away from the resonance line and near
$\theta_{\rm s}-\theta_{\rm s,res}=\pm\pi$ they bend towards the resonance. In
the end (Fig. 5, bottom panel), the orbits just above and below the separatrix
originate from around $J_{\rm s}-J_{\rm s,res}\approx\pm 120\,{\rm
kpc}^{2}\,{\rm Gyr}^{-1}$, respectively, which leads to the discontinuous step
in density at the resonance line when moving along $J_{\rm s}$ at $\theta_{\rm
s}-\theta_{\rm s,res}=\pm\pi$.

Figure 6: Guiding orbits of particles in the galactic plane ($\beta=0$)
trapped in the CR of the bar (gray ellipse). Resonantly trapped orbits at
$\theta_{\rm s}-\theta_{\rm s,res}>0~{}(x<0)$ are pulled forward by the bar
(gain angular momentum), while those at $\theta_{\rm s}-\theta_{\rm
s,res}<0~{}(x>0)$ are pulled back by the bar (lose angular momentum). The net
torque on the bar is determined by the phase distribution of these trapped
dark matter. As they slowly phase mix, the net torque attenuates.

It is instructive to relate the behaviour in the slow angle-action plane with
the Chandrasekhar’s classical picture of dynamical friction, i.e. as an over-
density lagging behind the moving object (bar). At the CR, the slow angle (18)
is $\theta_{\rm s}=2[\theta_{\psi}-(\varphi_{\rm b}-\theta_{\varphi})]$ where
$\theta_{\psi}$ is the azimuthal angle of the guiding orbit in its orbital
plane, and $\varphi_{\rm b}-\theta_{\varphi}$ is the azimuthal angle of the
bar in the galactic plane measured from the ascending node. The centre of
oscillation of the CR is at the bar’s minor axis $\theta_{\psi}-(\varphi_{\rm
b}-\theta_{\varphi})=\pi/2$, so $\theta_{\rm s,res}=\pi$, and thus the range
$\theta_{\rm s}-\theta_{\rm s,res}\in[-\pi,\pi]$ corresponds to
$\theta_{\psi}-(\varphi_{\rm b}-\theta_{\varphi})\in[0,\pi]$. Figure 6 shows
typical in-plane $(\beta=0)$ orbits in the bar’s corotating frame. For
$\theta_{\rm s}-\theta_{\rm s,res}>0$, the trapped orbits are lagging behind
the bar, so they exert a negative torque on the bar, while for $\theta_{\rm
s}-\theta_{\rm s,res}<0$, they are ahead of the bar, exerting a positive
torque on the bar. The net torque is then determined by the distribution in
$\theta_{\rm s}$: In the beginning, the overdensity at small $z$-angular
momentum (or $J_{\rm s}$) librates towards $\theta_{\rm s}-\theta_{\rm
s,res}>0$ thus giving rise to a net negative torque on the bar. Subsequently,
this overdensity librates back and forth between $\theta_{\rm s}-\theta_{\rm
s,res}>0$ and $\theta_{\rm s}-\theta_{\rm s,res}<0$, providing a systematic
oscillation of dynamical friction, which will fade with phase mixing.
Eventually, the distribution in $\theta_{\rm s}$ approaches equilibrium and
the net torque on a constantly rotating bar converges to zero as we saw in
Fig. 1.

If the external perturbation is transient (e.g. short-lived spiral arms), once
the perturbation fades away, the slow angle will recover it role as a true
angle variable and thus the halo will resume phase mixing in $\theta_{\rm s}$.
However, since the libration angle $\theta_{\ell}$ at which the orbits got
captured and were released are generally different, the final DF near
resonances will be flatter than the original distribution. This ‘scar’ in the
DF will affect the evolution of the renewed spiral arms (Sellwood & Carlberg,
2014; Sridhar, 2019), and can cause self-gravitational instabilities, exciting
further spiral modes (Sellwood & Carlberg, 2019).

###  4.3 Onset of chaos by overlap of resonances

(a) Resonances with $N_{\psi}=2,N_{\varphi}=2$. $N_{r}$ spans from 0 (left
bottom) to 9 (top right).  (b) Resonances with $N_{\psi}=0,N_{\varphi}=2$.
$N_{r}$ spans from 1 (left bottom) to 10 (top right).  (c) Resonances with
$N_{\psi}=-2,N_{\varphi}=2$. $N_{r}$ spans from 2 (left bottom) to 11 (top
right).

Figure 7: Bar-halo resonances ${\bm{N}}=(N_{r},N_{\psi},N_{\varphi})$ at three
different orbital planes $\beta=\cos^{-1}(L_{z}/L)=0,\pi/2,\pi$ from left to
right. The thick black curves mark the location of resonances, the dotted
curves mark the maximum extent of the separatrix, and the thin black lines
mark the direction of libration. The colors show the phase-space density in
log scale. Top panel: Resonances with $N_{\psi}=2$. Resonantly trapped volume
is largest at $\beta=0$ and vanishes at $\beta=\pi$. The black symbols, placed
on the outer Lindblad resonance ${\bm{N}}=(1,2,2)$, mark the initial actions
of orbits shown in Fig. 8. Middle panel: Resonances with $N_{\psi}=0$. Trapped
volume is largest at $\beta=\pi/2$ and decays towards $\beta=0,\pi$. Bottom
panel: Resonances with $N_{\psi}=-2$. Trapped volume increase towards
$\beta=\pi$.

By averaging over the fast motions, we have approximated the Hamiltonian near
a resonance with an integrable Hamiltonian $\bar{H}$ (23). We now ask whether
the system remains integrable when the neglected terms of the full Hamiltonian
$\delta H=H-\bar{H}$ are added back in. To consider the effect of $\delta H$,
we Fourier expand it in the new angles

$\displaystyle\delta
H({\bm{\theta}}^{\prime\prime},{\bm{J}}^{\prime\prime})=\sum_{{\mathbf{m}}}\hat{\Psi}_{{\mathbf{m}}}({\bm{J}}^{\prime\prime})\mathrm{e}^{i{\mathbf{m}}\cdot{\bm{\theta}}^{\prime\prime}}.$
(40)

The Kolmogorov–Arnold–Moser (KAM) theorem (Arnol’d, 1963) implies that these
terms will only slightly perturb the motion and the system will stay
integrable so long as they are small and sufficiently far from satisfying a
resonance condition ${\mathbf{m}}\cdot\bm{\Omega}^{\prime\prime}=0$ where
$\bm{\Omega}^{\prime\prime}=(\Omega_{{\rm f}_{1}},\Omega_{{\rm
f}_{2}},\Omega_{\ell/{\rm c}})$. However, if $\delta H$ contains a resonant
term, it may render the system non-integrable. These secondary resonances form
new resonant islands (either inside or outside the original separatrix) and
can in principle be modelled in much the same way as we have dealt with the
original resonance of the unperturbed system $H_{0}$ (e.g. Lichtenberg &
Lieberman, 1992; Malhotra, 1998; Wisdom, 2004). However, when the separatrices
of the secondary resonances come sufficiently close to the separatrix of the
original resonances, orbits may behave stochastically, moving from one domain
of resonance to another (Chirikov, 1979). A significant secondary resonance is
expected when two resonances of the unperturbed Hamiltonian $H_{0}$ reside
close to each other333Outside the original separatrix, the circulation
frequency approximates the slow frequency and thus the resonance condition
${\mathbf{m}}\cdot\bm{\Omega}^{\prime\prime}=0$ in the resonant system
$\bar{H}({\bm{J}}^{\prime\prime})$ approximates the original resonance
condition ${\mathbf{n}}\cdot\bm{\Omega}-{\rm n}_{\varphi}\Omega_{\rm p}=0$ in
the unperturbed system $H_{0}({\bm{J}})$.. The resonances in a bar-halo system
indeed overlap at particular regions of phase space and some orbits there turn
chaotic/irregular.

Figure 8: Motion of orbits trapped at the OLR. The symbol on the left top
corner represent the position of the initial actions marked in Fig. 7. The
black orbits are integrated numerically using the full Hamiltonian and the
blue orbits are obtained analytically from the averaged Hamiltonian. The orbit
at the top-right panel exhibits chaotic oscillations as it is overlapped by
the ${\bm{N}}=(2,0,2)$ resonance. All orbits take similar values of
$\varepsilon\sim 0.9$. Figure 9: Both panels show the same orbit in the angle-
action plane of the ${\bm{N}}=(1,2,2)$ (OLR, left-hand) and the
${\bm{N}}=(2,0,2)$ resonance (right-hand). Colours encode the time. The orbit
is initially on the OLR and then passes to the ${\bm{N}}=(2,0,2)$ resonance.

To illustrate the overlapping of resonances, we show in Fig. 7 the
configuration of bar resonances with $N_{\varphi}=2$ and $N_{\psi}=2,0,-2$
from top to bottom. A figure overlaying the three resonances is given in
Appendix F. In a spherical system, the resonant condition only depends on $L$
and $J_{r}$. However, the volume of the trapped zone depends in addition on
the orbital inclination $\beta\equiv\cos^{-1}(L_{z}/L)$. We thus show three
plots from left to right with $\beta=0$ (prograde), $\beta=\pi/2$
(perpendicular), and $\beta=\pi$ (retrograde). Following Binney (2018), we
draw the furthest excursions of trapped orbits with dotted curves and mark the
direction of libration with thin black lines. The colors show the phase-space
density of the unperturbed distribution function (43).

The top panels of Fig. 7 show resonances with $N_{\psi}=2$ which are familiar
in disc dynamics: the resonance at the bottom-left corner is the corotation
resonance (CR, $N_{r}=0$) and the neighboring one is the outer Lindblad
resonance (OLR, $N_{r}=1$). In our model, resonances with negative $N_{r}$
including the inner Lindblad resonance (ILR, $N_{r}=-1$) have negligibly small
volume444This contrasts with the models of W04 which are dominated by the ILR.
The discrepancy is most likely due to the difference in the adopted bar
pattern speed: our pattern speed $\Omega_{\rm p}\simeq 35.8\,{\rm Gyr}$
corresponds to $R_{\rm ILR}\sim 0.25\,{\rm kpc}$ in our Hernquist halo,
whereas their initial $\Omega_{\rm p}$ is much lower which implies larger
$R_{\rm ILR}$ and thus more mass trapped in the ILR.. The trapped volume is
largest at $\beta=0$ (left, in-plane prograde), and decreases towards
$\beta=\pi$ (right, in-plane retrograde) where it vanishes. The panels in the
middle row present resonances with $N_{\psi}=0$. Orbits trapped in these
resonances complete $N_{\varphi}$ radial oscillations in $N_{r}$ bar period
regardless of their azimuthal phase with respect to the bar. Weinberg & Katz
(2007) called these resonances the direct radial resonances (DRRs). The
strength/size of DRRs peaks at $\beta=\pi/2$. The lowest panels show
resonances with $N_{\psi}=-2$ which are strongest at $\beta=\pi$.

Comparing the different sets of resonances, we can see from Fig. 7 that the
resonances partly overlap especially near $\beta\approx\pi/2$, where all
resonances fill a finite phase-space volume (see Fig. 18 for a directly
superimposed plot). To see the impact of resonance overlap, we plot in Fig. 8
four orbits trapped in the OLR ${\bm{N}}=(1,2,2)$ with different initial
actions marked by black points in Fig. 7. As expected, the two orbits at
$\beta=0$ (circle and rhombus, left column in Fig. 8) are well modelled by the
averaged Hamiltonian (blue), since only resonances with $N_{\psi}=2$ has a
non-zero volume. The right-hand column of Fig. 8 (triangle and square) shows
two contrasting examples for $\beta=\pi/2$: the bottom right orbit (square) is
mildly distorted as it is near, but still outside the second resonance, while
the top right orbit (triangle) shows marked excursions as it is within the
region of the ${\bm{N}}=(2,0,2)$ resonance. We note that we ascribe only the
large-scale distortions to the presence of the second resonance, e.g. the
irregular dip near $\theta_{\rm s}-\theta_{\rm res}\approx 0$. The step-like
behaviour seen in the figure comes from the high eccentricity of the orbit
(large $J_{r}$ and small $L$), i.e. the fast change in $L_{z}$ when an orbit
rapidly passes near the bar at its pericentre.

When the domains of two resonances overlap, orbits may also pass from one
resonance to another. Figure 9 shows an example of such an orbit moving
between the OLR and the $(2,0,2)$ resonance. When viewed in the slow plane of
the OLR (left panel), the orbit is initially trapped and then escapes the
resonance after roughly one cycle of libration. In contrast, from the
viewpoint of the $(2,0,2)$ resonance (right panel), the orbit is at first
outside the resonance but then becomes trapped soon after it left the OLR at
around $9\,{\rm Gyr}$.

The important question here is how much does the overlap of resonances affect
the calculation of the net torque based on the averaged Hamiltonian? The first
thing to mention is that the three sets of resonances with $N_{\psi}\in
2,0,-2$ have their main territories at $\beta=0,\pi/2,\pi$ respectively, so
the majority of resonant phase-space is well isolated. Secondly, the
irregular/chaotic orbits found at the overlapped regions still librate quasi-
periodically around one of the two resonances over a finite time interval.
Thirdly, these oscillations on top of the libration behave mostly
stochastically without any apparent law, so it should not leave a systematic
bias when averaged over the whole set of orbits. Hence we do not expect
significant loss of accuracy in estimating the torque using the averaged
Hamiltonian. We will demonstrate this below with test particle simulations.

###  4.4 Total torque on the halo

Figure 10: Total torque on orbits trapped in each resonance divided by the
total halo mass. From top to bottom, the three strongest resonances are shown
(note the different y-axis scale in each figure). The analytical curves (blue)
agree well with that of the simulation (black).

After these preliminary considerations, we can now achieve our goal to
estimate the torque between the bar and the dark halo. Using equation (30),
(36), and (39), we obtain the total torque by summation over all resonances
${\bm{N}}$:

$\displaystyle\left\langle\frac{\mathrm{d}L_{z}}{\mathrm{d}t}\right\rangle$
$\displaystyle=\sum_{{\bm{N}}}\int\mathrm{d}^{3}{\bm{J}}^{\prime\prime}\mathrm{d}^{3}{\bm{\theta}}^{\prime\prime}~{}f({\bm{\theta}}^{\prime\prime},{\bm{J}}^{\prime\prime},t)\frac{\mathrm{d}L_{z}}{\mathrm{d}t}$
$\displaystyle=8(2\pi)^{2}\sum_{{\bm{N}}}N_{\varphi}\int\mathrm{d}^{2}{\bm{J}}_{\rm
f}\Psi J_{0}$
$\displaystyle\times\left[\int_{0}^{1}\\!\mathrm{d}\varepsilon\varepsilon^{2}\hat{K}\\!\int_{0}^{2\pi}\\!\\!\mathrm{d}\theta_{\ell}f(\theta_{\ell}-\Omega_{\ell}t,{\bm{J}}^{\prime\prime},0)~{}{\rm
sn}\left(\left.\hat{K}\theta_{\ell}\right|\varepsilon\right){\rm
dn}\left(\left.\hat{K}\theta_{\ell}\right|\varepsilon\right)\right.$
$\displaystyle\left.\pm\int_{1}^{\varepsilon_{\rm
cut}}\\!\mathrm{d}\varepsilon\tilde{K}\\!\int_{0}^{2\pi}\\!\\!\mathrm{d}\theta_{\rm
c}f(\theta_{\rm c}-\Omega_{\rm c}t,{\bm{J}}^{\prime\prime},0)~{}{\rm
sn}\\!\left(\\!\left.\tilde{K}\theta_{\rm
c}\right|\varepsilon^{-1}\\!\right){\rm
cn}\\!\left(\\!\left.\tilde{K}\theta_{\rm
c}\right|\varepsilon^{-1}\\!\right)\right],$ (41)

where $f$ is the distribution function (39) dependent on the
libration/circulation angle $\theta_{\ell/{\rm c}}$ and the resonant actions
${\bm{J}}^{\prime\prime}=(J_{{\rm f}_{1}},J_{{\rm f}_{2}},J_{\ell/{\rm c}})$,
$\Psi$ is the bar perturbation (55), $J_{0}$ measures the width of the
resonance (25), and $\varepsilon$ parametrizes the value of the Hamiltonian
$\bar{H}$ (25). sn, cn and dn are the Jacobi’s elliptic functions (65-67) and
we have defined $\hat{K}\equiv 2K(\varepsilon)/\pi$ and $\tilde{K}\equiv
K(\varepsilon^{-1})/\pi$ where $K$ is the complete elliptic integral of the
first kind (61). The first term in the square bracket describes the torque
applied on the trapped phase-space while the second term describes the torque
on the untrapped phase-space which we cut at $\varepsilon_{\rm cut}(>1)$ to
avoid duplicately integrating the phase space between neighbouring resonances.
However duplication cannot be fully avoided where resonances overlap. As
described above and demonstrated below, we expect the effects of overlapping
to be minor and defer detailed modelling for the overlap regions to future
study. Here, we cut the integration universally at $\varepsilon_{\rm cut}=2-3$
which is roughly the ratio of the width and interval between the neighbouring
resonances (Fig. 7). The $\pm$ sign in (41) denotes integration over the upper
and the lower circulating regions for which the mapping from $(\theta_{\rm
c},{\bm{J}}^{\prime\prime})$ to ${\bm{J}}=(J_{r},L,L_{z})$ is different.

Figure 10 compares the net torque on trapped orbits calculated with the
analytical approach (blue) to a test-particle simulation (black). From top to
bottom, we show the torques on the three strongest resonances:
$(N_{r},N_{\psi},N_{\varphi})=(0,2,2)$, $(1,0,2)$, and $(1,2,2)$. We select
the trapped orbits in the test-particle simulation by demanding that they
never pass $\theta_{\rm s}-\theta_{\rm s,res}=\pm\pi$. We then place the
particles back to their initial coordinates, re-run the simulation exclusively
for these particles, and sum up their torque. For all three resonances, the
analytical torque (blue) captures the principal features of the simulation
(black), i.e. the torque fluctuates and damps as phase mixing progresses
inside each resonance.

Figure 11: Top panel: Distribution of trapped orbits in libration period.
Bottom panel: Distribution weighted by the amplitude of the torque, which
provides the actual spectrum of the total torques shown in Fig. 10.

To better quantify the time dependence of the torque, we plot the spectra of
libration periods in the top panel of Fig. 11. The distribution of each
resonance peaks at a period which matches the main oscillating period of the
total torque shown in Fig. 10. Since the majority of orbits with long
libration periods are those trapped at small $\Psi$ (rather than those near
the separatrix) and thus have little impact on the total torque, we plot on
the lower panel of Fig. 11 the distribution when weighted by the amplitude of
the torque on each orbit. Note though that this is not a Fourier analysis of
Fig. 10 since the torque on individual orbits becomes increasingly non-
sinusoidal when approaching the separatrix (Fig. 3). The weighted distribution
clarifies that the dominant contribution to the torque comes from the
corotation resonance and in particular from those with libration period $\sim
2\,{\rm Gyr}$.

Figure 12: The libration period of the CR at the peak of its distribution
$F(T_{\ell})$ vs. bar pattern speed $\Omega_{\rm p}$ and bar strength $A$.

The libration period depends on both the bar pattern speed and the bar
strength which vary over time. A simple scaling suggests
$T_{\ell}\sim\omega_{0}^{-1}=|G\Psi|^{-1/2}\sim\Omega_{\rm p}^{-1}A^{-1/2}$.
Figure 12 shows the mode of the distribution of libration period for a given
pattern speed $\Omega_{\rm p}$ and bar strength $A$. As expected, the
libration period is shorter for higher pattern speed and for stronger bars.
Since the bar directly after formation is expected to be faster and weaker
than at present, the initial torque (and thus the initial pattern speed) in an
expected parameter range typical for the Milky Way should have fluctuated with
roughly a Gyr period.

Figure 13 compares the contribution to the total torque from trapped and
untrapped phase space with three cutoff $\varepsilon_{\rm cut}=\\{2,2.5,3\\}$.
As before, we summed over the resonances with $N_{r}\in[-10,10]$,
$N_{\psi}\in-2,0,2$, and $N_{\varphi}=2$. The torques on trapped and untrapped
orbits rise and decay on different timescales which simply reflects the
difference in their orbital periods, i.e. the period of circulation is
generally shorter than that of libration. As we increase the upper limit of
integration $\varepsilon_{\rm cut}$, the total torque on untrapped orbits gets
larger and slightly shifts towards early time because contributions from
orbits with short circulation periods far away from the resonances are added
into the calculation.

Figure 13: Comparison between the total torque exerted on resonantly trapped
(black) and untrapped (blue) orbits in the dark halo. We cut the integration
over the untrapped phase space at $\varepsilon_{\rm cut}$ to avoid duplication
by neighbouring resonances. Dynamical friction by trapped and untrapped orbits
rises and decays over different timescales. Figure 14: Total torque on the
dark halo. The theoretical curve (blue, equation 41) quantitatively explains
the behaviour of the simulation (black).

Figure 14 shows the sum of trapped and untrapped components (blue) which
quantitatively reproduces the total torque of the simulation (black). The
result assures us that the trapped orbits comprise the majority of the torque
at late times. In the very early times $<0.1\,{\rm Gyr}$, the torque of the
simulation rises sharply as high as $1.6{\,{\rm kpc}}^{2}{\,{\rm Gyr}}^{-2}$
(see Fig. 1 for the full range) but our model fails to predict this rapid
response for several reasons: (i) we have neglected the fast non-resonant
terms of the perturbation by averaging over the fast angles, (ii) the second
order Taylor approximation makes our model inaccurate in the region of rapid
circulation far from the resonances, (iii) there is a substantial volume of
untrapped phase-space that were not integrated due to the cutoff at
$\varepsilon_{\rm cut}$.

##  5 Slowdown of the bar

The natural progression of this work will be a self-consistent model for the
slowdown/growth of the bar. However, the bar’s moment of inertia is not a
constant, the halo’s self-gravity will amplify the torque, and the inflow of
gas towards the nuclear disc provides an additional angular momentum term.
Resolving these issues is beyond the scope of this paper, and so must be
deferred to future analysis. Instead, here we will discuss the effects of bar
slowdown on the angular momentum balance using the bar’s slowing rate directly
measured from local stellar kinematics (Chiba et al., 2021) which agrees with
N-body+smooth particle hydrodynamics simulation of the bar in a live dark halo
(Aumer & Schönrich, 2015).

Figure 15: Time evolution of the phase-space near the CR when the bar is
slowing down. The left column shows the phase-space density in log scale while
the right column is coloured by the initial $J_{\rm s}$. The phase-space
spiral continues to wind up as the resonance migrates upwards (note the
increasing $y$ axis from the top panel to the bottom). The bar slows down at
rate $\eta=-\dot{\Omega}_{\rm p}/\Omega_{\rm p}^{2}=0.004$ starting from
$\Omega_{\rm p}=81.8\,{\rm Gyr}^{-1}$ and ends up being $\Omega_{\rm
p}=31.0\,{\rm Gyr}^{-1}$ after $5\,{\rm Gyr}$.

Figure 15 shows the phase space perturbed by a bar that slows with
$\Omega_{\rm p}\propto 1/t$ from an initial value $\Omega_{\rm p0}=81.8\,{\rm
Gyr}^{-1}$ with rate $\eta\equiv-\dot{\Omega}_{\rm p}/\Omega_{\rm
p}^{2}=0.004$555A note on wording: This slowing rate implies $s\sim\eta/A=0.2$
which is in the ‘slow’ regime $(s<1)$. However, Chiba et al. (2021) classified
this model as a ‘rapidly’ decelerating bar in contrast to a yet slower model
that evolves almost adiabatically.. The bar strength is fixed to $A=0.02$. As
the bar slows, the resonantly trapped orbits get dragged with the resonance
towards larger $J_{\rm s}$ while they keep phase mixing internally. The dot-
dashed curves show the separatrix evaluated in a ‘time-frozen’ Hamiltonian
(when fixing the pattern speed) while the solid curves mark the approximate
boundary of orbits that remain trapped in the moving resonance (Appendix G).
As described in Chiba et al. (2021), the deceleration tilts the potential of
the resonance (76), making the effective trapped region shrink and shift
towards positive $\theta_{\rm s}$, which could also be ascribed to the Euler
force arising from the deceleration of the rotating frame of reference.

###  5.1 Capture and escape from resonance

Resonance sweeping leaves a significant change in the distribution through
capture and loss. When the resonance sweeps outwards, the untrapped orbits
swept over by the upper separatrix (dot-dashed curves in Fig. 15) will
temporarily get trapped, as adiabaticity is broken near the separatrix. Orbits
that crossed the upper separatrix near $\theta_{\rm s}-\theta_{\rm res}=\pi$
may reach the adiabatic region (solid curves), where the libration action is
conserved, and thus remain trapped and get dragged by the resonance. Other
orbits will escape the resonance before completing one cycle of libration,
proceeding to the lower circulating zone. Depending on the angle $\theta_{\rm
s}$ at which these uncaptured orbits enter the resonance, they will travel a
different amount of distance in $J_{\rm s}$: orbits crossing the separatrix
near $\theta_{\rm s}-\theta_{\rm res}=-\pi$ will immediately escape the
resonance, ending up at a relatively low $J_{\rm s}$, while those that could
almost but not quite reach the adiabatic region will be dragged a large
distance in $J_{\rm s}$ before they get damped out of the resonance. Thus,
after the passage of the resonance, orbits originally at the same action
$J_{\rm s0}$ (same color in the right column of Fig. 15) will spread over a
wide range of $J_{\rm s}$. These orbits then phase mix and form phase-space
stripes in the lower circulating region.

As seen in Fig. 15, the entire trapped volume (and the adiabatic region
within) generally grows as the bar slows. Inside the adiabatic region, trapped
orbits conserve their libration action $J_{\ell}$, so the newly captured
orbits accrete on the surface of the adiabatic region. As a result, the
resonance grows like a tree ring, where the $J_{\rm s}$ of origin increases
from the core towards the surface (Chiba & Schönrich, 2021).

Previous studies formalized the capturing and escaping processes in the
adiabatic limit ($s\ll 1$) by the transition probability averaged over the
libration/circulation angles $\theta_{\ell/{\rm c}}$ (e.g. Henrard, 1982;
Sridhar & Touma, 1996; Collett et al., 1997). However, the angle-averaged
probability has little use in modelling dynamical friction since it assumes
uniform distribution in $\theta_{\ell/{\rm c}}$. To model the transient
behaviour of dynamical friction, a phase dependent prescriptions for the
capture/loss is required. A straightforward way would be to numerically
integrate the equations of motion of the averaged Hamiltonian with time-
dependent pattern speed (Weinberg & Katz, 2007). Alternatively, one could
return to the Eulerian approach (section 3) and fully address the non-linear
terms of the CBE. This approach can in principle model the resonant structure
with arbitrary time-dependence (e.g. bar growth), although it has the
difficulty that increasing numbers of higher order modes will be required to
model the fine structures developed by phase mixing.

###  5.2 Transfer of angular momentum by resonant dragging

In addition to the torque caused by phase mixing, the moving resonances give
rise to another type of angular momentum transfer. The conservation of
libration action ensures that the adiabatic region preserves the phase-space
density at the time of capture. Since the resonances move towards higher
angular momentum, this implies that the population on resonances is much
denser than the surrounding phase space. When resonances move outward, the
freed phase space below them will be filled with orbits jumping across the
resonances, and just like lifting a heavy body in water, the halo gains net
angular momentum which further slows the bar (i.e. positive feedback).

If we accelerate the bar instead, the resonances sweep inwards. As long as
they are still denser than the surrounding phase space, the resonance frees
angular momentum, further accelerating the bar (i.e. again, positive
feedback). However, as the resonances proceed deep into the inner halo, the
surrounding density rises and so the amount of freed angular momentum
decreases up to even becoming negative when the resonances are underdense
(i.e. negative feedback).

TW84 and Weinberg (1985) called these types of angular momentum exchange the
‘dynamical feedback’ and presented a formula to predict its amount in the slow
(adiabatic) limit although they only figured in the change in angular momentum
due to orbits that jump across the resonance and separately discussed the
effect of resonant dragging. They gave an insightful interpretation of
dynamical feedback as a change in moment of inertia of the bar: since
dynamical feedback is to first order proportional to the change in pattern
speed, ($\tau_{\rm feedback}=C\dot{\Omega}_{\rm p},$ where $C>0$), the
equation of motion for the bar becomes (Weinberg, 1985)

$\displaystyle\dot{\Omega}_{\rm p}=\frac{\tau_{\rm friction}+\tau_{\rm
feedback}}{I_{\rm b}}=\frac{\tau_{\rm friction}}{I_{\rm b}-C}$ (42)

where the term $\tau_{\rm friction}$ represents dynamical friction and $I_{\rm
b}$ is the bar’s moment of inertia. This equation implies that dynamical
feedback translates as a reduction in the moment of inertia. However we
caution that (i) the equation is only valid in the slow limit where the orbits
adiabatically respond to the change in pattern speed. The response at
resonances is always delayed by an order of the libration period, meaning that
for a rapidly fluctuating torque caused by e.g. the bar-spiral interactions,
the bar will react with its bared moment of inertia. (ii) $C$ is not a
constant and depends both on the pattern speed $\Omega_{\rm p}$ as well as its
derivative $\dot{\Omega}_{\rm p}$. In particular, $C$ is a decreasing function
of $\dot{\Omega}_{\rm p}$ (i.e. the faster the bar slows, the smaller the
adiabatic/dragged region) which implies that dynamical feedback stabilizes the
bar’s slowing rate.

##  6 Conclusions

This work is a first attempt to push beyond the standard theory of dynamical
friction that quantifies the transfer of angular momentum between galactic
bars (or in fact any galactic substructures) and dark matter haloes. In
particular, we have:

  * •

shown how the standard linear perturbation theory of dynamical friction breaks
down at dynamical times for bars slowing at moderate rate (the slow regime) as
measured in the Milky Way,

  * •

developed an improved description for angular momentum transfer in this slow
regime by modelling the resonant dynamics using secular perturbation theory,

  * •

shown that dynamical friction oscillates with the libration period of the main
resonances, which would drive a significant fluctuation in bar pattern speed.

The standard treatment (; ) rests on the linearized collisionless Boltzmann
equation: the approach assumes that the bar perturbation interacts exclusively
with the pristine distribution of the halo. This assumption is valid while the
perturbed distribution remains small which, in the slow regime, only holds for
a few bar rotation periods. We have shown that the standard formalism
reasonably predicts the angular momentum transfer directly after formation of
the bar $(\lesssim 500\,{\rm Myr})$, but fails to predict the subsequent decay
of the torque driven by the non-linear response at resonances.

When the bar forms, orbits near resonances become resonantly trapped. Theses
trapped orbits periodically exchange $z$-angular momentum with the bar as they
slowly librate around the resonance. For instance, orbits trapped at the
corotation resonance around the stable Lagrange points have increasing angular
momentum while on the trailing side of the bar and are losing angular momentum
while on the leading side.

The net transfer of angular momentum to the trapped orbits is determined by
their libration angle distribution: for a typical halo in equilibrium, the
freshly trapped orbits are predominantly at lower angular momentum (libration
angle near $\theta_{\ell}\sim\pi$). This phase imbalance initially leads to a
net negative torque on the bar as this overdensity collectively librates up to
higher angular momentum. However, after half a libration period, the
overdensity starts librating back to lower angular momentum and can transfer
angular momentum back to the galactic bar. Since libration periods increase
towards the separatrix, this overdensity gradually winds up into a phase-space
spiral inside the resonance and the net torque follows secular damped
oscillations. The untrapped orbits outside the separatrix similarly exchange
angular momentum with the bar as they phase mix in their circulation angle,
although they have a much shorter mixing time.

We have shown that due to the near-harmonic effective potential around the
resonance centre, libration periods are quite similar across the trapped
volume. This provides enough coherence to cause several oscillations of the
torque in the libration period. A straightforward prediction from this is that
after bar formation, the pattern speed of the bar fluctuates with the typical
libration frequencies of the main resonances. And indeed, while we are not
aware of this having been noticed in previous literature, a closer inspection
of previous N-body simulations show these predicted bar pattern speed
oscillations (e.g. Aumer & Schönrich, 2015, Fig.2). We also note that these
oscillations in bar pattern speed take place on an order of magnitude longer
timescale $(>1\,{\rm Gyr})$ than the short-term oscillations caused by
alignments between bar and spiral patterns (Wu et al., 2016; Hilmi et al.,
2020).

The secular behaviour of orbits near resonances is modelled using the standard
method of averaging. If the bar rotates constantly, the averaged Hamiltonian
is integrable and thus admits an angle-action coordinate for each resonance.
With these new coordinates, the long-term friction can be modelled by an
analytical distribution function that simply evolves in the
libration/circulation angles at constant rates. The averaged Hamiltonian is,
however, not accurate when neighbouring resonances lie too close: overlap of
resonances leads to the onset of chaos where orbits move from one resonance to
another in a stochastic manner. In our current bar model, the limited regions
of chaotic phase-space did not significantly affect the estimation of the
total torque, although this may not be the case with a more realistic bar with
additional higher order modes (e.g. octopole).

This paper has shown qualitatively how dynamical friction may work in the slow
regime. To obtain quantitative prediction, a self-consistent model is
required. In particular we must model: (i) The growth of the bar which changes
its effective moment of inertia. (ii) The self-gravitational perturbation of
the halo which can amplify the torque on the bar by several factors (Weinberg,
1989; Dootson & Magorrian, 2021). (iii) The angular momentum supply from the
gas which sheds angular momentum to the bar as they sink into the inner bar-
region (e.g. Regan & Teuben, 2004; Athanassoula et al., 2013).

Finally, the phase-space spiral predicted in this paper is also expected to
form inside the resonances of the stellar disc. If we can identify this spiral
pattern in the observational data, it will allow us to directly constrain the
age of the Galactic bar.

## Acknowledgements

It is our pleasure to thank W. Dehnen for providing us with detailed comments
on our draft. We also thank J. Binney, J. Magorrian, E. Athanassoula, and M.
Semczuk for many helpful comments. We in particular thank M. Weinberg for a
fruitful and enlightening discussion. R.C. acknowledges financial support from
the Takenaka Scholarship Foundation and the Royal Society grant
RGF$\backslash$R1$\backslash$180095\. R.S. is supported by a Royal Society
University Research Fellowship. This work was performed using the Cambridge
Service for Data Driven Discovery (CSD3), part of which is operated by the
University of Cambridge Research Computing on behalf of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS
capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC
operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure.

## Data availability

The codes used to produce the results are available from the corresponding
author upon request.

## References

  * Arnol’d (1963) Arnol’d V. I., 1963, Russian Mathematical Surveys, 18, 9
  * Athanassoula (2003) Athanassoula E., 2003, MNRAS, 341, 1179
  * Athanassoula (2007) Athanassoula E., 2007, MNRAS, 377, 1569
  * Athanassoula et al. (2013) Athanassoula E., Machado R. E. G., Rodionov S. A., 2013, MNRAS, 429, 1949
  * Aumer & Schönrich (2015) Aumer M., Schönrich R., 2015, MNRAS, 454, 3166
  * Banik & van den Bosch (2021) Banik U., van den Bosch F. C., 2021, arXiv e-prints, p. arXiv:2103.05004
  * Binney (2016) Binney J., 2016, MNRAS, 462, 2792
  * Binney (2018) Binney J., 2018, MNRAS, 474, 2706
  * Binney (2020) Binney J., 2020, MNRAS, 495, 895
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Brizard (2013) Brizard A. J., 2013, Commun. Nonlinear Sci. Numer. Simul., 18, 511
  * Ceverino & Klypin (2007) Ceverino D., Klypin A., 2007, MNRAS, 379, 1155
  * Chandrasekhar (1943) Chandrasekhar S., 1943, ApJ, 97, 255
  * Chiba & Schönrich (2021) Chiba R., Schönrich R., 2021, MNRAS, 505, 2412–2426
  * Chiba et al. (2021) Chiba R., Friske J. K. S., Schönrich R., 2021, MNRAS, 500, 4710
  * Chirikov (1979) Chirikov B. V., 1979, Phys. Rep., 52, 263
  * Clarke & Gerhard (2021) Clarke J., Gerhard O., 2021, arXiv e-prints, p. arXiv:2107.10875
  * Collett et al. (1997) Collett J. L., Dutta S. N., Evans N. W., 1997, MNRAS, 285, 49
  * Collier & Madigan (2021) Collier A., Madigan A.-M., 2021, ApJ, 915, 23
  * Debattista & Sellwood (2000) Debattista V. P., Sellwood J. A., 2000, ApJ, 543, 704
  * Dehnen (2000) Dehnen W., 2000, AJ, 119, 800
  * Dootson & Magorrian (2021) Dootson D., Magorrian J., 2021, to be submitted
  * Halle et al. (2018) Halle A., Di Matteo P., Haywood M., Combes F., 2018, A&A, 616, A86
  * Henrard (1982) Henrard J., 1982, Celestial Mechanics, 27, 3
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hernquist & Weinberg (1992) Hernquist L., Weinberg M. D., 1992, ApJ, 400, 80
  * Hilmi et al. (2020) Hilmi T., et al., 2020, MNRAS, 497, 933
  * Hunt et al. (2019) Hunt J. A. S., Bub M. W., Bovy J., Mackereth J. T., Trick W. H., Kawata D., 2019, MNRAS, 490, 1026
  * Kaasalainen (1994) Kaasalainen M., 1994, MNRAS, 268, 1041
  * Katz et al. (2019) Katz D., et al., 2019, A&A, 622, A205
  * Lichtenberg & Lieberman (1992) Lichtenberg A., Lieberman M., 1992, Regular and Chaotic Dynamics. Springer-Verlag
  * Lynden-Bell & Kalnajs (1972) Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1
  * Malhotra (1998) Malhotra R., 1998, in Lazzaro D., Vieira Martins R., Ferraz-Mello S., Fernandez J., eds, Astronomical Society of the Pacific Conference Series Vol. 149, Solar System Formation and Evolution. p. 37
  * Monari et al. (2017) Monari G., Famaey B., Fouvry J.-B., Binney J., 2017, MNRAS, 471, 4314
  * Monari et al. (2019) Monari G., Famaey B., Siebert A., Wegg C., Gerhard O., 2019, A&A, 626, A41
  * Petersen et al. (2016) Petersen M. S., Weinberg M. D., Katz N., 2016, MNRAS, 463, 1952
  * Regan & Teuben (2004) Regan M. W., Teuben P. J., 2004, ApJ, 600, 595
  * Sellwood & Carlberg (2014) Sellwood J. A., Carlberg R. G., 2014, ApJ, 785, 137
  * Sellwood & Carlberg (2019) Sellwood J. A., Carlberg R. G., 2019, MNRAS, 489, 116
  * Sormani et al. (2015) Sormani M. C., Binney J., Magorrian J., 2015, MNRAS, 454, 1818
  * Sridhar (2019) Sridhar S., 2019, ApJ, 884, 3
  * Sridhar & Touma (1996) Sridhar S., Touma J., 1996, MNRAS, 279, 1263
  * Tremaine & Weinberg (1984) Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729
  * Trick et al. (2021) Trick W. H., Fragkoudi F., Hunt J. A. S., Mackereth J. T., White S. D. M., 2021, MNRAS, 500, 2645
  * Weinberg (1985) Weinberg M. D., 1985, MNRAS, 213, 451
  * Weinberg (1989) Weinberg M. D., 1989, MNRAS, 239, 549
  * Weinberg (1994) Weinberg M. D., 1994, ApJ, 420, 597
  * Weinberg (2004) Weinberg M. D., 2004, arXiv e-prints, pp astro–ph/0404169
  * Weinberg & Katz (2007) Weinberg M. D., Katz N., 2007, MNRAS, 375, 425
  * Wigner (1959) Wigner E. P., 1959, Group Theory and Its Application to the Quantum Mechanics of Atomic Spectra. Academic Press
  * Wisdom (2004) Wisdom J., 2004, AJ, 128, 484
  * Wu et al. (2016) Wu Y.-T., Pfenniger D., Taam R. E., 2016, ApJ, 830, 111

##  Appendix A Distribution function of Hernquist halo

The distribution function of an isotropic Hernquist halo normalized to 1 is
(Hernquist, 1990)

$\displaystyle f(E)$ $\displaystyle=\frac{1}{8\sqrt{2}\pi^{3}r_{\rm
s}^{3}v_{g}^{3}}$
$\displaystyle\times\frac{3\sin^{-1}\sqrt{\mathcal{E}}+\sqrt{\mathcal{E}(1-\mathcal{E})}(1-2\mathcal{E})(8\mathcal{E}^{2}-8\mathcal{E}-3)}{(1-\mathcal{E})^{\frac{5}{2}}}$
(43)

where $v_{g}=\sqrt{GM/r_{\rm s}}$, and $\mathcal{E}=-E/v_{g}^{2}$ is the
dimensionless binding energy. We draw initial condition for the test particles
by first sampling the energy from the differential energy distribution (Binney
& Tremaine, 2008)

$\displaystyle\mathcal{N}(E)$ $\displaystyle=f(E)g(E)$ (44)

where $g(E)$ is the density of states (i.e. phase-space volume per unit
energy)

$\displaystyle g(E)$ $\displaystyle=\frac{2\sqrt{2}\pi^{2}r_{\rm
s}^{3}v_{g}}{3\mathcal{E}^{\frac{5}{2}}}\bigl{[}3(8\mathcal{E}^{2}-4\mathcal{E}+1)\cos^{-1}\sqrt{\mathcal{E}}$
$\displaystyle\hskip
85.35826pt-\sqrt{\mathcal{E}(1-\mathcal{E})}(4\mathcal{E}-1)(2\mathcal{E}+3)\bigr{]}.$
(45)

We then choose the initial radius $r$ at fixed $E$ from the distribution

$\displaystyle P(r)$
$\displaystyle\propto\int\mathrm{d}^{3}{\bm{x}}\mathrm{d}^{3}{\bm{v}}\delta\left[\frac{1}{2}{\bm{v}}^{2}+\Phi({\bm{x}})-E\right]\delta\left[|{\bm{x}}|-r\right]$
$\displaystyle=(4\pi r)^{2}\sqrt{2\left[E-\Phi(r)\right]}.$ (46)

With $E$ and $r$, the initial speed $v=\sqrt{2\left[E-\Phi(r)\right]}$ is
determined. Finally, we pick angles randomly over a sphere in configuration
space (i.e. spherical) and in velocity space (i.e. isotropic). We confirmed
that the density and anisotropy distribution of test particles constructed in
this manner are unaltered after $10\,{\rm Gyr}$ of iteration without
perturbation by the bar.

##  Appendix B Fourier coefficients of bar potential

We follow the standard method of TW84, first expanding the bar’s potential in
spherical harmonics

$\displaystyle\Phi_{\rm
b}(r,\vartheta,\varphi,t)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\Phi_{lm}(r)Y_{lm}\left(\vartheta,\varphi-\Omega_{\rm
p}t\right).$ (47)

For our bar model (3), only terms with $l=2$ and $m=\pm 2$ are non-zero and
their coefficients are

$\displaystyle\Phi_{2,2}(r)=\Phi_{2,-2}(r)=\frac{\Phi_{\rm
b}(r)}{2Y_{2,2}(\pi/2,0)},$ (48)

where $\Phi_{\rm b}(r)$ is the radial profile of the bar potential (4). The
Fourier coefficients $\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)$ of the bar
potential expanded in the angles
${\bm{\theta}}=(\theta_{r},\theta_{\psi},\theta_{\varphi})$ can then be
expressed as (TW84)

$\displaystyle\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)=\int\frac{\mathrm{d}^{3}{\bm{\theta}}}{(2\pi)^{3}}\Phi_{\rm
b}(r,\vartheta,\varphi,t)\mathrm{e}^{-i{\mathbf{n}}\cdot{\bm{\theta}}}$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}i^{m-{\rm
n}_{\psi}}\delta_{m}^{{\rm n}_{\varphi}}Y_{l}^{{\rm
n}_{\psi}}\left(\frac{\pi}{2},0\right)d_{{\rm
n}_{\psi}m}^{l}(\beta)W_{lm}^{{\rm n}_{r}{\rm
n}_{\psi}}({\bm{J}})\mathrm{e}^{-im\Omega_{\rm p}t},$ (49)

where $d_{{\rm n}_{\psi}m}^{l}(\beta)$ is the Wigner’s small $d$-matrix (e.g.
Wigner, 1959)

$\displaystyle d_{{\rm
n}_{\psi}m}^{l}(\beta)=\sum_{t}(-1)^{t}\frac{\sqrt{(l+{\rm n}_{\psi})!(l-{\rm
n}_{\psi})!(l+m)!(l-m)!}}{(l-m-t)!(l+{\rm n}_{\psi}-t)!t!(t+m-{\rm
n}_{\psi})!}$ $\displaystyle\hskip
71.13188pt\times\left(\cos\frac{\beta}{2}\right)^{2l+{\rm
n}_{\psi}-m-2t}\left(\sin\frac{\beta}{2}\right)^{2t+m-{\rm n}_{\psi}}.$ (50)

and

$\displaystyle W_{lm}^{{\rm n}_{r}{\rm
n}_{\psi}}({\bm{J}})\equiv\frac{1}{\pi}\int_{0}^{\pi}\mathrm{d}\theta_{r}\Phi_{lm}(r)\cos\left[{\rm
n}_{r}\theta_{r}+{\rm n}_{\psi}(\theta_{\psi}-\psi)\right].$ (51)

Note that $\theta_{\psi}-\psi$ (the azimuthal deviation from the guiding
centre) only depends on $\theta_{r}$ for a given ${\bm{J}}$, so the integral
leaves no dependence on the angles. We may exert the $\delta_{m}^{{\rm
n}_{\varphi}}$ in (49) and express the coefficients in the form of
$\hat{\Phi}_{{\mathbf{n}}}({\bm{J}},t)=\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\mathrm{e}^{-i{\rm
n}_{\varphi}\Omega_{\rm p}t}$ where

$\displaystyle\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\equiv$
$\displaystyle\sum_{l=0}^{\infty}i^{{\rm n}_{\varphi}-{\rm
n}_{\psi}}Y_{l}^{{\rm n}_{\psi}}\left(\frac{\pi}{2},0\right)d_{{\rm
n}_{\psi}{\rm n}_{\varphi}}^{l}(\beta)W_{l{\rm n}_{\varphi}}^{{\rm n}_{r}{\rm
n}_{\psi}}({\bm{J}}).$ (52)

The Fourier coefficients of the bar potential expanded in the slow-fast angles
${\bm{\theta}}^{\prime}=(\theta_{{\rm f}_{1}},\theta_{{\rm f}_{2}},\theta_{\rm
s})$ at resonance ${\bm{N}}$ are

$\displaystyle\hat{\Psi}_{{\bm{k}}}({\bm{J}}^{\prime})$
$\displaystyle=\int\frac{\mathrm{d}^{3}{\bm{\theta}}^{\prime}}{(2\pi)^{3}}\Phi_{\rm
b}(r,\vartheta,\varphi,t)\mathrm{e}^{-i{\bm{k}}\cdot{\bm{\theta}}^{\prime}}$
$\displaystyle=\int\frac{\mathrm{d}^{3}{\bm{\theta}}^{\prime}}{(2\pi)^{3}}\left[\sum_{{\mathbf{n}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\mathrm{e}^{i\left({\mathbf{n}}\cdot{\bm{\theta}}-{\rm
n}_{\varphi}\Omega_{\rm
p}t\right)}\right]\mathrm{e}^{-i{\bm{k}}\cdot{\bm{\theta}}^{\prime}}$
$\displaystyle=\sum_{{\mathbf{n}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\int\frac{\mathrm{d}^{3}{\bm{\theta}}^{\prime}}{(2\pi)^{3}}$
$\displaystyle\times\mathrm{e}^{i\left[\left({\rm n}_{r}-N_{r}\frac{{\rm
n}_{\varphi}}{N_{\varphi}}-k_{{\rm f}_{1}}\right)\theta_{{\rm
f}_{1}}+\left({\rm n}_{\psi}-N_{\psi}\frac{{\rm
n}_{\varphi}}{N_{\varphi}}-k_{{\rm f}_{2}}\right)\theta_{{\rm
f}_{2}}+\left(\frac{{\rm n}_{\varphi}}{N_{\varphi}}-k_{\rm
s}\right)\theta_{\rm s}\right]},$ (53)

where the last line follows from the relation
$(\theta_{r},\theta_{\psi},\theta_{\varphi})=(\theta_{{\rm
f}_{1}},\theta_{{\rm f}_{2}},(\theta_{\rm s}-N_{r}\theta_{{\rm
f}_{1}}-N_{\psi}\theta_{{\rm f}_{2}})/N_{\varphi}+\Omega_{\rm p}t)$. Since
${\rm n}_{\varphi}$ is restricted to $m=\pm 2$, we have for resonances with
$N_{\varphi}=2$

$\displaystyle\hat{\Psi}_{{\bm{k}}}({\bm{J}}^{\prime})$
$\displaystyle=\sum_{{\mathbf{n}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\int\frac{\mathrm{d}^{3}{\bm{\theta}}^{\prime}}{(2\pi)^{3}}$
$\displaystyle\hskip 42.67912pt\times\mathrm{e}^{i\left[\left({\rm n}_{r}\mp
N_{r}-k_{{\rm f}_{1}}\right)\theta_{{\rm f}_{1}}+\left({\rm n}_{\psi}\mp
N_{\psi}-k_{{\rm f}_{2}}\right)\theta_{{\rm f}_{2}}+\left(\pm 1-k_{\rm
s}\right)\theta_{\rm s}\right]}$
$\displaystyle=\sum_{{\mathbf{n}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\delta_{{\rm
n}_{r},k_{{\rm f}_{1}}\pm N_{r}}\delta_{{\rm n}_{\psi},k_{{\rm f}_{2}}\pm
N_{\psi}}\delta_{k_{\rm s},\pm 1}$ $\displaystyle=\hat{\Phi}_{(k_{{\rm
f}_{1}}\pm N_{r},k_{{\rm f}_{2}}\pm N_{\psi},\pm
N_{\varphi})}({\bm{J}})\delta_{k_{\rm s},\pm 1}.$ (54)

Therefore, the resonant term $\Psi({\bm{J}}^{\prime})\equiv
2\left|\hat{\Psi}_{(0,0,1)}({\bm{J}}^{\prime})\right|$ in the averaged
Hamiltonian (22) is

$\displaystyle\Psi({\bm{J}}^{\prime})=2\left|\sum_{l=0}^{\infty}i^{N_{\varphi}-N_{\psi}}Y_{l}^{N_{\psi}}\left(\frac{\pi}{2},0\right)d_{N_{\psi}N_{\varphi}}^{l}(\beta)W_{lN_{\varphi}}^{N_{r}N_{\psi}}({\bm{J}})\right|.$
(55)

##  Appendix C Conservation of circulation action

Figure 16: Changes in circulation action $J_{\rm c}$ (blue, left $y$-axis) at
the CR as a function of $\varepsilon$ and bar growth time $t_{\rm growth}$.
The black dashed curve shows the circulation period (right $y$-axis).

Figure 16 demonstrates the conservation of circulation action $J_{\rm c}$ in
our growing bar simulation presented in Fig. 5. The blue curves show the
change in circulation action $\Delta J_{\rm c}/J_{\rm c0}=(J_{\rm c}-J_{\rm
c0})/J_{\rm c0}$ of untrapped orbits in the lower circulating region for a
given bar growth time $t_{\rm growth}$. $J_{\rm c0}$ and $J_{\rm c}$ are
obtained by numerically integrating the orbits’ $J_{\rm s}$ over $\theta_{\rm
s}$ (33) before and after bar formation. We plot the mean change in $J_{\rm
c}$ averaged over the final circulation angle $\theta_{\rm c}$. The $x$-axis
ranges from $\varepsilon=1$ (separatrix) to $1.4$ and the top $x$-axis
indicates the corresponding coordinate in $J_{\rm s}-J_{\rm s,res}$ at
$\theta_{\rm s}-\theta_{\rm s,res}=0$. The change in $J_{\rm c}$ becomes
larger towards the separatrix and with shorter bar growth time. With bar
growth time $t_{\rm growth}=4$ Gyr used in Fig. 5, most of the untrapped
orbits conserve $J_{\rm c}$ to high accuracy $(<0.1\%)$. For reference, we
also plot the circulation period calculated by orbit integration (black
dashed, right $y$-axis) which is several factors smaller than the adopted
$t_{\rm growth}$ apart from the vicinity of the separatrix.

##  Appendix D Linear response

We have shown in section 3 that standard linear perturbation theory fails to
predict the torque beyond few hundred Myrs. Here we take a closer look at this
issue by directly observing the linear response in the slow angle-action
plane. The Fourier coefficients of the linear response (13) to a constantly
rotating perturbation is

$\displaystyle\hat{f}_{{\mathbf{n}}}({\bm{J}},t)$
$\displaystyle=i{\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-i{\mathbf{n}}\cdot\bm{\Omega}(t-t^{\prime})}\left[\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\mathrm{e}^{-{\rm
n}_{\varphi}\Omega_{\rm p}t^{\prime}}\right]$
$\displaystyle=i{\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\mathrm{e}^{-i{\mathbf{n}}\cdot\bm{\Omega}t}\left[\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{i\left({\mathbf{n}}\cdot\bm{\Omega}-{\rm
n}_{\varphi}\Omega_{\rm p}\right)t^{\prime}}\right]$
$\displaystyle={\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\frac{\mathrm{e}^{-i{\rm
n}_{\varphi}\Omega_{\rm
p}t}-\mathrm{e}^{-i{\mathbf{n}}\cdot\bm{\Omega}t}}{{\mathbf{n}}\cdot\bm{\Omega}-{\rm
n}_{\varphi}\Omega_{\rm p}}.$ (56)

The total linear response is then

$\displaystyle f_{1}({\bm{\theta}},{\bm{J}},t)$
$\displaystyle=\sum_{{\mathbf{n}}}\hat{f}_{{\mathbf{n}}}({\bm{J}},t)\mathrm{e}^{i{\mathbf{n}}\cdot{\bm{\theta}}}$
$\displaystyle=\sum_{{\mathbf{n}}}{\mathbf{n}}\cdot\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\frac{\mathrm{e}^{i\left({\mathbf{n}}\cdot{\bm{\theta}}-{\rm
n}_{\varphi}\Omega_{\rm
p}t\right)}-\mathrm{e}^{i\left({\mathbf{n}}\cdot{\bm{\theta}}-{\mathbf{n}}\cdot\bm{\Omega}t\right)}}{{\mathbf{n}}\cdot\bm{\Omega}-{\rm
n}_{\varphi}\Omega_{\rm p}}.$ (57)

Note that this equation contains the resonant frequency at the denominator but
it does not diverge at the resonance because the numerator also vanishes
there. From the condition that the perturbed potential is real
$\hat{\Phi}_{-{\mathbf{n}}}({\bm{J}},t)=\hat{\Phi}^{\ast}_{{\mathbf{n}}}({\bm{J}},t)$,
one can express $f_{1}$ as

$\displaystyle f_{1}({\bm{\theta}},{\bm{J}},t)$
$\displaystyle\\!=\\!\sum_{{\mathbf{n}}}{\mathbf{n}}\\!\cdot\\!\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\left|\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\right|\frac{\cos{\left(\theta_{\rm
s}\\!-\\!\theta_{\rm s,res}\right)}\\!-\\!\cos{\left(\theta_{\rm
s}\\!-\\!\theta_{\rm s,res}\\!-\\!\Omega_{s}t\right)}}{\Omega_{s}}$
$\displaystyle\\!=\\!-\sum_{{\mathbf{n}}}{\mathbf{n}}\\!\cdot\\!\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\left|\hat{\Phi}_{{\mathbf{n}}}({\bm{J}})\right|\sin\left(\theta_{\rm
s}\\!-\\!\theta_{\rm
s,res}\\!-\\!\Omega_{s}t/2\right)\frac{\sin\left(\Omega_{s}t/2\right)}{\Omega_{s}/2},$
(58)

where $\theta_{\rm s}={\mathbf{n}}\cdot{\bm{\theta}}-{\rm
n}_{\varphi}\Omega_{\rm p}t$ and $\Omega_{s}={\mathbf{n}}\cdot\bm{\Omega}-{\rm
n}_{\varphi}\Omega_{\rm p}$. A similar equation is presented in Weinberg &
Katz (2007) (equation 15). At the resonance $\Omega_{s}\rightarrow 0$, the
term $\sin\left(\Omega_{s}t/2\right)/\left(\Omega_{s}/2\right)$ becomes $t$,
implying linear growth of perturbation $f_{1}$ in time. For linear theory to
hold, $f_{1}$ must remain sufficiently smaller than the variation of the
background distribution $\Delta f_{0}$ across the resonance $\Delta J_{\rm
s}\sim J_{0}=\sqrt{-\Psi/G}$ (25). This requires

$\displaystyle 1\gg\frac{f_{1}}{\Delta f_{0}}\sim\frac{\Psi}{\Delta J_{\rm
s}}t\sim\sqrt{-\Psi G}t=\omega_{0}t.$ (59)

Hence the timescale at which linear theory is valid in the slow limit is set
by the libration period $\omega_{0}^{-1}\gg t$.

Figure 17: Comparison between simulation (left), linear perturbation theory
(middle) and secular perturbation theory (right). As in Fig. 4, we show the
time evolution of the phase space density near the corotation resonance at
$(J_{{\rm f}_{1}},J_{{\rm f}_{2}})=(10,0)\,{\rm kpc}^{2}\,{\rm Gyr}^{-1}$.
Linear theory properly predicts the phase flow for the first few hundred Myrs
but fails to predict the subsequent behaviour.

Figure 17 compares the simulated phase flow near the corotation resonance
(left) with linear perturbation theory (middle) and secular perturbation
theory (right) described in section 4. Here, the libration period at the core
of the resonance is $T_{\ell}=2\pi/\omega_{0}\sim 1.2\,{\rm Gyr}$. For the top
five plots, the time since bar formation increases in steps of $0.25\,{\rm
Gyr}$, while the bottom plot shows the prediction after $10\,{\rm Gyr}$. The
linear response oscillates sinusoidally in the slow angle with phase
velocities that vanish towards the resonance, as described by the term
$\sin{\left(\theta_{\rm s}\\!-\\!\theta_{\rm
s,res}\\!-\\!\Omega_{s}t/2\right)}$ in (58). For the first few hundred Myrs
(top two rows), this shearing motion results in a change in distribution
similar to that caused by the motion of libration and circulation. However,
errors in linear theory become prominent beyond $t\gtrsim 0.2T_{\ell}$.

Despite the growing error in linear theory, past works (e.g. Lynden-Bell &
Kalnajs, 1972, TW84) have commonly assumed that the perturbation was switched
on in the distant past which is equivalent to taking $t\rightarrow\infty$. As
noted in Weinberg & Katz (2007), in this time-asymptotic limit, the function
$\sin\left(\Omega_{s}t/2\right)/\left(\Omega_{s}/2\right)$ in (58) approaches
a delta function $\pi\delta(\Omega_{s}/2)$666 In the sense that,
$\displaystyle\lim_{t\rightarrow\infty}\int\mathrm{d}^{3}{\bm{J}}g({\bm{J}})\frac{\sin\left(\Omega_{s}t/2\right)}{\Omega_{s}/2}=g({\bm{J}})\pi\delta(\Omega_{s}/2)$
(60) , i.e. the wavelength of the linear response in frequency (or action)
space gets narrower and narrower with time, and eventually an integration over
phase space will only leave contribution from the perfectly resonant orbits.
This is the prediction of the LBK formula (16). This time-asymptotic limit not
only ignores transient phenomena relevant in real galaxies, but also leads to
qualitatively wrong conclusions: in the limit $t\rightarrow\infty$, LBK
predicts that contribution from non-resonant orbits vanish by phase mixing and
only the perfectly resonant orbits give rise to a non-zero torque, whereas in
truth, the resonant orbits would also phase mix within the trapped zone
resulting in a zero torque as shown in Fig. 1.

##  Appendix E Elliptic functions

The complete elliptic integral of the first $K(k)$ and second $E(k)$ kinds are
defined as

$\displaystyle K(k)$
$\displaystyle=\int_{0}^{\frac{\pi}{2}}\frac{\mathrm{d}\theta}{\sqrt{1-k^{2}\sin^{2}\theta}},$
(61) $\displaystyle E(k)$
$\displaystyle=\int_{0}^{\frac{\pi}{2}}\mathrm{d}\theta\sqrt{1-k^{2}\sin^{2}\theta}.$
(62)

The incomplete elliptic integral of the first kind $F(\phi|k)$ is defined as

$\displaystyle
u=F(\phi|k)\equiv\int_{0}^{\phi}\frac{\mathrm{d}\theta}{\sqrt{1-k^{2}\sin^{2}\theta}}.$
(63)

and the Jacobi elliptic functions are defined via the inverse of $F(\phi|k)$:

$\displaystyle\phi=F^{-1}(u|k),$ (64) $\displaystyle{\rm
sn}(u|k)\equiv\sin\phi,$ (65) $\displaystyle{\rm cn}(u|k)\equiv\cos\phi,$ (66)
$\displaystyle{\rm dn}(u|k)\equiv\sqrt{1-k^{2}\sin^{2}\phi}.$ (67)

The derivatives of the Jacobi elliptic functions are

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}u}{\rm sn}(u|k)$ $\displaystyle={\rm
cn}(u|k){\rm dn}(u|k),$ (68) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}u}{\rm
cn}(u|k)$ $\displaystyle=-{\rm sn}(u|k){\rm dn}(u|k),$ (69)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}u}{\rm dn}(u|k)$
$\displaystyle=-k^{2}{\rm sn}(u|k){\rm cn}(u|k).$ (70)

##  Appendix F Overlap of resonances

Figure 18: Overlap of resonances with $N_{\psi}=2$ (black), $N_{\psi}=0$
(blue) and $N_{\psi}=-2$ (red). See text for detail.

Figure 18 overlays all three sets of resonances with $N_{\psi}=2$ (black),
$N_{\psi}=0$ (blue) and $N_{\psi}=-2$ (red). As in Fig. 7, $N_{\varphi}$ is 2
for all resonances and $N_{r}$ increase from bottom left to top right: 0 to 9
(black), 1 to 10 (blue), and 2 to 11 (red). The thick solid curves mark the
resonance, the dotted curve mark the maximum width of the trapped zone, and
the thin lines indicate the direction of libration. All three resonances
intersect near $L=0$ where $\Omega_{r}=2\Omega_{\psi}$ is satisfied. Overlap
of bar resonances is absent at orbital inclination $\beta=0$ (left) and
$\beta=\pi$ (right), implying that chaos due to resonance overlap is of less
concern in disk dynamics. At $\beta=\pi/2$ (middle), we see resonances
overlapping particularly at low $L$ and high $J_{r}$.

##  Appendix G Orbit subject to a slowing bar

We model an orbit trapped and dragged by a slowing bar using the Hamiltonian
(22) with a time-dependent pattern speed $\Omega_{\rm p}(t)$,

$\displaystyle\bar{H}(\theta_{\rm s},{\bm{J}}^{\prime},t)$
$\displaystyle=H_{0}({\bm{J}}^{\prime})-N_{\varphi}\Omega_{\rm p}(t)J_{\rm
s}+\Psi({\bm{J}}^{\prime})\cos\left(\theta_{\rm s}-\theta_{\rm s,res}\right).$
(71)

As before, we Taylor expand around $J_{\rm s}=J_{\rm s,res}$ which satisfies
the resonance condition at $t=t_{\rm res}$:

$\displaystyle\bar{H}(\theta_{\rm s},{\bm{J}}^{\prime},t)$
$\displaystyle\simeq\left[H_{0}({\bm{J}}^{\prime})-N_{\varphi}\Omega_{\rm
p}(t)J_{\rm s}\right]_{J_{\rm
s,res}}\\!+\left[\frac{{\partial}H_{0}}{{\partial}J_{\rm
s}}-N_{\varphi}\Omega_{\rm p}(t)\right]_{J_{\rm s,res}}\\!\\!\Delta$
$\displaystyle+\frac{1}{2}G({\bm{J}}_{\rm f},J_{\rm
s,res})\Delta^{2}+\Psi({\bm{J}}_{\rm f},J_{\rm s,res})\cos\theta,$ (72)

where we wrote $\Delta\equiv J_{\rm s}-J_{\rm s,res}$,
$\theta\equiv\theta_{\rm s}-\theta_{\rm s,res}$ and
$G\equiv\frac{{\partial}^{2}H_{0}}{{\partial}J_{\rm s}^{2}}$ for convenience.
Removing the first bracket term and substituting
$\frac{{\partial}H_{0}}{{\partial}J_{\rm s}}\bigr{|}_{J_{\rm
s,res}}={\bm{N}}\cdot\bm{\Omega}({\bm{J}}_{\rm f},J_{\rm
s,res})=N_{\varphi}\Omega_{\rm p}(t_{\rm res})$ in the second bracket term
yields (Weinberg, 1994)

$\displaystyle\bar{H}(\theta_{\rm s},{\bm{J}}^{\prime},t)$
$\displaystyle=-N_{\varphi}\left[\Omega_{\rm p}(t)-\Omega_{\rm p}(t_{\rm
res})\right]\Delta+\frac{1}{2}G\Delta^{2}+\Psi\cos\theta.$ (73)

From the Hamilton’s equations of motion:

$\displaystyle\dot{\theta}=-N_{\varphi}\left[\Omega_{\rm p}(t)-\Omega_{\rm
p}(t_{\rm res})\right]+G\Delta~{}~{}{\rm
and}~{}~{}\dot{\Delta}=\Psi\sin\theta,$ (74)

we obtain

$\displaystyle\ddot{\theta}+\omega_{0}^{2}\left(\sin\theta-s\right)=0,$ (75)

where $\omega_{0}^{2}\equiv-G\Psi$ and
$s\equiv-\frac{N_{\varphi}\dot{\Omega}_{\rm p}}{\omega_{0}^{2}}$ is the
dimensionless speed parameter (TW84). If we assume that $\dot{\Omega}_{\rm p}$
vary little during the typical libration period $2\pi/\omega_{0}$, $s$ can be
assumed constant and we obtain the following energy integral:

$\displaystyle E_{\rm p}=\frac{1}{2}{\dot{\theta}}^{2}+V(\theta)~{}~{}{\rm
where}~{}~{}V(\theta)=\omega_{0}^{2}\left(-\cos\theta-s\theta\right).$ (76)

Trapped orbits are confined in the potential well of $V(\theta)$ which is
linearly tilted due to the second term (see Fig.4 of Chiba et al., 2021). The
crest of the potential is at $\theta_{\rm
sep}=\sin^{-1}(s)~{}~{}(\pi/2\leq\theta_{\rm sep}\leq\pi)$ and the minimum
energy $E_{\rm p}$ required to reach the crest is

$\displaystyle E_{\rm{p,sep}}=\omega_{0}^{2}\left(-\cos\theta_{\rm
sep}-s\theta_{\rm sep}\right).$ (77)

Orbit with $E_{\rm p}=E_{\rm{p,sep}}$ approximately marks the boundary of
phase space that remains trapped in resonance. Hence, by a slight abuse of
language, we will refer to the orbit $E_{\rm p}=E_{\rm{p,sep}}$ as the
‘separatrix’. Inserting (77) and (74) to (76) yields an equation for the
separatrix:

$\displaystyle\Delta_{\pm}$ $\displaystyle=\frac{N_{\varphi}\left[\Omega_{\rm
p}(t)\\!-\\!\Omega_{\rm p}(t_{\rm
res})\right]}{G}\pm\frac{\omega_{0}}{G}\sqrt{2\left[\cos\theta\\!-\\!\cos\theta_{\rm
sep}\\!+\\!s\left(\theta\\!-\\!\theta_{\rm sep}\right)\right]}.$ (78)

The first term describes the movement of the resonance since $t_{\rm res}$ and
the second term describes the shape of the separatrix. Since the first term is
time dependent, the separatrix is not closed in phase space. However, in the
slow limit $s\ll 1$ the distance the resonance travels in a typical libration
period is smaller than the size of the resonance, so we may mark the
approximate area of trapped phase space by neglecting the first term. This is
drawn in Fig. 15. In comparison to the separatrix of the time-frozen
Hamiltonian (dot-dashed black), the separatrix of the time-dependent
Hamiltonian (solid black) is contracted and shifted towards $\theta$,
describing the trapped boundary of the simulation remarkably well. The first
term in fact breaks/opens the separatrix and thus allows orbits to enter or
escape the trapped region.

Setting $s=0$ in equation (78), we recover the equation for the separatrix of
a standard pendulum Hamiltonian (e.g. Lichtenberg & Lieberman, 1992):

$\displaystyle\Delta_{\pm}=\pm\frac{2\omega_{0}}{G}\cos\frac{\theta}{2}.$ (79)