Improving the accuracy of simulated chaotic $N$-body orbits using smoothness

TL;DR

This paper demonstrates that increasing the smoothness of symplectic integrators in N-body simulations significantly improves long-term accuracy, highlighting the importance of Hamiltonian smoothness in dynamical studies.

Contribution

It introduces the idea that preserving Hamiltonian smoothness in hybrid integrators enhances simulation accuracy without extra computational cost.

Findings

Smooth hybrid integrators reduce Jacobi constant error by about 5 orders of magnitude.

Increasing smoothness improves long-term stability of chaotic orbits.

Smoothness of the N-body Hamiltonian is a crucial property for accurate simulations.

Abstract

Symplectic integrators are a foundation to the study of dynamical $N$-body phenomena, at scales ranging from from planetary to cosmological. These integrators preserve the Poincar\'e invariants of Hamiltonian dynamics. The $N$-body Hamiltonian has another, perhaps overlooked, symmetry: it is smooth, or, in other words, it has infinite differentiability class order (DCO) for particle separations greater than $0$. Popular symplectic integrators, such as hybrid methods or block adaptive stepping methods do not come from smooth Hamiltonians and it is perhaps unclear whether they should. We investigate the importance of this symmetry by considering hybrid integrators, whose DCO can be tuned easily. Hybrid methods are smooth, except at a finite number of phase space points. We study chaotic planetary orbits in a test considered by Wisdom. We find that increasing smoothness, at negligible extra computational cost in particular tests, improves the Jacobi constant error of the orbits by about $5$ orders of magnitude in long-term simulations. The results from this work suggest that smoothness of the $N$-body Hamiltonian is a property worth preserving in simulations.

Equations

$$A$$

$$B$$

$$T_1$$

$$V_{p1}=-\sum _{0<i<j}\frac{G{m}_i{m}_j}{Q_{ij}}\left(1-K(Q_{ij})\right),$$

$$T_0=\frac{1}{2m_0}{\left(\sum _{i\ne 0}{\mathbf{P}}_i\right)}^2,\text{and}{V}_{p2}=-\sum _{0<i<j}\frac{G{m}_i{m}_j}{Q_{ij}}K(Q_{ij}).$$

$${\mathrm{e}}^{h\hat{\tilde{H}}}$$

$${\tilde{H}}^{\mathrm{BAB}}$$

$${\tilde{H}}^{\mathrm{ABA}}$$

$$\{\{A,B\},B\}$$

$$\{\{B,A\},A\}$$

$$-\{\{T_1,V_{\mathrm{pl}}\},V_{\odot }\}-\{\{V_{\odot },T_0\},T_1\}$$

$$\{\{A,B\},B\}$$

$$\{\{B,A\},A\}$$

$${\dot{\mathbf{Q}}}_i$$

$${\dot{\mathbf{V}}}_i$$

$$\frac{\partial K(Q_{ij})}{\partial {\mathbf{Q}}_{ij}}={\frac{{\mathbf{Q}}_{ij}}{Q_{ij}}\frac{\partial K(r)}{\partial r}|}_{r=Q_{ij}}=\frac{{\mathbf{Q}}_{ij}}{Q_{ij}}{K}^{\prime }(Q_{ij}).$$

$${\dot{\mathbf{Q}}}_i$$

$${\dot{\mathbf{V}}}_i$$

$${\mathrm{e}}^{h\hat{\tilde{H}}}={\mathrm{e}}^{h\hat{A}}\circ {\mathrm{e}}^{h\hat{B}}.$$

$$\tilde{H}=T_1+T_0+V_{\text{sun}}+V_{\mathrm{pl}}+V_{p2}+\frac{h}{2}\left(\{V_{p2},T_1\}+\{T_0,V_{\text{sun}}\}\right)+\mathcal{O}(h^2).$$

$${\dot{\mathbf{Q}}}_i$$

$${\dot{\mathbf{V}}}_i$$

$${\mathbf{a}}_{0i}$$

$${\mathbf{a}}_{1i}$$

$${\mathbf{b}}_{0i}$$

$${\mathbf{b}}_{1i}$$

$$+1/3Q_{ij}^2{K}^{\prime \prime }(Q_{ij}))-{\mathbf{V}}_{ij}{Q}_{ij}^2(K(Q_{ij})-Q_{ij}{K}^{\prime }(Q_{ij}))]$$

$$+\frac{G}{Q_i^5}\left[-3{\mathbf{Q}}_i({\mathbf{Q}}_i\cdot {\mathbf{P}}_{\text{sun}})+{\mathbf{P}}_{\text{sun}}{Q}_i^2\right]),$$

$$H=\frac{p^2+q^2}{2}.$$

$$H_1=\frac{q^{2}K(q)+p^2}{2}\text{and}{H}_2=\frac{q^2(1-K(q))}{2}.$$

$${\mathrm{e}}^{h\hat{\tilde{H}}}={\mathrm{e}}^{h\widehat{H_2}}\circ {\mathrm{e}}^{h\widehat{H_1}}.$$

$$\tilde{H}=H+\frac{h}{2}\left\{H_1,H_2\right\}+\mathcal{O}(h^2),$$

$$\{H_1,H_2\}=p\left(\frac{q^2}{2}\frac{dK}{dq}-q(1-K)\right).$$

$$\delta=\frac{\Delta H}{H_{0}}=\frac{-\frac{h}{2}\left(\{H_{1},H_{2}\}\bigg{\rvert}_{+}-\{H_{1},H_{2}\}\bigg{\rvert}_{-}\right)}{H_{0}}=sapq^{2}.$$

$$C_J=\frac{1}{2}{\mathbf{v}}^2+\Phi (\mathbf{x})-\mathbf{\omega }\cdot (\mathbf{x}\times \mathbf{v}),$$

$$\Phi (\mathbf{x})=-\frac{G{m}_1}{|\mathbf{x}-{\mathbf{x}}_1|}-\frac{G{m}_2}{|\mathbf{x}-{\mathbf{x}}_2|}.$$

$$\mathbf{x}$$

Full text

1 Introduction

Symplectic integrators have been successful at studying a wide range of $N$-body particle dynamics in astrophysics, from planetary systems to galaxies with dark matter. They are usually severely limited in that they require a timestep, but gravitation has no timescales associated with it. To combat this limitation, a number of techniques have been employed (Hernandez, 2019) such as block timestepping or hybrid integration, which have become the basis of popular dynamics codes. These modifications change the smoothness of the $N$-body integrator Hamiltonian. The $N$-body Hamiltonian is a $C^{\infty}$ function for the domain of separations greater than [math], meaning it can get differentiated infinite times with respect to phase space coordinates and still be continuous. It is also called smooth. Smoothness measures the number of continuous derivatives. Traditional symplectic integrators constructed through operator splitting (Yoshida, 1990; Wisdom & Holman, 1991) respect this symmetry of the $N$-body Hamiltonian and are obtained from $C^{\infty}$ Hamiltonians. However, hybrid integrators (Chambers, 1999; Rein et al., 2019; Hernandez, 2019; Wisdom, 2017) are derived from $C^{n}$ Hamiltonians, where $n$ is an integer satisfying $n\geq 0$. A $C^{n}$ function can be differentiated $n$ times and still be continuous. $n$ is denoted the differentiability class order (DCO). A $C^{1,1}$ function has a Lipschitz continuous, rather than continuous, derivative. Block time-stepping methods Farr & Bertschinger (2007) and integrator subcycling methods Springel (2005) are not derived from smooth Hamiltonians; they switch in analogy to the ‘modified Heaviside function’ from Rein et al. (2019). It remains to be seen whether Hamiltonian smoothness is a symmetry which integrators, whether they are symplectic or not, should mimic, just as they might mimic the conservation of Poincaré invariants, the existence of a global Hamiltonian (Hairer et al., 2006, Section IX.3.2), the time-reversibility (Hernandez & Bertschinger, 2018), or the angular and total momentum conservation. If an integrator is not Hamiltonian, we can still quantify its smoothness through its modified differential equations (e.g. Hairer et al., 2006; Hernandez & Bertschinger, 2018).

Hernandez (2019) has explored the question of integrator smoothness for regular orbits. The work considered hybrid integrators, because it is simple to modify the smoothness of the Hamiltonian from which they are derived, which we define as the smoothness of the integrator. Hybrid integrators are $C^{\infty}$ except at a finite number of position coordinates. Consistent with the predictions of Arnaud & Xue (2016), Hernandez (2019) found that, apart from caveats, smoothness of at least $C^{1,1}$ yielded periodic orbits. Increasing the smoothness past $C^{1,1}$ did not help bound the error.

The situation is different for chaotic orbits. Rein et al. (2019) studied the effect of Hamiltonian smoothness during a short-term chaotic simulation which integrates over discontinuities twice or four times, except in the modified-Heaviside case. The errors were studied using approximations to commutators which may or may not converge. Rein et al. (2019) argued that over this short-term simulation, Hamiltonian smoothness made no significant impact. However, the advantage of symplectic integrators is in their long-term behavior.

In this work, we study whether the smoothness of the $N$-body Hamiltonian matters for simulating long-term chaotic orbits. We consider hybrid integrators, as in Hernandez (2019), because their smoothness is easy to tune. We construct a hybrid $N$-body integrator combining ideas from Chambers (1999), Rein et al. (2019), and Wisdom (2017). According to the modified differential equations (MDEs) (e.g. Hernandez & Bertschinger, 2018; Hairer et al., 2006), we predict Hamiltonian smoothness should have a significant impact on the error of integrators. To test these ideas, we consider a restricted three-body problem in which the test particle executes a chaotic orbit in which it can visit both masses. We observe integrator performance improvement by increasing the DCO from $n=0$ to $n=4$. The error of the integrator apparently saturates for higher $n$. The MDEs suggest that the effects of nonsmoothness, for this particular hybrid integrator, should disappear as the timestep approaches 0. We verify this numerically by rerunning our test with a step $10$ times smaller. We also find that increasing the DCO improves the error in regular orbits. This work provides evidence that the smoothness of the $N$-body problem is a feature which should be followed as closely as possible by $N$-body integrators. This is not the case for the other methods described above such as block time-stepping.

In Section 2 we show how we construct the hybrid $N$-body integrators. We carry out error analysis of a similar hybrid integrator using MDEs in Section 3 and give an example of how a non-smooth function leads to secular error. Numerical experiments demonstrating the significance of Hamiltonian smoothness are carried out in Section 4. Section 5 concludes and discusses the impact of this work on $N$-body integrations.

2 Algorithm description

To construct a hybrid $N$-body integrator, we combine ideas from Chambers (1999) and Wisdom (2017). A coordinate system must be chosen. As pointed out by Chambers (1999), Jacobi coordinates are problematic if the planets have mass. Jacobi coordinates are appropriate if massive bodies cannot have close encounters, but massless test particles can have close encounters with the massive bodies. This was the situation studied by Wisdom (2017). We thus use Democratic Heliocentric coordinates (e.g. Hernandez & Dehnen, 2017), which accommodate both massive and massless possibly having close encounters. The Hamiltonian in these coordinates is split into two parts, $A$ and $B$, given by (Chambers, 1999):

[TABLE]

We defined,

[TABLE]

$K$ is a function such that $0\leq K\leq 1$. $A$ is Liouville integrable if all $K$ equal $1$. Then, a Kepler advancer such as that of Wisdom & Hernandez (2015) can be used to solve it. Otherwise, $A$, is not described by periodic orbits, and numerical approximation methods, such as those relying on computing Taylor series solutions to the ordinary differential equations (ODEs) to high powers in an expansion parameter must be used. Such a method is that of Rein & Spiegel (2015), which uses substeps to calculate series coefficients. $B$ is always integrable and simple to solve. $A$ and $B$ both conserve linear and angular momentum, which we verified numerically.

We construct an integrator using

[TABLE]

where $h$ is a timestep. Chambers (1999) and Rein et al. (2019) used the first version, $BAB$. We call the second version $ABA$. The $n=1$ planet problem (the Kepler problem) is not solved exactly by these integrators. To see this, let the maximum $i$ in Eqs. (1) and (2) be $1$. This problem has a solution, involving modifying the terms placed in $A$ and $B$. (e.g. Hernandez & Dehnen, 2017, Section 4.2) showed how to do this when all $K=0$, but the generalization to allow other $K$ is immediate. We tested this alternative splitting, but it did not perform better in the tests in this paper because solving the $n=1$ problem is not the only consideration for a good integrator (Hernandez & Dehnen, 2017, Figure 6, and error analysis). $\tilde{H}$ is given by the BCH formula (e.g. Hairer et al., 2006) as,

[TABLE]

$\{\}$ are Poisson brackets (Goldstein et al., 2002, Section 9.5). If all $K$ are 1, the Poisson bracket terms are given by Hernandez & Dehnen (2017):

[TABLE]

$\epsilon=m/m_{0}$, where the largest planetary mass $m$ gives upper error bounds. $\epsilon\ll 1$. If all $K=0$, the Poisson brackets are,

[TABLE]

Thus, the $ABA$ form has a lower error in the two limits. We performed numerical experiments of this prediction by integrating a system composed of the Sun, Jupiter, and Saturn. We tested three cases: $K$ was always 0, $K$ was always $1$, and $K$ allowed to enter $0<K<1$, using a $C^{2}$ (see Section 4) $K$. We used a variety of stepsizes and integration times. In all cases, the absolute magnitude of the energy error was about twice smaller in the $ABA$ case, as compared to $BAB$. This is consistent with (4) and (5), where the magnitude of the $\{\{B,A\},A\}$ coefficient is twice smaller for $ABA$. $ABA$ is approximately just as expensive as $BAB$ because evolving the final $A$ from a step can be combined with evolving the first $A$ of the next step. In this paper, we use the $BAB$ integrator form in our tests, but we note that the efficiency can be improved by using $ABA$. The goal of this paper is not to optimize efficiency. The equations of motion for Hamiltonian $A$ are (Rein et al., 2019),

[TABLE]

where $\boldsymbol{V}_{i}=\boldsymbol{P}_{i}/m_{i}$ is the barycentric velocity. Using velocities instead of momenta is one way to ensure we can describe test particles. We used (Rein et al., 2019),

[TABLE]

The equations for Hamiltonian $B$ are,

[TABLE]

Since $\ddot{\boldsymbol{Q}}_{i}$ is [math], $\dot{\boldsymbol{Q}}_{i}$ is a constant. $\dot{\boldsymbol{V}}_{i}$ is constant because $\dot{\boldsymbol{Q}}_{ij}$ is [math].

3 Modified differential equations

Although we have written the equations of motion for $A$ and $B$, the integrators of (3) can be described by sets of ODEs called modified differential equations (MDEs) (Hairer et al., 2006; Hernandez & Bertschinger, 2018). We derive these for the first-order integrator,

[TABLE]

The Section 3.1 example also uses a first order integrator. Here,

[TABLE]

The set of ODEs comes from $\tilde{H}$. For $i>0$,

[TABLE]

We defined,

[TABLE]

where $\boldsymbol{P}_{\sun}=\sum_{i>0}\boldsymbol{P}_{i}$ (Hernandez & Dehnen, 2017). Eqs. (15) are not Taylor series in $h$: $h$ is a constant parameter. The MDEs do not need to converge in $h$, they correctly describe the dynamics. Eqs. (16) are generally nonzero for massive and massless particles. $\boldsymbol{a}_{0i}$ is the heliocentric velocity and $\boldsymbol{b}_{0i}$ is the Newtonian force: in the limit $h\to 0$, these are the only surviving terms and the MDEs exactly solve the gravitational $N$-body equations. We can see that the integrator (13) suffers from the same problem as integrators (3), and any integrator split into the $A$ and $B$ from (1): the $n=1$ problem is not solved exactly, even when $K(Q_{ij})$ and its derivatives are 0. We see this because terms $\boldsymbol{a}_{1i}$ and $\boldsymbol{b}_{1i}$ are nonzero, although $\boldsymbol{a}_{1i}=0$ for test particles. The coefficients as the power of $h$ increases depend on higher derivatives of $K$. Up to first order in $h$, two derivatives of $K$ are needed. As we integrate these MDEs over discontinuities in those derivatives, unresolved timescales can occur (Hernandez, 2019), leading to large errors. As $h$ is increased, these discontinuities will become more important. The MDEs corresponding to integrator $BAB$, Eq. (3a), also the one considered in Chambers (1999) and Rein et al. (2019), do not have form (15); the leading error term is $\mathcal{O}(h^{2})$. However, Eqs. (15) and (16) let us see how discontinuities would enter any MDEs.

3.1 Example error analysis

Here we show and explain an example of a finite DCO $K$ giving rise to secular errors in orbits. This subsection can be omitted without loss of continuity. Consider the simple harmonic oscillator (SHO), with Hamiltonian,

[TABLE]

We are concerned in this work with chaotic orbits, but for this example, we consider a periodic orbit. We construct a hybrid integrator to solve this problem. This integrator was considered in (Hernandez, 2019, Section 4.1). The Hamiltonian is split into,

[TABLE]

The integrator is,

[TABLE]

$H_{1}$ is solved using Bulirsch–Stoer (as in Hernandez (2019)). Let $x=(q+2)/2$. If $x<0$, $K(q)=0$. If $x>0$, $K(q)=1$. Otherwise, $K(q)=x$. This $K$ is a $C^{0}$ function, so integrator (19) is not symplectic, according to Hernandez (2019). As in Hernandez (2019), initially, $q=1$ and $p=0$. $h=P/100$, where the exact period is $P=2\pi$. The runtime is $t=2P$. The energy error as a function of time is plotted in Fig. 1.

In the first period, four steps, shown with stars, integrate over a discontinuity in $dK/dr$. Stars locate the end of those steps. A green star means the region $0<K<1$ was left, while a blue star means $0<K<1$ was entered. The bottom subpanel zooms in on a portion of the top subpanel, and shows the energy error floor shifting. We now explain the magnitude and sign of this shift. A succession of such shifts lead to long-term error deterioration.

First, we calculate $\tilde{H}$, given by the BCH formula. If the series converges, we can ignore high powers in $h$ and assume $\tilde{H}$ is conserved. We have,

[TABLE]

The Poisson brackets are calculated as,

[TABLE]

Note that $\tilde{H}$ describes exact, continuous dynamics, explained by MDEs, although the integrator samples discrete points along the solution path. We calculate the error across a discontinuity. Using (20) and (21), and ignoring higher $h$ powers,

[TABLE]

$H_{0}$ is the initial energy $(1/2)$. \bigg{\rvert}_{+} indicates the limit from the right (forward time) of the discontinuous point while \bigg{\rvert}_{-} is the limit from the left. $a=h/(8H_{0})$. $s$ is $+1$ if $0<K<1$ was entered while it is $-1$ if $0<K<1$ was left. We calculate the first four $\delta$ values of Fig. 1 approximately, because only data at discrete values is available from the integrator. We use the coordinates after integrating over the discontinuity. Labelling each $\delta$, $\delta_{1}=-6.2\times 10^{-5}$, $\delta_{2}=0.027$, $\delta_{3}=0.039$, and $\delta_{4}=-5.8\times 10^{-5}$. The maximum $\delta$ and maximum $pq^{2}$ are, respectively, $0.039$ and $\approx 2.5$. The jump in the bottom subpanel of Fig. 1 is $+0.035$, so this error analysis provides an approximation for this jump. Reducing the step by 10 to $h=P/1000$, there are again four discontinuities over the first period. The maximum magnitude of $pq^{2}$ is similar, $\approx 2.3$, and $a$ is $10$ times smaller. The error floor jumps in this case by $+0.004$, so the error analysis result again matches approximately the data.

4 Numerical experiments

4.1 Chaotic orbits

We return attention to the $N$-body problem. We have yet to specify the functional form of $K$. To discuss $K$, first consider a function $x=(r-h_{1}R)/(h_{2}R)$, a function used by Wisdom (2017). $r$ is the separation between two non-stellar bodies, and $R$ is a constant. If $x<0$, $K=0$. If $x>1$, $K=1$. Otherwise, we use the following progressively smoother functions,

$C^{0}$: $K=x$ 2. 2.

$C^{1}$: $K=-2x^{3}+3x^{2}$ 3. 3.

$C^{2}$: $K=6x^{5}-15x^{4}+10x^{3}$ 4. 4.

$C^{3}$: $K=-20x^{7}+70x^{6}-84x^{5}+35x^{4}$. 5. 5.

$C^{4}$: $K=70x^{9}-315x^{8}+540x^{7}-420x^{6}+126x^{5}$ 6. 6.

$C^{5}$: $K=-252x^{11}+1386x^{10}-3080x^{9}+3465x^{8}-1980x^{7}+462x^{6}$

The functions increase monotonically and have odd-symmetry about $x=1/2$. They are displayed in Fig. 2.

Wisdom (2017) used a $C^{2}$ function and Chambers (1999) used a $C^{1}$ function. The code MERCURY, however, described in Chambers (1999), uses a $C^{2}$ function. Thus, the forces in MERCURY are smoother than those of Wisdom (2017), although the former has symplecticity and time-reversibility errors (Hernandez, 2016).

Following Wisdom (2017), if, for a given pair of non-stellar objects, $x>h_{3}R$ at the start of $A$, we use Kepler solvers to solve the Hamiltonian terms with their indices $i$ and $j$. $h_{3}R>1$ by construction. Otherwise, Bulirsch–Stoer is used for that pair in $A$, to account situations in which $K$ becomes less than $1$ during the evolution of the $A$ Hamiltonian.

We consider a simple chaotic system with two degrees of freedom, the restricted three-body problem. It has one constant of motion. In an inertial frame it is,

[TABLE]

where $\omega$ is the binary angular velocity, $\boldsymbol{\omega}=\sqrt{G(m_{1}+m_{2})/a^{3}}\hat{\boldsymbol{z}}$, and

[TABLE]

$a$ is the semi-major axis of the circular orbit, $m_{1}$ and $m_{2}$ are the masses, and $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$ are their positions. $m_{1}$ is the star. $\boldsymbol{x}$ and $\boldsymbol{v}$ are the position and velocity of the test particle. In units of au, day, and solar masses, the test particle has heliocentric position and velocity

[TABLE]

$m_{1}=1$ and $m_{2}$ is defined through $\mu=m_{2}/(m_{1}+m_{2})=0.01$. $a=5.2$. Initially, all bodies are aligned with the $x$-axis, with the star on the left. The center of mass is at [math]. This problem was analyzed by Wisdom (2017), so a comparison is possible in integrator performance. We let $R$ be the Hill radius, $R=a(m_{2}/(3m_{1}))^{1/3}$, $h_{1}=h_{2}=1.5$, and $h_{3}=4.0$, as in Wisdom (2017). For all tests in this paper, we checked that when Kepler solvers were used to solve $A$, a close encounter condition, $K<1$, was not detected after evolving the $A$ Hamiltonian. This justifies the choice of $h_{3}$.

Using function 3, $h=8$ days, and running for a total time $t=500$ yrs, we show the trajectory in the rotating frame in Fig. 3.

We have also indicated with green stars if a transition out of the region $0<K<1$ happened during the timestep yielding the coordinates. This means an integration over a discontinuity was performed. We also show with blue stars when a transition into $0<K<1$ happened. $54$ steps integrated over a discontinuity, out of a total of $22,830$ steps. Transitions should occur at separations of $r=3.0R$ and $r=1.5R$, respectively, which we verified numerically.

Next, we measure the Jacobi constant error for long-term integrations over $10,000$ yrs using various smoothing functions. We expect that smoother functions will have better error. We use $h=8$ days, but to reduce variations in the error, we plot the median error every 1000 steps, so that a total of 456 points are plotted. The result is in Fig. 4. We also plot data from Figure 9 of Wisdom (2017) (“Wisdom”). For these data, median absolute errors every $20$ steps are plotted, so that there are also a total of $456$ points in the range $t=0$ to $t=10,000$ yrs. We also plot corrected data from Fig. 10 of Wisdom (2017) as “Wisdomc.” These data behave similarly to the median errors from above, and the two curves lie on top of each other. We indicate by integers $n$ the DCO.

Integrations with smoother transition functions have smaller error, although there is no significant improvement from $n=4$ to $n=5$. We can improve on the Wisdom (2017) data simply by increasing smoothness at virtually no additional cost. Running to 500 yrs, the $C^{5}$ experiment took about $1\%$ longer in wall-clock time than the $C^{2}$ test. The $C^{5}$ test was faster than the $C^{0}$ test, but note different runs may have a different number of evaluations of the polynomial function. We used the MATLAB software for this timing test. But our goal is not to improve on the Wisdom (2017) result— if our focus were on accuracy, we could reduce the error on all our curves by using instead the $ABA$ forms. Our goal is to demonstrate the importance of the transition function.

Because the error as a function of time strongly depends on the initial conditions, we want to see if our results hold for a variety of orbits. Thus, we perturb the initial conditions (25), but not so much that the initial conditions’ Jacobi constant no longer satisfies the condition for a chaotic exchange oribt (Dehnen & Hernandez, 2017, Section 3.3). For the test of Fig. 3, the variation in either test particle position coordinate is $\approx 10$, while the variation in either velocity is $\approx 0.03$. Thus, we set a position and velocity scale as $x_{\text{s}}=10$ and $v_{\text{s}}=0.01$, respectively. Then $8$ initial conditions are generated by perturbing one of the positions or velocities at a time from the initial conditions above. The positions are perturbed by $\pm x_{\text{s}}\delta$ and velocities by $\pm v_{\text{s}}\delta$, where $\delta=10^{-3}$. For the long timescales we study in our tests, we can perturb initial conditions even by small amounts like $10^{-15}$. Assuming a Lyapunov time of $4$ yrs, if we perturb the initial conditions by $10^{-15}$, the trajectories will differ by unity after $\approx 136$ yrs. However, in one such test, significant differences between orbits were not apparent until after $400$ yrs. We repeat the test from Fig 4, but with $8$ initial conditions with each transition function. We plot their errors by again taking median values and using thin lines. A thick line represents an average in log space, or the geometric mean, of the eight curves. The result is shown in Fig. 5.

This plot confirms the improvement in error from increasing smoothness. As in Fig. 4, there is improvement in error until $n=4$. It is unclear whether this is an improvement from $n=4$ to $n=5$.

We verified the orbits from nearby initial conditions diverge quickly. For a certain integrator and the unperturbed initial conditions (25), Dehnen & Hernandez (2017) reported a Lyapunov time of $4$ yrs. For the $n=2$ case, we compared the orbit in which the $x$ coordinate is perturbed by $-x_{\text{s}}\delta$ and the orbit for which $y$ is perturbed by $+x_{\text{s}}\delta$. We calculated the Euclidean norm of the distance between the two orbits, which reached unity before $5$ yrs. The Euclidean norm in velocity space reached its saturation of approximately $10^{-2}$ at a similar time.

Based on the MDEs from (15), we expect that as $h$ is decreased, the smoothness of the transition will matter less, as its impact on the error becomes smaller. To show this, we rerun the initial conditions (25) and their $8$ perturbed orbits using a step $10$ times smaller, $h=0.8$ days. Median values are taken every $10,000$ steps for plotting, so that the same number of points are plotted as in Fig. 5. We compare the Jacobi constant errors in Fig. 6 between the two timestep runs, up to $t=1,000$ yrs. The top panel uses $h=0.8$ days. There is no clear difference once $n>1$ in this case, while for $h=8$ days case, differences are hard to see for $n>2$. Also note all thick curves have smaller error for the $h=0.8$ days case.

We want to demonstrate that large errors in the Jacobi constant are associated with integrations over discontinuities. In Fig. 7 we plot the error in the Jacobi constant for the $C^{2}$ function over the first $50$ yrs, using $h=8$ days. We also plot when an integration over a discontinuity occurred. We find the first discontinuity is associated with a secular increase in the Jacobi constant error. Further discontinuities are associated with other error spikes.

4.2 Revisiting regular orbits

This work has focused on the effects of smoothness on chaotic orbits. Hernandez (2019) showed that regular orbits were generally ensured with $C^{1,1}$ integrators. We explore here whether increasing smoothness past $C^{1,1}$ improves the regular orbits.

We consider the Kepler orbit in Hernandez (2019) and use map (9) in that paper. We run it for $t=10000P$, $h=P/100$, and eccentricity $e=0.7$. $P=2\pi$ is the period. The initial conditions are at apoapse ($r=1+e$, $p=0$, using the same notation as Hernandez (2019)). An error is calculated at each apocenter passage using a cubic polyynomial interpolation. Specifically, assume apocenter passage occurs between $t_{1}$ and $t_{1}+h$. Then, apocenter occurs at,

[TABLE]

where $a$, $b$, $c$, and $d$ are known constants and $t$ is found from $p=\dot{r}=0$. Fig. 8 plots the errors vs time for this experiment, where the smoothness of the map is varied. The amplitudes for the $C^{1}$, $C^{2}$, and $C^{\infty}$ runs are, approximately, $1.9\times 10^{-3}$, $9\times 10^{-4}$, and $7\times 10^{-4}$. Thus, our chaotic results apply to regular orbits; error amplitudes decrease as smoothness is increased. But a minimum smoothness is needed to guarantee the orbits are regular.

We verified that using a given function as a force, rather than Hamiltonian, switching function, gives smaller errors. This is expected as the force is obtained from differentiating the potential. We also used other conventional integrators besides Bulirsch–Stoer which are less sensitive to smoothness, but found no significant impact to our results.

5 Conclusion

Symplectic integration has been popular and successful at describing $N$-body dynamics in astrophysics. Symplectic integration is part of a broader category of integration called geometric integration (Hairer et al., 2006) in which solutions to ordinary differential equations (ODEs) are calculated through approximations which respect geometry of the ODEs. This is useful for any long-term dynamical studies. Symplectic integration is supposed to conserve the Poincaré invariants associated with Hamiltonian dynamics, leading to bounded errors over long time intervals. Popular symplectic integrators (Wisdom & Holman, 1991; Springel, 2005) also aim to preserve the other geometric features of the ODEs such as the time-symmetry, and the angular and linear momentum conservation.

Another geometric feature which is perhaps overlooked is the smoothness of an ODE. The $N$-body Hamiltonian is infinitely smooth, or $C^{\infty}$ for physical separations greater than [math]. Symplectic integrators are described by modified differential equations, which are also derived from a Hamiltonian with some smoothness. A number of popular symplectic methods, and perhaps the most versatile and useful ones (Chambers, 1999; Duncan et al., 1998; Springel, 2005) are only $C^{n}$ smooth, where $n$ is a small integer, or worse, one might argue a smoothness concept is not applicable to them. Any method with adaptive block stepping also has reduced smoothness. This work has explored the role of smoothness in the accuracy of $N$-body integrations. To explore this question, we have considered hybrid (Chambers, 1999; Rein et al., 2019; Wisdom, 2017; Hernandez, 2019) integrators, whose smoothness is easily tuned through a switching function. Hybrid integrators are $C^{\infty}$ in most of phase space, except for a few points in phase space which are rarely or never hit, in which there is a smaller differentiability class order (DCO) (except if an expensive $C^{\infty}$ transition function is used). However, Hernandez (2019) has shown that the performance of hybrid integrators degrades for regular orbits even when those points are never encountered in practice, and explained why.

This work has focused on the broad class of chaotic orbits, and how their errors might be affected by smoothness of the method. We have provided evidence that Hamiltonian smoothness has significant impacts on $N$-body performance. In a test first considered by Wisdom (2017) of a restricted three-body problem, the integrator performance significantly improved as $n$ was increased from [math] to $4$ by about $5$ orders of magnitude in the Jacobi constant error. We argued that this result can be explained through the MDEs. We also verified the expectation that smoothness should become less important as the timestep goes to [math]. We also showed regular orbits were improved by smoothness. Our results indicate that just as using symplectic integrators has become a pillar of dynamical astronomy, using integrators that are as smooth as possible is a consideration which is arguably just as important. Increasing smoothness, or the DCO, may not be easy to do in the case of block time-stepping methods, but for other methods, like hybrid integrators, increasing smoothness may be negligibly more computationally expensive, as we explored.

We remark on the relationship between smoothness and symplecticity: Hernandez (2019) has argued that a symplectic integrator should have at least $C^{1,1}$ smoothness ($C^{1}$ and Lipschitz continuous in the ODEs), because these appear to guarantee the existence of periodic orbits. This was verified numerically, and was predicted by Arnaud & Xue (2016), with some caveats described in Hernandez (2019). According to this criteria, we have explored both symplectic and non-symplectic integrators in this work, depending on the smoothness of the transition function.

Finally, we note that we have not investigated in detail the role smoothness plays in nonsymplectic integrators, but it is reasonable to assume smoothness can improve any integrator described by MDEs, symplectic or not.

6 Acknowledgements

I thank Jun Makino, Scott Tremaine, Matt Payne, Jack Wisdom, and the anonymous referees for comments that improved this work.