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

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 and 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 per cent of the total torque on the Milky Way’s bar.

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. Debattista & Sellwood 2000; Tiret & Combes 2007; Ghafourian, Roshan & Abbassi 2020), but also on its kinematic distribution (e.g. Athanassoula 2003; Long, Shlosman & Heller 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 disc 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 ≤ s < 1 ‘slow’ and 1 ≤ s ‘fast’.

Our current understanding on dynamical friction is limited to either of the two extreme evolutions: the fast limit (s ≫ 1) or the slow limit (s ≪ 1). Fast limit refers to a perturber rapidly changing its frequency such that the resonances sweep the orbits before any non-linearity can develop. In this limit, a time-dependent linear perturbation theory¹ 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).² 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 non-linearly and saturates due to resonant trapping. Trapping is a non-linear 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, Weinberg & Katz 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 phase-space position along each resonance.

The only study that has addressed dynamical friction in the general fast–slow regime is TW84. By solving the non-linear 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 seminumerical approach, i.e. we numerically solve the collisionless Boltzmann equation (CBE) using a Hamiltonian analytically averaged over the fast motions. This seminumerical 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 non-linear 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 non-linear torque formula derived by TW84.

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 (DF) 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 DF |$\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:

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 Ω_p(t):

where (r, |$\vartheta$|⁠, φ) 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, Shlosman & Heller 2006; Villa-Vargas, Shlosman & Heller 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_b and its corotation radius r_CR is kept constant, that is, b ≡ r_b/r_CR = const. This results in the following form (Chiba, Friske & Schönrich 2021):

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 and b = 0.28 in approximation to the hydrodynamical model of Sormani, Binney & Magorrian (2015) which is constrained by the kinematics of the gas in the inner Galactic disc.

Next we model the time-dependence of the bar’s pattern speed Ω_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 Ω_p. Thus, a better model would be to assume the following dimensionless quantity constant:

Integrating the above differential equation, we have

where Ω_p0 is the initial pattern speed. A constant η thus implies a pattern speed evolving in inverse proportion to time.³

Fig. 1 shows the pattern speed curves for a range of slowing rate η. The time is given in units of the bar’s initial rotation period T_b0 = 2|$\pi$|/Ω_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; Bovy et al. 2019; Monari et al. 2019; Sanders, Smith & Evans 2019; Asano et al. 2020; Binney 2020; Chiba & Schönrich 2021; Kawata et al. 2021; Clarke & Gerhard 2022; Gaia Collaboration 2023; Leung et al. 2022; Li et al. 2022; Lucey et al. 2023) 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}$|⁠. Fig. 1 then shows that if η = 0.001 the bar would have slowed to its present state over |$t \sim 160 T_{\rm b0} \sim 13.6 \, {\rm Gyr}$|⁠, while if η = 0.01 it merely took |$t \sim 16 T_{\rm b0} \sim 1.36 \, {\rm Gyr}$|⁠. Thus, considering the age of the Galactic bar, a reasonable range for the slowing rate would be η ∈ [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 η. Modelling of the local stellar kinematics suggests η = 0.0036 ± 0.0011 (Chiba et al. 2021).

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 semirestricted N-body simulation of Lin & Tremaine (1983), Weinberg (1985), and 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_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,

where α 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)

with current bar properties: bar semimajor 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, Gerhard & Portail 2015), resulting in |$I_{\rm rigid}\simeq 6 \times 10^{10} \, {\rm M}_\odot \, {\rm kpc}^2$|⁠. The parameter α 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 disc, not due to the 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 τ 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 α = 1, 4, 10, and 16. We choose an initial pattern speed |$\Omega _{\rm p0} = 72 \, {\rm Gyr}^{-1}$|⁠. For all α, our pattern speed model (blue) provides a reasonable fit to the simulations, where the best-fitting η is reported on the right. In agreement with Weinberg (1985), we find that a rigid bar (α = 1) slows significantly in a few bar rotation periods. Beane et al. (2022) report α = 4 from their N-body + hydrodynamical simulation. The evolution with α = 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, Combes & Semelin 2005; Villa-Vargas, Shlosman & Heller 2010; Hopkins & Quataert 2011; Athanassoula, Machado & Rodionov 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 disc (⁠|$4{-}16~{{\ \rm per\ cent}}$|⁠). We therefore explore dynamical friction over a range of slowing rate, η ∈ [0.001, 0.032], where it is to be understood that a fast bar evolution (η ∼ 0.01) is expected in gas-poor galaxies, while a slow bar evolution (η ∼ 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. TW84; Binney 2018), although we repeat them here for completeness.

To simplify the dynamics, we introduce angle-action variables |$({\boldsymbol \theta },{\boldsymbol 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({\boldsymbol J})$|⁠, and thus the equations of motion take the simple form: |$\dot{{\boldsymbol J}} = -{\partial }H_0/{\partial }{\boldsymbol \theta }= {\boldsymbol 0}$| and |$\dot{{\boldsymbol \theta }} = {\partial }H_0/{\partial }{\boldsymbol J}\equiv {\boldsymbol \Omega}({\boldsymbol 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 space.⁴ There are some freedom in the choice of the actions; we choose |${\boldsymbol 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 |${\boldsymbol \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

where |${\boldsymbol 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 _{\boldsymbol n}$| and coupled it to the exponential. |$\Phi _{\boldsymbol n}$| is given in equation (B3) of Chiba & Schönrich (2022). Near a resonance

where |${\boldsymbol N}\equiv (N_r,N_\psi ,N_\varphi)$|⁠,⁵ the equations of motion generated by H involve a term with |${\boldsymbol n} = {\boldsymbol N}$| that evolves slowly, while others (⁠|${\boldsymbol n} \ne {\boldsymbol N}$|⁠) rapidly alternate. This time-scale 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 time-scale disparity by making a canonical transformation to the fast–slow angle-action variables |$({\boldsymbol \theta }^{\prime },{\boldsymbol 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 |${\boldsymbol N}$| (TW84),

using the generating function of the second kind

We write the new Hamiltonian as

where the perturbation is expanded in the new angles |${\boldsymbol \theta }^{\prime }$| with new indices |${\boldsymbol 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

Averaging equation (13) over the fast angles |${\boldsymbol \theta }_{\rm f}\equiv (\theta _{{\rm f}_1},\theta _{{\rm f}_2})$| yields

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_φ = 2, it can be shown that terms other than k_s = ±1 vanish (Chiba & Schönrich 2022). Hence, we have

where Ψ ≡ 2|Ψ_{(0, 0, 1)}| (Chiba & Schönrich 2022, appendix B). For brevity, we have omitted the dependence on the constant fast actions |${\boldsymbol 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 Ω_p(t) and the associated growth of the bar Ψ(J_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.

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). Fig. 3 (left-hand column) shows the level curves of |$\bar{H}$| in the slow plane (θ_s, J_s) at some instance of time. We chose the corotation resonance |${\boldsymbol N}= (0,2,2)$| at |${\boldsymbol J}_{\rm f}= {\boldsymbol 0}$|⁶ as an example, but the discussion is general. The top left-hand figure shows the unperturbed Hamiltonian |$\bar{H}_0$|⁠. The trajectories are everywhere straight, and the Hamiltonian peaks exactly at the resonance J_s = J_{s, res} because the slope |${\partial }\bar{H}_0/{\partial }J_{\rm s}= {\boldsymbol N}\cdot {\boldsymbol \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_{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 θ_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_s but the phase-space area enclosed by the contours of the time-frozen |$\bar{H}$|

where we have divided the area by 2|$\pi$| to make it an action variable. The angles θ_ℓ conjugate to J_ℓ describe the phase of libration (or circulation). Adiabatic theory assumes that |$\bar{H}$| evolves slowly enough that actions J_ℓ 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_ℓ change discontinuously. The fate of an orbit crossing the separatrix is known to be extremely sensitive to its angle θ_ℓ and thus cannot be practically predicted. One can, however, calculate the probability of transition averaged over θ_ℓ (Neishtadt 1975; Henrard 1982), and this has been used to model the angle-averaged (coarse-grained) DF 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_ℓ to be sufficiently shorter than the time-scale at which the separatrices move by the width of the trapped phase-space. Let v_± be the speed of the separatrices above and below the resonance. To order of magnitude, v_± 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 δJ_res changes:

where the dimensionless parameter

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

As we shall see later (equations 34 and 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

For the problem at hand, γ is typically of the order of 10⁻¹ over the majority of phase space (Appendix  B), so the adiabatic condition is largely dictated by the ratio η/A. Models of the Galactic bar suggest A ∼ 10⁻² while η ∼ 10⁻³ (Section 2), so η/A ∼ 10⁻¹. 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

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

A suitable generating function for this transformation is

The new Hamiltonian is then

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

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 Δ.

The right-hand 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-hand 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-hand figure) no longer appear straight but move towards low Δ 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-hand panel) 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 θ = 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.

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

where the integral runs along constant time-frozen |$\widehat{H}$|⁠. Fig. 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 |$\boldsymbol {w}_i$| at each time step, we calculated the libration action by integrating a new orbit starting from |$\boldsymbol {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_ℓ 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

The adiabatic condition is then

which, to leading order, requires

As compared to the original criterion (21), the new criterion depends directly on the size of the parameter γ (equation 19). This implies that if the shape of the separatrix in the comoving frame varies little (γ ≪ 1), the adiabaticity does not necessary require a slowly moving resonance (η/A ≪ 1) as in the classical adiabatic assumption. In the example given above (Fig. 4), γ ∼ 0.1 so γη/A ∼ 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 Δ up to second order for the unperturbed term and to zeroth order for the perturbed term

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

where

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:

When the resonance is static (⁠|$\dot{J}_{\rm s, res}= 0$|⁠), ω₀ describes the libration frequency of orbits with infinitesimal libration amplitude and thus provides the time-scale of libration, while J₀ 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}\ne 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

where

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

Using G (equation 33), we obtain

Substituting this into equation (36), we have

which is equivalent to equation (89) of TW84. This expression indicates that, as opposed to the dimensionless slowing rate η (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 η is a global parameter, s is a local parameter that depends on the resonance |${\boldsymbol N}$| as well as on the position on each resonance set by the fast actions |${\boldsymbol J}_{\rm f}$|⁠.

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

where V is the dimensionless effective potential

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

Fig. 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 protoplanetary disc (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 Δ = 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 Δ and plugging it into equation (27)

where θ^± are the roots of the equation |$\widehat{H}+ F (\cos \theta + s \theta) = 0$|⁠. Fig. 6 illustrates the definition of θ^±.

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

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 time.⁷

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

where the function C(s) is the resonant-volume factor⁸

with domain 0 ≤ s ≤ 1 and range 0 ≤ C ≤ 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

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 θ_res as the angle at which the potential is locally minimum:

Fig. 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 θ_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 measures 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, φ_L − φ_b, which yields the speed of the resonance as s = sin [2(φ_L − φ_b)]. We will pursue this possibility in a separate paper.

4 EVOLUTION OF THE DF

Using the averaged Hamiltonian derived in the previous section, we now study the evolution of the 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({\boldsymbol \theta },{\boldsymbol J},t)$| is described by the CBE

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

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 |${\boldsymbol 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 |${\boldsymbol \theta }_{\rm f}$| at all times, since the averaged trajectories are independent of |${\boldsymbol \theta }_{\rm f}$|⁠. The problem then reduces to solving the trajectories in the 2D slow angle-action space

where, as with the Hamiltonian, we have omitted references to |${\boldsymbol 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.

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 |${\boldsymbol \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.

Fig. 8 shows the time-evolution of the phase-space density in the slow angle-action plane (θ_s, J_s). As before, we choose the corotation resonance |${\boldsymbol 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 η = 0.002. This results in the speed parameter increasing from s ≃ 0.134–0.184, which falls under the slow regime (s < 1).

The left-hand 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_s, where the phase-space density is high. As the bar slows, the resonance moves towards larger J_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-hand 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 |${\boldsymbol \theta }_{\rm f}$| effectively removes all Fourier series of perturbations other than the target resonance.

4.1 Friction and feedback in the slow non-linear regime

Fig. 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_s, and the negative gradient in f₀ 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 and 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 stripes⁹ 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}$|⁠.

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. Fig. 9 shows the variation of the phase-space distribution with slowing rate η. We chose η such that they are approximately evenly spaced in logarithm with width Δlog₁₀η ≈ 0.5. The top row plots the distribution after the bar slowed from Ω_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 η. For η = 0.001 and 0.0032, we observe new trapped layers because the trapped volume increases with time, as shown in the bottom row. For η = 0.01, however, the volume reduces, causing trapped orbits to escape the resonance and create a dense trailing wake.

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

where the prime denotes a time derivative. Since our bar potential has a fixed amplitude and elongates along with the expanding resonance, Ψ′/Ψ is minor. In contrast, −G′/G is large and positive since G rapidly approaches zero as the resonance moves towards larger J_s (Appendix  B). The third term, on the right-hand side, of equation (53) describes the change in the resonant-volume factor C (equation 47) and is typically negative since dC/ds < 0 while s′ > 0 (Fig. 9, middle row). For small s, and thus C ∼ 1, the second term dominates, resulting in an increasing trapped volume. However, in the transition regime (third left-hand 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-hand 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):

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

where |$\Omega _s({\boldsymbol J},t) \equiv {\boldsymbol N}\cdot {\boldsymbol \Omega}({\boldsymbol J})-n_\varphi \Omega _{\rm p}(t)$|⁠. We transform |${\boldsymbol \theta }$| to the fast–slow angles and average f₁ over the fast angles

Since |$\Phi _{{\boldsymbol N}}$| is non-zero only when n_φ = ±m = ±2 (Chiba & Schönrich 2022, appendix B), we have for resonances with N_φ = 2,

where δ_{n, N} is the Kronecker delta. Using the reality condition |$\Phi _{-{\boldsymbol n}} = \Phi _{{\boldsymbol n}}^{*}$| and the relation |${\boldsymbol N}\!\cdot \left({\partial }f_0/{\partial }{\boldsymbol J}\right) = {\partial }f_0/{\partial }J_{\rm s}$|⁠, we finally obtain

Fig. 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_lib since the non-linear response is precluded (TW84). Hence, while linear theory works in the slow regime (s < 1) only for a short period (t ≲ t_lib), it remains valid in the fast regime (s > 1) over an extended period (t ≳ t_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 DF over a Cartesian grid (96 × 96 × 96) in 3D velocity space:

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_max = Φ₀(r_max), where we set |$r_{\rm max} = 10 r_{\rm s}= 200 \, {\rm kpc}$|⁠. Calculations with r_max > 10r_s made little visible difference.

Fig. 11 shows the density response at the Galactic plane (z = 0) caused by the four strongest resonances: the corotation resonance |${\boldsymbol N}=(0,2,2)$|⁠, the outer Lindblad resonance |${\boldsymbol N}=(1,2,2)$|⁠, and the direct radial resonances |${\boldsymbol N}=(1,0,2)$| and |${\boldsymbol N}=(2,0,2)$|⁠. The bar lies on the x-axis and is rotating anticlockwise. As before, we slow the bar from Ω_p = 72–|$36\, {\rm Gyr}^{-1}$|⁠.

The density response varies significantly with slowing rate η as the evolutionary regime makes a transition from slow to fast. For η as small as 0.001 (Fig. 11 left-hand 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 θ_res shifts anticlockwise when the resonance is moving (see discussion surrounding Figs 5–7). As η is increased from the 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 |${\boldsymbol 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/non-linear. 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.

Fig. 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 η (left-hand 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 non-linear 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_fb is a coefficient that has the dimension of a moment of inertia, and denote the usual dynamical friction as τ_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

which implies that dynamical feedback changes the moment of inertia of the perturber (Weinberg 1985). If I_fb > 0, dynamical feedback helps the bar spin down (positive feedback); otherwise, trapping stabilizes the bar’s pattern speed (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 feedback from trapped and untrapped orbits are in fact closely related through the CBE, 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 DF any transient perturbations that phase mixes over time. This is effected by averaging the DF over the angle of libration

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

and integrate over the slow phase-space dθ_sdJ_s = dθ_ℓdJ_ℓ

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

where we have plugged in equation (38). We remind the reader that G is |${\partial }({\boldsymbol N}\cdot {\boldsymbol \Omega})/{\partial }J_{\rm s}$| evaluated at the resonance J_s = J_{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 |${\boldsymbol 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 |${\boldsymbol 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

where the function

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

using the identity ∑_i f(x_i)/|g′(x_i)| = ∫dx f(x)δ(g(x)) (e.g. Kanwal 1998), where δ 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

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

This equation describes the dynamical feedback from all the trapped orbits dragged by the moving resonances: the function B_ℓ measures the mass of trapped orbits, and the 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_ℓ ≥ 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.

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 DF 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 DF are of second order in the strength of perturbation A (e.g. Fouvry, Pichon & Prunet 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₀. 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₀ with |$-\dot{L}_z$| and integrating it over the slow phase-space d²S = dθ_sdJ_s, while excluding the region of libration S_ℓ from the domain of integration:

The first term in the second line is zero since the torque is periodic in the slow angle θ_s. We next Taylor expand f₀ in J_s to first order around the resonance J_{s, res}

where Δ = J_s − J_{s, res}. We similarly expand the Hamiltonian (equation 32) which makes the second term vanish since the domain of integration S_ℓ becomes symmetric about Δ = 0, and Δ is odd while |$\dot{L}_z$| is even in Δ. 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

where the function

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

Substituting B₀ to equation (72), and using equation (44), we obtain

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 and 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

This equation states that the strength of dynamical feedback is determined by the difference between the mass of the trapped orbits B_ℓ and the mass of the background orbits B₀ 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 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₀ − B_ℓ. Thus,

Our stability condition differs from that of TW84 (equation 104), which does not include feedback from trapped orbits, so B₀ − B_ℓ is always positive. Since we expect B₀ − B_ℓ < 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 DF 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 η/A → 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 ηγ/A → 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_{s, res} at the time of capture t_cap:

Fig. 13 compares this model (left-hand panel) with the original distribution (middle panel) for η = 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_{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-hand panel) exhibits the phase spiral with little offset, confirming that the model is reasonable.

Using this model for |$\tilde{f}_{\ell }$|⁠, we compute the dynamical feedback (75). Fig. 14 shows the dynamical feedback by the corotation resonance as a function of time for various slowing rate η, ranging from 0.001 (blue) to 0.01 (red). As before, we fix the bar’s initial and final pattern speeds (Ω_p = 72–|$36 \, {\rm Gyr}^{-1}$|⁠), so the time the bar takes to reach its final pattern speed (black dot) gets shorter with increasing η. The feedback torque is initially zero since B_ℓ = B₀, 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 η ∼ 0.0025. To understand this, we analyse the feedback |$\hat{\tau }_{\rm fb}$| from orbits near a single point of resonance at a specific |${\boldsymbol J}_{\rm f}$|⁠. From equation (65), we have

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_ℓ as |$J_{\ell, {\rm sep}} f_{\ell 0} = (8 J_0 / \pi) C(s) f_{\ell 0}$|⁠, where f_ℓ0 is the averaged density of the librating region, and write

The equation indicates that feedback from each point on the resonance scales as sC(s). Fig. 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 ∼ 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 η (Fig. 14) can be understood as the integration of this curve over the entire resonant surface.

Fig. 16 compares dynamical feedback with the total torque by the corotation resonance. The latter is computed by sampling particles from f₀ 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.

Fig. 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.

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 η of the order of 10⁻³ (black and blue), around 30 per cent of the torque is composed of dynamical feedback. For large η (≳ 10⁻²) 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 friction¹⁰, 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 |${\boldsymbol 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 \lt 0$| for all v_x > 0. Then, the phase-space density of the untrapped region at the current velocity of the perturber v_px, is larger than the phase-space density of the trapped region, which consists of orbits captured at v_x > v_px. Hence, the momentum gained by the untrapped orbits exceeds the momentum loss of the trapped orbits. Consequently, particles exert a net negative force (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 DF 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 ≫ 1). The mechanism of dynamical friction can be summarized as follows:

• 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 DF, but undergoes damped oscillations as the distribution phase mixes (Banik & van den Bosch 2022; Chiba & Schönrich 2022).

• 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 DF, 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 ∼ 0.4 and decays towards s = 1.

• 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.

• In the fast limit (s ≫ 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 DF, 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 per cent of the total torque on the Galactic bar, if the bar has been evolving with η ∼ 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 DF (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 2 (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 non-linear 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\R1\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.

Footnotes

References

APPENDIX A: NUMERICAL INTEGRATION

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

and the Hamilton’s equations are

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 Ω(J_r, L) and part of Ψ, namely the function |$W_{lm}^{N_rN_\psi }(J_r,L,t)$| (TW84; Chiba & Schönrich 2022, appendix B). The tables are then interpolated linearly.

APPENDIX B: THE PARAMETERS G AND γ

The left-hand column of Fig. B1 plots the resonant frequency |${\boldsymbol N}\!\cdot\!{\boldsymbol \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 |${\boldsymbol J}_{\rm f}$|⁠.

The middle column of Fig. B1 plots |$G \equiv {\partial }({\boldsymbol N} \!\cdot\! {\boldsymbol \Omega})/{\partial }J_{\rm s}|_\textrm {res}$| (equation 33), which is the derivative of the left-hand column with respect to J_s along constant |${\boldsymbol 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_{s, res}(= L_{z, res}/N_φ), with decreasing pattern speed |$N_\varphi \Omega _{\rm p}= {\boldsymbol N}\!\cdot \! {\boldsymbol \Omega}|_\textrm {res}$|⁠.

The right-hand column of Fig. B1 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, γ is of the order of 10⁻¹.