The main idea of the AC scheme is to reduce the effort of evaluating the force contribution from distant particles by combining two polynomials based on separate time-scales. Splitting the total force on a given particle into an irregular and a regular component by

we can replace the full N summation in (1.1) by a sum over the n nearest particles, together with a prediction of the distant contribution. Likewise, (3.1) may be differentiated once for the purpose of low-order prediction. This procedure can lead to a significant gain in efficiency, provided the respective time-scales are well separated and n ^ N. The two-time-scale scheme for one particle is illustrated in Fig. 3.1.

In order to predict the coordinates of neighbours at the irregular timesteps, we need to extrapolate the regular force and its first derivative. The former can be readily obtained from (2.3), without the term DR and using appropriate time arguments, t — Tk, where Tk refers to the regular force times. Similarly, differentiation of the regular part of (3.1) using Tk yields the corresponding derivative which should be added to F^.

To implement the AC scheme, we form a list for each particle containing all members inside a sphere of radius Rs. In addition, we include any particles within a surrounding shell of radius 21/3Rs satisfying the condition R • V < 0.1R'^/ATi, where again V = vi — v and ATi denotes the regular time-step. This ensures that fast approaching particles are selected from the buffer zone. A subsidiary search criterion based on m/r2 may be introduced for including more distant massive bodies.

The size of the neighbour sphere is modified at the end of each regular time-step when a total force summation is carried out. A selection criterion based on the local number density contrast has proved itself for a variety of problems, but may need modification for interacting subsystems. To sufficient approximation, the local density contrast is given by c=N (Rs )3 • (3.2)

where ni is the current membership and rh is the half-mass radius. In order to limit the range, we adopt a predicted membership np

subject to np being within [0.2nmax, 0.9nmax], with nmax denoting the maximum permitted value.* The new neighbour sphere radius is then adjusted using the corresponding volume ratio, which gives

An alternative and simpler strategy is to stabilize all neighbour memberships on the same constant value np = n0 [Ahmad & Cohen, 1973; Makino & Hut, 1988]. Since it is not clear whether this strategy is advantageous, the original density-related expression which has proved itself is still being used. A recent analysis by Spurzem, Baumgardt & Ibold  concluded that the Makino-Hut relation np a:

N3/4

may not be necessary for large systems with realistic density contrasts. Hence it appears that relatively small neighbour numbers of about 50 are sufficient to take advantage of the AC scheme even for N ~ 40 000. In any case, the actual code performance does not depend sensitively on the neighbour number because of compensating factors: reduction of ni speeds up the irregular force summation but leads to shorter regular time-steps. A number of refinements are also included as follows:

• To avoid a resonance oscillation in Rs, we use (np/ni)1/6 in (3.4) if the predicted membership lies between the old and new values

* As a recent refinement, np = np(rh/r) is now used outside the half-mass radius.

• Rs is modified by the radial velocity factor 1 + r, • v^ATi/r2 outside the core

• The volume factor is only allowed to change by 25%, subject to a time-step dependent cutoff if AT, < 0.01tcr

• If n < 3, r < 3rh and the neighbours are leaving the standard sphere, Rs is increased by 10%

Following the recalculation of the regular force, the gain or loss of particles is recorded when comparing the old and new neighbour list. Regular force differences are first evaluated, assuming there has been no change of neighbours. This gives rise to the provisional new regular force difference

where FRld denotes the old regular force, evaluated at time T0, and the net change of irregular force is contained in the middle brackets. In the subsequent discussions of the AC method, regular times and time-steps are distinguished using upper case characters. All current force components are obtained using the predicted coordinates (which must be saved), rather than the corrected values based on the irregular term Dj since otherwise (3.5) would contain a spurious force difference. The higher differences are formed in the standard way, whereupon the regular force corrector is applied (if desired).

A complication arises because any change in neighbours requires appropriate corrections of both the force polynomials, using the principle of successive differentiation of (3.1). The three respective Taylor series derivatives (2.4) and (2.5) are accumulated to yield the net change. Each force polynomial is modified by first adding or subtracting the correction terms to the corresponding Taylor series derivatives (2.3) (without the D4 term), followed by a conversion to standard differences by (2.7).

Implementation of the AC scheme requires the following additional set of regular variables: Fr, DR, DR, DR, AT, To, Tu T2, T3, as well as the neighbour sphere radius Rs and neighbour list of size nmax + 1. The corresponding code for softened potentials (nbody2) has been described in considerable detail elsewhere [Aarseth, 2001b]. In the point-mass case, close encounters are treated by KS regularization, giving rise to an analogous code for the AC scheme [Aarseth, 1985a], as well as a new Hermite AC code with full regularization. As an application of the AC method, we mention a scheme for studying the interaction between a star cluster and a molecular cloud in which the latter is treated by smoothed particle hydrodynamics [Theuns, 1992a].

The combination of the Hermite method with the AC scheme (HACS) gives rise to a very efficient treatment. Following the discussion in chapter

2, the implementation of these two methods is much facilitated. Here we summarize some of the main features for evaluating the derivatives. A full description with actual comparisons is given elsewhere for a softened interaction potential [Makino & Aarseth, 1992].

In HACS, we need to form the current total force for the predictions as well as for the next corrector step by the low-order extrapolation

Similarly, the total force derivative is just the sum

At the end of a regular time-step, we need the equivalent of the difference (3.5), namely the change in the regular force assuming there has been no loss or gain of neighbours, which yields

Note that the inner bracket of (2.20) and (2.21) requires the reverse sign of this expression. Likewise, the corresponding new regular force derivative based on the old neighbours is given by fR1) = FR1) - (f(1) - Fj1)), (3.9)

where the quantities F(1) denote the respective new derivatives and the terms in brackets contain the net change. Finally, we combine (3.9) with the previous regular force derivative to form the higher derivatives used for the corrector. If required, the corresponding corrected values at the end of the interval may be obtained by including contributions from the change of neighbours.

In the Hermite versions, we employ the code convention that the quantities Dk denote Taylor series derivatives, rather than differences as used in the old polynomial formulation. Consequently, the first derivatives are initialized by