Scattered Disk Dynamics: The Mapping Approach

We derive, and discuss the properties of, a symplectic map for the dynamics of bodies on nearly parabolic orbits. The orbits are perturbed by a planet on a circular, coplanar orbit interior to the pericenter of the parabolic orbit. The map shows excellent agreement with direct numerical integrations and elucidates how the dynamics depends on perturber mass and pericenter distance. We also use the map to explore the onset of chaos, statistical descriptions of chaotic transport, and sticking in mean-motion resonances. We discuss implications of our mapping model for the dynamical evolution of the solar system's scattered disk and other highly eccentric trans-Neptunian objects.


INTRODUCTION
Small bodies beyond the orbit of Neptune are understood to be the remnants of a formerly much larger planetesimal population that was dramatically sculpted by dynamical interactions with the outer planets during the solar system's early evolution (e.g., Gladman & Volk 2021).The orbits of these trans-Neptunian objects (TNOs) have been influenced to varying degrees by the subsequent billions of years of evolution due to various secular and resonant dynamical phenomena, the Galactic tide and kicks from passing stars, as well as chaos (see review by Saillenfest 2020).Understanding this long-term dynamical evolution is necessary for connecting the orbits of the TNOs we observe today to the solar system's formation and early evolution.Similar mechanisms drive the evolution of the debris disks found around many young stars (Wyatt 2020) and presumably govern the production of interstellar objects such as 'Oumuamua.
Many TNOs are on highly eccentric orbits.The most significant population of highly eccentric TNOs is the "scattered disk", comprised of objects for which chaotic dynamics induced by Neptune scattering leads to quasi-random semi-major axis evolution.The population of "detached" TNOs resides on large semi-major axis, high-eccentricity orbits similar to those of the scattered population, but these are not actively experiencing chaotic semi-major axis evolution because their perihelion distances are large enough that they do not make close approaches to Neptune or the other planets.How these bodies were placed on their current orbits, dynamically isolated from any significant influence by the solar system's giant planets, is a puzzle.Finally, Neptune's distant mean-motion resonances (MMRs) are inferred to host large small-body populations on high-eccentricity orbits ( ≳ 0.5; Pike et al. 2015;Volk et al. 2018;Crompvoets et al. 2022).It is not clear whether these populations are comprised entirely of scattered-disk objects experiencing transient periods of "resonance sticking" (e.g., Yu et al. 2018) or if some fraction of the populations are permanent captures from Neptune's outward migration early in the solar system's history.The high-eccentricity scattered, detached, and resonant populations are sometimes referred to collectively as "extreme TNOs" or eTNOs (e.g., de La Fuente Marcos & de La Fuente Marcos 2014).Drawing inferences about the solar system's dynamical history from the populations of eTNOs requires understanding the dynamics of highly eccentric orbits subject to planetary perturbations.
Analytic methods are desirable as a supplement to direct numerical simulations of eTNO dynamical evolution.Unfortunately, classical disturbing function expansions are generally of little use when considering the dynamics of the large eccentricities typical of eTNOs (though see Batygin et al. 2021).To overcome these shortcomings, numerous previous studies have turned to mapping approaches to model the behaviour of highly eccentric orbits (e.g., Petrosky 1986;Petrosky & Broucke 1988;Chirikov & Vecheslavov 1989;Malyshkin & Tremaine 1999;Pan & Sari 2004;Shevchenko 2011;Hadden et al. 2018).These mapping models approximate particles' orbital evolution by considering the effects of successive impulsive energy "kicks" experienced by the particles during their pericenter passage.Mapping models thus provide a computationally efficient alternative to direct -body simulations, which require time steps short enough to resolve Neptune's orbital period.
In addition to the benefits of computational efficiency, the mapping approach lends itself to the application of theoretical results on planar area-preserving maps.There exists a rich theory of chaos and transport in planar symplectic maps (see reviews by Mackay et al. 1984;Meiss 1992Meiss , 2015)).This body of theory can be applied, for example, to understand when large-scale chaos is expected to occur, to describe how chaotic trajectories diffuse, and to explicate the role of "resonance sticking" in producing power-law escape time distributions from regions of phase space.All of these aspects of Hamiltonian dynamics are relevant to understanding the long-term evolution of scattered-disk objects and other eTNOs.
In this paper, we develop a mapping approach to the dynamics of highly eccentric orbits subject to an inner, planetary-mass perturber.We focus on the dynamics of eTNOs perturbed by Neptune, though the results can readily be generalized to other contexts.This mapping provides an analytic model for the description of the mean-motion resonance dynamics of highly eccentric orbits.It also allows us to derive a statistical description of the evolution of scattered bodies using a Fokker-Planck approach.
We restrict ourselves to a solar system containing a single giant planet on a circular orbit (usually taken to have Neptune's mass and semimajor axis) and to bodies orbiting in the same plane as the planet.The latter restriction avoids complications such as Kozai-Lidov oscillations and should be a reasonable approximation for objects in the scattered and detached disks, which typically have only modest inclinations.
Our paper is structured as follows.In Section 2, we introduce our iterative map equations and also derive a "local" approximation to those equations, valid when particles' fractional variations in orbital energy are small.In Section 3, we use our map equations to derive an integrable Hamiltonian model for the dynamics of mean-motion resonances (MMRs) of arbitrary order in the high-eccentricity regime.We also apply our MMR model to predict the boundary between chaotic and regular motion as a function of both pericenter distance and semi-major axis.In Section 4, we present a statistical description of the evolution of chaotic orbits and derive a Fokker-Planck equation governing the evolution of scattered bodies' orbital energies.Section 5 summarizes our conclusions.

THE MAP
We approximate the dynamical evolution of test particles on highly eccentric orbits using a twist map.We consider a test particle with semi-major axis  and pericenter distance  on a highly eccentric orbit about a star of mass  * .The test particle is subject to perturbations from a planet of mass   ≪  * on a circular orbit of radius   ; we assume that  ≫   and  ∼   .The test particle is assumed to orbit in the same plane as the planet.We define the dimensionless quantities  =   / * ,  =   / and  =   /; we refer to  as the "energy" since both  and  = − 1 2   * / are proportional to the inverse semi-major axis.Conservation of the Jacobi constant implies that, far from pericenter, the particle's Tisserand parameter,  ∝ 1 2  +  −1/2 (2 − /) 1/2 , is fixed.Since we are considering highly eccentric orbits with  ≫   and  ∼   , a constant Tisserand parameter implies that during a pericenter passage the changes in pericenter, , and energy  are related by |/| ∼ || ≪ 1.Thus we may treat the particle's pericenter distance, , as fixed.
Each time the test particle passes through pericenter, its interactions with the perturber result in an energy "kick", moving it onto a new orbit with semi-major axis  ′ and energy  ′ =   / ′ .As we show in Appendix A, this kick can be described by a function,   (), such that  ′ =  − 2  ()/ where  is the angular separation between the test particle and perturber when the test particle is at pericenter.The test particle returns to pericenter after a time  ′ where 2/ ′ =    ′3/2 and   is the mean motion of the planet.At this return, the angular separation between the test particle and perturber is  ′ =  +    ′ .Therefore, we have the following iterative map describing the evolution of the test-particle orbit at successive pericenter passages: Equations ( 1) and ( 2) constitute a two-dimensional area-preserving twist map depending on two parameters,  and .
Previous studies have approximated the kick function,   /, appearing in Equation ( 1) by means of series expansions (e.g., Petrosky & Broucke 1988), the interpolation of numerical integrals (e.g.Pan & Sari 2004;Hadden et al. 2018), or fits to direct -body integrations (Malyshkin & Tremaine 1999).We present a novel derivation of the Fourier series expansion of   () in Appendix A. The advantage of this Fourier series representation of the kick function is that it enables us to predict the onset of dynamical chaos and chaotic diffusion rates as a function of pericenter distance ( §3.2).In particular, we show that, for non-crossing orbits ( < 1),   () can be written as the cosine series where exact expressions for the coefficients   () are given by Equation (A17) in Appendix A. Numerical implementations of the mapping in Equations ( 1) and (2) require truncating the sum over index  in Equation (3) at a finite value.Generally, this is not a serious limitation because the coefficients decay exponentially with increasing .Moreover, in Appendix A we show that they admit an asymptotic approximation   () ≈ () −() / for constants  and  that depend only on  (eq.A29).We take advantage of this asymptotic behaviour in our numerical implementation of the map by using the identity where we have dropped an unimportant constant term,  0 () +1 2 ()().We choose  max to achieve a 5% or better fractional error in the amplitudes   .The asymptotic approximation   () ≈ () −() /, derived in Appendix A, is valid for  < 8/9, which corresponds to the limiting pericenter distance, as  → 0, reached by particles that are initially on circular orbits co-orbital with the perturber.While our focus in this paper is on non-crossing orbits, we note that Equation (4) suggests that the functional form of the kick function for particles on crossing orbits ( > 1) can be obtained by replacing the cosh  term in Equation ( 4) with cos  int where  int () is the orbital phase at which the particle's orbit intersects the planet's (see Eq. A13).Fitting the overall amplitude of the kick function derived in this manner may provide a simpler approximation of -body results than the piecewise polynomial interpolation used by Malyshkin & Tremaine (1999) to represent the kick function of crossing orbits.

Local approximation of map
Equations ( 1) and ( 2) provide an efficient yet accurate means of modeling the dynamics of comet-like orbits.Valuable insight into these dynamics is obtained by approximating the map locally in the vicinity of an :1 mean-motion resonance, that is, when the planet orbits  times for each orbit of the test particle.In order to do so, we will define  0 =  −2/3 as the nominal energy of the resonance and assume that variations about this nominal value, , remain small (|/ 0 | ≪ 1).Writing  =  0 1 − 2 3  , the map can be approximated locally as (5) in which we have omitted O ( 2 ) terms in the second equation, under the assumption that ||/ 0 ≪ 1.These equations again constitute a two-dimensional area-preserving twist map, with the dependence on both perturber mass and resonance index, , entering through a single parameter, The map variable  in Equations ( 5) and ( 6) has the following interpretation: when  = 0 the angle  advances by exactly 2( + ) with each pericenter passage of the test particle.This means that at integer values, i.e.,  =  with integer , the perturber completes  +  orbits before the test particle returns to pericenter, corresponding to an  +  : 1 mean-motion resonance (MMR).More generally,  = / with  and  integers corresponds to ( + ) :  MMRs.We will refer to such a resonance as a th order MMR.This differs from the standard nomenclature of celestial mechanics but, as noted by Pan & Sari (2004), is more appropriate when considering the dynamics of highly eccentric orbits with large semimajor axes.
The dynamics of the local map, Equations ( 5) and ( 6), are unaltered by adding or subtracting integer values to the variable  so that we can consider the dynamics by taking the value of  modulo 1. 1 This invariance with respect to the addition of integer values to  reflects a repeating pattern of MMRs between successive :1 MMRs.These resonances occur at rational values of  between 0 and 1.In particular, all of the resonances up to a maximum order, , occur at values of  given by the Farey sequence,   , i.e., the sequence of reduced fractions between 0 and 1 having denominators less than or equal to .This repeating pattern of resonances is analogous to the repeating pattern of resonances between adjacent  :  − 1 MMRs that are dynamically dominant for closely-spaced planet pairs as described in Hadden & Lithwick (2018).
Figure 1 shows a comparison between direct -body simulations and the local approximation of the map given by Equations ( 5) and (6).
-body integrations were done with the ias15 integrator (Rein & Spiegel 2015) of the REBOUND code (Rein & Liu 2012).The initial conditions are chosen near the orbit of TNO 148209 (2000 CR105) with a perihelion distance of  = 44.1 AU and near the 20:1 MMR with Neptune, which we model as a circular and coplanar perturber at semimajor axis   = 30 AU. Figure 1 can be compared with the surfaces of section shown in Figure 1 of Volk & Malhotra (2022).The -body surface of section in Figure 1 is constructed from a set of initial conditions with fixed Jacobi constant.Section points are generated by recording the test particle's osculating orbital period at aphelion and then the angular separation between the test particle and the perturber at its subsequent perihelion passage.

Dynamics of Resonances
The twist map in Equations ( 5) and ( 6) can be derived from equations of motion where the independent variable, , is the particle's mean anomaly, and is a 2-periodic delta function.Equations ( 8)-( 9) can be derived from the "Hamiltonian" function where (, ) are a canonically conjugate coordinate-momentum pair and the particle's mean anomaly, , plays the role of the independent "time" variable.
Let  0 = (/) −  denote the nominal value of  at the : MMR.Introducing new canonical variables  =  −  0 and  =  − (/) , dropping an unimportant constant, and only retaining resonant harmonics, the transformed Hamiltonian is Note that in the case of first-order resonances (i.e.,  = 1), the sum over  in Equation ( 11) is simply   ().
Figure 2 compares trajectories computed directly via the mapping, Equations ( 5) and (6), to level curves of Hamiltonians for various resonances given by Equation (11).The phase space topology of :1 MMRs, shown in the bottom row of Figure 2, is distinct from that of the higher-order MMRs shown in the top row of Figure 2. First-order :1 MMRs possess two unstable fixed points, at (, ) = (0, 0) and (, ) = (, 0), along with two stable "asymmetric" libration islands situated at  = 0 and  = ± asym () where  asym satisfies > 0. The value of  along the outer separatrices of first-order MMRs is given, as a function of , by In contrast with :1 MMRs, higher-order MMRs exhibit simpler phase space topologies, with each  th order MMR possessing a chain of  pendulum-like libration islands.Since typically   ≫   for  > 1, the dynamics of these MMRs can be approximated by retaining only the lowest harmonic from the sum appearing in Equation ( 11).With this approximation F res becomes the Hamiltonian of a pendulum and we obtain a resonance half-width, measuring the distance from the elliptic fixed point to the separatrix, of Δ  = 2 √︁  |  ()|/(2) for an : MMR.In terms of semi-major axis, the full width of a  th order MMR is approximately in which we have approximated the width of a first-order MMR by its width at  = .Note that the maximal width of first-order MMRs actually occurs at  =  asym (), but this width differs little from the width at  = , which yields a slightly simpler expression in Equation ( 13).

Onset of Chaos and Destruction of KAM Barriers
To predict the onset of large-scale chaos, we use the resonance covering fraction criterion developed in Hadden & Lithwick (2018) (see also Quillen 2011).This criterion predicts that resonance overlap and chaos occurs when the total width of all resonances in a given region of phase space exceeds the size of that region.Specifically, we follow Hadden & Lithwick (2018) and define a "resonance optical depth",  res , as the covering fraction of resonances in a unit interval of the map variable , and predict that chaos occurs when  res = 1.The sum of resonance widths in a unit interval of the map variable  is given by where (), Euler's totient function, counts the number of resonances of order  falling in a unit interval of  (see Hadden & Lithwick 2018). 2n Figure 3 we plot our prediction of the onset of chaos, where  res = 1, on a plot of semi-major axis versus perihelion distance for particles subject to a Neptune-mass perturber on a circular orbit at 30 AU.We compare our prediction to a grid of Lyapunov times computed via -body simulations that indicate where chaotic behaviour occurs.Each -body simulation was integrated for a total of 50 test particle orbits using the WHFast integrator of the REBOUND code (Rein & Tamayo 2015), an implementation of the Wisdom-Holman symplectic mapping algorithm (Wisdom & Holman 1991) that includes variational equations.Variational equations were integrated along each orbit and used to compute Lyapunov times, taken to be the inverse of the Lyapunov characteristic exponent computed via the Simulation.calculate_lyapunovmethod in REBOUND.
We compare our prediction for resonance overlap and the onset of chaos to the resonance overlap boundary derived by Batygin et al. (2021), shown as a blue dashed curve in Figure 3.The criterion derived by Batygin et al. (2021) accounts only for the overlap of :1 and :2 MMRs and the disagreement between their prediction and the numerical results, especially for small , underscores the increasing dynamical importance of high-order  :  MMRs with  > 2 as  → 1. 3  The long-term evolution of chaotic trajectories will be qualitatively different depending on whether or not their motion is bounded by KAM tori, which serve as absolute barriers to transport in the phase space of our two-degree-of-freedom system.We conduct numerical experiments with the map to predict the location of the last KAM barrier in phase space as a function of pericenter distance.We adopt the following procedure to estimate the location of the last KAM curve: . Finally, we approximate the sum over the remaining terms as an integral using where Γ[ •, • ] is the incomplete Gamma function. 3The formulae that Batygin et al. ( 2021) derive for the widths of  :1 and  :2 MMRs also differ from those predicted by Equation ( 13).(i) We set a fixed pericenter distance and select an initial value  to determine the parameter  appearing in our local approximation of the map (Equation 7).
(iii) The trajectories are iterated up to a maximum of 3 × 10 7 / times (corresponding to a time of approximately 5 Gyr) or until any of the trajectories escapes the region 0 <  < 1.
(iv) If any trajectory escapes, we reduce  by 1 and repeat the experiment.Otherwise, we take  =   ×  2/3 as the semi-major axis of the last KAM curve.
The computed locations of the last KAM barriers are plotted in Figure 3.For large , the location of the last KAM barrier closely matches the boundary for the onset of chaos while at smaller pericenter distances there is an increasing gap between the onset of chaotic behaviour and the location of the last KAM curve, indicating a growing region of bounded chaos.
We plot the pericenter and semi-major axis distances of TNOs with 45 AU <  < 1000 AU and 30 AU <  < 65 AU obtained from a query of the Minor Planet Center (MPC) database. 4Figure 3 shows that essentially all TNOs with  ≳ 100 AU are within the predicted resonance overlap region.Furthermore, these objects are generally near or below the curve marking the last KAM barrier at a given pericenter distance and are therefore subject to large-scale migration in semi-major axis.(Additional dynamical effects arising from the other giant planets, Neptune's eccentricity, and inclination will generally render the KAM boundary plotted in Figure 3 a partial, rather than absolute, barrier to transport for objects slightly above the curve.)There is an apparent lack of low perihelion distance objects ( ≲ 40 AU) at large semi-major axes, despite a strong observational bias towards discovery of low perihelion orbits.We show below that the rate of chaotic diffusion exhibits a steep (i.e., exponential) dependence on pericenter distance (Fig. 6) and the relative excess of high-perihelion distance objects is likely attributable to the exponentially longer expected survival times at higher pericenter distance.Since quasi-circular orbits with  ≳ 40 AU are far from the resonance-overlap region and cannot reach large semi-major axes via chaotic scattering with Neptune the high-, large  objects presumably owe their present orbits to torques from stellar passages or the Galactic tide experienced while scattering to large  from lower values of  (Vokrouhlický et al. 2019), or perhaps to a rogue planet (Gladman & Chan 2006;Silsbee & Tremaine 2018).

STATISTICAL DESCRIPTION OF ENSEMBLE EVOLUTION
Here we consider an ensemble of chaotic particle orbits and seek a statistical description of their evolution.We begin in Section 4.1 with an analysis of chaotic diffusion in the local approximation of the map that was presented in Section 2.1.We derive the local map's chaotic diffusion rate as a function of the parameters  =   / and  (Equation 7).Then, in Section 4.2 we derive a Fokker-Planck description of the evolution of the energy of the scattering particles.While modeling comet energy evolution via a diffusion equation is not new (e.g., van Woerkom 1948;Yabushita 1980), past works have assumed that the rate of diffusion per pericenter passage is independent of energy.This approximation is valid as  → ∞ but, as we show in Section 4.1, breaks down for semi-major axes  ≲ 1000 AU where the majority of known TNOs reside.We also solve the Fokker-Planck equation with boundary conditions that more accurately capture the behaviour of scattered disk objects than those used in past works.We compare our Fokker-Planck description of orbital evolution to direct -body simulations in Section 4.3.Finally, in Section 4.4 we examine deviations from diffusive behaviour exhibited by the map that result from resonant sticking.

Diffusion in the local approximation
For chaotic trajectories, the evolution of the map variable , when viewed on sufficiently long time scales in regions of phase space where no KAM barriers are present, is well-approximated as a random walk and the mean square displacement of an initial ensemble of nearby trajectories will increase linearly in time.We define a normalized diffusion coefficient for this random walk (see Appendix B1 for intermediate steps) where  is the map parameter defined in Equation ( 7), angle brackets indicate averages over initial values of  and , i.e., ⟨  (, ) 2 is referred to as the "quasi-linear" approximation of the diffusion rate (see, e.g., Karney et al. 1982).
The terms ⟨ 0   ⟩ capture correlations between kicks  iterates apart and can be computed as follows: the functions   ( 0 ,  0 ) are periodic in both  0 and  0 and so can be written in terms of a Fourier expansion as it follows from the definition of the functions   ( 0 ,  0 ) that â (,) Therefore, each correction term can be expressed as the following sum over Fourier amplitudes: where ..denotes the complex conjugate of the preceding term.Cary & Meiss (1981) describe a procedure for expressing the correlation terms, ⟨ 0   ⟩, analytically for a broad class of maps that includes our Equations ( 5) and ( 6).Using their procedure, it is straightforward to show that ⟨ 0  1 ⟩ = 0.However, the analytic expressions for higherorder correlations are unduly cumbersome due to the fact that the function     () appearing in Equation ( 5) has a non-zero Fourier amplitude at every harmonic.Therefore, we implement a numerical procedure to compute correlation terms ⟨ 0   ⟩ up to some maximum order, , and thereby obtain an estimate of the diffusion coefficient in Equation ( 15), which we denote as To calculate D  (, ), we compute  iterates of the map on a grid of initial ( 0 ,  0 ) values, recording   ( 0 ,  0 ) for  = 1, . . .,  at each grid point.Using the grid of computed   ( 0 ,  0 ) values, amplitudes â (,)

𝑘
are computed by means of a fast Fourier transform (FFT) algorithm.The resulting amplitudes are then used to compute the sum in Equation ( 16).These are in turn used to compute D  (, ) according to Equation ( 15) with the sum over  limited to  ≤ .
Figure 4 illustrates the behaviour of D  (, ) as a function of  for orbits with pericenter distance  = 40 AU subject to a Neptune-like perturber, at four different semi-major axes.The semi-major axes range from  ∼ 240 AU, just after the last KAM barrier has disappeared, to  ∼ 900 AU, corresponding to a regime of strong resonance overlap (see Figure 3).As the semi-major axis decreases, regular trajectories associated with low-order MMRs occupy progressively more phase-space volume.At the 23:1 MMR ( = 243 AU), where regular resonant islands take up the largest fraction of the phase space, the diffusion coefficient estimates D  (, ) show large oscillations about the numerically measured value, even at large .These oscillations arise from the contributions of the regular islands to the ensemble-averaged correlations, ⟨ 0   ⟩, a generic feature of force correlations in maps with a mixed phase space (e.g., Karney et al. 1982).We are interested in the diffusion rate of trajectories in the connected chaotic component of the map, which we assume to be well-approximated by the mean value, ≈ 1  max  max =0 D  (, ), for sufficiently large  max such that the oscillating contributions to ⟨ 0   ⟩ from regular islands are averaged out.In Figure 5 we calculate mean diffusion coefficients over a range of pericenter and orbital distances, taking  max = 64.We compare these semi-analytic diffusion coefficients to numerical measurements derived by evolving ensembles of chaotic trajectories.The numerical measurements are obtained by initializing 1000 trajectories with  0 = 0 and  0 drawn from a Gaussian distribution centered on  0 = 1/2 with dispersion  = 10 −3 so that the initial ensemble resides near the unstable period-2 orbit associated with an 2 : 2 + 1 MMR.The trajectories of the ensemble are evolved for 5,000 iterations and the variance of the ensemble's  values is recorded at each iteration.The diffusion coefficient, D (, ), is then estimated from a linear fit of the variance versus iteration data.Both numerical and semi-analytic results show that the diffusion rate declines to zero as the location of the last KAM curve, indicated by stars for each pericenter distance in Figure 5, is approached.We also plot the values of D 4 (, ) for the different pericenter distances, shown as dashed lines in Figure 5, to highlight that good agreement with numerical results can be obtained with relatively few low-order corrections to the quasi-linear diffusion rate.However, we caution that the good agreement between D 4 (, ) and the numerical results is deceiving at semi-major axes near the last KAM curve: there, the values D  show large oscillations about the numerically-measured value for  > 4, similar to the case of the 23:1 MMR in Figure 4, so that the agreement is generally worse if more terms are added.Nonetheless, we find empirically that the numerically determined diffusion rates are nearly always well-approximated by D 4 (, ) except very near the semi-major axis of the last KAM curve.Thus, D 4 (, ) is useful as a relatively accurate estimate of the diffusion rate that can be calculated very efficiently.

Formulation of Fokker-Planck Equation
We now use our results on the diffusion rates of chaotic trajectories in the local version of the mapping to obtain a statistical description of the chaotic evolution of test particles in the restricted three-body problem.Consider the evolution of an ensemble of test particles.Let (, ) be the number density of particles with inverse semi-major axes in the range [,  + ) at time .We assume that the ensemble members are randomly distributed in orbital phase.Let  ( ′ |) ′ be the probability that a particle with inverse semi-major axis  transitions to an orbit with inverse semi-major axis in the range [ ′ ,  ′ +  ′ ) in one orbital period   of the perturber due to a gravitational interaction.Then (, ) In order to derive the Fokker-Planck equation, the integrand in the right-hand side of Equation ( 18) is written as ( − , ) ( +  − | − ), then expanded as a Taylor series in − to obtain Truncating the expansion at second order5 and integrating by parts, we obtain where are referred to as the drift and diffusion coefficients, respectively.In Section 4.1, we showed that the mean square kick in  = 3 2  (1 − / 0 ) per perihelion passage is given by  2 D (, ).Noting that  = 3 5/3 ,  0 =  −2/3 , and that a particle's orbital period is equal to    −3/2 , this translates to a mean square kick in  per perturber orbital period of  () = 4 −1   2  3/2 D (3 −5/2 , ), which is precisely the diffusion coefficient appearing in Equation ( 20).Deriving the drift coefficient, (), nominally requires working out the expected kick to second order in , whereas our map has only approximated the kick function to first order in .However, we can take advantage of the fact that the underlying equations governing the particles' dynamics are Hamiltonian to derive a relationship between () and  ().We provide such a derivation, working directly from Hamilton's equations, in Appendix B2.Here we provide a shorter alternative derivation based on the concept of detailed balance.
Since the Hamiltonian flow must preserve phase space volume, a distribution that is uniform in both the action Λ ∝ 1/ √  and its associated angle variable  is stationary.So (, ) ∝  −3/2  is a stationary distribution.Equation (20) says that the current in  is ()(, ) − 1 2    ()(, ).Since detailed balance means that for an equilibrium stationary distribution there are no net probability currents, we obtain The resulting Fokker-Planck equation can then be written as where   =   4 2 D QL and  () = 3 −5/2 .Figure 6 shows the diffusion timescale,   , as a function of pericenter distance for a Neptune-like perturber.The figure illustrates that the diffusion timescale is a steep function of pericenter distance, and for  ≳ 40 AU, is well approximated by the exponential relationship where we have normalized the scaling relationship Equation ( 25) such that the nominal parameters correspond to Neptune's mass and orbital distance.
When the diffusion coefficient in Equation ( 24) is approximated by the quasi-linear value, i.e., when the  dependence of D ( (), ) is ignored, then the Fokker-Planck equation reduces to the diffusion equation modeled by Yabushita (1980).He derives the analytic solution where  1 is a modified Bessel function, as the Green's function for Equation ( 24) assuming initial conditions (, 0;  0 ) = ( −  0 ) and subject to the boundary conditions  QL (0, ;  0 ) = lim →∞  QL (, ;  0 ) = 0.There are two corrections to the diffusion model of Yabushita (1980) that we apply in order to capture the statistical evolution of test particles using the Fokker-Planck model in Equation ( 24).First, the existence of a last KAM curve at some sufficiently large energy,  KAM , implies that a more appropriate boundary condition when solving Equation ( 24) is to impose a reflecting boundary satisfying    3/2 (, ) Second, as we showed in Section 4.1, the diffusion coefficient D (, ) can deviate significantly from the quasi-linear diffusion rate near this upper boundary.Below in Section 4.3, we will take these deviations into account by fitting a simple parametric form to the calculated local diffusion coefficient, D (, ).

Comparison to N-body integrations
We compare our Fokker-Planck model and the iterated mapping to direct -body integrations in Figures 7 and 8.The -body integrations follow an ensemble of 4 × 10 4 test particles with initial pericenter distance  = 35 AU and semi-major axes  = 400 AU, perturbed by a Neptune-mass planet on a circular orbit at 30 AU about a 1 ⊙ star.Integrations are done with the REBOUND WHFast integrator (Rein & Tamayo 2015) using a time step equal to 1/20th of Neptune's orbital period.We also evolve an ensemble of 10 4 map orbits, initialized with  = 30/400 and random initial  values.For these orbits, we augment the mapping Equations ( 1) and ( 2) with an auxillary time variable, , updated at each step according to  ′ =  +    −3/2 .
Figure 7 shows the normalized histograms of test particle energies, , for -body integrations and map orbits at different times.The histograms show mostly excellent agreement, except that a handful of -body particles reach  values between the 4:1 MMR (indicated by the dashed vertical line) and the 3:1 MMR, whereas the map orbits are excluded from this region.The disagreement between the mapping and -body integrations regarding the behaviour near the 4:1 MMR can be attributed to the mapping model's assumption of a constant pericenter distance.Conservation of the Tisserand parameter implies that particles initialized with a pericenter distance of  = 35 AU and  = 400 AU will have a slightly smaller pericenter distance of  ≈ 34.5 AU upon reaching the 4:1 MMR.We find by conducting the test to determine the last KAM barrier described in Section 3.2, that at this slightly lower pericenter distance, there is no longer a bounding KAM at the 4:1 MMR, though the flux of particles across the resonance is extremely slow.(It requires on the order of ∼ 5 × 10 4 map iterations, corresponding to ∼ 30 Myr, for the first particle from an ensemble of 10 3 particles initialized in the vicinity of (, ) = 0.5 to escape the region 0 <  < 1.)No -body particles reach period ratios interior to the 3:1 MMR, consistent with additional numerical tests with the mapping that indicate a bounding KAM curve exists near this MMR when  = 34.5 AU.
Figure 7 also compares the predictions of two Fokker-Planck models to the particle energy distributions obtained with the mapping and -body simulations.For both models, we impose an absorbing boundary at  = 0 and a reflecting boundary at  KAM = 4 −2/3 .The first Fokker-Planck model, which we refer to as the "quasi-linear" model, approximates the diffusion coefficient as  () = 4 3/2 D QL (/  ) and the corresponding density distribution is predicted by Equation (30).For the second model, which we will refer to as the "numerical" model, we adopt a simple parametric model for the diffusion coefficient given by where  () = 3 −5/2 and the parameter values ( 1 , ,  2 ,  * , Δ * ) = (0.136573, 0.971338, 0.286504, 0.269842, 0.170774) are obtained from a fit to the data shown in Figure 5. Equation ( 24) is then solved numerically to obtain the density distribution, (, ), at the times shown in Figure 7.The consequences of the reduced diffusion rate at values of  ≳ 0.2 are especially evident at 100 Myr, where the numerical model predicts a much smaller fraction of particles at  ≳ 0.2 than the quasi-linear model.The numerical model also predicts a significantly larger fraction of particles surviving in this region at late times compared to the quasi-linear model.Both Fokker-Planck models significantly under-predict the number of very loosely bound particles evident in the smallest histogram bins.In this region, the approximation of particle orbits as a time-continuous random process, which underlies the derivation of the Fokker-Planck equation, breaks down, as we now show.Figure 8 plots the fraction of surviving particles versus time for both the mapping and -body simulations and compares these to the prediction of Equation (31).The mapping and -body integrations show good agreement, while the Fokker-Planck model severely underpredicts the fraction of bound particles at late times.We find that the disagreement is due to the fact that at  ≪ 1 the approximation of trajectories as a time-continuous random walk breaks down.To demonstrate this we simulate particle trajectories, {( 1 ,  1 ), ( 2 ,  2 ), ...}, as a discrete random walk, where ( 1 ,  1 ) = (  (/  ) 3/2 ,   /) and with Δ  ∼ N (0, 1) a zero-mean Gaussian random variable with unit variance.The random walks are stopped at the first step, , for which   < 0. We simulate 2 17 random walks and plot the fraction of surviving particles as a function of time in Figure 8.The discrete random walk shows excellent agreement with the survival fractions of the mapping and -body simulations.

Resonance sticking
The phenomenon of "resonance sticking" can cause the evolution of chaotic trajectories to deviate from a random-walk process and lead to power-law tails in the distribution of times needed for trajectories to escape from regions of phase space containing a mixture of regular and chaotic trajectories (e.g., Meiss 2015).In their study of evolution of comets on Neptune-crossing orbits, Malyshkin & Tremaine (1999) find that the population of comets surviving up to ∼ 5 Gyr is dominated by orbits that experience prolonged periods of resonance sticking.By contrast, the agreement between the survival times of the discrete random-walk model, which cannot produce sticking, and the -body integrations presented in Figure 8 indicates that resonance sticking plays a minimal role in the survival of particles up to 5 Gyr for the initial semi-major axis and pericenter distance adopted in these experiments.The influence of sticking is not apparent in the simulations of Malyshkin & Tremaine (1999) until the comets are depleted by a factor of ∼ 10 −4 , whereas the total particle population in Figure 8 is only depleted by a factor of ∼ 0.05 after 5 Gyr, so it is not entirely surprising that sticking has a negligible influence on our survival statistics.Nonetheless, we find that resonance sticking does occur in our map and we briefly explore its effect on the escape times from local regions of phase space below.
In Figure 9, we use the local map to compute the fraction of orbits remaining in the vicinity of :1 MMRs with 11 ≤  ≤ 16 as a function of  31), while the solid blue curve shows the results of a numerical integration of the Fokker-Planck equation using the diffusion coefficient given by Equation ( 32).The red dotted line shows the results from the discrete random walk (Equation 34).
time for a range of pericenter distances.Here and throughout this section, all map parameters are reported assuming a Neptune-mass perturber located at an orbital distance of   = 30 AU.For each :1 MMR and pericenter distance, we initialize an ensemble of 655,360 orbits near the unstable fixed point of the first-order resonance at (, ) = (0, 0) by drawing their initial  values from a Gaussian centered at  = 0 with dispersion   = 10 −5 while setting their initial  = 0.Each orbit is then evolved until it escapes the region −1 <  < 1.The fraction of surviving orbits,  survive , is plotted against time in Figure 9 (here time is taken to be the number of map iterations times  ×   where  is set by the particular MMR).
For  survive ≳ 0.01, all of the curves in Figure 9 exhibit the exponentially decaying behaviour expected for a normal diffusive process.This is illustrated clearly in Figure 10, where we re-plot the survival fraction curves shown in Figure 9, but with the time coordinate re-scaled by a best-fit exponential decay time scale,   .(We determine the decay time scale,   , by means of a least-squares fit of each  survive () curve shown in Figure 9 using a simple exponential decay model,  −/  .)For  survive ≲ 0.01, many of the curves exhibit apparent algebraic decays at late times, indicative of resonance sticking (e.g., Karney 1983).As is evident in Figure 10, there does not appear to be a universal power-law slope for the decay of survival fractions, which instead ranges approximately between  −2 and  −1 behaviour.
Closer inspection reveals that changes in the late-time survival fraction behaviour at a given pericenter distance are associated with changes in the phase-space topology of the :1 MMRs as  varies.This is illustrated in Figure 11, where we plot a series of survival fraction curves for  = 35.5 AU at different :1 MMRs accompanied by two-dimensional histograms showing the densities of the trajectory points of the longest-surviving orbits.Specifically, the histograms are constructed by recording the trajectories of 32,768 map points initialized near (, ) = (0, 0) and binning the trajectory points of the 33 orbits (= 0.1% of 32, 768) that remain in the region −1 <  < 1 the longest.To highlight the contribution of sticking in :1 MMRs to the survival fraction behaviour, we compare survival fractions versus time in two different regions near each :1 resonance in the top panels of Figures 11.The first region we consider includes the range −1 <  < 1, as in Figure 9, with trajectories initialized near the :1 MMR located at  = 0.The second region includes only the range 0 <  < 1, and we initialize trajectories near the (2 + 1):2 MMR, with points near (, ) = (0, 1/2).Since neighboring :1 MMRs sit on the boundary of this second region, orbits that stick to the boundary of these MMRs quickly exit the region 0 <  < 1, and so sticking to :1 MMRs has negligible influence on the escape times for this second region.The resulting survival fractions for the first and second regions are plotted as solid and   that :1 MMRs produce the most dramatic sticking effects.Furthermore, it appears that the :1 MMRs are "stickiest" at or very near the value of  at which chaotic orbits can first traverse the region between the two asymmetric libration islands on either side of the hyperbolic fixed point at (, ) = (, 0).We observe similar behaviour across all pericenter distances in our numerical experiments: sticking in :1 MMRs is most significant near the values of  at which the last invariant symmetric libration curve (i.e., the last regular :1 resonant trajectory in which the resonance angle librates symmetrically about ) vanishes.

SUMMARY AND CONCLUSIONS
We have presented a new iterated map for modeling the dynamical evolution of highly eccentric test particles subject to a planetary mass perturber.Utilizing a new formula for the energy kick received by test particles at each perihelion passage, the map extends the validity of previous mapping treatments of this problem to semi-major axes and pericenter distances relevant to the bulk of the observed extreme TNO population.The map shows excellent agreement with direct -body integrations, at a dramatically reduced computational cost, while also elucidating how the dynamics depend on perturber mass and orbital period through the local approximation introduced in Section 2.1.In Section 3.1, we showed how the new kick function can be used to derive an integrable Hamiltonian model for the dynamics of distant mean-motion resonances (MMRs).These models are well-suited to study the resonant dynamics of highly eccentric, distant TNOs and do not require resorting to the usual expansions in powers of semi-major axis ratio or eccentricities.Our Hamiltonian model for MMRs also allows us to predict resonance widths and to determine, for a given pericenter distance and perturber mass, the semi-major axis at which resonance overlap leads to the onset of dynamical chaos.We provide a novel analytic criterion for the onset of chaos, given by Equation ( 14), that agrees well with numerical simulations.This new criterion improves upon the resonance overlap criterion of Batygin et al. (2021), which approximates the onset of chaos at large semi-major axis but performs poorly at smaller semi-major axes more typical of the observed TNO population.
In Section 4, we have presented a statistical description of the transport of chaotic orbits.We have shown how the mapping can be used to determine the local energy diffusion rate in Section 4.1.We demonstrate that, in regions where regular islands occupy a non-negligible fraction of phase space, long-time correlations can significantly reduce the diffusion rate relative to a naive prediction that treats successive energy kicks during pericenter passage as uncorrelated.Examining the diffusion rates shown in Figure 5 at semi-major axes and pericenter distances typical of the observed TNO population shown in Figure 3, it is clear that the effects of correlations can play a significant role in the chaotic transport of much of this population.
The local diffusion rates derived from the map in Section 4.1 serve as the basis for a Fokker-Planck description of the global transport of the two ensembles from the regions −1 <  < 1 and 0 <  < 1, respectively and found them to be nearly identical after multiplying the escape times of the second ensemble by a factor of 4.
particles in energy space, derived in Section 4.2.Our Fokker-Planck model reduces to the diffusion equation treated in previous works (e.g., van Woerkom 1948;Shteins & Riekstyn'sh 1961;Yabushita 1980) when one ignores deviations from the quasi-linear value of the diffusion coefficient introduced by correlations and approximates it as  () ∝  3/2 .Note that the mapping enables us to compute the value of the diffusion coefficient explicitly, whereas past works have either relied on order-of-magnitude estimates or simply worked in re-scaled time units so that the diffusion constant is unity.Additionally, the existence of KAM tori at short orbital periods for a given pericenter distance (or, more precisely, Jacobi constant) results in a reflecting boundary condition at some orbital energy,  KAM , rather than the absorbing boundary at  → ∞ used by Yabushita (1980).We derived a new analytic solution (Equation 30) to the Fokker-Planck problem under these boundary conditions, assuming the quasi-linear diffusion approximation.We demonstrated in Section 4.3 that our simple Fokker-Planck formulation provides good agreement with direct -body integrations.Our comparisons in Figure 7 reveal that deviations of the diffusion rate from its quasi-linear value can have a large influence on the number of particles at small semi-major axes (i.e., large ), where they are most observable.We also showed that the discrete nature of the random walk executed by particles, which maintain essentially constant orbital energies between successive pericenter passages, enhances the number of very loosely bound particles relative to the predictions of the Fokker-Planck model, which assumes a continuous-time random walk.
Finally, in Section 4.4, we showed that sticking in MMRs leads to late-time algebraic decay in the fraction of particles remaining in local regions of phase space.First-order :1 MMRs generally produce the most significant sticking effects, especially near the orbital period at which the last regular symmetric libration trajectory vanishes for a given pericenter distance.An interesting avenue of future investigation would be to explain the sticking behaviour around :1 MMRs and the late-time algebraic decay of the survival-time distributions in terms of a Markov tree model for transport (e.g., Meiss & Ott 1986;Alus et al. 2014Alus et al. , 2017)), especially near critical parameter values at which the last invariant symmetric libration curve breaks up.
There are a number of prospects for refining the dynamical model we have developed here.First, we have only considered a circular and coplanar perturber.It is relatively straightforward to extend the perturbative treatment in Appendix A to derive the increments in inclination, ascending node, and argument of pericenter induced by an eccentric, inclined perturber with each perihelion passage.From these, it should be possible to derive higher-dimensional symplectic mappings to describe the orbital evolution of inclined, eccentric test particles subject to an eccentric, massive perturber.Extending the mapping approach to include inclination dynamics may be particularly important for modelling the chaotic transport of eTNOs: the -body simulations of Gallardo et al. (2012) show that particles' chaotic semi-major axis diffusion can be intermittently interrupted by temporary capture into MMRs, during which time the inclination and pericenter distance can be raised by a Kozai-like resonance (see their Figure 14).Kozai cycling while sticking to MMRs has also been proposed as a mechanism for stranding inclined, detached TNOs at high pericenter distances during the late stages Neptune's outward migration (Gomes 2003).The map could also be augmented to capture the influence of additional interior perturbing planets, secular torques from the galactic tides, and stochastic perturbations from passing stars.and Equation (A24), for  > 0 we can write The integral on the right-hand side of Equation (A26) can be approximated for large  using the method of steepest descent.To do so, define In the limit  → 0, we get  → 1,  → 2 √ 2 3 3/2 − 1 + 1 2 log(/2), and  → 2 3/4  1/4 .The point,  * (), at which the function   exhibits a saddle occurs on the real axis for  = 8/9 so that the steepest descent method used to approximate the integral in Equation (A27) is no longer valid for  ≥ 8/9.Interestingly, this value of  is the limiting value, as  → 0, reached by particles with a Tisserand parameter of  = 3 (corresponding to an orbit that coincides with the planetary orbit; see, e.g., Pan & Sari 2004).In other words, particles with larger limiting values of  can achieve orbits that cross the perturber's for some positive value of .
Since we are orbit-averaging, we can integrate by parts over .Then the first and fourth term cancel, as do the second and third.This proves that The drift and diffusion coefficients   and   are the mean and mean-square changes in the action per unit time.If the action is related to some other variable  through  =  () then the mean and mean-square changes in ,  and , are given by   =  ′ () + 1 2  ′′ (),  2  = [  ′ ()] 2 , (B14) so Equation (B13) becomes If  is the inverse semi-major axis then  () ∝  −1/2 so  ′ () ∝  −3/2 and we recover Equation ( 23).
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Comparison of the local map equations, Equations (5) and (6), and direct  -body simulation.The map parameters are  = 5.15 × 10 −5 ,  = 20, and  = 0.68, corresponding to a perihelion distance  = 44.1 AU for a planet at Neptune's semi-major axis.The  -body simulations follow particles near the 20:1 MMR with a perturber having Neptune's mass and semi-major axis.The vertical axis is the ratio of the orbital period to Neptune's orbital period.

Figure 2 .
Figure 2. Comparison of resonant trajectories computed with the mapping to level curves of the analytic resonance Hamiltonian, Equation (11), for resonances up to order  = 5.The value of the analytic Hamiltonian is obtained as a function of  and  by substituting ( , ) = (  ,  − / ) in Equation (11).The parameters of the map are  = 0.005 and  = 0.75.For the case of particles subject to perturbations from Neptune, these parameters correspond to orbits with pericenter distances of  ≈ 40 AU near the 8 : 1 MMR.

Figure 3 .
Figure 3. Dynamical map in semi-major axis, , versus perihelion distance, , for test particles perturbed by Neptune (/ ⊙ = 5.15 × 10 −5 ,  = 30 AU).The color scale shows the value of the Lyapunov time computed on a grid of  -body simulations.The red line shows the analytic prediction of the resonance overlap criterion, setting  res = 1 in Equation (14).The blue line shows the analytic resonance overlap criterion ofBatygin et al. (2021).Observed low-inclination TNOs ( < 30 • ) are plotted as cyan points.See text for additional details.

Figure 4 .
Figure 4.The left column shows trajectories of individual chaotic orbits of the local map, Equations (5) and (6), in the vicinity of various  : 1 MMRs of a Neptune-mass perturber ( = 5.15 × 10 −5 ,  = 30 AU) for a pericenter distance of  = 40 AU (i.e.,  = / = 0.75).The right panel shows the fractional difference between the diffusion rate D  (  , ) defined in Equation (17) and the quasi-linear diffusion estimate D QL , as a function of the number, , of correlation terms, ⟨ 0   ⟩, included in the calculation of the diffusion coefficient.Dashed horizontal lines show diffusion rates measured numerically from ensembles of map trajectories, as described in the text.

Figure 5 .
Figure 5.Comparison of the numerically measured diffusion rate of a ball of initially close trajectories and our semi-analytic method for predicting diffusion rates using Equation (17) (see text for details).Diffusion rates normalised by the quasi-linear rate are plotted as a function of semi-major axis for different pericenter distances.An arbitrary offset is included for each pericenter distance to improve legibility.Numerically estimated values are plotted as open squares.Lines show our semi-analytic estimates of the diffusion coefficient, computed as either D (  , ) ≈ 1 64

Figure 7 .
Figure 7.Comparison between particle orbital energy distributions at different times predicted by direct  -body simulations, the iterated mapping, and the Fokker-Planck equation, for an ensemble of test particles subject to a Neptune-mass perturber at   = 30 AU. Particles are initialized with  = 400 AU and a pericenter distance of  = 35 AU.Histograms in each panel show the number of particles, normalized by the initial number of particles,  0 , in bins of  =   /.The predictions of the Fokker-Planck models are plotted as blue lines.The dashed blue line shows the analytic prediction of Equation (30), obtained by approximating the diffusion coefficient as  (  ) = 4 2  3/2 D  (  /), while the solid blue line shows a numerical solution to the Fokker-Planck equation adopting the diffusion coefficient given in Equation (32).The Fokker-Planck equations are solved with a reflecting boundary at    = 4 −2/3 , indicated by the vertical dashed line in each panel.

Figure 8 .
Figure 8. Fraction of surviving particles versus time for an ensemble of particles having initial semimajor axis  = 400 AU and pericenter distance  = 35 AU and subject to a Neptune-mass perturber on a circular orbit at 30 AU.The black squares show  -body results while cyan crosses indicate results from the mapping.The blue dashed line shows the prediction of the quasi-linear Fokker-Planck model, given by Equation (31), while the solid blue curve shows the results of a numerical integration of the Fokker-Planck equation using the diffusion coefficient given by Equation (32).The red dotted line shows the results from the discrete random walk (Equation34).

Figure 9 .
Figure 9. Survival fraction versus time for chaotic orbits escaping the vicinity of  :1 MMRs of Neptune at different pericenter distances.The fraction of orbits,  survive , from an initial ensemble of 655, 360 orbits, that remain within −1 <  < 1 is plotted against time (computed as    times the number of iterations).Different color curves show different pericenter distances.

Figure 10 .
Figure 10.Survival fraction versus time normalized by a best-fit decay time scale,   .Each curve shows the survival fraction for one of the parameter combinations of  and pericenter distance, , plotted in Figure 9. Power-law slopes of −1, −3/2 and −2 are indicated by colored lines for comparison.

Figure 11 .
Figure11.Top panels show the surviving fraction of orbits in the neighborhoods of various  :1 MMRs as a function of time.To highlight the contribution of resonance sticking in  :1 MMRs, the solid lines show the surviving fraction of orbits initialized near  = 0 that remain in the region −1 <  < 1, while the dashed lines show the surviving fraction of orbits initialized near  = 0.5 that remain in the region 0 <  < 1. (Exit time values for the latter curves are multiplied by a factor of 4 in order to account for the smaller region size.)The bottom panels show two-dimensional histograms constructed from the 100,000 iterations of the 0.1% of trajectories from an initial ensemble of 32, 768 initialized near  = 0 that remain in the region −1 <  < 1 the longest.The  values of these trajectories are plotted modulo 1 in the histograms.

Figure A1 .
Figure A1.Contour plot showing the real part of the function   () appearing in the integral defined by Eq. (A27) for  = 1/4.The red line indicates a branch discontinuity.The contour C is indicated by the blue line, which passes through the saddle point,  * (), indicated by a black 'x', along the path of steepest descent.