Subsections


Smoothed Particle Hydrodynamics

Smoothed particle hydrodynamics (SPH) is an N-body integration scheme that was introduced independently by Gingold and Monaghan (1977) and Lucy (1977) in order to model complicated non-symmetric phenomena in astrophysics. For this reason, it is particularly well suited for modelling the Capture Theory. This chapter is a presentation of the existing theory, along with several new extentions that were developed during the course of this investigation2.1.


Motivation

SPH has several advantages over other methods, that originate directly from its discrete particle nature. For example, in Figure 2.1 two centrally condensed spheres of gas are being modelled. The concentric circles represent contours of equally spaced constant density. Two dimensions are shown for clarity, but the following discussion is equally valid in three.

Figure 2.1: Comparison of simulation methods. (a) Grid based, (b) SPH.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/grida3.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 2.1(a) uses a traditional finite difference grid, whilst Figure 2.1(b) uses SPH. There are three important features to note. First, the size of a cell in the grid based method must be such that the central condensation is resolved adequately. This leads to a small cell size everywhere, even in regions near the surface of the sphere. SPH handles this contrast in density by distributing a higher concentration of particles in the central regions than at the surface. This leads directly to the second point:

In SPH, the simulated region is represented by particles which carry a certain mass2.2, thus computational run-time is not wasted following the properties of regions devoid of matter. In a simple grid based method, empty regions of space would exist in the memory of the computer, and their evolution in time would be calculated also, leading to longer simulations, or lower spatial resolution (achieved by reducing the fineness of the grid).

The final main difference is in the way the system evolves in time. If in this example the material is flowing towards the centre of each sphere, then two different concepts arise. In the grid based method the material flows through the cells, resulting in an Eulerian system, whereas in SPH the particles move, making it Lagrangian in nature. The advantage of a Lagrangian system in this example, is that as the density in the centre increases, more particles move into this region, resulting in the spatial resolution there increasing. In the grid based method, the density in the central cells of each sphere increases as the mass flows in, but the resolution in that region remains the same. Sophisticated grid based methods exist which subdivide grid cells containing more than a certain mass (Couchman, 1991; Couchman et al., 1995), thus increasing the resolution in that region, but when to divide, and perhaps subdivide again, is almost completely arbitrary. To summarize, SPH handles density gradients and voids quite naturally, and also the resolution in a system is essentially automatic.

So far, the picture built up of SPH is one of point mass particles, distributed according to the density of the material. The `Smoothed' part of SPH is now introduced, according to the desire to represent the particles as a continuum, rather than individual points in space. This is necessary, as the density, which is a continuum variable, is required in order that the hydrodynamic forces on a particle may be calculated, as shall be seen later. In Figure 2.1(a), the density at the centre of a grid cell can easily be calculated from knowing the volume of the grid cell, and the mass associated with it. The density at any point within a grid cell can then be interpolated from surrounding cells via a simple weighting function. In SPH, the procedure is not quite so straightforward, as in general a regular structure with which to associate a volume to a particular particle, is not present. However, at the heart of SPH exists an interpolating function, which is used to estimate the physical quantities and their derivatives at any point in space from disordered particles.


The Kernel and the Smoothing Length

The purpose of the kernel is to transform the representation of a particle from a point mass one, to a spatially spread one. More specifically, the kernel is a bell shaped function which has a smooth and continuous first (and usually second) derivative, and is smooth and continuous itself. A simple mathematical function which meets these requirements is the Gaussian. This function is used to spread the mass (and other properties) of the particles in space, and is normalised so as to conserve the physical quantities in the system. This results in Equation 2.1:

\begin{displaymath} W_{G}(r,h)={\frac{1}{h^{\nu}\pi^{\nu/2}}} \exp\Big(-\frac{r^{2}}{h^{2}}\Big) \end{displaymath} (2.1)

where $\nu$ is the number of dimensions, $r$ is the distance from the particle at which the quantity is being estimated, and $h$ is a scaling variable with dimensions of length. The scaling variable is referred to as the smoothing length, as it controls the degree to which the particle is spread in space. This normalised smoothing function is referred to as the kernel, and is denoted by $W$. The subscript $G$ on the kernel in Equation 2.1 indicates that it is the Gaussian. Given a distribution of particles, the kernel can be used to calculate the density at any intermediate position in space by using:
\begin{displaymath} \rho_{i}=\sum_{j=1}^{N}m_{j}W(r_{ij},h_{j}) \end{displaymath} (2.2)

where $r_{ij}$ is the distance between particle $j$ and the point in space being considered (position $i$), $m_{j}$ is the mass of particle $j$, $h_{j}$ is the smoothing length of particle $j$, and the sum is over all the particles present ($N$). As an example of what this may look like, consider the one dimensional scenarios presented in Figure 2.2.

Figure 2.2: One dimensional density profiles. Unit particle mass and spacing, variable smoothing length. (a) h=0.5, (b) h=1.0, (c) h=2.0.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/kernelsa2.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 2.2 shows three lines of equally spaced particles, all with the same unit mass and spacing. The theoretical corresponding continuum density profile for this distribution is a flat line of $\rho=1$. The density profiles shown are sampled from the particles using Equation 2.2.

Each particle in Figure 2.2(a) has a relatively small smoothing length associated with it, and thus the resulting density profile is significantly non-uniform. The other extreme, each particle having a relatively large smoothing length, is shown in Figure 2.2(c). In this figure, each particle has been smoothed over a large distance, resulting in an effective decrease in the spatial resolution present. If the resolution was actually required to be lowered, then in terms of computational expense, a fewer number of particles with a greater mass and interparticle distance, should have been employed to sample the region, and a smoothing length giving a density profile similar to Figure 2.2(b) should have been used. The above argument for the size of the smoothing length also applies to the force calculation, introduced in Section 2.4, i.e. the forces on a particle are undersampled in Figure 2.2(a), and oversampled in Figure 2.2(c).

Although the Gaussian kernel is simple, and interpolates with high accuracy, it has the practical disadvantage that it is not finite in extent. Thus all particles in a system contribute to the summation in Equation 2.2, even the ones that are at such a distance that their contribution is effectively zero. Usually when this kernel is utilised, the contribution is assumed to be zero when the distance from the particle is greater than $3h$. Nevertheless, a lot of particles still have to contribute to the calculation for the sampling to be adequate, which leads to an increased computational effort. To overcome this problem, many authors use a kernel which is more spatially compact, and finite in extent. A popular one, which is used in this work, is the spherically symmetric $\beta$ spline kernel proposed by Monaghan and Lattanzio (1985):

\begin{displaymath} W_{\beta}(r,h)=\frac{\sigma}{h^{\nu}} \left\{\begin{array}{c... ...\leq 2\\ 0&\textrm{\hspace{22pt}$u\geq 2$} \end{array}\right. \end{displaymath} (2.3)

where $u=r/h$, and $\sigma$ is a normalisation constant with the values $2/3$, $10/7\pi$, $1/\pi$, in one, two, and three dimensions respectively. This kernel has the practical advantage that it is exactly zero after two smoothing lengths.

To compare the two smoothing kernels they are equated at the origin, so that $h_G$ can be expressed in terms of $h_\beta$ (or visa versa), and renormalised for convenience. In three dimensions this requires that:

\begin{displaymath} \frac{1}{\pi^{\frac{3}{2}}h_{G}^3}=\frac{1}{\pi h_{\beta}^3} \end{displaymath} (2.4)

leading to:
\begin{displaymath} h_{G}=\pi^{-\frac{1}{6}}h_{\beta} \end{displaymath} (2.5)

The resulting graphs, along with first derivatives with respect to distance, are shown in Figure 2.3. Although this is a three dimensional calculation, the graphs are exactly the same in one and two dimensions, as it is the shape of the kernel that is being compared. As can be seen, the kernels are very similar, the $\beta$ spline being slightly more compact, as $(b)>(a)$ in the inner region. From hereafter when referring to the kernel, unless stated otherwise, the $\beta$ spline is inferred.

Figure 2.3: Comparison of the kernels. (a) $W_{G}$, (b) $W_{\beta}$, (c) $\mathbf{\nabla}W_{G}$, (d) $\mathbf{\nabla}W_{\beta}$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/Comp02.eps,angle=0,width=1.0\linewidth}} \end{figure}


Variable Smoothing Lengths

Up until now the smoothing lengths have all been assumed to be the same, and unchanging in time. This was the case during the pioneering times of SPH, where low numbers of particles, and small density gradients were all that could be handled. As computational power increased, it was realised that the resolution in a simulation could be improved by using individual particle smoothing lengths. How exactly to specify, and subsequently evolve the smoothing length of a particle, has been achieved in various ways. As an example of one way, consider the spheres of Figure 2.1(b) again. It is clear from Figure 2.2, that if the smoothing lengths are not reduced in size during the inflow of material, the potential for increasing the resolution in the central region is lost, due to the large smearing present. In order to reduce this effect, it is reasonable to suggest that the smoothing lengths should be inversely proportional to the cube root of the local density (in $3D$):
\begin{displaymath} h=h_{0}\bigg(\frac{\rho_0}{\rho}\bigg)^\frac{1}{3} \end{displaymath} (2.6)

But as shown earlier, Equation 2.2, the density is a function of the smoothing length, resulting in an inconsistency. One way of avoiding this inconsistency is to use Equation 2.2 to calculate the density, then use Equation 2.6 to update the smoothing lengths, and repeat the two calculations until the smoothing length or the density of a particle is unchanging between iterations. This can be a very slow procedure requiring many iterations, and is hardly ever put into practice in this form.

Constant Number of Neighbours

A more common procedure for specifying the smoothing length, is to require that the kernel contains a constant number of particles, i.e. $h$ is specified such that the kernel just overlaps onto the $N_{N}^{th}$ neighbour2.3, the distance to this particle being $2h$ for the $\beta$ spline kernel. Typical values for $N_{N}$ in one, two, and three dimensions are $5$, $15$, and $55$ respectively. Most implementations that use this method are quite relaxed about including exactly $N_{N}$ neighbours, claiming that to do so would require too much computational run-time. In these cases, the number of neighbours of a particle typically varies between $N_{N}/2$ and $2N_{N}$ during a simulation.

However, when a constant number of neighbours is not maintained, an additional time dependence is introduced into the equations of motion: Consider for the moment, forces that are a function of position and smoothing length only2.4. If a simulation has exactly the same distribution of particles at a time $t_{2}$ as it had at an earlier time $t_{1}$, then in theory the forces between the particles should also be exactly the same. In a simulation where the kernel extends out to the $N_{N}^{th}$ neighbour at all times, this will be the case. However, when this criteria is not met, the smoothing lengths at times $t_{1}$ and $t_{2}$ will in general differ, resulting in different interparticle forces also. Oscillatory motion is an example where the effects of this can be demonstrated. Here an effective lag on the size of the smoothing length, can produce either erroneous damping, or exponential growth in the velocity, depending on the exact implementation of the smoothing length update. When using Equation 2.6 to update the smoothing lengths, lag, and consequent damping of the velocity field can occur if insufficient iterations are performed.

The extent of this effect has not been investigated, and perhaps is negligible in most SPH simulations. Nevertheless it is present, and to avoid complicating analysis, an efficient method of finding the $N_{N}^{th}$ neighbour of a particle was developed:

Neighbour Search Algorithm

  1. To find the smoothing length of particle $i$, first collect all particles that are within $2ah_{i}$ of $i$, where $a\approx1.2$, and $h_{i}$ is the previous calculation of the smoothing length. If the number of particles collected is less than $N_{N}$, increase $a$ and repeat.
  2. It is now tempting to just pass the collected particles to a quicksort function, to sort in order of distance from $i$, and set $h_{i}$ to be half the distance to the $N_{N}^{th}$ particle. If this procedure was adopted, a great deal of redundant information would be calculated, reducing the overall efficiency of the method. Ideally, only the position of the $N_{N}^{th}$ particle is required, the order of all other particles being unimportant.
  3. The algorithm for achieving this is essentially the same as the quicksort algorithm, but with the difference that only one particle is sorted into order, i.e. the yet unknown $N_{N}^{th}$ neighbour of particle $i$ is inserted into the $N_{N}^{th}$ position in the list.
  4. In a quicksort algorithm, a pivot is used to split the list into two. The method is most efficient when the pivot is chosen such that the list is split in half. A new pivot is chosen for each list, and both are split into smaller lists. This continues recursively, until each list contains only one element, at which time the list is sorted.

    Figure 2.4: Determining the distance to the $N_{N}^{th}$ particle. For this example, P is the value of the pivot and $N_{N}=5$. (a) Traditional quicksort, (b) Reduced quicksort.
    \begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/qsorta1.eps,angle=0,width=0.85\linewidth}} \end{figure}

  5. The method is exactly the same here, except that following the division into two lists, only the list which contains the $N_{N}^{th}$ particle is considered. Knowing which list this resides in is achieved simply by keeping track of the number of particles inserted into each side. Here the pivot is taken to be the average of the values in a list. In practice, step 1 delivers a list of the distance squares from $i$, and to avoid using the rather slow square root function, these are used throughout the algorithm.
A comparison of the two schemes is shown in Figure 2.4. In practice, approximately $10$ times more particles are collected than shown here, so that the appeal of the new reduced quicksort algorithm is much greater. The computational run-time of the algorithm is limited by the time required to locate the particles before sorting. In an efficient data structure, similar to that used in Chapter 3, this is not a problem, the smoothing length determination taking typically $10\%$ of the run-time in the force calculation routine. This method also performs just as efficiently for variable mass particles, where the smoothing length is set to enclose a constant mass, rather than a constant number of particles. This involves keeping track of the total mass in each list, rather than the number of particles. See Appendix A.3 for the pseudo-code of the above algorithm in the context of a tree-code and constant enclosed mass.


Forces

In most computational techniques there are four types of forces: Attractive, repulsive, damping, and external forces. The first three forces can be considered to act within the system, such that when they are calculated, every piece of matter interacts with every other piece of matter in the system. If the system is split up into $N$ discrete bodies, as in SPH, this leads to a computational run-time proportional to $N^{2}$. External forces act on each part of the system individually, leading to run-times proportional to $N_{E}N$, where $N_{E}$ is the number of external forces.

The combined forces act on the matter of the system, and if not in balance cause it to accelerate according to Newton's second law of motion:

\begin{displaymath} \mathbf{F}=m\frac{d\mathbf{v}}{dt} \end{displaymath} (2.7)

where $\mathbf{F}$ is the force, $\frac{d\mathbf{v}}{dt}$ is the acceleration, $m$ is the mass of the material in question, and the convention is that bold fonts denote vector quantities. In SPH, the forces are centred on the particles, such that $m$ is the mass of a particle, and the associated acceleration of the particle is along the direction of the combined forces. One form of expressing this acceleration for the forces involved in this work, is via the momentum equation shown below for particle $i$:
\begin{displaymath} \frac{d\mathbf{v}_{i}}{dt}= -\frac{\mathbf{\nabla}P_{i}}{\rho_{i}} + \mathbf{a}_{i}^{visc} - \mathbf{\nabla}\phi_{i} \end{displaymath} (2.8)

The first and second terms are the contributions to the acceleration from the hydrodynamical gas pressure force (repulsive) and gas viscosity (damping), where $P_{i}$ and $\rho_{i}$ are the pressure and density at particle $i$ respectively, and $\mathbf{a}_{i}^{visc}$ is the acceleration due to the viscosity. The third term involves the gradient of the gravitational field at $i$, and in this work is a mixture of internal and external forces, as the material is self-gravitating, and is also influenced by external gravitational sources. Each of the three terms in Equation 2.8 is now examined.


Hydrodynamical Forces

The hydrodynamical component of the acceleration is easily derived by considering Figure 2.5.

Figure 2.5: Deriving the hydrodynamical acceleration of particle $i$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/Pboxa1.eps,angle=0,width=0.7\linewidth}} \end{figure}

In this figure, a small box of mass belonging to particle $i$ is shown. This is not necessarily all of the mass of particle $i$, but some small proportion centred on it. The lengths of each side of the box are denoted by $a$, $b$ and $c$. The acceleration along the x-axis is found by considering the associated pressures in the positive and negative directions ($P_{1}$ and $P_{2}$) of the surrounding gas, along the same axis. By expressing the force and mass in Newton's second law as:
\begin{displaymath} (F_{x})_{i}=bc(P_{1}-P_{2}) \end{displaymath} (2.9)


\begin{displaymath} m_{i}=abc\rho_{i} \end{displaymath} (2.10)

the acceleration along the $x$ axis at $i$ becomes:
\begin{displaymath} \frac{d(v_{x})_{i}}{dt}=\frac{1}{\rho_{i}}\frac{P_{1}-P_{2}}{a} \end{displaymath} (2.11)

and in the limit of a vanishingly small box:
\begin{displaymath} \nabla_{x}P_{i}=\frac{P_{2}-P_{1}}{a} \end{displaymath} (2.12)

Extending this to three dimensions introduces vector quantities, resulting in the required equation:
\begin{displaymath} \frac{d\mathbf{v}_{i}}{dt}= -\frac{\mathbf{\nabla}P_{i}}{\rho_{i}} \end{displaymath} (2.13)

This is the standard hydrodynamical continuum equation for compressible fluid flows.

How to calculate the density at a point has already been discussed in Section 2.2, both for SPH and grid based methods, but nothing has been said so far about how to calculate derivatives. In a grid based method, the acceleration of the material in a grid cell could be calculated by simply applying the approximation of Equation 2.11, where $P_{1}$ and $P_{2}$ would be the pressure from neighbouring grid cells of cell $i$. In this grid based approximation to the pressure gradient, more accuracy and stability is obtained by considering the pressure of the grid cell, along with grid cells further out along a given axis, to form an equation of the pressure field about the cell, along that axis. This approach results in a more smoothly varying pressure field than the simple two neighbour approximation. A very similar method to estimating the pressure gradient at particle $i$ is used in SPH, i.e. a weighted average of the estimated pressure gradient from particle $i$'s neighbours is used.

General SPH Equations and Derivatives

In SPH, the standard equation for estimating a physical quantity $A$, at position $i$, is given by:

\begin{displaymath} A_{i}=\sum_{j=1}^{N_{N}}m_{j}\frac{A_{j}}{\rho_{j}}W(r_{ij},h_{j}) \end{displaymath} (2.14)

where position $i$ is usually the position of a particle, and the sum is over all $i$'s neighbours, including the contribution from $i$. $A_{j}$ is the desired quantity $A$ associated with particle $j$. When $A$ is the density, Equation 2.2 results. A similar approach is adopted for calculating the gradient of $A$ at $i$:
\begin{displaymath} \mathbf{\nabla}{A_{i}}=\sum_{j=1}^{N_{N}}m_{j}\frac{A_{j}}{\rho_{j}} \mathbf{\nabla}W(r_{ij},h_{j}) \end{displaymath} (2.15)

where $\mathbf{\nabla} W$ is the first derivative of the kernel with respect to distance. For the Gaussian this is:
\begin{displaymath} \mathbf{\nabla} W_{G}(r,h)=-2\frac{r}{h^{2}}W_{G}(r,h) \mathbf{\hat{r}} \end{displaymath} (2.16)

and for the $\beta$ spline:
\begin{displaymath} \mathbf{\nabla} W_{\beta}(r,h)=\frac{\sigma}{h^{\nu+1}} \mat... ...textrm{\hspace{22pt}$u\geq 2$} \end{array}\right. \vspace{3mm} \end{displaymath} (2.17)

The notation is as described in Section 2.2, where both kernels and their derivatives were compared in Figure 2.3.

Although Equation 2.15 can be used as it stands, higher accuracy is obtained by rewriting the gradient of $A$ to include the density, usually as $\mathbf{\nabla}(\rho A)$ or $\mathbf{\nabla}(\frac{A}{\rho})$, and applying Equation 2.15 to the result. This requires that:

\begin{displaymath} \mathbf{\nabla} A=\frac{\mathbf{\nabla}(\rho A)-A\mathbf{\nabla}\rho}{\rho} \end{displaymath}

leading to:
\begin{displaymath} \mathbf{\nabla} A_{i}=\frac{\sum_{j=1}^{N_{N}}m_{j}(A_{j}-A_{i}) \mathbf{\nabla} W(r_{ij},h_{j})}{\rho_{i}} \end{displaymath} (2.18)

or that:

\begin{displaymath} \mathbf{\nabla} A=\rho\mathbf{\nabla}\Bigg(\frac{A}{\rho}\Bigg)+\frac{A}{\rho}\mathbf{\nabla}\rho \end{displaymath}

leading to:
\begin{displaymath} \mathbf{\nabla} A_{i}=\sum_{j=1}^{N_{N}}m_{j}\Bigg(\frac{\rh... ...}+ \frac{A_{i}}{\rho_{i}}\Bigg)\mathbf{\nabla} W(r_{ij},h_{j}) \end{displaymath} (2.19)

Equation 2.18 is usually used when calculating the divergence of the velocity field2.5. Equation 2.19 is preferred in this work, and in general for calculating the pressure gradient, as it results in a symmetrical force pair calculation, once the smoothing lengths have been averaged. The pressure component of the acceleration of particle $i$ can now be expressed as:
\begin{displaymath} \frac{d\mathbf{v}_{i}}{dt}=-\sum_{j=1}^{N_{N}}m_{j}\Bigg(\fr... ...rac{P_{i}}{\rho_{i}^{2}}\Bigg)\mathbf{\nabla} W(r_{ij},h_{ij}) \end{displaymath} (2.20)

where $h_{ij}$ has been introduced as a means of symmetrising the force between $i$ and $j$, and is either the arithmetic, or geometric mean of $h_{i}$ and $h_{j}$. Another approach to symmetrising the force is kernel averaging, $W(r_{ij},h_{ij})=\frac{1}{2}(W(r_{ij},h_{i})+W(r_{ij},h_{j}))$, and similarly for $\mathbf{\nabla} W$.

The simplest form is used in this work, $h_{ij}=\frac{1}{2}(h_{i}+h_{j})$, due to there being no discernible advantage of any one method over another. The smoothing length average also applies to the density calculation, Equation 2.2, to avoid a violation of the reciprocity principle (Serna et al., 1996; Evrard, 1988).

Artificial Viscosity

Artificial viscosity is included as a means of dissipating macro-scale translational energy in the system, into micro-scale energy of the molecules. In most astrophysical problems, the molecular viscosity of the gas is very small, such that the dissipation should ideally only occur in shocks. For example, without some form of energy dissipation, supersonic colliding streams of gas would penetrate each other quite readily, as the pressure present would never be sufficient to prevent this from happening, Figure 2.6(a).

Figure 2.6: Colliding streams of gas. The depth of hue is an indication of the temperature present. (a) Without viscosity, (b) With viscosity.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/viscosa1.eps,angle=0,width=1.0\linewidth}} \end{figure}

A more realistic situation, shown in Figure 2.6(b), is rapid particle deceleration in the shock interface, during which the energy of motion is transformed into the internal energy of the gas. This has the effect of increasing the density and temperature in the shock interface, and correspondingly the pressure, leading to the inflow of gas being further reduced. Preventing particle interpenetration is an important issue in SPH, as this can lead to physically unrealistic regions of space, containing particles with locally uncorrelated velocities and temperatures.

Many different forms of the artificial viscosity have emerged recently, each one being more suited to a particular type of problem (Watkins et al., 1996; Navarro and Steinmetz, 1997; Selhammar, 1997). It is likely that in the future one of these methods will become the standard prescription, but at present however, the standard equations of Monaghan and Gingold (1983) given below, are used here and by most workers in the field. This form of the artificial viscosity was designed with shocks in mind, and hence many of the terms are related to the shock phenomena.

\begin{displaymath} \mathbf{a}_{i}^{visc}=-\sum_{j=1}^{N_{N}}m_{j}\Pi_{ij}\mathbf{\nabla}W(r_{ij},h_{ij}) \end{displaymath} (2.21)

where
\begin{displaymath} \Pi_{ij}=\frac{-\alpha\xi_{ij}\overline c_{ij} + \beta\xi_{ij}^{2}}{\overline \rho_{ij}} \end{displaymath} (2.22)


\begin{displaymath} \xi_{ij}=\left\{\begin{array}{cl} \displaystyle \frac{\mathb... ...mathbf{v}_{ij}\cdot\mathbf{r}_{ij}\geq 0\\ \end{array}\right. \end{displaymath} (2.23)

and $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$, $\overline c_{ij}=(c_{i}+c_{j})/2$ is the average speed of sound of particles $i$ and $j$, $\overline \rho_{ij}=(\rho_{i}+\rho_{j})/2$, and $\eta^{2} \approx 0.01$ prevents numerical divergence for small particle separations. The force only acts when particles are approaching each other, hence the zero term in Equation 2.23.

The constant viscosity parameters, $\alpha$ and $\beta$, are of order unity, and control the strength of the viscosity. The term involving $\alpha$ is linear in $\xi_{ij}$, and is associated with the effective bulk and shear viscosity of the gas. The quadratic term is necessary where strong shocks are present to prevent interparticle penetration, and can be left out otherwise. In a typical simulation, $\alpha=0.5$ and $\beta=1.0$, although values a factor of five or so larger are often used to produce bodies in static equilibrium, on timescales much shorter than would have occurred otherwise.

Equation 2.21, is readily combined with Equation 2.20, to produce the total hydrodynamical acceleration at particle $i$:

\begin{displaymath} \frac{d\mathbf{v}_{i}}{dt}=-\sum_{j=1}^{N_{N}}m_{j}\Bigg(\fr... ...}}{\rho_{i}^{2}}+\Pi_{ij}\Bigg)\mathbf{\nabla}W(r_{ij},h_{ij}) \end{displaymath} (2.24)

This equation conserves total linear momentum, due to the imposed symmetry in both the pressure and viscosity components. In addition, angular momentum is also conserved, as the forces are directed along lines joining the particles, i.e. no torques are generated.


Gravitational Forces

A distribution of point masses contributes to the gravitational acceleration at a point $i$ by:

\begin{displaymath} -\mathbf{\nabla}\phi_{i}=-G\sum_{j=1}^{N}\frac{m_{j}}{r_{ij}^{2}} \frac{\mathbf{r}_{ij}}{r_{ij}} \end{displaymath} (2.25)

where $G$ is the gravitational constant. There are two reasons why this form of the acceleration is unacceptable in SPH:

Gravitational Softening

The first reason why a point mass force cannot be used, is due to the large acceleration which can occur when particles become very close. This situation can happen quite readily in SPH, due to quantities such as position and velocity not being so correlated on scales less than a smoothing length. That is to say, even though the smoothed mass and velocity of a group of particles, may suggest that they are all travelling in the same direction with uniform particle density; individual particles in the flow could have a small component of velocity in a different direction, or be much closer together than the average particle separation. This can lead to particles accelerating towards each other indefinitely, due to the pressure force tending to zero for small particle separations, and the gravitational force tending to infinity.

One method which is sometimes used to overcome this problem is gravitational softening via a Plummer law:

\begin{displaymath} -\mathbf{\nabla}\phi_{i}=-G\sum_{j=1}^{N}\frac{m_{j}}{(r_{ij}^{2}+\epsilon^{2}_{ij})} \frac{\mathbf{r}_{ij}}{r_{ij}} \end{displaymath} (2.26)

where $\epsilon_{ij}$ is approximately an order of magnitude less than the smoothing length, and is fixed or varying in time, depending on the implementation. If the latter is the case, then $\epsilon_{ij}=\frac{1}{2}(\epsilon_{i}+\epsilon_{j})$ is commonly used.

Consistent Representation of Matter

The second reason for a replacement gravitational term is not so important, but stems from how the distribution of matter is interpreted. When deriving Equation 2.25, the particles were assumed to be point masses, with a force between particle pairs of:
\begin{displaymath} {F}_{ij}=\frac{Gm_{i}m_{j}}{r^{2}_{ij}} \end{displaymath} (2.27)

along the line joining the particles. However, in following with the spirit of SPH, the particles can be considered to be part of a continuum, and the gravitational forces modified to correspond with this view. This leads to a more natural softening, and a force which is not so sensitive to particle disorder.

The method can be explained by considering a particle pair, and how the force between the pair is modified. In the SPH representation, there are three simple ways of visualising how the gravitational force acts between the pair: Particle-Particle interaction, Sphere-Sphere interaction, and Particle-Sphere interaction. The first, both particles thought of as point masses, has already been discussed, and found to be problematic in the absence of artificial softening. The second, both particles viewed as spread out distributions of matter, is perhaps the most realistic, but it is also the most difficult to implement.

Sphere-Sphere Gravitational Softening

In this method, the kernel is used to spread out the mass of each particle in space, via the normal SPH procedure. This leads to the gravitational force being calculated from the same physical representation as the hydrodynamical force, resulting in a more consistent force calculation. If the smoothed distributions are not overlapping, $r_{ij} > (2h_{i}+2h_{j})$, then the force is given by Equation 2.27. This is because the distributions of matter are spherically symmetric, allowing the force to be calculated as if they were point masses.

Figure 2.7: Gravitational force calculation: Sphere-Sphere interaction.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/overgrava1.eps,angle=0,width=0.8\linewidth}} \end{figure}

If the kernels associated with each sphere are overlapping, then the force given by Equation 2.27 is attenuated by an amount that depends on the degree of overlap. An example of this scenario is shown in Figure 2.7. The force between the spheres, which is subsequently applied to the particles, is still directed along the line joining the particles, but is reduced by an amount that depends on $r_{ij}$, $h_{i}$, and $h_{j}$. It is very difficult to find an analytical solution to this problem, due to the $\beta$ spline kernel being defined in three regions. However, a numerical solution was found, and reduced to a two dimensional lookup table. The procedure for doing this is not given here, as it was concluded that the final method given next, gave very similar results in typical experimental conditions, and was by far the simpler of the two.

Point-Sphere Gravitational Softening

The third method is a mixture of the first two: One of the particles is considered to be a point mass, and the other is considered to be spread out in space. This method is currently the most commonly used, due to it being straightforward to implement, and only mildly conflicting with the hydrodynamical interpretation of the mass distribution. For these reasons, it is the method used in this thesis.

When considering the acceleration on a point mass $i$, due to a spherically symmetric distribution of matter $j$, only the enclosed mass of $j$, contained in a concentric sphere of radius $r_{ij}$, need be considered. An example of this is shown in Figure 2.8, where the shaded region denotes the enclosed mass of sphere $j$, with respect to particle $i$.

Figure 2.8: Gravitational force calculation: Particle $i$-Sphere $j$ interaction.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/grava2.eps,angle=0,width=0.6\linewidth}} \end{figure}

The enclosed mass of a sphere distributed according to the kernel is given by:
\begin{displaymath} M_{j}(r,h)=\int^{r}_{0}4\pi r^{2}m_{j}W(r,h)dr \end{displaymath} (2.28)

where $r$ is the distance from the centre of the sphere. For the $\beta$ spline kernel in three dimensions, this results in:
\begin{displaymath} M_{j}(r,h)=m_{j} \left\{\begin{array}{cl} \frac{4}{3}u^{3}-\... ...\leq 2\\ 1&\textrm{\hspace{22pt}$u\geq 2$} \end{array}\right. \end{displaymath} (2.29)

where $u=r/h$. Substituting this for $m_{j}$ in Equation 2.25, the acceleration at $i$ becomes:
\begin{displaymath} -\mathbf{\nabla}\phi_{i}=-G\sum_{j=1}^{N}\frac{M_{j}(r_{ij},h_{ij})}{r_{ij}^{2}} \frac{\mathbf{r}_{ij}}{r_{ij}} \end{displaymath} (2.30)

where as usual, the smoothing lengths have been averaged to symmetrise the force. An alternative expression for this force is given by Hernquist and Katz (1989), where it is clear that singularities do not occur:
\begin{displaymath} -\mathbf{\nabla}\phi_{i}=-G\sum_{j=1}^{N}m_{j}g(r_{ij},h_{ij})\mathbf{r}_{ij} \end{displaymath} (2.31)

where $g$ is given by:
\begin{displaymath} g(r,h)= \left\{\begin{array}{cl} \frac{1}{h^{3}}\Big[\frac{4... ...c{1}{r^{3}}&\textrm{\hspace{22pt}$u\geq 2$} \end{array}\right. \end{displaymath} (2.32)

Also given is the gravitational potential at a point, which is used in Section 2.10 to calculate the instantaneous gravitational energy of the system:
\begin{displaymath} \phi_{i}=-G\sum_{j=1}^{N}m_{j}f(r_{ij},h_{ij}) \end{displaymath} (2.33)

where $f$ is given by:
\begin{displaymath} f(r,h)= \left\{\begin{array}{cl} -\frac{2}{h}\Big[\frac{1}{3... ...\frac{1}{r}&\textrm{\hspace{22pt}$u\geq 2$} \end{array}\right. \end{displaymath} (2.34)

As expected, differentiating Equation 2.34 with respect to $r$ and dividing by $-r$ results in Equation 2.32.

In this work, the $h$ in the above equations related to gravity is in general the same as in the hydrodynamical calculations. Some implementations of SPH prefer to set a lower limit on the gravitational smoothing length, or keep it constant in time, but not necessarily space. Setting a minimum value of $h$ is similar to setting a maximum value of the density and, indirectly, a minimum value for the time-step. This has advantages in simulations where material is flowing into a region of space, and it is the nature of the flow that is important, and not the resulting high density region. The reason for avoiding high density regions in this example, is because the simulation would grind to a halt due to small time-steps being required to model the high density region correctly. Another reason for keeping the gravitational smoothing length constant in time, is that it does not introduce complications when considering the energy of the system, as discussed in Section 2.10.


The Equation of State

When using the hydrodynamic acceleration derived earlier, Equation 2.20, in addition to a numerical value for the density of the gas, a numerical value for the pressure must also be available. This is achieved via the equation of state. Three different forms of equation of state are used in this thesis. The first is derived from the ideal gas law, and is expressed as:

\begin{displaymath} P_{i}=(\gamma -1)u_{i}\rho_{i} \end{displaymath} (2.35)

where $\gamma$ is the ratio of specific heat capacities of the gas, and $u_{i}$ is the internal energy per unit mass of particle $i$, also referred to as the specific internal energy. The second is a polytropic equation of state:
\begin{displaymath} P_{i}=K\rho_{i}^{\gamma} \end{displaymath} (2.36)

where $K$ is a measure of the entropy of the gas, and is globally constant when used in this work. In some implementations (Bate, 1995), $K$ is allowed to vary both globally and locally, via a $dK_{i}/dt$ term similar in form to the $du_{i}/dt$ term explained in Section 2.10.

These first two equations both assume that the ratio of specific heats, $\gamma$, is unchanging throughout the calculation. The third equation has no such restriction, as the temperature appears explicitly in the equation of state:

\begin{displaymath} P_{i}=\frac{k\rho_{i}T_{i}}{\mu_{i}} \end{displaymath} (2.37)

where $k$ is Boltzmann's constant, $T_{i}$ is the temperature, and $\mu_{i}$ is the mean molecular weight of the material associated with particle $i$ given by:
\begin{displaymath} \frac{1}{\mu_{i}}=\Bigg[\frac{X}{(2-a_{i})\mu_{H}}+\frac{Y}{\mu_{He}}\Bigg] \end{displaymath} (2.38)

where $a_{i}$ is the degree of dissociation of hydrogen molecules in the gas associated with $i$, and $X$ and $Y$ are the fractions by weight of hydrogen and helium respectively. In this work, globally constant values of $X=0.7$ and $Y=0.3$ are used. In theory, heavier elements should be included into this calculation, but their contribution is small, and are justifiably ignored for the pressure calculation.

The temperature of a gas is related to the translational component of the kinetic energy of its constituent molecules. At very low temperatures $\approx 10K$, only translational degrees of freedom are present, resulting in all the internal energy of the gas being associated with the thermal component. The ratio of specific heats is related to the number of degrees of freedom, $N_{D}$ by:

\begin{displaymath} \gamma=\frac{N_{D}+2}{N_{D}} \end{displaymath} (2.39)

leading to $\gamma=5/3$ at low temperatures. For a monatomic gas, such as helium, this is the appropriate value to use at all temperatures. However, for a molecular gas the value is a function of temperature, and to some extent pressure. For hydrogen in particular, as the temperature increases, the rotational degrees of freedom perpendicular to the line joining the atoms, slowly thaw. This accounts for a further two degrees of freedom, decreasing the ratio of specific heats to $\gamma=7/5$ once they are fully active when $T\approx300K$. As the temperature increases further, the energy of vibration along the line joining the atoms, slowly comes into equipartition, tending the ratio of specific heats to $\gamma=9/7$. Before these two2.6additional degrees of freedom have reached equipartition however, the hydrogen molecules begin to dissociate into their component atoms. At low pressures this occurs around $2000K$, and requires about $10$ times the total internal energy of the molecule at that temperature. Once all the hydrogen has dissociated, the ratio of specific heats is once again $\gamma=5/3$.

The standard SPH formalism provides for the specific internal energy of the gas to be followed in time, via a $du_{i}/dt$ term, as shall be seen in Section 2.10. A means of expressing the temperature and mean molecular weight of the gas, as a function of its specific internal energy is therefore required.

Experimental results are available for the specific heat capacity as a function of temperature at constant pressure for many gases. The source used here is a compilation of numerous published works, weighted by how reliable they were judged to be (Touloukian and Ho, 1970). Two sets of values are available: The first for a pressure of one atmosphere, and the second for an interpolated estimate of zero pressure. The latter is the appropriate one to use in these simulations, however for hydrogen both sets of results are within $0.1\%$ of each other. Two sets of data were supplied, because in general for other materials the correlation between $C_{P}^{0}$ and $C_{P}^{1}$ is not so close.

The specific heat capacity of molecular hydrogen at zero pressure is given by:

$\displaystyle C_{P}^{0}(J\,kg^{-1}K^{-1})_{100K\to400K}$ $\textstyle =$ $\displaystyle 6.1508\times10^{3}$ (2.40)
$\displaystyle +(6.7013\times10^{1})T$ $\textstyle -$ $\displaystyle (1.8591\times10^{-1})T^{2}+(1.7636\times10^{-4})T^{3}$  
$\displaystyle C_{P}^{0}(J\,kg^{-1}K^{-1})_{400K\to1500K}$ $\textstyle =$ $\displaystyle 1.4943\times10^{4}$ (2.41)
$\displaystyle -2.0498T$ $\textstyle +$ $\displaystyle (2.6065\times10^{-3})T^{2}-(5.0110\times10^{-7})T^{3}$  

A simple cubic function is fitted between $0K$ and $100K$, constructed from the theoretical value at $0K$: $\gamma=5/3\Rightarrow C_{P}=10.318kJ\,kg^{-1}K^{-1}$, and the value of Equation 2.40 at $100K$: $C_{P}=11.169kJ\,kg^{-1}K^{-1}$:
\begin{displaymath} C_{P}(J\,kg^{-1}K^{-1})_{0K\to100K}=1.0318\times10^{4}+(8.51\times10^{-4})T^{3} \end{displaymath} (2.42)

Equation 2.41 is used between $1500K$ and $2000K$, as the behaviour of it is realistic in this region, and by $2000K$ it has not yet exceeded the maximum theoretical value: $\gamma=9/7\Rightarrow C_{P}=18.572kJ\,kg^{-1}K^{-1}$.

Figure 2.9: Specific heat capacity as a function of temperature. (a) Hydrogen at constant pressure, (b) Hydrogen at constant volume, (c) Hydrogen and helium mix at constant volume. Beyond $2000K$ the specific heat capacities remain constant.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/heatcapa3.eps,angle=0,width=0.95\linewidth}} \end{figure}

For simplicity in this work, above $2000K$ all of the hydrogen is assumed to have dissociated, resulting in a specific heat capacity twice that at $0K$. This is because the number of degrees of freedom has reduced back to three, and there are now twice as many bodies per kilogram. The degree of dissociation is actually pressure dependent: At a pressure of one atmosphere, $50\%$ dissociation has occurred by $4000K$, and total dissociation by $6000K$; for a tenth of an atmosphere, the corresponding temperatures are $3500K$ and $5000K$ (Schofield, 1979). These values tend towards $2000K$ for zero pressure, and is the reason for using that temperature as a simplified cross over point in this work. The pressure present in hydrostatic gas bodies with constant central temperature scales as $M^{-2}$. A Jupiter mass protoplanet in quasi-static collapse, with a central temperature of $2000K$, has a central pressure of around a tenth of an atmosphere. This cross over approximation should therefore be applied with caution when studying the evolution of small mass bodies beyond $2000K$. However, for larger mass protoplanets, and for protostars, it is a reasonably valid assumption. The specific heat capacity as a function of temperature is shown in Figure 2.9.

The specific internal energy as a function of temperature, can be found by integrating the heat capacity equations with respect to temperature in the appropriate ranges. For example, the specific internal energy of the hydrogen-helium mix at $1000K$ is given by the area under the curve between $0K$ and $1000K$ in Figure 2.9(c). The specific heat capacity at constant volume is used in this calculation, due to the energy of expansion already being accounted for in the SPH energy equations:

\begin{displaymath} C_{V}=C_{P}-\mathcal{R} \end{displaymath} (2.43)

where $\mathcal{R}$ is the molar gas constant, and accounts for the difference in energy just mentioned.

In the above integration, the energy of dissociation of a hydrogen molecule, $4.476eV$, must be included. For the given mix of gases, and current dissociation approximation, this corresponds to a sink of energy at $2000K$ equivalent to $150MJkg^{-1}$. Whilst the molecules are absorbing this energy, the mean molecular weight gradually reduces in value. This is included in a linear fashion, i.e. when $75MJ$ has been supplied to a kilogram of gas, half of the hydrogen molecules present in that kilogram are assumed to have dissociated. The mean molecular weight and temperature can be expressed as a function of specific internal energy by solving the resulting equations numerically, Figure 2.10.

Figure 2.10: Dependence on specific internal energy. (a) Mean molecular weight, (b) Temperature.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/tasua2.eps,angle=0,width=1.0\linewidth}} \end{figure}

One advantage of this method, is that other sources of latent heat can be included in the same manner as the energy of dissociation. In fact the latent heat of evaporation of icy grains was included at low temperatures, but found to have negligible effect, and was omitted in this presentation for clarity2.7. The dependences of $\mu$ and $T$ on $u$, are stored in one dimensional look up tables for use during a simulation.

Some implementations of SPH use variations on the polytropic equation of state, Equation 2.36, to simulate both the dissociation of hydrogen, and the varying opacity of the material (Bate, 1998; Bonnell, 1994). Normally a piece-wise polytropic equation of state, divided according to density is used. Each density regime has an effective value of $\gamma$ that controls the behaviour of the gas in that regime, e.g. $\gamma=1.0$ results in the gas being isothermal, and is appropriate at low densities. At high densities, the gas is diatomic and opaque, and $\gamma=7/5$ is used. At higher densities still, a value of $\gamma\approx1.15$ accounts for the energy required to dissociate hydrogen molecules. Once the dissociation is complete, the last density range corresponds to an opaque atomic gas, where $\gamma=5/3$ is appropriate.

Procedures like this are required in simulations that contain both optically thin and optically thick regions of space. This is due to the standard SPH equations lacking a satisfactory radiative energy transfer scheme. More is said on this in Section 3.2, where a new general radiative energy transfer scheme is presented.


Energy Equations

So far, only the forces acting on the particles have been examined. How these forces affect the properties of the particles is the subject of this section and the next. In this section, the differential equations that govern how the properties of the particles evolve in time are derived. Also derived are the global energy equations that keep track of the different forms of energy in the system, and can be used to give an indication of how accurately the integration is performing. The method for integrating the particle properties in time is explained in Section 2.11.

In this implementation of SPH, there are three intrinsic properties of the particles that must be followed in time: Position, velocity and the internal energy of the particles if a non polytropic equation of state is used. Differential equations for updating the velocity have already been derived from the forces, Equation 2.8. The change in position with time is equivalent to the velocity, thus:

\begin{displaymath} \frac{d\mathbf{x}_{i}}{dt}=\mathbf{v}_{i} \end{displaymath} (2.44)

The change in internal energy per unit mass can be derived from the first law of thermodynamics, applied to an inviscid adiabatic process:
\begin{displaymath} \frac{du}{dt}=-\Bigg(\frac{P}{\rho}\Bigg)\mathbf{\nabla}\cdot\mathbf{v} \end{displaymath} (2.45)

This essentially says that the change in volume is proportional to the change in internal energy of the gas. Applying Equation 2.18 to Equation 2.45 results in:
\begin{displaymath} \frac{du_{i}}{dt}=\frac{P_{i}}{\rho_{i}^{2}}\sum_{j=1}^{N_{N}}m_{j}\mathbf{v}_{ij}\cdot\mathbf{\nabla}W(r_{ij},h_{ij}) \end{displaymath} (2.46)

Some implementations of SPH reformulate Equation 2.45 as:
\begin{displaymath} \frac{du}{dt}=-\mathbf{\nabla}\Bigg(\frac{P\mathbf{v}}{\rho}\Bigg)+\mathbf{v}\cdot\mathbf{\nabla}\Bigg(\frac{P}{\rho}\Bigg) \end{displaymath} (2.47)

and apply Equation 2.18, resulting in:
\begin{displaymath} \frac{du_{i}}{dt}=\sum_{j=1}^{N_{N}}m_{j}\Bigg(\frac{P_{j}}{... ...{2}}\Bigg) \mathbf{v}_{ij}\cdot\mathbf{\nabla}W(r_{ij},h_{ij}) \end{displaymath} (2.48)

and then average this with Equation 2.46. The reason for this approach is that it results in the same equation as found by applying energy conservation2.8 to the inviscid pressure acceleration of Equation 2.20:
\begin{displaymath} \frac{du_{i}}{dt}=\frac{1}{2}\sum_{j=1}^{N_{N}}m_{j}\Bigg(\f... ...^{2}}\Bigg)\mathbf{v}_{ij}\cdot\mathbf{\nabla}W(r_{ij},h_{ij}) \end{displaymath} (2.49)

This symmetrised form of the energy equation, although equally valid in principal to Equation 2.46, can produce negative particle internal energies. The reason for this was pointed out by Benz (1990), and is due to $du_{i}/dt$ being overestimated in some cases. This is because the pressure at particle $j$ enters into the $du_{i}/dt$ calculation of particle $i$, and if larger in magnitude can force $u_{i}$ negative. Although using Equation 2.49 with the integration scheme implemented here will not produce negative energies, it can lead to very small time-steps being predicted, and thus Equation 2.46 is preferred. Combined with the artificial viscosity, Equation 2.46 becomes:
\begin{displaymath} \frac{du_{i}}{dt}=\sum_{j=1}^{N_{N}}m_{j}\Bigg(\frac{P_{i}}{... ..._{ij}\Bigg)\mathbf{v}_{ij}\cdot\mathbf{\nabla}W(r_{ij},h_{ij}) \end{displaymath} (2.50)

The total energy of the system, $E_{S}$, can be monitored by combining the individual components together. The appropriate energies that require to be considered are: Kinetic energy, internal energy, and gravitational potential energy:
\begin{displaymath} E_{S}=E_{K}+E_{U}+E_{\Phi} \end{displaymath} (2.51)

The total kinetic energy is given by:
\begin{displaymath} E_{K}=\frac{1}{2}\sum_{i=1}^{N}m_{i}v_{i}^{2} \end{displaymath} (2.52)

The total internal energy when using $du/dt$ is given by:
\begin{displaymath} E_{U}=\sum_{i=1}^{N}m_{i}u_{i} \end{displaymath} (2.53)

The thermal component of the internal energy is $\frac{3}{2}(\gamma-1)$ times that of Equation 2.53, or can be calculated directly from the temperature:
\begin{displaymath} E_{T}=\frac{3}{2}k\sum_{i=1}^{N}\frac{m_{i}T_{i}}{\mu_{i}} \end{displaymath} (2.54)

When using a polytropic equation of state, the total internal energy is given by:
\begin{displaymath} E_{U}=\frac{K}{\gamma-1}\sum_{i=1}^{N}m_{i}\rho_{i}^{\gamma-1} \end{displaymath} (2.55)

Two ways of accounting for the gravitational energy in the system are used here. The first is an instantaneous measurement, obtained from the distribution of the particles. The second is similar in approach to how the internal energy is updated throughout the calculation.

Given that the mass of the particles are spread out by the $\beta$ spline kernel, and that the gravitational force is applied using this assumption, then the gravitational energy due to the separation of the particles is given by:

\begin{displaymath} E_{\Phi}=-G\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}m_{i}m_{j}f(r_{ij},h_{ij}) \end{displaymath} (2.56)

where $f$ was given in Equation 2.34. In addition to this potential energy of separation, each particle has an effective self gravitational potential energy associated with it. This arises from the spread nature of the particles, and is inversely proportional to the gravitational smoothing length. In calculations that maintain a constant gravitational smoothing length, the gravitational self energy of the particles remains constant, and can be ignored. However, if the smoothing lengths vary in time, then the change in self potential energy of the particles must be accounted for, not only in the energy equations, but also in the force equations for consistency.

The extra terms that are introduced, coined the $\nabla h$ terms, describe how the force between a pair of particles varies with $h_{ij}$. In theory, $\nabla h$ terms should also be included in the hydrodynamic force and energy equations. Nelson and Papaloizou (1994) have investigated the inclusion of the hydrodynamic $\nabla h$'s, and concluded that they are correction terms to the standard equations, and that omitting them results in non conservation of energy at the $10\%$ level over the course of a typical simulation2.9.

Even though the $\nabla h$ terms have not been included in this work, the total measured energy of the system remains constant in the absence of gravitational forces. This is due to the internal energy being updated consistently from the hydrodynamic force. The gravitational energy calculated from Equation 2.56 is only consistent with the gravitational force equation when the smoothing lengths are kept constant in time. If the total energy is required to remain constant in the presence of gravitational forces, as a check on the performance of the integration say, then a similar approach to updating the internal energy can be adopted for updating the gravitational energy. Initially the gravitational energy of the system is calculated using Equation 2.56. It is then subsequently integrated forward in time via the differential equation: