197 7MNRAS.181. .

375G

Mon. Not. R. astr. Soc. (1977) 181, 375-389

Smoothed particle hydrodynamics: theory and

application to non-spherical stars

R. A. Gingold and J. J. Monaghan* Institute of Astronomy,

Madingley Road. Cambridge, CB3 OHA

Received 1977 May 5, in original form February 2

Summary. A new hydrodynamic code applicable to a space of an arbitrary

number of dimensions is discussed and applied to a variety of poly tropic
stellar models. The principal feature of the method is the use of statistical
techniques to recover analytical expressions for the physical variables from a
known distribution of fluid elements. The equations of motion take the form
of Newtonian equations for particles. Starting with a non-axisymmetric distri-
bution of approximately 80 particles in three dimensions, the method is found
to reproduce the structure of uniformly rotating and magnetic polytropes to
within a few per cent. The method may be easily extended to deal with more
complicated physical models.

1 Introduction

Many of the most interesting problems in astrophysics involve systems with large departures
from spherical symmetry. This may occur either because the initial state lacks spherical
symmetry, as in the case of a protostar forming from a dense interstellar cloud, or because
non-spherical forces arising from rotation or magnetic fields, as in the case of the fission of a
rotating star, play an important part in the dynamics. Frequently these sources of non-
spherical symmetry will be found combined.
Because of the complexity of these systems numerical methods are required to follow
their evolution. However, the standard finite difference representations of the continuum
equations are of limited use, because of the very large number of grid points required to
treat each coordinate on an equal footing. If, for example, 20 points along the radial direc-
tion give adequate accuracy for a spherical polytrope, we may require (20)3 such points to
give the same accuracy for a highly distorted polytrope. This difficulty is mirrored in the
evaluation of multiple integrals.
For the astrophysical problems a numerical method which allows reasonable accuracy for
a small number of points is required. Ideally it should also be simple to program and robust.
An early attempt to provide such an alternative to the standard finite difference method was
made by Pasta & Ulam (1959). They replaced the continuous fluid by a fictitious set of
* Permanent address: Mathematics Department, Monash University, Clayton, Victoria 3168, Australia.

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

376 R. A. Gingold and J. J. Monaghan

particles with inter-particle forces designed to mimic the true pressure and other body
forces. The weakness in this method is that transport processes are difficult to include
correctly.
A better method is to make use of the Lagrangian description of fluid flow which auto-
matically focuses attention on fluid elements. In the discrete version, parcels of fluid move
according to the Newtonian equations with forces due to the pressure gradient and other
body forces: gravity, rotation and magnetic. The central feature of our analysis^ is the
method we use to determine the forces from the current positions of the fluid elements.
For fluid elements of equal mass, the number per unit volume must be proportional to
the density. In addition, unless special symmetry is introduced from the start, the positions
of the elements will be random because of the complicated motion which is inevitable for
large 7V-body systems. We therefore make the assumption that, at any time, the positions of
the fluid elements are randomly distributed according to the density. To recover the density
from the known distribution of elements is then equivalent to recovering a probability distri-
bution from a sample. Statisticians have given two methods for doing this which are well
suited to the fluid problem. The first is the smoothing kernel method (Bartlett 1963; Parzen
1962), and the second is the delta spline technique (Boneva, Kendall & Stepanov 1971). Both
methods may be thought of as an approximation to an integral determined according to the
Monte Carlo procedure. Since the Monte Carlo method is known to give reasonable estimates
of multiple integrals with fewer points than finite difference methods often require, it is
plausible to expect a reduction in work if the statistical smoothing methods are used. We call
this method smoothed particle hydrodynamics (SPH).
In this paper we first give a detailed description of the smoothing method and establish
conditions which guide the choice of the smoothing kernel. Static spherical polytropes are
then studied by relaxing from an initial non-spherical configuration with a damping term in
the equations of motion. The free non-spherical oscillations of polytropes are then examined.
Finally the departures from spherical symmetry produced by uniform rotation and magnetic
fields in polytropes are determined and compared with results from perturbation theory.

2 Recovering distributions and body forces

2.1 THE DENSITY DISTRIBUTION

The equation of motion of the /th element of fluid with volume Auy, centre of massry and
density py is

d2 tj
PyAUy— =-AUyVP+PyAUyFy, (2.1)
dt

or

d i
2ii Vp + Fy, (2.2)
dt2 Pi

where Fy is the body force acting on the element of fluid and VP is the pressure gradient at
ry. Since in our approximation the element of fluid is described dynamically by a point, we
shall call it a particle, and (2.2) the equation of motion of the /th particle.
It is convenient to begin our analysis by considering the calculation of a smoothed
*Leon Lucy has proposed and experimented with a similar method. See the acknowledgment.

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

Smoothed particle hydrodynamics 377

density from a set of points, the various r;- distributed according to the density. Following
Parzen (1962), we consider a smoothed density ps(r) defined by

Ps(r) =
J^-OpCO*', (2.3)

where IP is a function satisfying the condition

Jw(r) dr=l, (2.4)

where the integration is over all space.

If p(r ) is unknown, (2.3) cannot be evaluated, but if we have a set of N points (i*!, r2,... ,
r^v) distributed according to p, the integral can be evaluated by the Monte Carlo method
(Hammersley & Handscomb 1964). Thus, defining Pjv(r) by

M ¿V
PN(r)=-Z W(r~tj), (2.5)

where

M=jp(t)dr, (2.6)

we find, with E denoting the expectation

I f f N
E [Piv(r)] S^J...J Pjv(r) II PO-,)*,- = ps(r). (2.7)

In our numerical procedure only one sample distribution is produced each time. The
equality (2.7) is therefore to be understood as implying that if we were to create an
ensemble of models, each starting with a different array of points consistent with the initial
conditions, then the ensemble average of Pjv(r) would be ps(r).
The error involved in replacing ps(r) by Pjv(r) is ±a, where o is defined by

a^Kp^-PsOr))2]

.M2 _ , 1 \M ^ 12
(2.8)

To complete the chain of analysis it is necessary to show that a PP(r) can always be
chosen so that, as ^increases, ps(r) becomes a better approximation to p(r). We establish
this result in the next section.

2.2 CHOOSING THE KERNEL ^(r)

Intuitively it seems reasonable to expect that H^r) can be made more like 0(r) as TV becomes
larger. If this is the case then

PsW^PO*) as TV-*«*.

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

378 R. A. Gingold and J. J. Monaghan

To make this result more precise it is convenient to write h/(r) in the form

H/(r) = ¿£(r//z), (2.9)

h
where /z is a parameter with the dimensions of length and the space is assumed to be three
dimensional. By an easy generalization of a theorem due to Parzen (1962) we find that if

/z->0 asTV-*00

and if Kin) is a Borel function satisfying

AT(u)du-l, |u2£(u)|->0 as lui-*00, |/f(u)|du<co

f

where the integrals are over all space, then

PN(r)~>p(r) as/V-»°o.

It proves convenient to choose h'(r) to be an even function. Typical examples of kernel

functions in three dimensions are

/ 1 f2 3rt(l-|rl//Q r..,S(\r\/h)
exp(-r2/h2) (ii) (2.10)
3 M
®fe) 4 nh ~¡r

where H is the Heaviside step function and S is the spherical delta spline discussed in
Appendix 1. Each of the functions in (2.10) is a member of a sequence of functions which
represents the delta function.
In addition to requiring pN{t) -> p(r) we require that o should be as small as possible. To
satisfy these conditions we choose h by minimizing the functional

Z,(r) = {^[p^(r)]-p(r)}2 + ^[(pw(r)-ps(r))2]

(2.11)
= £’[(pA)(r)-p(r))2].

Using (2.5) we find

¿(r) =~j^2(r - r') p(r')dx' + (l -d)[JlV(r - r') p(r') *'] '

+
P2(r) — 2p(r) JJV(r — r') p(r') dr . (2.12)

Since IV(r-r') is strongly peaked at r = r', we can expand p(r') about r. Keeping only the
dominant terms, we find

M f „ , Í V2p r )2
2 2
Z,(r)~ —p(r)jlk (r)c?r +|—JfV(r')r' (ir'| . (2.13)

Using (2.9) the minimum of L(r) is found to occur at

27 Mp
(2.14)
N (V2p)2

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

Smoothed particle hydrodynamics 379

Since p is unknown (2.14) cannot be used directly except to infer

/z oc l/Af^andimin oc l/Af4/7.

We find an appropriate choice of h to be given by

h = b{(t2)-{x)2Ÿl\ (2.15)

where b is adjustable and

Ipcdr j
<c(r)> = = — Zcir,-).
M Nj’

In the problems we consider, the derivative of the smoothed function is required to be

continuous. For this reason the second of the kernel functions in (2.10) is not useful. Of the
other possible kernel functions we have concentrated on the Gaussian and spline functions.
To decide between the Gaussian and the spline kernels we took a known distribution and
distributed a set of points. The goodness of fit [p;v(r) -p(r)]2 was then evaluated for various
values of b. For forty points there is negligible difference between the two kernels, but for
80 points the Gaussian was much more accurate. For this reason we prefer the Gaussian
kernel and use if for the results reported here.
For our stellar models we choose b by requiring the smoothed particle model to fit the
known density of the spherically symmetric hydrostatic model. More elaborate procedures
could be used but we have found those described to be successful for the models considered.

2.3 THE GRAVITATIONAL POTENTIAL

We use the gravitational potential 0 defined by

ßN(r')cti
(2.16)
|r —r'|

Using (2.5)

GM N rW(r’ - r^dr
0= (2.17)
N pW ¡r-r'|

rW(r'-rj)dr'
(2.18)
//=
J Ir-r'| ’

can be evaluated easily noting that

V2/y = — 4 7T — iy). (2.19)

We find

GM N 47t ru¡
v0 = ~— E W(u) u2du)Vuj, (2.20)
«/Jo

where
u = r
/ -r/

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

380 R. A. Gingold and J. J. Monaghan

For the Gaussian W defined by (2.10) (i) with/= l/h2

GM N 2 /A1/2r . , 1 n 1
2
N r L exp (-/«/)--
W/J o exp (—fu ) du\J Vu¡ (2.21)

and

// = “(“) exp (—fu2) du.

Uj\TÏ/ J0

The equivalent formulae for the delta spline W involve polynomials, and are easier to
evaluate.

2.4 GENERAL DISTRIBUTIONS

To find the smoothed version of any other scalar (or vector) field A(r), we define the
smoothed field >ls(r) by

A W{t-r)A{x)dx\ (2.22)

where, in general, the kernel differs from that in (2.3). Then an estimate of ^s(r) is

M N A{ri)
^N(r)=- I ^(r-iy) (2.23)
N j~ i p(iy)

The error in this estimate is ± a where now

^20/)
W\x r
/) (2.24)
P2(r/)

The approximations involved become better when ^4(r) is distributed similarly to the
density. This is the case for temperature and entropy, but for the magnetic field it is not in
general true. To deal with this case importance sampling is useful and we discuss its applica-
tion in the next subsection. Where the field has known symmetry properties antithetic
variables can be used to improve the accuracy.

2.5 THE MAGNETIC FIELD AND CURRENT

According to the prescription given in Section 2.4 an estimate of the magnetic field is given
by

M N B,
B
Ar(r)=- E ^(r-r/)-, (2.25)
N
i~ i Pi
and an estimate of the current by

, M T- B.
=
•Mr) eoC — E Wx —. (2.26)
N j Pi

However it is usually the case that there is field inside and outside the star and (2.25) is
then a poor approximation to Bs(r), and (2.26) is an even poorer approximation to Js(r). To

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

Smoothed particle hydrodynamics 381

improve the approximation we use importance sampling (Hammersley & Handscomb 1964)
in the form

Bs(r) = B0(r)+jH/(r-r') [B(r')-B0(r')]íír', (2.27)

where B0(r) is an approximation to B(r). To obtain this approximation we solve

e<)C2Vx Bo =-Mr), (2-28)

in the form

JatO-)x (r — r)dt
Bn=B ext i I (2.29)
FnC2 JÍ |r-r/13

where Bext is any superimposed external field. It could, for example, be the field permeating
an interstellar cloud from which a star is forming. Substituting (2.26) into (2.29) we find

M " r/B, n^ 4n
Bo(r) Bext - B, lV(r - r/) (2.30)
47TJV ^ llpy ’ 9r/ dr,- ]■

where Ij is defined by (2.18).

For the Gaussian kernel (2.10)(i) with/= l//z2, the approximate field becomes

where u = Uy = r — iy.
The field we use is

M W(r — r)
®iv(r)= ®o(r) +"77 S {®(r/) —
(2.32)
N
j Pi
and the current is obtained from the curl of (2.32). Thus

ec\C2M W
r = x —
J/v( ) Tl X {2 BO*/) Bo O'/)} + JextC1*)? (2.33)
N
i Pi
where Jext is the current associated with Bext. This procedure, as we show later, gives a
good fit to the current and the field.

3 Equations of motion

The equations of motion of the fluid particles for a uniformly rotating polytrope of index n,
with an internal magnetic field B are

d2tf dtj ^l+fl) — dtj JxB

—' = -r—' -Kp¡in~l -—Vp V<t>-fix(ßx r) -2fix — + 1,2,...,7V
2 dt
dt dt n
(3.1)

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

382 R. A. Gingold and J. J. Monaghan

where the pressure P = Kßl + Yln, SI is the angular velocity, B the magnetic field and J the
current. The damping term Tdtj/dt has been introduced to allow static models to be
calculated. For the rotating models considered here only the static structure is required,
and the Coriolis term can be dropped. Using the dimensionless variables Xy, D, t, b defined
by

p = XD, Tj = axj, t = ßr, B = 5b (3.2)

where

j (n + l)K\1ln-1 a2n
(3.3)
4nG ’ KX^il+n)’

ftnd X and B can be chosen for convenience. (3.1) becomes, on dropping the Coriolis term,

d2 Xj yd\ • n (Vxb)xb
' - D)ln~l VD;— — - w2i x (i x x.) + v (3.4)
2
dr dr 4n ’

where 7 is a new damping constant

e0c2nB2
ß= co 1, (3.5)
(=?)' 4TrGX2a.2’
and is the scaled gravitational potential where

I <t>Q
3
GXa GM
and we have chosen the value of ß = /Z>(x)dx such that were the representations of integrals
by sums in Section 2 to be exact, then D(0) would equal unity and X would be the central
density. Our scaled variables are therefore similar to the usual poly tropic variables
(Chandrasekhar 1939). However, since in our models D(0) + 1, our length scale is related to
the poly tropic variable I by £= £)(0)(” ~ |x |.
For the models considered here the magnetic field variation is calculated in the flux
freezing approximation

(3.6)

where d/dt is a derivative following the motion. To integrate this equation forward we
replace u by the smoothed velocity field. Equation (3.6) has the advantage that it automatic-
ally generates the quantity B/p required at each fluid element to produce the smoothed
field.

4 Numerical tests - spherical models

4.1 CONSTRUCTION OF STATIC MODELS

To construct a static model we follow the damped motion of a set of particles from some
initial distribution of position and velocity until the system comes to rest. Typically the
particles were initially at rest, distributed in space either according to a random Gaussian
distribution or alternatively on a spherically symmetric cubic lattice. In the former case the

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

Smoothed particle hydrodynamics 383

Figure 1. The central density D(0) as a function of time r for two damped hydrodynamic sequences.
The initial configurations are given in the text.

initial coordinates of the particles were adjusted so that the centre of mass was at the centre
of the coordinate system. As a check, the position and velocity of the centre of mass were
monitored throughout the calculations.
The approach to equilibrium for two initial configurations with different degrees of
damping is illustrated in Fig. 1 for a polytrope of index 1. The solid line represents the
behaviour of £>(0) as a function of the scaled time r, in a sequence that commences with 33
particles on a cubic lattice and 7 = 0.05. The broken curve shows a sequence, with 7 = 0.15,
commencing with 40 particles distributed normally about the origin.
The models finally obtained are found to be nearly independent of both the damping,
the initial configuration, and the number of particles. These model sequences commence, of
course, with a good deal of spherical symmetry. Quite irregular initial distributions can also
be successfully treated. Fig. 2(a) shows the density profile in the (A, 7) plane of an initially
non-spherically symmetric distribution which leads to the symmetric distribution of Fig.
2(b) representing a polytrope of index 1.

4.2 STATIC STRUCTURE

Polytropes of index 1 and 1.5, constructed using about 40 particles and a wide range of
smoothing constant b (defined by equation 2.15) were found to have density profiles which
matched the true density to within a few per cent over the bulk of the star. The density
profile in the outermost 10 per cent of the polytropic radius for the polytrope of index 1,
however, reflects the nature of the smoothing function rather more than it does the actual
distribution of matter. The size of this region can be decreased by employing a larger
number of particles.
Sequences that commenced with a high degree of spherical symmetry yielded similarly
highly symmetric polytropes. The departure from spherical symmetry in other cases was less
than 2 per cent.
Polytropes of index 2.5 were also constructed employing both a range in the number of
particles and in the value of the smoothing constant b. In Fig. 3 we display the appropriately
scaled density profiles of two such models. The curve represents the true density for n - 2.5
while the filled circles show the density profile of a model constructed with7V=80. Since
this model is highly symmetrical the density profile along only one axis is shown. Also
13

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

384 R. A. Gingold and J. J. Monaghan

Figure 2. An example of initial and final smoothed density in the xvy plane for a polytrope of index 1.
The initial state was a superposition of two Gaussian density distributions.

Figure 3. Density profiles for a polytrope of index 2.5. The Emden density is shown: ; the
80-particle SPH is shown: • • • •. The variation in density for a given X along the x,y,z axes is indicated
by the size of the filled circle. The analogous variation for 40 particles is shown by the bar.

indicated in Fig. 3 is the less symmetrical density profile of a model constructed with one
half as many particles. The range in density at points on each of the three coordinate axes
for this model is indicated by the vertical bars. The density profiles in this case fit the true
density more closely than might appear from the figure. Along each of the coordinate axes

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

Smoothed particle hydrodynamics 385

Table 1. Parameters of polytropes of index 2.5.

N 40 40 40 40 40 80 200 200 200

lib2 0.5 1.0 1.5 2.0 2.5 2.0 1.6 2.0 3.0
l/h2 0.35 0.60 1.05 1.69 2.81 1.01 0.67 0.79 1.19
D{0) 0.72 1.22 1.97 2.97 4.92 1.52 1.06 1.08 1.34
<£2> 1.17 1.88 2.15 2.27 2.31 2.55 2.47 2.65 3.01

the density profile is quite good, but the peak value is offset from the origin. Nevertheless,
the improvement achieved by increasing the number of particles from 40 to 80 is remarkable
and surpasses the yjN improvement we would expect in Monte Carlo integrations. This is
probably due to 40 particles being intrinsically too few.
Some brief details of various models for n—2.5 are given in Table 1. In each case the
sequence commenced with a Gaussian distribution of particles about the origin. Although
the final values of Z)(0) vary with both the smoothing parameter and the number of
particles, in each case the density profile after dividing by D(0) is similar. Also displayed in
Table 1 are the mean squared radial position of the particles given by

These values are smaller than the value of 4.8 obtained by performing the integrations using
the density profile in the above expression. This discrepancy is not surprising since the
integrand p£4 has a sharp maximum beyond the position of the bulk of our particles.

4.3 UNDAMPED OSCILLATIONS

Several hydrodynamic sequences were followed with damping excluded. These were found
to oscillate in a mixture of modes reflecting the initial state of the model. In each case a
dominant period of oscillation of the central density was manifest. This matched the periods
of oscillation of polytropes of index n in the range 1-2.5 given by Kopal (1938) to within
10 per cent for 7V~40. This error can be reduced by using more particles. During extended
runs over many cycles of large amplitude oscillations (ôZ)(0)/Z)(0)~ 0.3) the total energy#
of systems with 40 was found to oscillate with |ô#/#|^0.1. This error is consistent with
replacing integrals by sums according to the Monte Carlo procedure.

5 Numerical tests — non-spherical models

5.1 UNIFORMLY ROTATING POLYTROPES

Uniformly rotating polytropes were studied to determine the accuracy with which the
technique reproduced a non-spherical structure.
In Fig. 4 the polar and equatorial density profiles are shown for a rapidly rotating poly-
trope of index 1.5 for which co2 = 0.024. The figure also shows the density profiles for the
same model obtained using the approximation technique of Monaghan & Roxburgh (1965).
It is clear that the agreement is good. All models were found to be symmetric about the
rotation axis and the equator to within 5 per cent.
Because our models do not have D(0)= 1, the parameter a of Monaghan & Roxburgh is
related to our co2 by a = 2cj2/nD(0). The model shown is therefore on the verge of breakup.
Since our method does not produce fluid particles near the edge, the critical go corre-

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

386 R. A. Gingold and J. J. Monaghan

Figure 4. The density profiles for a uniformly rotating polytrope of index 1.5. The SPH results are shown
thus: polar density: lower curve; equatorial density: upper curve. Perturbation analysis (Monaghan &
Roxburgh 1965) shown by and a a a a .

sponding to breakup cannot be determined accurately. Of course, with more particles, and
therefore a smaller /z, the critical co can be calculated as accurately as desired. Alternatively
test particles could be introduced.

5.2 MAGNETIC POLYTROPES

The static structure of polytropes with both poloidal and toroidal fields was studied by
starting with a static, non-rotating, polytrope and then superimposing the field. The poly-
trope was then allowed to relax to a static structure. Because the main purpose of this study
was to explore the numerical method we chose initial fields which were known solutions of

Figure 5. Poloidal magnetic field and current in a polytrope of index 1. Perturbation analysis (Monaghan
1965) shown thus: . The smoothed initial field and current shown thus: The final static
field and current shown thus: .

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

Smoo thed particle hydrodynamics 387

Figure 6. Density profiles for a polytrope of index 1 with a dipole poloidal field. Perturbation analysis
(Monaghan 1965) shown thus: polar density: • • • •; equatorial density: ■ ■ ■. The SPH density shown
thus: . The full field smoothing method has been used.

the first-order perturbation equations. The poloidal field was taken from Monaghan (1965)
and the toroidal field from Roxburgh (1966).
In Fig. 5 we show the initial field and current on the x axis calculated from the analytical
expression for a dipole field in a polytrope of index 1. Also shown is the initial and final
smoothed field and current calculated according to the procedure of Section 2.5. The agree-
ment between the initial field and its smoothed equivalent is very good. We believe it could
be further improved by adjusting the smoothing parameter or by adopting a different
value of this parameter for each component of the magnetic field.
The analytical equilibrium field is based on a first-order perturbation analysis which
assumes the field can be constructed from a non-perturbed density. Since we find density
perturbations of ~10 per cent, we expect the final field to differ from the first-order pertur-
bation results by quantities of this order. The difference between the initial and final field
and current shown in Fig. 5 is therefore not unexpected.
In Fig. 6 the equatorial and polar density profiles are shown for both the present
numerical calculations, and for the first-order perturbation results. Since our models have
£>(0) # 1, and a field and current which differ from the analytical one by approximately a
scale factor, the relation between r? and the factors to and k of Monaghan (1965) is approxi-
mately
V /computed #2(0)\ 2
nk2D{0)lJtlfn \ analytical5Z(0)/
The agreement between the first-order perturbation results and our numerical results is
very good. Small changes, of the order of 10 per cent of the deviations from the unperturbed
density, are to be expected because of errors in the perturbation method, but this has a
negligible effect on the density profiles.
The toroidal field investigated is zero outside the polytrope and the smoothed field can
be obtained satisfactorily without using the importance sampling device of Section 2.5.
With 40 particles in a polytrope of index 1 the initial smoothed field reproduced the
analytical field to within < 5 per cent.
During the calculations the constancy of the magnetic flux was monitored and found to
remain constant to within 2 per cent.

6 Computational requirements
All of the sequences described in this paper were stored in 60—90K bytes of core storage in
an IBM 360/165. A typical 7V=40 sequence with no magnetic field requires about 0.25 s per

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

388 R. A. Gingold and J. J. Monaghan

time step and about 200 time steps (50s) per model. The time step is fixed by calculating
the minimum of h/vm^ (^/^max)172 and h/D1'2” where umax is the maximum particle
velocity, Fmax the maximum force and the last criterion is based on the speed of sound. The
minimum is then multiplied by a constant. If this constant is in the range 0.1—0.5 the
integration is stable without requiring excessive time. The time can be halved by storing the
current forces and densities at the particle positions. We are currently performing calculations
for dynamical sequences leading to fission. These require approximately four minutes of
computing time for an 80-particle configuration.

Conclusions
The results of this study show that the smoothed particle method is a simple technique
which gives satisfactory results for oscillating polytropes, and for polytropes which relax
from a non-spherical initial state to a spherical final state. Rotation and magnetic fields may
be included without difficulty, and the comparison with the perturbation results shows that
moderate distortion can be reproduced accurately. Structure on a finer scale or greater
accuracy can always be obtained by increasing the number of particles and by using the
devices known to improve Monte Carlo integration methods.

Acknowledgment
In a lecture give at the Institute of Astronomy, Cambridge in 1976 Leon Lucy discussed the
use of smoothing techniques for hydrodynamic codes. His ideas were adumbrated to us by
our colleagues, but the mathematical development in this paper is independent of his work.

References
Bartlett, M. S., 1963. Sankhya (A), 25, 245.
Boneva, L. I., Kendall, D. & Stepanov, I., 1971. J. R. stat. Soc., 33,1.
Chandrasekhar, S., 1939. An introduction to the study of stellar structure, p.87. University of Chicago
Press.
Hammersley, J. M. & Handscomb, D. C., 1964. Monte Carlo methods, Methuen, London.
Kopal, Z., 1938. Mon. Not. R. astr. Soc., 99,33.
Monaghan, J. J., 1965. Mon. Not. R. astr. Soc., 131,105.
Monaghan, J. J. & Roxburgh, I. W., [Link]. Not. R. astr. Soc., 131,13.
?aiZQn,E., 1962. Ann. Math. Stat., 33, 1065.
Pasta, J. R. & Ulam, S., 1959. Mathematical tables and other aids to computation, 13,1.
Roxburgh, I. W., [Link]. Not. R. astr. Soc., 132, 347.

Appendix
In one dimension the simplest representation of a sample is by a histogram. To smooth the
histogram the constraints of minimizing the slope, while retaining reproducibility of the
data, can be used. The resulting smoothing function is the delta spline of Boneva, Kendall
& Stepanov (1971).
In three dimensions there are various possible generalizations. Our experiments have been
based on the following.
Around a sample point construct the unit ball, i.e. the sphere of unit radius. This is one
generalization of the unit histogram. Surround the ball by concentric shells of radius 77=/.
Now construct the spherical delta spline S(r) by the rules
r n+i
Min 47t J j r2 dr with 47r Sr2dr = 8ft 1,2,3.
J r¿

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System

197 7MNRAS.181. .375G

Smoothed particle hydrodynamics 389

These rules ensure minimization of the slope with the constraint that the integral over all
space is equal to the contribution from within the unit ball. The resulting set of equations is
easily solved and the spherical delta spline is found to oscillate with an exponentially
decreasing amplitude.
An alternative we haven’t experimented with is based on the subdivision of the space into
cubes. Then the function to be minimized is (33iS/3x3>,3z)2 and the delta spline becomes just
a product of one-dimensional delta splines.

© Royal Astronomical Society • Provided by the NASA Astrophysics Data System