Dynamical Friction, Buoyancy, and Core-stalling. I. A Nonperturbative Orbit-based Analysis

In the perturbative framework of TW84 and KS18, dynamical friction arises from the LBK torque, which only has a nonzero contribution from pure resonances, i.e., orbits that obey a commensurability condition between the (circular) frequency of the perturber, Ω_P, and the frequencies of the field particles in the unperturbed potential. Even the more general, self-consistent torque introduced by BB21 is formulated in terms of these frequencies.

In the nonperturbative framework adopted in this paper, in which we consider fully perturbed orbits ² in the galaxy+perturber potential to arbitrary order, the frequencies of the individual field particles vary with time due to energy and angular momentum exchanges with the perturber; the original actions of the unperturbed galaxy are no longer conserved, and neither are the frequencies associated with the corresponding angles (Tremaine & Weinberg 1984; Fouvry & Bar-Or 2018). Hence, a field particle will not satisfy a commensurability condition throughout its orbital evolution, but rather will find itself "trapped," librating around resonance(s) with the perturber. In fact, this is what happens when the field particle along the horseshoe orbit in Figure 1 moves from B to C and from D to E; its azimuthal frequency, Ω_ϕ , is swept back and forth through a near-corotation resonance with the circular frequency of the perturber, Ω_P. This same principle also underlies the physics of radial migration in disks due to interactions with transient spirals (e.g., Carlberg & Sellwood 1985; Sellwood & Binney 2002; Daniel & Wyse 2015). Dynamical friction arises from an imbalance between the number of field particles that "sweep up" versus "sweep down" in frequency space, and this imbalance itself arises from gradients in the distribution function.

The geometry of our dynamical model is illustrated in Figure 3. It depicts the galaxy (large, shaded circle), the perturber (solid black dot), and the COM in the corotating (x,y) frame that we will adopt throughout. For convenience, we define the following mass ratios: q ≡ M_P/M_G is the mass ratio of the infalling perturber and the host galaxy, while q_enc(R) ≡ M_P/M_G(R) is the mass ratio of the perturber and the galaxy enclosed within R. The distances between the COM and the galactic center and between the COM and the perturber are given by q_G R and q_P R, respectively, where

$$\begin{eqnarray}\begin{array}{rcl}{q}_{{\rm{G}}} & = & \displaystyle \frac{{M}_{{\rm{P}}}}{{M}_{{\rm{P}}}+{M}_{{\rm{G}}}(R)}=\displaystyle \frac{{q}_{\mathrm{enc}}(R)}{1+{q}_{\mathrm{enc}}(R)},\\ {q}_{{\rm{P}}} & = & \displaystyle \frac{{M}_{{\rm{G}}}(R)}{{M}_{{\rm{P}}}+{M}_{{\rm{G}}}(R)}=\displaystyle \frac{1}{1+{q}_{\mathrm{enc}}(R)}.\end{array}\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

Throughout this paper, we adopt dimensionless units to describe our dynamical system. All length scales are expressed in units of r_s, the scale radius of the galaxy, masses are expressed in units of the mass of the galaxy, M_G, and velocities are expressed in units of $\sigma ={({{GM}}_{{\rm{G}}}/{r}_{{\rm{s}}})}^{1/2}$. The corresponding, characteristic timescale is r_s/σ.

For convenience, we consider the perturber to be a point mass, but we emphasize that the analysis that follows can be easily extended to accommodate any other (spherically symmetric) perturber potential. In our dimensionless units, we then have that the perturber potential

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{P}}}=-q/{r}_{{\rm{P}}}.\end{eqnarray} \tag{ 2 }$$

Throughout this work, we adopt q = 0.004 (i.e., the mass of the perturber is only 0.4% of that of the galaxy). Unlike the perturbative treatments in TW84 and KS18, though, which require q to be small, our analysis is also valid for more massive perturbers.

In order to contrast dynamical friction in cored and cuspy density profiles, we consider two different density profiles for the galaxy: a Plummer sphere, which has a central constant-density core with central logarithmic density gradient, $\gamma \equiv {\mathrm{lim}}_{r\to 0}d\mathrm{log}\rho /d\mathrm{log}r=0$ (Plummer 1911), and a Hernquist sphere, which has a central γ = −1 cusp (Hernquist 1990). Both have the advantage that the density and potential are given by simple, analytical expressions. For the Plummer sphere, the density and potential (in our dimensionless units) are given by

$$\begin{eqnarray}&&{\rho }_{{\rm{G}}}(r)=\displaystyle \frac{3}{4\pi }\displaystyle \frac{1}{{\left(1+{r}^{2}\right)}^{5/2}},\,\,\,\,\,\,{{\rm{\Phi }}}_{{\rm{G}}}(r)=-\displaystyle \frac{1}{\sqrt{1+{r}^{2}}},\end{eqnarray} \tag{ 3 }$$

while for the Hernquist sphere we have that

$$\begin{eqnarray}&&{\rho }_{{\rm{G}}}(r)=\displaystyle \frac{1}{2\pi }\,\displaystyle \frac{1}{r\,{\left(1+r\right)}^{3}},\,\,\,\,\,{{\rm{\Phi }}}_{{\rm{G}}}(r)=-\displaystyle \frac{1}{1+r}.\end{eqnarray} \tag{ 4 }$$

Figure 4 plots these density profiles (left-hand panel) and corresponding logarithmic density gradients, $d\mathrm{log}\rho /d\mathrm{log}r$ (right-hand panel), as functions of radius. The magenta and black vertical dashed lines indicate R = 0.2 and 0.5, respectively. These are the orbital radii of the perturber considered in this paper. As we demonstrate below, in the case of the Plummer host, these radii bracket the bifurcation radius, R_bif (≈0.39 for our fiducial case), at which the orbital makeup of the Plummer sphere undergoes a drastic change due to a bifurcation of some of the Lagrange points, which in turn impacts the nature (retarding versus enhancing) of the torque on the perturber. In the case of the Hernquist sphere, no such bifurcation occurs.

Zoom In Zoom Out Reset image size

Throughout, we assume that the galaxies have isotropic velocity distributions, such that their distribution functions are ergodic (i.e., depend only on energy). In the case of the Plummer sphere, we have

$$\begin{eqnarray}&&{f}_{0}(\varepsilon )=\displaystyle \frac{3}{7{\pi }^{3}}{\left(2\varepsilon \right)}^{7/2},\end{eqnarray} \tag{ 5 }$$

while for the Hernquist sphere,

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0}(\varepsilon ) & = & \displaystyle \frac{1}{8\sqrt{2}{\pi }^{3}}\\ & & \times \,\displaystyle \frac{3{\sin }^{-1}\sqrt{\varepsilon }+\sqrt{\varepsilon (1-\varepsilon )}(1-2\varepsilon )(8{\varepsilon }^{2}-8\varepsilon -3)}{{\left(1-\varepsilon \right)}^{5/2}}.\end{array}\end{eqnarray} \tag{ 6 }$$

Here, ε = −E_0G (E_0G is the unperturbed galactocentric energy), and the subscript "0" indicates that these distribution functions correspond to the unperturbed galaxies. Both distribution functions have been normalized such that

$$\begin{eqnarray}&&{\rho }_{{\rm{G}}}(r)=4\pi {\int }_{0}^{{{\rm{\Psi }}}_{{\rm{G}}}}\sqrt{2({{\rm{\Psi }}}_{{\rm{G}}}-\varepsilon )}\,{f}_{0}(\varepsilon )\,d\varepsilon ,\end{eqnarray} \tag{ 7 }$$

with Ψ_G = − Φ_G.

Since the gravitational potential and hence the Hamiltonian in the restricted three-body problem is time-variable, energy is not a conserved quantity. And due to the lack of spherical symmetry, neither is angular momentum. However, as is well-known (see, e.g., Binney & Tremaine 2008), the Jacobi integral,

$$\begin{eqnarray}&&{E}_{{\rm{J}}}=E-{{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\cdot {\boldsymbol{L}}=\displaystyle \frac{1}{2}{\dot{{\boldsymbol{r}}}}^{2}+{{\rm{\Phi }}}_{\mathrm{eff}}\left({\boldsymbol{r}}\right),\end{eqnarray} \tag{ 8 }$$

is a conserved quantity. Here, r is the position vector of the field particle with respect to the COM (see Figure 3), and Ω _P = (0, 0, Ω_P) with

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{P}}}=\sqrt{\displaystyle \frac{G\,[{M}_{{\rm{G}}}(R)+{M}_{{\rm{P}}}]}{{R}^{3}}},\end{eqnarray} \tag{ 9 }$$

the angular frequency of the perturber with respect to the COM, which, in our dimensionless units, is given by

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{P}}}={\left(\displaystyle \frac{1}{R}{\left.\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{G}}}}{\partial {r}_{{\rm{G}}}}\right|}_{{r}_{{\rm{G}}}=R}+\displaystyle \frac{q}{{R}^{3}}\right)}^{1/2}.\end{eqnarray} \tag{ 10 }$$

E and L are, respectively, the perturbed energy and angular momentum (per unit mass) of the field particle in the nonrotating, inertial frame, given by

$$\begin{eqnarray}&&E={E}_{0}+{{\rm{\Phi }}}_{{\rm{P}}}=\displaystyle \frac{1}{2}| \dot{{\boldsymbol{r}}}+{{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times {\boldsymbol{r}}{| }^{2}+{{\rm{\Phi }}}_{{\rm{G}}}\left({\boldsymbol{r}}\right)+{{\rm{\Phi }}}_{{\rm{P}}}\left({\boldsymbol{r}}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\boldsymbol{L}}={\boldsymbol{r}}\times \left(\dot{{\bf{r}}}+{{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times {\boldsymbol{r}}\right).\end{eqnarray} \tag{ 12 }$$

Here, E₀ is the unperturbed energy, i.e., the part of the Hamiltonian without the perturber potential, and Φ_G and Φ_P are the gravitational potentials due to the galaxy and the perturber, respectively. The effective potential in Equation (8) is defined as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{eff}}\left({\boldsymbol{r}}\right)={{\rm{\Phi }}}_{{\rm{G}}}({r}_{{\rm{G}}})+{{\rm{\Phi }}}_{{\rm{P}}}({r}_{{\rm{P}}})-\displaystyle \frac{1}{2}| {{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times {\boldsymbol{r}}{| }^{2},\end{eqnarray} \tag{ 13 }$$

where r_G and r_P are the distances to the field particle from the galactic center and the perturber, respectively, and are given by

$$\begin{eqnarray}\begin{array}{rcl}{r}_{{\rm{G}}}^{2} & = & {r}^{2}+{q}_{{\rm{G}}}^{2}{R}^{2}+2\,{q}_{{\rm{G}}}R\,r\cos \phi ,\\ {r}_{{\rm{P}}}^{2} & = & {r}^{2}+{q}_{{\rm{P}}}^{2}{R}^{2}-2\,{q}_{{\rm{P}}}R\,r\cos \phi .\end{array}\end{eqnarray} \tag{ 14 }$$

Here, r = ∣ r ∣, ϕ is the counterclockwise angle between r and the line connecting the COM and the perturber positioned along the positive x-axis (see Figure 3), and q_G and q_P are the mass ratios given by Equation (1). The third term in Equation (13) is the potential due to the centrifugal force. Plugging in the expression for Ω_P, and using the fact that ∂Φ_G/∂r =GM_G(r)/r² and q_enc(R) = M_P/M_G(R), the effective potential reduces to

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{eff}}({\boldsymbol{r}})={{\rm{\Phi }}}_{{\rm{G}}}({r}_{{\rm{G}}})-\displaystyle \frac{q}{{r}_{{\rm{P}}}}-\displaystyle \frac{1+{q}_{\mathrm{enc}}(R)}{2}\,\displaystyle \frac{{r}^{2}}{R}{\left.\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{G}}}}{\partial {r}_{{\rm{G}}}}\right|}_{{r}_{{\rm{G}}}=R}.\end{eqnarray} \tag{ 15 }$$

As already mentioned above, in the perturbed potential, energy and angular momentum are no longer constants of motion. Instead, the only conserved quantity in the restricted three-body case considered here is the Jacobi energy, E_J. A field particle therefore gains and loses energy and angular momentum (which is exchanged with the perturber) as it traverses its orbit. In order to compute the rates at which the energy and angular momentum of a field particle change as function of time, we integrate its orbit using the equation of motion in the corotating frame (Binney & Tremaine 1987), which is given by

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}=-{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{G}}}-{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{P}}}-2\left({{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times \dot{{\boldsymbol{r}}}\right)-{{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times \left({{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times {\boldsymbol{r}}\right).\end{eqnarray} \tag{ 16 }$$

Here, the first and second terms on the RHS denote the gravitational accelerations due to the galaxy and the perturber, respectively, while the third and the fourth terms correspond to accelerations due to the Coriolis and centrifugal forces, respectively. In cylindrical coordinates, the above reduces to the following radial and azimuthal equations of motion:

$$\begin{eqnarray}\begin{array}{rcl}\ddot{r}-r\,{\dot{\phi }}^{2} & = & -\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{G}}}}{\partial r}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial r}+2\,{{\rm{\Omega }}}_{{\rm{P}}}\,r\,\dot{\phi }+{{\rm{\Omega }}}_{{\rm{P}}}^{2}\,r,\\ r\,\ddot{\phi }+2\dot{r}\,\dot{\phi } & = & -\displaystyle \frac{1}{r}\,\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{G}}}}{\partial \phi }-\displaystyle \frac{1}{r}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial \phi }-2\,{{\rm{\Omega }}}_{{\rm{P}}}\,\dot{r}.\end{array}\end{eqnarray} \tag{ 17 }$$

The latter can be combined with Equations (14) to yield an expression for the torque,

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal T } & = & \displaystyle \frac{{dL}}{{dt}}=-\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{G}}}}{\partial \phi }-\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial \phi }\\ & = & \displaystyle \frac{{rR}\,\sin \,\phi }{1+{q}_{\mathrm{enc}}(R)}\left[\displaystyle \frac{{q}_{\mathrm{enc}}(R)}{{r}_{{\rm{G}}}}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{G}}}}{\partial {r}_{{\rm{G}}}}-\displaystyle \frac{1}{{r}_{{\rm{P}}}}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial {r}_{{\rm{P}}}}\right],\end{array}\end{eqnarray} \tag{ 18 }$$

where $L={r}^{2}(\dot{\phi }+{{\rm{\Omega }}}_{{\rm{P}}})$ is the total angular momentum of the field particle in the inertial frame. Equation (18) is an expression for the combined torque, exerted by both the perturber and the galaxy on the field particle. For a slowly evolving circular orbit of the perturber, i.e., nearly constant Ω_P, as considered in this paper, E_J = E − Ω _P · L is a conserved quantity. Hence, the corresponding rate of energy change of the field particle is simply given by

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}={{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\cdot \displaystyle \frac{d{\boldsymbol{L}}}{{dt}}.\end{eqnarray} \tag{ 19 }$$

Because of this equality, throughout this paper we will talk about ΔE and ΔL interchangeably. Note that, depending on the sign of the torque ${ \mathcal T }={dL}/{dt}$, the perturber can either lose (dynamical friction) or gain energy (dynamical buoyancy). Also note that dynamical friction or buoyancy results in a nonzero time-derivative of Ω _P , which, following TW84 and KS18, has been ignored in the above equations. Since we are mainly interested in examining dynamical friction near the core-stalling radius, where ∣d Ω _P /dt∣ vanishes, this is justified. In fact, it is justified as long as the timescale for dynamical friction is sufficiently long, i.e., we are in what TW84 refer to as the "slow" regime.

Throughout this paper, all orbit integrations are performed using an exactly Hamiltonian-conserving algorithm proposed by Kotovych & Bowman (2002) for simulating general N − body systems. It ensures that the Jacobi Hamiltonian is conserved up to machine precision for all the orbits we have integrated.

To get a better understanding of dynamical friction, it is instructive to study the different kinds of stellar orbits that arise in presence of the perturber. Using Equation (16), we numerically integrate stellar orbits in the corotating frame under the combined gravitational potential of the perturber plus galaxy. Along each orbit, we then register the time evolution of the orbital energy and angular momentum. We emphasize that, in doing so, the perturber is fully accounted for (i.e., is not treated as a small perturbation). For the sake of simplicity, though, we restrict ourselves to 2D, and only study the dynamics in the orbital plane of the perturber.

One can gain valuable insight regarding the orbital families by examining the system's equipotential contours, which can be parameterized by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{eff}}\left({\boldsymbol{r}}\right)={E}_{{\rm{J}}}.\end{eqnarray} \tag{ 20 }$$

These contours are zero-velocity curves (ZVCs) since they map out the locations in the corotating frame where the field particles of a given Jacobi energy E_J have zero velocity (in the corotating frame). Therefore, field particles along an orbit can only occasionally touch its ZVC and can only access regions on the side of its ZVC where its Jacobi energy E_J > Φ_eff( r ).

Of particular relevance are the fixed points, also known as the Lagrange points, where the effective force in the corotating frame vanishes. These are given by the roots of

$$\begin{eqnarray}&&{\rm{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}=0.\end{eqnarray} \tag{ 21 }$$

As we discuss below, and in detail in Paper II, the number of Lagrange points depends on the inner logarithmic slope γ of the galaxy density profile and the galactocentric distance R of the perturber.

All orbits in the restricted three-body problem have some sense of circulation, either around the galactic center, around the perturber, around the COM, or around a specific Lagrange point. ³ We can discriminate between these different cases by considering the circular part of their Jacobi energy, E_Jc, evaluated in the neighborhood of a center of perturbation (COP) (either the location of the perturber or a stable Lagrange point such as L4, L5, or L0),

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{Jc}} & = & {E}_{{\rm{J}}}-\displaystyle \frac{{\left({\rm{\Delta }}{{\rm{\Omega }}}_{0}\right)}^{2}}{4\left|{c}_{0}\right|}\\ & & -\,\left({\kappa }_{0}+\displaystyle \frac{{b}_{0}}{\left|{c}_{0}\right|}{\rm{\Delta }}{{\rm{\Omega }}}_{0}\right){J}_{r}-\left({a}_{0}+\displaystyle \frac{{b}_{0}^{2}}{\left|{c}_{0}\right|}\right){J}_{r}^{2}.\end{array}\end{eqnarray} \tag{ 22 }$$

Here, J_r is the radial action, ΔΩ₀ = Ω₀ − Ω_P, Ω₀ and κ₀ are the azimuthal and radial epicyclic frequencies, respectively, and a₀, b₀, and c₀ are constants that depend on the galaxy potential, evaluated at the COP (see Appendix A for details). The family of an orbit is dictated by the values of E_Jc computed in the neighborhood of L4/L5 (${E}_{\mathrm{Jc}}^{\left(4\right)}$), L0 (${E}_{\mathrm{Jc}}^{\left(0\right)}$) and the perturber (${E}_{\mathrm{Jc}}^{{\rm{P}}}$), respectively, relative to the values of the effective potential, Φ_eff, at the various Lagrange points and the location of the perturber. In what follows, we use ${E}_{{\rm{J}}}^{\left(k\right)}$, with k = 0, 1, ..., 5 to (approximately) indicate the value of Φ_eff at the k^th Lagrange point (e.g., ${E}_{{\rm{J}}}^{(3)}$ indicates the Φ_eff value corresponding to the equipotential/zero-velocity contour that passes through L3), and ${E}_{{\rm{J}}}^{{\rm{P}}}$ to indicate the value of the effective potential at the location of the perturber (see Appendix A for details). For nearly circular orbits with J_r ≈ 0, E_Jc ≈ E_J and the orbital families are roughly dictated by the equipotential contours.

We start our census of the orbital families by considering a Plummer galaxy with a massive perturber (as always, assumed to be a point mass with q = 0.004) orbiting at R = 0.5, which is outside of the bifurcation radius (see Figure 4). The equipotential contours (Φ_eff = E_J) in this case are as depicted in Figure 5(a). The system has six Lagrange points (L0, L1, L2, L3, L4, and L5) as indicated. Of these, L0 (which coincides with the galactic center), L1, L2, and L3 all lie along the x-axis, while L4 and L5 are located symmetrically on both sides of it, each forming an equilateral triangle with L0 and the perturber. As discussed in detail in Paper II, the Lagrange points L0, L4, and L5 are stable fixed points (centers), while L1, L2 and L3 are unstable fixed points (saddles).

Zoom In Zoom Out Reset image size

We identify different orbital families based on the circular part of the Jacobi energy, E_Jc, as specified in Table 1; for nearly circular orbits, this amounts to considering only the Jacobi energy since they lie close to their ZVCs. All such near-circular orbits with ZVCs inside the same shaded region in Figure 5(a) have similar morphology and are taken to belong to the same orbital family. These families are separated by the ZVCs passing through the saddle points, known as separatrices. Note, though, that since J_r can vary along an orbit, certain orbits (especially those with higher eccentricities in the inertial frame) can transition between different orbital families by undergoing separatrix crossings. We shall address these special kinds of orbits separately toward the end of this section, and shall proceed with the delineation of orbital families using E_Jc for now.

Table 1. Different Orbital Families in the Corotating Frame of Our Restricted Three-body Framework

Orbit type	${E}_{\mathrm{Jc}}^{\left(4\right)}$	${E}_{\mathrm{Jc}}^{\left(0\right)}$	${E}_{\mathrm{Jc}}^{{\rm{P}}}$	L	COC	Friction (F)/
						Buoyancy (B)/
						Negligible (N)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
Horseshoe (γ = 0)	${E}_{\mathrm{Jc}}^{\left(4\right)}\lt {E}_{{\rm{J}}}^{(3)}$	${E}_{\mathrm{Jc}}^{\left(0\right)}\gt \max \left[{E}_{{\rm{J}}}^{(1)},{E}_{{\rm{J}}}^{(2)}\right]$	⋯	⋯	L3	F
Horseshoe (γ < 0)	${E}_{\mathrm{Jc}}^{\left(4\right)}\lt {E}_{{\rm{J}}}^{(3)}$	⋯	${E}_{\mathrm{Jc}}^{{\rm{P}}}\gt \max \left[{E}_{{\rm{J}}}^{(1)},{E}_{{\rm{J}}}^{(2)}\right]$	⋯	L3	F
Pac-Man (γ = 0)	⋯	${E}_{\mathrm{Jc}}^{\left(0\right)}\lt {E}_{{\rm{J}}}^{(1)}$	${E}_{\mathrm{Jc}}^{{\rm{P}}}\gt {E}_{{\rm{J}}}^{(2)}$	L(1) < L < L(2)	L0	F/B
Pac-Man (γ < 0)	⋯	⋯	${E}_{{\rm{J}}}^{(2)}\lt {E}_{\mathrm{Jc}}^{{\rm{P}}}\lt {E}_{{\rm{J}}}^{(1)}$	L(1) < L < L(2)	L0	F/B
Tadpole (R > Rbif)	${E}_{{\rm{J}}}^{(3)}\lt {E}_{\mathrm{Jc}}^{\left(4\right)}\lt {E}_{{\rm{J}}}^{\left(4\right)}$	⋯	⋯	⋯	L4/L5	F/B
Tadpole (R ≤ Rbif)	${E}_{{\rm{J}}}^{\left(0\right)}\lt {E}_{\mathrm{Jc}}^{\left(4\right)}\lt {E}_{{\rm{J}}}^{\left(4\right)}$	⋯	⋯	⋯	L4/L5	F/B
Center-phylic (γ = 0)	⋯	${E}_{{\rm{J}}}^{\left(0\right)}\lt {E}_{\mathrm{Jc}}^{\left(0\right)}\lt {E}_{{\rm{J}}}^{(1)}$	⋯	L < L(1)	L0	N
Center-phylic (γ < 0)	⋯	⋯	${E}_{\mathrm{Jc}}^{{\rm{P}}}\lt {E}_{{\rm{J}}}^{(1)}$	L < L(1)	L0/cusp	N
Perturber-phylic	⋯	⋯	${E}_{{\rm{J}}}^{{\rm{P}}}\lt {E}_{\mathrm{Jc}}^{{\rm{P}}}\lt \min \left[{E}_{{\rm{J}}}^{(1)},{E}_{{\rm{J}}}^{(2)}\right]$	L(1) < L < L(2)	P	N
COM-phylic	⋯	⋯	${E}_{\mathrm{Jc}}^{{\rm{P}}}\lt {E}_{{\rm{J}}}^{(2)}$	L > L(2)	COM	N

Note. Column (1) indicates the name of the orbital family used throughout this paper. Columns (2), (3), and (4) indicate the bounds on the circular part of the Jacobi energy, E_Jc ≈ E_J − κ₀ J_r (κ₀ is the value of the radial epicyclic frequency evaluated at the center of perturbation and J_r is the radial action; see Appendix A for details), evaluated in the neighborhood of L4/L5, L0, and the perturber, i.e., ${E}_{\mathrm{Jc}}^{\left(4\right)}$, ${E}_{\mathrm{Jc}}^{\left(0\right)}$, and ${E}_{\mathrm{Jc}}^{{\rm{P}}}$, respectively. Column (5) indicates the angular momentum, L. Column (6) indicates the center of circulation (COC), where "P" refers to the perturber, and column (7) indicates whether these orbits contribute significantly to dynamical friction (F) or buoyancy (B) or negligibly to either of the two (N). ${E}_{{\rm{J}}}^{\left(k\right)}$, where k = 0, 1, ..., 5, approximately denotes the value of Φ_eff at the k^th Lagrange point (see Appendix A for details), while ${E}_{{\rm{J}}}^{{\rm{P}}}$ denotes that at the location of the perturber (${E}_{{\rm{J}}}^{{\rm{P}}}=-\infty$ for a point mass). ${L}^{\left(k\right)}$, with k = 1, 2, denotes the value of the angular momentum at the k^th Lagrange point. Note that Pac-Man orbits are absent when ${E}_{{\rm{J}}}^{(2)}\gt {E}_{{\rm{J}}}^{(1)}$, which is always the case if the galaxy has a central cusp or the perturber is at large R in a cored galaxy. Orbits that are further away from corotation resonance can cross the separatrix corresponding to L1, L2, or L3, due to changes in J_r , thereby taking on the morphology of a different orbital family, constituting what we call "Chimera orbits" (see Section 4.2 and Appendix B for details).

Download table as:  ASCIITypeset image

Let us start with the yellow-shaded region in Figure 5(a). These are orbits that circulate the galactic center (which coincides with the stable Lagrange point L0 for γ > −1 and the central cusp for γ ≤ −1). These are characterized by ${E}_{{\rm{J}}}^{\left(0\right)}\lt {E}_{\mathrm{Jc}}^{\left(0\right)}\lt {E}_{{\rm{J}}}^{(1)}$ for central cores (γ = 0) and ${E}_{\mathrm{Jc}}^{{\rm{P}}}\lt {E}_{{\rm{J}}}^{(1)}$ for steeper profiles (γ < 0). Additionally, they have lower angular momentum than that at L1, L⁽¹⁾, i.e., they have L < L⁽¹⁾. Their orbital frequency is typically much larger than that of the perturber, and particles on these orbits are thus far from corotation resonance. In what follows, we shall refer to such orbits as "center-phylic." An example is shown in the top row of Figure 6. As is evident from the right-hand panel, the orbital energy varies very little with orbital phase. As a consequence, field particles on these center-phylic orbits exchange very little energy with the perturber, and thus do not contribute significantly to dynamical friction.

Zoom In Zoom Out Reset image size

There is a similar family of nonresonant orbits, with ${E}_{{\rm{J}}}^{{\rm{P}}}\lt {E}_{\mathrm{Jc}}^{{\rm{P}}}\lt \min [{E}_{{\rm{J}}}^{(1)},{E}_{{\rm{J}}}^{(2)}]$, that, in the corotating frame, only circulate the perturber. These orbits, which we call "perturber-phylic," are restricted to the Roche lobe centered on the perturber (shaded light blue in Figure 5(a)). Their angular momentum is higher than that at L1, L⁽¹⁾, but smaller than that at L2, L⁽²⁾, i.e., they have L⁽¹⁾ < L < L⁽²⁾. An example is shown in the middle row of Figure 6. Note that, due to the proximity to the perturber, the orbital energy along this orbit changes drastically and rapidly. Because of the rapid oscillations of orbital energy, the net energy exchange from all field particles on these perturber-phylic orbits is negligible, and this orbital family therefore is also not a significant contributor to dynamical friction.

Next, there is a family of low-E_J orbits with ${E}_{\mathrm{Jc}}^{{\rm{P}}}\lt {E}_{{\rm{J}}}^{(2)}$, that circulate the COM of the combined galaxy+perturber system. Their ZVCs (for near-circular orbits) fall in the unshaded region of Figure 5(a) (outside of the equipotential contour that passes through L2), as their angular momentum prevents them from entering the "central" (shaded) regions, i.e., they have L > L⁽²⁾. An example of such a "COM-phylic" orbit is shown in the bottom row of Figure 6. It reveals small fluctuations in orbital energy on a relatively short timescale. Since there are roughly equal numbers of field particles along each phase of these COM-phylic orbits, they also have a negligible, net contribution to dynamical friction (i.e., at each point in time, these orbits contribute roughly equal numbers of energy gainers as energy losers).

Next, we discuss the three families that are the dominant contributors to dynamical friction. They all have azimuthal frequencies that are comparable to that of the perturber, i.e., Ω_ϕ ≈ Ω_P, such that their libration time in the corotating frame is long. In fact, along these orbits, Ω_ϕ − Ω_P oscillates back and forth about ± Ω_r /N, where Ω_r is the radial frequency and the integer N is the number of radial excursions or epicycles for every libration. The typical range of N is [2, ∞) for realistic galaxy profiles, with N → ∞ marking the corotation resonance, i.e., N is larger the closer the orbit is to corotation. Therefore, they are "near-corotation-resonant" (NCRR), i.e., they librate about the near-corotation resonances, Ω_ϕ − Ω_P = ± Ω_r /N. When the perturber is farther out, M_G(R) ≫ M_P, implying ${{\rm{\Omega }}}_{{\rm{P}}}\approx \sqrt{{{GM}}_{{\rm{G}}}(R)/{R}^{3}}\approx {{\rm{\Omega }}}_{\phi }$ in the vicinity of L4 and L5 (since these two Lagrange points are both at a distance, R, from the galactic center). Therefore, N is large, i.e., N ≫ 1, and the orbits librating about L4 and L5 are close to corotation resonance. As the perturber penetrates deeper into the core region, M_P becomes comparable to M_G(R), and Ω_P significantly exceeds Ω_ϕ near L4 and L5, thereby pushing the orbits farther away from corotation resonance (smaller N), as pointed out by KS18.

The first, and probably most well-known, among the NCRR orbits is the family of so-called "horseshoe orbits," which we already encountered in Section 2. These have ${E}_{\mathrm{Jc}}^{\left(4\right)}\lt {E}_{{\rm{J}}}^{(3)}$ and ${E}_{\mathrm{Jc}}^{\left(0\right)}\gt \max [{E}_{{\rm{J}}}^{(1)},{E}_{{\rm{J}}}^{(2)}]$ for central cores (γ = 0), while ${E}_{\mathrm{Jc}}^{\left(4\right)}\,\lt {E}_{{\rm{J}}}^{(3)}$ and ${E}_{\mathrm{Jc}}^{{\rm{P}}}\gt \max [{E}_{{\rm{J}}}^{(1)},{E}_{{\rm{J}}}^{(2)}]$ for steeper profiles (γ < 0). The ZVCs of near-circular horseshoe orbits fall within the dark blue-shaded region in Figure 5(a) and can only cross the x-axis at the side of L0 opposite to the perturber; the Lagrange point L1 acts as a barrier, forcing the particle to take a long "detour" around the center of the galaxy. They have a net sense of circulation around L3, with a libration frequency ∣Ω_lib∣ ≪ Ω_P. As is evident from the top row of Figure 7 (see also Figure 1), the orbital energy can vary drastically along the orbit, undergoing rapid changes when close to the perturber, where the perturber's force pulls the field particle either inward or outward.

Zoom In Zoom Out Reset image size

Somewhat similar to the horseshoe orbits is a family of orbits that we call "Pac-Man" orbits. These are characterized by ${E}_{\mathrm{Jc}}^{\left(0\right)}\lt {E}_{{\rm{J}}}^{(1)}$ and ${E}_{\mathrm{Jc}}^{{\rm{P}}}\gt {E}_{{\rm{J}}}^{(2)}$ for central cores (γ = 0), while ${E}_{{\rm{J}}}^{(2)}\lt {E}_{\mathrm{Jc}}^{{\rm{P}}}\lt {E}_{{\rm{J}}}^{(1)}$ for steeper profiles (γ < 0). Additionally, they have L⁽¹⁾ < L < L⁽²⁾. They differ from the horseshoe orbits in that they have a net sense of circulation around L0. The Jacobi energy of the near-circular Pac-Man orbits is less than that of L1, which allows their ZVCs to cross the x-axis at the side of L0 that coincides with the perturber, and fall within the green-shaded region of Figure 5(a). Rather than taking a "detour," these orbits can therefore take a "shortcut," which changes their characteristic shape such that they resemble the iconic flashing-dot-eating character of the popular 1980s computer game Pac-Man (see middle row of Figure 7). We emphasize that Pac-Man orbits are only present when ${E}_{{\rm{J}}}^{(1)}\gt {E}_{{\rm{J}}}^{(2)}$. For a given galaxy potential and mass of the perturber, this puts a constraint on the galactocentric distance of the perturber, R. For the Plummer potential and our fiducial mass ratio q = 0.004, Pac-Man orbits are only present when the perturber is located at R ≲ 1.23. When further out, Pac-Man orbits are absent such that the equipotential contours and orbital families are similar for both cored and cuspy galaxy profiles.

The final family of NCRR orbits are known as "tadpole" orbits, a name that again relates to their characteristic shape in the corotating frame (see bottom row of Figure 7). These are characterized by ${E}_{{\rm{J}}}^{(3)}({E}_{{\rm{J}}}^{\left(0\right)})\lt {E}_{\mathrm{Jc}}^{\left(4\right)}\lt {E}_{{\rm{J}}}^{\left(4\right)}={E}_{{\rm{J}}}^{(5)}$ for R > R_bif (R ≤ R_bif), and have a net sense of circulation around either L4 or L5. Their ZVCs fall within the red-shaded region of Figure 5(a).

Along all NCRR orbits (horseshoe, Pac-Man, and tadpole), the energy and angular momentum oscillate with a large amplitude and long period, and the star is up/downscattered through near-corotation resonances by interactions with the perturber. This can be understood in terms of slow and fast action-angle variables, which exist in the neighborhood of a resonance and are related to the radial and azimuthal action-angle variables by a canonical transformation (e.g., Tremaine & Weinberg 1984; Lichtenberg & Lieberman 1992; Chiba & Schönrich 2021). The NCRR orbits librate about the commensurability condition Ω_ϕ − Ω_P ∓ Ω_r /N = 0. The corresponding angle, θ_s = θ_ϕ ∓ θ_r /N − Ω_P t, is called the slow angle, and the action conjugate to it is called the slow action, J_s , which is proportional to the angular momentum. Note that, close to the commensurability condition, d θ_s /dt = Ω_ϕ − Ω_P ∓ Ω_r /N ≃ 0, indicating that θ_s indeed varies slowly. And while it does, the corresponding slow action undergoes large changes. Both J_s and θ_s librate about the near-corotation resonances with a time period, T_lib, that is much larger than the orbital time of the perturber (see Contopoulos 1973; Chiba & Schönrich 2021, for detailed derivations using perturbative expansions of the Hamiltonian around resonances). In fact, for orbits that come arbitrarily close to the separatrices, T_lib approaches infinity.

Contrary to the slow angle, the fast angle, which is nothing but the radial angle, θ_r , varies rapidly along an orbit, while its conjugate action, the fast action, J_f = J_r ± L/N, is nearly invariant. In general, the faster the angle changes, the closer its corresponding fast action is to an adiabatic invariant. Therefore, the NCRR orbits have two integrals of motion, the Jacobi Hamiltonian, E_J (which is exactly conserved), and the fast action, J_f (which is very nearly conserved), and are nearly integrable. ⁴ For the very nearly corotation resonant orbits, N ≫ 1, and therefore J_f ≈ J_r , i.e., the orbital eccentricity (in the inertial frame) remains nearly constant. This is, however, not the case for orbits farther away from corotation resonance, which can show very interesting dynamics, as we shall see shortly.

The relative abundances of the different orbital families depend on the orbital radius R of the perturber. For example, Figure 5(b) shows the equipotential contours of the same Plummer galaxy as in Figure 5(a), but with the perturber orbiting inside the central core, at R = 0.2. Now only four Lagrange points are present; both L1 and L3 have disappeared. As the perturber approaches the galactic center, the Roche lobes around the galactic center and the perturber coalesce to form a single lobe surrounding the perturber. As we show in Paper II, this is associated with the merging, or "bifurcation" of L3, L0, and L1 at a critical bifurcation radius, R_bif, which leaves only L0, L2, L4, and L5, and changes the stability of L0 from being a center to a saddle. As a consequence, neither horseshoe nor center-phylic orbits survive. In addition, the contribution of the tadpole orbits is also significantly diminished. Instead, the dominant orbital families in the central core region are the perturber-phylic orbits and the Pac-Man orbits. As we will see, this has profound implications for dynamical friction.

The orbital configuration is particularly sensitive to the density profile of the galaxy. The lower two panels of Figure 5 show the equipotential contours of a Hernquist galaxy with a perturber at R = 0.5 (Figure 5(c)) and R = 0.2 (Figure 5(d)). In such a cuspy galaxy, there is no L0 (L0 is replaced by the cusp), and the five Lagrange points (L1, L2, L3, L4, and L5) survive throughout, for any value of the orbital radius of the perturber, R, without the occurrence of any bifurcation. As a consequence, in this galaxy potential, there are never any Pac-Man orbits and the relative abundances of different orbital families show a much weaker dependence on R than in the case of the Plummer sphere. How all of this relates to dynamical friction will be discussed in more detail in Sections 5 and 6.

Before proceeding with the computation of the dynamical friction torque from the various orbits, we first discuss a potential complication. We have defined orbital families on the basis of E_Jc, but family is not an invariant property for all orbits. In fact, an orbit can change its family in the course of its evolution. This is because the orbit determinant, E_Jc, as expressed in Equation (22), is not an invariant quantity. It not only involves E_J, which is an integral of motion and thus conserved, but also the radial action, J_r , which is typically not constant along an orbit. In particular, J_r can undergo significant changes along orbits that are farther away from corotation resonance, since only a linear combination of J_r and L, and not J_r alone, is the fast action in this case. Therefore, the value of E_Jc can potentially cross over from that corresponding to one orbital family to another, which corresponds to the orbit undergoing separatrix crossing due to a change in the radial action enabled by the perturber, altering its morphological appearance. We call such orbits "Chimera orbits." ⁵ These Chimera-like transitions occur between trapped regions of neighboring resonances on either side of a separatrix (see Appendix A) or a chaotic island formed by the overlap of resonances (see Chiba & Schönrich 2021, for a detailed discussion in the context of bar-like perturbations). For example, the metamorphosis between horseshoes and tadpoles occurs near L3, while that between horseshoes, Pac-Mans, and center-phylic orbits happens near L1. Finally, the transition between Pac-Man, COM-phylic, and perturber-phylic orbits occurs in the neighborhood of L2. We show several examples of such Chimera orbits in Appendix B. Not all orbits show this Chimera behavior. The very nearly corotation resonant orbits are nearly circular and thus have small J_r . Since J_r is a fast action along such orbits, it remains almost constant, i.e., the orbits remain nearly circular and do not exhibit Chimera characteristics.

When the separatrix crossing along a Chimera orbit results in a perturber-phylic phase, we speak of resonant capture (Henrard 1982), which as pointed out in Tremaine & Weinberg (1984), can "dress" the perturber with a cloud of captured stars. Note, though, that in the "slow" regime considered here, in which the orbital radius of the perturber is taken to be invariant, these stars can undergo separatrix crossing again, transitioning back to a Pac-Man or a COM-phylic orbit. Similarly, when a separatrix crossing results in a transition from a "trapped" NCRR state to an "untrapped" COM-phylic state, the transition is sometimes called "scattering" (e.g., Daniel & Wyse 2015).

Chimera orbits are difficult to account for in our treatment because they do not have a clear periodic behavior, i.e., they do not have a well-defined libration time. However, we find that most of them typically behave as an archetypal orbit of their family for many orbital periods before revealing their Chimera nature, i.e., they are "semi-ergodic" (similar to the semi-ergodic orbits identified by Athanassoula et al. 1983 in their study of barred galaxies). This is akin to how Arnold diffusion in KAM theory can cause chaotic orbits to behave quasi-regularly for extended periods (e.g., Lichtenberg & Lieberman 1992). Hence, we conjecture that their relevance to dynamical friction is captured, at least to leading order, by our following treatment of the NCRR orbital families.

In order to compute the torque on the perturber due to a single orbit, we proceed as follows. We numerically integrate the orbit of a massless field particle in the presence of the perturber, registering its position r , velocity $\dot{{\boldsymbol{r}}}$, energy E, and angular momentum L , as a function of time $t^{\prime}$. We use $t^{\prime}$ to indicate the phase of a particle along this orbit. We have seen in Sections 2 and 4.2 that, as a particle moves along the perturbed orbit, it undergoes changes in energy and angular momentum due to exchanges with the perturber. Hence, after some time Δt, a particle starting from phase $t^{\prime}$ has transferred a net amount of energy ${\rm{\Delta }}E({\rm{\Delta }}t)=E(t^{\prime} +{\rm{\Delta }}t)-E(t^{\prime} )$ to the perturber. Here, $E(t^{\prime} )$ is the perturbed energy of a particle at phase $t^{\prime}$, given by Equation (11). To work out the total energy exchanged with the perturber by all stars associated with the orbit in question, we need to integrate ΔE(Δt) along the orbit, weighted by the relative number of stars at each point along the orbit. This weight is given by ${f}_{0}({E}_{0{\rm{G}}}(t^{\prime} ))$, with f₀ the unperturbed DF, and

$$\begin{eqnarray}&&{E}_{0{\rm{G}}}(t^{\prime} )=\displaystyle \frac{1}{2}{\left|\dot{{\boldsymbol{r}}}+{{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times {\boldsymbol{r}}-{{\boldsymbol{v}}}_{{\bf{G}}}\right|}^{2}+{{\rm{\Phi }}}_{{\rm{G}}}\end{eqnarray} \tag{ 23 }$$

the galactocentric energy of the star at phase $t^{\prime}$ in absence of the perturber, where ${{\boldsymbol{v}}}_{{\bf{G}}}=-{{\rm{\Omega }}}_{{\rm{P}}}{q}_{{\rm{G}}}R\,\hat{y}$ is the circular velocity of the galactic center about the COM. If we use $s(t^{\prime} )$ to parameterize the path length along the phase-space trajectory traced out by the orbit, then the total energy exchanged with the perturber along this orbit, some time Δt after the perturber was introduced, is given by the following line integral:

$$\begin{eqnarray}&&{\rm{\Delta }}E({\rm{\Delta }}t)=\displaystyle \frac{1}{{ \mathcal A }}{\int }_{s}{ds}(t^{\prime} )\left[E(t^{\prime} +{\rm{\Delta }}t)-E(t^{\prime} )\right]\,{f}_{0}({E}_{0{\rm{G}}}(t^{\prime} )),\end{eqnarray} \tag{ 24 }$$

where ${ \mathcal A }$ is a normalization factor (see below).

Typically, an orbit in the corotating frame will not be exactly closed and the integration limit therefore will have no boundaries. However, for the NCRR orbits discussed in Section 4.2, the orbit is approximately periodic in the corotating frame, with a period T_lib set by the time it takes the particle to librate about its Lagrange point (the COC in column 6 of Table 1), which we compute by a Fourier analysis of the orbit in the corotating frame. In the vicinity of the stable Lagrange points, L4 and L5, T_lib can be analytically computed using a perturbative method, as discussed in Paper II.

The line integral in Equation (24) has to be performed along the phase-space trajectory, and therefore the differential line element ${ds}(t^{\prime} )$ is given by ${ds}=\sqrt{{\left|d{{\boldsymbol{r}}}_{i}\right|}^{2}+{\left|d{\dot{{\boldsymbol{r}}}}_{i}\right|}^{2}}$. Knowing that the Jacobian for the transformation from $t^{\prime}$ to the arc length $s(t^{\prime} )$ is given by

$$\begin{eqnarray}&&\displaystyle \frac{{ds}}{{dt}^{\prime} }=\sqrt{{\left|{\dot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}+{\left|{\ddot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}},\end{eqnarray} \tag{ 25 }$$

with ${\dot{{\boldsymbol{r}}}}_{{\bf{i}}}$ and ${\ddot{{\boldsymbol{r}}}}_{{\bf{i}}}$ being the velocity and acceleration in the inertial frame, respectively, we can approximate the line integral as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}E({\rm{\Delta }}t) & \approx & \displaystyle \frac{1}{{ \mathcal A }}{\displaystyle \int }_{0}^{{T}_{\mathrm{lib}}}{dt}^{\prime} \,\sqrt{{\left|{\dot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}+{\left|{\ddot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}}\\ & & \times \,\left[E(t^{\prime} +{\rm{\Delta }}t)-E(t^{\prime} )\right]\,{f}_{0}({E}_{0{\rm{G}}}(t^{\prime} )),\end{array}\end{eqnarray} \tag{ 26 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal A } & = & {\displaystyle \int }_{s}{ds}(t^{\prime} ){f}_{0}({E}_{0{\rm{G}}}(t^{\prime} ))\\ & = & {\displaystyle \int }_{0}^{{T}_{\mathrm{lib}}}{dt}^{\prime} \,\sqrt{{\left|{\dot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}+{\left|{\ddot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}}\,{f}_{0}({E}_{0{\rm{G}}}(t^{\prime} )).\end{array}\end{eqnarray} \tag{ 27 }$$

Note that, with this normalization, ΔE(Δt) is the average energy per star exchanged with the perturber in a time Δt along the orbit in question.

The inertial acceleration vector is given by

$$\begin{eqnarray}&&{\ddot{{\boldsymbol{r}}}}_{{\bf{i}}}=-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 28 }$$

where Φ = Φ_P + Φ_G is the total potential, while the velocity vector in the inertial frame is related to that in the corotating frame, $\dot{{\boldsymbol{r}}}$, by

$$\begin{eqnarray}&&{\dot{{\boldsymbol{r}}}}_{{\bf{i}}}=\dot{{\boldsymbol{r}}}+{{\boldsymbol{\Omega }}}_{{\boldsymbol{P}}}\times {\boldsymbol{r}}.\end{eqnarray} \tag{ 29 }$$

We perform this line integral for the three NCRR orbits (horseshoe, Pac-Man, and tadpole) shown in Figure 7. All three orbits correspond to our fiducial q = 0.004 point-mass perturber in a Plummer potential at R = 0.5. The solid blue, dotted–dashed green, and dotted red lines in Figure 8 show the resulting ΔE for the horseshoe, Pac-Man, and tadpole orbits, respectively, as function of Δt. Note that ΔE(Δt = T_lib) = 0; as discussed in Section 2, along each NCRR orbit, particles both gain and lose energy, and the net effect for a single particle over a full libration period is zero. However, due to the nonuniform phase distribution along each orbit, which arises from the unperturbed phase-space distribution, f₀(E_0G), we see that ΔE is positive for all 0 < Δt < T_lib. A positive ΔE indicates that the field particles along these orbits gain net energy from the perturber, and thus that the perturber experiences dynamical friction. As the field particles gain energy, their Ω_ϕ decreases. The perturber in turn loses energy and falls in, with increasing Ω_P. This puts the original NCRR orbits out of near-corotation resonance. Therefore, ΔE(Δt) is only relevant for the dynamics of the system for relatively small Δt. The exact choice of Δt to consider is somewhat ambiguous; it should be indicative of the timescale over which the perturber moves through the resonances, which in turn depends on the strength of dynamical friction. In what follows, we take Δt = T_orb, the orbital time of the perturber. None of our qualitative conclusions are sensitive to this particular choice.

Zoom In Zoom Out Reset image size

The solid blue, dotted–dashed green, and dotted red curves in Figure 8 correspond to NCRR orbits in the case where the perturber is orbiting at R = 0.5, just outside the core of the Plummer sphere. For comparison, the dashed, green curve in Figure 8 indicates the ΔE(Δt) for a Pac-Man orbit in the case where the perturber is at R = 0.2, well inside the core of the Plummer galaxy. In this case, ΔE is negative, indicating that this orbit contributes a positive, enhancing torque. Note that, since the torque on the field particle is given by ${dL}/{dt}={{\rm{\Omega }}}_{{\rm{P}}}^{-1}({dE}/{dt})$ (see Equation (19)), the average torque on the perturber due to an orbit between t = 0 and t = T_orb is equal to $-{\rm{\Delta }}E/\left({{\rm{\Omega }}}_{{\rm{P}}}{T}_{\mathrm{orb}}\right)=-{\rm{\Delta }}E/(2\pi )$, i.e., $\mathrm{sign}({ \mathcal T })=-\mathrm{sign}({\rm{\Delta }}E)$. We thus see that some of the NCRR orbits can give rise to dynamical buoyancy, rather than friction. An explanation of the latter is discussed in Section 6.

Having demonstrated how to compute the contribution to dynamical friction from individual orbits, in the form of ΔE(T_orb), one can in principle obtain the total torque by summing over all orbits, properly weighted by their relative contribution to the distribution function. In practice, though, this is far from trivial. First of all, sampling all orbits numerically is tedious to the point that one is better off just running an N-body simulation. Second, some orbits are difficult to integrate accurately, especially some Chimera orbits that reveal semi-ergodic behavior, and the perturber-phylic orbits along which the energy varies rapidly with time. Hence, the nonperturbative method adopted in this paper is not well-suited to an accurate computation of the total dynamical friction torque. Nonetheless, it gives valuable insight as to the inner workings, in an orbit-based sense, of dynamical friction and buoyancy.

As an example, we now proceed to investigate the contribution to the torque, in terms of ΔE(T_orb), from a modest subsample of orbits. In what follows, we continue to treat the dynamics in 2D (i.e., we only consider orbits in the x–y plane depicted in Figure 3). We densely sample the part of the orbital parameter space corresponding to the NCRR orbits, which is most relevant for dynamical friction. We first sample the starting point (x₀, y₀) by setting y₀ = 0 and sampling x₀ uniformly over the range dominated by the NCRR horseshoe and Pac-Man orbits (roughly the region inside the ${E}_{{\rm{J}}}^{(2)}$ separatrix marked by the solid line in Figure 5). Note that, by sampling orbits that intersect the x-axis, we exclude tadpole orbits with large E_J that librate in small regions around L4 and L5. After sampling x₀, we uniformly sample E_J over the range $[{{\rm{\Phi }}}_{\mathrm{eff}}({x}_{0},0),{E}_{{\rm{J}}}^{\left(4\right)}]$. Although orbits with ${E}_{{\rm{J}}}\gg {E}_{{\rm{J}}}^{\left(4\right)}$ are far from corotation resonance, thereby contributing less to dynamical friction, those with small, positive values of ${E}_{{\rm{J}}}-{E}_{{\rm{J}}}^{\left(4\right)}$ are NCRR and have similar contribution to the torque as those with ${E}_{{\rm{J}}}\lesssim {E}_{{\rm{J}}}^{\left(4\right)}$. Therefore, we consider ${E}_{{\rm{J}}}^{\left(4\right)}$ to be only an approximate rather than a hard cutoff for the NCRR orbits. Finally, we sample the initial velocities v_x,0 and v_y,0 under the constraint that

$$\begin{eqnarray}&&{v}_{0}=\sqrt{{v}_{{\rm{x}},0}^{2}+{v}_{{\rm{y}},0}^{2}}=\sqrt{2[{E}_{{\rm{J}}}-{{\rm{\Phi }}}_{\mathrm{eff}}({x}_{0},0)]}.\end{eqnarray} \tag{ 30 }$$

Note that both x ₀ and v ₀ are defined in the corotating frame. We numerically integrate the orbits for 100 T_orb, with T_orb being the orbital time of the perturber, after which we estimate the libration time, T_lib, by noting the consecutive time stamps at which each orbit crosses the abscissa of its center of circulation (see Table 1) after making a 2π circulation about it. Finally, we compute ΔE ≡ ΔE(T_orb∣x₀, E_J, v_x,0, v_y,0) using Equation (26).

In order to allow for a meaningful comparison of the torque contribution from each of these orbits, we weight the ΔE per star, given by Equations (26)–(27), by the average phase-space density associated with that orbit. This yields the total energy exchange per unit phase space from an orbit, given by

$$\begin{eqnarray}&&{\left({\rm{\Delta }}E\right)}_{{\rm{w}}}\equiv \displaystyle \frac{{\int }_{s}{ds}(t^{\prime} ){f}_{0}({E}_{0{\rm{G}}}(t^{\prime} ))}{{\int }_{s}{ds}(t^{\prime} )}\,{\rm{\Delta }}E.\end{eqnarray} \tag{ 31 }$$

Given that the time-averaged torque (per unit phase space) on the perturber contributed by an individual orbit is yielded by

$$\begin{eqnarray}&&{{ \mathcal T }}_{{\rm{w}}}=-\displaystyle \frac{1}{{{\rm{\Omega }}}_{{\rm{P}}}}\,\displaystyle \frac{{\left({\rm{\Delta }}E\right)}_{{\rm{w}}}}{{\rm{\Delta }}t},\end{eqnarray} \tag{ 32 }$$

(see Equation (19)), we know that the torque per unit phase space contributed by the orbit can be expressed as

$$\begin{eqnarray}&&{{ \mathcal T }}_{{\rm{w}}}=-\displaystyle \frac{1}{2\pi }\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&\times \displaystyle \frac{{\int }_{0}^{{T}_{\mathrm{lib}}}{dt}^{\prime} \,\sqrt{{\left|{\dot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}+{\left|{\ddot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}}\left[E(t^{\prime} +{\rm{\Delta }}t)-E(t^{\prime} )\right]{f}_{0}({E}_{0{\rm{G}}}(t^{\prime} ))}{{\int }_{0}^{{T}_{\mathrm{lib}}}{dt}^{\prime} \,\sqrt{{\left|{\dot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}+{\left|{\ddot{{\boldsymbol{r}}}}_{{\bf{i}}}\right|}^{2}}},\end{eqnarray} \tag{ 34 }$$

where we have used the fact that we adopt Δt = T_orb = 2π/Ω_P, and we have rewritten ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}$ using Equations (26) and (27).

Figure 9 plots ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}$ for the Plummer sphere as a function of x₀ and E_J for $\left|{v}_{{\rm{y}},0}\right|=\tfrac{1}{3}{v}_{0}$. Results for other values of $\left|{v}_{{\rm{y}},0}\right|$ are very similar, but with the overall amplitudes in ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}$ decreasing as $\left|{v}_{{\rm{y}},0}\right|\to {v}_{0}$. For each $\left({x}_{0},{E}_{{\rm{J}}},\left|{v}_{{\rm{y}},0}\right|\right)$, there are four combinations of (v_x,0, v_y,0), given by $(\pm \sqrt{{v}_{0}^{2}-{v}_{{\rm{y}},0}^{2}},\pm \left|{v}_{{\rm{y}},0}\right|)$. The values of ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}$ shown are the sums of these four cases combined.

Zoom In Zoom Out Reset image size

The left- and right-hand panels correspond to R = 0.5 and R = 0.2, respectively. They show the results for a total of 4 × 2500 different orbits. Redder colors denote more positive values of ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}$ (i.e., stronger dynamical friction), while bluer colors indicate more negative values (i.e., more pronounced dynamical buoyancy). Note that, for R = 0.5, i.e., when the perturber is outside the core, ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}$ is predominantly positive, indicating that nearly all the NCRR orbits (horseshoes, Pac-Mans, and some tadpoles, with x₀ on either side of L3 and L1) exert a retarding torque (i.e., dynamical friction). However, when R = 0.2 and the perturber is orbiting inside the core, almost the entire orbital parameter space (dominated by the NCRR Pac-Mans and tadpoles) contributes to dynamical buoyancy (i.e., ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}\lt 0$). Clearly, there is a profound transition in the total torque once the perturber enters the core.

When the perturber is outside the core (left-hand panel), the contribution from the center-phylic orbits, which occupy the range of x₀ on either side of the galactic center (L0, marked by the vertical, dashed line) is completely negligible. The same holds for the COM-phylic orbits near the leftmost edge of the plot. When the perturber is inside the core (right-hand panel), one again sees that orbits with starting positions close to L₀ contribute a negligible torque. Unlike in the left-hand panel, though, these are not center-phylic orbits. After all, those vanish when the perturber crosses the bifurcation radius. Rather, these are predominantly Pac-Man and tadpole orbits, but unlike their counterparts with starting positions a bit further away from the (unstable) L₀, they happen to exert negligible torque. Note that some of the COM-phylic orbits with x₀/R ≲ − 0.75 also contribute a (positive) torque. Their net contribution, though, is significantly smaller than that from the NCRR Pac-Man orbits, and rapidly weakens when x₀/R becomes smaller (i.e., further away from the galactic center).

Figure 10 is the same as Figure 9, but for our Hernquist galaxy. For both R = 0.5 (left-hand panel) and R = 0.2 (right-hand panel), it is clear that the total torque is negative (retarding) and dominated by the NCRR orbits. Most importantly, there is no transition in the sign of the total torque as one approaches the center, consistent with the notion that buoyancy and core-stalling are absent if the central density profile is cuspy. Another difference with respect to the Plummer sphere is that, while there is no significant contribution to the torque from the COM-phylic orbits for R = 0.5 nor for R = 0.2, the center-phylic orbits now make a significant contribution to the total torque. Although each of these orbits has a very small ΔE(T_orb), the steepness of the distribution function toward the galactic center means that they are abundant, thus receiving a large weight. When the perturber is at R = 0.5, there are roughly equal numbers of center-phylic orbits with positive and negative ${\left({\rm{\Delta }}E\right)}_{{\rm{w}}}$ (note the alternating red and blue stripes on either side of the galactic center). As a consequence, the net torque contribution from the entire population of center-phylic orbits is small.

Zoom In Zoom Out Reset image size

Finally, we emphasize that the above inventory of the torque from individual orbits is incomplete. First of all, we have restricted the range of x₀ such that it does not include any perturber-phylic orbits. The reason is that they are difficult to integrate, while their contribution to the torque is negligible for reasons discussed in Section 4.2. Second, by only picking starting points along the x-axis, we have selected against tadpole orbits with large E_J, which are typically confined to small regions centered on L4 or L5. We have examined several such orbits and found their behavior to be very similar to that of the horseshoe and Pac-Man orbits in terms of their contribution to the torque. Third, we have restricted the E_J values of the orbits up to ${E}_{{\rm{J}}}^{\left(4\right)}={E}_{{\rm{J}}}^{(5)}$. This is because orbits with ${E}_{{\rm{J}}}\gg {E}_{{\rm{J}}}^{\left(4\right)}$ are far from corotation resonance (with drift time steeply falling with increasing E_J) and consequently less important for dynamical friction. However, orbits with small, positive values of ${E}_{{\rm{J}}}-{E}_{{\rm{J}}}^{\left(4\right)}$ have similar contribution to the torque as those with ${E}_{{\rm{J}}}\lesssim {E}_{{\rm{J}}}^{\left(4\right)}$. Thus, we use ${E}_{{\rm{J}}}^{\left(4\right)}$ only as an approximate cutoff for the NCRR orbits. Finally, and most significantly, we have only considered orbits of field particles confined to the orbital plane of the perturber, i.e., those with z = 0 and v_z = 0. We presume that this does not significantly impact any of our conclusions regarding the contributions of the NCRR horseshoe and Pac-Man orbits, as the third dimension merely allows for an additional vertical oscillation not accounted for in our 2D planar treatment (in particular, no new orbital families are introduced by allowing motion in the z-direction, since there exist no Lagrange points off the orbital plane). However, the relative contributions of the different NCRR orbits to the total torque may be significantly different from what emerges from the 2D analysis presented here. In particular, the tadpole orbits would dominate the phase space and therefore might contribute more significantly to the overall torque in 3D. This is a caveat of our approach that we leave for future work.

L4 and L5 are located at a distance r_G = R from the galactic center, and at an angle φ = ± π/3. The nearest saddle point to L4/L5 is L3, at a distance of r_G = r_sep = r₃ from the galactic center, and with φ = φ_sep = π. In the region centered on L4/L5 and bounded by the L3 separatrix, the disturbing potential Φ₁, given by Equation (A4), is bounded by

$$\begin{eqnarray}&&\begin{array}{l}-q\left[\displaystyle \frac{1}{{r}_{3}+R}+\displaystyle \frac{{r}_{3}}{{R}^{2}}\right]={{\rm{\Phi }}}_{1}({r}_{3},\pi )\\ \quad \lt {{\rm{\Phi }}}_{1}\lt {{\rm{\Phi }}}_{1}(R,\pm \pi /3)=-\displaystyle \frac{q}{2R}.\end{array}\end{eqnarray} \tag{ A9 }$$

Therefore, in the vicinity of L4, L5, and L3, J_ϕ has real roots for φ that are restricted to 0 < φ < π (or − π < φ < 0) when ${h}_{0}^{{\prime} (4)}+{{\rm{\Phi }}}_{1}({r}_{3},\pi )\lt {E}_{\mathrm{Jc}}^{{\prime} (4)}\lt {h}_{0}^{{\prime} (4)}+{{\rm{\Phi }}}_{1}(R,\pm \pi /3)$, i.e., when the "trapping criterion,"

$$\begin{eqnarray}&&\begin{array}{l}{h}_{0}^{{\prime} (4)}-q\left[\displaystyle \frac{1}{{r}_{3}+R}+\displaystyle \frac{{r}_{3}}{{R}^{2}}\right]={E}_{{\rm{J}}}^{{\prime} (3)}\\ \quad \lt {E}_{\mathrm{Jc}}^{{\prime} (4)}\lt {E}_{{\rm{J}}}^{{\prime} (4)}={h}_{0}^{{\prime} (4)}-\displaystyle \frac{q}{2R},\end{array}\end{eqnarray} \tag{ A10 }$$

is satisfied. Here, ${E}_{\mathrm{Jc}}^{{\prime} (4)}$ and ${h}_{0}^{{\prime} (4)}$ are, respectively, the circular part of the Jacobi energy given by Equation (A8) and the unperturbed Jacobi energy, both evaluated at the COP, L4/L5. These "trapped" orbits that lie inside the L3 separatrix are the tadpole orbits. The "untrapped" orbits that lie beyond the L3 separatrix (φ > φ_sep = π), i.e., that satisfy the condition,

$$\begin{eqnarray}&&{E}_{\mathrm{Jc}}^{{\prime} (4)}\lt {E}_{{\rm{J}}}^{{\prime} (3)},\end{eqnarray} \tag{ A11 }$$

are the horseshoe orbits. The tadpole orbits are trapped, librating around L4/L5, while the horseshoe orbits are untrapped w.r.t. L4/L5. However, as we shall see shortly, the horseshoe orbits are trapped inside the L1/L2 separatrix and are therefore still librating about near-corotation resonances. Since J_r can vary along an orbit, certain orbits with high J_r can cross the L3 separatrix, and the trapped tadpoles and untrapped horseshoes can metamorphose into each other, showing Chimera behavior (see Appendix B).

A perturbative analysis around L0 can only be performed for a perfect central core, i.e., $\gamma \equiv {\mathrm{lim}}_{r\to 0}d\mathrm{log}\rho /d\mathrm{log}r=0$. For density profiles with γ < 0, Ω₀ and κ₀ diverge like ${r}_{{\rm{G}}}^{\gamma /2}$ toward the galactic center, around which the perturbative expansion of ${E}_{{\rm{J}}}^{{\prime} }$ given in Equation (A1) is thus not defined.

For a cored galaxy, L0 is a stable Lagrange point located at the galactic center, i.e., r_G = 0. The nearest saddle point to L0 is L1, located along the x-axis (φ = φ_sep = 0) at a distance r_G = r_sep = r₁ from the galactic center. In the region centered on L0 and bounded by the L1 separatrix, Φ₁ is bounded by

$$\begin{eqnarray}&&q\left[-\displaystyle \frac{1}{\left|{r}_{1}-R\right|}+\displaystyle \frac{{r}_{1}}{{R}^{2}}\right]\lt {{\rm{\Phi }}}_{1}\lt q\left[-\displaystyle \frac{1}{R}+\displaystyle \frac{{r}_{1}}{2{R}^{2}}\right].\end{eqnarray} \tag{ A12 }$$

Hence, in the vicinity of L1, the roots of J_φ are real for restricted values of φ when ${E}_{\mathrm{Jc}}^{{\prime} (0)}\gt {h}_{0}^{{\prime} (0)}+{{\rm{\Phi }}}_{1}({r}_{1},0)$, i.e., when

$$\begin{eqnarray}&&{E}_{\mathrm{Jc}}^{{\prime} (0)}\gt {E}_{{\rm{J}}}^{{\prime} (1)}={h}_{0}^{{\prime} (0)}+q\left[-\displaystyle \frac{1}{\left|{r}_{1}-R\right|}+\displaystyle \frac{{r}_{1}}{{R}^{2}}\right].\end{eqnarray} \tag{ A13 }$$

Here, ${E}_{\mathrm{Jc}}^{{\prime} (0)}$ and ${h}_{0}^{{\prime} (0)}$ are, respectively, the circular part of the Jacobi energy given by Equation (A8) and the unperturbed Jacobi energy, both evaluated at the COP, L0. These trapped orbits that lie inside the L1 separatrix are the horseshoe orbits. These orbits are therefore in a librating state around near-corotation resonances (despite being untrapped w.r.t. L4/L5). The untrapped orbits that lie outside the L1 separatrix, and that satisfy the condition,

$$\begin{eqnarray}&&{E}_{\mathrm{Jc}}^{{\prime} (0)}\lt {E}_{{\rm{J}}}^{{\prime} (1)},\end{eqnarray} \tag{ A14 }$$

are the Pac-Man and the center-phylic orbits. The Pac-Mans have higher angular momentum (L⁽¹⁾ < L < L⁽²⁾) than the center-phylic orbits (L < L⁽¹⁾), where ${L}^{\left(k\right)}$, with k = 1, 2, denotes the value of L at the k^th Lagrange point. Although the Pac-Mans are beyond the L1 separatrix, they are still trapped inside the L2 separatrix (as we shall see shortly) and therefore librating about near-corotation resonances. The center-phylic orbits on the other hand rotate about the galactic center. Since J_r can vary along an orbit, certain orbits with high J_r can cross the L1 separatrix, resulting in Chimera-like metamorphosis between the trapped horseshoe and untrapped Pac-Man and center-phylic orbital families (see Appendix B).

In the vicinity of the perturber (i.e., the region centered on the perturber and bounded by the L2 separatrix), for a given r_G, Φ₁ varies in the range

$$\begin{eqnarray}&&q\left[-\displaystyle \frac{1}{\left|{r}_{{\rm{G}}}-R\right|}+\displaystyle \frac{{r}_{{\rm{G}}}}{{R}^{2}}\right]\lt {{\rm{\Phi }}}_{1}\lt q\left[-\displaystyle \frac{1}{R}+\displaystyle \frac{{r}_{{\rm{G}}}}{2{R}^{2}}\right].\end{eqnarray} \tag{ A15 }$$

The nearest saddle point to the perturber is L2, located along the x-axis (φ = φ_sep = 0) at a distance r_G = r_sep = r₂ from the galactic center. Hence, in the neighborhood of L2, J_φ has real roots for restricted values of φ when ${E}_{\mathrm{Jc}}^{{\prime} {\rm{P}}}\gt {h}_{0}^{{\prime} {\rm{P}}}+{{\rm{\Phi }}}_{1}({r}_{2},0)$, i.e., when

$$\begin{eqnarray}&&{E}_{\mathrm{Jc}}^{{\prime} {\rm{P}}}\gt {E}_{{\rm{J}}}^{{\prime} (2)}={h}_{0}^{{\prime} {\rm{P}}}+q\left[-\displaystyle \frac{1}{\left|{r}_{2}-R\right|}+\displaystyle \frac{{r}_{2}}{{R}^{2}}\right].\end{eqnarray} \tag{ A16 }$$

Here, ${E}_{\mathrm{Jc}}^{{\prime} {\rm{P}}}$ and ${h}_{0}^{{\prime} {\rm{P}}}$ are, respectively, the circular part of the Jacobi energy given by Equation (A8) and the unperturbed Jacobi energy, both evaluated at the COP, which in this case is the perturber. These trapped orbits that lie inside the L2 separatrix are Pac-Mans when ${E}_{{\rm{J}}}^{{\prime} (1)}\gt {E}_{{\rm{J}}}^{{\prime} (2)}$ and horseshoes otherwise. The Pac-Mans are therefore in a librating state about near-corotation resonances, even if they are untrapped w.r.t. L0. The untrapped orbits that lie outside the L2 separatrix, and that satisfy the condition,

$$\begin{eqnarray}&&{E}_{\mathrm{Jc}}^{{\prime} {\rm{P}}}\lt {E}_{{\rm{J}}}^{{\prime} (2)},\end{eqnarray} \tag{ A17 }$$

are the COM-phylic, perturber-phylic, and center-phylic (for γ < 0 profiles) orbits. The COM-phylic orbits have higher angular momentum (L > L⁽²⁾) than the perturber-phylic orbits (L⁽¹⁾ < L < L⁽²⁾). While the COM-phylic orbits rotate about the COM of the galaxy-perturber system, the perturber-phylic orbits rotate about the perturber. Due to variation of J_r along an orbit, Chimera-like transitions can occur between the trapped Pac-Man and untrapped COM-phylic and perturber-phylic orbital families (see Appendix B).