# Dynamical friction and feedback on galactic bars in the general fast-slow
regime

Rimpei Chiba1,2,3
1Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George Street, Toronto, ON M5S 3H8, Canada
2Mullard Space Science Laboratory, University College London, Holmbury St.
Mary, Dorking, Surrey, RH5 6NT, UK
3Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Parks
Road, Oxford, OX1 3PU, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Current theories of dynamical friction on galactic bars are based either on
linear perturbation theory, which is valid only in the fast limit where the
bar changes its pattern speed rapidly, or on adiabatic theory, which is
applicable only in the slow limit where the bar’s pattern speed is near-
constant. In this paper, we study dynamical friction on galactic bars spinning
down at an arbitrary speed, seamlessly connecting the fast and slow limits. We
treat the bar-halo interaction as a restricted $N$-body problem and solve the
collisionless Boltzmann equation using the fast-angle-averaged Hamiltonian.
The phase-space distribution and density wakes predicted by our averaged model
are in excellent agreement with full 3D simulations. In the slow regime where
resonant trapping occurs, we show that, in addition to the frictional torque,
angular momentum is transferred directly due to the migration of the trapped
phase-space: trapped orbits comoving with the resonance typically gain angular
momentum, while untrapped orbits leaping over the trapped island lose angular
momentum. Due to the negative gradient in the distribution function, gainers
typically outnumber the losers, resulting in a net negative torque on the
perturber. Part of this torque due to the untrapped orbits was already
identified by Tremaine & Weinberg who named the phenomenon dynamical feedback.
Here, we derive the complete formula for dynamical feedback, accounting for
both trapped and untrapped orbits. Using our revised formula, we show that
dynamical feedback can account for up to $30\%$ of the total torque on the
Milky Way’s bar.

###### keywords:

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

††pubyear: 2023††pagerange: Dynamical friction and feedback on galactic bars
in the general fast-slow regime–B

##  1 Introduction

Dynamical friction is a process by which a heavy object moving in space
gravitationally transfers net momentum to the surrounding matter
(Chandrasekhar, 1943). It is one of the key processes that drives the secular
evolution of galaxies, most notably the spin-down of galactic bars (e.g.
Weinberg, 1985; Hernquist & Weinberg, 1992; Athanassoula, 1996) and the
orbital decay of satellite galaxies (e.g. Lin & Tremaine, 1983; White, 1983;
Weinberg, 1986). Dynamical friction is important for probing dark matter, as
it depends not only on the existence of dark matter (e.g. Tiret & Combes,
2007; Debattista & Sellwood, 2000; Ghafourian et al., 2020), but also on its
kinematic distribution (e.g. Athanassoula, 2003; Long et al., 2014; Sellwood,
2016) and even on its particle nature (e.g. Hui et al., 2017; Lancaster et
al., 2020; Wang & Easther, 2022).

The original study of dynamical friction by Chandrasekhar (1943) considered
the drag force on a point-mass perturber moving through an infinite
homogeneous background of particles. Under the assumption of homogeneity (i.e.
a flat background potential), field particles undergo a single hyperbolic
encounter with the perturber. However, in realistic stellar systems with a
finite inhomogeneous distribution, most particles are gravitationally bound to
the system, meaning that they encounter the perturber periodically. Due to
this periodic interaction, the exchange of momentum between the particle and
the perturber becomes most effective when the two are in orbital resonance.

Tremaine & Weinberg (1984, hereafter TW84) were the first to describe
dynamical friction in terms of resonances. Using perturbation theory, they
derived an analytic expression for dynamical friction in spherical systems
which states that friction arises exclusively from resonant orbits. Their
dynamical friction formula is referred to as the LBK formula, named after
Lynden-Bell & Kalnajs (1972) who derived a similar equation in the context of
disk dynamics. TW84 further showed that the quality of dynamical friction
changes depending on the speed at which the resonance sweeps the phase space.
They classified the evolution based on the dimensionless speed of the
resonance $s$, where they called $0\leq s<1$ ‘slow’ and $1\leq s$ ‘fast’.

Our current understanding on dynamical friction is limited to either of the
two extreme evolutions: the fast limit ($s\gg 1$) or the slow limit ($s\ll
1$). Fast limit refers to a perturber rapidly changing its frequency such that
the resonances sweep the orbits before any nonlinearity can develop. In this
limit, a time-dependent linear perturbation theory111For a time shorter than
the libration period, linear theory works regardless of the speed of the
resonance (e.g. Magorrian, 2021) (see also Section 4.4). Since the libration
period scales as the inverse square root of the perturbation amplitude, linear
theory holds if the perturbation is weak (e.g., Poisson shot noise, Fouvry et
al., 2015b) or short lived (e.g., transient spiral arms, Carlberg & Sellwood,
1985). gives a decent approximation (Weinberg, 2004; Banik & van den Bosch,
2021). It is worth noting that the time-dependent linear-torque formula
reduces to the classical LBK formula if one assumes that (i) an infinite time
has passed since the perturber emerged (the time-asymptotic limit) and (ii)
the perturber rotates at a fixed pattern speed (the slow limit)222Some authors
refer to the time-asymptotic limit as the adiabatic approximation and the slow
limit as the secular approximation (Banik & van den Bosch, 2021; Kaur & Stone,
2022). The linear approximation, however, breaks down precisely at the
intersection of these two limits (e.g. Chiba & Schönrich, 2022). Thus the LBK
formula is at best qualitative. The structure of the density wakes
corresponding to the LBK formula was studied in detail by Kaur & Stone (2022).

The slow limit is the opposite extreme where the perturber evolves at a near-
constant frequency. In this limit, the orbital response near resonances grows
nonlinearly and saturates due to resonant trapping. Trapping is a nonlinear
effect which cannot be treated with linear theory. In the slow limit, however,
the change in the perturbing frequency is assumed to be slow enough that the
evolution of the trapped orbits can be modelled by identifying their adiabatic
invariants by the method of averaging (e.g. Neishtadt, 1975; Henrard, 1982;
Sridhar & Touma, 1996). Chiba & Schönrich (2022) formulated dynamical friction
on galactic bars in this limit and showed that the phase mixing of the trapped
orbits gives rise to an oscillating friction that damps over several libration
periods. Banik & van den Bosch (2022) proposed a similar mechanism to explain
the oscillation and subsequent stalling of sinking satellites at the core of
galaxies (e.g. Read et al., 2006).

Unfortunately, neither the fast nor slow limit is appropriate for describing
the evolution of real galaxies. The perturber inevitably changes its frequency
in response to dynamical friction, yet the rate of change is not always fast
enough to justify the linear approximation. For example, both $N$-body
simulations (e.g. Petersen et al., 2016; Halle et al., 2018) and observations
(e.g. Monari et al., 2019; Binney, 2020; Chiba & Schönrich, 2021) indicate
that important evolutions of barred galaxies occur in the general slow regime
where substantial amount of stars and dark matter are captured into resonance.
To complicate things further, some orbits may evolve in the fast regime while
others evolve in the slow regime, since $s$ depends not only on the rate of
change of the pattern speed but also on the libration frequency of individual
orbits, which varies with resonances and with the position along each
resonance.

The only study that has addressed dynamical friction in the general fast-slow
regime is TW84. By solving the nonlinear orbital motions, they formulated the
net change in the angular momentum of orbits swept by resonances in the time-
asymptotic limit. They found that untrapped orbits swept by a slowly moving
resonance provide a torque that directly depends on the rate of change of the
pattern speed and therefore modifies the effective moment of inertia of the
perturber. TW84 termed this torque dynamical feedback to distinguish it from
dynamical friction. There are, however, two issues in TW84’s formalism: First,
because they take the time-asymptotic limit, their theory fails to capture the
transient behaviour of dynamical friction which lasts for several libration
periods and can therefore have a significant impact on the perturber’s
evolution. Secondly, they did not take into account (though discussed) the
dynamical feedback exerted by trapped orbits. In view of the conservation of
phase-space volume, the loss (gain) of angular momentum by untrapped orbits
leaping over the resonance and the gain (loss) by the trapped orbits dragged
by the resonance occur simultaneously and should be considered together.

In this paper, we address the shortcomings of TW84 in the context of bar-halo
interaction. First, we study dynamical friction and feedback using a semi-
numerical approach, i.e., we numerically solve the collisionless Boltzmann
equation using a Hamiltonian analytically averaged over the fast motions. This
semi-numerical approach, inspired by Weinberg & Katz (2007), allows us to
understand the mechanism of dynamical friction in the general fast-slow regime
from the fundamental viewpoint of nonlinear dynamics in a fully time-dependent
manner. We then derive a complete analytic expression of dynamical feedback
that incorporates contributions from both trapped and untrapped orbits. The
latter is shown to be equivalent to the nonlinear-torque formula derived by
TW84. We discuss its application to the Milky Way’s bar.

This paper is structured as follows. In Section 2, we introduce our bar-halo
model. In Section 3, we describe the fundamental phase-space dynamics around a
moving resonance using an angle-averaged Hamiltonian with an action variable
comoving with the resonance. In Section 4, we compute the time evolution of
the distribution function in the slow angle-action plane and discuss the
mechanism of dynamical friction and feedback in the general fast-slow regime.
In Section 5, we calculate the density response in real space and connect the
phase-space dynamics to the generation of density wakes. In Section 6, we
formulate dynamical feedback from both trapped and untrapped orbits using the
angle-averaged distribution function $\tilde{f}$. We then model the evolution
of $\tilde{f}$ using an extension of adiabatic theory and quantify dynamical
feedback on the Galactic bar. We summarize our results in Section 7.

##  2 Model

We restrict our analysis to the behaviour of test particles moving under the
combined potential of the bar and the dark halo, i.e., we ignore the self-
gravity of the halo’s response. Following our previous work (Chiba &
Schönrich, 2022), we adopt a spherical isotropic halo with a Hernquist (1990)
profile:

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

where $G$ is the gravitational constant, $M=1.5\times 10^{12}\,{\rm
M}_{\odot}$ is the total halo mass, and $r_{\rm s}=20\,{\rm kpc}$ is the halo
scale radius.

We perturb this halo with a quadrupole bar rotating at a time-dependent
pattern speed $\Omega_{\rm p}(t)$:

$\displaystyle\Phi_{1}(r,\vartheta,\varphi,t)=\Phi_{\rm
b}(r,t)\sin^{2}\vartheta\cos
2\left[\varphi-\int_{0}^{t}\mathrm{d}t^{\prime}\Omega_{\rm
p}(t^{\prime})\right],$ (2)

where $(r,\vartheta,\varphi)$ are the spherical coordinates. The polar angle
$\vartheta$ is measured from the pole, so the Galactic mid-plane is at
$\vartheta=\pi/2$. Since the bar typically grows as it slows (e.g. Debattista
& Sellwood, 2000; Martinez-Valpuesta et al., 2006; Villa-Vargas et al., 2009;
Athanassoula, 2014; Aumer & Schönrich, 2015; Fujii et al., 2019), we model the
radial profile of the bar’s potential such that the ratio between the bar’s
characteristic length $r_{\rm b}$ and its corotation radius $r_{\rm CR}$ is
kept constant, that is, $b\equiv r_{\rm b}/r_{\rm CR}={\rm const}$. This
results in the following form (Chiba et al., 2021):

$\displaystyle\Phi_{\rm b}(r,t)=-\dfrac{Av_{\rm
c}^{2}}{2}\left[\dfrac{r}{r_{\rm
CR}(t)}\right]^{2}\left[\dfrac{b+1}{b+r/r_{\rm CR}(t)}\right]^{5},$ (3)

where $A$ is the dimensionless strength of the bar and $v_{\rm c}=235\,{\rm
km}\,{\rm s}^{-1}$ is the local circular velocity. We adopt $A=0.02,b=0.28$ in
approximation to the hydrodynamical model of Sormani et al. (2015) which is
constrained by the kinematics of the gas in the inner Galactic disk.

Next we model the time-dependence of the bar’s pattern speed $\Omega_{\rm
p}(t)$. The simplest model one could think of is $\dot{\Omega}_{\rm p}={\rm
constant}$. However, $N$-body simulations typically find that the bar’s
slowing rate $-\dot{\Omega}_{\rm p}$ declines with decreasing $\Omega_{\rm
p}$. Thus a better model would be to assume the following dimensionless
quantity constant:

$\displaystyle\eta\equiv-\dfrac{\dot{\Omega}_{\rm p}(t)}{\Omega_{\rm
p}^{2}(t)}={\rm const}.$ (4)

Integrating the above differential equation, we have

$\displaystyle\Omega_{\rm p}(t)=\dfrac{\Omega_{\rm p0}}{1+\eta\Omega_{\rm
p0}t},$ (5)

where $\Omega_{\rm p0}$ is the initial pattern speed. A constant $\eta$ thus
implies a pattern speed evolving in inverse proportion to time333An
alternative model would be to assume $\dot{\Omega}_{\rm p}/\Omega_{\rm p}$
constant, in which case the pattern speed declines exponentially..

Figure 1: Model of the bar’s pattern speed $\Omega_{\rm p}$ (equation 5) as a
function of time for various dimensionless slowing rates $\eta$. The axes are
given in units of the bar’s initial pattern speed $\Omega_{\rm p0}$ and
rotation period $T_{\rm b0}=2\pi/\Omega_{\rm p0}$.

Figure 1 shows the pattern speed curves for a range of slowing rate $\eta$.
The time is given in units of the bar’s initial rotation period $T_{\rm
b0}=2\pi/\Omega_{\rm p0}$. For the Milky Way’s bar, its current pattern speed
is measured to be around $\Omega_{\rm p1}=33-41\,{\rm Gyr}^{-1}$ (Portail et
al., 2017; Monari et al., 2019; Sanders et al., 2019; Bovy et al., 2019; Asano
et al., 2020; Binney, 2020; Chiba & Schönrich, 2021; Kawata et al., 2021;
Clarke & Gerhard, 2022; Li et al., 2022; Leung et al., 2022; Lucey et al.,
2022; Gaia Collaboration et al., 2022) and thus $T_{\rm b1}\sim 170\pm
20\,{\rm Myr}$. Suppose that the bar was born at a pattern speed twice as fast
as its current value, so $T_{\rm b0}=T_{\rm b1}/2=85\,{\rm Myr}$. Figure 1
then shows that if $\eta=0.001$ the bar would have slowed to its present state
over $t\sim 160T_{\rm b0}\sim 13.6\,{\rm Gyr}$, while if $\eta=0.01$ it merely
took $t\sim 16T_{\rm b0}\sim 1.36\,{\rm Gyr}$. Thus, considering the age of
the Galactic bar, a reasonable range for the slowing rate would be
$\eta\in[0.001,0.01]$. Adopting a larger initial pattern speed will not change
this estimate significantly, while a smaller initial pattern speed in
principle allows for an arbitrary small $\eta$. Modelling of the local stellar
kinematics suggests $\eta=0.0036\pm 0.0011$ (Chiba et al., 2021).

Figure 2: Evolution of the bar’s pattern speed in semi-restricted $N$-body
simulations (black) fitted with equation (5) (blue).

To examine the adequacy of our model (5), and to clarify the conditions under
which fast and slow evolutions are expected, we calculate the pattern speed by
evolving the bar in accordance with the angular momentum transferred to the
halo. The method is similar to the semi-restricted $N$-body simulation of Lin
& Tremaine (1983); Weinberg (1985); Sellwood (2006), but we make the bar grow,
which requires a time-dependent model for the bar’s moment of inertia $I$.
Following Beane et al. (2022), we let $I$ vary with the bar’s corotation
radius $r_{\rm CR}$, but we choose the scaling $I\propto r_{\rm CR}^{3}$
instead of $r_{\rm CR}^{2}$ as the bar grows in mass as well as in length.
Thus

$\displaystyle\frac{I(t)}{I_{\rm rigid}}=\alpha\left[\frac{r_{\rm
CR}(t)}{r_{\rm CR1}}\right]^{3},$ (6)

where $\alpha$ is a free parameter and $r_{\rm CR1}=6.6\,{\rm kpc}$ is the
current corotation radius (Chiba & Schönrich, 2021; Clarke & Gerhard, 2022).
We have written $I$ in units of the $z$-moment of inertia of a rigid
homogeneous ellipsoid (e.g. Weinberg, 1985)

$\displaystyle I_{\rm rigid}=\frac{M_{\rm b}}{5}(a_{1}^{2}+a_{2}^{2}),$ (7)

with current bar properties: bar semi-major axis $a_{1}=5\,{\rm kpc}$, minor
axis $a_{2}=2\,{\rm kpc}$, and bar mass $M_{\rm b}=10^{10}\,{\rm M}_{\odot}$
(Wegg et al., 2015), resulting in $I_{\rm rigid}\simeq 6\times 10^{10}\,{\rm
M}_{\odot}\,{\rm kpc}^{2}$. The parameter $\alpha$ thus roughly measures the
degree of deviation from a rigid bar. We assume that the change in the bar’s
moment of inertia is due to the gain of mass from the stellar disk, not due to
dynamical friction. Hence, the term $\dot{I}\Omega_{\rm p}$ is excluded from
the bar’s equation of motion, $I\dot{\Omega}_{\rm p}=\tau$, where $\tau$ is
dynamical friction. While dynamical friction may deform the bar and change
$I$, we ignore this effect here.

The black curves in Fig. 2 show the evolution of the bar pattern speed with
$\alpha=1,4,10,16$. We choose an initial pattern speed $\Omega_{\rm
p0}=72\,{\rm Gyr}^{-1}$. For all $\alpha$, our pattern speed model (blue)
provides a reasonable fit to the simulations, where the best-fit $\eta$ is
reported on the right. In agreement with Weinberg (1985), we find that a rigid
bar ($\alpha=1$) slows significantly in a few bar rotation periods. Beane et
al. (2022) report $\alpha=4$ from their $N$-body+hydrodynamical simulation.
The evolution with $\alpha=4$ is still fast, which is in line with $N$-body
simulations without gas. A number of theories have been proposed to prevent
such a rapid slowdown: (i) the dark halo has a low central density (Debattista
& Sellwood, 2000) or is non-existent (Ghafourian et al., 2020), (ii) the halo
has a net spin (Long et al., 2014; Fujii et al., 2019), or (iii) the gas
provides angular momentum to the bar (Bournaud et al., 2005; Villa-Vargas et
al., 2010; Hopkins & Quataert, 2011; Athanassoula et al., 2013). The last
point was recently revisited by Beane et al. (2022) who showed that bars can
spin down at an arbitrary low rate within a realistic fraction of gas in the
disk ($4-16\%$). We therefore explore dynamical friction over a range of
slowing rate, $\eta\in[0.001,0.032]$, where it is to be understood that a fast
bar evolution $(\eta\sim 0.01)$ is expected in gas-poor galaxies, while a slow
bar evolution $(\eta\sim 0.001)$ requires a certain supply of angular momentum
from the gas.

##  3 Phase-space dynamics around a moving resonance

In this section, we analyse the phase-space dynamics in the neighbourhood of a
resonance moving at an arbitrary speed. The derivation up to equation (16) can
be found elsewhere (e.g. Tremaine & Weinberg, 1984; Binney, 2018), although we
repeat them here for completeness.

To simplify the dynamics, we introduce angle-action variables
$({\bm{\theta}},{\bm{J}})$, which are the natural coordinates to describe
near-integrable systems (e.g. Binney & Tremaine, 2008). By construction, the
Hamiltonian of an integrable orbit is a function of the actions alone
$H_{0}=H_{0}({\bm{J}})$, and thus the equations of motion take the simple
form: $\dot{{\bm{J}}}=-{\partial}H_{0}/{\partial}{\bm{\theta}}={\bm{0}}$ and
$\dot{{\bm{\theta}}}={\partial}H_{0}/{\partial}{\bm{J}}\equiv\mathbf{\Omega}({\bm{J}})={\rm
constant}$. Actions are not only integrals of motion, but have the important
property of being adiabatic invariants, making them especially valuable in
describing the long-term evolution of galaxies.

In our model, the spherical symmetry and time symmetry of the unperturbed
potential (equation 1) guarantee the existence of a single set of angle-action
coordinates that covers the entire phase space444The Hamiltonian of a general
axisymmetric potential is not guaranteed to possess a global system of angle-
action coordinates. However, an integrable Hamiltonian that closely
approximates the true Hamiltonian can be constructed with the torus mapping
technique introduced by McGill & Binney (1990). Binney (2018) applied this
technique to model the bar’s perturbation on the stellar disk in a realistic
three-dimensional galactic potential.. There are some freedom in the choice of
the actions; we choose ${\bm{J}}=(J_{r},L,L_{z})$ where $J_{r}$ is the radial
action, $L$ is the norm of the angular momentum vector, and $L_{z}$ is its $z$
component. The conjugate angles
${\bm{\theta}}=(\theta_{r},\theta_{\psi},\theta_{\varphi})$ describe the phase
of radial oscillation, the phase of azimuthal motion in the orbital plane, and
the longitude of the ascending node, which is a constant of motion.

Since the angle variables are $2\pi$ periodic, we can expand the perturbed
potential (equation 2) into a Fourier series

$\displaystyle H({\bm{\theta}},{\bm{J}},t)$
$\displaystyle=H_{0}({\bm{J}})+\sum_{\bm{n}}\Phi_{\bm{n}}({\bm{J}},t)~{}\mathrm{e}^{i\left[{\bm{n}}\cdot{\bm{\theta}}-n_{\varphi}\int\\!\mathrm{d}t^{\prime}\Omega_{\rm
p}(t^{\prime})\right]},$ (8)

where ${\bm{n}}\equiv(n_{r},n_{\psi},n_{\varphi})$ is a triple of integers,
and we have extracted the rotational time-dependence of the perturbation from
the coefficients $\Phi_{\bm{n}}$ and coupled it to the exponential.
$\Phi_{\bm{n}}$ is given in equation (B3) of Chiba & Schönrich (2022). Near a
resonance

$\displaystyle{\bm{N}}\cdot\mathbf{\Omega}({\bm{J}})-N_{\varphi}\Omega_{\rm
p}(t)=0,$ (9)

where ${\bm{N}}\equiv(N_{r},N_{\psi},N_{\varphi})$555Without loss of
generality, we assume $N_{\varphi}>0$., the equations of motion generated by
$H$ involve a term with ${\bm{n}}={\bm{N}}$ that evolves slowly, while others
(${\bm{n}}\neq{\bm{N}}$) rapidly alternate. This timescale separation suggests
that we may capture the secular dynamics near resonances by averaging $H$ over
the fast motions and retaining only the resonant term. Following the usual
procedure, we exploit this timescale disparity by making a canonical
transformation to the fast-slow 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})$, specific
to each resonance ${\bm{N}}$ (TW84),

$\displaystyle\theta_{{\rm f}_{1}}=\theta_{r},\hskip 19.91692pt\theta_{{\rm
f}_{2}}=\theta_{\psi},\hskip 19.91692pt\theta_{\rm
s}={\bm{N}}\cdot{\bm{\theta}}-N_{\varphi}\int\\!\mathrm{d}t~{}\Omega_{\rm
p}(t),$ (10) $\displaystyle J_{{\rm
f}_{1}}=J_{r}-\frac{N_{r}}{N_{\varphi}}L_{z},\hskip 14.22636ptJ_{{\rm
f}_{2}}=L-\frac{N_{\psi}}{N_{\varphi}}L_{z},\hskip 14.22636ptJ_{\rm
s}=\frac{L_{z}}{N_{\varphi}},$ (11)

using the generating function of the second kind

$\displaystyle S({\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}\int\\!\mathrm{d}t~{}\Omega_{\rm
p}(t)\right]J_{\rm s}.$ (12)

We write the new Hamiltonian as

$\displaystyle H^{\prime}({\bm{\theta}}^{\prime},{\bm{J}}^{\prime},t)$
$\displaystyle=H^{\prime}_{0}({\bm{J}}^{\prime},t)+\sum_{{\bm{k}}}\Psi_{{\bm{k}}}({\bm{J}}^{\prime},t)~{}\mathrm{e}^{i{\bm{k}}\cdot{\bm{\theta}}^{\prime}},$
(13)

where the perturbation is expanded in the new angles ${\bm{\theta}}^{\prime}$
with new indices ${\bm{k}}\equiv(k_{{\rm f}_{1}},k_{{\rm f}_{2}},k_{\rm s})$.
The unperturbed Hamiltonian acquires an explicit time dependence due to the
varying pattern speed

$\displaystyle
H^{\prime}_{0}({\bm{J}}^{\prime},t)=H_{0}({\bm{J}}^{\prime})+\frac{{\partial}S}{{\partial}t}=H_{0}({\bm{J}}^{\prime})-N_{\varphi}\Omega_{\rm
p}(t)J_{\rm s}.$ (14)

Averaging equation (13) over the fast angles ${\bm{\theta}}_{\rm
f}\equiv(\theta_{{\rm f}_{1}},\theta_{{\rm f}_{2}})$ yields

$\displaystyle\bar{H}(\theta_{\rm s},{\bm{J}}^{\prime},t)$
$\displaystyle=\frac{1}{(2\pi)^{2}}\int\mathrm{d}^{2}{\bm{\theta}}_{\rm
f}H^{\prime}({\bm{\theta}}^{\prime},{\bm{J}}^{\prime},t)$
$\displaystyle=\bar{H}_{0}({\bm{J}}^{\prime},t)+\sum_{k_{\rm s}}\Psi_{k_{\rm
s}}({\bm{J}}^{\prime},t)~{}\mathrm{e}^{ik_{\rm s}\theta_{\rm s}},$ (15)

where $\bar{H}_{0}=H^{\prime}_{0}$ and $\Psi_{k_{\rm
s}}\equiv\Psi_{(0,0,k_{\rm s})}$. For the dominant resonances with
$N_{\varphi}=2$, it can be shown that terms other than $k_{\rm s}=\pm 1$
vanish (Chiba & Schönrich, 2022). Hence, we have

$\displaystyle\bar{H}(\theta_{\rm s},J_{\rm s},t)=\bar{H}_{0}(J_{\rm
s},t)+\Psi(J_{\rm s},t)\cos\left(\theta_{\rm s}-\theta_{\rm s,res}\right),$
(16)

where $\Psi\equiv 2|\Psi_{(0,0,1)}|$ (Chiba & Schönrich, 2022, Appendix B).
For brevity, we have omitted the dependence on the constant fast actions
${\bm{J}}_{\rm f}\equiv(J_{{\rm f}_{1}},J_{{\rm f}_{2}})$, which may be
regarded as parameters rather than dynamical variables.

Our focus here is on the impact of the time dependence of $\Omega_{\rm p}(t)$
and the associated growth of the bar $\Psi(J_{\rm s},t)$. The evolution of
resonant systems in the classical adiabatic limit is described in detail by
Neishtadt (1975) and Henrard (1982). In the following, we summarize their key
results, which we further develop in the subsequent Section 3.2.

Figure 3: Level curves of the averaged Hamiltonian (bottom) and its
unperturbed component (top). Left column: Hamiltonian $\bar{H}$ (equation 16)
in the static action coordinate $J_{\rm s}$. Right column: Hamiltonian
$\widehat{H}$ (equation 25) in the comoving action coordinate $\Delta=J_{\rm
s}-J_{\rm s,res}(t)$ when $\eta=0.002$. $J_{0}$ is the characteristic width of
the resonance defined in equation (34).

###  3.1 Adiabatic evolution near a resonance

The slow evolution of mechanical systems is best described in terms of
adiabatic invariants. In the classical studies of resonant systems, adiabatic
invariants are found by observing the phase-space trajectories in the ‘time-
frozen’ Hamiltonian of form (16). Figure 3 left column shows the level curves
of $\bar{H}$ in the slow plane $(\theta_{\rm s},J_{\rm s})$ at some instance
of time. We chose the corotation resonance ${\bm{N}}=(0,2,2)$ at
${\bm{J}}_{\rm f}={\bm{0}}$666For the corotation resonance, ${\bm{J}}_{\rm
f}={\bm{0}}$ corresponds to unperturbed orbits that are circular ($J_{r}=0$)
with zero inclination ($L_{z}=L$). as an example, but the discussion is
general. The top-left figure shows the unperturbed Hamiltonian $\bar{H}_{0}$.
The trajectories are everywhere straight, and the Hamiltonian peaks exactly at
the resonance $J_{\rm s}=J_{\rm s,res}$ because the slope
${\partial}\bar{H}_{0}/{\partial}J_{\rm
s}={\bm{N}}\cdot\mathbf{\Omega}-N_{\varphi}\Omega_{\rm p}$ vanishes there. The
bottom-left figure shows the contours of the perturbed Hamiltonian $\bar{H}$.
Because $\bar{H}_{0}$ takes an extremum at $J_{\rm s,res}$, the presence of
even the smallest perturbation can change the topology of the phase-space
trajectories near the resonance: the level curves of $\bar{H}$ connect phase
space above and below the resonance, resulting in trajectories that oscillate
in $\theta_{\rm s}$. This is the regime of libration. Librating orbits are
referred to as being trapped in resonance. Untrapped orbits, in contrast,
remain on either side of the resonance and freely circulate, i.e., their slow
angle evolves monotonically (albeit not linearly in time).

The adiabatic invariant in the perturbed system is no longer $J_{\rm s}$ but
the phase-space area enclosed by the contours of the time-frozen $\bar{H}$

$\displaystyle J_{\ell}=\frac{1}{2\pi}\oint\mathrm{d}\theta_{\rm s}J_{\rm s},$
(17)

where we have divided the area by $2\pi$ to make it an action variable. The
angles $\theta_{\ell}$ conjugate to $J_{\ell}$ describe the phase of libration
(or circulation). Adiabatic theory assumes that $\bar{H}$ evolves slowly
enough that actions $J_{\ell}$ of most orbits are approximately preserved. The
adiabaticity, however, inevitably breaks down at the separatrix where the
orbital period becomes infinite. There, orbits may transit from one regime of
phase space to another and the actions $J_{\ell}$ change discontinuously. The
fate of an orbit crossing the separatrix is known to be extremely sensitive to
its angle $\theta_{\ell}$ and thus cannot be practically predicted. One can,
however, calculate the probability of transition averaged over $\theta_{\ell}$
(Neishtadt, 1975; Henrard, 1982), and this has been used to model the angle-
averaged (coarse-grained) distribution function of stellar systems near
resonances (Sridhar & Touma, 1996). We will come back to this point in Section
6 when we model dynamical feedback.

The adiabatic assumption requires the typical libration period $T_{\ell}$ to
be sufficiently shorter than the timescale at which the separatrices move by
the width of the trapped phase-space. Let $v_{\pm}$ be the speed of the
separatrices above and below the resonance. To order of magnitude, $v_{\pm}$
can be written as the sum of the rate at which the resonance moves
$\dot{J}_{\rm s,res}$ and the rate at which the width of the trapped region
$\delta J_{\rm res}$ changes:

$\displaystyle v_{\pm}\sim\dot{J}_{\rm s,res}\pm\delta\dot{J}_{\rm
res}=\left(1\pm\gamma\right)\dot{J}_{\rm s,res},$ (18)

where the dimensionless parameter

$\displaystyle\gamma\equiv\frac{\delta\dot{J}_{\rm res}}{\dot{J}_{\rm s,res}}$
(19)

describes the change in resonant width with respect to the change in resonant
position. The adiabatic condition requires

$\displaystyle T_{\ell}\ll\frac{\delta J_{\rm res}}{{\rm
max}(v_{+},v_{-})}=\frac{\delta J_{\rm
res}}{\left(1+|\gamma|\right)|\dot{J}_{\rm s,res}|}.$ (20)

As we shall see later (equations 34, 38), the libration period generally
scales as $T_{\ell}\sim 1/\sqrt{A}$, where $A$ is the order of perturbation
(equation 3), the resonance width scales as $\delta J_{\rm res}\sim\sqrt{A}$,
while the speed of the resonance as $|\dot{J}_{\rm s,res}|\sim\eta$ (see also
Chiba et al., 2021, Section 2.4). Therefore, to leading order, the adiabatic
assumption requires

$\displaystyle\left(1+|\gamma|\right)\frac{\eta}{A}\ll 1.$ (21)

For the problem at hand, $\gamma$ is typically of order $10^{-1}$ over the
majority of phase space (Appendix B), so the adiabatic condition is largely
dictated by the ratio $\eta/A$. Models of the Galactic bar suggest $A\sim
10^{-2}$ while $\eta\sim 10^{-3}$ (Section 2), so $\eta/A\sim 10^{-1}$. Hence,
in real bar-halo systems, the adiabatic criterion (21) is unlikely to be
sufficiently satisfied.

###  3.2 Adiabatic theory in the comoving frame

In this section, we show that the range of applicability of the adiabatic
theory can be extended – the adiabatic condition can be eased – by making a
simple canonical transformation.

As discussed in the previous Section, the adiabatic criterion depends on the
speed with which the separatrix moves and is determined predominantly by the
change in the position of the resonance rather than the accompanying change in
the shape of the separatrix. This suggests that we may find a phase-space
coordinate in which the separatrix appears more stationary by making a
canonical transformation to a slow action variable comoving with the sweeping
resonance

$\displaystyle\Delta\equiv J_{\rm s}-J_{\rm s,res}(t).$ (22)

The slow angle need not be changed, but for brevity, we redefine

$\displaystyle\theta\equiv\theta_{\rm s}-\theta_{\rm s,res}.$ (23)

A suitable generating function for this transformation is

$\displaystyle S(\theta_{\rm s},\Delta,t)=\left(\theta_{\rm s}-\theta_{\rm
s,res}\right)\left[\Delta+J_{\rm s,res}(t)\right].$ (24)

The new Hamiltonian is then

$\displaystyle\widehat{H}(\theta,\Delta,t)=\widehat{H}_{0}(\theta,\Delta,t)+\Psi(\Delta,t)\cos\theta,$
(25)

where the unperturbed part now has an additional term that is linear in the
slow angle $\theta$:

$\displaystyle\widehat{H}_{0}(\theta,\Delta,t)=\bar{H}_{0}+\frac{{\partial}S}{{\partial}t}=\bar{H}_{0}[\Delta+J_{\rm
s,res}(t),t]+\dot{J}_{\rm s,res}(t)\theta.$ (26)

This new unperturbed Hamiltonian $\widehat{H}_{0}$ has lost periodicity in the
slow angle because it returns the value of the old Hamiltonian $\bar{H}_{0}$
along lines of constant comoving action $\Delta$.

The right column of Fig. 3 shows the change in the level curves of the time-
frozen Hamiltonian caused by this canonical transformation. Note that the axes
of the left and right columns remain the same, although different labels are
used to highlight the change in canonical coordinates. In the new comoving
coordinates, the unperturbed trajectories (top-right figure) no longer appear
straight but move towards low $\Delta$ because the resonance is moving up. The
sign of $\dot{\theta}$ switches from negative to positive as the orbits pass
the resonance. As with the case in the static coordinate, the presence of the
perturbation (bottom-right) changes the topology near the resonance, but
crucially, the phase space of libration has diminished in size. In addition,
the libration is no longer symmetric around $\theta=0$: the angles of the
stable fixed points (the local extrema of $\widehat{H}$) have shifted
positively, while that of the unstable fixed points (the saddle points) have
shifted negatively. Unlike the trapped orbits, the phase curves of the
untrapped orbits are open: they approach the resonance from above, make a
U-turn at the resonance, and transition to the lower circulating regime.

Figure 4: Time variation of the libration action of an orbit trapped and
dragged by a moving resonance. The bar slowed from $\Omega_{\rm p}=60$ to
$30\,{\rm Gyr}^{-1}$ in approximately $7.5\,{\rm Gyr}$ (slowing rate
$\eta=0.002$) and the orbit was initially placed at $J_{\rm s}=550\,{\rm
kpc}^{2}\,{\rm Gyr}^{-1}$. The libration action calculated along the contours
of the time-frozen Hamiltonian in the comoving frame $\widehat{H}$ (blue) is
better preserved than that in the static frame $\bar{H}$ (black).

The phase-space area enclosed by the trapped orbits in the comoving frame
defines a new action of libration

$\displaystyle J_{\ell}\equiv\frac{1}{2\pi}\oint\mathrm{d}\theta\Delta,$ (27)

where the integral runs along constant time-frozen $\widehat{H}$. Figure 4
compares the variation of the action of libration defined by the time-frozen
$\bar{H}$ (static action coordinate) with that defined by the time-frozen
$\widehat{H}$ (comoving action coordinate). Specifically, we integrated a
single trapped orbit in the time-dependent Hamiltonian (here, static or
comoving makes no difference), and for each phase point $\bm{w}_{i}$ at each
time step, we calculated the libration action by integrating a new orbit
starting from $\bm{w}_{i}$ but in the time-frozen Hamiltonian (here, static or
comoving is relevant). Clearly, the libration action in the comoving frame is
better preserved. The dotted lines mark the time when the orbit enters the
region of libration (defined by the time-frozen $\bar{H}$ and $\widehat{H}$
respectively), which is slightly earlier in the static case because the
trapped region is larger (cf. Fig. 3).

The reason for the improvement in the conservation of $J_{\ell}$ can be
understood by re-evaluating the adiabatic condition. In the comoving frame,
the resonance is fixed, so the speed of the separatrix is given solely by the
change in shape of the trapped region

$\displaystyle v_{\pm}\sim\pm\delta\dot{J}_{\rm res}=\pm\gamma\dot{J}_{\rm
s,res}.$ (28)

The adiabatic condition is then

$\displaystyle T_{\ell}\ll\frac{\delta J_{\rm res}}{{\rm
max}(v_{+},v_{-})}=\frac{\delta J_{\rm res}}{|\gamma\dot{J}_{\rm s,res}|},$
(29)

which, to leading order, requires

$\displaystyle|\gamma|\frac{\eta}{A}\ll 1.$ (30)

As compared to the original criterion (21), the new criterion depends directly
on the size of the parameter $\gamma$ (equation 19). This implies that if the
shape of the separatrix in the comoving frame varies little ($\gamma\ll 1$),
the adiabaticity does not necessary require a slowly moving resonance
$(\eta/A\ll 1)$ as in the classical adiabatic assumption. In the example given
above (Fig. 4), $\gamma\sim 0.1$ so $\gamma\eta/A\sim 0.01$, which satisfies
the adiabatic criterion better than in the static coordinate by an order of
magnitude. While the degree of improvement depends on the chosen orbit, the
evolution is always more adiabatic in the comoving frame. This implies that
the instantaneous phase flow is better represented by the level curves of
$\widehat{H}$ rather than that of $\bar{H}$. We will demonstrate in Section 4
that the phase-space flow identified in the comoving action coordinate indeed
describes the numerical simulations remarkably well.

To gain further insights into the phase-space structure of $\widehat{H}$, let
us employ the standard pendulum approximation (e.g. Lichtenberg & Lieberman,
1992). That is, we Taylor expand $\widehat{H}$ in $\Delta$ up to second order
for the unperturbed term and to zeroth order for the perturbed term

$\displaystyle\widehat{H}(\theta,\Delta,t)$
$\displaystyle\simeq\left.\bar{H}_{0}\right|_{\Delta=0}+\left.\frac{{\partial}\bar{H}_{0}}{{\partial}\Delta}\right|_{\Delta=0}\Delta+\frac{1}{2}\left.\frac{{\partial}^{2}\bar{H}_{0}}{{\partial}\Delta^{2}}\right|_{\Delta=0}\Delta^{2}$
$\displaystyle+\left.\Psi\right|_{\Delta=0}\cos\theta+\dot{J}_{\rm
s,res}\theta.$ (31)

The first term can be removed as it does not affect the slow dynamics, and the
second term is zero at the resonance. Hence

$\displaystyle\widehat{H}(\theta,\Delta,t)$
$\displaystyle=\frac{1}{2}G\Delta^{2}-F\cos\theta+\dot{J}_{\rm s,res}\theta,$
(32)

where

$\displaystyle
G\equiv\left.\frac{{\partial}^{2}\bar{H}_{0}}{{\partial}\Delta^{2}}\right|_{\Delta=0}=\left.\frac{{\partial}\left({\bm{N}}\cdot\mathbf{\Omega}\right)}{{\partial}J_{\rm
s}}\right|_{J_{\rm s}=J_{\rm s,res}}~{}~{}~{}{\rm
and}~{}~{}~{}F\equiv-\left.\Psi\right|_{\Delta=0}.$ (33)

The Hamiltonian (32) is equivalent to that of a pendulum subject to an
external torque. For the main bar resonances, the quantity $G$ is negative in
much of the phase space (Appendix B). $F$ is defined to be negative such that
the Hamiltonian resembles that of a real pendulum. In the following, we thus
assume $FG>0$. Cases with $FG<0$ can be dealt similarly by shifting the angle
as appropriate. The two parameters together define two important scales for
the motion near a resonance:

$\displaystyle\omega_{0}=\sqrt{FG},~{}J_{0}=\sqrt{F/G}~{}~{}~{}\leftrightarrow~{}~{}~{}F=-\omega_{0}J_{0},~{}G=-\omega_{0}/J_{0}.$
(34)

When the resonance is static ($\dot{J}_{\rm s,res}=0$), $\omega_{0}$ describes
the libration frequency of orbits with infinitesimal libration amplitude and
thus provides the timescale of libration, while $J_{0}$ is one-half of the
maximum libration amplitude and thus characterizes the size of the trapped
region (Lichtenberg & Lieberman, 1992). Note that when $\dot{J}_{\rm
s,res}\neq 0$, both $G$ and $F$, as well as the two scale parameters, are time
dependent because they are evaluated on the moving resonance.

Using equation (34), we may rewrite the Hamiltonian (32) as

$\displaystyle\widehat{H}(\theta,\Delta,t)=\frac{1}{2}G\Delta^{2}-F\left(\cos\theta+s\theta\right),$
(35)

where

$\displaystyle s\equiv\frac{\dot{J}_{\rm
s,res}}{\omega_{0}J_{0}}=\frac{T_{\ell}}{T_{\rm d}}$ (36)

is a dimensionless parameter that describes the ratio of the characteristic
libration time $T_{\ell}=\omega_{0}^{-1}$ to the time the resonance takes to
drift by its width $T_{\rm d}=J_{0}/\dot{J}_{\rm s,res}$. This parameter $s$
is in fact the dimensionless speed parameter introduced in TW84. To
demonstrate this, we derive $\dot{J}_{\rm s,res}$ by differentiating the
resonance condition (9) with time

$\displaystyle 0$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}\left[{\bm{N}}\cdot\mathbf{\Omega}\left(J_{\rm
s,res}\right)-N_{\varphi}\Omega_{\rm p}\right],$ $\displaystyle=\dot{J}_{\rm
s,res}\left.\frac{{\partial}\left({\bm{N}}\cdot\mathbf{\Omega}\right)}{{\partial}J_{\rm
s}}\right|_{J_{\rm s}=J_{\rm s,res}}-N_{\varphi}\dot{\Omega}_{\rm p}.$ (37)

Using $G$ (equation 33), we obtain

$\displaystyle\dot{J}_{\rm s,res}=\frac{N_{\varphi}\dot{\Omega}_{\rm p}}{G}.$
(38)

Substituting this into equation (36), we have

$\displaystyle s=-\frac{N_{\varphi}\dot{\Omega}_{\rm p}}{GF},$ (39)

which is equivalent to equation (89) of TW84. This expression indicates that,
as opposed to the dimensionless slowing rate $\eta$ (equation 4), which
measures $\dot{\Omega}_{\rm p}$ with respect to the bar’s dynamical time
$\Omega_{\rm p}^{-1}$, the speed parameter $s$ measures $\dot{\Omega}_{\rm p}$
with respect to the characteristic libration time
$\omega_{0}^{-1}=(GF)^{\\-1/2}$. Hence, while $\eta$ is a global parameter,
$s$ is a local parameter that depends on the resonance ${\bm{N}}$ as well as
on the position on each resonance set by the fast actions ${\bm{J}}_{\rm f}$.

To clarify the factors featuring the Hamiltonian (35), we further rewrite it
in a dimensionless form:

$\displaystyle\widehat{H}=F\left[\frac{1}{2}\left(\frac{\Delta}{J_{0}}\right)^{2}+V(\theta)\right],$
(40)

where $V$ is the dimensionless effective potential

$\displaystyle V(\theta)=-\cos\theta-s\theta.$ (41)

The phase-space structure of $\widehat{H}$ is thus characterized solely by
$s$.

Figure 5: The effective potential of the resonance (41) and the contours of
the Hamiltonian (40) for various speed $s$ increasing from top to bottom.

Figure 5 shows the potential and the contours of the (time-frozen) Hamiltonian
(40) in the slow plane for various $s$. A similar phase portrait has been
observed in a variety of physical phenomena that involve a sweeping resonance:
the radial migration of planet in proto-planetary disk (Ogilvie & Lubow,
2006), the propagation of waves in an inhomogeneous plasma (e.g. Shklyar &
Matsumoto, 2009; Artemyev et al., 2018), and the acceleration of
particles/molecules by chirped lasers (e.g. Kalyakin, 2008; Armon & Friedland,
2016).

The top panel of Fig. 5 depicts the phase flow in the slow limit ($s=0$,
constant pattern speed) which resembles that of a standard pendulum. The phase
flow streams anticlockwise around the fixed points since $F<0$. Unlike the
original Hamiltonian (Fig. 3), this Hamiltonian is symmetric about the
resonance $\Delta=0$ since we have neglected the high-order Taylor series.

The middle two panels of Fig. 5 show the phase flow in the slow regime
($0<s<1$, resonance moving slowly). The potential gets tilted and the trapped
phase-space (grey area) accordingly shrinks, allowing untrapped orbits to
transit across the resonance. Only the trajectories of the trapped orbits are
closed. The action of libration (the enclosed phase-space area) can be
calculated by solving equation (40) for $\Delta$ and plugging it into equation
(27)

$\displaystyle J_{\ell}$
$\displaystyle=\frac{1}{2\pi}\oint\mathrm{d}\theta\sqrt{2J_{0}^{2}\left(\widehat{H}/F+\cos\theta+s\theta\right)}$
$\displaystyle=\frac{\sqrt{2}J_{0}}{\pi}\int_{\theta^{-}}^{\theta^{+}}\mathrm{d}\theta\sqrt{\cos\theta-\cos\theta^{+}+s\left(\theta-\theta^{+}\right)}$
(42)

where $\theta^{\pm}$ are the roots of the equation
$\widehat{H}+F(\cos\theta+s\theta)=0$. Figure 6 illustrates the definition of
$\theta^{\pm}$.

Figure 6: Schematic diagram of the effective potential of the resonance.

The bottom panel of Fig. 5 presents the fast regime $(s>1)$. The effective
potential no longer forms a local minima since

$\displaystyle\frac{\mathrm{d}V}{\mathrm{d}\theta}=\sin\theta-s<0~{}~{}\forall~{}~{}\theta,$
(43)

meaning that no orbits can stay trapped – the resonance sweeps past the orbits
before they can complete a full cycle of libration. As described in TW84 and
Chiba & Schönrich (2022) (Appendix D), this is the only regime where linear
perturbation theory works near resonances beyond the libration time777Note
though that other sources of time-dependent evolution in galaxies may disrupt
motions of libration and hence allow the system to remain in the linear
regime..

Figure 7: Black curve: The resonant-volume factor $C$ (equation 47) as a
function of the dimensionless speed $s$ (equation 36). Blue and red curves:
The local minimum $\theta_{\rm res}$ and maximum $\theta_{\rm sep}^{+}$ of the
resonant potential $V$ in the range $0<\theta_{\rm res}<\pi/2<\theta_{\rm
sep}^{+}<\pi$. Green curve: The root $\theta_{\rm sep}^{-}$ of the equation
$V(\theta)=V(\theta_{\rm sep}^{+})$ closest to $\theta_{\rm sep}^{+}$ (see
Fig. 6). The difference $\Delta\theta=\theta_{\rm sep}^{+}-\theta_{\rm
sep}^{-}$ is the angular width of the trapped phase-space.

None of the trajectories presented in Fig. 5 escape or get captured into the
resonance since we have frozen $\widehat{H}$ in time. For orbits to cross the
separatrix, the volume of trapping must change. The trapped volume is given by
the action of libration (equation 42) evaluated along the separatrix

$\displaystyle J_{\ell,{\rm sep}}=\frac{8J_{0}}{\pi}C(s),$ (44)

where the function $C(s)$ is the resonant-volume factor888The function $C(s)$
is equivalent to  $\displaystyle Q(s)=\frac{\sqrt{2}}{\pi
s}\int_{\theta^{-}_{\rm sep}}^{\theta^{+}_{\rm
sep}}\mathrm{d}\theta\sin\theta\sqrt{\cos\theta-\cos\theta^{+}_{\rm
sep}+s\left(\theta-\theta^{+}_{\rm sep}\right)}$ (45) derived in TW84
(equation 99) up to a constant factor because  $\displaystyle
4\sqrt{2}C(s)-\frac{\pi}{\sqrt{2}}Q(s)$
$\displaystyle=\frac{1}{s}\\!\int_{\theta^{-}_{\rm sep}}^{\theta^{+}_{\rm
sep}}\\!\mathrm{d}\theta\left(s\\!-\\!\sin\theta\right)\sqrt{\cos\theta\\!-\\!\cos\theta^{+}_{\rm
sep}\\!+\\!s\\!\left(\theta\\!-\\!\theta^{+}_{\rm sep}\right)}$
$\displaystyle=\frac{1}{s}\left[\frac{2}{3}\left\\{\cos\theta-\cos\theta^{+}_{\rm
sep}+s\left(\theta-\theta^{+}_{\rm
sep}\right)\right\\}^{\frac{3}{2}}\right]_{\theta^{-}_{\rm
sep}}^{\theta^{+}_{\rm sep}}$ $\displaystyle=0,$ (46) where the last line
holds since $V(\theta^{+}_{\rm sep})=V(\theta^{-}_{\rm sep})$ (see Fig. 6).

$\displaystyle C(s)=\dfrac{1}{4\sqrt{2}}\int_{\theta^{-}_{\rm
sep}}^{\theta^{+}_{\rm
sep}}\mathrm{d}\theta\sqrt{\cos\theta-\cos\theta^{+}_{\rm
sep}+s\left(\theta-\theta^{+}_{\rm sep}\right)},$ (47)

with domain $0\leq s\leq 1$ and range $0\leq C\leq 1$. The upper bound of the
integral $\theta_{\rm sep}^{+}$ is at the local maximum of the effective
potential (Fig. 6), which corresponds to the saddle point of the Hamiltonian.
Hence, $V^{\prime}(\theta_{\rm sep}^{+})=\sin\theta^{+}_{\rm sep}-s=0$, which
gives

$\displaystyle\theta_{\rm
sep}^{+}=\sin^{-1}s,~{}~{}~{}~{}\left(\frac{\pi}{2}\leq\theta_{\rm
sep}^{+}\leq\pi\right).$ (48)

The lower bound of the integral $\theta^{-}_{\rm sep}$ is given by the root of
the equation $V(\theta)=V(\theta_{\rm sep}^{+})$ closest to $\theta_{\rm
sep}^{+}$. For subsequent use, we also define $\theta_{\rm res}$ as the angle
at which the potential is locally minimum:

$\displaystyle\theta_{\rm res}=\sin^{-1}s,~{}~{}~{}~{}\left(0\leq\theta_{\rm
res}\leq\frac{\pi}{2}\right).$ (49)

Figure 7 plots the resonant-volume factor $C$ (black) as a function of the
speed parameter $s$. The function is unity at $s=0$ and monotonically declines
with increasing $s$ until it vanishes at $s=1$. We also plot the angles
$\theta_{\rm sep}^{+}$ (red) and $\theta_{\rm sep}^{-}$ (green), which
demarcate the trapped region. At $s=0$, these angles are separated by $2\pi$.
As $s$ increases, they converge and eventually coincide at $s=1$, in line with
the vanishing of $C$.

It is worth noting that as $s$ approaches unity, the bounding angles
$\theta_{\rm sep}^{\pm}$ and the centre of libration $\theta_{\rm res}$
converge to $\pi/2$, not zero. This has an interesting observational
implication: detecting a shift in the centre of libration allows us to
directly measure $s$ and thus the bar’s slowing rate $\dot{\Omega}_{\rm p}$.
For the corotation resonance, this amounts to measuring the angular deviation
of the Lagrange points from the bar’s axes, $\varphi_{\rm L}-\varphi_{\rm b}$,
which yields the speed of the resonance as $s=\sin\left[2\left(\varphi_{\rm
L}-\varphi_{\rm b}\right)\right]$. We will pursue this possibility in a
separate paper.

##  4 Evolution of the distribution function

Using the averaged Hamiltonian derived in the previous section, we now study
the evolution of the distribution function (DF) around a moving resonance. The
modelled DF can then be used to derive integrated quantities such as the
density response in real space (Section 5).

The evolution of the DF $f({\bm{\theta}},{\bm{J}},t)$ 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.$
(50)

The CBE states that, unlike in a collisional model (e.g. Hamilton et al.,
2022), the phase-space density along an orbit is conserved. Thus, we
immediately have

$\displaystyle
f({\bm{\theta}},{\bm{J}},t)=f({\bm{\theta}}_{0},{\bm{J}}_{0},0),$ (51)

and all we need is to calculate the phase-space trajectories under the
modelled Hamiltonian. Since we do not model the system self-consistently,
i.e., $H$ does not depend on $f$, determination of the trajectories can be
done independently for each point in phase space.

In the averaged system $\bar{H}$, the two fast actions ${\bm{J}}_{\rm f}$ are
conserved, so the trajectories only have four degrees of freedom. We further
assume that the DF is initially fully phase-mixed in the angles. The DF then
remains uniform in the fast angles ${\bm{\theta}}_{\rm f}$ at all times, since
the averaged trajectories are independent of ${\bm{\theta}}_{\rm f}$. The
problem then reduces to solving the trajectories in the 2D slow angle-action
space

$\displaystyle f(\theta_{\rm s},J_{\rm s},t)=f(\theta_{\rm s0},J_{\rm s0},0),$
(52)

where, as with the Hamiltonian, we have omitted references to ${\bm{J}}_{\rm
f}$.

In analysing the dynamics, it is often sufficient to observe only a fraction
of the entire phase space at a few snapshots. We take advantage of this by
using the backward-integration technique (e.g. Vauterin & Dejonghe, 1997;
Dehnen, 2000). Integrating the trajectories backward in time saves significant
numerical cost because we need only to integrate trajectories that pass the
phase space of interest at the time of taking the snapshot.

Figure 8: Time evolution of the phase-space distribution around a moving
resonance in the slow regime ($\eta=0.002,s\sim 0.16$). Left column: evolution
with the 1D averaged Hamiltonian $\bar{H}$ using the CBE (equation 52). Right
column: evolution with the 3D full Hamiltonian $H$ using a standard test-
particle simulation.

In contrast to the study in the slow limit (Chiba & Schönrich, 2022), where
analytical solutions to the equations of motion were available, here we solve
the Hamilton’s equations numerically following the method of Weinberg & Katz
(2007). This allows us to abandon the pendulum approximation, making the
averaging in ${\bm{\theta}}_{\rm f}$ the only approximation we rely on. Note
that the result of the numerical calculation does not hinge on whether we use
the static or comoving action coordinate. For the sake of computational
simplicity, we use the former. The latter will be used to demarcate the phase
space of trapping and to model dynamical feedback in Section 6. We summarize
the computational method in Appendix A.

Figure 8 shows the time-evolution of the phase-space density in the slow
angle-action plane $(\theta_{\rm s},J_{\rm s})$. As before, we choose 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}$. The bar forms at $t=0$ with
pattern speed $\Omega_{\rm p}=72\,{\rm Gyr}^{-1}$ and spins down to $36\,{\rm
Gyr}^{-1}$ in approximately $6.9\,{\rm Gyr}$ with slowing rate $\eta=0.002$.
This results in the speed parameter increasing from $s\simeq 0.134$ to
$0.184$, which falls under the slow regime ($s<1$).

The left column of Fig. 8 shows the DF calculated using the 1D averaged
Hamiltonian $\bar{H}$. The black closed curves mark the separatrix of the
resonance defined by the time-frozen Hamiltonian in the comoving frame
$\widehat{H}$ (equation 25). The resonance is initially at small $J_{\rm s}$,
where the phase-space density is high. As the bar slows, the resonance moves
towards larger $J_{\rm s}$, carrying the trapped orbits and capturing new
orbits as it grows in volume. The decrease in background density may cause an
optical illusion that the trapped region is getting denser, but in fact the
phase-space density of the trapped region is conserved, as ensured by the CBE.

The right column of Fig. 8 shows the snapshots obtained from a standard test-
particle simulation using the 3D full Hamiltonian $H$. We integrated orbits in
the galactic plane ($J_{{\rm f}_{2}}=0$) forward in time and measured the
density in the range $J_{{\rm f}_{1}}\in[0,20]\,{\rm kpc}^{2}\,{\rm
Gyr}^{-1}$. The 1D averaged model describes the full 3D simulation remarkably
well. The only structure the 1D model fails to reproduce is the groove seen at
$J_{\rm s}\approx 100\,{\rm kpc}^{2}\,{\rm Gyr}^{-1}$ which is created by the
inner 1:4 resonance. This is expected since averaging the Hamiltonian over the
fast angles ${\bm{\theta}}_{\rm f}$ effectively removes all Fourier series of
perturbations other than the target resonance.

###  4.1 Friction and feedback in the slow-nonlinear regime

Figure 8 clarifies the mechanism by which angular momentum is transferred from
the bar to the halo in the general slow regime. We describe this based on two
classifications: (i) trapped or untrapped, and (ii) dynamical friction or
dynamical feedback.

####  4.1.1 Dynamical friction by trapped orbits

In the trapped region, the motion of libration generates a winding phase-space
spiral, which gives rise to an oscillating dynamical friction on the
perturber. This was the main point of discussion in Chiba & Schönrich (2022)
and Banik & van den Bosch (2022). The phase spiral emerges because the
libration frequency drops towards the separatrix and the density along the
motion of libration is non-uniform due to the initial negative gradient in
$f$. Phase mixing acts to erase this gradient and, during this process, net
angular momentum is transferred to the halo. As mixing proceeds, the
frictional torque undergoes damped oscillations. However, because the
libration region is growing in size and capturing new orbits, fresh
inhomogeneity is constantly introduced into the trapped region, preventing
dynamical friction from vanishing entirely.

####  4.1.2 Dynamical feedback by trapped orbits

When the bar spins down due to dynamical friction, the resonance accordingly
moves, typically towards larger angular momentum (Appendix B). Since the
trapped orbits comove with the resonance, they gain angular momentum, meaning
that an additional negative torque is applied on the bar. Following TW84, we
refer to this type of torque as dynamical feedback, although the original
definition differs (see Sections 4.1.4 and 6). Dynamical feedback is
fundamentally different from friction in that it arises only when the
frequency of the perturber changes. Since dynamical feedback is proportional
to the rate at which the bar slows down, it modifies the bar’s effective
moment of inertia, typically reducing it. This allows the bar to spin down
more rapidly in response to dynamical friction. We will formulate and quantify
dynamical feedback in Section 6.

####  4.1.3 Dynamical friction by untrapped orbits

The untrapped orbits circulating on either side of the resonance also exert
dynamical friction in much the same manner as the trapped orbits: the motion
of circulation mixes orbits with different initial $J_{\rm s}$, and the
negative gradient in $f_{0}$ assures that this results in a net transfer of
angular momentum to the halo. This dynamical friction works regardless of the
value of $s$.

####  4.1.4 Dynamical feedback by untrapped orbits

In the slow regime, untrapped orbits transiting from one side of the resonance
to the other experience a net change (typically a loss) in angular momentum
that is proportional to the volume of the trapped phase-space. This torque is
caused only when the resonance shifts in position and is therefore classified
as dynamical feedback. In fact, it was this feedback torque by the untrapped
orbits that TW84 first identified and dubbed ‘dynamical feedback’. As with the
feedback from trapped orbits, this torque acts to modify the bar’s moment of
inertia, albeit in the opposite direction. The total dynamical feedback is
given by the difference between the two and is typically negative, hence
promoting the bar’s spin-down (see Section 6).

###  4.2 Wave form of the perturbation

A close inspection of Fig. 8 reveals that the phase spirals in the trapped
region have two different wave forms: (i) At the core of the resonance, the
phase spiral is continuous, as it stems from the smooth initial distribution
of orbits captured instantly at the time of bar formation. (ii) At the outer
region, which formed later, the density varies discontinuously between the
neighbouring layers. This discontinuity arises from orbital capture at the
saddle point of the Hamiltonian (cf. Figs. 5,6). At this point, two separate
phase-space trajectories approach from above and below and become adjacent
upon entering the trapped region. Since these newly trapped orbits were not
adjacent prior to capture, the phase-space density in the outer trapped region
becomes discontinuous.

Similarly, in the untrapped region, the phase stripes999It should be noted
that phase ‘spirals’ and ‘stripes’ are coordinate-dependent patterns: the
spirals in the trapped phase-space appear as stripes if plotted in, e.g., the
angle-action of libration $(\theta_{\ell},J_{\ell})$, while the stripes in the
untrapped phase-space appear as spirals if plotted in, e.g., Cartesian
canonical coordinates $(\sqrt{2J_{\rm s}}\cos\theta_{\rm s},\sqrt{2J_{\rm
s}}\sin\theta_{\rm s})$. exhibit two distinct wave forms: (i) The stripes
below the initial position of the resonance $(J_{\rm s}\lesssim 200\,{\rm
kpc}^{2}\,{\rm Gyr}^{-1})$ appear smooth as they originate from the
instantaneous formation of the bar. (ii) In contrast, perturbations left by
the passage of the resonance have a discontinuous form again due to the saddle
point of $\widehat{H}$.

Figure 9: Dependence of the phase-space distribution near a resonance on the
bar’s slowing rate $\eta$. Top row: the distribution after the bar has slowed
from $\Omega_{\rm p}=72$ to $36\,{\rm Gyr}^{-1}$. Middle row: time evolution
of the speed parameter $s$ (equation 36). The dotted lines mark the critical
value $s=1$ above which trapping is absent. Bottom row: time evolution of the
phase-space area of trapping represented by $J_{\ell,{\rm sep}}$ (equation
44). Decrease in $J_{\ell,{\rm sep}}$ at $\eta=0.01$ (third column) implies
resonant escape.

###  4.3 Variation with bar slowing rate

Since our method is applicable to any speed $s$, it allows us to seamlessly
connect the slow and fast limits previously explored via adiabatic theory and
linear theory, respectively. Figure 9 shows the variation of the phase-space
distribution with slowing rate $\eta$. We chose $\eta$ such that they are
approximately evenly spaced in logarithm with width
$\Delta\log_{10}\eta\approx 0.5$. The top row plots the distribution after the
bar slowed from $\Omega_{\rm p}=72$ to $36\,{\rm Gyr}^{-1}$. The left-three
figures show evolution in the slow regime, where the speed parameter $s$
remains below unity as reported in the middle row. As described in Section 3,
the size of the trapped region diminishes with increasing $\eta$. For
$\eta=0.001$ and 0.0032, we observe new trapped layers because the trapped
volume increases with time, as shown in the bottom row. For $\eta=0.01$,
however, the volume reduces, causing trapped orbits to escape the resonance
and create a dense trailing wake.

Figure 10: Phase-space distribution near a moving resonance predicted by
linear perturbation theory (equation 58). Dotted lines mark the current
position of the resonance. Linear theory qualitatively describes the near-
resonant dynamics (Fig. 9) only in the fast regime (right figure) where
libration is absent.

Whether the trapped volume increases or not depends on a number of factors, in
particular, from equation (44),

$\displaystyle\frac{\mathrm{d}J_{\ell,{\rm
sep}}}{\mathrm{d}t}=\frac{8J_{0}}{\pi}\left[\frac{1}{2}\left(\frac{\Psi^{\prime}}{\Psi}-\frac{G^{\prime}}{G}\right)C+\frac{\mathrm{d}C}{\mathrm{d}s}s^{\prime}\right],$
(53)

where the prime denotes a time derivative. Since our bar potential has a fixed
amplitude and elongates along with the expanding resonance,
$\Psi^{\prime}/\Psi$ is minor. In contrast, $-G^{\prime}/G$ is large and
positive since $G$ rapidly approaches zero as the resonance moves towards
larger $J_{\rm s}$ (Appendix B). The third term on the r.h.s of equation (53)
describes the change in the resonant-volume factor $C$ (equation 47) and is
typically negative since $\mathrm{d}C/\mathrm{d}s<0$ while $s^{\prime}>0$
(Fig. 9 middle row). For small $s$, and thus $C\sim 1$, the second term
dominates, resulting in an increasing trapped volume. However, in the
transition regime (left-third column of Fig. 9), where $s$ is close to unity
and thus $C$ is small (Fig.7), the third term becomes dominant, causing the
resonance to shrink.

The right column of Fig. 9 shows evolution in the fast regime where trapping
is absent. The striated perturbation remains, but it now has a smooth
distribution since neighbouring orbits no longer split at the resonance due to
the disappearance of the saddle points of $\widehat{H}$ (cf. Fig. 5).

###  4.4 Comparison with linear perturbation theory

The perturbed distribution in the fast regime should align with predictions
from linear perturbation theory since librating orbits are absent. To confirm
this, we derive the linear response in the slow angle-action space, allowing
for the general time-dependence of the perturber. The linearized CBE has the
following general solution (Weinberg, 2004):

$\displaystyle f_{1}({\bm{\theta}},{\bm{J}},t)$
$\displaystyle=\sum_{\bm{n}}\left[i{\bm{n}}\\!\cdot\\!\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-i{\bm{n}}\cdot\mathbf{\Omega}(t-t^{\prime})}\widehat{\Phi}_{{\bm{n}}}({\bm{J}},t^{\prime})\right]\mathrm{e}^{i{\bm{n}}\cdot{\bm{\theta}}}.$
(54)

We separate the rotational time-dependence of the perturbation
$\widehat{\Phi}_{{\bm{n}}}({\bm{J}},t)=\Phi_{{\bm{n}}}({\bm{J}},t)\mathrm{e}^{-in_{\varphi}\int^{t}_{0}\mathrm{d}t^{\prime}\Omega_{\rm
p}(t^{\prime})}$ and rewrite

$\displaystyle f_{1}({\bm{\theta}},{\bm{J}},t)$
$\displaystyle=i\sum_{\bm{n}}{\bm{n}}\\!\cdot\\!\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\mathrm{e}^{i\left[{\bm{n}}\cdot{\bm{\theta}}-n_{\varphi}\int^{t}_{0}\mathrm{d}t^{\prime}\Omega_{\rm
p}(t^{\prime})\right]}$ $\displaystyle\hskip
8.53581pt\times\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-i\int_{t^{\prime}}^{t}\mathrm{d}t^{\prime\prime}\Omega_{s}({\bm{J}},t^{\prime\prime})}\Phi_{{\bm{n}}}({\bm{J}},t^{\prime}),$
(55)

where
$\Omega_{s}({\bm{J}},t)\equiv{\bm{n}}\cdot\mathbf{\Omega}({\bm{J}})-n_{\varphi}\Omega_{\rm
p}(t)$. We transform ${\bm{\theta}}$ to the fast-slow angles and average
$f_{1}$ over the fast angles

$\displaystyle\bar{f}_{1}(\theta_{\rm s},{\bm{J}},t)$
$\displaystyle=\frac{1}{(2\pi)^{2}}\int\mathrm{d}^{2}{\bm{\theta}}_{\rm
f}~{}i\sum_{\bm{n}}{\bm{n}}\\!\cdot\\!\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}$
$\displaystyle\hskip
8.53581pt\times\mathrm{e}^{i\left[\left(n_{r}-N_{r}\frac{n_{\varphi}}{N_{\varphi}}\right)\theta_{{\rm
f}_{1}}+\left(n_{\psi}-N_{\psi}\frac{n_{\varphi}}{N_{\varphi}}\right)\theta_{{\rm
f}_{2}}+\frac{n_{\varphi}}{N_{\varphi}}\theta_{\rm s}\right]}$
$\displaystyle\hskip
8.53581pt\times\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-i\int_{t^{\prime}}^{t}\mathrm{d}t^{\prime\prime}\Omega_{s}({\bm{J}},t^{\prime\prime})}\Phi_{{\bm{n}}}({\bm{J}},t^{\prime}).$
(56)

Since $\Phi_{{\bm{n}}}$ is non-zero only when $n_{\varphi}=\pm m=\pm 2$ (Chiba
& Schönrich, 2022, Appendix B), we have for resonances with $N_{\varphi}=2$,

$\displaystyle\bar{f}_{1}(\theta_{\rm s},{\bm{J}},t)$
$\displaystyle=i\sum_{\bm{n}}{\bm{n}}\\!\cdot\\!\frac{{\partial}f_{0}}{{\partial}{\bm{J}}}\delta_{n_{r},\pm
N_{r}}\delta_{n_{\psi},\pm N_{\psi}}\delta_{n_{\varphi},\pm N_{\varphi}}$
$\displaystyle\hskip
8.53581pt\times\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{i\left[\frac{n_{\varphi}}{N_{\varphi}}\theta_{\rm
s}-\int_{t^{\prime}}^{t}\mathrm{d}t^{\prime\prime}\Omega_{s}({\bm{J}},t^{\prime\prime})\right]}\Phi_{{\bm{n}}}({\bm{J}},t^{\prime}),$
(57)

where $\delta_{n,N}$ is the Kronecker delta. Using the reality condition
$\Phi_{-{\bm{n}}}=\Phi_{{\bm{n}}}^{*}$ and the relation
${\bm{N}}\\!\cdot\left({\partial}f_{0}/{\partial}{\bm{J}}\right)={\partial}f_{0}/{\partial}J_{\rm
s}$, we finally obtain

$\displaystyle\bar{f}_{1}(\theta_{\rm s},{\bm{J}},t)$
$\displaystyle=-\frac{{\partial}f_{0}}{{\partial}J_{\rm
s}}\int_{0}^{t}\\!\mathrm{d}t^{\prime}\sin\left[\theta_{\rm s}-\theta_{\rm
s,res}-\int_{t^{\prime}}^{t}\\!\mathrm{d}t^{\prime\prime}\Omega_{s}({\bm{J}},t^{\prime\prime})\right]$
$\displaystyle\hskip 8.53581pt\times|\Phi_{{\bm{N}}}({\bm{J}},t^{\prime})|.$
(58)

Figure 10 plots the perturbed distribution $f=f_{0}+\bar{f}_{1}$ predicted by
equation (58). Comparison with Fig. 9 confirms that linear theory is
qualitatively correct in the fast regime (right figure), although the errors
are catastrophic in the slow regime (left three figures) where resonant orbits
librate. Note that the accuracy of linear theory also depends on the evolution
time, which shortens as the slowing rate increases. We have confirmed,
however, that linear theory continues to succeed in the fast regime even
beyond the typical libration period $t_{\rm lib}$ since the nonlinear response
is precluded (TW84). Hence, while linear theory works in the slow regime
($s<1$) only for a short period ($t\lesssim t_{\rm lib}$), it remains valid in
the fast regime ($s>1$) over an extended period ($t\gtrsim t_{\rm lib}$).

##  5 Density response

Although the phase-space distribution contains the full information of the
halo’s response, the density response in real space will help us understand
the nature of dynamical friction and feedback.

We calculate the density response by integrating the distribution function
over a Cartesian grid ($96\times 96\times 96$) in 3-dimensional velocity
space:

$\displaystyle\rho({\bm{x}},t)=\int\mathrm{d}^{3}{\bm{v}}f(\theta_{\rm
s},{\bm{J}},t).$ (59)

For each grid point, we compute $f$ by transforming the coordinates to the
fast-slow angle-action variables and evolving only the slow variables backward
in time until $t=0$ (equation 52). In practice, we truncate the velocity
integral at a fixed energy $E_{\rm max}=\Phi_{0}(r_{\rm max})$, where we set
$r_{\rm max}=10r_{\rm s}=200\,{\rm kpc}$. Calculations with $r_{\rm
max}>10r_{\rm s}$ made little visible difference.

Figure 11: Density wakes at the equatorial plane ($z=0$) caused by (a) the
corotation resonance ${\bm{N}}=(0,2,2)$, (b) the outer Lindblad resonance
${\bm{N}}=(1,2,2)$, and (c,d) the direct radial resonances ${\bm{N}}=(1,0,2)$
and ${\bm{N}}=(2,0,2)$. The density is calculated from the distribution
function evolved with the fast-angle averaged Hamiltonian. The bar lies on the
$x$-axis and is rotating anti-clockwise. The black circles mark the resonant
radius for circular orbits ($J_{r}=0$).

(a) Superposition of density wakes generated by the four dominant resonances
using the 1D averaged Hamiltonian.  (b) Standard test-particle simulation
using the 3D Hamiltonian in Cartesian coordinate.  (c) Residual between the 1D
model and the 3D simulation.

Figure 12: Density wake at the galactic plane calculated using the 1D averaged
Hamiltonian (top) and the 3D full Hamiltonian (middle). Bottom row shows the
residual between the two.

Figure 11 shows the density response at the galactic plane ($z=0$) caused by
the four strongest resonances: the corotation resonance ${\bm{N}}=(0,2,2)$,
the outer Lindblad resonance ${\bm{N}}=(1,2,2)$, and the direct radial
resonances ${\bm{N}}=(1,0,2)$ and ${\bm{N}}=(2,0,2)$. The bar lies on the
$x$-axis and is rotating anticlockwise. As before, we slow the bar from
$\Omega_{\rm p}=72$ to $36\,{\rm Gyr}^{-1}$.

The density response varies significantly with slowing rate $\eta$ as the
evolutionary regime makes a transition from slow to fast. For $\eta$ as small
as $0.001$ (Fig. 11 left column), the resonances mostly evolve in the slow
regime, so the trapped orbits dragged from the inner halo form a density
island migrating outwards with the resonance. These density islands are
misaligned with the bar because the center of libration $\theta_{\rm res}$
shifts anti-clockwise when the resonance is moving (see discussion surrounding
Figs. 5-7). As $\eta$ is increased from left to right, the density island
rotates further. Eventually, the evolution shifts to the fast regime, and the
density response deforms into a striated pattern – a density wake – that no
longer drifts with the resonance. Note that slow and fast regimes typically
coexist since orbits crossing the galactic plane take a wide range of
${\bm{J}}_{\rm f}$ and thus pass the resonance at different speed $s$.

It is worth noting that the density response in the fast/linear regime is
mostly in phase with the bar, implying that the halo’s self-gravity will
amplify the bar’s perturbation (Dootson & Magorrian, 2022). However, the
density response gradually becomes out of phase with the bar as the evolution
becomes slow/nonlinear. This implies that, in the slow regime, the self-
gravity of the halo may weaken the bar’s perturbation, causing the bar to
evolve even slower.

Figure 12 compares the superposition of the density response induced by the
four strongest resonances (top row) with the 3D test-particle simulation
(middle row). The superposition of 1D averaged models reproduces almost all
features of the 3D simulation. This result confirms that the density response
for each resonance has been computed correctly and that the assumption behind
the method of averaging (i.e. resonances are well separated) is acceptable.
The bottom row shows the residual between the 1D model and the 3D simulation.
At small $\eta$ (left column), the residual shows weakly negative (blue)
features in the trapped region which indicate that the trapped mass in the
full simulation is slightly smaller than in the averaged model. This is likely
due to the partial overlap of resonances (Fig.8 and Appendix F of Chiba &
Schönrich, 2022), where the chaotic behavior of orbits can hinder them from
remaining trapped on the moving resonances. The residual is also non-
negligible near the galactic centre, where all resonances induce coherent
quadrupole perturbations, implying that adding contributions from other minor
resonances may further improve the model.

##  6 Dynamical feedback

In this section, we formulate dynamical feedback introduced by TW84. As
described in Section 4.1.2, dynamical feedback is the transfer of momentum
that arises directly from the migration of the trapped phase-space. Suppose,
for instance, that the bar slowed down by dynamical friction and, as a result,
the resonance shifted towards larger angular momentum. The orbits trapped in
resonance will then be dragged along with the resonance, carrying further
angular momentum away from the bar. Dynamical feedback refers to this feedback
torque. Dynamical feedback only occurs in the slow nonlinear regime, where
there is a finite volume of trapping, and is thus not predicted by standard
linear perturbation theory (e.g. Lynden-Bell & Kalnajs, 1972; Weinberg, 2004).

Dynamical feedback directly depends on the rate of change of the bar’s pattern
speed which has an interesting consequence. To see this, let us write
dynamical feedback as $\tau_{\rm fb}=I_{\rm fb}\dot{\Omega}_{\rm p}$, where
$I_{\rm fb}$ is a coefficient that has the dimension of a moment of inertia,
and denote the usual dynamical friction as $\tau_{\rm fr}$. The equation of
motion of the bar $I\dot{\Omega}_{\rm p}=\tau_{\rm fr}+\tau_{\rm fb}$ (Section
2) can then be written as

$\displaystyle(I-I_{\rm fb})\dot{\Omega}_{\rm p}=\tau_{\rm fr},$ (60)

which implies that dynamical feedback changes the moment of inertia of the
perturber (Weinberg, 1985). If $I_{\rm fb}>0$, dynamical feedback helps the
bar spin down (a positive feedback); otherwise, trapping stabilizes the bar’s
pattern speed (a negative feedback). As we shall see, dynamical feedback by
the trapped orbits is typically positive.

The original study by TW84 in fact did not concern dynamical feedback by the
trapped orbits. Instead, they studied feedback from the untrapped orbits
passed by a slowly moving resonance. The feedbacks from trapped and untrapped
orbits are in fact closely related through the collisionless Boltzmann
equation, which states that phase-space flow is incompressible: When the
trapped island moves, the untrapped orbits in its path are forced to skirt
around it and flow to the far side of the resonance, as illustrated in Fig. 5.
Therefore, an acceleration (deceleration) of the trapped orbits must be
accompanied by a deceleration (acceleration) of the surrounding untrapped
orbits. The net dynamical feedback is determined by the imbalance between
these two effects.

Here, we provide the full formula of dynamical feedback, incorporating
contributions from both trapped and untrapped orbits.

###  6.1 Dynamical feedback by trapped orbits

We begin with the dynamical feedback by the trapped orbits. Since dynamical
feedback is the torque that arises purely from the movement of the trapped
phase-space, we must first remove from the distribution function any transient
perturbations that phase mixes over time. This is effected by averaging the
distribution function over the angle of libration

$\displaystyle\tilde{f}_{\ell}({\bm{J}}_{\rm
f},J_{\ell},t)=\frac{1}{2\pi}\int\mathrm{d}\theta_{\ell}f_{\ell}(\theta_{\ell},{\bm{J}}_{\rm
f},J_{\ell},t).$ (61)

We multiply $\tilde{f}_{\ell}$ with minus the torque on individual halo orbits

$\displaystyle\dot{L}_{z}=N_{\varphi}\left(\dot{J}_{\rm
s,res}+\dot{\Delta}\right)$ (62)

and integrate over the slow phase-space $\mathrm{d}\theta_{\rm
s}\mathrm{d}J_{\rm s}=\mathrm{d}\theta_{\ell}\mathrm{d}J_{\ell}$

$\displaystyle\hat{\tau}_{\ell}\left({\bm{J}}_{\rm
f},t\right)=-N_{\varphi}\int_{0}^{J_{\ell,{\rm
sep}}}\\!\mathrm{d}J_{\ell}\int_{0}^{2\pi}\\!\mathrm{d}\theta_{\ell}\tilde{f}_{\ell}\left({\bm{J}}_{\rm
f},J_{\ell},t\right)\left(\dot{J}_{\rm s,res}+\dot{\Delta}\right).$ (63)

This is the torque felt by the bar due to orbits trapped in resonance
${\bm{N}}$ at a specific ${\bm{J}}_{\rm f}$. Since $\dot{\Delta}$ is $2\pi$
periodic in $\theta_{\ell}$, the second term vanishes, leaving

$\displaystyle\hat{\tau}_{\ell}\left({\bm{J}}_{\rm f},t\right)=-\frac{2\pi
N_{\varphi}^{2}\dot{\Omega}_{\rm p}}{G({\bm{J}}_{\rm f},J_{\rm
s,res})}\int_{0}^{J_{\ell,{\rm
sep}}}\\!\mathrm{d}J_{\ell}\tilde{f}_{\ell}\left({\bm{J}}_{\rm
f},J_{\ell},t\right),$ (64)

where we have plugged in equation (38). We remind the reader that $G$ is
${\partial}({\bm{N}}\cdot\mathbf{\Omega})/{\partial}J_{\rm s}$ evaluated at
the resonance $J_{\rm s}=J_{\rm s,res}$ (equation 33), and is not the
gravitational constant.

Depending on the geometry of the resonance, there could exist more than one
trapped region at a fixed ${\bm{J}}_{\rm f}$ satisfying the same resonance
condition: the resonance curve on the $(L,J_{r})$ plane may have multiple
intersections with the planes of constant ${\bm{J}}_{\rm f}$. Such cases are
limited but found, for instance, at the corotation resonance at small $L$ and
large $J_{r}$ (Appendix B). Hence we write

$\displaystyle\hat{\tau}_{\ell}\left({\bm{J}}_{\rm f},t\right)=-2\pi
N_{\varphi}^{2}\dot{\Omega}_{\rm p}\sum_{i}\frac{B_{\ell}({\bm{J}}_{\rm
f},J_{\rm s,res}^{i},t)}{G({\bm{J}}_{\rm f},J_{\rm s,res}^{i})},$ (65)

where the function

$\displaystyle B_{\ell}({\bm{J}}_{\rm f},J_{\rm
s,res}^{i},t)=\int_{0}^{J_{\ell,{\rm
sep}}}\\!\mathrm{d}J_{\ell}\tilde{f}_{\ell}\left({\bm{J}}_{\rm
f},J_{\ell},t\right)$ (66)

quantifies the amount of trapped orbits near each resonance $J_{\rm
s,res}^{i}$. Equation (65) can be formally rewritten as

$\displaystyle\hat{\tau}_{\ell}\left({\bm{J}}_{\rm f},t\right)=-2\pi
N_{\varphi}^{2}\dot{\Omega}_{\rm p}\\!\int\\!\mathrm{d}J_{\rm
s}\mathrm{sgn}(G)B_{\ell}({\bm{J}}_{\rm f},J_{\rm
s},t)\delta\\!\left({\bm{N}}\\!\cdot\\!\mathbf{\Omega}-N_{\varphi}\Omega_{\rm
p}\right),$ (67)

using the identity
$\sum_{i}f(x_{i})/|g^{\prime}(x_{i})|=\int\mathrm{d}xf(x)\delta(g(x))$ (e.g.
Kanwal, 1998), where $\delta$ is the Dirac’s delta function.

To get the total dynamical feedback by trapped orbits, we integrate
$\hat{\tau}_{\ell}$ over the fast angle-actions and sum over all resonances

$\displaystyle\tau_{\ell}$
$\displaystyle=\sum_{{\bm{N}}}\int\mathrm{d}^{2}{\bm{J}}_{\rm
f}~{}\mathrm{d}^{2}{\bm{\theta}}_{\rm
f}~{}\hat{\tau}_{\ell}\left({\bm{J}}_{\rm f},t\right)$
$\displaystyle=-(2\pi)^{3}\dot{\Omega}_{\rm
p}\\!\sum_{{\bm{N}}}N_{\varphi}^{2}\\!\int\\!\mathrm{d}^{3}{\bm{J}}^{\prime}\mathrm{sgn}(G)B_{\ell}({\bm{J}}^{\prime}\\!,t)\delta\\!\left({\bm{N}}\\!\cdot\\!\mathbf{\Omega}-N_{\varphi}\Omega_{\rm
p}\right),$ (68)

where we used the notation ${\bm{J}}^{\prime}=({\bm{J}}_{\rm f},J_{\rm s})$.
Switching variables from ${\bm{J}}^{\prime}$ to ${\bm{J}}=(J_{r},L,L_{z})$,
where the Jacobian determinant is
$\left|{\partial}({\bm{J}}^{\prime})/{\partial}({\bm{J}})\right|=1/N_{\varphi}$,
we finally obtain

$\displaystyle\tau_{\ell}=-(2\pi)^{3}\dot{\Omega}_{\rm
p}\\!\sum_{{\bm{N}}}N_{\varphi}\\!\\!\int\\!\\!\mathrm{d}^{3}{\bm{J}}\mathrm{sgn}(G)B_{\ell}({\bm{J}},t)\delta\\!\left({\bm{N}}\\!\cdot\\!\mathbf{\Omega}-N_{\varphi}\Omega_{\rm
p}\right).$ (69)

This equation describes the dynamical feedback from all the trapped orbits
dragged by the moving resonances: the function $B_{\ell}$ measures the mass of
trapped orbits, and the $\mathrm{sgn}(G)$ encodes the direction at which they
move in $L_{z}$ in response to the change in pattern speed. In the absence of
resonance migration ($\dot{\Omega}_{\rm p}=0$), dynamical feedback vanishes as
expected. Since $B_{\ell}\geq 0$ by definition and $G<0$ for most cases
(Appendix B), trapped orbits generally provide a net positive feedback on the
perturber, thereby decreasing the perturber’s moment of inertia.

Figure 13: Left: model for the phase-mixed distribution function (equation
77). Middle: original distribution function (equation 52). Right: residual
between the mixed model and full distribution. For all panels, the bar slowed
from $\Omega_{\rm p}=72$ to $36\,{\rm Gyr}^{-1}$ with slowing rate
$\eta=0.001$.

###  6.2 Dynamical feedback by untrapped orbits: the TW formula

We next formulate dynamical feedback by the untrapped orbits. As in the
previous section, we first aim at removing the transient perturbations from
the distribution function by averaging it over the angle. Unlike the trapped
orbits, however, the angle of circulation for the untrapped orbits cannot be
defined from their trajectories in the comoving frame since the orbits are not
closed (cf. Fig. 5). The motion of circulation can instead be described by
expanding perturbation series about the unperturbed orbit in the static frame.
Perturbation theory then predicts that changes to the angle-averaged
distribution function are of second order in the strength of perturbation $A$
(e.g. Fouvry et al., 2015a). Hence, to first order in $A$, the phase-mixed
distribution of the untrapped (circulating) orbits $\tilde{f}_{\rm c}$ can be
approximated by the unperturbed distribution $f_{0}$. We will demonstrate
below that this approximation leads to the formula derived by TW84.

With the aforementioned approximation, $\tilde{f}_{\rm
c}=f_{0}+\mathcal{O}(A^{2})$, the net torque on the bar by the untrapped halo
orbits can be obtained by multiplying $f_{0}$ with $-\dot{L}_{z}$ and
integrating it over the slow phase-space
$\mathrm{d}^{2}S=\mathrm{d}\theta_{\rm s}\mathrm{d}J_{\rm s}$, while excluding
the region of libration $S_{\ell}$ from the domain of integration:

$\displaystyle\hat{\tau}_{\rm c}\left({\bm{J}}_{\rm f},t\right)=$
$\displaystyle-\int_{S-S_{\ell}}\mathrm{d}^{2}Sf_{0}({\bm{J}}_{\rm f},J_{\rm
s})\dot{L}_{z}$ $\displaystyle=$
$\displaystyle-\int_{S}\mathrm{d}^{2}Sf_{0}({\bm{J}}_{\rm f},J_{\rm
s})\dot{L}_{z}+\int_{S_{\ell}}\mathrm{d}^{2}Sf_{0}({\bm{J}}_{\rm f},J_{\rm
s})\dot{L}_{z}.$ (70)

The first term in the second line is zero since the torque is periodic in the
slow angle $\theta_{\rm s}$. We next Taylor expand $f_{0}$ in $J_{\rm s}$ to
first order around the resonance $J_{\rm s,res}$

$\displaystyle\hat{\tau}_{\rm c}\left({\bm{J}}_{\rm
f},t\right)=\int_{S_{\ell}}\mathrm{d}^{2}S\left[f_{0}({\bm{J}}_{\rm f},J_{\rm
s,res})+\left.\frac{{\partial}f_{0}}{{\partial}J_{\rm s}}\right|_{J_{\rm
s,res}}\Delta\right]\dot{L}_{z},$ (71)

where $\Delta=J_{\rm s}-J_{\rm s,res}$. We similarly expand the Hamiltonian
(equation 32) which makes the second term vanish since the domain of
integration $S_{\ell}$ becomes symmetric about $\Delta=0$, and $\Delta$ is odd
while $\dot{L}_{z}$ is even in $\Delta$. For the remaining first term, we can
follow the same procedure as in the previous section beginning with equation
(63). The dynamical feedback by all the untrapped orbits is then

$\displaystyle\tau_{\rm c}=(2\pi)^{3}\dot{\Omega}_{\rm
p}\sum_{{\bm{N}}}N_{\varphi}\\!\int\\!\mathrm{d}^{3}{\bm{J}}\mathrm{sgn}(G)B_{0}({\bm{J}},t)\delta\\!\left({\bm{N}}\cdot\mathbf{\Omega}-N_{\varphi}\Omega_{\rm
p}\right),$ (72)

where the function

$\displaystyle B_{0}({\bm{J}},t)=J_{\ell,{\rm
sep}}({\bm{J}},t)f_{0}({\bm{J}})$ (73)

describes the mass of untrapped orbits displaced by the trapped orbits.

Substituting $B_{0}$ to equation (72), and using equation (44), we obtain

$\displaystyle\tau_{\rm c}=(2\pi)^{3}\dot{\Omega}_{\rm
p}\\!\sum_{{\bm{N}}}N_{\varphi}\\!\\!\int\\!\\!\mathrm{d}{\bm{J}}\mathrm{sgn}(G)J_{0}Q(s)f_{0}({\bm{J}})\delta\\!\left({\bm{N}}\\!\cdot\\!\mathbf{\Omega}\\!-\\!N_{\varphi}\Omega_{\rm
p}\right),$ (74)

where $Q(s)=8C(s)/\pi$ (equation 45) and $J_{0}=\sqrt{F/G}=\sqrt{\Psi/|G|}$
(equation 34). Equation (74) is the formula of dynamical feedback derived by
TW84 (equations 86,100) and summarized in Weinberg (1985) (equation 39).

###  6.3 Total dynamical feedback

The formulae derived in the previous sections can be combined to give the
total dynamical feedback

$\displaystyle\tau_{\rm fb}$ $\displaystyle=\tau_{\ell}+\tau_{\rm c}$
$\displaystyle=(2\pi)^{3}\dot{\Omega}_{\rm
p}\sum_{{\bm{N}}}N_{\varphi}\\!\int\\!\mathrm{d}^{3}{\bm{J}}\mathrm{sgn}(G)\left(B_{0}-B_{\ell}\right)\delta\\!\left({\bm{N}}\\!\cdot\\!\mathbf{\Omega}-N_{\varphi}\Omega_{\rm
p}\right).$ (75)

This equation states that the strength of dynamical feedback is determined by
the difference between the mass of the trapped orbits $B_{\ell}$ and the mass
of the background orbits $B_{0}$ that would otherwise occupy the trapped
volume. This is a consequence of the incompressibility of phase space and is
similar to the physics of buoyancy – the trapped region is ‘submerged’ in the
phase space of untrapped orbits, and the torque required to move the trapped
region is proportional to the density difference between the two regions.

The phase-space density of a dark halo in equilibrium typically declines with
increasing angular momentum. Hence, so long as the resonance sweeps the phase
space monotonically, the gainers of angular momentum always surpass the
losers, resulting in a net negative torque on the bar. This negative torque
provides a positive feedback when the bar decelerates while it provides a
negative feedback when the bar accelerates. Thus dynamical feedback makes the
bar easier to spin down than to spin up. If the bar switched between
acceleration and deceleration in mid-course, dynamical feedback may provide a
positive torque, resulting in the opposite conclusion. Generally, the
stability condition for the bar’s rotation depends on (i) the sign of the
moment of inertia of the bar $I$, which was assumed positive in the discussion
above, (ii) the sign of $G$ (equation 34), and (iii) whether the trapped
phase-space is overdense or not, $B_{0}-B_{\ell}$. Thus

$\displaystyle\begin{cases}\mathrm{sgn}(I)=\mathrm{sgn}\left[G\left(B_{0}-B_{\ell}\right)\right]~{}~{}~{}{\rm
destabilizing~{}(pos.~{}feedback)}\\\
\mathrm{sgn}(I)=-\mathrm{sgn}\left[G\left(B_{0}-B_{\ell}\right)\right]~{}~{}~{}{\rm
stabilizing~{}(neg.~{}feedback)}\end{cases}$ (76)

Our stability condition differs from that of TW84 (equation 104), which does
not include feedback from trapped orbits, so $B_{0}-B_{\ell}$ is always
positive. Since we expect $B_{0}-B_{\ell}<0$ for a slowing bar, and $I>0$
while $G<0$ in most cases (Appendix B), our revised stability condition
implies that dynamical feedback will destabilize the bar, as we demonstrate
below.

To calculate equation (75), we need a model for the phase-mixed distribution
function of the trapped region $\tilde{f}_{\ell}$ (61). For this, we extend
the model of Sridhar & Touma (1996) which provides a prescription for
$\tilde{f}_{\ell}$ in the extreme adiabatic limit $\eta/A\rightarrow 0$
(Section 3.1). In essence, the evolution of $\tilde{f}_{\ell}(J_{\ell},t)$ is
assumed adiabatic everywhere except at the separatrix. If the separatrix
expands, orbits are captured into resonance, and the value of
$\tilde{f}_{\ell}$ at the expanding separatrix is given by the average density
of the newly captured orbits. If the separatrix shrinks, orbits at the
separatrix escape the resonance, so the domain of $\tilde{f}_{\ell}$ simply
diminishes but the distribution over the remaining domain is unaltered.

We extend this model to the adiabatic evolution in the comoving frame
$\eta\gamma/A\rightarrow 0$ (Section 3.2) where the separatrix is deformed. In
the extreme adiabatic limit, $\tilde{f}_{\ell}$ at the expanding separatrix
can be determined from the angle-averaged density of the untrapped orbits
$\tilde{f}_{\rm c}$ at the separatrices (see Sridhar & Touma, 1996, Table 1).
However, as mentioned in Section 6.2, $\tilde{f}_{\rm c}$ is not available in
the comoving frame since untrapped orbits are not closed and thus the angles
cannot be defined. As a crude approximation, we take $\tilde{f}_{\ell}$ at the
expanding separatrix to be the unperturbed density at the position of the
resonance $J_{\rm s,res}$ at the time of capture $t_{\rm cap}$:

$\displaystyle\tilde{f}_{\ell}({\bm{J}}_{\rm f},J_{\ell})\simeq
f_{0}({\bm{J}}_{\rm f},J_{\rm s,res}[t_{\rm cap}(J_{\ell})]).$ (77)

Figure. 13 compares this model (left) with the original distribution (middle)
for $\eta=0.001$. Our model $\tilde{f}_{\ell}$ predicts a flat distribution at
the core because capture at the time of bar formation occurs while $J_{\rm
s,res}$ is fixed. The real phase-mixed distribution will not be flat unless
${\partial}f/{\partial}J_{\rm s}$ is constant across the trapped region.
Outside the core, our model $\tilde{f}_{\ell}$ declines smoothly towards the
surface. The residual between the model and the original $f$ (Fig. 13, right)
exhibits the phase spiral with little offset, confirming that the model is
reasonable.

Figure 14: Time evolution of the dynamical feedback (75) by the corotation
resonance. The slowing rate $\eta$ increases from 0.001 (blue) to 0.01 (red)
with equal spacing 0.00025.

Using this model for $\tilde{f}_{\ell}$, we compute the dynamical feedback
(75). Figure 14 shows the dynamical feedback by the corotation resonance as a
function of time for various slowing rate $\eta$, ranging from $0.001$ (blue)
to $0.01$ (red). As before, we fix the bar’s initial and final pattern speeds
($\Omega_{\rm p}=72$ to $36\,{\rm Gyr}^{-1}$), so the time the bar takes to
reach its final pattern speed (black dot) gets shorter with increasing $\eta$.
The feedback torque is initially zero since $B_{\ell}=B_{0}$, but decreases
with time, becoming more negative, as the density contrast between the trapped
and untrapped region gets higher (cf. Fig. 8). Interestingly, feedback takes a
peak at an intermediate slowing rate $\eta\sim 0.0025$.

Figure 15: Dynamical feedback $\hat{\tau}$ from orbits with a specific
${\bm{J}}_{\rm f}$ as a function of the speed of the resonance $s$.

To understand this, we analyse the feedback $\hat{\tau}_{\rm fb}$ from orbits
near a single point of resonance at a specific ${\bm{J}}_{\rm f}$. From
equation (65), we have

$\displaystyle\hat{\tau}_{\rm fb}$
$\displaystyle=\hat{\tau}_{\ell}+\hat{\tau}_{\rm c}$ $\displaystyle=2\pi
N_{\varphi}^{2}\dot{\Omega}_{\rm
p}(B_{0}-B_{\ell})/G=\mathrm{sgn}(\dot{\Omega}_{\rm
p}G)N_{\varphi}\Psi(B_{0}-B_{\ell})s,$ (78)

where $s=|N_{\varphi}\dot{\Omega}_{\rm p}/(G\Psi)|$ is the speed parameter
(equation 36), allowing for arbitrary signs in $\dot{\Omega}_{\rm p}$ and $G$.
Let us approximate $B_{\ell}$ as $J_{\ell,{\rm sep}}f_{\ell
0}=(8J_{0}/\pi)C(s)f_{\ell 0}$, where $f_{\ell 0}$ is the averaged density of
the librating region, and write

$\displaystyle\hat{\tau}_{\rm fb}=\hat{\tau}_{0}sC(s),~{}~{}{\rm
where}~{}~{}\hat{\tau}_{0}=\mathrm{sgn}(\dot{\Omega}_{\rm
p}G)N_{\varphi}\Psi(8J_{0}/\pi)(f_{0}-f_{\ell 0}).$ (79)

The equation indicates that feedback from each point on the resonance scales
as $sC(s)$. Figure 15 shows this torque $\hat{\tau}_{\rm fb}$ as a function of
$s$. At small speed $s$, the torque rises with increasing $s$ because the
resonance shifts in angular momentum at a faster rate. As $s$ increases
further, the torque saturates at $s\sim 0.4$ and starts declining because the
trapped phase-space shrinks with $s$. The feedback drops to zero at $s=1$
where trapping vanishes entirely. The non-monotonic dependence of the total
feedback on $\eta$ (Fig. 14) can be understood as the integration of this
curve over the entire resonant surface.

Figure 16: Comparison between dynamical feedback (smooth line) and total
torque (jagged line) by the corotation resonance for $\eta=0.001$.

Figure 16 compares dynamical feedback with the total torque by the corotation
resonance. The latter is computed by sampling particles from $f_{0}$ and
evolving their slow variables forward in time. The total torque oscillates at
early times as orbits captured at the time of bar formation collectively
librate and phase mix (Chiba & Schönrich, 2022). As the bar slows and grows in
length, both the total torque and the dynamical feedback get stronger.

Figure 17: Dynamical feedback by the four strongest resonances for
$\eta=0.001$.

Figure 17 compares dynamical feedback from different resonances. The
corotation resonance is by far the dominant source of dynamical feedback,
overwhelming others by more than a factor of 4.

Figure 18: Comparison between the total dynamical feedback (smooth lines) and
the total torque (jagged lines) computed from 3D test-particle simulations.

Finally, Fig. 18 compares the total dynamical feedback (sum of the four
strongest resonances) with the total torque calculated from 3D test-particle
simulations. The result suggests that, for small $\eta$ of order $10^{-3}$
(black and blue), around 30% of the torque is composed of dynamical feedback.
For large $\eta$ $(\gtrsim 10^{-2})$ in which the fast regime prevails (red
and green), dynamical feedback is absent and most of the torque is due to
dynamical friction.

###  6.4 Dynamical feedback in infinite homogeneous system

Dynamical feedback also applies to a point-mass perturber moving through an
infinite homogeneous distribution of lighter particles. Here, trapped orbits
correspond to those gravitationally bound to the moving perturber. When the
perturber decelerates due to dynamical friction101010The phase-space flow in a
decelerating Kepler potential is studied by Namouni & Guzzo (2007), it drags
the trapped particles along with it. Hence, the trapped particles act to
increase the effective inertia of the perturber. However, as described in the
previous sections, conservation of phase-space volume requires that the
deceleration of the trapped particles be accompanied by an acceleration of the
surrounding untrapped particles, and it is the imbalance between the two that
determines the net dynamical feedback. Let us consider the case in which the
perturber is decelerating monotonically in its direction of motion along the
$x$-axis, so ${\bm{v}}_{\rm p}=(v_{{\rm p}x},0,0)$, and the particles
initially have an isotropic velocity distribution with a negative gradient
${\partial}f_{0}/{\partial}v_{x}<0$ for all $v_{x}>0$. Then, the phase-space
density of the untrapped region at the current velocity of the perturber
$v_{{\rm p}x}$, is larger than the phase-space density of the trapped region,
which consists of orbits captured at $v_{x}>v_{{\rm p}x}$. Hence, the momentum
gained by the untrapped orbits exceeds the momentum loss of the trapped
orbits. Consequently, particles exert a net negative force (a positive
feedback) on the perturber, meaning that the effective inertia of the
perturber in fact decreases. If the perturber accelerates instead, its
effective inertia will increase. The conclusion is the same as in the bar-halo
system. The classical formula of dynamical friction by Chandrasekhar (1943)
fails to capture this phenomenon because it is based upon the assumptions that
the perturber has a constant velocity and that all field particles move pass
the perturber on Keplerian hyperbolae, that is, all orbits are assumed to be
unbound (untrapped).

##  7 Discussion and conclusion

We studied dynamical friction in the general fast-slow regime, where neither
linear perturbation (fast limit) theory nor classical adiabatic (slow limit)
theory is applicable. This intermediate regime is expected in realistic galaxy
evolution, but has thus far eluded proper description.

We first showed that the phase-space dynamics around a drifting resonance is
naturally described in the comoving frame of the resonance since the adiabatic
condition there is better fulfilled than in the static frame. The phase-space
flow in the comoving frame is determined by a single dimensionless parameter
$s$, which is the ratio of the characteristic libration time to the
resonance’s migration time. We then solved the evolution of the distribution
function numerically using the angle-averaged Hamiltonian, and provided an
explanation to dynamical friction for a range of $s$, linking the theories
developed earlier in the slow ($s=0$) and fast limits ($s\gg 1$). The
mechanism of dynamical friction can be summarized as follows:

(i) In the slow limit $(s=0)$, the phase-space flow is similar to that of a
pendulum, where trapped orbits librate inside the separatrix, while untrapped
orbits circulate outside. As particles move along these orbits, they exchange
angular momentum with the perturber. Dynamical friction arises from the
imbalance between the gainers and losers, and is initially non-zero due to the
negative gradient in the distribution function, but undergoes damped
oscillations as the distribution phase mixes (Chiba & Schönrich, 2022; Banik &
van den Bosch, 2022).

(ii) In the general slow regime $(0<s<1)$, the trapped phase-space comove with
the resonance, while its volume contracts. This contraction in the trapped
volume allows the untrapped orbits to become open in phase space, transiting
from one side of the resonance to the other. In this regime, in addition to
the dynamical friction described above, the moving resonances give rise to
dynamical feedback, which refers to the transfer of angular momentum directly
caused by the movement of the resonances. For example, a decelerating bar
typically drives the resonances towards larger angular momentum $L_{z}$, thus
yielding $L_{z}$ to the trapped orbits while receiving $L_{z}$ from the
untrapped orbits that leap over the drifting resonance. Due to the negative
gradient in the distribution function, gainers of $L_{z}$ generally outnumber
the losers, resulting in a net negative torque on the perturber. This argument
holds irrespective of the sign of the change in pattern speed
$\dot{\Omega}_{\rm p}$. Dynamical feedback therefore facilitates the
deceleration of the perturber (positive feedback) but hinders its acceleration
(negative feedback). Since the phase-space volume diminishes in size with
increasing $s$, dynamical feedback peaks at $s\sim 0.4$ and decays towards
$s=1$.

(iii) In the general fast regime $(s>1)$, the trapped region vanishes, and so
does dynamical feedback. Here, angular momentum is transferred only through
dynamical friction by the untrapped orbits. These untrapped orbits freely flow
across the resonance and generate a striated perturbation downstream of the
resonance. Since non-linear libration is absent in this regime, the perturbed
distribution can be reproduced qualitatively by standard linear perturbation
theory.

(iv) In the fast limit $(s\gg 1)$, the error of the linear solution diminishes
since $s$ scales inversely with the perturbation amplitude. This is in line
with the conclusion of TW84.

Dynamical feedback was first predicted by TW84 although their description was
incomplete as they only indicated feedback by the untrapped orbits. By
identifying dynamical feedback as the torque arising from the phase-mixed
distribution function, we derived the full analytic expression for dynamical
feedback, accounting for both trapped and untrapped orbits. The untrapped
component of our formula is shown to be equivalent to the equation derived in
TW84. Our full formula highlights that dynamical feedback from each vicinity
of the resonance is proportional to the difference between the mass of the
trapped orbits and that of the untrapped orbits displaced by the trapped
orbits. This follows from the conservation of phase-space volume and is
reminiscent of the Archimedes’ principle in fluid dynamics – just like the
force required to lift an object submerged in a fluid, the torque required to
shift the trapped island in phase space scales with the difference between the
density of the trapped island and that of the background. Using our formula,
we showed that dynamical feedback can comprise up to 30% of the total torque
on the Galactic bar, if the bar has been evolving with $\eta\sim 0.003$ as
measured by Chiba et al. (2021).

Real galaxies are far more complex than the simple bar-halo system studied
here: there are various sources of time-dependent perturbations (e.g. spiral
arms, molecular clouds, orbiting substructures) that may stochastically
scatter orbits and disrupt the resonant phase-space. The resulting diffusion
in phase space is likely to sustain dynamical friction by recovering the
gradient in the distribution function (Hamilton et al., 2022), but is expected
to suppress dynamical feedback as it acts to reduce the density difference
between the trapped and untrapped regions. On the other hand, the self-gravity
of the halo’s response, which is neglected in this study, may coherently
strengthen or weaken the bar’s perturbation (e.g. Weinberg, 1989; Chavanis,
2012). Linear theory predicts that self-gravity enhances dynamical friction on
bars by almost a factor of two Dootson & Magorrian (2022). This is in line
with the halo’s density response (Fig. 11) which is largely in phase with the
bar in the fast-linear regime. In the slow-nonlinear regime, however, self-
gravity may weaken both friction and feedback as the density response becomes
out of phase with the bar. Work is in progress to quantify these effects.

## Acknowledgements

I thank J-B. Fouvry, N. Frankel, S. Hadden, C. Hamilton, J. Hunt, N. Murray,
S. Tremaine and members of the NAOJ for valuable discussions. I am also
grateful to my mentors R. Schönrich, J. Magorrian and J. Binney for their
helpful advice and encouragement. This work was supported by the Royal Society
grant RGF$\backslash$R1$\backslash$180095 and the Natural Sciences and
Engineering Research Council of Canada (NSERC), [funding reference
#DIS-2022-568580].

## Data availability

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

## References

  * Armon & Friedland (2016) Armon T., Friedland L., 2016, Phys. Rev. A, 93, 043406
  * Artemyev et al. (2018) Artemyev A. V., Neishtadt A. I., Vainchtein D. L., Vasiliev A. A., Vasko I. Y., Zelenyi L. M., 2018, Communications in Nonlinear Science and Numerical Simulations, 65, 111
  * Asano et al. (2020) Asano T., Fujii M. S., Baba J., Bédorf J., Sellentin E., Portegies Zwart S., 2020, MNRAS, 499, 2416
  * Athanassoula (1996) Athanassoula E., 1996, in Buta R., Crocker D. A., Elmegreen B. G., eds, Astronomical Society of the Pacific Conference Series Vol. 91, IAU Colloq. 157: Barred Galaxies. p. 309
  * Athanassoula (2003) Athanassoula E., 2003, MNRAS, 341, 1179
  * Athanassoula (2014) Athanassoula E., 2014, MNRAS, 438, L81
  * 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, ApJ, 912, 43
  * Banik & van den Bosch (2022) Banik U., van den Bosch F. C., 2022, ApJ, 926, 215
  * Beane et al. (2022) Beane A., et al., 2022, arXiv e-prints, p. arXiv:2209.03364
  * 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
  * Bournaud et al. (2005) Bournaud F., Combes F., Semelin B., 2005, MNRAS, 364, L18
  * Bovy et al. (2019) Bovy J., Leung H. W., Hunt J. A. S., Mackereth J. T., García-Hernández D. A., Roman-Lopes A., 2019, MNRAS, 490, 4740
  * Carlberg & Sellwood (1985) Carlberg R. G., Sellwood J. A., 1985, ApJ, 292, 79
  * Chandrasekhar (1943) Chandrasekhar S., 1943, ApJ, 97, 255
  * Chavanis (2012) Chavanis P.-H., 2012, Physica A Statistical Mechanics and its Applications, 391, 3680
  * Chiba & Schönrich (2021) Chiba R., Schönrich R., 2021, MNRAS, 505, 2412–2426
  * Chiba & Schönrich (2022) Chiba R., Schönrich R., 2022, MNRAS, 513, 768
  * Chiba et al. (2021) Chiba R., Friske J. K. S., Schönrich R., 2021, MNRAS, 500, 4710
  * Clarke & Gerhard (2022) Clarke J. P., Gerhard O., 2022, MNRAS, 512, 2171
  * Debattista & Sellwood (2000) Debattista V. P., Sellwood J. A., 2000, ApJ, 543, 704
  * Dehnen (2000) Dehnen W., 2000, AJ, 119, 800
  * Dootson & Magorrian (2022) Dootson D., Magorrian J., 2022, to be submitted
  * Fouvry et al. (2015a) Fouvry J.-B., Pichon C., Prunet S., 2015a, MNRAS, 449, 1967
  * Fouvry et al. (2015b) Fouvry J. B., Pichon C., Magorrian J., Chavanis P. H., 2015b, A&A, 584, A129
  * Fujii et al. (2019) Fujii M. S., Bédorf J., Baba J., Portegies Zwart S., 2019, MNRAS, 482, 1983
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2206.06207
  * Ghafourian et al. (2020) Ghafourian N., Roshan M., Abbassi S., 2020, ApJ, 895, 13
  * Hairer et al. (2006) Hairer E., Hochbruck M., Iserles A., Lubich C., 2006, Oberwolfach Reports, 3, 805
  * Halle et al. (2018) Halle A., Di Matteo P., Haywood M., Combes F., 2018, A&A, 616, A86
  * Hamilton et al. (2022) Hamilton C., Tolman E. A., Arzamasskiy L., Duarte V. N., 2022, arXiv e-prints, p. arXiv:2208.03855
  * 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
  * Hopkins & Quataert (2011) Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027
  * Hui et al. (2017) Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, Phys. Rev. D, 95, 043541
  * Kalyakin (2008) Kalyakin L. A., 2008, Russian Mathematical Surveys, 63, 791
  * Kanwal (1998) Kanwal R. P., 1998, Generalized functions theory and technique: Theory and technique. Springer Science & Business Media
  * Kaur & Stone (2022) Kaur K., Stone N. C., 2022, MNRAS, 515, 407
  * Kawata et al. (2021) Kawata D., Baba J., Hunt J. A. S., Schönrich R., Ciucă I., Friske J., Seabroke G., Cropper M., 2021, MNRAS, 508, 728
  * Lancaster et al. (2020) Lancaster L., Giovanetti C., Mocz P., Kahn Y., Lisanti M., Spergel D. N., 2020, J. Cosmology Astropart. Phys., 2020, 001
  * Leung et al. (2022) Leung H. W., Bovy J., Mackereth J. T., Hunt J. A. S., Lane R. R., Wilson J. C., 2022, arXiv e-prints, p. arXiv:2204.12551
  * Li et al. (2022) Li Z., Shen J., Gerhard O., Clarke J. P., 2022, ApJ, 925, 71
  * Lichtenberg & Lieberman (1992) Lichtenberg A., Lieberman M., 1992, Regular and Chaotic Dynamics. Springer-Verlag
  * Lin & Tremaine (1983) Lin D. N. C., Tremaine S., 1983, ApJ, 264, 364
  * Long et al. (2014) Long S., Shlosman I., Heller C., 2014, ApJ, 783, L18
  * Lucey et al. (2022) Lucey M., Pearson S., Hunt J. A. S., Hawkins K., Ness M., Petersen M. S., Price-Whelan A. M., Weinberg M. D., 2022, arXiv e-prints, p. arXiv:2206.01798
  * Lynden-Bell & Kalnajs (1972) Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1
  * Magorrian (2021) Magorrian J., 2021, MNRAS, 507, 4840
  * Martinez-Valpuesta et al. (2006) Martinez-Valpuesta I., Shlosman I., Heller C., 2006, ApJ, 637, 214
  * McGill & Binney (1990) McGill C., Binney J., 1990, MNRAS, 244, 634
  * Monari et al. (2019) Monari G., Famaey B., Siebert A., Wegg C., Gerhard O., 2019, A&A, 626, A41
  * Namouni & Guzzo (2007) Namouni F., Guzzo M., 2007, Celestial Mechanics and Dynamical Astronomy, 99, 31
  * Neishtadt (1975) Neishtadt A., 1975, Journal of Applied Mathematics and Mechanics, 39, 594
  * Ogilvie & Lubow (2006) Ogilvie G. I., Lubow S. H., 2006, MNRAS, 370, 784
  * Petersen et al. (2016) Petersen M. S., Weinberg M. D., Katz N., 2016, MNRAS, 463, 1952
  * Portail et al. (2017) Portail M., Gerhard O., Wegg C., Ness M., 2017, MNRAS, 465, 1621
  * Read et al. (2006) Read J. I., Goerdt T., Moore B., Pontzen A. P., Stadel J., Lake G., 2006, MNRAS, 373, 1451
  * Sanders et al. (2019) Sanders J. L., Smith L., Evans N. W., 2019, MNRAS, 488, 4552
  * Sellwood (2006) Sellwood J. A., 2006, ApJ, 637, 567
  * Sellwood (2016) Sellwood J. A., 2016, ApJ, 819, 92
  * Shklyar & Matsumoto (2009) Shklyar D., Matsumoto H., 2009, Surveys in Geophysics, 30, 55
  * Sormani et al. (2015) Sormani M. C., Binney J., Magorrian J., 2015, MNRAS, 454, 1818
  * Sridhar & Touma (1996) Sridhar S., Touma J., 1996, MNRAS, 279, 1263
  * Tiret & Combes (2007) Tiret O., Combes F., 2007, A&A, 464, 517
  * Tremaine & Weinberg (1984) Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729
  * Vauterin & Dejonghe (1997) Vauterin P., Dejonghe H., 1997, MNRAS, 286, 812
  * Villa-Vargas et al. (2009) Villa-Vargas J., Shlosman I., Heller C., 2009, ApJ, 707, 218
  * Villa-Vargas et al. (2010) Villa-Vargas J., Shlosman I., Heller C., 2010, ApJ, 719, 1470
  * Wang & Easther (2022) Wang Y., Easther R., 2022, Phys. Rev. D, 105, 063523
  * Wegg et al. (2015) Wegg C., Gerhard O., Portail M., 2015, MNRAS, 450, 4050
  * Weinberg (1985) Weinberg M. D., 1985, MNRAS, 213, 451
  * Weinberg (1986) Weinberg M. D., 1986, ApJ, 300, 93
  * Weinberg (1989) Weinberg M. D., 1989, MNRAS, 239, 549
  * 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
  * White (1983) White S. D. M., 1983, ApJ, 274, 53

##  Appendix A Numerical integration

In Sections 4 and 5, we solve the collisionless Boltzmann equation (52) by
numerical integration of the equations of motion averaged over the fast
angles. The averaged Hamiltonian reads

$\displaystyle\bar{H}(\theta_{\rm
s},{\bm{J}},t)=\bar{H}_{0}({\bm{J}})+\Psi({\bm{J}},t)\cos\left(\theta_{\rm
s}-\theta_{\rm s,res}\right),$ (80)

and the Hamilton’s equations are

$\displaystyle\dot{\theta}_{\rm s}$
$\displaystyle=\frac{{\partial}\bar{H}}{{\partial}J_{\rm
s}}={\bm{N}}\\!\cdot\\!\mathbf{\Omega}({\bm{J}})-N_{\varphi}\Omega_{\rm
p}(t)+\frac{{\partial}\Psi({\bm{J}},t)}{{\partial}J_{\rm
s}}\cos\left(\theta_{\rm s}-\theta_{\rm s,res}\right),$ (81)
$\displaystyle\dot{J_{\rm s}}$
$\displaystyle=-\frac{{\partial}\bar{H}}{{\partial}\theta_{\rm
s}}=\Psi({\bm{J}},t)\sin\left(\theta_{\rm s}-\theta_{\rm s,res}\right).$ (82)

Since $\bar{H}$ is non-separable, explicit symplectic integrators like leap-
frog cannot be used. Following Weinberg & Katz (2007), we employ the implicit
midpoint rule, which is a second-order symplectic integrator (Hairer et al.,
2006). In practice, we tabulate the orbital frequencies
$\mathbf{\Omega}(J_{r},L)$ and part of $\Psi$, namely the function
$W_{lm}^{N_{r}N_{\psi}}(J_{r},L,t)$ (Tremaine & Weinberg, 1984; Chiba &
Schönrich, 2022, Appendix B). The tables are then interpolated linearly.

##  Appendix B The parameters $G$ and $\gamma$

Figure 19: The resonant frequency ${\bm{N}}\\!\cdot\\!\mathbf{\Omega}$ (left),
its derivative
$G={\partial}({\bm{N}}\\!\cdot\\!\mathbf{\Omega})/{\partial}J_{\rm s}$
(middle), and the parameter $\gamma=\delta\dot{J}_{\rm res}/\dot{J}_{\rm
s,res}$ (right). The dotted lines mark contours of constant ${\bm{J}}_{\rm
f}$. The scales for $G$ and $\gamma$ are linear within $[10^{-2},10^{2}]$ and
$[10^{-1},10^{1}]$, respectively, and logarithmic elsewhere.

The left column of Fig. 19 plots the resonant frequency
${\bm{N}}\\!\cdot\\!\mathbf{\Omega}|_{\textrm{res}}=N_{\varphi}\Omega_{\rm p}$
for the major bar-halo resonances. Each contour represents the resonance curve
at a given pattern speed. The dotted lines indicate the contours of constant
fast actions ${\bm{J}}_{\rm f}$.

The middle column of Fig. 19 plots
$G\equiv{\partial}({\bm{N}}\\!\cdot\\!\mathbf{\Omega})/{\partial}J_{\rm
s}|_{\textrm{res}}$ (equation 33), which is the derivative of the left column
with respect to $J_{\rm s}$ along constant ${\bm{J}}_{\rm f}$ (dotted lines).
$G$ is mostly negative apart from very small $L$ (see also TW84, Fig. 5).
Negative $G$ implies an increase in the position of resonance in $z$-angular
momentum, $J_{\rm s,res}(=L_{z,\textrm{res}}/N_{\varphi})$, with decreasing
pattern speed $N_{\varphi}\Omega_{\rm
p}={\bm{N}}\\!\cdot\\!\mathbf{\Omega}|_{\textrm{res}}$.

The right column of Fig. 19 plots the dimensionless parameter
$\gamma=\delta\dot{J}_{\rm res}/\dot{J}_{\rm s,res}$ (equation 19), which
measures how much the trapped volume changes as the resonance moves. Apart
from regions where $G$ is close to zero, $\gamma$ is of order $10^{-1}$.