{'fraction_non_alphanumeric': 0.03472893350375659, 'fraction_numerical': 0.013670231006640271, 'mean_word_length': 4.653917334978408, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gauss's method for secular dynamics, softened

We show that the algorithm proposed by Gauss to compute the secular evolution of gravitationally interacting Keplerian rings extends naturally to softened gravitational interactions. The resulting tool is ideal for the study of the secular dynamical evolution of nearly Keplerian systems such as stellar clusters surrounding black holes in galactic nuclei, cometary clouds, or planetesimal discs. We illustrate its accuracy, efficiency and versatility on a variety of configurations. In particular, we examine a secularly unstable unstable system of counter-rotating disks, and follow the unfolding and saturation of the instability into a global, uniformly precessing, lopsided (m=1) mode.


INTRODUCTION
The theory of secular variations of the planetary orbits was introduced by Lagrange over two centuries ago. This approach to celestial mechanics is based on the observation that the orientation of eccentric orbits in a Keplerian potential is fixed: in a coordinate representation, the argument of periapsis ω and the longitude of the node Ω are stationary. Because of these conserved quantities, the time-averaged torque between two planets on fixed orbits is generally non-zero if the orbits are both eccentric and/or mutually inclined. Thus, if the typical planetary mass is m and the mass of the host star is M⋆, the interplanetary gravitational forces lead to both short-term oscillations, having periods of order the orbital period P and fractional amplitudes of order m/M⋆, and long-term or secular oscillations, having periods of order (M⋆/m)P , fractional amplitudes in eccentricity and inclinations of order unity, and constant semi-major axis. The large amplitude of the secular oscillations means that in most cases they dominate the evolution of the planetary system. In the solar system, secular variations are usually computed as power series in the eccentricities e, inclinations i, and mass ratios m/M⋆ of the planets, which are all small. Current secular theories employ differential equations of motion containing ∼ 10 5 terms up to order (m/M⋆) 2 , e 5 , and i 5 (Laskar 1986).
Expansions of the equations of motion in powers of eccentricity and inclination work well for the orbits of planets in the solar system (median eccentricity and inclination 0.05 and 1.8 • ), but not so well for orbits with larger eccentricities and inclinations, such as those of asteroids (median eccentricity and inclination 0.14 and 7 • ), extrasolar planets (median eccentricity 0.2), the stars that surround the central black holes found in many galaxies, and planets or binary stars executing Kozai oscillations excited by a companion star (Innanen et al. 1997;Mazeh et al. 1997;Holman et al. 1997;Wu & Murray 2003;Fabrycky & Tremaine 2007).
The primary goal of this paper is to develop numerical methods for secular dynamics that are valid for large eccentricities and inclinations. In so doing we shall restrict ourselves to equations of motion that are first order in the ratio m/M⋆. This approximation does not work well in the solar system, mostly because of near-resonances between the outer planets (e.g., the 'great inequality' between Jupiter and Saturn; see also Murray & Holman 1999), but should be more useful in systems where the mass ratio m/M⋆ is smaller and/or we are seeking a qualitative understanding of the secular evolution rather than an accurate ephemeris.
Although the results we describe in this paper are applicable to any near-Keplerian system, for conciseness we shall refer to the central body as a 'star' and the satellites in orbit around it as 'planets', except in §4 where we examine the performance of the algorithm with parameters appropriate for stars orbiting a black hole. To first order in the planetary masses, the equations of motion for secular dynamics are derived from a governing Hamiltonian that represents the sum of the time-averaged interaction potential energies between all pairs of planets in the system. Therefore the central ingredient in any numerical method for secular dynamics is the efficient evaluation of the average interaction energy between two fixed Kepler orbits, each specified by its semi-major axis a and shape parameters e (eccentricity), ω (argument of periapsis), I (inclination), and Ω (longitude of node). The averaging can be thought of as smearing each planet into a 'ring' in which the mass per unit length is inversely proportional to the speed of the planet in its orbit. We call a code that follows the secular gravitational dynamics of N planets an 'N-ring code' by analogy with N-body codes that follow the full gravitational dynamics. The simplest way to calculate the interaction energy between rings labelled by α and β is to choose K points equally spaced in mean anomaly on each ring and then compute − K i,j=1 Gmαm β /(K 2 ∆ij ), where mα and m β are the masses of the two bodies and ∆ij is the distance between point i on ring α and point j on ring β (e.g., Gürkan & Hopman 2007). This approach is rather inefficient, both because it requires O(K 2 ) computations per ring pair and because it converges slowly, especially if the orbits have a close approach (min ∆ij ≪ max ∆ij ). For example, Gürkan & Hopman (2007) typically find that K ∼ 10 3.5 is needed for 1 per cent relative accuracy in the potential.
A better approach, which we pursue in this paper, is to compute the potential from ring α at an arbitrary point r by integrating analytically over the mass elements of ring α, and then to numerically integrate this potential over the locus of points r that lie on ring β. The rather involved analysis needed to do this was accomplished by Gauss and Hill in the nineteenth century 1 . We shall call this approach 'Gauss's method', and the use of Gauss's method in an N-ring code the 'Gaussian ring algorithm'.
When he wrote his classic memoir on the force exerted by an elliptic ring on a mass point, Gauss (1818) achieved at least three lasting outcomes: he provided dynamics with a vivid illustration of the averaging principle; he introduced the arithmeticgeometric mean as an efficient method for computing Legendre's elliptic functions; and he showed how the evaluation of the first-order average of the planetary equations can be dramatically accelerated with the help of such functions. Hagihara (1971) gives a thorough review of Gauss's method, its algorithmic manifestations and its applications. Its main attraction of course is the efficient averaging of the equations of motion over the fast, orbital, time scale. The resulting differential equations, in the case of the solar system, can be evolved with a timestep of at least 500 yr-over 10 4 times larger than the timestep needed to follow the unaveraged equations of motion. The main drawback of Gauss's method, as already noted, is that it neglects perturbations that are second-order or higher in the planetary masses, some of which are near-resonant terms that, in the solar system, can be almost as strong as the secular terms that are first order in the masses.
An implementation of Gauss's method by one of us (J.T.) was used to explore Kozai (1962) oscillations of the then newly discovered planet 16 Cygni Bb (Holman et al. 1997). In most applications, Kozai oscillations arise from the influence of a distant companion body. In the case of 16 Cygni Bb, this is a stellar companion to the host star with a projected separation of ∼ 850 AU, 500 times the planetary semi-major axis. At such large separations the perturbing potential from the stellar companion is dominated by its quadrupole moment, and in this case the averaged interaction potential between the companion and the planet can be computed analytically. Nevertheless, Gauss's method offers important insight into this problem. The reason is that in the quadrupole approximation the averaged interaction potential is axisymmetric in the reference frame whose equator is the orbital plane of the companion, even if the companion orbit is eccentric. Thus both the energy and the component of angular momentum along the axis of the companion orbit are integrals of motion, and the secular dynamics is described by a Hamiltonian with only one degree of freedom. The existence of this second integral is an accidental consequence of the properties of the quadrupole potential, and without the quadrupole approximation the integral is not present. Gauss's method allows us to explore the two degree-of-freedom problem defined by the full secular dynamics, which for example shows that in many cases initially circular planetary orbits are chaotic and allows us to estimate the width of the chaotic zone 2 .
What drives much of the current interest in Gauss's method is the wish to study the evolution of a distribution of stars within the sphere of influence of the central black hole in a galaxy. In this situation the central object enforces Keplerian dynamics, perturbed by the mutual gravitational attraction between the stars. Here then is a natural context for Gauss's method. However, in such problems, one immediately faces close encounters between particles, which translate into ring crossings. At a crossing, the force between rings changes discontinuously 3 , a jump that is difficult for numerical integrators to follow. Such was the experience of Rauch & Tremaine (1996), who used an unsoftened version of Gauss's method to simulate resonant relaxation around black holes and found that the speed-up resulting from averaging is virtually cancelled by the slow-down caused by ring crossings.
In N-body simulations, one often resorts to softened gravitational interactions to alleviate the cost of resolving nearcollisions between particles (and to better approximate collisionless behavior). It is natural to ask whether Keplerian rings that perturb each other with softened gravitational forces can be evolved with the help of a softened variant of Gauss's method. We were pleased to learn that the answer is yes, and we will proceed to tell you how. We start in §2 by reviewing the underpinnings of Gauss's algorithm, and how they are modified for the softened case. We thereby derive a formulation of Gauss's method for softened gravitational interactions. We discuss how to implement the method numerically in §3. In §4 we describe numerical examples that illustrate its accuracy, its strengths, and its limitations (limitations of secular averaging in general), by direct comparison with simulations of the analogous unaveraged system. These tests lead us to §5, where we present simulations of a (linearly) unstable configuration of two counter-rotating annuli, which track the instability to its (nonlinear) saturation as a global, uniformly precessing, lopsided (m = 1) mode. We conclude in §6 with a discussion of the main features, potential applications, and future improvements of our algorithm.

GAUSS'S METHOD REVISITED
We argued in the preceding section that to compute the interaction potential between two Keplerian rings using K equally spaced points per ring required O(K 2 ) calculations. Gauss (1818) showed that the work can be reduced to O(K) calculations by integrating analytically over one ring with the help of Legendre's elliptic functions. Hill (1882) elaborated the arguments in Gauss's paper, and described an algorithm for the exact calculation of secular evolution which is accurate to first order in the masses and all orders in eccentricity and inclination. Halphen (1888), and Goriachev (1937) after him, followed an alternative, geometric, route to elliptic functions, which was later described in tensorial language by Musen (1963). While Gauss's exposition is purely algebraic in flavor, Halphen exploits the geometry of the cone (with base on the disturbing ring and apex at the point where the force is evaluated), to motivate a transformation to the cone's principal axes that makes Gauss's insight transparent. While the algebraic manipulations work through effortlessly for softened interactions, it is less obvious how to soften things geometrically, since the softening scale breaks a crucial symmetry. We choose to cover the algebraic derivation in the body of the text, relegating the geometrical treatment to Appendix A.
The algorithm, softened or not, has two independent components: on one hand, the expression of the averaged equations of motion in a flexible, singularity-free coordinate system, while presuming a means, efficient or not, for computing averages of interest (the subject of this section); on the other, an efficient method for computing these averages (the subject of §3). Taken together, these maneuvers produce a regularized vector field which one can evolve with giant strides using an appropriate numerical integrator.
The calculations below largely follow Hill (1882); they are repeated here partly to include the generalization to softening, partly to update Hill's derivation to modern notation and style, and partly to express the results in a vectorial, coordinate-free notation that is much more efficient for numerical integration algorithms.

Regularized equations
We consider a particle (called the 'planet', though it could be a star, moon, planetesimal, comet, etc.) bound to a massive central body (called the 'star', though it could be a black hole, planet, etc.). The planet follows a Keplerian orbit and is subject to relatively weak external forces, which cause time variations in the (osculating) Keplerian orbital elements. We express the equations governing these variations in terms of non-singular vectorial variables, then recover their secular average for the case of interest, namely when particles perturb each other's Keplerian orbits with their mutual gravitational attraction (Newtonian or softened).
We denote by M⋆ the mass of the star, and adopt the following notation for properties pertaining to the planet: • m mass of the planet, • r radius vector from the central mass, with r = |r|, • v ≡ṙ velocity vector, • a the semi-major axis, • e the eccentricity, • n = G(M⋆ + m)/a 3 the mean motion, with P = 2π/n the orbital period, • l, E, and f the mean, eccentric and true anomalies, respectively.
When considering perturbations between planets, unprimed variables refer to the perturbing planet, and the primed counterparts to the perturbed planet.
In the planetary context, it is common to describe the perturbations to a Keplerian orbit using formulae that describe the variation of the orbital elements eccentricity, inclination, etc. (e.g., Murray & Dermott 1999). These equations suffer from singularities at zero eccentricity and inclination, as well as at unit eccentricity. The singularities at zero eccentricity and inclination are traditionally handled by a canonical transformation to appropriate Cartesian coordinates for the orbital elements (the familiar trick for dealing with the singularity of cylindrical polar coordinates at the origin). Such transformations work reasonably well in 'cold' settings where the eccentricity and inclination are small. Near unit eccentricity (nearly radial orbits) the singularities can also be removed by an appropriate canonical transformation (Tremaine 2001). However, in hot environments such as the stellar cusp around a black hole the eccentricities and inclinations take on all values so one must regularly switch between coordinate charts if one insists on working with elements such as eccentricity and inclination or even their Cartesian analogs. We initially implemented such switches, and found them costly and cumbersome. Instead, we chose to parametrize an orbit with the coordinate-free (though redundant) angular momentum and eccentricity (or Laplace, or Runge-Lenz) vectors (see Hagihara 1971or Murray & Dermott 1999. We associate with an orbit three orthogonal unit vectors:ẑ, along the rescaled angular momentum vector L = √ 1 − e 2ẑ (normal to the orbital plane, scaled by mna 2 ),x, pointing toward periapsis, andŷ =ẑ×x, perpendicular to both. The eccentricity vector is The position vector of the planet is r = xx + yŷ where x = r cos f = a(cos E − e), and y = r sin f = a √ 1 − e 2 sin E, with the radius r = a(1 − e cos E). We introducer, a radial unit vector, and the tangential vectort =ẑ×r; explicitlŷ We will also need the velocity vector v =ẋx +ẏŷ, whereẋ = −na 2 sin E/r andẏ = na 2 √ 1 − e 2 cos E/r (usingĖ = na/r). The semi-major axis, angular momentum vector, and eccentricity vector of the disturbed (primed) planet respond to a perturbing acceleration f ′ on that planet as The first equation expresses the change in the Keplerian energy resulting from the perturbing force, the second the change in angular momentum resulting from the associated torque, and with it the change in orientation of the orbital plane; the third governs the eccentricity and periapsis direction of the orbit, hence dictating shape and orientation in the orbital plane. We shall work with the radial (R), tangential (S) and normal components (W ) of the perturbing acceleration: Substitute into equation (3) to get where anḋ These equations govern the secular evolution of an osculating Keplerian orbit in the vectorial representation. They are general enough to accommodate a great variety of forces, conservative (external perturbers, such as satellites, planets, stars, or black holes, on asteroids, planets, or stars) or not (drag, tides, dynamical friction), that may affect the motion of a particle whose motion is otherwise dominated by a massive central body. In this work, we are concerned with analytic averaging over softened mutual gravitational interactions in a near-Keplerian cluster of particles, to which we shall restrict our attention in the following.

Averaged equations
The softened gravitational force on the primed planet, m ′ f ′ = −∇ ′ Φ, derives from the disturbing function and b is the softening length (primed and unprimed refer to perturbed and perturbing respectively). The first term is Plummer's potential (which reduces to Newton's for b = 0) coupling the two particles, the second (the indirect term) is a non-inertial contribution resulting from pinning the reference frame to the central mass. Note that we soften only the interplanetary potential, not the potential between each planet and the central star. We are interested in the secular evolution of the osculating Keplerian orbital elements, which we obtain by averaging the equations of motion (7) over the mean anomalies l and l ′ of the disturbing and disturbed planets. In other words we are after Before proceeding further, we make a couple of simplifying remarks. (i) The indirect term in Φ contributes a force r/r 3 that averages to zero over l. Therefore we may drop the indirect term. (ii) Another classical result of secular dynamics is that da ′ /dt averages to zero over l ′ , since the average is proportional to the work of a conservative force over a closed path. Thus we may treat the semi-major axis a ′ and mean motion n ′ of the disturbed planet as constants of motion when examining conservative perturbations in secular dynamics (in the case of unsoftened gravitational forces this is Poisson's theorem, but the result requires only that the perturbing force is conservative and time-independent and hence holds for softened forces as well).
Below, we shall demonstrate that the average of the direct force over l is analytic and expressible in terms of elliptic functions. For now, we denote by [R] l , [S] l , and [W ] l the l-averaged radial, tangential and normal components of the direct force, and note that the l-averaged equations can be directly translated from their unaveraged counterparts in equations (6) and (7), by replacing R, S and W by their respective l averages. Now consider the average over the orbit of the disturbed planet. For any function X(E ′ ), Inspecting equations (6) and (7) we find that the average over l ′ will at most involve the second-order term of the Fourier expansion of [R] l , [S] l and [W ] l in their periodic dependence on the eccentric anomaly E ′ , leading to: and Here with similar definitions for S n c , W n c , etc. For conservative forces, [da ′ /dt] ll ′ = 0 and this gives us the identity The cost of non-singular coverage of the full range in eccentricity and inclination using the vectors L ′ and A ′ has been to expand the phase space of an N-ring system from 4N (reduced from 6N as a result of averaging over the mean anomaly l ′ and conservation of the semi-major axis a ′ ), to the redundant 6N -dimensional phase space spanned by the angular momentum and eccentricity vector of each ring; N of the redundant components are removed by noting that L ′ and A ′ are orthogonal, while the other N reflect the relation between the magnitudes of the two vectors, (11)-(13) if the perturbing forces are conservative.

Analytic averaging
We shall go over the crux of our exercise, namely the demonstration that Gauss's averaging works out perfectly for softened potentials.
The average of the gravitational potential Φ over the orbit of the perturbing (unprimed) planet can be written where with Note the inequalities We pause briefly to emphasize the seemingly innocuous absorption of the softening b 2 in A b . As far as Gauss's method is concerned, this is the only direct change that is brought about by softening.
The average perturbing acceleration is where At this point, we write introduce a new angular variable T (the perspective anomaly) defined by and demand that where Q is a 3 × 3 matrix with elements Qij , i, j = 0, 1, 2. The geometric motivation for this change of variable is described more fully in Appendix A. Briefly, the orbit average [·] l can be thought of as a line integral around the orbit, with r and E the polar variables in the orbit plane for this line. We show in the Appendix that the orbit average is unchanged if we shift the line integral to a new contour, defined by the intersection of any plane with the cone that contains the orbit and has its apex at the disturbed planet. The perspective anomaly is a polar coordinate in this new plane. The plane is later chosen to make the line integral as simple as possible.
Since sin 2 E + cos 2 E = sin 2 T + cos 2 T = 1 we must have The first of these simply states that x is confined to a cone; the second restricts Q to be a matrix that conserves the quadratic form x T Mx, where x T is the transpose of the column vector x and M is the diagonal matrix diag (1, −1, −1). This requirement can be written Thus Q is a pseudo-orthogonal matrix; in the language of special relativity, Q is a Lorentz transformation in two spatial dimensions, M is the Minkowski tensor, and the vectors x and y lie on the light cone passing through the origin. We define a diagonal matrix with complex coefficients, C ≡ diag (1, i, i). Then where * denotes complex conjugation and I is the unit matrix. Let L ≡ C * QC * . Then using equations (26) and (27) we have Thus L is orthogonal if Q is pseudo-orthogonal, and any pseudo-orthogonal Q can be written in the form where L is orthogonal. Note that even though L may have complex coefficients, it is orthogonal (L −1 ij = Lji) rather than Hermitian (L −1 ij = L * ji ). Because we allow complex coefficients the inner product of a vector with itself, is not necessarily positive or even real. We may now re-write equation (17) as where Gauss now asks for the pseudo-orthogonal matrix Q that diagonalizes the quadratic form for where D is some diagonal matrix. Using equation (29) we re-write the quadratic form as where Since R is symmetric, (i) its eigenvectors e j can be chosen to be orthonormal, e j ·e k = δ jk ; (ii) the matrix L defined by Lij = e j i is orthogonal; (iii) L T RL = diag (λi), where {λi} are the eigenvalues of R. Thus equation (33) can be written as which is the desired diagonalization (32), with Q related to L by (29) and The eigenvalues of R are the solutions of the cubic equation Using the relations (18) it is straightforward to show that Furthermore from (19) a 2 (1 − e 2 ) < A b so y(λ) must have three real roots, one between −C and 0, one between 0 and a 2 (1 − e 2 ), and one between a 2 (1 − e 2 ) and A b . It proves useful to label these roots so that λ0 > λ1 > λ2; thus Relations between these roots include Explicit expressions are where The eigenvectors of R are where α k , which may be complex, is chosen so that e k ·e k = 1 or Using the relations (39) this can be manipulated into a more compact form, where l, m are the indices other than k from {0, 1, 2}. With the ordering (38), We choose the phases of the eigenvectors so that Then The determinant of this matrix is We now examine the relation between the eccentric anomaly E and the perspective anomaly T in more detail. From the (0, 0) component of the equations Q T MQ = M and QMQ T = M we have it is straightforward to show that Q00 > 1 so ∆ is real and positive. Then from equations (22) and (23) ∆x0 From equation (24) we also have Using the (0, 0), (0, 1) and (0, 2) components of the equations Q T MQ = M we can simplify this to so from equations (51) we have Similarly, and Using equations (39) and (47) it is straightforward to show that where s = ±1 is defined by equation (48). Thus Combining this result with (54), This result relates the perspective and eccentric anomalies and shows that the perspective anomaly circulates through 2π (if s = 1) or −2π (if s = −1) when the eccentric anomaly circulates through 2π. Moreover the differentials are related by where the last equality employs (59).
We may now combine equations (20), (22), (23), (30), and (35) into an expression for the orbit-averaged perturbing acceleration at position r: Replacing y1/y0 by sin T and y2/y0 by cos T , and dropping terms that integrate to zero, we have where the F j are defined by equation (21), and The integrals can be expressed in terms of the complete elliptic integrals Thus The averaged radial, tangential and normal components of the force ([R] l , [S] l , [W ] l ) are found by replacing f ′ in equations (4) by its averaged value (67). Thus, one of the averages is done, the force of a softened ring at a point is calculated. To get the fully averaged secular equations (12), a second average over the perturbed ring itself is required, as given by equation (14). The second average is accomplished by numerical quadrature as described in the following section.

REMARKS CONCERNING IMPLEMENTATION
To be sure, we are not yet at the point where we have an optimal N-ring code, partly because the improvements afforded by averaging lessen the urgency to optimize, partly because we prefer to explore novel dynamical phenomena rather than novel algorithms. Obviously, a careful optimization will eventually be necessary, but for now we simply outline the steps involved in the efficient implementation of the Gaussian ring algorithm, emphasizing difficulties encountered, and possible ways around them. We are engaged in optimizing and parallelizing the code, with a view to simulations of stellar clusters surrounding massive black holes.

Numerical averaging
In order to perform the numerical average over the orbit of the perturbed planet, which gives us equations (11)-(13) for the evolution of its angular momentum and eccentricity vectors, we need to compute the first three Fourier coefficients in the eccentric anomaly E ′ of the force vector (R 0 c , R 1 c , . . . , W 2 s ; see eq. 14). Direct summation at points equally spaced in E ′ (i.e., a discrete Fourier transform) will do in most cases, and we have adopted this approach in the calculations presented below.
The required number of points K for the evaluation of these Fourier coefficients depends both on the error one is willing to tolerate and on the eccentricities and separations of the rings. The Euler-Maclaurin summation formula tells us that the accuracy of the evaluation will improve with increasing K faster than any power of 1/K, typically as exp(−αK) for some constant α > 0. Of course, during close approaches between rings this desirable asymptotic behavior may only occur at rather large K, and the constant α may be disappointingly small.
Hill (1882), who was mainly concerned with planetary motion (nearly circular, nearly coplanar, well-separated rings), was happy with K = 12, even 8. More generally, he pointed out that 'if the number of these values be even, the order of the error. . .will be the same as that of the power of the eccentricities or mutual inclination of orbits, whose exponent is one less than the number of these values'. We can re-state Hill's conclusion-which we have confirmed using expansions of the disturbing functions in eccentricity and inclination-more precisely as follows: Consider two rings with eccentricities and semi-major axes ej, aj, j = 1, 2, and mutual inclination I. Let δ = |a2 − a1|/a where a = 1 2 (a1 + a2), and let e = max(e1, e2). Then if e ≪ δ, I ≪ δ, the fractional error in the time derivatives of the eccentricity and angular-momentum vectors is O[max(e/δ, I/δ) K−1 ] as e, I → 0. 4 This result also implies that the implementation of Gauss's method using numerical quadrature with K equally spaced points is formally equivalent to a secular perturbation theory in which the Hamiltonian is expanded to first order in the planetary masses and K th order in the eccentricities and inclinations (e.g., Libert & Henrard 2008, for K = 12).
Going by this estimate, Hill's use of K = 12 would yield a relative error of 10 −10 in the time derivatives of the eccentricity and angular-momentum vectors of a two-planet system consisting of Jupiter and Mercury (eccentricity 0.2); numerical experiments on this system, with varying orientations, yielded relative errors in the range 10 −9 -10 −12 in reasonable agreement with this estimate. To achieve a similar accuracy for the system consisting of Jupiter and Halley's comet our experiments indicate that one needs K ∼ 200-400, and this again depending on argument of node and perihelion, while keeping the other elements fixed. A softening length b = 0.01a or b = 0.1a can cut this number by about a factor of two or four, respectively. The stellar systems which we will be dealing with have higher eccentricities and inclinations than anything Gauss or Hill envisioned, as well as frequent ring crossings. One can get away with K = 16-32 for nearly circular rings, while achieving a relative error of 10 −10 or smaller in the evolution equations. For nearly radial orbits, we found that we needed up to K = 512 for the same level of accuracy.
It is useful to have a criterion for the accuracy of the numerical average, which can be used to provide an adaptive choice for K for each ring pair. Hill's criterion works when e ≪ δ, I ≪ δ, but not for near-radial or crossing rings. We devised an alternative criterion which is equivalent to Hill's in the limit where Hill's works. We rely on the observation that, had the averages been exact, the equation governing the evolution of a ring's semi-major axis would be [da ′ /dt] ll ′ = 0 (see discussion following eq. 9), which in turn implies the identity (15). Therefore we may simply choose K adaptively so that the numerical evaluation of the absolute value of the left side of (15), divided by n ′2 a ′ to make a dimensionless quantity, is less than some specified tolerance ǫ quad ('quad' for quadrature) for each ring pair. We start with K = 16 and proceed through doublings. Note that Richardsonian extrapolation is not useful here, because the asymptotic convergence is faster than any power of 1/K. What is an appropriate tolerance? A tolerance ǫ quad translates into a tolerance on the relative error in semi-major axes per secular or precession period that is no worse than (N − 1)ǫ quad . This in turn translates into into a relative error over M secular periods no worse than N M ǫ quad . Presuming one is interested in simulating a system of 10 5 particles over 100 precession periods, and one can live with cumulative relative error in the dominant Keplerian energy of 10 −4 , then ǫ quad ∼ 10 −11 should be amply sufficient for such a bound. For most ring-ring interactions, such a tolerance is achieved with 64 to 128 sectors per ring at a softening of 0.01 times the typical size of a ring. There are of course likely to be some ring pairs for which sampling of 16 to 32 times these values may be necessary, and in these cases adaptive sampling proves essential. So far we have described how to carry out numerical averaging over the perturbed ring using equally spaced points in eccentric anomaly. However, we are aware that more sophisticated quadrature routines may be better, particularly at ring crossings when the softening length b is small. We have conducted limited experiments with the QAGS adaptive integrator from GSL (intended to handle integrable singularities by adaptive application of the Gauss-Kronrod [Gauss again] rule). In the rather extreme case of two coplanar, intersecting rings, with a softening of only 0.001 times the orbital radius, QAGS required 1500 judiciously placed integration points to achieve a relative error of 10 −11 ; ten times more points were required by uniform sampling to achieve the same tolerance. The factor of 10 drops to 2 for a softening length which is ten times larger.

Numerical integration
The rings are evolved in time using a Bulirsch-Stoer (BS) integrator. The BS integrator is not tailored to this problem, but for a reliable, brute force exploration, it is efficient and simple to implement. The ideal of course would be to exploit the Lie-Poisson structure of the equations to construct a geometric integrator with an adaptive timestep that would handle the broad range of dynamical time-scales in typical stellar and planetary systems. This is certainly a noble undertaking, which we put off to a later stage, and we will be more than happy if others beat us to it. In any case, in our BS integrations, we set a tolerance ǫint for the error in the integrated vector which is similar to the tolerance ǫ quad that we demand in the numerical quadratures. In most cases, we are quite happy with a tolerance of ǫint = 10 −8 (at which we live with step sizes on the order of 0.1 secular periods); we do push things further in the simulations below, and evolve some configurations with a tolerance of ǫint = 10 −12 , and consequently smaller timesteps on the order of 0.01 secular periods.
Highly eccentric orbits are likely to be common in many nearly Keplerian systems (e.g., stellar clusters around black holes, cometary clouds). In the plane one has the lens orbits of Sridhar & Touma (1999); related families exist in three dimensions and were studied in oblate and prolate potentials by Sambhus & Sridhar (2000). The averaging process we describe here is well-behaved for nearly radial orbits (e → 1), as are the equations of motion for the eccentricity and angular-momentum vectors-which is one reason why we chose to work with vectors rather than orbital elements, which usually give singular equations of motion as e → 1 (Tremaine 2001). Nevertheless we have found that our BS integrator struggles with nearly radial configurations, slowing the integration by a factor of 10 to 1000 in some extreme cases.
Part of the slow down is dynamical, and will be faced by any adaptive integration scheme. The other, and more critical, part is a consequence of the violation of the constraint L 2 + A 2 = 1 in the course of the numerical evolution; it manifests as a square-root singularity in the equations of motion (11)-(13), and forces a slowing down in the integration, as steps are reduced to avoid that singularity. One can imagine a number of solutions to this problem. As suggested above, the best would be to design a geometric integrator which is adapted to the Lie-Poisson structure of the equations at hand, and hence preserves the integral in question automatically. Less elegant, but still effective, would be to enforce the constraint L 2 + A 2 = 1 by hand, renormalizing one or the other vector at regular intervals to eliminate the small errors introduced by the integration scheme. Without modifying the BS integrator, one could impose a reflecting boundary condition at high eccentricity. There are also physical processes that mitigate the problem in many systems: these include general-relativistic precession (Appendix B) and tidal disruption, collisions, or swallowing of stars by the central object. A careful examination of the proper implementation of these and related alternatives would take us too far afield. In the ring simulations described below, which include generalrelativistic precession, we encountered an occasional near-radial ring that slowed down the integration of the whole cluster. We found that this problem could be removed either by brute-force removal of highly eccentric orbits (e > 0.999) or by exploiting the redundancy in the present vectorial formulation. In the second solution, we simply used the eccentricity vector A to solve for the eccentricity at high eccentricities and the angular-momentum vector L to solve for the eccentricity otherwise. The resulting algorithm still slowed down at high eccentricity passages, but only by that dynamical factor which reflects the higher precession rates at these eccentricities, boiling down to a factor of 10 decrease in performance during these phases. We found this intermittent slowing down-about a factor of ten when a ring (or a population of rings) developed mean eccentricity greater than 0.8-tolerable in the cases we studied, and believe it manageable with larger systems, when the code is eventually parallelized.

Evaluation of elliptic integrals
To compute the elliptic integrals one could either iterate using the arithmetic-geometric mean (Gauss once again), interpolate in a table, or rely on Chebyshev polynomial expansions. We tried all three methods and found that smart searches through tables were twice as costly as Chebyshev expansions, which were marginally better than the arithmetic-geometric mean (for an absolute error goal of 10 −14 ). In our current implementation, elliptic integrals are calculated with the arithmetic-geometric mean.

Additional physics
Nearly Keplerian stellar systems are often subject to additional effects with substantial secular contributions. We can easily accommodate the average contributions of tides (Galactic for cometary clouds, bulge/halo for stellar clusters around a central black hole). We can also account for dynamical friction contributions, in a manner analogous to related effects which one of us has handled in studies of the history of the lunar orbit (Touma & Wisdom 1998). Finally, we can also account for the secular effects of general-relativistic corrections (Appendix B).

THE ALGORITHM TESTED
In this section, we test the Gaussian ring algorithm by comparing its results to the behavior of the equivalent unaveraged N-body system, integrated with the same Bulirsch-Stoer (BS) algorithm at equivalent tolerance. Alternatively, one may match numerically averaged solutions with results of analytical averaging, an exercise which we leave to a forthcoming discussion of analytical theories of softened systems of coupled rings (Mroué & Touma, in preparation).
The desired accuracy of the calculation is specified by two parameters: the first is ǫ quad , the fractional accuracy of the numerical averaging, which was discussed in §3; and the second is ǫint, the tolerance specified in the numerical integrator for the ordinary differential equations.
A technical point is that when comparing secular computations to unaveraged calculations, one should choose the initial conditions of the unaveraged system to ensure that the comparison is not affected by high-frequency oscillations (Laskar & Simon 1988). This can be done by experimenting with various initial conditions for the unaveraged calculations to find one at which the high-frequency oscillation is near a node, or by techniques such as 'warming up' the unaveraged integration (Saha & Tremaine 1992). In the present paper, however, we simply use the same initial conditions for both N-body and averaged calculations.
The simulations below involve stars orbiting a 10 7 M⊙ black hole. Distances are measured in parsecs, time in years, and mass in units of the black-hole mass. In these units, Kepler's law relating the period P and semi-major axis a is P = 2.96 × 10 4 yr (a/ pc) 3/2 . Orbital inclinations are measured with respect to the x-y plane of a right-handed reference frame, with origin at the central black hole; longitudes of node and periapsis are referred to the x-axis of this frame. The dominant relativistic correction (apsidal precession due to the central mass) is sometimes included, as described in Appendix B.

Softened Kozai oscillations
To start out, consider a system that represents a star orbiting one of the components of a binary black hole at the centre of a galaxy. The black-hole companion is represented by a fixed ring having mass Mo = 1, semi-major axis ao = 10, and eccentricity eo = 0.5 (recall that semi-major axis is in units of parsecs and mass in units of the mass of the primary black hole, 10 7 M⊙). The fixed orbital plane of the black hole binary serves as the x-y plane of the reference frame, with the x-axis taken perpendicular to the line of apsides, the y-axis pointing to the periapsis, and the z-axis along the orbit normal. The star is represented by a ring with Mi = 10 −7 , ai = 0.1, and initial eccentricity ei = 0.01. Its orbital plane is initially inclined by Ii = 60 • , cuts the plane of the binary along the x-axis, and has a periapsis which is 90 • from the node. The softening is 0.1ai; relativistic precession is switched off. We first carried out a direct (restricted) three-body integration (BS tolerance ǫint = 10 −12 ) with these initial conditions, supplemented by randomly chosen mean anomalies. Then we performed a number of ring simulations for comparison. At a fixed sampling of 16 sectors per ring, a maximum error of 10 −12 in the Keplerian energy was comfortably achieved. The Gaussian ring algorithm follows the secular evolution with giant strides ranging from 1.4 × 10 7 (ǫint = 10 −14 ) to 4.5 × 10 7 yr (ǫint = 10 −6 ).
Below, we compare results of the secular ring computation at ǫint = 10 −12 with the benchmark particle simulation. In this rather conservative run, the timestep in the averaged simulation (∼ 1.8 × 10 7 yr, 1/40 of the Kozai period) was 2 × 10 5 times larger than in the direct three-body integration (∼ 90 yr, 1/10 of the orbital period).
Neither energy nor angular momentum is preserved in the restricted, eccentric, three-body calculations. In the averaged dynamics, the star's energy is conserved but its angular momentum is not. The energy conservation provides a check on the accuracy of the Gaussian ring algorithm: the ring's secular energy was conserved to a maximum fractional error of 4 × 10 −10 , over a hundred Kozai cycles 5 .
The orbital elements of the star from the three-body and secular codes are shown superimposed in Fig. 1. The eccentricity and inclination variations are due to Kozai oscillations. To the naked eye, the runs are practically indistinguishable in this figure. A closer look reveals two types of (minor) differences between the two codes, as shown in Fig. 2: (i) small-amplitude variations over the binary's orbital period, which average out in the secular dynamics; (ii) relatively larger variations with the same period as the Kozai cycle. The latter are mainly due to differences in the initial conditions between the three-body and secular codes; we were able to reduce them significantly by using averaged initial conditions for the secular dynamics as described above. While the differences in initial conditions are caused by second-order terms in the masses, the time evolution of differences is not; in fact, ring simulations with slight variations in initial conditions revealed similar (qualitative and quantitative) time variations in the differences of their orbital elements.

Counter-rotating rings
The linearized dynamics of softened, planar, counter-rotating, nearly Keplerian systems of stars and gas was considered by Touma (2002). There it was shown that (i) such configurations are prone to violent m = 1 (lopsided) instabilities; (ii) the bifurcation to instability, the growth rate and the pattern speed are functions of softening; (iii) the instability is associated with interaction of modes of opposite energy.
The secular theory employed by Touma (2002) was based on a Hamiltonian that was truncated at O(e 2 ). Mroué & Touma (in preparation) extend the analysis to O(e 4 ), and provide a description of the unfolding of the bifurcation, the appearance of uniformly precessing ring configurations, and the secular dynamics around them. The Gaussian ring algorithm offers an ideal tool for exploring this secular dynamics to all orders in inclination and eccentricity. Results on the planar two-ring problem will be presented in Mroué & Touma. Here, we further test our method by following the development of a nearly coplanar, counter-rotating system of two rings.
Specifically, we follow two equal-mass rings, in initially nearly circular (0.01 eccentricity), nearly coplanar (1 • inclination), counter-rotating orbits. We consider a case in which Mo = Mi = 10 −6 , ao = 1 and ai = 0.5, with softening b = 0.3ao and relativistic precession switched off. This is a linearly unstable configuration, in which the eccentricity of both rings is expected to grow (note that reducing the softening length below 0.26ao switches off the instability, at this ratio of semi-major axes). The evolution is shown in Fig. 3, where one observes oscillations in the eccentricity from near-circular to near-radial and back. Even more dramatic is the inclination dynamics: the rings experience repeated flips between prograde and retrograde near the peaks of the eccentricity oscillations. The cause of this out-of-plane instability is not fully diagnosed, though it smells like Kozai cycles. To be sure, this behavior, which may stimulate bending instabilities in near-Keplerian discs, is markedly 360 Ω (degrees) Figure 1. Osculating elements of a star in a black-hole binary, as determined from a direct three-body integration (solid curve) and from the Gaussian ring algorithm (crosses) over a single Kozai cycle. different from the behavior of counter-rotating configurations studied in the context of the three-body problem. There, it has been understood for decades that counter-rotating configurations are rather stable, even more so than their prograde counterparts. It appears that softening (or heat) modifies the secular frequency spectrum to a point that angular momentum exchange between rings leads to an instability, which first promotes eccentricity growth to a near-radial configuration, then inclination growth to flipped orientations, and back.
In the course of these acrobatics, the secular algorithm satisfied an averaging tolerance ǫ quad = 10 −14 , by working its way up to 128 subdivisions per ring (from a minimum of 16) in the high eccentricity/inclination phases of the evolution. The dynamics was followed with a mean timestep of about 1/100 of the oscillation period, dictated by an integrator tolerance ǫint = 10 −12 . At this tolerance, the relative error in the secular energy was no more than a few times 10 −9 over 100 oscillations; the fractional error in the total angular momentum was bounded by 3 × 10 −11 over that same period.
The reader may wonder how well we do against particle integrations in these vigorously unstable configurations. The situation is illustrated in Fig. 4: there is good agreement over a cycle, giving way to quantitative disagreement but similar qualitative trends. This discrepancy has little to do with proper adjustment of initial conditions. Rather, second-order effects  Figure 2. Difference between the osculating elements of a star in a black-hole binary, as determined from a direct three-body integration and from the Gaussian ring algorithm, over ten Kozai cycles. in the masses, which are neglected in ring simulations, are amplified in close encounters that result in the large eccentricity and inclination phases of this system. This is particularly apparent in Fig. 5, where jumps in particle semi-major axis accumulate to a few percent in less than five cycles. These jumps, which occur at close encounters/crossings between particles, are not accounted for in the secular representation of the dynamics. They bring about shifts in the secular frequencies which are responsible for the discrepancies observed in Fig. 4 6 . 6 In comparisons between ring and particle simulations it should be borne in mind that there are several possible particle analogs to a ring. Here, we associate a single particle with a ring, to highlight the discrepancies that the practitioner should keep in mind when studying secular few-body dynamics, say in the context of planetary systems. Alternatively, one could think of a ring as a guiding-center trajectory that provides the mean phase-space position of a distribution of particles, with total mass equal to the mass of the ring and spatial extent set by the ring's softening; this point of view is more natural to stellar clusters in galactic nuclei. Finally, the ring could represent a point in the five-dimensional space of actions and resonant angles, i.e., a component of a phase-space distribution function that is uniformly distributed in mean anomaly. The error associated with Gauss's method depends on the analogy being used. Finally, while softening/heat can induce instability in a counter-rotating system of rings, general-relativistic precession can suppress the instability. The importance of normally negligible relativistic effects is already appreciated in studies of Kozai oscillations, (e.g., Holman et al. 1997), and we demonstrate that the same applies in this case by switching on relativistic precession in the same configuration (Fig. 6).

COUNTER-ROTATING SYSTEMS OF RINGS
We conclude our battery of tests with a modest, yet challenging, look at systems of counter-rotating rings. For some time now, it has been recognized that counter-rotation is likely to be associated with lopsidedness in self-gravitating systems of gas and stars. Touma (2002) pointed out that counter-rotating instabilities are at work in nearly Keplerian systems and speculated that non-linear evolution of unstable discs may lead to the formation of eccentric discs such as the one observed in the 'double' nucleus of M31 (Tremaine 1995;Peiris & Tremaine 2003). We limit our discussion to a couple of runs, which are illustrative of both the challenges faced by the algorithm, and the typical dynamical fate of such systems.
We follow the dynamics of a system of 100 rings, half of them prograde with a total mass of 10 −3 relative to the central point mass, and the other half retrograde, with the same semi-major axis distribution as the prograde rings, and 20 per cent of their mass. In this illustrative example, we picked the rings from uniform distributions in semi-major axis (over the interval points per ring and a maximum of K = 512, enforcing an averaging error ǫ quad = 10 −10 . The system was followed for 30 Myr. During this period, the secular energy was conserved to within a fractional error of 10 −9 , with ten times smaller error in the total angular momentum. The run took about two weeks on a single 3.2 Ghz processor. An order of magnitude improvement in performance can be achieved if one is willing to live with a tolerance ǫint = 10 −6 . The evolution of the system in eccentricity and inclination is shown in Fig. 7. In the first 5 Myr the system undergoes significant growth in eccentricity, at near constant inclination; by 5 Myr the mean eccentricity of the prograde population is around 0.6 and the mean eccentricity of the less massive retrograde population is larger at 0.85. This result is not new: a closely related planar system is known to be linearly unstable (Touma 2002), with the lighter/retrograde component experiencing a more significant growth of eccentricity than the heavier/prograde component. It is the later fully non-linear behavior that we can confidently follow with the Gaussian ring algorithm for the first time. Following the initial growth of eccentricity, and in a manner reminiscent of the behavior witnessed in the previous section for two rings, the rings in the retrograde/lighter component experience a dramatic growth and spread in inclination, practically reversing orientation by 7.8 Myr. The mean eccentricity of both populations increases slightly during this flipover phase, with the retrograde rings  As the time increases, from top left to bottom right, one witnesses the growth of eccentricity at near constant inclination, followed by the flip-over of the retrograde rings, then their attempted return to their initially coplanar configuration. The more massive prograde rings experience a significant net growth in their mean eccentricity, while remaining clustered close to their initial inclination. Figure 8. Evolution of the projected and smoothed surface density of the full system. The projection is on the reference x-y plane, lumping together contributions from prograde and retrograde orbits. The central point mass is at the origin, marked by a cross. Color refers to density levels. An m = 1 mode grows as a result of the counter-rotation instability, capturing rings in its way, with the prograde population maintaining the precessing single-lobed pattern, all the while the initially retrograde population prepares for disintegration. Figure 9. The projected surface density evolution of the rings over longer times, when viewed in negative color. The highly eccentric and dispersed population of initially retrograde rings projects onto a lightly hued halo, around the heavier lopsided configuration of prograde rings, captured in a precessing m = 1 mode. reaching a mean of 0.9 around 6 Myr. Then, the initially retrograde rings try to regain their home-base; they almost succeed around 10.3 Myr, but suffer significant losses along the way. At this point their mean inclination is about 140 • , and their mean eccentricity has decreased to 0.6; the prograde rings are still in low-inclination orbits, with a mean eccentricity of 0.45. The prograde rings are also bunched in the conjugate angles, the longitude of periapsis in particular. This is apparent in the associated spatial densities, shown in Fig. 8. The non-linear saturation of the instability leads to an oscillating, eccentric ring composed of prograde particles, with a mean inclination of 7.5 • , a mean eccentricity of around 0.6, and an increasingly steady prograde precession with a frequency 0.8 rad Myr −1 (Fig. 9). Meanwhile the culprits behind the instability, namely the initially retrograde rings, are dispersed into a population with a wide spread in inclination, node and argument of periapsis, and a mean eccentricity of 0.7 (Fig. 10).
In short, the mode, which is ignited by the instability, grows through resonant capture of both prograde and retrograde populations, then saturates as the retrograde population (which has a lower angular momentum) escapes from resonance, and disperses through (partially stochastic) phase mixing; this leaves a captured prograde population which maintains a lopsided  m = 1 mode that is increasingly steady in its (prograde) precession. Just for fun, and once the mode was fully established, we removed all initially retrograde rings, and left the prograde population to itself. Perhaps not surprisingly, the captured prograde population, shown in Fig. 11, maintained its cohesive state, now precessing uniformly with a slightly higher frequency of 0.9 rad Myr −1 for at least four precession periods. At a mean timestep of 1.1 × 10 4 yr, the integration kept fractional errors in energy and angular momentum below 10 −9 . One could think of the resulting mode as an analog of the lopsided mode which was artificially excited by Jacobs & Sellwood (2001). The main differences are: (i) the Jacobs-Sellwood simulation was a two-dimensional N-body simulation, while ours is a three-dimensional N-ring simulation; (ii) the Jacobs-Sellwood simulation used N =100,000 while ours used N = 100 (each of our rings should be regarded as representing many stars); (iii) the Jacobs-Sellwood simulation formed an eccentric disc using eccentric initial conditions, while we relied on counter-rotation to generate an eccentric equilibrium through the non-linear evolution of an instability.
A careful treatment of the growth and saturation of the counter-rotating instability would take us too far afield. It is in fact the subject of an upcoming study by Touma and Kazandjian (in preparation), who perform sufficiently resolved N-body simulations to permit a detailed phase-space dissection of the dynamics in question. Their results, which are based on long Figure 11. The evolution of the prograde rings after the retrograde rings were removed, leaving a slightly thickened, lopsided, and uniformly precessing annulus.
term N-body simulations of up to 500,000 counter-rotating particles around a dominant central body, are in remarkable qualitative agreement with the salient features of the 100 ring simulations reported here. Still, the full fledged confrontation of our algorithm with a conventional N-body code will have to wait for ring simulations that are sufficiently resolved (10 5 rings, pursued over a hundred precession periods) to permit detailed phase-space analysis, and comparison of the dynamics. These simulations are surely within the reach of a parallel version of our algorithm, running on a state-of-the-art cluster, and we hope to report on such undertakings in the near future.

DISCUSSION
The algorithm we have described computes the secular or orbit-averaged evolution of systems of particles in a Keplerian force field and interacting via softened gravity. It is accurate to first order in the masses of the particles, and to all orders in eccentricity and inclination. The order for a system of N particles, which is 6N in a conventional N -body code (assuming that momentum conservation is used to solve for the motion of the central body), is reduced by 2N as a result of averaging and the consequent conservation of semi-major axis. However, computing with angular momentum and eccentricity vectors forces a redundancy as it resolves a degeneracy. In other words, in order to facilitate the traversal of high eccentricity and inclination regimes, we simulate a redundant 6N -dimensional system of equations. The algorithm uses analytic averaging over the orbit of the perturbing particle (Gauss's method) and numerical averaging over the orbit of the perturbed particle. The timesteps are now fractions of secular time-scales, which are at least a factor of M⋆/(N m) longer than steps taken in unaveraged simulations of equivalent accuracy. The price for this is that our method cannot capture the effects of mean-motion resonances, which appear in the secular dynamics at second order in the masses.
The Gaussian ring algorithm provides us with a powerful numerical tool for the study of secular dynamics of nearly Keplerian N-body systems, especially those with large eccentricities and inclinations (so that perturbation expansions converge slowly or not at all) and/or crossing orbits. Our method resolves the dominant, secular evolution of these particles with much greater accuracy and speed than conventional N-body codes: in the examples we have examined the typical timesteps are two to three orders of magnitude larger, depending on the mass of the N-body cluster relative to the central mass (although the calculations per timestep are more lengthy). Moreover a single ring should be thought of as containing many stars, so the N required for reliable results should be much smaller.
We have discussed the implementation of the method ( §3). The accuracy goals of the algorithm are encapsulated in two parameters, ǫ quad and ǫint, the first expressing the desired fractional accuracy of the numerical quadrature of the forces over the perturbed ring, and the second expressing the desired accuracy of the integration of the equations of motion. The accuracy of the quadrature can be assessed by evaluating the rate of change of the semi-major axis ( §3.1), which should be zero in secular dynamics. Provided one can live with moderate softening (b 0.01a), between 16 and 128 force evaluations per ring should be sufficient to yield ǫ quad ∼ 10 −11 ; 10 to 100 times that number of force evaluations may be required for extreme configurations or very small softening.
Much can be done to enhance the efficiency and accuracy of the method: (i) adaptive quadrature routines can reduce the number of force evaluations required for ring pairs experiencing a close approach; (ii) geometric integrators for this most illustrious of Lie-Poisson systems should enhance the long-term accuracy of the numerical integration; (iii) the cost of quadratures makes the method eminently parallelizable in a distributed environment.
JT acknowledges the generous hospitality of Stephane Colombi at the Institute for Astrophysics in Paris, and a visiting position at the Institute for Advanced Study. We also acknowledge support from NSF grants AST-0206038, AST-0507401, and AST-0807432, and NASA grant NNX08AH24G. We thank Abdel Hussein Mroué, S. Sridhar, and Michael Balabane for discussions.

APPENDIX A: HALPHEN'S CONE: A DELAYED GEOMETRIC INTERLUDE
where ds is an element of surface area and the integral is over the surface area bounded by the curve C, i.e., the area swept out by the vector r − r⋆ in one orbit of the perturber. Now consider the cone having vertex at the origin (the attracted particle) that passes through the curve C, i.e., the set of all vectors parallel to r(l), 0 l < 2π. The curve C lies on this cone, and if it is smoothly distorted to some other curve on the cone the additional contribution to Fx will be zero, because the normal to the cone ds at r is perpendicular to r. We conclude that the integral in equation (A2) is independent of the curve C so long as it lies on the cone. We may therefore choose a curve C on the cone that simplifies the line integrals in (A2).
We note that y(0) = (1 − e 2 ) 2 a 2 h 4 > 0, y(h 2 ) = −e 2 y⋆ 2 h 4 < 0, which implies that the roots are real, one negative and two positive, and we may therefore assume that µ1 > µ2 > 0 > µ3. The eigenvectors of W may be written where θ k is real, chosen so that h k ·h k = 1. We assume that θ1 and θ2 are positive and then choose the sign of θ3 so thatê1, e2,ê3 form a right-handed triad, that is,ê3 =ê1×ê2. Then the orthogonal matrix O is defined by Oij = e j i and det O = +1. The equation of the cone in the new coordinates is We choose the curve C to be the intersection of the plane u3 = const with the cone. A parametric equation for C is then u1 = αu3 cos ψ, u2 = βu3 sin ψ, α = |µ3|/µ1, β = |µ3|/µ2.
The curve C as defined in equation (A2) is oriented counter-clockwise as viewed by a distant observer on the positive angular-momentum axis (the z axis). This is the same orientation as the curve defined by equation (A13) if and only if h and u3 (the distance of the integration plane from the origin) have the same sign. To account for this we assume u3 > 0 and replace the pre-factor 2πh in equation (A2)  We may correct Newton's equations with the dominant secular contribution from general-relativistic effects. These effects are mainly felt by nearly radial orbits on their close encounters with the central black hole. When averaged over the orbital period these corrections contribute a Hamiltonian The gradients are ∇ L Hgr = 0 and The secular equations of motion are (Tremaine et al. 2008) To this order, the main effect is on precession of periapsis, leaving eccentricity, orientation and semi-major axis of the orbit unaffected dA dt = 6πGM⋆ c 2 aP (1 − e 2 ) 3/2 L×A.