Hybrid Symplectic Integrators for Planetary Dynamics

TL;DR

This paper reevaluates the switching function criteria in hybrid symplectic integrators like MERCURY, proposing simpler alternatives that maintain accuracy in planetary dynamics simulations.

Contribution

It clarifies the criteria for choosing switching functions and introduces a simpler scheme using a Heaviside function without sacrificing accuracy.

Findings

Heaviside switching function performs as well as polynomial functions in short-term accuracy.

Original polynomial switching function choice in MERCURY was based on incorrect motivation.

Simpler switching schemes can be used without loss of performance.

Abstract

Hybrid symplectic integrators such as MERCURY are widely used to simulate complex dynamical phenomena in planetary dynamics that could otherwise not be investigated. A hybrid integrator achieves high accuracy during close encounters by using a high order integration scheme for the duration of the encounter while otherwise using a standard 2nd order Wisdom-Holman scheme, thereby optimizing both speed and accuracy. In this paper we reassess the criteria for choosing the switching function that determines which parts of the Hamiltonian are integrated with the high order integrator. We show that the original motivation for choosing a polynomial switching function in MERCURY is not correct. We explain the nevertheless excellent performance of the MERCURY integrator and then explore a wide range of different switching functions including an infinitely differentiable function and a Heaviside function. We find that using a Heaviside function leads to a significantly simpler scheme compared to MERCURY, while maintaining the same accuracy in short term simulations.

Equations

$$H$$

$$H_0$$

$$H_K$$

$$H_I$$

$$H_J$$

$${\hat{H}}_{WH}(dt)$$

$$\frac{1}{12}[[{\hat{H}}_J,{\hat{H}}_K],{\hat{H}}_K],\frac{1}{24}[[{\hat{H}}_J,{\hat{H}}_K],{\hat{H}}_J],$$

$$\frac{1}{12}[[{\hat{H}}_I,{\hat{H}}_K],{\hat{H}}_K],\frac{1}{24}[[{\hat{H}}_I,{\hat{H}}_K],{\hat{H}}_I],$$

$$\text{and}\frac{1}{12}[[{\hat{H}}_J,{\hat{H}}_K],{\hat{H}}_I].$$

$$\underset{\widehat{=}[{\hat{H}}_J,{\hat{H}}_K{]}^{-1}(-dt)}{\underbrace{{\hat{H}}_K(-dt)\circ {\hat{H}}_J(-dt)\circ {\hat{H}}_K(dt)\circ {\hat{H}}_J(dt)}}\circ {\hat{H}}_K(-dt)$$

$$\circ \underset{\widehat{=}[{\hat{H}}_J,{\hat{H}}_K](dt)}{\underbrace{{\hat{H}}_J(-dt)\circ {\hat{H}}_K(-dt)\circ {\hat{H}}_J(dt)\circ {\hat{H}}_K(dt)}}\circ {\hat{H}}_K(dt).$$

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

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

$${\mathbf{V}}_i\equiv \frac{{\mathbf{P}}_i}{m_i}\text{and therefore}{\dot{\mathbf{V}}}_i=\frac{{\dot{\mathbf{P}}}_i}{m_i}.$$

$$\frac{\partial K(Q_{ij})}{\partial 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$$

$$L(r)$$

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

$$L_{\mathrm{merc}}(r)=\{\begin{array}{ll}0& y\le 0\\ 10y^3-15y^4+6y^5& 0<y<1\\ 1& y\ge 1,\end{array}$$

$$y=\frac{r-0.1r_{\mathrm{crit}}}{0.9r_{\mathrm{crit}}}$$

$$L_{\inf }(r)$$

$$f(y)=\{\begin{array}{ll}0& y\le 0\\ e^{-1/x}& y>0\end{array}.$$

$${\mathbf{Q}}_i=\{\begin{array}{ll}{\mathbf{q}}_i-{\mathbf{q}}_0& i\ne 0\\ \frac{1}{M}{\sum }_{j=0}^{N-1}{m}_j{\mathbf{q}}_j& i=0.\end{array}$$

$${\mathbf{P}}_i=\{\begin{array}{ll}{\mathbf{p}}_i-\frac{m_i}{M}{\sum }_{j=0}^{N-1}{\mathbf{p}}_j& i\ne 0\\ {\sum }_{j=0}^{N-1}{\mathbf{p}}_j& i=0.\end{array}$$

Full text

1 Introduction

Since ancient times astronomers have predicted the locations of planets in the night sky. But only the advent of computers has made it possible to calculate the orbital evolution of planetary systems over millions and even billions of years (Laskar & Gastineau, 2009). One major breakthrough was the development of mixed variable symplectic maps for the $N$-body problem which we will refer to as Wisdom-Holman integrators (Wisdom & Holman, 1992). The idea behind the Wisdom-Holman integrator is to split the motion of a planet into two steps: a dominant Keplerian motion around the star (the Kepler step), and the motion due to small interactions between the planets (the interaction step).

The Wisdom-Holman scheme works well as long as the interactions between planets can be considered perturbations to their predominant Keplerian motion around the star. However, if two planets come close to each other, the planet-planet interactions can become the dominant part of the motion. In that case the Wisdom-Holman integrator becomes inaccurate and might require excessively small timesteps to achieve an acceptable solution.

Duncan et al. (1998) and Chambers (1999) provide a solution to this problem in the form of hybrid symplectic integrators. The idea is to move terms from the interaction part to the Keplerian part during close encounters, to ensure that one part always remains a perturbation. Although solving the Keplerian part during a close encounter now becomes more difficult than simply solving Kepler’s equation, it can be done relatively efficiently with a high order integrator such as a Burlish-Stoer or Gauß-Radau scheme. As long as close encounters happen infrequently, the high order integrator is rarely used and has a negligible effect on the runtime. The scheme is therefore almost as fast as a Wisdom-Holman integrator. Furthermore, the scheme is also symplectic, ensuring good long term conservation of energy and angular momentum.

In this paper we look specifically at the MERCURY integrator. MERCURY is freely available online and has become a widely used tool in running many otherwise impossible types of simulations. We discuss inconsistencies in the original paper by Chambers (1999) which describes the MERCURY algorithm111A similar discussion in Duncan et al. (1998) is correct, but it describes a slightly different and more complex scheme.. We note that Wisdom (2017) also looked at these inconsistencies and came independently of us to the same conclusion. Despite these inconsistencies, the resulting scheme has very good characteristics. We provide mathematical and physical explanations for them in this paper. With this new understanding, it is easy to construct new hybrid integrators. We describe one new integrator that is much simpler than MERCURY, but achieves the same accuracy.

2 Hybrid Symplectic Integrators

We begin by introducing the notation and coordinate system we use. We have one star and an additional $N$ planetary bodies with Cartesian coordinates $\mathbf{q}_{i}$ and canonical momenta $\mathbf{p}_{i}$ in some inertial frame. We will identify the stellar object with $i=0$ and the planets with $i>0$. The $N$-body Hamiltonian in these coordinates is

[TABLE]

To evolve this system forward in time we could calculate the equations of motion for this Hamiltonian. However, in general the resulting differential equations cannot be solved analytically and are usually stiff and thus hard to solve numerically.

The hybrid integration schemes discussed in this paper make use of democratic heliocentric coordinates $Q_{i}$ and corresponding canonical momenta $P_{i}$ (see Appendix A). The advantage of these coordinates is that the Hamiltonian can be written as a sum of four terms:

[TABLE]

where

[TABLE]

In the above splitting, we’ve introduced the arbitrary scalar function $K(r)$ which depends only on the distance between pairs of particles. Note that the terms involving $K$ cancel out when adding $H_{K}$ and $H_{I}$ and thus do not change the evolution of the system.

In this paper, we will use the letter $K$ to label switching functions that appear in the Hamiltonian and refer to them as Hamiltonian switching functions. Later, we will also encounter force switching functions that first appear in the equations of motions. To avoid any confusion, we will use the letter $L$ to label force switching functions. There is a relation between $K$ and $L$ which we derive in Sec. 3.

Each of the above terms has a physical interpretation. $H_{0}$ describes the motion of the centre of mass. The term $H_{K}$ describes the Keplerian motion of the planets, ignoring planet-planet interactions (while $K=1$). This term is separable and can be solved for each planet independently. The term $H_{I}$ corresponds to the planet-planet interactions. The term $H_{J}$ is related to the barycentric motion of the star and is referred to as the jump term (it changes the positions of particles instantaneously but keeps the momenta constant).

We can calculate the equations of motion for each of these Hamiltonians. The resulting differential equations for $H_{0}$, $H_{I}$, and $H_{J}$ are trivial to solve. Only the equations for $H_{K}$ require more work, i.e. solving Kepler’s equation (while $K=1$). Evolving the system under the influence of a Hamiltonian can be thought of as an operator acting on a state corresponding to the initial conditions. The idea of an operator splitting integration method is to approximate the evolution of the full system by consecutively evolving the system under the partial Hamiltonians. We can construct arbitrarily high order integrators by applying these operators in the right order for different amounts of time (Yoshida, 1990). In particular, a second-order integrator can be constructed from the above splitting by using the following chain of operators

[TABLE]

where an operator $\widehat{H}_{A}(dt)$ corresponds to evolving the system under the Hamiltonian $H_{A}$ for a time $dt$. Note that we recover the standard Wisdom-Holman integrator in democratic heliocentric coordinates if we set $K=1$.

For any finite timestep $dt$, $\widehat{H}_{WH}$ is only approximately equal to $\widehat{H}$. We can analyze the error of this splitting scheme with the help of the Baker-Campbell-Hausdorff (BCH) formula. The splitting scheme error arises because some pairs of operators do not commute. The error consists of an infinite series of commutators between operators. In our case, the action of an operator corresponds to solving the differential equations resulting from its respective Hamiltonian. As a result, the splitting scheme error can be expressed as a series of nested Poisson brackets (see e.g. Hernandez & Dehnen, 2017). One caveat to keep in mind when working with the BCH formula is that it is only a formal asymptotic series and does not necessarily converge everywhere in phase space. We come back to this issue later.

The commutators that appear in the BCH formula can be grouped together by powers of the timestep $dt$. The specific splitting scheme in Eq. (7) is time symmetric and therefore only even orders of $dt$ appear (Yoshida, 1990). Further note that $\widehat{H}_{J}$ commutes with $\widehat{H}_{I}$. Thus, the lowest order non-zero error terms that appear are:

[TABLE]

In this paper, we do not attempt to calculate approximations to these expressions analytically using Poisson brackets, but do so numerically. This allows us to easily extend the calculations to arbitrarily high order (7th order in our case).

There are two way to approach this. Consider a commutator of the form $[\widehat{A},\widehat{B}]$. Starting from a state $z$, we can calculate $\delta z\equiv(\widehat{A}\widehat{B})\,z-(\widehat{B}\widehat{A})\,z$, and similarly for higher order commutators. If $\delta z=0$, then the commutators commute. To quantify the statement, we can use any norm on $\delta z$. The measure we are typically interested in is the energy, or more specifically the energy error. Note that calculating the energy of $\delta z$ does not make sense. However, we can calculate the energy error by comparing the energy of the state $z$ to that of the state $z+\delta z$. Alternatively, we can calculate the state $z^{\prime}\equiv(\widehat{B}^{-1}\widehat{A}^{-1}\widehat{B}\widehat{A})\,z$. Here, $z^{\prime}$ is a valid state and we can directly compare the energy of the states $z$ and $z^{\prime}$.

In this paper we use the latter approach. Note that it is a simple chain of operators and no additions or subtractions of state vectors are needed. This approach avoids some issues with finite floating point precision compared to the standard approach where one needs to subtract nearly equal numbers multiple times when calculating the energy difference of the states $z$ and $z+\delta z$. As a specific example, consider the first (third order) commutator that appears in the BCH formula of the hybrid symplectic intergators: $\small[\small[\widehat{H}_{J},\widehat{H}_{K}\small],\widehat{H}_{K}\small]$. The effect of this commutator can be calculated by the following chain of operators:

[TABLE]

The two approaches are not exactly equal, only to within order $dt$. Our approach is only an approximation of the commutator appearing in the BCH formula and is accurate to within $\mathcal{O}(dt)$ only. We could apply the Zassenhaus formula (the dual of the BCH formula) to calculate higher order corrections, which are again chains of the same operators. However, the approximations are good enough for our purpose: they will allow us to demonstrate which are the dominant contributions to the error terms in the BCH formula and understand their origins, especially for higher order commutators which are otherwise not easily tractable.

3 The MERCURY Code

The Wisdom-Holman map (setting $K=1$ in the above notation) performs poorly during close encounters because some of the commutators in the error term become large. To construct an integrator that can more accurately resolve close encounters, Chambers (1999) suggests to use a non-constant Hamiltonian splitting function $K$ as a basis for the MERCURY scheme.

The idea, as outlined in Chambers (1999), is that the switching function keeps the interaction Hamiltonian small even during close encounters between planets. As a result the switching scheme errors from the BCH formula should remain small as well.

Let us use Hamilton’s equations and derive the equations of motion for the Hamiltonian $H_{K}$:

[TABLE]

where we have introduced the democratic heliocentric velocities

[TABLE]

We also used the fact that the gradient of $K$ satisfies

[TABLE]

The equations of motion above are not those of Keplerian motion but now also include interaction terms. Nevertheless they can be solved with a high order integrator such as a Bulirsch-Stoer algorithm or a high order Gauß-Radau integrator (Everhart, 1985; Rein & Spiegel, 2015).

Similarly, the equations of motion for $H_{I}$ can be derived as

[TABLE]

These equations remain trivial to solve even in the case of $K\neq 1$.

Somewhat surprisingly, the equations of motion derived above are not the equations of motion implemented in the MERCURY code. Specifically, the terms involving $K^{\prime}(Q_{ij})$ are not present in the publicly available version of MERCURY. Therefore, the MERCURY code does also not correspond to the Hamiltonian splitting in Eq. (2).

This then begs the question222 This inconsistency was discussed at the Toronto Meeting on Numerical Integration Methods in Planetary Science in 2017 and an explanation was first published by Wisdom (2017). We independently came to the same conclusion. : What is MERCURY solving? To find the answer to this question, we define a new switching function

[TABLE]

With this definition, we can rewrite the equations of motion in terms of $L$. For example, Eq. (10) becomes

[TABLE]

Note that no derivatives of $L$ appear in the equations of motion. We can think of $L$, which is changing between [math] and $1$, as a weight in the inter-particle forces. We therefor refer to it as the force switching function. If we set $L$ to the piecewise polynomial function

[TABLE]

with

[TABLE]

then we recover exactly the integration scheme as implemented in the MERCURY code333Note that this is not the function described in the paper which has roughly the same shape but is slightly different from the one implemented in the publicly available version of MERCURY. We here use the function implemented in the MERCURY code.. In the above equation, $r_{\rm crit}$ is some critical switching distance. The value depends on the specific application but it typically a small multiple of the mutual Hill radius of the particles.

To find out what the corresponding Hamiltonians are that this scheme is integrating, we need to solve the differential equation in Eq. (15) to find the Hamiltonian switching function $K$, given the force switching function $L$. The solution is a one parameter family due to a gauge symmetry. We choose the solution that is continuous by matching $K_{\rm merc}(r_{\rm crit})=1$. We plot the functions $K_{\rm merc}$ and $L_{\rm merc}$ in the top panel of Fig. 1. Note that the function $K_{\rm merc}$ does not reach exactly [math] at $r=0.1r_{\rm crit}$. This means that there is a finite contribution to the interaction Hamiltonian $H_{\rm I}$ for any choice of $r>0$ (see function $K_{\rm merc}(r)/r$ in the bottom panel of Fig. 1).

This is inconsistent with the discussion of Chambers (1999) which argues that the switching function keeps the interaction Hamiltonian small. That argument would only be correct if the force switching function $L_{\rm merc}$ were used as a Hamiltonian switching function (in the MERCURY code $L_{\rm merc}$ is used as a force switching function). As an illustration, we also plot the function $L_{\rm merc}(r)/r$ in the bottom panel of Fig. 1. This is a potential which now reaches [math] at $r=0.1r_{\rm crit}$. However, whereas using $L_{\rm merc}$ as the Hamiltonian switching function would lead to a valid symplectic switching integrator, it will not perform as well as an integrator using $K_{\rm merc}$. Looking at the bottom panel of Fig. 1, this is because $L_{\rm merc}(r)/r$ has a local minima around $0.8r_{\rm crit}$ leading to an unphysical repulsive force between particles in some intermediate regime. One consequence of the force changing sign is that the differential equations corresponding to $L_{\rm merc}(r)/r$ are stiffer during the close encounter. Note that in both cases, particles will feel no force from the potentials in the interaction Hamiltonian while $r<0.1r_{\rm crit}$.

It is worth pointing out that for small $r$ the shape of the effective potential $K_{\rm merc}(r)/r$ is similar to that of potentials used in other $N$-body systems where one simply ignores the close range interactions by having a finite smoothing length (similar to a kernel in smooth particle hydrodynamics).

We would like to stress that despite this inconsistency in the derivation of the equations of motion, the MERCURY algorithm is solving a Hamiltonian system that converges to the original system in the limit of $dt\rightarrow 0$.

For the remainder of this paper, we will refer to the force switching function $L(r)$ simply as the switching function. The corresponding Hamiltonian switching function $K(r)$ can always be calculated by solving Eq. (15).

4 Switching functions

We now explore different force switching functions and their effects on the accuracy of the integration. Our motivation comes from the fact that the function used by MERCURY has been determined heuristically and, as we have shown above, the original justification is not correct. It is thus not clear if one can improve the integration scheme by simply choosing a better switching function. In this paper, we will present and discuss results for the following four switching functions. These are extreme cases which will allow us to explain the general ideas behind choosing switching functions, which can then be applied to any arbitrary switching function.

(i) We set $L(r)=1$. The integrator simply becomes the (non-switching) Wisdom-Holman integrator in democratic heliocentric coordinates.

(ii) We set $L(r)$ to the discontinuous Heaviside function with a jump at $r=r_{\rm crit}$.

(iii) Same as (ii) but the value of $L$ is kept constant for the duration of the entire timestep. To do that we approximate the trajectories along straight lines for the duration of the timestep and set $L=0$ if the trajectories come closer than $r_{\rm crit}$, and $L=1$ otherwise. There are cases where this breaks time reversibility. We could avoid this by doing this procedure iteratively. However, in all tests that we have performed this does not seem to be an issue. Further, note that this procedure formally makes the integrator phase space dependent in a similar way than adaptive timesteps do. We comment on whether this has any real world consequences in Sec. 7.

(iv) We use the switching function $L_{\rm merc}(r)$ given in Eq.(17). Thus the integrator becomes the MERCURY integrator.

(v) Finally, we consider an infinitely differentiable switching function given by

[TABLE]

where $y$ is given by Eq. (18) and $f(y)$ is

[TABLE]

The shape of the function $L_{\rm inf}$ is very similar to that of $L_{\rm merc}$ (shown in Fig. 1). However the function $L_{\rm inf}(r)$ and all its derivatives are smooth and continuous for all $r$, even at the boundaries. Note however, that just as $L_{\rm merc}(r)$, $L_{\rm inf}(r)$ is not analytic, i.e. a Taylor series expansion will not converge globally.

5 Test Setup and Mercurius

We now look at one representative test case in which a close encounter occurs between two giant planets. We work in units where $G=1$ and the stellar mass is $m_{0}=1$. Both planets have a mass similar to that of Jupiter, $m_{1}=m_{2}=10^{-3}$. The inner planet is on an orbit with semi-major axis $a_{1}=1$ and eccentricity $e_{1}=0.1$. The outer planet is on a nearby orbit with $a_{2}=1.1$ and $e=0.2$. We use a timestep corresponding to 0.5% of the inner planet’s period. We adjust the phases such that a close encounter occurs at $t\approx 0$ with a close approach distance of approximately $5\cdot 10^{-5}$.

To run these simulations we use MERCURIUS, our own implementation of the MERCURY integrator444In reference to and appreciation of MERCURY, we call our implementation MERCURIUS (the Latin word for Mercury).. We here only give a short overview of the implementation as we plan to publish a detailed code description paper. The MERCURIUS integrator is freely available as part of the REBOUND integrator package. Internally it uses the Kepler solver of the WHFast integrator (Rein & Tamayo, 2015) and the IAS15 integrator (Rein & Spiegel, 2015) to evolve $H_{\rm K}$. Using MERCURIUS allows us to easily experiment with different switching functions. We can also call the individual operators manually which allows us to calculate the commutators numerically as described above. We also ran tests with the original MERCURY code and implemented wrapper functions to make the fortran code callable from python. These functions allow us to call MERCURY’s time-stepping functions directly from python, avoiding any potential issues that may arise due to coordinate transformations, unit conversions, time-stepping logic, or input/output files. This might have been an issue in an earlier study looking at the symplecticity of MERCURY (see Appendix B). Aside from small differences at the floating point precision limit, if we use the same switching function as implemented in MERCURY, we find perfect agreement between MERCURY and MERCURIUS in all our results and therefore only show the results using MERCURIUS.

6 Results

Figure 2 summarizes our results for the test case and switching functions introduced above. The columns correspond to the different switching functions. The horizontal axis of every panel corresponds to time, with the close encounter occurring at $t\approx 0$. The particles enter the critical distance $r_{\rm crit}$ at $t\approx-1$ and exit at it again at $t\approx 1$. The top row shows the value of the switching function $L$ as a function of time. The second row shows the absolute value of the gravitational force between the planets as calculated during the interaction step. We also show the first four derivatives of the force. The third, fourth, and fifth rows show approximations of the third, fifth, and seventh order commutators, respectively. We colour commutators that only include $\widehat{H}_{K}$ and $\widehat{H}_{I}$ (i.e. not $\widehat{H}_{J}$) red, and all other commutators orange. We do not show commutators that are always zero (i.e. those only involving $\widehat{H}_{I}$ and $\widehat{H}_{J}$). As a measure of the size of the commutators, we use the energy error as described above. Note that the energy error does not account for phase errors. Thus an energy error of zero does not imply a perfect solution. The plots for the commutators also show the actual energy error occurred over one timestep step, as well as the integrated energy error, i.e. the energy error measured relative to the beginning of the simulation.

The first column shows the evolution for the standard WH integrator. The forces, its derivatives, and therefore all commutators get very large during the close encounter. As a result, the integrated energy error goes off the chart. As expected, this integrator cannot accurately resolve the close encounter. Note that the higher order commutators become comparable to the lower order commutators during the close encounter. The physical interpretation for this is that timestep is too large to resolve the characteristic timescale (periastron passage of the close encounter). Also note that the dominant commutators are those involving only $\widehat{H}_{K}$ and $\widehat{H}_{I}$.

The second column shows the evolution of using the Heaviside function as the switching function. Note that the commutators which only involve $\widehat{H}_{K}$ and $\widehat{H}_{I}$ vanish when $L=0$. This happens because $\widehat{H}_{I}$ becomes the identity operator, i.e. not having any effect, and it thus commutes with any other operator. Commutators involving $\widehat{H}_{K}$ and $\widehat{H}_{J}$ do not vanish when $L=0$ and become the dominant contribution to the energy error during that time. However, these are several orders of magnitude smaller and do not diverge during the close encounter unless the planets are very eccentric and close to periastron. Note that there is a spike in most commutators exactly when $L$ transitions at $r_{\rm crit}$. This can be understood by imagining a scenario where the particles start out just outside of $r_{\rm crit}$. To calculate the commutator $\small[\widehat{H}_{K},\widehat{H}_{I}\small]$, let us initially apply $\widehat{H}_{I}(dt)$. When calculating the effect of this operator the value of $L$ will be $1$. If we then apply $\widehat{H}_{K}(dt)$ to the result, it might take the particles inside of $r_{\rm crit}$. Thus, if we then apply $\widehat{H}_{I}(-dt)$, the value of $L$ will be [math]. Let us finally apply $\widehat{H}_{K}(-dt)$ again which might take the particles outside of $r_{\rm crit}$. This sequence of operators corresponds to the commutator $\small[\widehat{H}_{K},\widehat{H}_{I}\small]$ which is one of the nested commutators that appears at all orders in the BCH formula. The above argument makes it clear that the commutator will be particularly large when the switching function changes a lot during a timestep. This is the case for the Heaviside function, but it is important to note that for any other switching function that changes significantly during a timestep this is also the case, no matter if the switching function is smooth or has a discontinuity. The consequence is that the actual energy error (both the single step error and the integrated error) jump up when the particles transition $r_{\rm crit}$. Despite this undesirable jump, the integrator is at least somewhat better at resolving the close encounter than the standard Wisdom-Holman integrator. A physical argument can be given in terms of timescales again: if the switching function changes a lot during one timestep, we’re effectively introducing a very short timescale that the integrator is unable to resolve.

We can avoid large changes of the switching function during the timestep, by simply calculating the value of $L$ once, at the beginning of the timestep, and then keeping it fixed. The results of this modification to the Heaviside switching function are shown in the third column. Several features are identical to the standard way of using the Heaviside switching function (second column). In particular, all commutators that only involve $\widehat{H}_{K}$ and $\widehat{H}_{I}$ still vanish when $L=0$. However, we see that we have successfully removed the spikes in the commutators when $L$ changes. The commutators involving $\widehat{H}_{K}$ and $\widehat{H}_{I}$ simply go to zero. As a consequence the energy error does not increase when particles pass over the $r_{\rm crit}$ transition, despite the sudden change in the value of the switching function $L$. This can be explained following the same scenario described in the last paragraph. This is a remarkable result: This integration scheme, effectively using a binary switching function, is significantly simpler to implement than those using a continuous switching function. However, as one can see in the figure, it provides equal accuracy for one close encounter, at reduced computational complexity. We further comment on this case below.

The fourth column shows the evolution using the standard MERCURY algorithm and its polynomial switching function. Once again all commutators that only involve $\widehat{H}_{K}$ and $\widehat{H}_{I}$ vanish when $L=0$. However, in contrast to the integrators using the Heaviside function, this only happens for $r<0.1r_{\rm crit}$. Although these commutators are finite for $r>0.1r_{\rm crit}$ they are significantly smaller than in the case of the Wisdom-Holman integrator. Note that the third order commutators are smooth. This is because the specific choice of switching function ensures that the forces and its first two derivatives are continuous. These derivatives appear in the Poisson bracket and therefore in the commutator. However looking at the higher order commutators, one can see the effect of non-continuous higher derivates due to the finite differentiable switching function $L_{\rm merc}$. Fortunately, these commutators are significantly smaller than the third order commutators and therefore do not contribute much to the total energy error. Remember that the plots only show approximations of the commutators. This becomes apparent in the 7th order commutators where the spikes near the discontinuities of the force derivatives at $r_{\rm crit}$ are somewhat smeared out.

The fifth column shows the evolution using an infinitely differentiable switching function. Here, all derivatives of the force, and therefore all Poisson brackets, are continuous and smooth. Because of the smoothness requirement, the higher order derivatives of the force become highly oscillatory. This is a generic feature of any infinitely differentiable switching function and not specific to the function we chose. The commutators look similar to those of the MERCURY algorithm but one significant difference is that the higher order commutators become larger and more oscillatory. As in the previous case this is because these commutators are directly related to the higher order derivatives of the forces. For moderate timesteps like the one we’ve chosen in this example, the energy error of the simulation is still dominated by the third order commutators. However, for larger timesteps, there can be cases where the higher order commutators will dominate the error. We can think of this in physical terms once again. The highly oscillatory derivatives of the force corresponds to very short timescales. Because we are using a finite and fixed timestep, these timescales can at some point be no longer accurately resolved, leading to larger errors. The appearance of these small timescales is a purely numerical artefact (somewhat related to timestep resonances discussed by Rauch & Holman, 1999). Mathematically speaking, the differential equations corresponding to the higher order commutators become very stiff. Note that in the standard Wisdom-Holman integrator (left panel), the forces get large, but do not oscillate.

7 Conclusions

We revisited the motivation behind hybrid symplectic integrators and the choice of switching functions. We found that the derivation of equations of motion by Chambers (1999) omits derivative terms. As we showed, it turns out that despite this inconsistency the MERCURY integrator is a switching integrator. However, it effectively uses a very different switching function than the one described by Chambers (1999).

Motivated by the somewhat arbitrary choice of switching function in MERCURY, we explored a wide range of different switching functions to see if it is possible to improve the accuracy of this kind of integrators any further. In particular, we presented results of two extreme cases: one infinitely differentiable function, and the Heaviside function. In summary, we found that the smoothness alone is not a good criterion for choosing a switching function.

Our results show that an infinitely differentiable switching function does not perform much better over one close encounter than the polynomial function used by MERCURY. We attribute this to two reasons. First, the beneficial effects of a smoother function only show up in higher order commutators which are several orders of magnitude smaller than the dominant third order commutators. Second, the higher order derivates of an infinitely differentiable switching function necessarily become highly oscillatory. The integrator cannot resolve the associated small timescales and thus behaves no better than an integrator involving finite differentiable functions.

We show that using a Heaviside switching function can have surprisingly good properties if the value of the switching function does not change during the timestep. The resulting scheme is significantly simpler than the MERCURY scheme but achieves the same accuracy in the test case of a single close encounter presented. What still needs to be studied in greater detail is the long term evolution in systems with repeated close encounters. We ran additional tests of 1000 highly collisional planetary systems and found that the integrator using the Heaviside function performs on average as well as MERCURY. We will discuss the implementation of this new integrator, as well as the implementation of the MERCURIUS integrator, in an upcoming code description paper.

In practice one is mostly interested in how accurate an integrator can resolve a close encounter. In any numerical simulation, one can calculate a precise answer to this question using any arbitrary metric one thinks is approriate. A mathematical (or maybe even somehwat philosophical) question is whether these integrators we discussed in this paper are symplectic or not. The authors of this paper were not able to agree on the answer to this question. We will therefore leave the final answer up to the reader, but comment on a few aspects of this question.

First, it is important to keep in mind that we very quickly reach the limits of finite floating point precision in systems with close encounters555Imagine a scenario where two planets on orbits with $a\approx 1$ have a close encounter with a minimal encounter distance of one Earth radius, $\approx 10^{-5}$. Using double floating point precision, we only have $\approx 11$ decimal digits to work with when representing the planets’ positions in heliocentric coordinates. When calculating any derivative, which is required for checking the symplecticity, we further loose roughly half of the significant digits.. One could argue that we should not talk about concepts of differential geometry such as symplecticity at all if we work with floating point numbers because derivatives formally do not even exist. We run simulations on a discrete phase space where only finite differences exist.

Second, note that all switching functions we used are non-analytic. That is not a coincidence but a requirement if we want the integrators to fall back to the Wisdom-Holman integrator whenever particles are far away from each other. As a consequence, the Taylor series of the switching function, and therefore the Poisson Brackets in the BCH formula will not converge to the correct solution globally.

Third, when choosing a switching function, differentiability is not the right criteria to consider. If the mathematical concept of differtiability were an important requirement for a switching function, we could simply come up with a new infinitely differentiable function that is arbitrarily close to the original non-smooth function. We have effectivly done that for $L_{\rm merc}$ by introducing $L_{\rm inf}$, which did not lead to improvements in any of our tests. What matters is how much the switching function changes during a (finite) timestep. If it changes a lot, then the corresponding differential equations become stiff and the errors large. What also matters is that the switching function is continous to ensure the existance and uniquness of a solution666This only matters when calculating the solution of the $H_{K}$ step. It does not matter for the solution of the $H_{I}$ step because the switching function is constant. (Hernandez, 2019). Note that integrator with the Heaviside function completely changes the nature of the algorithm and can among other things break time reversibility. This might be important for some long term integrations, but not for simulations where encounters are rare (for example for planetary system which go unstable), or if physical collisions break time reversibility anyway.

Fourth, all integrators will converge to the correct solution in the limit of the timestep going to zero. We thus argue that the simplest and most reliable way to test that a given integrator gives reliable answers is to change the timestep and watch for changes in the results, i.e. a classical convergence test.

Acknowledgments

We thank the referee John Chambers for providing us with helpful comments which allowed us to improve and clarify the manuscript. We thank Scott Tremaine for helpful discussions at various stages of this project. This research has been supported by the NSERC Discovery Grant RGPIN-2014-04553, the 2018 Fields Undergraduate Summer Research Program at the Fields Institute for Research in Mathematical Sciences in Toronto, and the Centre for Planetary Sciences at the University of Toronto Scarborough. Support for this work was also provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51423.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This research was made possible by the open-source projects Jupyter (Kluyver et al., 2016), iPython (Pérez & Granger, 2007), and matplotlib (Hunter, 2007; Droettboom et al., 2016).

Appendix A Democratic Heliocentric Coordinates

The hybrid integration schemes discussed in this paper make use of democratic heliocentric coordinates (Duncan et al., 1998). Let us assume we have one star and an additional $N$ planetary bodies with Cartesian coordinates $\mathbf{q}_{i}$ and canonical momenta $\mathbf{p}_{i}$ in some inertial frame. Then, the democratic heliocentric coordinates are defines as:

[TABLE]

Here, $M=\sum_{j=0}^{N-1}m_{j}$ is the total mass of the system. The corresponding canonical momenta are

[TABLE]

Appendix B Symplecticity of Mercury

While working with MERCURY and MERCURIUS we also looked at measuring the symplecticity error. To do that we follow Hernandez (2016) in calculating the Jacobian using a finite difference approach. This is non-trivial because one quickly runs into issues of finite floating point precision. The results of Hernandez (2016) suggest that the MERCURY algorithm might be non-symplectic for their binary planet test case. To investigate this issue, we implemented wrapper functions to directly call the MERCURY subroutines which evolve the state by one timestep, but ignore all the input, output and other bookkeeping login in the MERCURY package. Our results with this method differ from those obtained by Hernandez (2016) and indicate that MERCURY is symplectic (whether we should really call the scheme symplectic or not, becomes a somewhat philosophical rather than practical question, see discussion).

We attribute this discrepancy to some non-symplectic operations that the MERCURY package applies to the state vector before and after the actual integration (e.g. unit conversions and coordinate transformation), and a loss of precision when using ASCII files as input and output files. We also ran tests with our own implementation of the algorithm, MERCURIUS, and confirm the results obtained with MERCURY.