Subsections
Conclusion
This thesis has been concerned with developing the tools necessary to simulate the
Capture Theory for the origin of planetary systems. The main numerical achievement of
this thesis is the development of a robust radiation transport algorithm, suitable
for use with SPH, particularly spatial tessellation tree implementations. With respect
to the Capture Theory, it has been demonstrated that the original premise of a
quasi-statically collapsing protostar-Sun interaction, is not conducive to planetary
formation, whereas a new type of encounter, ``tidal induced fragmentation followed by
capture'', involving a more freely collapsing protostar is.
The Capture Theory proposes that planetary material is captured from a nearby forming
star, and that the direct collapse of this material results in the formation of large
planets. The principal difference between the Capture Theory and its predecessors, is
that the roles of the interacting stars are reversed. Early tidal theories proposed
that the planetary material was torn from the Sun during a close encounter with a more
massive star. It was subsequently shown that the direct collapse of solar material to
form planets was not possible, Spitzer (1939), and that insufficient angular
momentum would be imparted to the material to account for that contained in the
orbital motion of the planets, Russell (1935). In addition, the spectrographic
evidence indicates a cold origin of planets.
A major problem which nebula theories of planetary formation must address, is the
distribution of angular momentum in the Solar System: The Sun with 
of the
mass, has only 
of the total angular momentum contained in its spin, the
remaining 
residing in the orbits of the planets. Nebula theories which
address this problem directly, by assuming that the Sun captured the nebula material,
or that the density distribution in the forming nebula is centrally peaked, are unable
to provide planets as a result of the direct collapse of nebula material, due to the
large tidal forces from the central body preventing collapse, Jeans (1919).
Numerical examples of tidal disruption and gravitational fragmentation, in the context
of early theories, but equally applicable to the Capture Theory, were presented in
Chapter 1.
Early simulations of the Capture Theory demonstrated that both the proposed capture
of material, and the subsequent condensations into protoplanets was possible,
Woolfson (1964); Dormand and Woolfson (1971). However, the assumptions that were necessary in order to
numerically contract the problem sufficiently for simulation purposes, were
subsequently shown, in Chapter 5 of this thesis, not to be valid for
the type of protostars which had originally been envisaged. A review of the Capture
Theory, including the features of the Solar System that all planetary formation
theories should address, almost all of which are explicable in the context of the
Capture Theory, concluded Chapter 1.
The SPH technique was primarily chosen because of its ability to simulate large
density gradients and voids, both of which are present at various stages in the
Capture Theory. SPH has matured considerably since its introduction by
Gingold and Monaghan (1977) and Lucy (1977), despite it still being referred to as a new
method. Unlike the first numerical methods devised to simulate the Capture Theory, SPH
does not require the introduction of pseudo-forces, and thus the physical processes
are represented more accurately than previously possible.
One outstanding question in SPH, is with regards to how the smoothing lengths of
particles should be defined. The consensus is tending towards a definition, whereby a
constant number of a particle's neighbours are contained in twice its smoothing length
at all times. The computational expense in achieving this, usually results in most
implementations obeying this definition only approximately. When this is so, the
smoothing length ceases to be a function of position only, and consequently an
additional time dependence is introduced into the equations of motion. A fast method
for determining the distance to the 
neighbour of a particle was
developed. This involves using the previous calculation of the smoothing length to
collect 
particles via one tree traversal per particle. A reduced
quicksort algorithm sorts only the 
particle into the correct position in
the list. This procedure typically consumes 
of the total computational run-time
of the force calculation routine.
In the simulation domain of the Capture Theory, the number of degrees of freedom of a
hydrogen molecule is a function of temperature. As the temperature rises between
and 
, the number of degrees of freedom increases. With relation to the
early stages of planetary formation, this increase could potentially absorb
substantial quantities of energy, which would have otherwise contributed to supporting
the protoplanet hydrodynamically, via translational (thermal) molecular motion, i.e. a
filaments ability to fragment, and the rate of collapse of protoplanets, could be
significantly enhanced. A novel approach of including this temperature dependence into
the usual SPH equations was developed. This method also allows for the energy absorbed
during changes of state, and the resulting change in mean molecular mass, if any, to
be included, and is particularly useful in accounting for the energy of dissociation
of hydrogen molecules, and the consequent change of state, incurred at
.
Most implementations of SPH employ a low order integration scheme, with time-step set
by local physical attributes of the system, e.g. the Courant time. These methods are
derived from their grid based counterparts. However, in particle based methods, the
opportunity of controlling the time-step implicitly is available. Most integration
schemes which employ automatic time-step control are used to solve high accuracy
problems, and are therefore usually of high order. Those which are of low order are
also in general concerned with high accuracy problems, and tend only to become
efficient when relatively small time-steps are used. These methods are therefore not
applicable to SPH, where relatively large time-steps, and large local errors are
acceptable, due to the ensemble of particles being of more interest than the
individual particles. This led to the development of a new embedded second order
integration scheme, with automatic time-step control suitable for large time-step SPH
applications. This method is based on a combination of the forward Euler scheme and
the improved Euler scheme, and when auxiliary equations are not solved, FSAL (first
step as last) can be implemented to reduce the number of force calculations per
time-step from two to one.
When particle methods are used to simulate problems containing long-range forces,
the number of calculations, and hence the computational run-time, scales as 
.
If the particles are grouped together in a hierarchical tree structure,
Barnes and Hut (1986), the computational run-time can be reduced to scale as 
log. In
hierarchical tree based methods, the forces between neighbouring particles are still
calculated directly, but particles which are sufficiently far away are grouped
together, and an average force applied from their centre of mass. Spatial grouping of
particles also results in an efficient method of locating a particle's neighbours,
which are required for the smoothing length, density, and hydrodynamical force
calculations. Other than the unique implementation of this method, the pseudo-code of
which is given in Appendix A, little advancement was made in this
area during the course of this thesis.
The Barnes-Hut type tree referred to above was implemented in order to reduce the
computational run-time of a simulation. As a result of its inclusion, it was realised
that the same principles of space division and data storage could be used in a new
radiation transport algorithm. Prior to this realisation, an attempt at including
radiation transport into SPH, yet to be satisfactorily achieved elsewhere, had focused
on an method, in which energy was transported directly between particles,
attenuated by a uniform interparticle grid surrounding the simulation domain. This
method had several serious drawbacks. The first and most obvious was that the
computational run-time scaled as , and was therefore unsuitable for tree based
implementations of SPH. The second drawback was that the method involved a uniform
grid, and was therefore unsuitable for solving problems with large density gradients
and voids, exactly the sort of problems which SPH is well suited to. Disregarding the
previous two drawbacks, the final and most important drawback was that the method
only performed well in the low opacity regime: erroneous absorption of energy by the
grid, led to an incorrect limiting behaviour in opaque regimes.
Following the insights into radiation transport, in the context of SPH, gained from
the first `preliminary' investigation, a more robust transport scheme was developed.
This scheme is based on a Barnes-Hut tree, and therefore resolves space according to
particle density. In order to reduce the unnecessary resolution incurred when
particles become uncharacteristically close, a method of closing, ``tree pruning'',
leafs that are determined to resolve space too finely is employed. The particle
properties are spread onto the tree using a modified interpolation procedure, and the
optical properties of leafs, and subsequently all levels of the tree, are calculated.
A Monte-Carlo type approach is adopted for transporting radiation energy, whereby each
leaf emits packets of energy per unit time, ``power packets'', on constrained random
walks throughout the tree. The power packets are either absorbed in other leafs, or
escape the system. The change in energy of the leafs during the time-step is then
interpolated back to the particles.
This scheme has several advantages over direct methods of radiation transport: The
spatial resolution, and hence the computational run-time, can be controlled during
tree pruning; the number of power packets emitted by the leafs, can be used to control
the accuracy and also the computational run-time; point sources can be included in the
simulation domain, without degrading the resolution afforded to the SPH particles. The
method performs equally well in all opacity regimes. There is however room for
improvement. For example, currently the dynamic time-step is divided into 
radiation time-steps of equal length, and a simple forward Euler integration scheme is
used to calculate the change in energy of leafs during each of these. A procedure for
determining the optimum number of radiation time-steps, perhaps not even of equal
size, along with the inclusion of a higher order, perhaps semi-implicit scheme, would
increase the physical accuracy, without increasing the computational run-time.
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. The gravitational force
calculation was tested in isolation to see how well it could reproduce the
theoretical behaviour of a uniform-density spherical distribution of material,
collapsing under self gravity alone. Several types of simulations were investigated
in order to gain insight into the SPH method. The deviations from the theoretical
results in each of these simulations could be explained in terms of `features' of
SPH. For example, when softened gravitational forces were employed, the additional
softening at the surface of the sphere resulted in a shell of reduced density
propagating radially inwards.
Once a component of the code has been tested, and has proven to operate as expected,
it is reasonable to include and test other parts of the code alongside it, rather than
testing them in isolation first. A good test of both gravitational and hydrodynamic
forces is that of the equilibrium structure, and oscillation of polytropic spheres of
gas - polytropes. When comparing the equilibrium radial profile against theory, the
experimental radius of the polytrope is required. The distance to the furthest
particle could be used for this purpose, or even the distance to that particle plus
twice its smoothing length, but both of these give unsatisfactory results. A novel
approach was investigated, in which the radius of the polytrope is computed from the
experimental calculation of the gravitational potential energy of the polytrope. With
the implementation of this procedure, the experimental and theoretical profiles
coincided almost exactly. The periods of oscillating polytropes were determined to be
slightly larger than those predicted by theory. The large amplitude of oscillation was
concluded to be the cause of this difference.
Each of the four radiation tests were designed to test the conditions present in the
Capture Theory. The first test was concerned with the luminosity of an isolated
uniform-density sphere of gas. For small time-steps, the results from this experiment
were within of the theoretical values for all opacity regimes. In the large
opacity regime, larger time-steps resulted in the luminosity reducing significantly,
due to the surface leafs cooling. It was concluded that if the theoretical results
were calculated to take account of this cooling, then the experimental results would
match much closer.
The second test was concerned with the rate of heating in a sphere of gas from a point
source radiator. As expected, the sampling errors were large for small opacities, due
to a relatively small number of power packets emitted by the point source being
absorbed in the sphere. Attention was drawn to the non spatially spherical
distribution of power packets, an issue to be addressed in future work. This
non-uniformity leads to the absorption being between a factor of 
and 
different to that predicted by theory. However, even these `large' values result in
equilibrium temperatures no more than 
in error, which is still relatively small
where radiation transport is concerned.
The third test was a combination of the first two, and was concerned with the
equilibrium temperature reached by the sphere. In this test, the same surface of the
sphere always faced the point source, resulting in a temperature profile across it
developing. The scatter in these profiles were deemed acceptable, and the mean
temperatures were observed to be comparable with those predicted from theory.
The final test investigated the rate at which the sphere reached an equilibrium
temperature profile. The sphere was initially at a temperature much cooler than its
mean equilibrium temperature. Due to the sphere being opaque, a thermal wave of
energy, originating from the point source radiator, propagated through it with finite
velocity, and a definite boundary between heated and unheated material was clearly
visible. It was determined that the velocity of this wave was of the order predicted
by theory.
The capture and subsequent fragmentation mechanisms central to the Capture Theory,
were investigated in Chapter 5. This chapter began with an
investigation into the contraction rate of quasi-statically collapsing bodies. This
investigation concluded that it was not until the mass of protoplanets exceeded around
Jupiter masses, that temperatures in the central regions could exceed
within 
years. The significance of this conclusion is that newly formed Jupiter
mass protoplanets, almost certainly remain in an extended state for many orbits.
Previous research, Schofield and Woolfson (1982a), had concluded that a fast hydrodynamical
collapse of such bodies could occur in a relatively short period of time,
years, such that the core was stable to disruption at perihelion.
However, the formation of a core stable to disruption is not so critical for the new
type of encounter, as the perihelion distance of the orbits is larger than previously,
and the time the protoplanets spend at aphelion longer.
During the investigation into quasi-statically collapsing bodies, it was concluded
that protostars could not exist in an extended state of contraction, which was the
underlying assumption of the original type of encounters. Even if it is assumed that
such bodies do exist, it was discussed why material drawn from them is not conducive
to planetary formation, the main reason being that the material never leaves the Roche
limit of the Sun. A new type of encounter, ``tidal induced fragmentation followed by
capture'', which is conducive to planetary formation, was suggested. In this new type
of encounter, the whole of the approaching protostar is tidally drawn into a long thin
filament, with properties such that it is Jeans critical, but within the Roche limit,
at perihelion. As the filament leaves perihelion, its orbit takes it outside of the
Roche limit, and it is able to fragment into several protoplanets, around half of
which are captured. Although this simulation domain has yet to be thoroughly
investigated, it was demonstrated in an example encounter, that heating of the
filament by solar radiation was important. The protoplanets formed in this sample
simulation were of mass around times that of Jupiter, such that their cores
reached the hydrogen dissociation point before first perihelion with the Sun, and
would have therefore been quite stable to tidal disruption.
There are many paths in which future work may take. Clearly a complete investigation
into the new type of capture mechanism is one. Following the development of this code,
and further increases in computational power, other areas of the Capture Theory are
now open to investigation. For example, satellite formation is an area which needs
reassessing in the light of current knowledge. In order that later stages of planetary
formation, and consequent satellite formation may be investigated, a more robust
treatment of artificial viscosity would be necessary, due to rotation, and
consequently angular momentum transport being of importance. Moving to larger scales,
a greater knowledge of the star formation process would remove some of the
arbitrarily chosen starting conditions from the new type of capture encounters.
Simulations connecting the star formation process to the capture process, and
ultimately to the planet and satellite formation process, will soon be possible. The
most useful improvement to the radiation transport algorithm, other than those already
discussed, would be the inclusion of frequency dependent opacities. This is particularly
relevant to the Capture Theory, due to the large difference in temperatures present
when solar radiation is included.
In its modified form, the Capture Theory provides a plausible mechanism for the origin
of planetary systems.
Front Page
Stephen Oxley
2002-01-19