{'fraction_non_alphanumeric': 0.026296501943364796, 'fraction_numerical': 0.011771238200999444, 'mean_word_length': 4.676500252143217, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An energy-conserving formalism for adaptive gravitational force softening in SPH and N-body codes

In this paper we describe an adaptive softening length formalism for collisionless N-body and self-gravitating Smoothed Particle Hydrodynamics (SPH) calculations which conserves momentum and energy exactly. This means that spatially variable softening lengths can be used in N-body calculations without secular increases in energy. The formalism requires the calculation of a small additional term to the gravitational force related to the gradient of the softening length. The extra term is similar in form to the usual SPH pressure force (although opposite in direction) and is therefore straightforward to implement in any SPH code at almost no extra cost. For N-body codes some additional cost is involved as the formalism requires the computation of the density via a summation over neighbouring particles using the smoothing kernel. The results of numerical tests demonstrate that, for homogeneous mass distributions, the use of adaptive softening lengths gives a softening which is always close to the `optimal' choice of fixed softening parameter, removing the need for fine-tuning. For a heterogeneous mass distribution (as may be found in any large scale N-body simulation) we find that the errors on the least-dense component are lowered by an order of magnitude compared to the use of a fixed softening length tuned to the densest component. For SPH codes our method presents a natural and elegant choice of softening formalism which makes a small improvement to both the force resolution and the total energy conservation at almost zero additional cost.


INTRODUCTION
This paper is concerned with the question of how best to represent the gravitational force when simulating self gravitating systems using particles. The simplest of such systems is a collection of stars which is usually replaced by a very much smaller number of computational particles. A more complicated example is the simulation of self gravitating gas with or without a stellar component, using the particle method SPH (Monaghan 2005).
Provided the number of computational particles is sufficient to resolve the important dynamical scales the simulation can give satisfactory results for most quantities though the slow relaxation of a galaxy is number dependent. In the case of N −body simulations it is, however, necessary to soften or smooth the forces between pairs of particles so that the binary collisions of the computational particles will not cause numerical artifacts.
The simple Plummer form of the softening where the force F (r) between a particle pair with masses ma and m b separated by distance r is and h is the softening length. Dehnen (2001), amongst others (e.g. Dyer & Ip 1993) has shown that a better choice is to use Kernel smoothing with kernels W (r, h) that have compact support. The softened force at a due to b then takes the form Provided h is constant Poison's equation shows that the softening is equivalent to calculating the local gravitational force on a point particle a due to a density or, when there is a collection of particles contributing to the force on particle a, the density is which is identical to the SPH density estimate. Kernels with compact support are zero beyond some specified distance proportional to the length scale h, and the pair force then has the correct value for the two sets of real particles represented by two computational particles.
The usual practice is to use a fixed value of h for all of the N −body particles. A key issue which arises in this context, and which has been the subject of a number of studies (Merritt 1996;Romeo 1998;Athanassoula et al. 2000;Dehnen 2001;Rodionov & Sotnikova 2005), is the 'optimal' choice of softening length, for too small a softening length will result in noisy force estimates, whilst too large a value will systematically 'bias' the force in an unphysical manner. In general, however, the 'optimal' choice depends on particular system under investigation and may not be known a priori (Athanassoula et al. 2000). Dehnen (2001) quantifies the errors arising from both Plummer and kernel softening of the above form. In all cases he finds the accuracy is improved if h is allowed to vary according to the local particle number density n in such a way that h is smaller when n is larger. Typically h ∝ 1/n 1/3 . For self-gravitating SPH calculations using a softening length which differs from the smoothing length can lead to unphysical results (Bate & Burkert 1997).
To retain conservation of linear and angular momentum it is necessary to use a symmetric form of F so that each particle in a pair interaction experiences an equal but opposite force. This can be achieved by using, for examplē h = 1 2 (hi+hj) in place of h in (2). However, because the softening length then varies in space the total energy of the system will not be conserved. The errors are often not large but can lead to secular increases in the total energy of the system, destroying the phase-space conservation which is crucial for accurate N −body simulation (Hernquist & Barnes 1990;Dehnen 2001;Rodionov & Sotnikova 2005).
In this paper we show how a Lagrangian for a self gravitating gas can be devised which has the softening of the force and the variation of h built in. The advantage of using a Lagrangian is that, provided it is constructed correctly, the conservation laws are automatically satisfied. In particular the conservation of energy and momentum is exact though, in practice, the accuracy is determined by the time stepping algorithm. The new equations of motion have an extra term in addition to the standard SPH and gravity terms. It is this term which guarantees energy conservation. We apply our algorithm to both static and dynamic problems. In some cases the new equations give results which are very similar to results obtained previously, but in some cases the results are very much improved.

KERNEL SOFTENING
A general formulation for force softening was given by Dehnen (2001) and we use a similar formulation here. The modified gravitational potential per unit mass may be written in the form where φ is the softening kernel which is a function of the particle separation and the softening length h (we use h to denote the softening length since it corresponds with the smoothing length used in the SPH density estimate). The kernel determines the functional form of the modified gravitational force law. For example, in the case of Plummer softening the softening kernel is given by Neglecting the spatial variation of h the force estimate corresponding to (5) is given bŷ where φ ′ = ∂φ/∂|r − r b |. The underlying smooth density field can be obtained from Poisson's equation giving where the density kernel is related to the softening kernel according to The kernel density given by (9) corresponds to the mass distribution of each particle being smoothed. Readers familiar with SPH will notice that (9) corresponds to the density estimate used in SPH calculations, where W is the usual SPH smoothing kernel. In general, the functional form of the softening kernel may be specified for either the potential term φ, the force evaluation term φ ′ or W . In each case the corresponding kernel for the other cases may be determined by a straightforward integration or differentiation. For example, in Nbody codes, it may be preferable to choose a kernel primarily for the force evaluation, from which the functional form of the potential and density kernel can be derived. In SPH the kernel is primarily used for the density estimate, where the most commonly used form is the cubic spline of Monaghan & Lattanzio (1985) where q = r/h. The corresponding force kernel is given by the functional form of which is given for the cubic spline in Appendix A. The softening kernel for the gravitational potential may be calculated from the force kernel using the form of which is also given in Appendix A for the cubic spline. For general kernels (13) combined with (12) can be integrated by parts to give where the last term is the constant of integration, determined by the requirement that φ → 0 as r → ∞ (note that we have also assumed a kernel with compact support of size 2h). The modified potential, force functions and the density kernel are shown in Figure 1 for the cubic spline. The reader should note that whilst we use the cubic spline as an example throughout this paper, the results derived in the following sections are quite general and any smoothing kernel may be used (including any of those suggested by Dehnen 2001).

LAGRANGIAN FORMULATION
The Lagrangian describing a self-gravitating gas is given by where Φ is the gravitational potential (5) and u is the thermal energy per unit mass. The equations of motion may be obtained via the Euler-Lagrange equations giving The advantage of using a Lagrangian to derive the equations of motion is that, provided the Lagrangian is symmetrised appropriately, momentum and energy conservation are guaranteed. Variational principles have been used extensively to derive conservative SPH formalisms for relativistic fluid dynamics (Monaghan & Price 2001), magnetohydrodynamics (Price & Monaghan 2004) and in the case of a spatially variable smoothing length (Springel & Hernquist 2002;Monaghan 2002). An adaptive softening length formalism may be derived by writing the gravitational part of the Lagrangian in the form where φ bc refers to φ(|r b − rc|). Swapping indices in the double summation shows that (17) is equivalent to averaging the softening kernels in the form The derivative of (17) is given by where We relate the smoothing length to the particle co-ordinates, assuming h = h(ρ), using where ρ is the density calculated by a summation over neighbouring particles in the form where W is the density kernel. The relationship between h and ρ means that this is a non-linear equation for both ha and ρa which can be solved self-consistently for each particle. The iterative method we use for doing so is described in detail in §4.2. The spatial derivative of (22) is where W is the density kernel and Ω is a term accounting for the gradient of the smoothing length given by Using (20), (21) and (23) in (19) and simplifying, we have The quantity ζ is defined as where ∂φ/∂h can be tabulated (or calculated) directly for the particular smoothing kernel used. For the cubic spline the expression is given in Appendix A. The derivation of the SPH pressure force from the thermal energy term in the Lagrangian (14) in the case of a spatially variable smoothing length has been described in detail elsewhere (e.g. Springel & Hernquist 2002;Monaghan 2002;Price & Monaghan 2004) and we simply use the result here. The final equations of motion take the form The first term in (27) corresponds to the softened gravitational force. The second term is present only in the case of adaptive softening lengths and it is the incorporation of this term which restores the energy conservation. The third term is the usual SPH pressure force allowing for a spatially variable smoothing length. The terms Ω and ζ required in the adaptive softening term (and for Ω also in the pressure force) are easily calculated alongside the density summation. The additional adaptive softening length term can be seen to have the same form as the pressure force, with the quantity P/ρ 2 replaced by ζ. Notice however that ζ is, for positive kernels, a negative definite quantity and therefore that the adaptive softening term acts in opposition to the usual pressure term (that is, in the direction of increasing the gravitational force). This is in line with the recent findings of Hubber et al. (2006); that SPH always underestimates the gravitational force at low resolution. In fact they suggest adding an additional contribution to the gravitational force based on the 'self-gravity' of an SPH particle. The new term derived above provides a similar contribution without the need for ad-hoc prescriptions.
Alternative formulations of the adaptive softening formalism given above are possible by symmetrising the Lagrangian in different ways. We use the formulation given above since it is simple and efficient to implement. As an example, we derive an alternative version based on the average softening length in Appendix B. Whilst the force derived using the average softening length is very similar to (27) but in keeping with the average h used in the variable smoothing length SPH formalism of Benz (1990), it has the practical disadvantage that we do not use the average smoothing length elsewhere in the calculations, neither in the density summation (since the density for particle a is calculated using only ha) nor in the SPH force term (see 27). Thus the calculation of the quantityζ (B6) for each particle requires the calculation of the kernel not only using ha (for the density and Ω) but additionally usingh ab which is not only inefficient but also rather inelegant in the numerical code. Thus, we do not use the average h formalism in this paper.
A further possibility, not examined in this paper, would be to use the product of the softening kernels in the force evaluation. The situation is complicated slightly in this case as the product form must either be used in the potential or force but not both. In any case the differences in the force evaluated using different symmetrised forms are very small. The suggestion put forward by Dyer & Ip (1993) that the force should be symmetrised by considering two overlapping spheres may also be used in a similar manner to derive an energy conserving formalism, but it is not clear that there is any advantage to be gained by doing so (see Dehnen 2001).
For reference the consistent forms of the continuity and internal energy equations for SPH simulations are given by and where v is the particle velocity. The continuity equation can be used to make a starting guess for the h iteration procedure used to determine the density (described in the following section). An alternative to using the internal energy equation is to evolve the entropy as an independent variable (Springel & Hernquist 2002) which is possible for ideal equations of state.

NUMERICAL TESTS
We test the adaptive softening length formalism derived in the previous section using three examples. The first ( §4.3) is a series of static tests used by Dehnen (2001) and Athanassoula et al. (2000) in order to estimate the force errors associated with softening formulations. We also consider a dynamic version of one of these tests in order to study the energy conservation properties of our new method ( §4.3.3). The second example ( §4.4) involves self-gravitating SPH and the static structure and dynamical oscillation of a polytrope.

Errors
In the static halo tests, we calculate the Average Square Error (ASE) in the gravitational force according to where fi is the force on particle i, N is the particle number and C is a normalisation constant. Unless otherwise specified we use C = 1/f 2 max where fmax is the maximum value of the exact solution. The Mean Average Square Error (MASE) is then the mean over all realisations, We choose this quantity rather than the Mean Integrated Square Error (MISE) used by Merritt (1996) and Dehnen (2001), given by where ρ(x) is the true density at a point x, f (x) is the force calculated at that point from the N −body distribution and M is the total mass. Calculation of the MISE is complicated by the need to integrate along radial grid points (involving calculation of the force at positions other than particle positions) and Athanassoula et al. (2000) find little difference between their results using MASE or MISE error measures. In our case the correction terms derived in §3 depend on a particle's own density estimate, so it makes sense to calculate errors only at particle positions (that is, using the MASE estimate). The reader should bear in mind that, using either the MASE or MISE as defined above, the total error tends to be dominated by the regions containing the largest forces. This can be somewhat misleading in comparing adaptive softening with fixed softening, as the fixed softening length is generally chosen to minimise the error in the densest regions, where the adaptive softening will not show a large difference. An example is given in §4.3.2 where a two-halo system is setup and we explicitly show the contribution to the MASE from each halo, even though the total error is dominated by the densest component.

Setting the softening length
The method we use for setting the softening length is identical to the method used by Price (2004) (see Price & Monaghan 2004) for setting the smoothing length in SPH calculations. A similar method is also used by Springel & Hernquist (2002) and hence also in the publicly available GADGET-2 code for N -body and SPH (Springel 2005). The idea is to regard the smoothing length as a function of density via the relation or more specifically where m is the particle mass, ρ is the mass density and η is a dimensionless parameter which specifies the size of the smoothing length in terms of the average particle spacing (similar to the parameter ǫ used by Dehnen 2001). The derivative is given by An equivalent interpretation of (34) is that the mass contained within a smoothing sphere is held constant (Springel & Hernquist 2002), that is where σ is the compact support radius of the kernel (= 2 for the cubic spline) and N neigh ≡ 4 3 π(ση) 3 may be used as an approximate measure of the number of neighbours contained within a smoothing sphere. Unless otherwise specified we use η = 1.2 in the variable smoothing/softening length formulations used throughout this paper, which in three dimensions is equivalent to ∼ 60 neighbours.
For a pure N-body simulation using unequal mass bodies, it may be advantageous to use a number density rather than the mass density for setting the softening length. The resulting gravitational force in that case is identical to the first two terms in (27) with mass density replaced by number density. In this paper we assume a mass density dependence consistent with SPH simulations.
The density is calculated by a direct summation over the particles in the form (22) which, via the relation (34) becomes a non-linear equation to be solved for both h and ρ. Dehnen (2001) suggests that even a rough approximation of the (number) density is sufficient for the purpose of adapting the softening length in N −body calculations. A similar argument may be made for setting the smoothing length in SPH calculations. In both cases, however, the situation changes once the gradient terms are incorporated into the equations of motion (as in this paper) since these terms are calculated on the basis of the h(ρ) (or h(n)) relation and may therefore introduce substantial inaccuracies into the solution if the density and smoothing (or softening) length are far from being consistent with (34).
From a practical point of view obtaining a selfconsistent solution to (34) and (22) is a relatively straightforward root-finding problem. The function to be solved may be written in terms of either h or ρ. Written in terms of h we have in the form where ρa is the density consistent with the current smoothing length ha calculated from the relation (34) and ρsum is the density calculated using ha from the summation over neighbouring particles (22). We use a Newton-Raphson iteration method, ie.
where the derivative of (38) is given by We find this method to be efficient and cost effective, particularly since the quantity Ω (defined in (24)) is already calculated alongside the density summation for use in the equations of motion. Convergence is determined for each particle individually according to the criterion |hnew − h|/h0 < ǫ where h0 is the smoothing (or softening) length at the start of the iteration procedure and typically we use ǫ = 10 −3 . We find that it is more efficient to perform the iterations by looping over the particles as the outer loop and iterating each particle at a time to convergence. We also find that it is no longer efficient to store a global neighbour list for all particles but rather to perform a neighbour search on-the-fly (e.g. using a treecode), recalculating where necessary and being stored only for the particle being iterated. This also represents a significant reduction in memory requirements for SPH calculations.
The Newton-Raphson iterations work extremely well provided that the initial estimate of h is reasonably close to the actual solution. This is almost always the case in the calculations since there is only a small change in ρ between timesteps. The only problems which may arise are in the first iterations on the initial conditions where h and ρ may be far from the relation (34). For this reason it is useful to revert to a bisection scheme (which is guaranteed to converge) in the case where the Newton-Raphson iterations do not converge (we set the limit for this as > 20 iterations).
In terms of cost the density iterations add only a small amount of extra work to SPH calculations. The exact work required depends on the nature of the simulation since more iterations are required when the density is changing rapidly. However the scheme is very efficient since iterations are only performed on the subset of particles whose densities are changing. The scheme can be made still more efficient by predicting an initial guess for h in the time evolution scheme using the time derivative where for the summation (22) the time derivative is given by (28). Using a prediction step we find that in general (although dependent on the dynamics of a particular simulation) only a small fraction of the particles require extra density calculations and that these particles then converge rapidly (in ∼ 2 − 3 iterations).

Isolated haloes
The first test we perform is to compare the (softened) gravitational force to the exact force given an analytic density profile (corresponding to typical structures formed in cosmological simulations or used as initial conditions in galaxy models). Following Dehnen (2001) we consider two different density profiles -Plummer spheres and Hernquist models. The density profile for the Plummer spheres are given by where M is the mass and rs is a parameter determining the concentration of the halo. The corresponding gravitational potential and force are given by and The cumulative mass profile for the Plummer sphere is given by The Hernquist (1990) model consists of a density profile given by The density and gravitational forces in this model are more difficult to resolve since the density profile is cusped near the origin. The potential and force are given by whilst the mass profile is given by The density profiles in each case are set up in the usual manner choosing three random deviates (x1, x2, x3) uniformly on (0, 1). The first is used as a position in the mass profile from which the radial coordinate is determined by rearranging M (r) to give r(m) where m is the mass fraction. In practice we use only mass fractions < 0.99 in order to prevent isolated particles being placed at extremely low densities. The second random number x2 is used to give a random azimuthal angle ϕ = π(2x2 − 1) whilst the third is used to give a spherical angle θ via the transformation θ = cos −1 (2x3 − 1) (necessary to prevent the distribution from clumping towards the poles). The result is a particle distribution which closely mirrors the analytic density profile although with errors decreasing like 1/ √ N due to the Monte Carlo nature of the distribution.
In the numerical simulations we use units of mass [M ] = 1, length [R] = 1 and time [τ ] = (GM/R 3 ) −1/2 . In these units GM = 1 such that the gravitational constant does not appear in the numerical equations. Correspondingly, force and energy (both per unit mass) are measured in units of GM/R 2 and GM/R respectively. For calculation of the mean error we compute 3 × 10 6 /N realisations for a halo of N particles.
The MASE calculated for Plummer haloes of N = 10 2 , 10 3 , 10 4 and 10 5 particles with M = 1 and rs = 1 are shown in Figure 2. The left panel shows the results using a fixed softening length, comparing both Plummer (solid lines) and cubic spline (dashed lines) softening kernels. The right panel shows the results using (cubic spline) adaptive softening, with (dashed line) and without (solid line) the energy-conserving term, also showing the variation with the "Number of Neighbours", by which we mean the parameter N neigh defined in (36).
Some general features are worth pointing out. Firstly, using a fixed softening length, there is a large variation in the total error depending on the choice of softening length (left panel). For softening lengths too small, the error is dominated by noise, reaching a maximum value once the softening length is smaller than the smallest particle separation. For softening lengths too large, the error is dominated by the bias in the force introduced by the softening procedure. For some intermediate choice of softening length, there is a balance between noise and bias which produces a minimum error. This gives rise to the concept of 'optimal softening' introduced by Merritt (1996) and examined in detail by Athanassoula et al. (2000) and Dehnen (2001), whereby the softening length can be 'fine tuned' for a particular simulation for minimum error. In principle this means that, for every N −body calculation, there is an optimal choice of softening length. The problem, demonstrated in the left panel of Figure 2, is that this 'optimal' choice not only depends on the parameters of the problem (for example, Figure 2  shows that the optimal choice changes with resolution, but Athanassoula et al. (2000) also discuss the dependence on the shape and degree of central concentration of the halo), but may also be different for different components of the same simulation.
By contrast, use of adaptive softening (right panel) shows only a weak dependence on the choice of the (adaptive) softening parameter (provided that the neighbour number is small compared to the total number of particles). To guide the reader, our fiducial choice of N neigh ≃ 60 (given by η = 1.2 in (34)) is marked by the open circle in each case. Making this reasonable choice in all cases gives a softening which (according to the MASE estimate) is close to the optimal choice of fixed softening. The exception is perhaps for the 10 4 particle halo, where the optimal choice of fixed softening gives a MASE ∼ 2.5× lower than that using adaptive softening with energy conservation. However, changing the fixed softening length up or down by a factor of two in either direction (ie. not a large range if the optimal choice is not known a priori ) means that even in this case that adaptive softening, even with the energy conservation term added, wins.
Comparison of the dashed (Plummer kernel) and solid lines (cubic spline kernel) in the left panel of Figure 2 confirms the conclusion reached by several authors that use of the cubic spline kernel is advantageous over the standard Plummer softening. In particular the optimal error for each halo is reached at a higher softening length for the cubic spline kernel and the slope in the error curve at high softening lengths is steepened, demonstrating that the cubic spline reduces the bias in the force estimate. This is a result of the compact support of the cubic spline kernel, which gives a force with zero bias outside of the kernel radius (see discussion in §1).
The right panel of Figure 2 also shows the influence of the new energy-conserving term in (27) on the force errors in a static configuration. This term appears to increase the noise but also lower the bias in the (adaptively softened) force estimate, meaning that the total error is greater for smaller neighbour numbers but lower for larger neighbour numbers. We attribute this to the fact that the extra term is related to the gradients in softening length: Where these gradients are spurious (due to noise), the extra term may increase the total error. Where the gradients are due to actual gradients in the density, the extra term correspondingly leads to a more accurate force estimate. This conclusion is also borne out by the results using Hernquist models (right panel of Figure 3). Here the extra term leads to a smaller MASE (compare the dashed and solid lines in the right hand panel of Figure 3) at a lower N neigh than in the Plummer case (dashed vs solid lines in the right hand panel of Figure 2). The density profile in the Hernquist model is strongly cusped near the origin, meaning that any improvement in the resolution of density gradients (e.g. from the new term) tends to improve the error estimate.
The Hernquist model was computed using M = 1, rs = 1 and N = 10 2 , 10 3 , 10 4 and 10 5 particles. The results using fixed Plummer (solid lines) and cubic spline (dashed lines) softening on the Hernquist model are shown in the left hand panel of Figure 3. The differences between Plummer and cubic spline softening are much smaller in this case than for the Plummer spheres (Figure 2), apart from a factor of ∼ 2 difference in the optimal choice of softening length for each kernel (that is, the optimal softening length for Plummer softening is approximately ∼ 1/2 of the optimal value using cubic spline softening).
These tests demonstrate that, for an isolated halo, the use of adaptive softening gives force errors which are close to optimal. Whilst there is not a significant improvement in the MASE compared to the use of an optimally-chosen fixed softening length, the use of adaptive softening removes the need for such fine-tuning. Furthermore it may not be possible to find a softening which is 'optimal' for all components of a simulation. In the following section we consider such an example, where the use of adaptive softening shows a clear improvement.

Two Plummer spheres
Next we consider two Plummer spheres placed at a fixed distance from each other, of equal mass but where one halo is much denser than the other. This situation may be representative of two haloes present in a typical cosmological N −body simulation or in a simulation of galaxy dynamics where more than one galaxy is present. A similar test was considered by Athanassoula et al. (2000) where a variety of mass ratios were also examined. Here we simply choose one representative case.
Both haloes are Plummer spheres, setup as described in §4.3.1. We use equal mass spheres with M = 0.5. The first sphere is placed at the origin, with concentration parameter rs = 1 whilst a second sphere with rs = 0.1 (ie. much denser) is placed some distance away at [x, y, z] = [10, 0, 0]. The MASE is calculated un-normalised in this case, that is with the normalisation factor C = 1 in order to make a meaningful comparison between the errors in each component.
The MASE resulting from 300 realisations of this configuration using a total of 10,000 particles (5000 per sphere) is shown in Figure 4 (solid line), using fixed cubic spline softening (left panel) and adaptive softening with the energy conservation term included (right panel), showing the variation with softening length in the former case and N neigh in the latter (where again the open circles correspond to our fiducial choice of N neigh ≃ 60). The total MASE (solid line) is completely dominated by the densest component, with results comparable to those shown in Figure 2. However, we also plot the contribution to the total MASE from the lessdense component (that is, the rs = 1 sphere at the origin) (dashed line).
The problem with the use of a fixed softening length in a general N −body simulation is evident from Figure 4, namely that the 'optimal' choice of softening length differs for each component. The choice which minimises the errors in the densest component produces errors in the least dense component that are over 1 1 2 orders of magnitude larger than the optimal choice of softening for that component. Conversely choosing the softening which is optimal for the least dense component produces a bias in the force estimate in the densest component leading to a MASE ∼ 2 orders of magnitude larger than the optimal choice of softening length in the dense component. Usual practise is therefore to choose the softening which minimises the softening in the densest component(s) (since this represents the largest contribution to the total error). However frequently one is interested in the properties of both (or all) components in an N −body simulation. This leads naturally to a need for adaptive force softening.
The results using our adaptive softening length formalism (including the energy conservation term) are shown in the right hand panel of Figure 4. For both components the resulting error is, as previously, close to the optimal choice of fixed softening, but here the softening is close to optimal for both components. This means that the force errors in the least dense component are ∼ 1 order of magnitude smaller than using a fixed softening length tuned to the densest component.

Halo relaxation
An extension to the static halo test is to examine the dynamic influence of the energy-conserving term. The initial conditions for this test are a Plummer sphere with an isotropic velocity distribution corresponding to a (dynamic) steady state. The particle velocities are setup in the manner described by Aarseth, Henon & Wielen (1974): The distribution function is where f (r, v, t)drdv is the total mass of particles with position r and velocity v at time t and E is the energy per unit mass of a body: The distribution function is sampled by scaling the velocities in terms of the maximum velocity at r, ie. the escape velocity Writing q = v/vesc, from (51) and (50) the probability distribution for q is proportional to where |q| < 1. This distribution is sampled using the Von Neumann rejection technique (Press et al. 1992): two uniform random deviates x4 and x5 are drawn. Noting that g(q) is always less than 0.1 (since |q| < 1) we adopt q = x4 if 0.1x5 < g(q), otherwise a new pair of random numbers is tried until the inequality is satisfied. The velocity modulus v is obtained using (52) and, using two more uniform random deviates x6 and x7, the velocities are given by The halo is evolved forwards in time using a standard second order leapfrog integrator with a global timestep controlled by the condition where h is the softening length, f is the force per unit mass and the minimum over all particles is used. The energy conservation during the evolution of the equilibrium halo is shown in Figure 5. Using adaptive softening without the additional term (solid line), fluctuations in the energy are observed from the changes in softening length which, although small, dominate over the errors due to timestepping. With the energy conserving term added (dashed line), only a small non-conservation of energy remains which can be shown to decrease as the timestep is made shorter.
As a slightly more demanding test, we also consider the relaxation of a perturbed Plummer sphere -that is, with the velocities drawn from the equilibrium distribution function as described above, but then multiplied by a factor of 1.2. This means that the halo initially expands before slowly relaxing into a dynamical equilibrium state. The evolution of the total energy in this case is shown in Figure 6. Using adaptive softening lengths, the change in softening lengths corresponding to the initial expansion are reflected as a secular increase in the total energy (solid line). Using the new energy conserving formalism (dashed line), this secular increase is not present and total energy is conserved to timestepping accuracy.

Static structure of a Polytrope
A simple test of self-gravitating gas dynamics is to verify the static structure of a polytrope by allowing an initial arrangement of gas to settle into hydrostatic equilibrium. In order to do so we setup ∼ 1000 SPH particles in a quasi-uniform spherical distribution and damp them into an equilibrium state using a polytropic equation of state P = Kρ γ with γ = 5/3. The low resolution is chosen in order to highlight the differences between various softening formalisms.
The exact manner in which the particles are initially setup is not particularly important, although a perfectly uniform arrangement tends to produce numerical artifacts in the collapsed particle configuration whilst a clumpy initial setup takes longer to settle to equilibrium. In this paper we use a quasi-uniform distribution achieved by placing particles initially on a uniform square lattice, cropped to ensure that r < 1 and with a small, random perturbation of amplitude 0.2∆ (where ∆ is the lattice spacing). The particle configuration is shifted slightly to ensure that the centre of mass is placed at the origin. Using a lattice spacing of ∆ = 0.15 results in a total of 1086 particles in the calculations.  . Total energy conservation during the dynamical relaxation of the perturbed 1000-particle Plummer sphere. In this case the initial velocities were multiplied by a factor of 1.2. Using adaptive softening lengths but without the new term (solid line) the change in the softening lengths caused by the initial expansion can be seen to cause a secular increase in the total energy. Adding this term (dashed line), the total energy is conserved to timestepping accuracy.
The exact solution for the polytrope static structure is computed by solving the equation numerically using a simple finite difference scheme. The solution is then scaled to give a polytrope of radius unity. In code units (discussed in §4.3) a polytrope of radius unity is obtained by choosing K = 0.4246 in P = Kρ γ . In all simulations the density and SPH smoothing length are calculated by direct summation using the iterative method described in §4.2. Also the variable smoothing length terms in the SPH equations are used throughout. In order to isolate the effects of the softening formalisms we Note that the SPH smoothing length is adaptive in all cases calculate the gravitational term by a direct summation over the particles (rather than using a treecode) with a standard second-order leapfrog scheme for time integration using a timestep controlled by a Courant condition based on the signal velocity (Monaghan 2005). The particles are damped to an equilibrium using a standard form of the SPH artificial viscosity (Monaghan 1997) together with a damping in the force equation which is independent of particle number, given by where f is the force per unit mass. Note that the polytropic equation of state means that the kinetic energy removed by the artificial viscosity and damping terms is not deposited as thermal energy but rather allowed to escape from the system. The equilibrium configuration of the γ = 5/3 polytrope with various softening formulations are shown in Figure 7 and may be compared in each case to the exact solution given by the solid line. The fractional errors (fi − fexact)/fexact are also shown in an inset plot in each panel. The top two panels show the results using fixed Plummer (top left) and cubic spline (top right) softening, where, not knowing the 'optimal' choice a priori we have used the rule-of-thumb given by Springel (2005) whereby the softening length is chosen to be ∼ 1/40 of the average particle spacing in the initial conditions. Thus guided we choose h = 0.06 for the cubic spline softening, using half of this value, h = 0.03, in the Plummer softening (see discussion in §4.3.1).
Using adaptive softening lengths without the energy conservation term (bottom left panel) shows a small improvement over the fixed softening results, mainly in the outer regions where the force estimate is much less noisy. The density resolution in the centre is slightly lower in this case, but this is substantially improved when the energy conserving term is incorporated (lower right panel). The error in the outer regions is also improved by the energyconservation term. The more compact distribution produced in this case is consistent with the additional term being always in the direction of increasing the gravitational force (see §3).

Polytrope oscillations
Having obtained the static structure, studying the radial oscillations of the polytrope provides a test of the energy conservation properties of the code. In order to do so we apply a radial velocity perturbation of vr = 0.2r to the static solutions obtained in the previous section. In order to distinguish effects due to the softening formulation from effects due to the timestepping algorithm we use a very low Courant number of Ccour = 0.05 for this test. In general, however, non-conservation effects from the softening formalism are much larger than effects due to timestepping. No artificial viscosity or damping is applied for this problem.
The evolution of the total energy of the system is shown in Figure 8 using cubic spline softening with fixed and adaptive softening lengths. Using a fixed softening length (solid line) the total energy is conserved exactly (ie. to timestepping accuracy). Adapting the softening length using the method described in §4.2 results in non-conservation of energy (dashed line). Incorporating the additional pseudopressure term into the adaptive softening formulation restores the total energy conservation (dot-dashed line).

SUMMARY
In this paper we have described an algorithm for using adaptive softening lengths in both SPH and N-body codes which retains the conservation of both momentum and energy. The formalism requires the computation of an additional gravitational force term which is similar in form to the SPH pressure force and is therefore straightforward to implement in any SPH code at almost no added cost. For pure N-body codes calculation of the additional term requires some extra work since quantities such as the density must be evaluated using the smoothing kernel. However even in this case the cost is small compared to the evaluation of the long-range gravitational forces using a treecode.
The softened gravitational force can be symmetrised either by using an average of the softening lengths or alternatively an average of the softening kernels, where the latter is preferred because of the manner in which the density is calculated. The choice of softening kernel is completely arbitrary, with calculations in this paper were made using the standard SPH cubic spline kernel (although any of the kernels proposed by Dehnen 2001 could be used).
Use of spatially variable ('adaptive') softening lengths is found to provide near-optimal softening for arbitrary mass distributions using a single, fiducial choice of the adaptive softening parameter N neigh . This contrasts to the results of Athanassoula et al. (2000) where the optimal (fixed) softening length was found to depend strongly on the number of particles and parameters such as the central concentration and shape of the mass distribution. For a mass distribution where more than one component is present, we find that the use of our adaptive softening length formalism can give more than an order of magnitude improvement in the errors on the least dense component.
The main advantage of the formalism presented here is that adaptive softening lengths can be used whilst maintaining energy conservation to timestepping accuracy. This was found to be particularly important in the case of collisionless N −body simulations where secular increases in the total energy were found to result from the use of adaptive softening lengths without the energy-conserving term. For self-gravitating SPH simulations the new formalism is a natural and self-consistent choice which is found to give a small improvement in resolution and energy conservation over traditional ad-hoc formulations for essentially zero additional cost.
q ≥ 2. (A3) Alternatively ∂φ/∂h can be evaluated from the potential and force functions according to where K(q) = hφ and K ′ (q) = h 2 φ ′ are the functional forms of the potential and force kernels.

APPENDIX B: ADAPTIVE SOFTENING LENGTH FORMALISM USING AVERAGED SOFTENING LENGTHS
A alternative way of symmetrising the gravitational potential is to use an average of the particle softening lengths. This is similar to the approach taken in the adaptive smoothing/softening length formalism used by Benz (1990). The main difference is that we symmetrise the gravitational potential rather than the force and are therefore able to account for the spatial variation of softening length in the equations of motion, leading to the conservation of both momentum and energy. Note that this is only made possible because of the self-consistent relationship between the density and the smoothing length described in §4.2.
Using an average of the softening lengths, the gravitational part of the Lagrangian can be written in the form