Subsections


Testing the Code

Before a numerical method can be confidently used to obtain new results, it must be tested to see how well it can reproduce known results. This chapter is a presentation of the tests that have been performed after each new major addition to the code, given in roughly the order of development.

Gravity

A standard test of gravitational forces is the free-fall collapse of a cloud of material due to self gravity alone (e.g., Golanski, 1999; Turner et al., 1995). In this test, a uniform-density spherical distribution of matter with mass $M$ and radius $R$, collapses homogeneously to a point in the free-fall time:

\begin{displaymath} t_{ff}=\frac{\pi}{2}\Bigg(\frac{R^{3}}{2GM}\Bigg)^{\frac{1}{2}} \end{displaymath} (4.1)

During the collapse, any element of the cloud which at time $t=0$ has radius $r_{0}$, at later times in the interval $0\le t \le t_{ff}$, has radius $r$ given by:
\begin{displaymath} \frac{t}{t_{ff}}=1-\frac{2}{\pi}\Bigg\{\sin^{-1}\Bigg[\bigg(... ...\frac{1}{2}}\bigg[1-\frac{r}{r_{0}}\bigg]^{\frac{1}{2}}\Bigg\} \end{displaymath} (4.2)

These equations are derived by considering how the motion of the surface varies with time, see Appendix B.1.

Initial Configuration

Several types of simulation conditions are possible for this test, allowing for the opportunity to examine various aspects of the SPH technique. The first, using a uniform grid and point mass force, is the most popular as it reproduces the free-fall theoretical results more closely than the others. In this test, particles of equal mass are positioned on the lattice points of one of the common crystal structures: Hexagonal closest packed (HCP), body centred cubic (BCC), or simple cubic (SC). An initial BCC configuration is used here, and throughout this thesis. This choice is based on two factors: The degree to which the particles resemble a spherical body after the grid has been truncated at a given radius from the central particle, and the size of the average interparticle separation - quantified by the packing fraction4.1 of the crystal structure.

The first of the reasons above, comes from the fact that most of the tests and simulations in this thesis have initial spherically symmetric starting configurations. The second is related to minimising the potential energy associated with the separation of the particles in equilibrium bodies (when pressure is present), i.e. a perturbation to an equilibrium body possessing a small packing fraction, sufficient to overcome the `potential barrier' associated with the large interparticle separation, could result in the particles adopting a new more stable configuration with a globally lower potential energy, and smaller interparticle separations. Large packing fractions are therefore desirable when initial stable equilibrium configurations are required.

Figure 4.1: Stereo representations of HCP and BCC spherical structures. Left and middle column - wide eyed view; Middle and right column - cross eyed view. From top to bottom: (a) $675$ particle HCP sphere, x-y plane, (b) Arbitrary plane, (c) $531$ particle BCC sphere, x-y plane, (d) Arbitrary plane.
\begin{figure}\begin{center} \begin{tabular}{@{}c@{}c@{}c@{}} \vspace{-5pt} \eps... ...carbr.eps,angle=0,width=0.33\linewidth}\\ \end{tabular}\end{center}\end{figure}

Figure 4.1 shows the structure present for a relatively small number of particles in an HCP and BCC sphere.

Although an HCP structure has the largest packing fraction ($0.74$), it lacks a desirable symmetry along Cartesian axes which complicates analysis of results, and when truncated to form a sphere it produces bodies which are noticeably non-spherical, especially when small numbers of particles are employed. An SC structure has the desirable symmetry, but a low packing fraction ($0.52$) means it is to be avoided. The BCC structure has both the desirable symmetry and a relatively high packing fraction ($0.68$). Other types of crystallographic structures were not investigated: It is suspected that FCC (Face centred cubic) and CCP (Cubic closest packing) may have the desirable features just discussed.


Uniform Grid - Point Mass Force

In these tests, the distances from the centre of the sphere of three particles are followed. Progressing outwards from the centre of the sphere, the first particle which encloses $10\%$ of the spheres mass is chosen. The same procedure is used to choose particles which enclose $50\%$ and $90\%$ of the spheres mass. These varied values allows the dependence on initial radius, and hence surface effects, to be investigated. Two simulations are performed for this uniform grid - point mass force scenario: A low resolution simulation containing $531$ particles, and a high resolution simulation, with approximately twice the spatial resolution, containing $4279$ particles. For simplicity, the simulated region is equal to one solar mass, and has a radius of one astronomical unit. This leads to a theoretical free-fall time $t_{ff}=1/\sqrt{32} \approx 0.177\,yr$, and an initial theoretical density $\rho=1.42\times 10^{-4}kg\,m^{-3}$.

Figure 4.2 compares the theoretical radii with the experimental radii as a function of time.

Figure 4.2: Main free-fall curves for $4279$ particles arranged on a uniform grid interacting via gravitational point forces alone. The solid lines represent the theoretical radius as a function of time for the experimental $10\%$, $50\%$ and $90\%$ mass radii shown. The experimental results are plotted every $5$ time-steps.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/freefall01a1.eps,angle=0,width=0.87\linewidth}} \end{figure}

The theoretical curves are calculated using the initial radial distance of the particles in Equation 4.2, so that the theoretical and experimental curves coincide at $t=0$. Over the range of times shown, it can be seen that the experimental results closely follow the theoretical curves. A magnified view near the end of the simulation, where the correspondence is not so close, is shown later in Figure 4.6(d) in Section 4.1.5. The spatial variance of the particles in this experiment is shown in Figure 4.3.

Figure 4.3: Spatial particle variance during the free-fall test: $4279$ particles placed on a uniform grid - point mass force. The scale marked on each frame corresponds to $0.1AU$
\begin{figure}\begin{center} \begin{tabular}{c@{\quad}c} \epsfig{figure=/home/st... ... (d)&0.40&0.97&0.31&$\infty$&380\\ \hline \end{tabular}\end{center}\end{figure}

Due to the symmetry in the simulation, approximately half of the circles in this figure correspond to more than one particle. The table at the bottom of the figure summarizes the simulation properties in each frame: In this simulation, a point mass force is produced by setting the number of neighbours of a particle to $1$, in order to remove the softening effect in the gravitational acceleration, Equation 2.30. The effect of this is to give initial particle densities which are approximately twice the theoretical initial density4.2. This also leads to the minimum interparticle separation during this experiment being directly related to the maximum density present, and likewise the maximum interparticle distance and the minimum density are directly related. Thus it is clear from the infinity in the table relating to Figure 4.3(d) that the experiment was halted due to two particles becoming very close.

In these simulations, a small tolerance for time-step prediction is employed, $\epsilon_{r}=0.001$, so as to retain any particle symmetry that may be present. This is necessary for point force simulations, as any small deviation from symmetry will result in particle-particle singularities forming early in the collapse.

The evolution of the particle positions with respect to each other can be examined by considering the frames in Figure 4.3. In a grid of infinite extent, it is clear that the initial uniform distribution of particles on a grid will be maintained throughout the collapse. In this simulation, it is the inability of discrete particles to model a spherical surface, which leads to the behaviour shown, uniform collapse, during which a shell of non-homogeneity gradually propagates radially inwards, corrupting the grid structure. The simulation halts when two particles form a singularity somewhere in the region behind the shell. Inside the propagation shell, particles are still positioned approximately on grid points. This propagating shell effect is more clearly demonstrated in the next experiment, where the use of softened gravitational forces allows the collapse to be followed for longer.

Uniform Grid - Softened Force

Figure 4.4: Spatial particle variance during the free-fall test: $4279$ particles placed on a uniform grid - softened force. The scale marked on each frame corresponds to $0.1AU$. The value in brackets corresponds to the approximate total size of the system at that time.
\begin{figure}\begin{center} \begin{tabular}{c@{\quad}c} \epsfig{figure=/home/st... ...)&0.15(0.4)&0.995&0.001&50.0&260\\ \hline \end{tabular}\end{center}\end{figure}

This scenario tests the code in the context of gravitationally softened forces. The scenario is exactly the same as the previous one, except that the number of neighbours of a particle is set to $55$. Two new features are introduced, one with negative aspects, and one with positive aspects. The positive feature is that the simulation is not halted due to singularities forming. Unfortunately, the softening that prevents the singularities from forming, is also responsible for additional, more dominant, surface effects: It is widely accepted that in SPH, surfaces require additional considerations in order that they may be modelled correctly. This is because the SPH representation of matter cannot describe arbitrarily high density and force gradients. In this example, the density is uniform across the majority of the sphere, gradually reducing to zero approximately two smoothing lengths beyond the surface, in a manner similar to that shown in Figure 2.2(b) in Chapter 2. Using a constant number of neighbours amplifies this effect, as the smoothing lengths of the surface particles are larger than the rest of the system. In a grid of infinite extent, the correct free-fall time and uniform grid structure would be observed even though the force between individual particles is softened.

The spatial variance of the particles in this experiment is shown in Figure 4.4. In the previous experiment, the corrupting shell effectively only displaced the particles from the grid points. Here, not only is this displacement effect present, but due to the force being softened at the surface of the sphere, the shell also propagates a region of reduced density radially inwards. As the shell moves inwards it becomes cube-like, Figure 4.4(d), verified by examining this figure from arbitrary angles. Due to the small tolerance, the symmetry throughout the simulation is preserved, i.e. no chaotic motion is introduced: An arbitrarily high precision integrator would have produced the same results. Thus, as in the previous experiment, approximately half of the circles correspond to more than one particle.

See also Figure 4.6(e) in Section 4.1.5, for further results from this experiment.

Displaced Grid - Softened Force

Figure 4.5: Spatial particle variance during the free-fall test: $4279$ particles placed on a displaced grid - softened force. The scale marked on each frame corresponds to $0.1AU$. The value in brackets corresponds to the approximate total size of the system at that time.
\begin{figure}\begin{center} \begin{tabular}{c@{\quad}c} \epsfig{figure=/home/st... ...&0.20(0.65)&0.99&0.0001&48.4&370\\ \hline \end{tabular}\end{center}\end{figure}

This is the final free-fall test scenario. The other possible scenario, displaced grid - point mass force, is omitted as it is clearly not feasible, due to singularities occurring almost immediately after the simulation has commenced. The test presented here is exactly the same as the previous one, except that the particles are displaced randomly from the initial grid positions, by up to half the distance to their nearest neighbour. This procedure produces a body which is more uniform and spherical in nature, than a body formed by positioning particles at random within the sphere. These initial conditions are the most commonly used in this thesis, as it is felt that symmetrical structures can produce misleading results, e.g. although the results obtained from the previous two scenarios are in good accord with the theoretical expectations, in general the initial distribution of material, even in initially uniform spherical configurations, would naturally be slightly inhomogeneous.

The spatial variance of the particles in this experiment is shown in Figure 4.5. In this simulation, only the inwardly propagating shell of reduced density is observable. As seen, the initial random fluctuations in the original density field are magnified at later times, resulting in fragmentation. Collapse within a fragment containing more than $55$ particles eventually halts the simulation. The type of structure formed in this simulation is similar to typical structures formed when rapid cooling of gas, or large initial density fluctuations, result in gravitational instability and dynamical fragmentation. The main structural difference between simulations of that type and this simulation, is that an inward radial flow of material is not present, due to the simulation domain normally being a part of an extended medium - usually accomplished by applying boundary conditions to the domain. Two example processes exhibiting this type of behaviour are the formation of large scale structure in the universe (e.g., Salmon et al., 1994), and the formation of stars in molecular clouds (e.g., Klessen et al., 1998).

See also Figure 4.6(f) in the following section, for further results from this experiment.


Free-Fall curves

Figure 4.6: Comparison of experimental free-fall deviations. Left column, $531$ particles: (a) Uniform-point, (b) Uniform-softened, (c) Displaced-softened. Right column, $4279$ particles: (d) Uniform-point, (e) Uniform-softened, (f) Displaced-softened.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/freefalla1.eps,angle=0,width=0.90\linewidth}} \end{figure}

Figure 4.2 in Section 4.1.2 showed the experimental and theoretical radii as a function of time for $4279$ particles arranged on a uniform grid, interacting via point mass gravitational forces alone. As already stated, this scenario reproduces the free-fall theoretical results more closely than any of the other scenarios. Initially every scenario follows the theoretical curves closely. It is not until later times that the deviation is noticeable. The deviations for each of the three simulations presented in the previous sections are shown in Figure 4.6. Also shown for comparison is the deviation for low resolution versions of these simulations containing $531$ particles, which is approximately half the spatial resolution of the high resolution simulations.


Hydrodynamics

There are several standard experiments available for testing the hydrodynamic component of the code. In simulations where strong shocks are present, it is the viscosity which needs to be handled correctly, and a good test of this is in the shock tube experiment (e.g., Watkins et al., 1996; Lombardi et al., 1999; Turner et al., 1995). The shock tube experiment is a purely hydrodynamic test that requires boundary conditions. Since boundary conditions have not been implemented in this code, and the opportunity of using a hydrodynamical and gravity based test is available, the shock tube test was not performed.

Once a component has been tested, and has proven to have been included correctly, it is reasonable to include and test other parts of the code alongside it, rather than testing them in isolation first. A good hydrodynamic-gravity test, which in essence provides its own boundary conditions, is the equilibrium structure and oscillation of polytropes (e.g., Golanski, 1999; Bate, 1995; Coates, 1980). For the purpose of these experiments, a polytrope is defined as being a model star, whose structure and evolution is governed by the polytropic equation of state, introduced in Equation 2.36 in Chapter 2:

\begin{displaymath} P=K\rho^{1+(1/n)} \end{displaymath} (4.3)

where $n$ is a dimensionless quantity called the polytropic index. Comparing Equation 4.3 with Equation 2.36 results in $\gamma=1+(1/n)$. If the model star is assumed to behave adiabatically, then $\gamma$ is the real ratio of specific heats as defined by the number of degrees of freedom, Equation 2.39, otherwise it is an effective ratio of specific heats. In Equation 4.3, $P$ is the pressure of the gas, and in these tests, $K$ is a scaling constant. Given that $K$ and $n$ are constant during the course of a simulation, both the experimental structure of equilibrium bodies, and the behaviour of perturbed equilibrium bodies can be compared with theory.

The Structure of Polytropes

The method used here for creating polytropes is as follows: Particles of equal mass are placed on a spherical BCC displaced grid of total radius $R$ and mass $M$. The gravitational potential energy of this arrangement is calculated from the particles using Equation 2.56 in Chapter 2 or from theory:

\begin{displaymath} \Omega=-\frac{3}{5}\bigg(\frac{GM^{2}}{R}\bigg) \end{displaymath} (4.4)

The values calculated from experiment and those calculated from theory differ by about $2.6\%$ for $1000$ particles, by about $0.59\%$ for $8000$ particles, and by about $0.11\%$ for $64000$ particles. It is desirable to use a value of $K$ close to the value that gives virial equilibrium, i.e. a value where the thermal energy present is equal to minus half the gravitational binding energy. When this is so, the thermal energy per unit mass is given by:
\begin{displaymath} u_{T}=-\frac{1}{2}\frac{\Omega}{M} \end{displaymath} (4.5)

If it is assumed that the model polytropes behave adiabatically, then a value for $K$ can be found by equating the ideal gas equation of state, Equation 2.35 in Chapter 2, with the polytropic equation of state, Equation 4.3.
\begin{displaymath} K=(\gamma-1)u\rho^{(1-\gamma)} \end{displaymath} (4.6)

Substituting for the internal energy per unit mass
\begin{displaymath} u=\frac{2}{3}\frac{u_{T}}{(\gamma-1)} \end{displaymath} (4.7)

leads to:
\begin{displaymath} K=\frac{GM}{5R\rho^{\frac{1}{n}}} \end{displaymath} (4.8)

The particles are then evolved into the desired polytropic equilibrium configuration. In doing so, the particles at the edge initially move outwards and those close to the centre move inwards. Whilst the particles are relaxing to form a polytrope, strong artificial viscosity parameters, $\alpha\approx 3$ and $\beta\approx 6$, are employed to quickly damp the ensuing oscillations.

A simple measure of how well the particles represent a polytrope can be obtained by comparing the radial variation of particle properties with theory. For a polytrope, the density as a function of radial distance is the appropriate quantity to compare. The theoretical density profile of a polytrope in equilibrium can be obtained by solving the governing Lane-Emden equations numerically, (e.g., Jeans, 1961; Coates, 1980). In the comparison procedure used here, the theoretical ratio of central density to mean density ( $\rho_{c}/\bar{\rho}$) is also required.

Three different values of the polytropic index were used for this test: $n=1$, $n=1.5$ and $n=2.5$. Although $n=1$ corresponds to an unrealistic ratio of specific heats, it serves to demonstrate the trends in behaviour. The theoretical density profiles, and the theoretical ratio of densities, were calculated simultaneously using the MAPLE programming language, and were seen to be the same as other results calculated elsewhere in the literature, Table 4.1.

Table: Ratio of central density to mean density as a function of $n$.
n $\rho_{c}/\bar{\rho}$
1.0 3.2899
1.5 5.9907
2.5 23.406


In order that the density profiles can be compared, the theoretical and experimental axes are required to be normalised to the same units. A natural choice of units is to set the theoretical central density ($\rho_{c}$) and the theoretical radius ($R$) to unity. The mass of the theoretical polytrope is then given by:
\begin{displaymath} M_{T}=\frac{4\pi}{3}\bigg(\frac{\bar{\rho}}{\rho_{c}}\bigg) \end{displaymath} (4.9)

Given that the experimental polytrope has a known mass $M$, and that the density at each SPH point is also known, the radius of the polytrope is the only other quantity required for comparison. The technique used here for determining this is novel in approach, and has the desirable side-effect of testing the gravitational potential energy calculation also. It is clear that the surface particles do not correspond to the effective surface of the polytrope, as they give finite density further out. Perhaps a better approximation, not investigated here, would be to take the distance to the surface particles plus twice (or some other factor of) their smoothing lengths.

The effective radius of a polytrope in static equilibrium can be derived by noting that its theoretical gravitational potential energy is equal to:

\begin{displaymath} \Omega=-\frac{3}{5-n}\bigg(\frac{GM^{2}}{R}\bigg) \end{displaymath} (4.10)

If $\Omega$ is calculated experimentally, the effective radius of the polytrope can be associated with the distance determined from Equation 4.10 above, referred to as $R_{\Omega}$. The degree to which $R_{\Omega}$ is larger than the distance to the surface described by the particles, depends on the polytropic index and the number of particles used to simulate the polytrope. The radial distances of the particles are matched to the theoretical results by dividing by $R_{\Omega}$. This scaling requires that the density of the particles be multiplied by $R_{\Omega}^{3}$. Finally the difference in masses have to be accounted for by multiplying the density of the particles by $M_{T}/M$.

Figure 4.7: Comparison of experimental (dots) with theoretical (solid line) density profiles. Many particles are not visible because of the close coincidence. The arrows indicate the extent of the experimental data. Left column, $531$ particles: (a) n=1, (b) n=1.5, (c) n=2.5. Right column, $4279$ particles: (d) n=1, (e) n=1.5, (f) n=2.5.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/staticpola2.eps,angle=0,width=0.90\linewidth}} \end{figure}

Figure 4.7 is a comparison of the experimental and theoretical density profiles of each of the three types of polytropes formed. The column on the left contains low resolution polytropes with $531$ particles, and the column on the right contains high resolution polytropes with $4279$ particles. The experimental data fits the theoretical curves remarkably well, despite the claim by Bate (1995) that for high polytropic indices there is significant deviation. It is certainly true that, as $n$ approaches $5$, the SPH method becomes less capable of modelling polytropes. This is because a finite number of particles cannot model an arbitrarily high density gradient. Presumably Bate was using a different method for determining the radius of the polytrope. Indeed, if he had chosen a larger effective radius, the density profiles would have coincided more closely.

So long as a polytrope in static equilibrium can be formed, and $R_\Omega$ is used as the effective radius, the experimental profile will closely coincide with the theoretical one. However, if a polytrope cannot be formed in static equilibrium, the number of particles used in the simulation must be increased, or the polytropic index decreased. For example, an $n=3$ polytrope could not be formed using $4279$ particles, despite several attempts at doing so. This was presumably due to the required density profiles being unattainable, exemplified by the central regions collapsing indefinitely. Table 4.2 shows the ratio of thermal energy to gravitational potential energy present in each of the polytropes formed.

Table 4.2: The ratio of thermal energy to gravitational potential energy of each polytrope formed. The theoretical value is equal to $0.5$.
Polytropic Index, n Number of Particles $-E_{T}/\Omega $
1.0 531 0.5049
  4279 0.4976
1.5 531 0.5004
  1067 0.4986
  1807 0.4986
  4279 0.4967
2.5 531 0.4941
  4279 0.4943


The value predicted by the virial theorem is $0.5$.

The Oscillation of Polytropes

This experiment tests several parts of the code simultaneously. In particular, it tests the dynamical interaction between the different forces, as well as testing the viscous damping, energy calculations, and the integration scheme. The tests presented here are a small subset of those possible, and indeed of those actually performed whilst developing the code. For example, the behaviour of the new integration scheme, along with each of the different equations of state described in Section 2.9, were tested using experiments similar to those described here.

If a polytrope is radially perturbed, it will oscillate about its equilibrium configuration with a period equal to the fundamental radial pulsation period. This period is of the order of the free-fall time of the polytrope:

\begin{displaymath} t_{ff}\sim(G\bar{\rho})^{-1/2} \end{displaymath} (4.11)

Given that the mass and radius of a polytrope are known, the period can be solved numerically from theory as a function of $n$, by solving the linear adiabatic wave equation (e.g., Cox, 1980; Coates, 1980). A polytrope with index $n=1.5$ is the only type of polytrope considered in these oscillation experiments, and the fundamental pulsation period calculated by Cox is used: An $n=1.5$ polytrope with a solar mass and radius pulsates with a period equal to $0.0703$ days. Using the scaling properties of Equation 4.11, the fundamental oscillation period of a polytrope with arbitrary mass and radius can easily be derived. In these simulations, the experimental equilibrium value of $R_\Omega$, before perturbation, is used as the effective radius of the polytrope, and thus the period is given by:
\begin{displaymath} Q=0.0703\Bigg(\frac{M_{\odot}R^{3}_{\Omega}}{MR^{3}_{\odot}}\Bigg)^{\frac{1}{2}} days \end{displaymath} (4.12)

In each of these tests, the originally static polytrope is instantaneously uniformly expanded radially by $20\%$, keeping the $K$ in the equation of state constant. This degree of radial expansion was chosen so that the results could be compared to those by Bate (1995). The subsequent behaviour depends on the number of particles used to model the polytrope, and more obviously, the amount of viscous damping present. The kinetic energy lost by damping would usually serve to heat the gas, but if this was implemented here, the final equilibrium configuration would not equal the original equilibrium configuration, as the value of $K$ would be larger, and thus the period of pulsation would vary over time.

Figure 4.8: Oscillation of polytropes. All energies are in units of the thermal equilibrium energy. Time is in units of the theoretical oscillation period. $n=1.5$. Constant viscosity: $\alpha=1$, $\beta=2$. Solid line, $531$ particles; Dotted line, $1067$ particles; Dashed line, $1807$ particles; Long-dashed line, $4279$ particles.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/oscildiffparta2.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 4.9: Oscillation of polytropes. All energies are in units of the thermal equilibrium energy. Time is in units of the theoretical oscillation period. $n=1.5$. Constant number of particles: $4279$. $\beta=2\alpha$. Solid line, $\alpha=4.0$; Dotted line, $\alpha=2.0$; Dashed line, $\alpha=1.0$; Long-dashed line, $\alpha=0.5$; Dot-dashed line, $\alpha=0.1$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/oscilpolytropesa1.eps,angle=0,width=1.0\linewidth}} \end{figure}

The results from these experiments are shown in Figure 4.8 and Figure 4.9. In Figure 4.8 the effect of varying the number of particles in the simulation is investigated, whilst Figure 4.9 demonstrates the effect of varying the viscosity parameters. The top frame in each of these figures displays the total kinetic energy as a function of time. The second frame displays the thermal energy as a function of time, which for an $n=1.5$ polytrope is equal to the total internal energy present. The third frame displays the ratio of thermal to gravitational potential energies as a function of time. The unit of time is the fundamental oscillation period, and the unit of energy is the initial thermal energy of the equilibrium configurations. In these experiments, typical values for the number of neighbours, $N_{N}=55$, the tolerance, $\epsilon_{r}=0.05$, and the opening angle, $\theta=0.57$, are used.

The graphs follow the trends demonstrated elsewhere in the literature, namely that decreasing the number of particles used to model the polytrope, effectively increases the amount of damping present, and that increasing the viscosity has the expected effect of decreasing the oscillation amplitudes. A possible explanation of the former behaviour was put forward by Bate (1995): For low resolution simulations, the particles are smoothed over a larger volume, and hence the viscous forces between particles act over larger distances and over larger velocity dispersions, increasing the bulk viscosity. However the claim by Bate that the oscillation periods are not dependent on the number of particles used to model the polytrope is clearly not the case here. The reason for this disagreement is probably related to how the effective radius is defined.

Although it may seem that reducing the number of particles produces a closer match to the theoretical period, it must be remembered that the theoretical period is only strictly valid for small perturbations from the equilibrium state. It can be seen from the first few cycles of the kinetic energy curves, that the polytropes spend significantly more time in an extended state, than in a compressed state, implying that a radial perturbation of $20\%$ does not correspond to the linear definition of a small perturbation. The theoretical period for this experiment could thus in fact be larger than the fundamental radial pulsation period, in which case, increasing the number of particles results in the usual expected increase in accuracy.


Radiation Transport

The radiation component of the code was developed in isolation from the gravitational and hydrodynamical components, and was subsequently incorporated exterior to these for reasons discussed in Section 3.2.11. Consequently, the testing during the development of this component was performed in isolation also. However, even if all components had been developed simultaneously, the differing nature and timescales present in the radiative and dynamical components, would have warranted separate testing. Certainly, had the dynamical components of the code been present during the development and testing of the radiation component, many of the errors detected would have gone unnoticed.

Although radiation transport is a relatively new area of research in SPH, it is fairly well developed for regular grid based methods. For this reason, standard tests of radiation transport exist (e.g., Stone et al., 1992). However, the tests developed and presented here are more relevant, as they have been designed to test the conditions present during typical simulations of the Capture Theory. The first two tests measure the `instantaneous' power radiated and received by a spherical body, whilst the second two tests examine the time evolution of the body.

Luminosity of an Isolated Sphere

This test is concerned with the rate at which energy is lost from the surface of an isolated uniform-density sphere of gas. The analytical solution for such a body is given by:

\begin{displaymath} L=4\pi\sigma R^{2}T^{4}\Big\{1-exp(-\tau)\Big\} \end{displaymath} (4.13)

where $\sigma$ is the Stefan-Boltzmann constant, $R$ is the radius of the sphere, $T$ is the temperature of the sphere, and $\tau$ is the optical depth across its diameter, given by:
\begin{displaymath} \tau=2R\rho\kappa \end{displaymath} (4.14)

where $\rho$ is the density, and $\kappa$ is the opacity of the gas. Since this test is completely scalable, the optical depth across the sphere is the only physical quantity of interest. The test is to compare the luminosity determined from experiment, with that given in Equation 4.13 for a range of optical depths. The experimental luminosity is calculated from the difference in energy of the particles before and after a small time-step. This procedure is preferred over the instantaneous rate of energy calculated on the grid during the first radiation time-step, as the scope of the code tested is greater. Obviously in theory, the two different methods should agree in the limit of a vanishingly small time-step. In practice they do so for sufficiently small time-steps, but as the time-step reduces to zero, the results diverge. This is because of the finite precision used to store the particle energies.

The experimental setup consists of $4279$ particles of equal temperature and mass, arranged on a displaced spherical grid. For radiation tests, the effective radius ($R_{e}$) of the body is taken as being equal to the distance to the furthest particle from the centre, plus twice its smoothing length. This is approximately $1.3$ times the distance to the furthest particle, for this number of particles and number of neighbours ($N_{N}=55$).

Given an optical depth across the sphere ($\tau$), the optical depth of leafs ($\tau_{g}$) of size $d$ are calculated using:

\begin{displaymath} \tau_{g}=\tau\frac{d}{2R_{e}} \end{displaymath} (4.15)

or alternatively, the opacity of the leafs which lead to the optical depth is calculated using:
\begin{displaymath} \kappa_{g}=\frac{\tau}{2\rho_{g}R_{e}} \end{displaymath} (4.16)

Figure 4.10 compares the experimental luminosity with that derived from theory, for a large range of optical depths, and for several different sized time-steps.

Figure 4.10: Comparison of experimental with theoretical luminosity. The solid curve shows the theoretical fraction of the black body luminosity as a function of optical depth. The dotted curves show the experimental luminosity divided by the theoretical luminosity for six different sized experimental time-steps.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/radtesta4.eps,angle=0,width=1.0\linewidth}} \end{figure}

In this figure, the solid curve corresponds to the theoretical fraction of the black body luminosity as a function of optical depth, and each of the dotted curves represent a different sized experimental time-step, the luminosity of these being expressed in units of the theoretical luminosity at that optical depth. The behaviour demonstrated is mostly consistent with that expected:

The reason why the large time-step curves deviate significantly from the desired value of unity, is that the body simply has sufficient time to cool. The finite precision storage of the energy is responsible for the corruption in the results at low optical depths for the smallest time-step. For large time-steps and optical depths above $1000$, the deviations are constant but considerably greater than for smaller optical depths. At these high optical depths, if the temperatures of the surface particles at the end of the time-step are inspected, they are found to be significantly higher than the temperatures required to produce this large discrepancy in the luminosity. This is because, at these high optical depths, radiation only escapes the system from surface leafs, which do cool to the required temperatures. Surface leafs are able to cool because adjacent leafs are unable to heat them efficiently, due to the large penetration depth attenuation factors involved, see Section 3.2.12.

Overall, the results are very encouraging, with most discrepancies being due to the deviation of the experimental data from the test scenario, i.e. if these experimental results were compared to the theoretical results of a cooling body, much closer correspondence would be obtained. Earlier radiation methods, and indeed earlier versions of this tree radiation method, could give deviations of many orders of magnitude in the high opacity regimes, even for small time-step experiments. It is unknown why the limiting experimental behaviour produces curves which lie a few tens of percent above unity around optical depths of $10$. Perhaps it is related to the radiation transfer within the tree being most efficient at these optical depths, due to the optical depths of the leafs being of order unity.

Although the experimental variables used in this test are not strictly necessary for future reproduction and verification of these results, they are given for completeness: $M=2\times 10^{27}\,kg$, $R=1.061\times 10^{11}\,m$, $R_{e}=1.375\times 10^{11}\,m$ and $T=426.2\,K$. The unit of time is a year, and the usual solar mix of gases with a ratio of specific heats equal to $\gamma=5/3$ was used. The pruning factor for tree construction was $P_{F}=2.0$, and the number of radiation time-steps was equal to $50$. Increasing either of these latter two tree variables had little effect on the results. These conditions are almost completely arbitrary, have little significance, and were used only for convenience.


Radiation Received from a Point Source

This experiment is particularly relevant to the Capture Theory, as it investigates how well the code can transport energy from a point source radiator, the primary star, to material modelled using SPH, the secondary star and protoplanets. The test is to compare the radiation absorbed in a spherical secondary body from a point source primary body, with that derived from theory. In similarity to the previous test, this is a scalable problem, where in theory, only the optical depth dependence need be investigated, the primary-secondary distance and radius of the secondary body being unimportant. However, the nature of the random walk sampling central to this radiation scheme, introduces additional considerations. In particular, if the secondary body is placed on an axis, it absorbs more radiation than if it is placed off axis at the same distance from the primary.

Although the random walk procedure has been rigorously tested during its development, and behaves as intended, it is apparent that the distribution of power packets emerging from a point source, does not have the desired spherical symmetry. The extent of this nonhomogeneity, which leads to the aforementioned erroneous absorption, is demonstrated in these tests.

In this experiment, and in general when point sources are present, a minimum spatial resolution is specified for the radiation tree structure. This is achieved by creating an additional set of virtual particles during the tree construction stage. These particles are arranged on a uniform grid encompassing the system, and have smoothing lengths of the order of the size of a grid cell, such that the leafs they create are not subsequently pruned. Typical high resolution experiments employ a $16^{3}$ grid for this purpose, whilst in low resolution experiments, an $8^{3}$ grid is used. This minimum resolution is intended to improve the random walk sampling procedure.

The experimental setup consists of $4279$ particles of equal temperature and mass, arranged on a displaced spherical grid. The effective radius, $R_{e}$, of the body is taken as being equal to the distance to the furthest particle from the centre, plus twice that particle's smoothing length. The distance between the centre of the primary body and the centre of the secondary body is denoted by $d$. A more relevant quantity of interest in these tests is the ratio of these distances:

\begin{displaymath} K=\frac{d}{R_{e}} \end{displaymath} (4.17)

Given that $K>1$, and that the optical depth, $\tau$, across the diameter of the secondary body is known, the theoretical fraction of radiation emanating from the primary that is absorbed in the secondary, can be calculated numerically, see Appendix B.2 for the derivation of this. In this experiment, a sufficiently small time-step is chosen to prevent significant heating or cooling occurring during the time-step. Although, this change in energy is not as critical as in the previous experiment, as the primary body is artificially kept at a constant temperature, and hence has constant luminosity throughout the time-step. The primary body emits between $10^{4}$ and $10^{5}$ power packets. Increasing the number emitted, increases the sampling accuracy, but also increases the computational run-time.

Two different values of $K$ are examined, $K=2$ and $K=5$. In each of these scenarios, the dependence on the angular position of the secondary body in the x-y plane, with respect to the primary body, is investigated. Figure 4.11 shows the results from these experiments.

Figure 4.11: Point source absorption test. (a) $K=2$. (b) $K=5$. Solid line, $\alpha=0^{\circ}$; Dotted line, $\alpha=15^{\circ}$; Dashed line, $\alpha=30^{\circ}$; Long-dashed line, $\alpha=45^{\circ}$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/rectesta1.eps,angle=0,width=1.0\linewidth}} \end{figure}

As the secondary body is rotated from the x-axis to the y-axis, the energy received reduces, reaching a maximum reduction to around a fifth of that received on an axis, at the midway point of $\alpha=45^{\circ}$. This is because the power packet density distribution on a spherical surface centred on a point source, is biased towards the axes. Indeed, if the secondary body is positioned on the corner of a cube centred on the primary, the radiation absorbed is reduced further still. This effect could conceivably be negated, perhaps by introducing additional weighting considerations into the constrained random walk algorithm. An approach such as this was only briefly investigated, as other more important discrepancies took precedence.

The behaviour of an individual curve is fairly consistent with that expected. At low optical depths, large fluctuations are present. This is because only a small proportion of the power packets incident on the surface of the secondary body are absorbed. Increasing the number of power packets emitted by the primary, reduces these sampling errors. The penetration depth attenuation of the energy received does not apply to power packets emitted from point sources, so that the constant limiting behaviour demonstrated by the curves at high optical depths is consistent with what should occur4.3.

Equilibrium Temperatures of Illuminated Bodies

This test is an extention to that presented in the previous section, as it is concerned with the equilibrium temperatures reached by bodies illuminated by a point source. The actual time evolution of the body, from a given starting temperature to the final equilibrium temperature, is investigated in the next section and is shown to be consistent with the theoretical scaling properties of the system.

Except where stated, the experimental setup is the same as in the previous section. If the illuminated sphere is assumed to be completely opaque, then the rate of heating of the sphere, due to a point source of luminosity $L$, is given by:

\begin{displaymath} \frac{dE_{BB}}{dt}=\frac{L}{2}\Bigg\{1-\bigg(1-\frac{1}{K^{2}}\bigg)^\frac{1}{2}\Bigg\} \end{displaymath} (4.18)

where $dE_{BB}/{dt}$ is the black body absorption rate, and $K$ is as defined in the previous section, see Appendix B.2 for the details. If the equilibrium temperature of the sphere is assumed to be the same throughout the body, then this temperature can be found by equating the rate of energy received in Equation 4.18, with the black body energy radiated by the sphere, which leads to:
\begin{displaymath} T_{eq}=\Bigg(\frac{L\Big\{1-\big(1-\frac{1}{K^{2}}\big)^\frac{1}{2}\Big\}}{8\pi\sigma R_{e}^{2}}\Bigg)^{\frac{1}{4}} \end{displaymath} (4.19)

Equation 4.19 also applies to semi-opaque bodies, due to the rate of energy received from the point source, being attenuated by the same amount as the rate of energy radiated by the sphere. Unless the sphere is rotating in such a way that each unit area of its surface shares an equal time facing the point source, the equilibrium temperature derived in Equation 4.19 will not be uniform throughout the sphere. In this experiment, the body is stationary at a given distance and orientation from the point source, so that the surface closest to the point source will be hotter than the surface furthest away. The magnitude of the difference in temperature depends weakly on how opaque the sphere is. This dependence is not investigated here, other than to say that it appears that the more opaque the body is, the larger the temperature difference. An optical depth across the sphere of $\tau=10$ was used in these experiments.

Figure 4.12: Equilibrium distribution of particle temperatures for $K=2$. A thin slice in the x-y plane of the $4279$ particles is shown. The particles are hued according to temperature, white being around $2.6$ times hotter than black. The radius of the circle associated with each particle is equal to half that particle's smoothing length.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/K2equila1.eps,angle=0,width=0.9\linewidth}} \end{figure}

In this test, the four experimental extremes of the last section are evolved in time until an equilibrium distribution of temperature is reached. Figure 4.12 shows the equilibrium distribution of temperature for $K=2$ and the sphere positioned on the x-axis. A thin slice in the x-y plane of the $4279$ particles is shown. The depth of hue is an indication to the temperature of the particles, white being around $2.6$ times hotter than black. The radius of the circle associated with each particle is equal to half that particle's smoothing length ($N_{N}=55$).

Figure 4.13: Point source equilibrium temperature test. Plotted on the y-axis is the temperatures of the SPH particles comprising the sphere, divided by the theoretical equilibrium temperature as determined from Equation 4.19. Along the x-axis is the distance of the particles from the point source, divided by the distance from the centre of the sphere to the point source. (a) $K=2$, $\alpha=0^{\circ}$, (b) $K=2$, $\alpha=45^{\circ}$, (c) $K=5$, $\alpha=0^{\circ}$, (d) $K=5$, $\alpha=45^{\circ}$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/equatesta2.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 4.13 shows the results for the experimental setup shown in Figure 4.12 and for the other three scenarios. Plotted on the y-axis is the temperatures of the SPH particles comprising the sphere, divided by the theoretical equilibrium temperature as determined from Equation 4.19. Along the x-axis is the distance of the particles from the point source, divided by the distance from the centre of the sphere to the point source.

In theory, the temperature of a particle should only depend on the distance it is from the point source, i.e. there should be no deviation in the temperature of particles at the same distance from the point source. The degree of scatter demonstrated in these graphs is deemed acceptable, and is probably due to the small difference in temperature between interior and surface particles. As seen in the previous section, the rate of energy absorbed by spheres placed between the x-axis and y-axis ( $\alpha=45^{\circ}$) is reduced significantly by almost an order of magnitude, compared to that of spheres placed on axis. Since the luminosity depends on the fourth power of the temperature, the reduction in the equilibrium temperatures demonstrated in Figure 4.13(b) and Figure 4.13(d) is insignificant in comparison to the deviations demonstrated in Figure 4.11. This is the reason why little effort was made in developing a solution to this discrepancy.

Radiation Transport in the Opaque Regime

This test is concerned with the rate at which energy is transported through an opaque medium. As a continuation to the experiments in the previous sections, the conditions of Figure 4.12 are used for this experiment. The test is to see if the body reaches an equilibrium temperature distribution, on a timescale which is consistent with that predicted by theory. In the opaque regimes present in the Capture Theory, it is particularly important that the rate at which energy is transported through the system is simulated correctly. For example, during the early stages of protoplanet formation, at a time when the planet is exhibiting damped oscillatory motion tending towards quasi-static collapse, if the energy generated in the centre during a compression cycle is transported to the surface too quickly, then this could lead to an erroneously large collapse rate. It is suspected that this is true of the work presented by Schofield (1979) on the contraction of newly formed protoplanets. In that work, a one dimensional finite difference code was used to solve the ensuing dynamical and radiative evolution of opaque protoplanets. An iterative integration scheme was used to solve the heat flow equations between shells of the protoplanet, and iteration ceased, when, to within a given tolerance, the temperature of shells between successive iterations were unchanging. An iterative approach like this can be applied successfully to linear heat flow problems, such as ordinary heat conduction, but is not necessarily automatically valid for non-linear heat flow problems, such as radiation diffusion.

Although the test is almost completely scalable, for clarity, a specific experimental setup is considered. The point source has a luminosity equal to that of the Sun ( $L_{\odot}=3.92\times10^{26}\,Watts$), the body has a mass equal to a Jupiter mass ( $M=2\times 10^{27}\,kg$), and a density equal to $\rho=4\times10^{-7}\,kg\,m^{-3}$. This results in an approximate effective radius of $R_{e}=0.92\,AU$ when $4279$ particles are used, and thus for $K=2$, a distance between centres equal to $1.84\,AU$. Initially the sphere is at a uniform temperature, a tenth that of the theoretical equilibrium temperature as given by Equation 4.19. For this setup, $T_{eq}=255\,K$ and thus the starting temperature is equal to $25.5\,K$. An optical depth across the sphere of $\tau=100$ is used, and this leads to a time-step equal to around half a year for the evolution to be adequately resolved in time4.4.

Figure 4.14 shows the distribution of temperature at the end of the first time-step.

Figure 4.14: Evolved temperature distribution after one time-step for the experimental setup shown in Figure 4.12. The particles are hued according to temperature, white being equal to 300K, and black equal to 25K. The radius of the circle associated with each particle is equal to half that particle's smoothing length. (a) All $4279$ particles projected onto the x-y plane, (b) A thin slice in the x-y plane, (b) All particles projected onto the y-z plane.
\begin{figure}\begin{center} \begin{tabular}{c@{\,\,\,}c@{\,\,\,}c} \epsfig{figu... ...idth=0.31\linewidth}\\ (a) & (b) & (c)\\ \end{tabular}\end{center}\end{figure}

As seen in Figure 4.14(c) the distribution is fairly symmetric about the line joining the centre of the bodies.

Figure 4.15: Temperature evolution during the course of the simulation. The time elapsed between frames is equal to 4 years (8 time-steps). A thin slice in the x-y plane of the $4279$ particles is shown. The particles are hued according to temperature, white being equal to $420\,K$ and black equal to $25\,K$.
\begin{figure}\begin{center} \begin{tabular}{c@{\,\,}c@{\,\,}c} \epsfig{figure=/... ...ps/12.eps,angle=0,width=0.31\linewidth}\\ \end{tabular}\end{center}\end{figure}

Figure 4.15 shows a sequence of frames demonstrating the evolution of the body during the course of the simulation. The time elapsed between frames is equal to $4$ years. As predicted by theory, a wave of energy, with a definite boundary between heated and unheated material, penetrates into the body with a finite speed. A complete description of the physics of this thermal, or `Marshak', wave is not relevant here, see (e.g., Zel'dovich and Raizer, 1966). What is relevant however, is the rate at which the wave penetrates into the material, and the governing scaling equations of the problem.

As demonstrated, the body reaches an equilibrium distribution of temperature in approximately $50$ years. The diffusion timescale for photons travelling through an optically thick medium is given by:

\begin{displaymath} t_{D}=\frac{x^{2}}{\lambda c} \end{displaymath} (4.20)

where $t_{D}$ is the time it takes photons to travel a distance $x$, $\lambda$ is the mean free path between collisions of the photons travelling through the material, and $c$ is the speed of light. This can also be expressed as:
\begin{displaymath} t_{D}=\frac{\tau x}{c} \end{displaymath} (4.21)

where $\tau$ is the optical depth of the material for photons travelling a distance $x$. If the radiation energy density present is much greater than the thermal energy density due to the translational motion of the molecules, the radiation will propagate through the body in a time given by Equation 4.20 or Equation 4.21. Indeed, a common experiment performed when testing radiation transport codes, is the evolution of a Gaussian pulse of radiation energy placed in a medium of known optical properties (e.g., Stone et al., 1992). If the radiation energy density is much greater than the thermal energy density, or only the radiation field is evolved in time, the pulse diffuses outwards on the diffusion timescale, maintaining its original Gaussian profile throughout. Its amplitude as a function of time being
\begin{displaymath} A(t)=\frac{A_{0}}{(t_{D}+1)^{3}} \end{displaymath} (4.22)

in three dimensions. The code presented here can reproduce this behaviour only approximately, and then only under a small range of conditions. This is because radiation is transported via the thermal energy field, so that if the energy of the radiation field is much larger than the energy of the thermal field, the implicit assumptions in several parts of the code become invalid. For example, because this code was not designed to handle these types of problems, it is possible for radiation to propagate faster than the speed of light. This problem could probably be eliminated by introducing a flux limiter, using analysis similar to that presented in Section 3.2.12. However, it is not until extremely large temperatures are reached that the radiation energy density becomes comparable to the thermal energy density, and thus this scenario is of no concern here.

The radiation energy per unit volume is defined as:

\begin{displaymath} U_{radiation}=\frac{4\sigma T^{4}}{c} \end{displaymath} (4.23)

where $\sigma$ is the Stefan-Boltzmann constant. The thermal energy per unit volume is defined as:
\begin{displaymath} U_{thermal}=\frac{3}{2}\frac{k\rho T}{\mu} \end{displaymath} (4.24)

where $k$ is Boltzmann's constant and $\mu$ is the mean molecular weight of the material. In this test, $U_{thermal}$ is approximately $5\times 10^{4}$ times larger than $U_{radiation}$, at $T=400\,K$. Since the radiation component is responsible for transporting energy through the system, and this is always coupled to the thermal component, it follows that if the thermal energy density is greater than the radiation energy density, the transport timescale is reduced by the ratio of these two:
\begin{displaymath} t_{transport}=\frac{3k\rho\tau x}{8\sigma\mu T^{3}} \end{displaymath} (4.25)

Substituting the values used in this experiment: $\mu=3.93\times10^{-27}\,kg$, $\rho=4\times10^{-7}\,kg\,m^{-3}$, $\tau=100$, $x=1.84\,AU$, and $T=400\,K$, results in a transport time approximately equal to $100$ years. Since this value is derived from an approximate scaling calculation, it is considered consistent with the approximate value of $50$ years obtained from experiment.

The scaling properties of Equation 4.25 were verified to be correct: Other than the optical depth, varying each of the quantities separately, or at the same time, and modifying the time-step appropriately, produces exactly the same numerical results to those shown in Figure 4.15. Aside from the obvious deviation in results as the optical depth reduces to zero, the correct limiting behaviour as the optical depth approaches infinity is not reproduced. The reason for this is related to the observation in the previous section, that as the optical depth increases, so does the difference in temperature between the front and back surfaces of the sphere. For example, an optical depth across the sphere of $\tau=100000$, results in an equilibrium distribution of temperatures similar to usual, but with the front surface at a temperature of around $1500\,K$, and the back surface at a temperature of around $300K$. If $1500\,K$ is used in Equation 4.25, along with the new optical depth, then a timescale consistent with experiment is obtained. The reason for an erroneously high equilibrium temperature distribution at high optical depths, is possibly related to the behaviour demonstrated in Figure 4.10. That is, for very high optical depths, the sphere can receive energy more efficiently than it can radiate it. In a real simulation, the effect is not so great, as the opacity of low density surface leafs is not artifically increased to maintain a constant optical depth per unit length across the sphere, using Equation 4.15 or Equation 4.16. If the real opacity of the material as a function of temperature is used, a temperature on the front surface of the sphere of around $700\,K$, and a temperature on the back surface of the sphere of around $30\,K$ is produced for these experimental conditions.

Front Page
Stephen Oxley 2002-01-19