Towards time symmetric N-body integration
Walter Dehnen
astro-ph.IM

[NO \title GIVEN]
    [NO \author GIVEN]
    December 30, 2023
======================


	Computational efficiency demands discretised, hierarchically organised, and individually adaptive time-step sizes (known as the block-step scheme) for the time integration of N-body models. However, most existing N-body codes adapt individual step sizes in a way that violates time symmetry (and symplecticity), resulting in artificial secular dissipation (and often secular growth of energy errors). Using single-orbit integrations, I investigate various possibilities to reduce or eliminate irreversibility from the time stepping scheme. Significant improvements over the standard approach are possible at little extra effort. However, in order to reduce irreversible step-size changes to negligible amounts, such as suitable for long-term integrations of planetary systems, more computational effort is needed, while exact time reversibility appears elusive for discretised individual step sizes.



	gravitation —
	methods: numerical —
	celestial mechanics





§ INTRODUCTION

Astrophysical N-body problems typically have a large range of dynamical time scales with factors of 10^2-4 between the shortest and longest orbital times. Consequently, instead of using time steps of (fixed or varying) global size, contemporary N-body algorithms advance each particle with time steps of a size which is individually adapted along its trajectory <cit.>. There are two components to such an individual time-stepping method: a time-step function and a time-stepping scheme. The time-step function (or criterion) returns an appropriate step size given the instantaneous state of a particle's trajectory. The time-stepping scheme, on the other hand, is a method that adapts the individual particle step sizes to follow these time-step criteria as best as possible.

This study is concerned solely with the second ingredient, the time-stepping scheme. There are two important conditions for such a scheme: (1) it should not hinder computational efficiency and (2) it should be time reversible (and/or support symplectic time integration[A symplectic integrator advances the 
trajectories by a canonical map. An equivalent statement is that the Jacobian

    J⃗ = ∂ξ(t+h)/∂ξ(t)

between the initial state ξ={x⃗,p⃗} of the system and that advanced by step size h satisfies J⃗^T·Ω⃗·J⃗=J⃗ with the symplectic matrix

    Ω⃗ =[  0⃗ -I⃗;  I⃗ -0⃗ ].

As a consequence, the geometric structure of phase space, most first integrals, and the Poincaré invariants are preserved and the energy error tends to be bounded, but see footnote <ref>.]). This latter condition is important to avoid artificial numerical dissipation <cit.>. In order to meet the first condition, all contemporary astrophysical N-body methods for large N employ the block-step method, where particle time-step sizes are discretised and hierarchically synchronised <cit.>. The original motivation for this scheme was the reduction in the number of predictions of particle positions, which are required for the computation of the forces on other particles. Moreover, with modern gravity solvers the simultaneous computation of all  gravitational force between N particles requires only 𝒪(Nln N) <cit.> or (fewer than) 𝒪(N) operations <cit.>, instead of 
𝒪(N^2) for a direct force summation, which allows considerable efficiency saving from the synchronisation. It appears that the block-step is the only possibility to achieve these savings and yet allow for individual time-step sizes. It is therefore mandatory to use this method.




However, (as far as I am aware) none of the contemporary N-body codes employing the block-step is time reversible or symplectic, and only little effort has been made towards that goal <cit.>. One problem is that integrating the mutual forces between any two particles with different step sizes for either particle cannot be reconciled with a canonical map and hence  symplecticity. Note, however, that the block-step method itself does not destroy symplecticity[Any dependence of h on ξ alters the Jacobian (<ref>) and destroys symplecticity <cit.>. However, when using discrete step sizes (as with the block-step scheme), the function h(ξ) is piece-wise constant with jump discontinuities. Thus, the neighbourhood of almost all trajectories use the same h and the integration remains symplectic (Tremaine, private communication, 2016). The exceptions (trajectories hitting the jumps in h) occupy a volume of measure zero in phase space and are irrelevant.].  (, see their Fig. 5) demonstrate this with an example where the Poincaré invariant is conserved to within round-off error using the block-step scheme.

Symplecticity can be restored by integrating the force between any two particles with the same time-step size (either that of the faster or slower of the two particles), resulting in exact conservation of total momentum <cit.>. However, such schemes are only reversible as long as the time-step adaption process (the application of the jumps discussed in footnote <ref>) is. 

Since the N-body dynamics is strictly reversible, irreversibility of the numerical integration method tends to result in artificial dissipation. As a consequence, the energy error (for example) is not guaranteed to be bounded, irrespective of whether the integrator is symplectic[Symplecticity by itself does not imply time reversibility. A simple counter example is the second-order accurate leapfrog integrator followed by a positional offset proportional to h^4. This is a symplectic, second-order accurate integrator, yet is not reversible and suffers from significant energy drift.] or not. This consideration suggests that reversibility of the integration scheme is more important than symplecticity. Unfortunately, so far no reversible yet efficient block-step-based time-stepping scheme is known. In fact, most practitioners determine the step size by the time-step function evaluated at the start of the time step <cit.>. This simple forward method violates time symmetry whenever the step size is adapted and is well known to destroy long-term stability <cit.>.
Thus, the most wanted ingredient for reversible N-body integration is a reversible method for adapting the individual particle step sizes given some time-step function.

The goal of this study is to consider ways to improve this situation. If h is not discretised, but continuous, the situation is much simpler and various adaptive time stepping methods have been proposed <cit.>. Here, I adapt several of these to discrete h and the block-step and study them in the context of single-orbit integrations. This constitutes a first test that any such method must pass before being considered for the full N-body problem. In particular, I do not toy with the time-step function or the underlying time integrator.

The paper is organised as follows. In section <ref> the problem is stated more formally, while sections <ref> and <ref> present two different approaches of improving reversibility. Section <ref> presents a third approach, which is more accurate but also more computationally expensive and hence more suitable for long-term integration of planetary systems. Finally, sections <ref> and <ref>discuss the findings and conclude.









§ STATEMENT OF THE PROBLEM


Let ξ denote the state of a particle orbit and ϕ_h a self-adjoint (time-symmetric) integrator used for individual time steps (i.e. the composition ϕ_-h∘ϕ_h obtains the identity map). Then the state of the integrated orbit after n steps of sizes h_i is

    ξ_n = ϕ_h_n-1/2 ξ_n-1
    		    = ϕ_h_n-1/2∘…∘ϕ_h_1/2 ξ_0,

where h_n+1/2=t_n+1-t_n. Thus, the step sizes naturally have half-integer indices to satisfy time symmetry. Since ϕ_h is reversible, so is the combined map (<ref>) iff the individual step sizes h are reversibly adapted. With the block-step scheme individual step sizes are discretised as

    h = 2^-r h_max

with integer rung r≥0 and arranged hierarchically, as shown in Fig. <ref>. As a consequence of the hierarchical structure, a change to higher rung (δ r>0: shorter step size) is always possible, but a change to lower rung (δ r<0: longer steps size) only if the synchronisation requirement is met (which is every 2^-δ r steps). Furthermore, the change in rung is usually limited to |δ r|≤1, corresponding to changing h by a factor two either way. These block-step synchronisation constraints imply that given a desired step size τ, there exist at any time t_n a unique block-step maximum step size

    h_block(τ,t_n) ≤τ

of the form (<ref>). The dependence of h_block on the actual simulation time t_n originates only from the block-step synchronisation constraint, i.e. the fact that doubling the step size (reducing the rung) is not always possible. In view of Noether's theorems, this dependence of the integration method on absolute time suggests that the energy error may never be fully bounded with the block-step method. However, one would expect this time dependence to affect the energy not in a systematic but rather pseudo-random way. In fact, when experimentally dropping the block-step synchronisation constraint in single-orbit integrations, I found if anything a deterioration of the energy conservation. In the notation used in the remainder of this paper, the time-dependence of h_block is suppressed for brevity.

If the integrator ϕ_h requires a force at the end (and start) of each time step, then the force computations are synchronised not only between particles with the same rung but between particles of all rungs larger than the smallest active rung at the given time (see also Fig. <ref>). This is the main reason why most astrophysical N-body methods for collision-less dynamics use the kick-drift-kick (rather than the drift-kick-drift) version of the leapfrog integrator

    ϕ_h = ϕ^K_h/2∘ϕ^D_h∘ϕ^K_h/2

with kick and drift operators

    ϕ^K_h∘(t, x⃗, v⃗)   =   (t, x⃗, v⃗-h∇⃗Φ(x⃗,t)),
    
    
    	
    ϕ^D_h∘(t, x⃗, v⃗)   =   (t+h, x⃗+hv⃗, v⃗),

where x⃗ and v⃗ denote position and velocity, respectively. All tests in this study are performed using this integrator, but most conclusions regarding the effect of the time-step adaptation schemes are also valid for other integrators, including those of higher order.

The statement of the problem then is:

    ∙  adapt h reversibly and such that h≈ T(ξ),

where T(ξ) denotes the time-step function. This problem is non-trivial because the states ξ, and hence T(ξ), are not synchronised with the step sizes h, as indicated by the half-integer indices for the latter. Thus, in order to solve the problem, some form of synchronisation of relation (<ref>) is required, and different synchronisation attempts result in different stepping schemes.




 §.§ What is wrong with the simple forward method?

As mentioned in the introduction, the state of the art for adapting individual step sizes is simply

    h_n+1/2 = h_block(T_n)
    	  with   T_n≡ T(ξ_n),

when the step size matches (block-step allowing) its optimal value at the start of the step. This simply ignores the synchronisation problem in equation (<ref>). As a consequence, any time dependence of T result in an O(h) synchronisation error in equation (<ref>). The forward method gives the adaptation scheme

    δ r
    	= {[                     -1                     if                h_n-1/2                      ≤                 12 T_n and block-step allows,;                      0                     if                h_n-1/2                      ≤                   T_n,;                     +1           4lotherwise. ].

Here, I have limited rung changes to |δ r|≤1, which is common practice and which is assumed in the remainder of this study unless otherwise stated.

For two orbits integrated with this scheme, Fig. <ref> plots the time evolution of radius, h, T, and the energy over one orbital period. The time-step function used is

    T(ξ) = η√(r^3/GM(r))

with η=0.02 and M(r) the mass enclosed at radius r (of course, such a time-step function is not directly available in N-body simulations, but reasonable approximations are, e.g. ). Evidently, the energy varies considerably over one orbit. Such oscillations of energy are characteristic of symplectic integrators. If using a (sufficiently small) fixed step size, the error made on approach to peri-centre is exactly undone on the way out again.




With the adaptive method, shorter time steps near peri-centre reduce the energy error there and improve the overall trajectory accuracy, while avoiding unnecessarily short steps for most of the orbit. However, the error is not exactly symmetric w.r.t. peri-centre. Instead, with the forward adaptation method, the integration accuracy is on average higher on the outward than on the inward part of the orbit and a residual error remains. For
an integrator of order n (n=2 for the leapfrog used here), the residual error is proportional to h^n+1∝η^n+1 per irreversible step-size change. As a consequence, the energy error grows over the long term.

This can be seen in Fig. <ref>, which plots  for an integration over 10^4 periods the ratio of the long-term energy error over the short-term energy error (maximum error over a single orbit as shown in Fig. <ref>) for η=0.02 (black) and 0.01 (red). Evidently, δ E_long/|δ E_short| grows roughly linearly at a rate ∝η. Thus, unless η is chosen exceedingly small the long-term error (for these eccentric orbits) is soon dominated by the adverse effects of irreversible step-size changes. This result is independent of the order of the underlying integrator ϕ_h.




 §.§ Counting stepping errors


The departure of the energy error from exactly symmetrical behaviour over each orbit, and hence the net energy error and long-term drift, are a consequence of each step-size change being irreversible, which results in deviations of the step size h from the condition (<ref>). While the forward scheme does not avoid such deviations, one can detect them a posteriori.

There are four types of deviations, depending on the sign of Ṫ and whether h was too large or too small. One may estimate the ideal step size as T_n+1/2 if known (i.e. for a longer step when in fact two short steps have been taken), otherwise as √(T_nT_n+1). Then a step is too long if h_n+1/2>√(T_nT_n+1) and two adjacent steps of the same size h are too short if a larger step was possible and 2h<T_middle.

The forward scheme (<ref>) only makes two types of errors: too long steps at Ṫ<0 and too short steps at Ṫ>0. Both tend to result in an energy error of the same sign (which is different for the two orbits considered in Fig. <ref>).

Errors of the other types (too long steps at Ṫ>0 or too short steps at Ṫ<0) did not occur, but would have resulted in energy errors of the opposite sign. I define the net rate of stepping errors as difference in the number of errors of different sign divided by the total number of step changes:

    [                               R_err≡                                    [          number of too long steps at                                  Ṫ<0;                                                                         +         number of too short steps at                                  Ṫ>0;                                                                         -          number of too long steps at                                  Ṫ>0;                                                                         -         number of too short steps at                                  Ṫ<0                                  ] /;                                                                           2ltotal number of step-size changes, ]

where a pair of adjacent too short steps is only counted once. For the forward scheme, the chance for a stepping error is only 50%, in agreement with the results from the orbit integrations (indicated in Fig. <ref>). Of course, R_err is a rough measure, as it only accounts for the typical sign but not the actual size of the errors. In reversed time, Ṫ changes sign and hence also R_err. Therefore, a reversible scheme should obtain R_err=0, even though it may not always chose the ideal step size. Such a scheme should have long-term energy errors comparable to the short-term energy error.


Since the short-term error is ∝ h^n∝η^n for an integrator of order n, the ratio of long-term errors owed to irreversible step-size changes to the short-term error scales like η R_err. For the forward method, R_err∼0.5 independent of η, and therefore δ E_long/|δ E_short|∝η, in agreement with Fig. <ref>. Of the two orbits presented, the Kepler orbit has three times as many step-size changes than the cusp orbit, but owing to the bimodal structure of the error (bottom left panel of Fig. <ref>), some of the resulting energy errors cancel, unlike the situation for the cusp orbit. Therefore, the net effect on the ratio δ E_long/|δ E_short| happens to be similar for these two orbits.





 §.§ How can the problem be solved?





  §.§.§ Reducing or eliminating stepping errors

If one can reduce or even eliminate stepping errors, or at least their net rate such that errors of opposite sign largely cancel, the integration will be near-reversible, even if no attempt is made at constructing an exactly reversible scheme.

One possibility is to extrapolate the time-step function into the future. In section <ref>, I consider such a method which has synchronisation error 𝒪(h^3), i.e. as good as the leapfrog's trajectory error, but still incurs stepping errors.

A more rigorous approach is to attempt to eliminate stepping errors by trying different h and then chosing those satisfying the time-step criteria. Applying this approach to each particle separately ignores the interdependence of particle orbits and time-step functions, but is close enough to reversibility for most practical purposes (see section <ref>). However, this approach requires twice as many force computations as actually used. This overhead can be reduced (but not eliminated) when combined with the extrapolation method, see section <ref>.




  §.§.§ Aiming for reversibility


When instead trying to obtain truly reversible schemes, past and future (step sizes) must be treated symmetrically, which severely restricts our hands in how to use our knowledge about the past. In order to avoid computationally expensive iterations, a reversible method must be explicit, i.e. the future step size must be obtained from the current value of the time-step function and the past step size in a reversible way. Since h is discretised, the choice for the future is always between a long step (either equal to or twice the previous step size) and two shorter steps. If the decision for the next step size is based solely on the first of these shorter steps, time symmetry is violated, while the second shorter step is beyond the horizon of prediction (at the moment of the decision). Therefore, in order to maintain time symmetry, the decision must be based solely on the merit of the long step: an explicit scheme must prefer longer steps. Attempts to obtain such schemes are presented in sections <ref> and <ref>.




§ EXTRAPOLATING AND ITERATING


For the situation of a continuous step size h, formally exact time symmetry can be obtained by matching h to some mean of the time-step function T at the start and end of the step:
[
	Matching h instead to the time-step function evaluated in the middle of the step 
	
    h_n+1/2 = T_n+1/2≡ T(ϕ_1/2h_n+1/2ξ_n)

	fails to obtain exact time symmetry, because in general
	ϕ_h/2≠ϕ_-h/2∘ϕ_h, such that T_n+1/2 obtained in the forward and backward directions differ.]

    h_n+1/2=μ(T_n,T_n+1)
    	  with  
    	T_n+1≡ T(ϕ_h_n+1/2∘ξ_n).

Here, μ(x,y) denotes a general mean, satisfying μ(x,y)=μ(y,x) and min{x,y}≤μ(x,y)≤max{x,y}, for example the arithmetic mean

    h_n+1/2=12(T_n+T_n+1)

<cit.>, the geometric mean

    h_n+1/2=√(T_nT_n+1),

or the harmonic mean. Equation (<ref>) is a non-trivial implicit relation for h_n+1/2, which requires an iterative approach for an exact solution.

For h discretised with the block-step, equation (<ref>) is naturally replaced by 

    h_n+1/2=h_block(μ(T_n,T_n+1)),

when the step size can only take a few discrete values, such that a finite number of iterations suffices for full convergence. However, owing to the interdependence of the particle trajectories, an exact solution requires iterating not just over the trajectory of each particle individually, but over the N-body trajectory of all particles combined. While this appears to give a reversible scheme <cit.>, it requires an enormous computational effort both in time and memory, and is completely unpractical.







 §.§ Extrapolating the time-step function


One possibility to avoid such iterations is to estimate T_n+1 and obtain an approximately reversible scheme. The lowest-order approximation is simply T_n+1≈ T_n, when equation (<ref>) results in the forward scheme (<ref>). At the next order, one can use the previous value of T and estimate

    √(T_nT_n+1)≈ 
    	T_n (T_n/T_n-1)^2 h_n+1/2/2 h_n-1/2,

which amounts to linear extrapolation of ln(T) in time. When inserted into equation (<ref>) this still gives an implicit relation for h_n+1/2, but with a trivial dependence. For the block-step with |δ r|≤1, h_n+1/2 can take only three allowed values and one can easily solve this relation to obtain the scheme

    δ r = {[                                 -1                 if   h_n-1/2 T_n-1                                  ≤ 12T_n^2
    		  and block-step allows,;                                  0               if   h_n-1/2^2 T_n-1                                  ≤                             T_n^3,;                                 +1                         otherwise. ].

Strictly, these conditions are not unique if h_n-1/2≥2T_n. While this should never occur for appropriate time-step functions, it can be easily resolved by testing for longer steps first, say.

Equation (<ref>) has an error ∝ h^2 (^2ln(T)/ t^2), but since even-order time derivatives are time symmetric, the synchronisation error of scheme (<ref>) is 𝒪(h^3) as opposed to 𝒪(h) for the forward scheme.

Fig. <ref> shows the ratio of long- to short-term energy error over 10^4 periods with this scheme for exactly the same two orbits already considered above with the naive forward scheme. The net rate R_err of stepping errors is much smaller than for the forward scheme (<ref>), representing a significant improvement.
As a result, the long-term and orbital energy variations are comparable in stark contrast to the situation for the forward scheme (see Fig. <ref>). Another significant improvement over the forward method is that R_err decreases with η. As a consequence, the ratio of irreversibility-induced long-term energy errors to the orbital energy error decreases faster than ∝η.

In this context, it is worth noting that not all of the measured long-term energy error is due to irreversible step-size changes (it also occurs for integrations without any violation of time symmetry, such as for the Kepler orbit integrated with η=0.01 in Fig. <ref> below). As already discussed in section <ref>, Noether's theorems and the dependence of h_block(τ) on absolute time suggest such pseudo-random behaviour of the long-term energy error. Alternatively, one may interpret this as a random walk caused by an incommensurability between period and time stepping (the cusp orbit exhibits some quasi-periodic behaviour of the error, clearly visible at η=0.01 in Figs. <ref>-<ref>, <ref> & <ref>, indicating a resonance between period and time stepping).





 §.§ A try-and-reject scheme



The extrapolation of the time-step function is only advisable if it is well-behaved, i.e. does not change much on time scales of itself, implying |Ṫ|≪1. This is not necessarily the case in N-body simulations, where fluctuations may occur due to close and insufficiently softened encounters as well as force approximation errors.

Therefore, I now consider an alternative approach to deal with the implicit relations. If one ignores the violation of time symmetry resulting from the mutual dependence of the particle trajectories, one can consider each particle trajectory separately. Since there are only few discrete values allowed for h_n+1/2, one can, instead of iterating, simply try them one after the other, starting with the shortest allowed step.

To this end, the position at t_n+12h_n-1/2, corresponding to halving the step size, is predicted first and the time-step function T at that point computed. If h_n-1/2>T the step size h_n+1/2=h_n-1/2 is too long (in the sense of equation <ref>) and instead the shorter step must be used, the first part of which has already be done by predicting the position. Otherwise, the step is continued to size h_n+1/2=h_n-1/2 at the end of which the procedure is repeated (subject to block-step synchronisation constraints).

In order to continue the step, the kick-drift-kick leapfrog can be re-phrased as a predictor-corrector scheme. The predictor simply predicts position and velocity into the future assuming a constant acceleration a⃗=a⃗_beg equal to that at the start  of the step. The prediction operation is associative and second-order accurate in position. A full time step consists of any number of predictions followed by a force computation (obtaining a⃗_end) and the correction step

    v⃗←v⃗ + 12 h 
    			(a⃗_end-a⃗_beg).

In practice, this scheme requires storage for a⃗_beg and h, the accumulated length of the step so far (which is incremented by the predictor). The scheme for a full time step for one particle is then as follows. 

	
  (i) Increment the rung (r← r+1) to try a shorter step.	
  (ii) 
		Predict ξ until the rung r is block-step synchronised.
	
  (iii) Calculate T and the force (obtaining a⃗_end)[Strictly speaking, a⃗_end is only required in (v), but since T typically depends on the gravitational potential and/or acceleration, the efficiency saving by evaluating a⃗_end only when needed is not available.].
	
  (iv) If 2h≤ T and if rung r-1 is block-step synchronised at time t+h, reject the step (discarding the force computation just made), decrement the rung (r← r-1), and recurse with point (ii).
	
  (v) Otherwise, finish the step by applying the correction (<ref>), then reset a⃗_beg←a⃗_end and h←0. 





This scheme is similar in spirit to that proposed by <cit.>, with the main difference that their scheme tests for longer steps first, while I test for shorter steps first. I should also note that this scheme violates time symmetry not only because of the mutual dependence of particle step sizes, but even when applied to a single orbit in a fixed potential. This is because the predictor is only time symmetric between the start and end of a step, but not with respect to intermediate times. The position predicted for the middle of a step (at unchanged step size) differs between forward and backward steps by 18h^2(a⃗_end-a⃗_beg), see also footnote <ref>. Consequently, the time-step function T evaluated at these positions  in (iii) differs too and with it possibly the decision about the step size in (iv).

However, these violations of time symmetry are not obviously systematic and do not introduce a coherent arrow of time. Hence, no significant systematic energy drift is to be expected. Fig. <ref> shows that for the same two orbits used before the ratio of long- to short-term energy error is slightly smaller than for the extrapolation scheme form the previous sub-section.




 §.§ Combining extrapolation with try-and-reject





The try-and-reject scheme is rather wasteful, as it requires (on average) twice as many force computations as it actually uses for the orbit integration. However, by combining extrapolation of the time-step function with the rejection idea, one may avoid to always try for a shorter step. One can parameterise the ignorance about the precise future values for T by a parameter λ≤1 and assume that

    T_n+1/2≥ T_n λ^2 h_n+1/2/2 h_n-1/2,

which corresponds to the assumption that Ṫ≳ln(λ). With this  
supposition, one can replace step (i) in the algorithm given in the previous sub-section by r ← r + δ r with

    δ r = {[                                           -1 if   h_n-1/2≤12T_n λ  and block-step allows,;                                           -0                     if   h_n-1/2^2 ≤ T_n^2λ,;                                           -1                                   otherwise. ].

For λ=0.8, Fig. <ref> shows the evolution of the ratio of long- to short-term energy error for the same orbits as considered before. The long-term error is similar to that for the try-and-reject scheme before (Fig. <ref>). For these orbits, the ratio of force calculations to time steps is down from 2 for the pure try-and-reject scheme to 1.19 and 1.36 for Kepler and cusp orbit, respectively.




§ ATTEMPTING EXPLICIT REVERSIBILITY


The schemes discussed in the previous section attempt to obtain T_n+1/2, i.e. synchronise the right-hand side of equation (<ref>). Instead, the schemes considered in this section attempt to synchronise the left-hand side of equation (<ref>). For continuous step sizes (unconstrained by the blockstep scheme), this can be accomplished by swapping the roles of h and T in the implicit relation (<ref>), yielding <cit.>

    μ(h_n-1/2,h_n+1/2) = T_n,

which is an explicit relation for h_n+1/2, i.e. requires no iterations for an exact solution. Since in reversed time h_n-1/2 and h_n+1/2 are swapped while T_n remains unaltered, this scheme is reversible.  proposed the harmonic mean, 2/μ=1/x+1/y. I will use the geometric mean μ^2=xy, which of the one-parameter family of averages

    [           μ^α             =   12[x^α+y^α]      for α≠0,;           lnμ             = 12[ln x+ln y]       for α=0 ]

is the only one that always obtains a unique positive h_n+1/2 for any given T_n, h_n-1/2>0. In the remainder of this study, I stick to the geometric mean, but I experimented with various α without finding α=0 inferior.

The resulting integration schemes for continuous h give long-term energy errors consistent with zero for the two test orbits used before (not shown). It seems expedient, therefore, to try to transfer this idea to the block-step scheme.

For the further quantification of (the lack of) reversibility, I introduce the net rate of irreversible step-size changes in close analogy to the net rate of stepping errors:

    [                         R_irr≡                              [           number of steps with                       h>h_back                         at Ṫ<0;                                                             +           number of steps with                       h<h_back                         at Ṫ>0;                                                             -           number of steps with                       h>h_back                         at Ṫ>0;                                                             -           number of steps with                       h<h_back                         at Ṫ<0                            ] /;                                                               4lnumber of step-size changes. ]

Here, h_back is the step size which the same scheme would have taken when applied after the step of size h in the backward direction. Like stepping errors, irreversibilities come in four flavours, depending on the sign of Ṫ and the relation between h and h_back. Note that R_irr is an exact measure, unlike R_err of equation (<ref>), which involved an estimation for the correct step size.









 §.§ Reversible adaptation of a continuous time step variable



A direct porting of equation (<ref>) to discrete step sizes obtains a scheme only slightly better than the forward scheme, as described in the appendix. The problem is that for certain intermediate values of T the discrete step size h flips between two values bracketing T. These flips are reversible only if step-size changes are unlimited, which allows the possibility of arbitrary long steps.

To avoid too long steps, one must not allow this flipping of h. This can be achieved by following the initial scheme (<ref>) with an auxiliary continuous variable τ instead of h and set

    h_n+1/2 = h_block(τ_n+1/2).

However, implementing this directly via 

    τ_n+1/2=T_n^2/τ_n-1/2

(for μ the geometric mean) fails. The problem is that τ oscillates around T with amplitude increasing at each change in step size[
	This increasing oscillation amplitude seems at odds with the reversible nature of the method. However, this can be understood in analogy to linear instability where time symmetry implies that each growing mode is accompanied by a decaying mode of no practical significance.]. If h remains unchanged, the relation (<ref>) places ln(τ_n-1/2), ln(T_n) and ln(τ_n+1/2) on a line when plotted against time (rather than step index). So, to avoid the amplification of oscillations when h changes, one may simply retain that linear relation, which gives

    τ_n+1/2=T_n(T_n/τ_n-1/2)^h_n+1/2/h_n-1/2

instead of (<ref>). Together with equation (<ref>), this is an implicit relation for h_n+1/2 and τ_n+1/2 given the previous values and T_n. Again, solving these relations is straightforward, obtaining the scheme

    δ r = {[                     -1                     if       h_n-1/2τ_n-1/2^2                      ≤               12 T_n^3 and block-step allows,;                      0                     if         h_n-1/2τ_n-1/2                      ≤                 T_n^2,;                     +1           4lotherwise. ].

The resulting ratio of long- to short-term energy error for the two test orbits is shown in Fig. <ref> when started at apo-centre and Fig. <ref> when started half way between apo- and peri-centre. In the latter case, one sees a significant difference between an orbit integration started with τ_1/2=T_0 or τ_1/2=√(T_0T_1) (and h_1/2=h_block(τ_1/2) in both cases). This difference can be understood by a significant number of irreversible steps (as reported via R_irr).

Fig. <ref> shows the mechanics of irreversible step-size changes with this scheme. They an occur for both signs of Ṫ and originate from an integration/extrapolation error of τ, which causes the value for τ in the middle of a long step to be different in the backward direction, if that step was not taken. The chance for such integration errors and hence irreversible step-size changes is greatly enhanced if the auxiliary quantity τ oscillates around the time-step function (as it does to some extent in Fig. <ref>). If τ_n-1/2=T_n-1/2, relation (<ref>) corresponds to a linear extrapolation of ln(τ)=ln(T) in time. But if τ_n-1/2≠ T_n-1/2 or if ln(T) has curvature, this extrapolation goes slightly wrong and causes step-to-step oscillations of τ around T, like those detected by  for non-discrete h. This explains why carefully choosing the initial τ greatly reduces the number of irreversible steps (and with it the energy error) reported in the bottom panels of Fig. <ref>.








 §.§ Can the situation be improved?



However, this reduction of irreversible steps through careful choice of the initial τ works only as long as the time-step function is sufficiently smooth. This is not necessarily the case in N-body simulations, where T(ξ) may fluctuate. 


In order to reduce the number of irreversible step-size changes, one may devise techniques to suppress these oscillations, even though this inevitably introduces another type of irreversibility. One technique is to re-align τ=T when log_2τ is close to a half-integer, i.e. as far away from any step-size changes as possible. The resulting test orbit integrations give comparable results to those obtained with a careful choice of the first step size (Fig. <ref>, bottom panels).



Another idea is to avoid irreversible step-size changes by means of the try-and-reject method of section <ref> (which violates time reversibility in a more subtle and less harmful way than irreversible steps).
Because one cannot detect the violation of reversibility in time to avoid it, the only possibility is to try avoiding stepping errors. Then, instead to prefer the longer step in case of ambiguity (as per the argument of section <ref>), one tries the shorter step but rejects it if the longer step is found acceptable. However, this introduces a new type of irreversible rung change, because not all stepping errors of the scheme (<ref>) are irreversible.




§ INTEGRATING THE STEP SIZE


The time stepping schemes considered in the previous two sections attempt to synchronise the right- and left-hand sides of h=T (equation <ref>), respectively. I now consider schemes which instead synchronise a differential form, like 


    Δ h^-1   =   Δ T^-1    or
    Δln(h)   =   Δln(T).


The left-hand and right-hand sides of these relations can be expressed as finite difference and differential, respectively. For non-discrete step sizes, this obtains the schemes

	

    1/h_n+1/2 - 1/h_n-1/2
    	= - Ṫ_n/T_n

<cit.> or

    ln h_n+1/2 - ln h_n-1/2
    	= Ṫ_n.


According to  and my own experiments, these schemes perform very well for continuous step sizes (no block-step): even for chaotic orbits h closely follows T. Most importantly, this scheme avoids step-by-step oscillations of h, and hence is unlikely to suffer from any significant number of irreversible step-size changes, when adapted to the block-step.

The price to pay is the need to compute the time derivative Ṫ of the time-step function. This is usually not efficiently possible in large-N simulations, but is a reasonable option for small-N methods.




Since the block-step does not allow incremental changes to h, one cannot use equations (<ref>) directly, but must use an auxiliary variable τ. One must then also account for the discrete step over which τ is integrated, giving
 
	
    1/τ_n+1/2   =   1/τ_n-1/2
    	 - 	h_n-1/2+h_n+1/2/2Ṫ_n/T_n^2,
    	
    lnτ_n+1/2   =   lnτ_n-1/2  + h_n-1/2+h_n+1/2/2Ṫ_n/T_n

in place of equations (<ref>).

Combined with equation (<ref>), one again has trivially solvable implicit relations for h_n+1/2. When restricting rung changes to |δ r|≤1 as usual, this obtains the stepping schemes
 

    δ r    =   {[                     -1                     if      2/τ_n-1/2 -3d/T_n                      ≤              1/h_n-1/2 and block-step allows,;                      0                     if       1/τ_n-1/2 -d/T_n                      ≤              1/h_n-1/2;                     +1           4lotherwise, ]. 
    δ r    =   {[                     -1                     if     ρ_n-1/2-3/2(d/ln2)                      ≤              r_n-1/2-1 and block-step allows,;                     -0                     if     ρ_n-1/2-3/2(d/ln2)                      ≤               r_n-1/2,;                     +1           4lotherwise. ].


where d ≡ h_n-1/2(Ṫ_n/T_n), while ρ≡-log_2(τ/h_max) is a continuous rung variable replacing τ. To initialise τ, one must integrate the first half time step, i.e.
 

    1/τ_1/2   =   1/T_0 - h_1/2Ṫ_0/2 T_0^2,
    	
    τ_1/2   =    T_0 exp(h_1/2Ṫ_0/2T_0).


Fig. <ref> shows the run of the ratio of long- to short-term energy error for this scheme and the two test orbits. There are still very few irreversible step-size changes, which can occur for either sign of Ṫ. These are caused by the same basic mechanism that is also responsible for irreversible step-size changes with the scheme (<ref>), see Fig. <ref>. However, because the time integration of τ is more accurate and avoids oscillations, the schemes (<ref>) only suffer very occasionally from  irreversible rung changes. For the e=0.9 Kepler orbit, for example, no irreversible step-size changed occurred for η=0.01 over 10^4 orbits (the small-amplitude trend of the long-term energy for this orbit must be owed to other sources, see also the last paragraph of Section <ref>).

One can differentiate the relation between h and T once more and numerically integrate a second-order differential equation for τ, for example using the kick-drift-kick method. This obtains τ_n+1  reversibly predicted from τ_n and τ̇_n and enables the time-step condition h_n+1/2≤μ(τ_n,τ_n+1). Such an approach is conceptually very similar to a zonal time-step function T=T(x⃗), discussed in section <ref> below. It is quiestionable, however, whether this significantly reduces the chances of irreversible step-size changes, since the integration schemes (<ref>) are already second-order accurate.








§ DISCUSSION


The usage of individual adaptive time-step sizes in N-body simulations can result in enormous efficiency savings. Without these, hardly any of the many computer simulations of stellar dynamics, large-scale structure formation, and galaxy formation would have been possible. However, this technique still suffers from fundamental limitations in that almost all contemporary implementations violate time symmetry at the basic level, by setting the step size equal to its desired value at the beginning of each step (the forward method). Such violations are well known to affect the long-term stability of the simulations, resulting in artificial dissipative behaviour, which is often (but not necessarily) revealed by a secular drift of the total energy.




 §.§ Can exact reversibility be achieved?


Adapting the size h of each time step reversibly is non-trivial because the 
time-step function T(ξ), which provides the ideal value for h given the state ξ of the trajectory, can only be evaluated at the start of each step to inform the choice of its size h. Thus, T and h are not naturally synchronised and to achieve time reversibility some form of synchronisation is necessary.

This problem of synchronising h and T also exists in the case of a single continuous step size (unconstrained by the block-step scheme). This simpler situation has been well studied and several solutions have been proposed and demonstrated to give good results <cit.>. These methods essentially integrate the step size itself in a reversible way.

However, porting these methods to the block-step scheme and at the same time retaining reversibility appears impossible (at least when efficiency is retained. So how does the block step render the problem so much more difficult?







  §.§.§ Discreteness and the horizon of predictability


The effect of the block-step scheme is to discretise the step sizes h. Without such discretisation synchronisation between T and h can be achieved by exactly solving an explicit reversible adaption condition, such as equations (<ref>) or (<ref>). When discretising h, this approach is no longer viable, as none of the discrete values will solve the reversibility condition exactly.

Instead, one has to follow an explicit and reversible adaption condition with an auxiliary continuous variable τ, which in turn is then used to inform the choice for the discretised h. Complications arise because the time evolution of τ depends on the actual discrete step sizes h used, which creates a mutual dependency between h and τ and hence an implicit relation for h. Because only a finite number of discrete values are allowed for h, this implicit relation can be solved without iterations, using the longest allowed step in case of ambiguity.

The reason for preferring the longest allowed step is that the merit of the alternative, a pair of two shorter steps, cannot be assessed without violating time symmetry, since the second short step is still in the future and beyond the horizon of predictability.

The integration of the continuous step-size variable τ is inevitably subject to errors and, as a consequence, different values are obtained by different sequences of steps. Consider the situation depicted in Fig <ref>. If two short steps are taken between t_0 and t_2. then the value τ_for for τ(t_1) predicted at time t_0 in the forward direction differs (in general) from τ_back, the value predicted for the same quantity at time t_2 in the backward direction.

If this difference is that between taking the long step from t_0 to t_2 or two shorter steps, an irreversible step-size change results. Thus, if two short steps are taken instead of one long one, then the value τ_back for the alternative longer step is beyond the horizon of predictability at the moment τ_for is computed. This mechanism only produces irreversibly step-size changes with h_for < h_back but of either sign of Ṫ. This is exactly what I find for the schemes (<ref>) and (<ref>): all irreversible step-size changes are of this kind and their frequency is substantially larger when removing the block-step's hierarchical synchronisation (but still requiring discretised step sizes). 







  §.§.§ The case of zonal time-step functions


A zonal time-step function T=T(x⃗) depends only on the particle position. This is useful in simulations of near-equilibrium galactic dynamics, when the local orbital time is well described by a simple function of radius. The discrete step sizes allowed by the block-step then correspond to radial zones <cit.>.

For zonal time-step functions (and using the kick-drift-kick leapfrog) it seems that exact reversibility can be obtained with the block-step, because the value T_n+1 can be readily computed from the information known at time t_n.
, for example, uses for h_n+1/2 the minimum of the discrete step sizes required at t_n and t_n+1. However, the trajectory is subject to integration errors and hence so is the step size, such that the problem discussed in the previous section still pertains, as explained in Fig. <ref>, and even for zonal time-step functions exact reversibility is elusive.




 §.§ Applicability to contemporary N-body methods

The importance of the adverse effects of irreversible step-size adaptation varies considerably between different applications of the N-body method and so does the necessity to and potential benefit of applying any of the methods considered in this study. One potentially serious problem is the noise level, as quantified by |Ṫ|, of any practical time-step-function implementation, which for most of the adaptation schemes will result in poorer reversibility and hence long-term stability than for single-orbit integrations.




  §.§.§ Simulations of collision-less dynamics

In contemporary N-body simulations of collision-less stellar dynamics (galaxy formation and interactions, as well as large-scale structure formation) the adverse effects of irreversible step-size adaptation from the naive forward scheme are presumably tolerable, i.e. the resulting errors remain small over the duration of the simulations (though this has not been rigorously validated), for two reasons. First, such simulations only cover ∼10^2-3 dynamical times, many fewer than the test orbit integrations presented here.

Second, the time-step function most commonly used in such simulations, T=η√(ϵ/|a⃗|), where ϵ is the softening length and a⃗ the acceleration, varies only weakly over typical orbits. (In fact, this time-step function is sub-optimal, since for the vast majority of orbits it gives unnecessarily short time steps, see  for a better alternative.) Moreover, in situations where ϵ is smaller than the inter-particle separation close encounters will result in |Ṫ|≪̸1. For most of the stepping schemes discussed, this inevitably leads to a degradation of the performance compared to the simple test-orbit integrations considered in this study. Thus, if a scheme with considerably better reversibility than the state-of-the-art forward method is required, some form of adaptive force softening appears desirable.




  §.§.§ Simulations of collisional dynamics

Simulations of star clusters, on the other hand, span usually ∼10^4-6 dynamical times, when adverse effects from violating time symmetry may well become important. However, such simulations are usually done very carefully, typically by controlling or monitoring the total energy error and hence guarding against artificial secular dissipation. 

In this case, the total net energy error may still be dominated by the adverse effects of the step-size adaptation. If this is so, then a more careful scheme for adapting the step sizes may improve the efficiency of collisional N-body simulations. The problem of fluctuating time-step functions is likely much less of a problem than for collision-less dynamics, simply because shorter time steps are chosen for enhanced accuracy. This directly reduces |Ṫ| and with it the rate of irreversible step changes for all of the adaptation schemes considered in this study and in contrast to the forward method.








  §.§.§ Simulations of planetary systems

Finally, N-body simulations of planetary systems span ∼10^8-12 dynamical times, enough for even mild secular dissipation to accumulate. In order to avoid this, such integrations are typically done with a single global time-step size (which may be adapted).

Of all the block-step based adaptation schemes considered in this study, only those of Section <ref> appear at all suitable for integrations of planetary systems. These techniques avoid oscillations of the auxiliary continuous step-size variable τ and only suffer from the apparently unavoidable irreversibilities of the type described in Figs. <ref>&<ref> and Section <ref>.

In order to assess how the rate of these irreversible step-size changes depends on the intended integration accuracy, i.e. on the parameter η of the time-step function (<ref>), I ran some integrations of the e=0.9 Kepler orbit for 10^6 orbits (100 times longer than in Fig. <ref>), see Fig. <ref>. Encouragingly, the net-rate R_irr of irreversible step-size changes (but also their total number) decrease like η^2 (the same scaling as for the energy error). This is because decreasing η increases the integration accuracy not only of the trajectory, but also of the continuous step-size variable τ, and hence reduces the chances of irreversible step-size changes.

As a consequence, the ratio of the irreversibility-driven long-term energy error over the short-term energy error decreases like ∼η^3 with decreasing η, independent of the order of the integrator (which determines the short-term energy error). As a consequence, the effects of irreversibilities on the long-term energy error are hardly relevant for sufficiently small η. Indeed, there is no visible correlation between the instances of irreversible step-size changes in the bottom panel of Fig. <ref> (indicated by thin red lines on top and bottom) and the run of the long-term energy error. This suggests that the methods of Section <ref> may well be useful in long-term integrations of planetary system dynamics.




§ CONCLUSION


My attempts to improve on the simple forward method for adapting individual particle step sizes were met with varied success. When merely trying to reduce deviations form reversibility (section <ref>), progress is possible either by extrapolating the step size or by a try-and-reject approach, which however comes at increased computational costs. Porting these methods to N-body simulations requires a well-behaved time-step function (without short-term fluctuations), which may require adaptive force softening.

The situation is more complicated when attempting explicit reversibility. Adapting the scheme of <cit.> to the block-scheme step (section <ref>) obtains methods that suffer from oscillations of the step size or of an auxiliary continuous step-size variable. These oscillations are already present in the original method of , but are much more problematic with the block-step scheme, where they cause irreversible step-size changes. More promising is the idea, pursued in section <ref>, to integrate a continuous step-size variable, analogous to the method of <cit.> for non-discrete h. This avoids the oscillations and obtains near-reversibity, but requires the computation of the time derivative of the time-step function. Such a method may well be useful in long-term integration of planetary systems, when reversibility is much more important than in any other N-body model.

It appears that no practical (explicit or nearly explicit) scheme exists to adapt individual block-step discretised particle step sizes exactly reversibly. As outlined in sub-section <ref>, the ultimate reason appears to be the discretisation of the step size itself.




§ ACKNOWLEDGEMENTS

I thank David M. Hernandez for a critical reading and Scott Tremaine for useful conversations on symplecticity and time symmetry. WD is partly supported by STFC grant ST/N000757/1.

mnras





§ A FLIPPING SCHEME




The obvious way to adapt the adaptation method (<ref>) to discrete step sizes is (with μ the geometric mean) to set

    h_n+1/2 = h_block(T_n^2/h_n-1/2),

which gives the scheme

    δ r = {[                     -1                     if                h_n-1/2                      ≤              1√(2) T_n and block-step allows,;                      0                     if                h_n-1/2                      ≤                   T_n,;                     +1           4lotherwise. ].

The only difference to the forward scheme (<ref>) is the factor 1/√(2) instead of 1/2 in the first clause.




Fig. <ref> shows the ratio of long- to short-term energy error for this scheme. Obviously, this is only a little better than the forward scheme (<ref>). So obviously constraining the time steps to follow the block-step breaks the reversibility of this method, but how?

Fig. <ref> plots the run of radius, h, T, and energy over the first orbit in the cusp model. One sees first that the step sizes flip between two values if 2^1/2-r≤(T/h_max)≤2^1-r for some rung r. By itself this flipping is reversible. However, there are ten such flips between steps 59 and 88 on the way to peri-centre, but only nine between steps 158 and 185 on the way out again – a violation of time symmetry. Upon closer inspection, one finds that the last flip on the in-falling part of the orbit (at step 88) is not reversible: in reversed time the scheme (<ref>) would not chose to change h=2^-7 to h=2^-6 at that moment.

This is a consequence of the restriction to changes in h by at most a factor two. If allowing larger rung changes, the algorithm jumps from h=2^-6 to h=2^-8 at step 88 and in reversed time jumps back. Unfortunately, enabling |δ r|>1 allows the occasional much too long time step with disastrous consequences for the integration accuracy. Using different limits for δ r and/or functions μ(x,y) can give some improvement over the reported attempt, but did not yield a reliable method significantly better than the simple forward scheme (<ref>). Moreover, the flipping of step sizes is inefficient for N-body force solvers, since in some places half the particles will have rung r and the other r+1, switching every other step.

The scheme (<ref>) only produces irreversible step-size adaptations of the first two types in equation (<ref>), which give energy errors of the same sign but opposite to the second two types. 


