Subsections


Additional Computational Techniques

In the previous chapter, the basic equations that govern the evolution of a self-gravitating gas were presented. This chapter is concerned with the additional methods that are required, in order that the Capture Theory may be modelled with as high a spatial resolution and physical accuracy as possible. The addition of radiative energy transfer, Section 3.2, is developed in order to improve physical accuracy, whilst a spatial tessellation tree, Section 3.1, is included to reduce the computational run-time.


Tree Gravity

Principles and Methods

In this section, the general procedure for implementing a hierarchical tree data structure (tree-code) is presented in the context of the gravitational force calculation. The tree-code is subsequently developed in Section 3.2 for the use of radiative energy transfer. The reason for it being referred to as a tree-code, is that the data structure is conceptually similar in form to a tree.

Figure 3.1: An example of a tree data structure corresponding to $64$ particles arranged on a two dimensional $8$ by $8$ grid, see Section 3.1.3.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/datatreea2.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 3.1 demonstrates the structure present in a typical tree-code. This structure is explored below, and explained in more detail in Section 3.1.2, where a more realistic example is given.

An alternative method to that given in the previous chapter for calculating the gravitational force between particles, is clearly necessary in order that the computational run-time may be reduced. This is because the long-range nature of the gravitational force, requires that every particle in the system interacts with every other particle. This leads to the computational run-time scaling as $O(N^{2})$, where the $O$ corresponds to an approximately constant computational overhead. The best any scheme can hope for is a computational run-time that scales as $O(N)$. Even if this could be exceeded, not much would be gained, as the hydrodynamical component of the force calculation scales as $O(N_{N}N)$, which reduces to $O(N)$.

The choice of alternative method is restrained by certain factors. The main factor is that the benefits associated with using the SPH method should not be undermined. Specifically, the system should still be Lagrangian, and the spatial resolution not compromised. One method could be to extrapolate the mass onto a uniform grid, and solve for the potential using a Fast Fourier Transform algorithm. The computational run-time could then easily be controlled by the number of cells in the grid. Unfortunately this method as given, could severely compromise the resolution in a highly clustered simulation. A way of avoiding this could be to recognise where the resolution is compromised, and subdivide the offending cells into a similar grid structure, but on a smaller scale. One such method that uses this approach with SPH is being developed by Couchman et al. (1995) . Applying this approach to the example of Section 2.1, results in Figure 3.2.

Figure 3.2: An example of grid based resolution refinement.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/ppma1.eps,angle=0,width=1.0\linewidth}} \end{figure}

Here the centres of both spheres have been re-sampled with a finer grid, resulting in an increased resolution in those regions. In general, the depth of subdivision and the number of regions subdivided is greater than shown here.

The method of dividing space according to resolution requirements is central to most efficient gravitational potential solvers. Usually when such an approach is adopted, the computational run-time is reduced at the expense of the accuracy in the computed force. In general, particles still incur a direct summation with their closest neighbours, but particles far away are grouped together, and an averaged force applied. The motivation behind this can be understood by the following example:

Consider the forces between a distant galaxy and our own. In theory the gravitational force between every single subatomic particle in both galaxies should be computed. Obviously this is an overwhelming task in practice, and in general, when the galaxies are well separated, it is sufficient enough to calculate just the force between the centre of masses. Again, an overwhelming task is confronted when calculating the relative motion of the stars within the galaxies, even if the stars are considered to be single bodies. In this case, grouping the stars into clusters, and calculating the force between the centre of mass of those, greatly reduces the computational load. Further grouping of the stars into subclusters, and similar subcluster-subcluster force calculations within a cluster, can be performed if necessary. Finally, the force between stars in the same (sub)cluster can be calculated directly. A simplified schematic overview of this method is shown in Figure 3.3.

Figure 3.3: Schematic overview of a cluster-cluster force calculation.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/cellcella3.eps,angle=0,width=1.0\linewidth}} \end{figure}

The diagram shows the force evaluations required to follow the motion of the two particles on the bottom left of the diagram. The double arrowed lines represent centre of mass force calculations between clusters at the same level in the tree hierarchy.

The method outlined above belongs to the cluster-cluster family of tree-codes. An overview of the procedures used in most implementations of this type of tree-code is now given. It may be useful to refer to Figure 3.3.

Cluster-cluster interactions are computed prior to calculating the force on the particles. This involves interacting every cluster with every other cluster at the same level in the hierarchy, and belonging to the same region of space, i.e. clusters with the same `parent' cluster. When calculating these interactions, extra accuracy is obtained by considering higher moments of the mass distribution than just the total enclosed mass and centre of mass. This is normally necessary in a cluster-cluster tree-code, and although not necessary, it is common practice in a particle-cluster tree-code, explained later. The potential field in a cluster, calculated from all other clusters, can then be determined; the higher the moments used, the more accurate the field. Once the cluster-cluster interactions have been performed, the force on a single particle can be calculated: First a direct force calculation with all particles belonging to the same cluster is carried out. From the previous cluster-cluster interactions a potential field gradient at the position of the particle is then estimated and used to approximate the force from the remaining particles.

An alternative approach to a cluster-cluster tree-code is a particle-cluster tree-code. The difference in approach between the cluster-cluster tree-code, and the particle-cluster tree-code is demonstrated in Figure 3.4.

Figure 3.4: Schematic overview of a particle-cluster force calculation.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/pcella3.eps,angle=0,width=1.0\linewidth}} \end{figure}

The diagram shows the force evaluations required to follow the motion of the same two sample particles as examined previously. The details of the actual particle-cluster tree-code used in this thesis, are presented in Section 3.1.2. An overview of the method is given now:

As the name implies, particles not only interact directly with neighbouring particles in the same cluster, but also directly with other clusters of particles. The centre of mass in every cluster is computed prior to calculating the force on the particles. As before, higher multipole moments describe the mass distribution more accurately. Once the mass distribution of every cluster in the tree hierarchy has been computed, the force on a single particle can be calculated: As before, a direct force calculation with all particles belonging to the same cluster is performed. The particle then interacts with the clusters. As seen in Figure 3.4, the further away a particle is from a collection of clusters, the higher up in the cluster hierarchy the particle interacts with.

In the past, cluster-cluster type tree-codes have been claimed to scale as $O(N)$. Recently this claim has been disputed, and a more conservative $O(N$log$N)$ estimation has been demonstrated in practice (Capuzzo-Dolcetta and Miocchi, 1998). Figure 3.3 is a very simple example of a cluster-cluster simulation, containing just two layers of clusters, not counting the particle layer. In general, many layers of clusters are present, so that the computational overhead hidden in the $O$ can be very large. Cluster-cluster tree-codes are generally more accurate than particle-cluster codes. Normally this is a good reason for using them, however in astrophysical simulations using SPH, the additional accuracy is often not required, due to reasons discussed in Section 2.11.

Particle-cluster tree-codes also scale as $O(N$log$N)$, where the $O$ is usually smaller than in cluster-cluster tree-codes. Capuzzo-Dolcetta and Miocchi concluded that their optimised implementation of the Barnes-Hut tree-code (a particle-cluster tree-code, see Section 3.1.2.) became more efficient than direct summing, for $N\approx 10^{3}$, and was approximately $5$ times faster than their optimised version of the Fast Multipole Method (a cluster-cluster tree-code) in the simulation domain $N\le2\times10^{5}$. It was not until the number of particles exceeded $\approx 10^{4}$, that the FMM outperformed direct summing. Both tree-codes were constrained to produce an error of no greater than $0.1\%$ in the force on a particle.

The method used in this thesis is a particle-cluster tree-code. This was chosen because of the reasons above, but also because this type of tree-code is simpler to understand and program. Having implemented the tree-code, it was realised that the same principles of space division and data storage, could be used for a new radiative energy transfer technique, see Section 3.2. Since other methods can be more spatially abstract in nature, it is uncertain whether it would be possible to implement the new radiation transport scheme, without explicitly introducing a Barnes-Hut type tree structure.


Tree implementation

The tree-code used in this thesis is most similar in form to that proposed by Barnes and Hut (1986). This is a top-down particle-cluster tree-code, where the construction of the tree is performed by recursive subdivision of the simulated region. Conversely, a bottom-up tree-code begins with the particles, and recursively groups them together into larger clusters to form a tree structure. The Barnes-Hut algorithm is split into three stages:
  1. Tree construction.
  2. Centre of mass calculation of all clusters via $1$ tree traversal.
  3. Force calculation on the particles via $N$ tree traversals.
Each step in the algorithm is now examined:


Tree Construction

For the time being it is simpler to consider only two dimensions. The methods are readily extended into three dimensions. The purpose of the tree construction stage is to subdivide the simulated region into squares, each containing no more than one particle. This is achieved with the following algorithm:
  1. Define a square region of space that just encompasses all of the particles.
  2. Insert a new particle into the largest square encompassing all of the particles.
  3. If the current square is subdivided, move the current particle down to the appropriate sub-square and goto 3.
  4. If the current square contains only the current particle goto 2.
  5. Divide the current square into four equal sub-squares, and descend both the current particle and the original particle into the appropriate sub-squares.
  6. If both particles still reside in the same square goto 5, otherwise goto 2.
Normally a tree algorithm such as this is implemented recursively, but in a programming language where this is not possible, an iterative implementation is an acceptable alternative. Obviously, if programming style or efficiency is unimportant, it can be used exactly as given. Since FORTRAN77 is the programming language used throughout this thesis, iteration is used. The pseudo-code for inserting a particle into the tree is given in Appendix A.1.

Figure 3.5 shows an example quad-tree after construction. Figure 3.5(b) is the tree representation of the spatial distribution, and was formed by starting at the top-left square in Figure 3.5(a), and rotating clockwise through the hierarchy. The dotted lines in Figure 3.5(b) correspond to empty squares, and are subsequently removed to reduce storage space and improve tree-traversal efficiency.

Figure 3.5: Representations of a quad-tree. (a) Spatial, (b) Tree.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/quadtreea4.eps,angle=0,width=1.0\linewidth}} \end{figure}

Applying this procedure to the example of Section 2.1, results in Figure 3.6.

Figure 3.6: Spatial tree structure for the example in Section 2.1.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/treegravexa1.eps,angle=0,width=0.6\linewidth}} \end{figure}


Tree Storage in Memory

In this implementation of the tree-code, the squares of the tree are arranged in memory in the same order in which they were created. When a square is subdivided, the created sub-squares are positioned at the end of the tree-list, one after another. The square that was subdivided contains a link forward to its first sub-square, and all sub-squares contain a link back to their `parent' square. The empty squares are removed by deleting their entry in the tree list, and shuffling every subsequent entry up a place, making sure that all links are updated.

The resultant data structure is arranged almost completely randomly in memory, due to the effectively random location of the particles in space. A more logical approach could greatly improve tree access efficiency. Salmon et al. (1994) have progressed some way to developing such a structure. In their implementation of the Barnes-Hut tree, the distribution of the tree in memory is closely related to the spatial distribution of the particles. This has advantages in parallel processor implementations, as a portion of space, along with the corresponding tree fragment, can be allocated to an individual processor. Using the data structure implemented here, the complete tree structure would have to be available to each processor. Their method has advantages in single processor machines also, as the tree data structure results in mostly linear memory access during tree traversals (the most heavily used tree function, see later.)

In the radiative transfer of energy, Section 3.2, the tree is most heavily used during the random walks of power packets (explained in that section) through the tree, resulting in the data structure of Salmon et al. being less attractive. However, that is not to say that a more logical memory arrangement is not possible. In hindsight, since the tree construction routine takes practically zero computational run-time, almost any additional routine to arrange the tree in memory more logically would be beneficial. In fact, the more logical approach of Section 3.2.6 for removing unnecessary squares, results in almost linear memory access for tree traversals.


Centre of Mass Calculation

This is achieved by traversing the tree once, progressively accumulating the mass and centre of mass at each level in the tree. A tree traversal involves moving from left to right through the tree structure, and is shown schematically in Figure 3.7 for the previous example.

Figure 3.7: Tree traversal for the example in Figure 3.5.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/quadtrava2.eps,angle=0,width=1.0\linewidth}} \end{figure}

The pseudo-code for calculating the mass and centre of mass of every square in the tree is given in Appendix A.2.

Force Calculation

Assuming a purely gravitational point mass system for the moment, the force on a single particle is found via one tree traversal. Whilst traversing the tree, the decision to open a square and consider the force from the sub-squares it contains, rather than use the current mass and centre of mass approximation, is decided upon by a MAC (Multipole Acceptability Criteria). Several MACs have been proposed in the past, see Salmon and Warren (1994) for a review of the most commonly used ones. In this implementation of the Barnes-Hut tree, only the zeroth and first order moments are calculated, so that the aforementioned multipole, refers to the centre of mass.

The MAC used in this work is the original one proposed by Barnes and Hut (referred to as the BH-MAC hereafter). It asserts that the current square can only be used, if the ratio of the size of the square to the distance between the particle and the centre of mass in the square, is less than a tunable parameter, $\theta$, called the opening angle. This leads to approximately the same percentage error in the force acting on a particle, regardless of the distance to the region it is interacting with. In the example of Figure 3.8 this criteria corresponds to $s/r<\theta$.

Figure 3.8: Testing the multipole acceptability criteria: The BH-MAC.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/maca1.eps,angle=0,width=1.0\linewidth}} \end{figure}

For three dimensional astrophysical simulations, fairly large values of the opening angle are commonly used, $0.7\le\theta\le 1.0$. However, Salmon and Warren warn that using values as large as this can potentially cause problems. They suggest a more conservative value of $\theta<1/\sqrt{3}$. This value of the opening angle corresponds to the limit where a particle could potentially fail to open the cube it belongs to. Figure 3.9(a) shows the interactions required to calculate the force on a sample particle, highlighted in Figure 3.9(b), for the two dimensional equivalent limit, $\theta=1/\sqrt{2}$.

Figure 3.9: Gravitational force interactions required for a sample particle with $\theta=1/\sqrt{2}$, (a) Spatial representation, (b) Tree representation.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/quadgrava1.eps,angle=0,width=1.0\linewidth}} \end{figure}


TreeSPH

In the explanation of tree-codes just given, only point mass gravitational forces were considered. The tree method also has advantages to offer the hydrodynamical force calculation. One such advantage is the efficient determination of the neighbour list of a particle via a tree traversal. The algorithm used to achieve this involves traversing the tree, only opening squares that the particle overlaps into. For efficiency, the $2h$ sphere of the particle is replaced by a box of sides $4h$ encompassing the sphere, and it is this which is tested for overlap. Once the particle layer is reached, the $2h$ sphere is used.

The approach of collecting the neighbours of a particle, results in an array of size $N$ by $N_{N}$ integers. Since $N_{N}$ can potentially be rather large3.1, most memory conservative codes prefer to re-find the neighbours of a particle for subsequent use, rather than store them. This is the method used here, but in hindsight, storage may have been possible, since an array of integers approximately that size is used during the transfer of radiative energy, see Section 3.2.

The complete tree-code algorithm for calculating the force is now given:

  1. Tree construction.
  2. Centre of mass calculation via $1$ tree traversal.
  3. Smoothing length calculation via $N$ tree traversals.
  4. Density calculation via $N$ tree traversals.
  5. Both gravitational and hydrodynamical force calculation in the same routine via $N$ tree traversals.
Stage 1 and Stage 2 have already been discussed, and pseudo-code given in Appendix A.1 and Appendix A.2 respectively. The smoothing length calculation is given in Appendix A.3 with respect to variable particle mass, where it is proposed that a constant enclosed mass is appropriate, rather than a constant number of neighbours. When the particles have equal mass, this algorithm reduces to a constant number of neighbours search. The density calculation, Stage 4, is given in Appendix A.4.

A few points need to be made concerning Stage 3 and Stage 4. Whilst traversing the tree to collect the required particles, the normal MAC is unnecessary. Instead, the current particle is enclosed in a box of sides $4h$, and squares in the tree are only opened if this box overlaps into them. This additional SPH-MAC is combined with the usual BH-MAC in Stage 5, where both the gravitational force and hydrodynamical force are calculated via one tree traversal per particle. When both MACs are used, the SPH-MAC takes precedence, i.e. a square must be opened if the $4h$ box surrounding the particle overlaps into it.

When using a tree-code for calculating particle interactions, symmetry in the particle quantities must be considered more carefully. In Stage 4 and Stage 5, whilst traversing the tree, quantities are calculated only at the position of the current particle. This is unlike an $N^{2}$ code, where quantities at both particles are incremented via a single interaction. The consequence of this is that ensuring symmetrisation in the density and force calculation is not so straightforward. Essentially the problem arises from the spread nature of the particles. In an $N^{2}$ code, a particle interacts with another particle if $r_{ij}<h_{ij}$. In a tree-code, particle $j$ may not interact with particle $i$, even though $r_{ij}<h_{ij}$. This would happen when $h_{j}$ is sufficiently smaller than $h_{i}$, and is unavoidable given the current MACs.

Symmetry is ensured by keeping track of particle interactions. Given the above hydrodynamic case of $i$ interacting with $j$, but not vice versa, symmetry is ensured by also incrementing the quantity at $j$ due to $i$, even though it is particle $i$ that is currently under consideration. This is fairly straightforward to implement for the hydrodynamic quantities of density, pressure force, and internal energy. The gravitational interaction is slightly more complicated, as the point-point force that will be calculated on $j$ due to $i$, when $j$ is under consideration, must be counteracted whilst particle $i$ is still under consideration.

The pseudo-code for Stage 5 is given in Appendix A.5, where the above methods concerning symmetrisation are presented more completely. The pseudo-code is for the general case of differing gravitational ($\epsilon$) and hydrodynamical ($h$) smoothing lengths.


Tree Radiation

In the previous section, a procedure for dividing space such that the spatial resolution closely followed the distribution of the particles was presented. This led to an efficient method for calculating the forces between particles. In this section, the aforementioned tree-code is extended for the use of radiative energy transfer.

Motivation

In the simulation domain explored in this thesis, the ability of the gas to radiate energy could potentially play an important role. This can be demonstrated by considering the timescales and optical depths present in a typical simulation. There are essentially three stages leading to planetary formation in the Capture Theory: protostar disruption, filament capture and condensation, and planetary contraction. Each of these stages can be classed into one of three regimes according to the opacity of the material present: transparent, semi-opaque, or opaque. Each of the stages can also be broadly classed according to the evolutionary timescale: dynamical, i.e. of the order of the free-fall time, or radiative, i.e. of the order of the time to lose energy via radiation.

Although this classification is fairly straightforward at the end of a simulation, it is difficult in this model to classify each stage a priori to simulation. This is because different simulation parameters produce differently classed stages, i.e. by varying the initial conditions, a stage determined to be opaque and controlled by radiative processes, could subsequently be determined to be opaque, but evolving on a dynamical timescale. As a relevant example of this, consider the planetary contraction stage. In this stage, the timescale for collapse is usually governed by how effectively the planet can radiate energy, the radiative timescale. However, given conditions which produce a sufficiently cool protoplanet such that the thermal energy is much less than the gravitational potential energy, the timescale for contraction will be of the order of the free-fall time.

This lack of a priori knowledge is one reason which prevents simplified radiative schemes from being applied. One such scheme which relies on the above classification was outlined in Section 2.9. In that scheme, the knowledge of the optical depth as a function of density is required. In implementations of SPH which use schemes similar to this, the gas is assumed to reside in either completely transparent or completely opaque regimes. It is also assumed that the dynamical timescale is much shorter than the radiative timescale in the opaque regime, and that the radiative timescale is much shorter than the dynamical timescale in the transparent regime. The latter assumptions are required in order that an explicit radiative energy transfer mechanism not be necessary, i.e. in the transparent regime, a compressed gas has sufficient time to cool to the same temperature as surrounding gas, and in the opaque regime, the change in energy due to radiative processes is much less than that due to dynamical ones, and thus can be ignored. In these methods, only the consequent behaviour of the gas due to it residing in a particular regime is important, the actual transfer of energy being unimportant. As outlined in Section 2.9, this behaviour is controlled via the effective ratio of specific heats in the equation of state.

As demonstrated above, the simulation domain in this model does not afford simplifications to the radiative energy transfer problem, firstly due to a lack of a priori knowledge of the conditions present, and secondly due to no desirable coupling between optical depth and timescale regimes. Other methods that deal explicitly with radiative energy transfer exist, but these cannot be satisfactorily applied here, as they usually cover only one optical regime, and normally require boundary conditions in order that energy can be lost at surfaces.

The new radiative energy transfer method presented in this section, was developed in order that the conditions present in the Capture Theory could be modelled as accurately as possible. Since all optical depth and timescale regimes are present in almost any combination, this led to a general SPH radiative energy transfer scheme being developed. Before the scheme is presented, radiative principles are discussed, and a first attempt at a general radiative scheme, including reasons for failure and what was learned from this attempt, is examined.

Radiative Principles and Assumptions

Both of the methods considered in this work essentially split the radiative transfer mechanism into two components: radiating and receiving. The simulated region is considered in terms of sub-regions of order of size $h$. To simplify the calculation, the gas is assumed to be in local thermal equilibrium (LTE). This implies that the temperatures of the components that constitute the gas are coupled together, i.e. hydrogen, helium, the metals, the ices, and the dust grains are all at the same local temperature. Certainly in the simulation domain under consideration here, the dust and gas temperatures are strongly coupled together, see Whitworth et al. (1998) for a fairly complete analysis of this coupling for a wide range of simulation domains. Given this assumption, it is possible to estimate the local Rosseland mean opacity, $\kappa$, from the temperature and density of the gas, as demonstrated later in Section 3.2.3. With the LTE approximation, the amount of energy emitted, and the amount of energy received, can both be related to the optical depth in the region under consideration:

A completely opaque region radiates energy at a rate according to a black body law:

\begin{displaymath} \frac{dE}{dt}=\sigma AT^{4} \end{displaymath} (3.1)

where $\sigma$ is the Stefan-Boltzmann constant, $A$ is the surface area of the radiating region, and $T$ is the surface temperature of that region. When the region is completely opaque, it also absorbs all radiation incident on it.

The optical depth of a region is defined as:

\begin{displaymath} \tau=\kappa\rho d \end{displaymath} (3.2)

where $\kappa$ and $\rho$ are the average values of the Rosseland mean opacity and the density in the region under consideration, and d is the depth of the region. The optical depth is used to attenuate the radiative energy emitted and received when the gas is not completely opaque. Consider the cubical region of space in Figure 3.10.

Figure 3.10: Radiation incident perpendicular to one face of a cuboid.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/cuberada1.eps,angle=0,width=0.66\linewidth}} \end{figure}

Here, radiation denoted by arrowed lines is incident perpendicular to a face. The proportion of the radiation that is absorbed in the region is given simply by $1-exp(-\tau)$, where $d$ is the depth of the region perpendicular to the absorbing face. Given LTE, the rate of radiative energy emitted perpendicular to a face is attenuated by this same amount:
\begin{displaymath} \frac{dE}{dt}=\sigma AT^{4}\Big\{1-exp(-\tau)\Big\} \end{displaymath} (3.3)

The above example is rather basic, but only a few modifications are required for it to be applicable to the tree scheme used in Section 3.2.5. In the first scheme (Section 3.2.4) more complicated geometry is employed, so that additional considerations need to be taken into account.


Opacity

The opacity of an astrophysical gas depends on the state and concentration of the individual components which compose the gas. Since some astrophysical applications are particularly sensitive to opacity3.2, a great deal of effort has been spent in recent years in investigating this dependence. This has led to accurate results, covering the whole of the Capture Theory's simulation domain, becoming available. Therefore in theory, every plausible combination of gas components and conditions can be accounted for. In practice, this is an impossible task, given the myriad of combinations that are possible. Thus some assumptions concerning the state and composition of the material throughout the simulation domain must be made.

A practical method is to use a set of opacities for one particular mix of materials, and apply a multiplication parameter to these opacities to simulate varying the composition. Incorporating a multiplication parameter is also useful for investigating how sensitive the simulation is to opacity, and consequently to radiative energy transfer in general. The mix used here is closest to what the solar-nebula theorists believe would have been present at various stages in their model. This is estimated from what is observed in star forming regions, the planets, and the surface of the Sun today. An in-depth description of the radiative processes present in this mix, which subsequently lead to an effective opacity, is not relevant here, instead see Alexander and Ferguson (1994); Gaustad (1963); Henning and Stognienko (1996); Alexander et al. (1983) for detailed discussions of the physics involved. However a brief qualitative description is in order:

At low temperatures the opacity is dominated by solid grains. Above $\approx 150K$ the volatile ices evaporate, but the residual refractory compounds in the grains still provide the major absorption. When the grains have completely evaporated, $\approx 1500K$, hydrogen molecules and $H^{-}$ are the major contributors to the opacity.

Two sources of opacity data are used in this thesis. The first is for temperatures up to $630K$, where only temperature dependence is necessary (Henning and Stognienko, 1996). The results from this paper were conveniently provided in a parameterised form:

\begin{displaymath} \kappa=aT^{b} \end{displaymath} (3.4)

where the values of $a$ and $b$ as a function of $T$ are given in Table 3.1. These values are plotted in Figure 3.11.

Table 3.1: Fit parameters for low temperature Rosseland mean opacities.
T(K) a b
10-135 1x$10^{-5}$ 2.1
135-155 1x$10^{0}$ -0.2
155-375 3x$10^{-3}$ 0.9
375-400 8x$10^{+4}$ -2.0
400-575 4x$10^{-2}$ 0.5
575-630 9x$10^{+18}$ -6.9



Table 3.2: Master grid of Rosseland mean opacities.
(-7,-4)(0,0)log T                             log R                                
   -4.0  -3.5  -3.0  -2.5  -2.0  -1.5  -1.0  -0.5  0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  
4.10  -1.493  -1.479  -1.439  -1.342  -1.147  -0.840  -0.430  0.023  0.476  0.887  1.217  1.453  1.618  1.742  1.851  1.966  2.096  2.228  2.374  2.527  
4.05  -1.509  -1.497  -1.447  -1.325  -1.097  -0.767  -0.350  0.078  0.467  0.773  0.987  1.135  1.247  1.344  1.447  1.570  1.717  1.885  2.056  2.237  
4.00  -1.529  -1.498  -1.429  -1.285  -1.054  -0.729  -0.354  -0.009  0.257  0.441  0.571  0.674  0.768  0.871  0.991  1.137  1.311  1.517  1.721  1.927  
3.95  -1.521  -1.482  -1.405  -1.267  -1.062  -0.802  -0.553  -0.357  -0.216  -0.107  -0.010  0.092  0.205  0.337  0.490  0.671  0.880  1.115  1.351  1.589  
3.90  -1.510  -1.475  -1.414  -1.326  -1.223  -1.126  -1.046  -0.971  -0.892  -0.799  -0.687  -0.554  -0.402  -0.230  -0.037  0.182  0.425  0.676  0.948  1.222  
3.85  -1.524  -1.527  -1.546  -1.593  -1.657  -1.713  -1.742  -1.725  -1.655  -1.539  -1.387  -1.210  -1.016  -0.804  -0.575  -0.326  -0.048  0.231  0.544  0.857  
3.80  -1.669  -1.790  -1.943  -2.111  -2.275  -2.412  -2.496  -2.500  -2.418  -2.268  -2.075  -1.862  -1.631  -1.386  -1.120  -0.826  -0.500  -0.154  0.207  0.541  
3.75  -2.112  -2.321  -2.537  -2.751  -2.952  -3.122  -3.229  -3.241  -3.147  -2.978  -2.757  -2.509  -2.235  -1.933  -1.593  -1.213  -0.819  -0.437  -0.087  0.225  
3.70  -2.765  -2.997  -3.224  -3.443  -3.644  -3.812  -3.921  -3.932  -3.836  -3.647  -3.387  -3.067  -2.690  -2.273  -1.851  -1.450  -1.078  -0.751  -0.441  -0.160  
3.65  -3.522  -3.742  -3.951  -4.147  -4.323  -4.465  -4.546  -4.530  -4.387  -4.113  -3.735  -3.307  -2.875  -2.469  -2.098  -1.759  -1.446  -1.176  -0.892  -0.628  
3.60  -4.297  -4.478  -4.642  -4.783  -4.889  -4.939  -4.912  -4.783  -4.545  -4.220  -3.853  -3.485  -3.136  -2.814  -2.509  -2.214  -1.926  -1.699  -1.440  -1.218  
3.55  -4.940  -5.033  -5.086  -5.100  -5.080  -5.030  -4.945  -4.817  -4.640  -4.420  -4.167  -3.896  -3.617  -3.330  -3.027  -2.718  -2.420  -2.333  -2.127  -1.958  
3.50  -5.182  -5.169  -5.147  -5.121  -5.102  -5.094  -5.086  -5.059  -4.996  -4.883  -4.695  -4.413  -4.019  -3.521  -2.979  -2.498  -2.169  -2.169  -2.169  -2.169  
3.45  -5.203  -5.207  -5.231  -5.274  -5.333  -5.404  -5.469  -5.503  -5.465  -5.279  -4.870  -4.270  -3.608  -3.014  -2.596  -2.363  -2.245  -2.245  -2.245  -2.245  
3.40  -5.387  -5.460  -5.550  -5.650  -5.745  -5.812  -5.818  -5.688  -5.319  -4.700  -3.990  -3.351  -2.885  -2.625  -2.507  -2.454  -2.428  -2.428  -2.428  -2.428  
3.35  -5.796  -5.896  -5.971  -6.004  -5.974  -5.832  -5.479  -4.881  -4.174  -3.529  -3.055  -2.802  -2.698  -2.658  -2.641  -2.631  -2.623  -2.623  -2.623  -2.623  
3.30  -6.133  -6.117  -6.060  -5.894  -5.515  -4.903  -4.194  -3.559  -3.134  -2.939  -2.869  -2.844  -2.834  -2.828  -2.825  -2.822  -2.820  -2.820  -2.820  -2.820  
3.25  -6.131  -5.888  -5.409  -4.724  -4.012  -3.463  -3.186  -3.088  -3.057  -3.045  -3.040  -3.037  -3.036  -3.034  -3.033  -3.032  -3.029  -3.029  -3.029  -3.029  
3.20  -5.011  -4.281  -3.702  -3.420  -3.329  -3.301  -3.291  -3.288  -3.285  -3.284  -3.284  -3.285  -3.285  -3.284  -3.284  -3.280  -1.881  -1.881  -1.881  -1.881  
3.15  -3.496  -3.448  -3.434  -3.429  -3.428  -3.427  -3.427  -3.427  -3.428  -3.428  -3.429  -2.133  -1.936  -1.756  -1.653  -0.974  -0.587  -0.587  -0.587  -0.587  
3.10  -3.423  -3.423  -3.423  -3.424  -3.424  -3.425  -2.346  -2.092  -1.862  -1.269  -0.940  -0.857  -0.646  -0.495  -0.407  -0.374  -0.343  -0.343  -0.343  -0.343  
3.05  -3.427  -3.427  -2.203  -1.957  -1.479  -1.045  -0.844  -0.743  -0.649  -0.551  -0.468  -0.431  -0.406  -0.401  -0.400  -0.399  -0.399  -0.399  -0.399  -0.399  
3.00  -1.263  -1.091  -0.889  -0.755  -0.671  -0.586  -0.512  -0.471  -0.454  -0.453  -0.452  -0.452  -0.452  -0.452  -0.452  -0.452  -0.452  -0.452  -0.452  -0.452  
2.95  -0.746  -0.680  -0.602  -0.534  -0.501  -0.499  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  -0.498  
2.90  -0.577  -0.542  -0.540  -0.539  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  -0.538  
2.85  -0.571  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  -0.570  
2.80  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  -0.590  


Figure 3.11: Low temperature opacities.
Figure 3.12: Opacity as a function of temperature at constant density. (a) $10^{-10}$, (b) $10^{-8}$, (c) $10^{-6}$, (d) $10^{-4}kg\,m^{-3}$.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/lowt... ...ve/tmp/hold/thesis/graphs/alltemp.eps,angle=0,width=1.0\linewidth}} \end{figure}

The second source covers $630K$ to $12500K$ (Alexander and Ferguson, 1994). Since the gas component of the medium is the major contribution to opacity at these temperatures, the opacity depends on density as well as temperature. It is usually the convention with temperature-density dependent opacities, to introduce a new variable, $R$, in order to facilitate rectangular tabulation of results:

\begin{displaymath} R=\frac{\rho}{T_{6}^{3}} \end{displaymath} (3.5)

where $T_{6}$ is the temperature in millions of degrees, and $\rho$ is normally expressed in c.g.s. units. For consistency in this thesis, $\rho$ is expressed in S.I. units.

These values of opacity (log $\kappa$), expressed as a function of log $T$ and log $R$, are given in Table 3.2. The opacities covered by this table range from log $R=-4.0$ to log $R=5.5$. These values of $R$ cover the density range $10^{-13}$ to $3$x $10^{-4}kg\,m^{-3}$ at $1000K$, and the density range $10^{-10}$ to $3$x $10^{-1}kg\,m^{-3}$ at $10000K$. The last three columns and the bottom row are extensions to the data of Alexander and Ferguson:

The bottom row, log $T=2.8$, was added so as to ensure a smooth transition between the two opacity data sets, all values being the same as the low temperature opacity data set at that temperature. The last three columns on the right (log $R=4.5$ to log $R=5.5$), were added so as to allow the results to be applicable at higher densities. For low temperatures in this region (below log $T=3.5$), the values are all equal to the values at log $R=4.0$ in the same row. Analysis of this table and other sources, show that this is a reasonably valid assumption. The values of the opacity above log $T=3.5$ (and log $R=4.0$) are estimated from other sources of opacity data available for mixes similar to this, and are seen to be consistent with the pattern of increase demonstrated by preceding opacities in the same row.

Simple linear interpolation is applied for intermediate values. If the material has properties which result in opacities outside of the table being required, the boundary values are used, and a warning flag set. This set of opacities, along with the first set, are plotted in Figure 3.12 for four sample densities.


A First Idea

A brief overview of a first attempt at including a general radiative energy transfer scheme into SPH is now given. The method is essentially a particle-particle transfer scheme, incorporating a uniform interparticle grid for absorption purposes.

Radiators

When considering the amount of energy radiated by the particles, each SPH particle is replaced with a sphere of uniform density centred on it. By assuming that this density corresponds to the value determined during the SPH force calculation routine, a radius for this sphere can be found:
\begin{displaymath} r_{i}=\Bigg(\frac{3m_{i}}{4\pi\rho_{i}}\Bigg)^\frac{1}{3} \end{displaymath} (3.6)

The sphere is assumed to radiate like a black body, attenuated by an amount similar to what was given in Equation 3.3, appropriate for a sphere:
\begin{displaymath} \frac{dE_{i}}{dt}=4 \pi \sigma r_{i}^{2}T_{i}^{4}\Big\{1-exp(-2\kappa_{i}\rho_{i}r_{i})\Big\} \end{displaymath} (3.7)

Receivers

The corresponding receiving element associated with each particle is taken as a uniform density truncated segment of a sphere centred on the radiating particle.

Figure 3.13: Radiators and receivers.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/firstrada3.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 3.13 shows this arrangement for two dimensions. In three dimensions, the choice of using a square or circular cross section is available. The exact definition of the segment used in practice is unimportant here, suffice to say many geometrical relationships are available.

Given the arrangement of Figure 3.13, the proportion of radiation absorbed from particle $i$, incident on the front face of the receiving element, is given by $1-exp(-\kappa_{j}\rho_{j} d)$. This rather simple result is the reason for representing the mass and energy of a particle as a truncated segment whilst it is receiving radiation, rather than as a cube or sphere, where calculating the amount absorbed from a point source is rather involved.

Interparticle Absorption Grid

A grid is required in order that absorption between particles can be accounted for, i.e. placing a particle between $i$ and $j$ in Figure 3.13, reduces the amount of radiative energy incident on $j$ due to $i$. A simple regular grid was employed for this purpose. In order that the absorption properties of a grid cell can be estimated, a temperature and density sampled from surrounding particles is calculated at the centre of the cell, and taken to apply throughout the cell. This can be achieved in various ways, but the modified SPH interpolation procedure given later in Section 3.2.7 was used.

When considering the radiative energy absorbed in $j$ due to $i$, several interparticle points are considered along the line joining the particles, between the surfaces of the radiating and receiving elements, the absorbing properties of which are estimated from the grid3.3. The amount of energy falling on the front face of the receiving element, from the radiating element, is attenuated by a factor determined using these points. The energy lost due to the interparticle absorption is not stored, as it is assumed that when the particles that contributed to the grid properties interact with $i$, this energy will be absorbed by them then.

Discussion

The merits of the method can be discussed in terms of physical accuracy and computational expense.

On the computational side, the method is essentially an $N^{2}$ code, as every particle interacts with every other, and is thus relatively easy to implement compared to tree based methods. However, this leads to large computational run-times when large numbers of particles are employed, and would not be computationally viable for tree based implementations of the force calculation. Moreover, a large computational overhead is present in the $N^{2}$ proportionality constant, due to the interparticle absorption calculation. Thus even an $N^{2}$ based force calculation routine would suffer considerably longer computational run-times, with the inclusion of this radiative energy transfer scheme. One way of circumventing these concerns is to apply the radiative energy transfer routine intermittently, e.g. every $10$ time-steps for an $N^{2}$ force calculation routine, and every $100$ time-steps say, for a tree based force calculation routine. Obviously this is not the most desirable of solutions.

Continuing with the computational side of the discussion, the use of a regular grid introduces resolution considerations. Specifically the number of cells in the grid must be such that the smallest detail in the simulation be resolved. In three dimensional clustered simulations, this can be an unachievable task, requiring many billions of grid cells, most of which will resolve space to an unnecessarily high degree of accuracy. Obviously the solution to this problem is to use a tree type grid structure, as demonstrated in earlier sections, and applied later in Section 3.2.5.

Physically the method performs well in the almost transparent and semi-opaque opacity regimes, producing errors at the $1\%$ level per time-step in test scenarios. However the method is beset with problems in the opaque regime. The major symptom is that during radiation tests, the amount radiated by the particles can be severely less than expected. Likewise, during absorption tests, the amount absorbed by the particles is less than expected. The reason for these problems stems from the relationship between the radiating, receiving, and grid geometries:

Figure 3.14: Radiation and absorption test scenario.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/firstproba2.eps,angle=0,width=1.0\linewidth}} \end{figure}

A simple demonstration of this can be given by considering the example scenario in Figure 3.14. Here the group of small circles on the right represents the radiating elements of the particles in a sphere of uniform density, represented approximately by the inner large circle. The outer large circle is the approximate extent to which the particles contribute to the grid properties, i.e. $2h$ from the edge particles. Beyond this, grid cells have zero density and temperature, and thus are completely transparent. The cross on the left of the diagram is a point source radiator, used to test the absorption properties of the particles.

Assuming that the particles are all at the same temperature initially, a theoretical rate of absorption in the sphere of particles from the point source can be determined. The rate of energy lost by the particles can also be found analytically. These details are not important here, but can be found in Section 4.3, where this scenario is used to test the tree radiation transport code. Only an instant in time is being considered in these tests, so force and velocity calculations are unnecessary. As already mentioned, the radiative energy transfer algorithm under investigation, performs well in the low opacity regime, where optical depths of the order of the size of the modelled region or greater are present. However, as the optical depth across the spheres increases significantly beyond unity, the simulation diverges severely from the analytical solution.

Considering the absorption test first: The reason for very little absorption in the sphere from the point source radiator, is due to the grid absorbing most of the radiation erroneously, as demonstrated in Figure 3.15(a), i.e. it is falsely assumed that the energy absorbed by the shaded region of the grid, will be absorbed by particles at a later stage.

Figure 3.15: Demonstrating erroneous grid absorption. (a) Absorption test for the point source shown in Figure 3.14, (b) Radiation test.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/firstradreca3.eps,angle=0,width=1.0\linewidth}} \end{figure}

The same is true for the radiation test, Figure 3.15(b), i.e. the surface spheres attempt to radiate to the outside world, but the grid intercepts the energy. This latter behaviour only manifests itself when sample points are placed at the edge of the simulation domain, and the energy received by them measured. This is because grid cells are only invoked during the absorption process. That does not mean that this behaviour is not important though, as it affects the energy transfered within the simulated sphere, and to other bodies that may be present in a simulation.

A possible way of alleviating these effects, could be to estimate the grid properties from the particles using the same uniform density spheres as determined for the radiating elements. However, even with this approach the erroneous absorption can still occur. This happens markedly when a particle's radiation sphere just overlaps onto the centre of a grid cell. This is because the properties at the centre of the cell are taken to apply throughout the whole cell, and thus a part of the cell could be beyond the desired radiation surface. A similar argument is applicable for the absorption elements.

For this effect to be completely suppressed, the geometry of the radiating and receiving elements must be the same as the grid geometry, i.e. square. It is also necessary to coincide the centre and size of the emitting and receiving elements, with the centre and size of the grid cells, or some multiple of them.

Conclusion

The conclusions to be drawn from this first attempt at including a generalised radiative energy transfer algorithm into SPH, can be summarised as:
  1. An $N^{2}$ radiation algorithm should be avoided for tree based force calculations.
  2. The use of a regular grid should be avoided for resolution reasons.
  3. In the large opacity regime, the radiating and receiving elements should be the same shape and size as the grid elements, and centred on them.


Tree Radiation Method

This tree implementation of radiative energy transfer does not violate any of the above conclusions. In fact, as already mentioned, it was developed as a result of the preceding investigation. By being tree based, both point 1 and point 2 are satisfied. Point 3 is automatically satisfied, as the grid is used for radiating, receiving and attenuating energy. Essentially the method uses a Monte-Carlo approach for transferring energy throughout the tree structure: The particle properties are spread onto the tree; energy is transfered within the tree via a random walk Monte-Carlo method; the resultant change in energy over the time-step is returned to the particles. Below is an expanded outline of this procedure:

  1. A tree is constructed in a similar fashion to that given in the force calculation section, but with a few minor alterations.
  2. The properties of the particles are `spread' onto the smallest squares in the tree (leafs), using the normal SPH technique adapted for this purpose.
  3. The optical properties of the leafs are calculated, and the optical properties of larger squares calculated from the sub-squares and leafs they contain, via a tree traversal.
  4. The rate of radiative energy loss from every leaf in the tree is calculated.
  5. This rate of energy loss (power) is divided equally to form approximately $100$ packets of power for each leaf, which are then sent on constrained random walks in all directions through the tree. The power packets either leave the tree or are absorbed in other leafs.
  6. The resulting change in leaf temperature, due to this transfer of energy, is calculated over the dynamical time-step, and this change is returned to the particles in an appropriately weighted manner.
Each step in the tree is now examined in detail.


Defining the Tree

Tree Construction

The tree is essentially a Barnes-Hut tree, but with a few modifications. Firstly, since the properties of the particles will be spread onto the tree at a later stage, the definition of the tree should be such that space is sufficiently resolved in all directions within $2h$ of a particle. The usual procedure for tree construction, see Section 3.1.3, can result in situations as shown in Figure 3.16(a), especially at surfaces.

Figure 3.16: Radiative tree construction. (a) Using N particles, (b) Using 5N (the 2D equivalent of 7N) particles.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/tree7na1.eps,angle=0,width=1.0\linewidth}} \end{figure}

This example shows a small region of a simulation, where a large fraction of the properties of the particles on the left, will not be sampled sufficiently by the tree. The procedure applied here to prevent this from happening, is to create a list of particles $7N$ large, which is used only for the purpose of constructing the radiation tree. The list consists of the positions and smoothing lengths of the normal set of particles, plus $6N$ particles offset by $2h_{i}$ in the $-x_{i}$, $+x_{i}$, $-y_{i}$, $+y_{i}$, $-z_{i}$, $+z_{i}$ directions for each particle (in $3D$). A smoothing length equal to the central particle is associated with each of these `virtual' particles, which is used for pruning the tree in the following section. Constructing the tree with this set of particles for the example in Figure 3.16(a), results in Figure 3.16(b). Leafs which do not contain a particle are not deleted, as they are used in later stages.

Tree Pruning

The Barnes-Hut tree construction algorithm can result in an unnecessarily high resolution around some particles. This occurs when particles are very close to each other, and is amplified by using $7N$ particles in the construction process. Clearly the resolution around particles should be of the order of a smoothing length, much higher resolution being an unnecessary waste. This is accomplished by a process coined ``tree pruning''. A square is pruned, meaning that all of its sub-squares are deleted, if all of its sub-squares are deemed to resolve space too finely. This condition is necessary in order that the leaf level of the tree covers space completely. A square is defined as resolving space too finely, if the average smoothing length of the particles contained within it, is greater than some factor of the size of the square:

\begin{displaymath} \bar{h}_{d}>P_{F}d \end{displaymath} (3.8)

where $\bar{h}_{d}$ is the average smoothing lengths contained in a square of sides $d$, and $P_{F}$ is a parameter which controls the degree of pruning and hence the resolution of the tree, called the pruning factor. Typical values of the pruning factor range from $1.5$ to $3$; higher values correspond to higher spatial resolutions. Unless otherwise stated $P_{F}=2$ is used here. This corresponds to approximately $10$ times more leafs than particles, and approximately $200$ leafs sampling each particle, in typical three dimensional simulations.

The tree pruning procedure is implemented in the opposite sense to what has been explained above, i.e. all squares are marked as `prunable', and a tree traversal determines which ones should not be pruned:

  1. Initially every square in the tree is marked as prunable.
  2. The average smoothing length of the particles contained in every square of the tree is calculated via a progressive accumulation during one tree traversal, similar to the centre of mass calculation in Section 3.1.5.
  3. Whilst traversing the tree, if Equation 3.8 for a square is false, the square's parent is marked as being unprunable. All squares above an unprunable square, are also marked unprunable during the same tree traversal, regardless of Equation 3.8 for their sub-squares.
  4. The tree is then copied to a more compact tree, only copying the parts of the tree with parent squares that are marked unprunable. The procedure for this results in almost linear memory access during tree traversals, as the copying process is performed using a tree traversal. See Section 3.1.4 for a discussion of the advantages of this.
The pseudo-code for the tree pruning and subsequent compact tree creation is given in Appendix A.6 and Appendix A.7 respectively.

Figure 3.17: Pruned radiation trees. (a) For the example in Figure 3.16, (b) For the example in Section 2.1.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/treepruneda2.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 3.17(a) shows the tree structure which may result following pruning for the example in Figure 3.16. Figure 3.17(b) shows the tree structure which may be generated for the example in Section 2.1.


Interpolating the Particle Properties to the Tree

Essentially the same tree procedure as given in Section 3.1.7 for calculating the density, is used to estimate the properties of the leafs from surrounding particles. The major difference being that the density at the centre of a leaf is calculated using:

\begin{displaymath} \rho_{g}=\sum_{i}^{N}m_{i}W(r_{ig},h_{i}) \end{displaymath} (3.9)

and that the internal energy per unit mass is calculated using:
\begin{displaymath} u_{g}=\frac{\sum_{i}^{N}u_{i}m_{i}W(r_{ig},h_{i})}{\rho_{g}} \end{displaymath} (3.10)

where the subscript $g$ denotes the centre of a leaf. The internal energy per unit mass can be calculated during the same routine as the density, by dividing by the $\rho_{g}$ in Equation 3.10 once both summations are complete. Since grids are scarcely used in SPH, Equation 3.10 is an extension to the usual SPH equations. Note that the internal energy is defined per unit mass, thus it, and hence the temperature, is uniform across the sphere of an isolated SPH particle, unlike the density which is bell shaped. Thus the internal energy per unit mass calculation of Equation 3.10 is a linear weighting from surrounding particles, the fraction contributed to a leaf by a particle, being proportional to the fraction of density contributed to the leaf by the particle in Equation 3.9.

By assuming that the properties calculated at the centre of the leaf, apply throughout the whole leaf, a mass and consequently an internal energy can be derived:

\begin{displaymath} m_{g}=\rho_{g}V_{g} \end{displaymath} (3.11)


\begin{displaymath} U_{g}=m_{g}u_{g} \end{displaymath} (3.12)

where $V_{g}$ is the volume of the leaf, and $U_{g}$ and $m_{g}$ correspond to the total internal energy and the mass associated with leaf $g$ respectively.

When adopting this procedure for calculating the properties on the grid, the total mass and internal energy is typically conserved to within $1\%$ or so, depending on the pruning factor. If the change in internal energy were to be returned to the particles in an absolute manner, this discrepancy may cause problems related to energy conservation, however the change is returned in such a way that these problems are avoided. See Section 3.2.13 for the details of this procedure.


Tree Absorption Properties

This section is concerned with calculating the absorption properties at every level in the tree.

Given the average density and temperature of a leaf, the opacity can be calculated as demonstrated earlier in Section 3.2.3, and consequently the optical depth, as given by Equation 3.2. The optical depth is used to calculate the proportion of radiative energy incident on a region that is absorbed. Since the radiative energy transfer process is discretised, as described in 3.2.9, the proportion of energy absorbed is more appropriately referred to as the probability of absorption of a region. Furthermore, in practice it is more convenient to calculate the transmission probability, denoted by $Z$, rather than the absorption probability of the leafs in the tree. Obviously the absorption probability is simply $1-Z$.

As will be seen later in Section 3.2.9, radiative energy only travels along the Cartesian axes in this method. This means that radiative energy is always incident perpendicular to a face of a square. Thus for the leafs, the transmission probability of radiative energy incident on a face is given by:

\begin{displaymath} Z=exp(-\tau) \end{displaymath} (3.13)

where $\tau$ is the optical depth of the leaf. This value is valid for all faces, as the properties calculated at the centre of a leaf are assumed to apply throughout the leaf uniformly in all directions. When considering the transmission probability of squares above the leafs, more sophistication than using only the average properties of the leafs they contain is available. This is achieved by weighting the transmission probabilities of the sub-squares a square contains, to form a transmission probability along each Cartesian axis (3 values of $Z$ in 3D) for each square. The appropriate weighting for this procedure can be demonstrated by considering the example in Figure 3.18.

Figure 3.18: Calculating the tree absorption properties.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/opacitycona1.eps,angle=0,width=0.56\linewidth}} \end{figure}

Squares $1$ to $4$ and squares $6$ to $8$ are leafs in a two dimensional tree structure. Transmission probabilities are known for these, the task is to find the transmission probabilities in the $x$ and $y$ directions of squares $5$ and $9$. Essentially this is trivial probability analysis, where multiplication for squares in series, and addition for squares in parallel is appropriate. The transmission probabilities for this example are:
\begin{displaymath} {\setlength\arraycolsep{2pt} \begin{array}{cl} Z_{5x}=&\frac... ...})\\ Z_{9y}=&\frac{1}{2}(Z_{5y}Z_{7}+Z_{6}Z_{8}) \end{array}} \end{displaymath} (3.14)

A tree traversal is used to accumulate the transmission probabilities of squares at all levels in the radiation tree.


Radiative Energy Transfer Process

The purpose of this and the following three sections is to estimate the change in energy at the leaf level of the tree, due to radiative processes during the dynamical time-step. Instead of attempting to consider every leaf with every other leaf, or perhaps considering a leaf with a group of leafs sufficiently far away (as in the gravitational tree force calculation), the radiative energy emanating from a leaf is absorbed in a sample distribution of surrounding leafs.

The first stage in the method involves calculating a rate of radiative energy loss from the surface of the leafs:

\begin{displaymath} \frac{dU_{g}}{dt}=6\sigma d^{2}_{g}T_{g}^{4}(1-Z_{g}) \end{displaymath} (3.15)

This rate of radiative energy loss (luminosity, or power) is divided equally into approximately $100$ units of power for each leaf. These units of power are termed ``power packets''. The actual number of these depends on the optical depth of the emitting leaf. Fewer power packets are used, both for very low optical depths, where power packets mostly leave the system, and for very high optical depths, where they are absorbed in neighbouring leafs. The power packets are transmitted perpendicular to the surface of the leafs, on constrained random walks through the tree structure. Using a constrained random walk reduces the number of steps required for power packets to leave the system in semi-opaque and transparent opacity regimes. It also produces a smoother sampling distribution when small numbers of power packets are employed. The final reason for using a constrained random walk is to soften the problem discussed in Section 3.2.10, where it is seen that using the random walk sampling procedure, artificially enhances the absorption probability of surrounding leafs.

When a power packet passes through a leaf, it is either absorbed in the leaf, or `scattered' from it. In a true tree random walk process, the scattering could be in any of the six axial directions, however, here the power packet cannot be scattered in the direction it came from, or back along the axis it was originally transmitted. For a uniform grid, i.e. all leafs at the same depth in the tree (and hence the same size), this results in power packet walks similar to those demonstrated in Figure 3.19.

Figure 3.19: Example constrained random walks for uniform leaf depth in the three opacity regimes. (a) Opaque, (b) Semi-opaque, (c) Almost transparent.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/randomwalka1.eps,angle=0,width=0.7\linewidth}} \end{figure}

In this example, sample power packets transmitted from leafs in each of the three opacity regimes are shown. The random walk method is of most interest in the semi-opaque regime, as in the other two, the power packets either leave the system, or are absorbed after their first step.

This uniform example is simple to implement: When a power packet moves into a leaf, a random number between $0$ and $1$ is generated. If the random number is less than $Z_{g}$ for that leaf, the power packet is scattered from the leaf, otherwise it is absorbed in it. However, the non-uniform general tree structure is more complicated, requiring additional considerations.

Figure 3.20 shows four scenarios which can occur when different levels of leafs are present. The scenarios are ordered in increasing complexity.

Figure 3.20: Power packet walk scenarios. (a) Leaf same size in the same square, (b) Leaf same size in a different square, (c) Larger leaf, (d) Subdivided square.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/randomlevela1.eps,angle=0,width=0.55\linewidth}} \end{figure}

Leaf same size in the same square, Figure 3.20(a)

In this scenario, a power packet is simply absorbed or scattered from a neighbouring leaf, in the manner just explained.

Leaf same size different square, Figure 3.20(b)

This scenario is similar to the previous one, the only difference being that a tree boundary is crossed. In order to find the appropriate leaf, a tree neighbour search is performed, see Section 3.2.10 for the details. This procedure involves specifying a leaf, and a direction of travel. The tree is searched upwards until it is possible to move in the specified direction to an adjacent square of the same size, without crossing a tree boundary (the found squares will share the same parent square). The tree is then searched downwards, descending to squares which are adjacent to the original leaf, until the required leaf on the other side of the boundary is found. The power packet is then either absorbed or scattered in the usual manner. If the neighbouring leaf is of a different size to the current leaf, then the following two scenarios are appropriate.

Larger leaf, Figure 3.20(c)

This scenario is exactly the same as the previous one, except that whilst descending the tree, if the leaf level is encountered prematurely, the power packet is either scattered or absorbed from the leaf currently being examined in the descent.

Subdivided square, Figure 3.20(d)

If necessary, the same neighbour searching procedure as used previously is performed. So as not to introduce additional complexity, or reduce the resolution unduly in the tree walking procedure, and hence the simulation in general, an additional criteria is introduced: A power packet can only be absorbed or scattered at the leaf level of the tree. This has already been implied to some extent in the previous scenarios, but here the definition is required explicitly, as the choice of using a subdivided square for these purposes is clearly possible.

The approach adopted to ensure that this criteria is met in a computational efficient and physically realistic manner, involves generating a random number between $0$ and $1$ as before, and comparing it against the $Z_{g}$ of the subdivided square, for the current axis the power packet is travelling along. If the random number is less than this transmission probability, the power packet is descended down the tree to completely random sub-squares, until a leaf is reached, at which point the power packet is scattered in the usual constrained manner from it.

If the random number is greater than the transmission probability of the square in question, the power packet is absorbed in one of the leafs contained in it. The leaf it is absorbed in is determined by consideration of the transmission probabilities of the square and leafs it contains. Again, trivial probability analysis is appropriate, but this time the direction along an axis is important. Consider the example of Figure 3.18 again, where it is assumed that a power packet from leaf $7$ has just been absorbed somewhere in square $5$, Figure 3.21.

Figure 3.21: Sub-square absorption example.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/opacitycontrana1.eps,angle=0,width=0.6\linewidth}} \end{figure}

The leafs of square $5$ have the following probabilities of absorbing the power packet, the sum of which totals $1$:
\begin{displaymath} {\setlength\arraycolsep{2pt} \begin{array}{cl} \vspace{6pt} ... ...laystyle{\frac{1}{2}}\frac{(1-Z_{4})}{(1-Z_{5y})} \end{array}} \end{displaymath} (3.16)

If the square which absorbs the power packet is subdivided, similar absorption probabilities are calculated for its sub-squares, and the tree is descended until the power packet is absorbed in a leaf. In three dimensions the factor of half in Equations 3.16 is replaced with a quarter. The absorption probability of the parent square is required in the divisor, so that the probability that the power packet is absorbed in one of its sub-squares is $1$.

Figure 3.22: Example constrained random walks for the multiple leaf depth example of Section 2.1.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/randomexamplea1.eps,angle=0,width=1.0\linewidth}} \end{figure}

In practice, all of the above scenarios are combined together in the same routine, as it is not known a priori which of the above scenarios will be encountered during the random walk. Figure 3.22 shows some example constrained random walks for the multiple leaf depth example of Section 2.1. See Appendix A.8 for the pseudo-code of a tree walk for one power packet.

Additional Tree Walking Considerations


Neighbour Search

This procedure is required when a power packet crosses a boundary, i.e. it moves to a square which does not share the same parent as it. The square is found by searching upwards in the tree, until it is possible to move in the required direction without crossing a tree boundary. Whilst ascending the tree, a list containing the relative positions of the ascended squares with respect to other squares that share the same parent is constructed. The list is read in reverse order during the descent, i.e. it is a stack. To determine which squares should be descended into, the current list value is compared with the direction of travel.

As an example of this, consider Figure 3.23.

Figure 3.23: Tree neighbour search. (a) Labelling notation, (b) Example search.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/treeneigha1.eps,angle=0,width=1.0\linewidth}} \end{figure}

Figure 3.23(a) shows the labelling notation used in this example. The stack generated for the power packet example of Figure 3.23(b) is `$331$'. Since the direction of travel is to the right, a value of $-1$ is added to each value of the stack to form `$220$'. With this labelling notation, $+1$, $-2$, and $+2$ are added to each value of the stack, for directions of travel leftwards, upwards, and downwards respectively. In this example, the tree is ascended until the top left square labelled $2$ is reached, at which point it is determined that a movement to the right is possible without crossing a boundary. The tree is then descended using the stack, i.e. first to sub-square $0$ of square $3$, then to square $2$, and finally to the required leaf, labelled $2$.

This procedure results in a full tree climb when a power packet crosses the central axes of the largest square in the tree.


Absorption Softening

In using random walks to sample the radiative energy emanating from a leaf, the probability of absorption of surrounding leafs is artificially enhanced. This can be demonstrated by considering Figure 3.24.

Figure 3.24: Example random walks demonstrating artificially enhanced absorption.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/treeradsofta1.eps,angle=0,width=0.57\linewidth}} \end{figure}

Here two sample random walks are shown. The long arrowed lines indicate the depth of material the radiative energy would have travelled through for a more physically realistic direct interaction between the leafs. As can be seen, the sampling process has artificially increased the absorption depth by just over a factor of two in these examples.

The theoretical basis of the procedure used here for softening this effect, involves keeping track of the total depth of material which has already been considered for absorption during the random walk of the power packet. The absorption depth of the current square is reduced to the difference in distance between the total absorption depth already considered, and the distance between the centre of the current square and the centre of the radiating leaf. The total absorption depth is then incremented by this amount.

The above procedure works well for uniform grids, almost completely counteracting the surplus absorption, however it cannot be applied directly to a tree structure. Instead an approximation to this method, which uses the distance from the radiating leaf to both the previous leaf and the current square, is used to reduce the absorption depth.

Figure 3.25: Example absorption softening.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/treeradexama1.eps,angle=0,width=0.62\linewidth}} \end{figure}

Figure 3.25 shows the distances required for a typical multi-depth random walk. The transmission probabilities at each stage in the walk are reduced as follows:
\begin{displaymath} Z\leftarrow Z^{\epsilon} \end{displaymath} (3.17)

where $\epsilon$ is the difference in distances divided by the size of the current square, and therefore $\epsilon\le 1$. For the example in Figure 3.25 these are:
\begin{displaymath} \epsilon_{1}=\frac{r_{1}}{d_{1}}: \epsilon_{2}=\frac{r_{2}-r... ...r_{3}}{d_{4}}: \epsilon_{5}=\frac{r_{5}-r^{\prime}_{4}}{d_{5}} \end{displaymath} (3.18)

Dividing by the size of the cell, $d$, removes the absorption depth which was used in the original transmission probability calculation. In this method, if $\epsilon$ is less than $0$, as is $\epsilon_{5}$ in this example, a value of $\epsilon=0$ is taken, otherwise a transmission probability of greater than $1$ is obtained. In theory, the negative distance this corresponds to should be stored, so that it can be subtracted from subsequent positive distances. Accounting for these negative distances for random walks through a uniform grid, is equivalent to the theoretical method of storing the total absorption depth discussed earlier.


Time integration

During the random walk radiative energy transfer routine, a list summarising each power packet transfer is compiled. The list is referred to as the power packet transfer list. Each value in the list corresponds to a power packet, and contains a link to the radiating leaf and a link to the absorbing leaf, or zero if the power packet left the system. The list is used to calculate the change in energy of the leafs during the time-step.

Due to the strong dependence on temperature, the radiative energy transfer equations could potentially require a much shorter time-step for explicit integration than the dynamical time-step. Usually in simulations where some part of the physics requires a much shorter time-step than the rest of the system, an implicit time integration scheme is applied. One such scheme commonly used in radiation hydrodynamics is the semi-implicit Crank-Nicholson integration scheme, for example see Stone et al. (1992). In that work, Newton-Raphson iteration is used to solve both the hydrodynamical equations of motion, and the radiation transport equations simultaneously over a dynamical time-step.

A simple, but effective method which requires neither implicit integration or iteration is implemented here. More complicated methods were investigated, but all proved to have apparently intractable problems associated with them. For example, whilst incorporating an implicit iterative integration scheme, problems related to the rapid transfer of energy through the system were encountered, the source of which appeared to be inherent to iterative schemes in general.

In this simple method, the particle properties excluding radiative energy transfer, are integrated over the dynamical time-step as described in Section 2.11. The average of the particle properties at the beginning of the dynamical time-step and at the end of the time-step, are used to construct the radiation tree and calculate the properties on it. The dynamical time-step is divided into approximately $100$ smaller radiation time-steps, and a forward Euler integration scheme is used to calculate the change in energy of the leafs over these smaller time-steps. The change in energy is returned relatively to the particles from the leafs. The difference in energies at the average particle positions, before and after radiative energy transfer, is then applied absolutely to the particle energies at the end of the time-step.

When an approach like this is adopted, i.e. additional physical processes are solved outside of the dynamical equations, FSAL, see Section 2.11, cannot be implemented. This is because the extra change in particle properties at the end of the time-step, now means that the most recent derivative calculation of the dynamical properties at the end of the time-step, cannot be approximated to correspond to the first derivative calculation at the beginning of the next time-step. The main reason for not solving the radiative energy transfer equations simultaneous to the dynamical ones, stems from the random nature of the radiative energy transfer process, and the corruption in time-step prediction it would cause: The time-step prediction algorithm relies on the difference in particle properties calculated between the forward Euler stage and the improved Euler stage. If this difference is swamped with a random component from the radiative energy transfer equations, the algorithm will cease to perform efficiently, or even not at all, i.e. every time-step being failed, due to the tolerance criteria never being met.

Implementation

The procedure for calculating the change in energy of a leaf over a radiation time-step, using the power packet transfer list, is now given. At the beginning of the first radiation time-step, the rate of heating in each leaf is calculated by summing together the rate of heating associated with each power packet absorbed in the leaf, via one loop over the power packet transfer list. This is combined with the rate of energy radiated from the leafs given in Equation 3.15, to form the total rate of change in energy of the leafs at the beginning of the current time-step:

\begin{displaymath} \frac{dU^{T}}{dt}=\frac{dU^{A}}{dt}-\frac{dU^{R}}{dt} \end{displaymath} (3.19)

where the superscripts $T$, $A$, and $R$ are abbreviations for the $T$otal rate of change of energy of a leaf, the rate of change due to $A$bsorption, and the rate of change due to $R$adiation respectively. The new leaf energies at the beginning of the next radiation time-step are calculated according to:
\begin{displaymath} U_{1}=U_{0}+\Delta t_{R}\frac{dU^{T}_{0}}{dt} \end{displaymath} (3.20)

where $\Delta t_{R}$ is the radiation time-step, and $\frac{dU^{T}_{0}}{dt}$ is the total rate of change of energy at the beginning of the time-step.

In theory after this first time-step, new tree absorption properties should be calculated, and the power packet transfer routine performed again. However, although calculating the tree absorption properties takes inconsequential computational run-time, the power packet transfer routine is very computationally demanding. Instead, the approximation that the opacities are constant over the dynamical time-step is applied. With this approximation, both the recalculation of the tree absorption properties, and the power packet transfer routine can be avoided.

The modified rate of energy radiated by a leaf at the beginning of radiation time-step $n$ is given by:

\begin{displaymath} \frac{dU^{R}_{n}}{dt}=\frac{dU^{R}_{0}}{dt}\Bigg(\frac{U_{n}}{U_{0}}\Bigg)^{4} \end{displaymath} (3.21)

where although the zero subscripts correspond to the values at the beginning of the first radiation time-step, using values at any previous radiation time-step is equally valid. The power packet transfer list is used to calculate the new rate of absorption in the leafs, following which new total rates of change of energies are calculated. The leaf energies are then integrated to the beginning of the next radiation time-step using the generalised form of Equation 3.20:
\begin{displaymath} U_{n+1}=U_{n}+\Delta t_{R}\frac{dU^{T}_{n}}{dt} \end{displaymath} (3.22)


Rate of Transfer Considerations

In the semi-opaque and almost transparent regime, the speed of light is the limiting factor for the rate at which energy can travel through the system during a time-step. When the simulations in this work are in this regime, large time-steps of the order of years are present. In one year, light travels $63200 AU$, which is around 3 orders of magnitude greater than the dimensions of the largest simulation. Even if the dynamical time-step is divided into hundreds of radiation time-steps, this factor is of no concern here, however the rate of transfer in opaque regimes does require consideration.

It is apparent that in completely opaque regimes, the rate of energy flow through the system, using the methods presented in earlier sections, is governed by the size of the leafs and the number of radiation time-steps. This is because the energy absorbed by one face of a leaf during a radiation time-step, is assumed to have been absorbed uniformly within it, i.e. all faces of a leaf have the same effective surface temperature at all times. Thus radiation incident on one face of a leaf, effectively travels to the opposite face at a rate proportional to the size of the leaf. The following analysis was used to derive a term which could be used to reduce the flow of energy to a more realistic rate in opaque regimes.

In Figure 3.26,

Figure 3.26: Radiation incident perpendicular to one face of a leaf.
\begin{figure}\centering {\epsfig{figure=/home/steve/tmp/hold/thesis/graphs/rateatta1.eps,angle=0,width=0.7\linewidth}} \end{figure}

radiation denoted by arrowed lines is incident perpendicular to a face of a leaf of size $d$. If the leaf is opaque, then in theory most of the radiation should be absorbed close to its surface. However, short of making major modifications to the code, the energy absorbed is actually distributed evenly within the leaf. To approximate the desired surface absorption, the theoretical average penetration depth is used to reduce the magnitude of the radiative energy which is absorbed throughout the leaf. The average depth the radiation penetrates into the leaf can be calculated by considering how the intensity of the radiation varies as a function of penetration depth:
\begin{displaymath} I_{x}=I_{0}exp(-\kappa_{g}\rho_{g}x) \end{displaymath} (3.23)

where $I_{0}$ is the intensity of the radiation incident on the surface of the leaf, $I_{x}$ is the intensity at a penetration depth $x$ into the leaf, and $\kappa_{g}$ and $\rho_{g}$ are the average opacity and density of the leaf respectively. The average absorption depth is thus:
\begin{displaymath} \bar{x}=\frac{\int_{0}^{d}x\,exp(-\kappa_{g}\rho_{g}x) dx}{\int_{0}^{d}exp(-\kappa_{g}\rho_{g}x) dx} \end{displaymath} (3.24)

Performing the integration and substituting for $\tau_{g}=\kappa_{g}\rho_{g}d$, and $Z_{g}=exp(-\tau_{g})$ results in:
\begin{displaymath} \frac{\bar{x}}{d}=\frac{\Big\{(\tau_{g}+1)Z_{g}\Big\}-1}{\tau_{g}(Z_{g}-1)} \end{displaymath} (3.25)

Inspection of Equation 3.25 shows that it has the correct limiting behaviour, i.e. as the optical depth reduces to zero (less than $\approx 0.1$), $\bar{x}/d$ tends to $1/2$, and as the optical depth tends to infinity (greater than $\approx 5$), $\bar{x}/d$ tends to $1/\tau_{g}$. The former limiting behaviour is naturally correct, whilst the latter results in a transfer of energy, similar to what would have been obtained if the diffusion approximation for radiative energy transfer would have been invoked, as demonstrated later.

It is not until after the power packet transfer list is compiled, that the penetration depth calculation of Equation 3.25 is exploited: Whilst looping over the power packet transfer list, collecting the rate of change of energy in a leaf due to other leafs, the total power absorbed in the leaf is reduced by an amount that depends on that leaf's average penetration depth, $1-(2\bar{x}/d)$. In order to conserve energy, the fraction of power that is not absorbed by the leaf is returned to the leafs where the power packets responsible for the heating originated.

In opaque regimes, the rate of energy flow between two adjacent leafs (denoted by $1$ and $2$) of equal size, is thus approximately given by:

\begin{displaymath} \frac{dU_{12}}{dt}=\frac{2d^{2}\sigma}{\bar{\tau}}(T^{4}_{1}-T^{4}_{2}) \end{displaymath} (3.26)

This equation can be compared with that obtained from the diffusion approximation: Using Fick's first law, the diffusion approximation for photons can be expressed as:
\begin{displaymath} S=-\frac{\lambda c}{3}\nabla U_{V} \end{displaymath} (3.27)

where $S$ is the energy flux, $\lambda$ is the mean free path of a photon, $c$ is the speed of light, and $U_{V}$ is the total energy of photons per unit volume. For an opaque gas, the mean free path is related to the opacity and density by:
\begin{displaymath} \lambda=\frac{1}{\rho\kappa} \end{displaymath} (3.28)

and the photon energy density is given by:
\begin{displaymath} U_{V}=\frac{4\sigma T^{4}}{c} \end{displaymath} (3.29)

where $\sigma$ is the Stefan-Boltzmann constant. The rate of energy flow between adjacent leafs of equal size, as given by the diffusion approximation, is thus:
\begin{displaymath} \frac{dU_{12}}{dt}=\frac{4d^{2}\sigma}{3\bar{\tau}}(T^{4}_{1}-T^{4}_{2}) \end{displaymath} (3.30)

which is to within a factor of two of the result obtained in Equation 3.26.

In practice, when calculating the rate of energy transfer between leafs, the average of the attenuation factors derived from Equation 3.25 is used. This is necessary to prevent super-heating occurring around $1500K$, where the opacity drops several orders of magnitude, Figure 3.12. Super-heating occurs because of the severe imbalance of the radiation flow. For example, if the centre of an opaque body were to reach a sufficiently hot temperature, such that the optical depth across the centre fell by several orders of magnitude, then the material present would be subject to extreme heating from the surrounding opaque regions. This is because the surrounding opaque regions are able to transfer energy to the central region, but the central region is unable to transfer energy back effectively to the opaque regions, due to the large difference in the attenuation factors.


Interpolating the Tree Energy to the Particles

This section is concerned with returning the change in energy at the leafs, as the result of the transfer of radiative energy during a dynamical time-step, back to the SPH particles. There are several possible methods which could be used for achieving this. For example, one method could be to perform the inverse of the procedure set out in Section 3.2.7 for calculating the leaf properties, i.e. the smoothing length of a particle could be used to collect energy from the leafs. This method, although relatively straightforward, has several problems associated with it. One problem is related to the conservation of energy. In general, the total energy of all the leafs is not equal to the total energy of all the particles, due to the finite resolution of the tree. This implies that although the energy of the leafs may remain constant during a time-step, the energy determined at the particles may not be so.

In order to relieve these energy conservation concerns, the change in energies of the leafs are collected by the particles. The proportion of the change in energy of leaf $g$, which is collected by particle $i$, depends on how much of leaf $g$'s original energy was contributed to by that particle:

\begin{displaymath} U^{\prime}_{i}=U_{i}+\sum_{g}(U^{\prime}_{g}-U_{g})F_{ig} \end{displaymath} (3.31)

where $U^{\prime}_{i}$ is the new total energy of particle $i$ following the transfer of radiative energy, $U_{i}$ is the original total energy, $U^{\prime}_{g}$ is the new total energy of leaf $g$, $U_{g}$ is the original total energy, and $F_{ig}$ is the fraction of leaf $g$'s original energy that was due to particle $i$. There are several approaches for calculating $F_{ig}$, all of which give a similar, if not identical result:
\begin{displaymath} F_{ig}=\frac{U_{ig}}{U_{g}} \end{displaymath} (3.32)

Using Equation 3.12 to expand the contribution from particle $i$ to the total energy of leaf $g$:
\begin{displaymath} U_{ig}=m_{ig}u_{ig} \end{displaymath} (3.33)

Using Equation 3.11 to expand the contribution from particle $i$ to the mass of leaf $g$, in terms of the density contributed from $i$ and the volume of the leaf:
\begin{displaymath} m_{ig}=\rho_{ig}V_{g} \end{displaymath} (3.34)

Using Equation 3.10 to expand the contribution from particle $i$ to the specific internal energy of leaf $g$3.4.
\begin{displaymath} u_{ig}=\frac{u_{i}m_{i}W(r_{ig},h_{i})}{\rho_{ig}} \end{displaymath} (3.35)

Combining the above equations leads to:
\begin{displaymath} F_{ig}=\frac{U_{i}V_{g}W(r_{ig},h_{i})}{U_{g}} \end{displaymath} (3.36)

Inserting this into Equation 3.31, and rearranging into a more convenient form, results in the equation implemented here:
\begin{displaymath} u^{\prime}_{i}=u_{i}\bigg(1+\sum_{g}\Bigg\{\frac{U^{\prime}_{g}}{U_{g}}-1\Bigg\}V_{g}W(r_{ig},h_{i})\bigg) \end{displaymath} (3.37)

Since $u_{i}$ was calculated using the average values at the beginning and end of the dynamical time-step, the difference $(u^{\prime}_{i}-u_{i})$ is added to the specific internal energies of the particles at the end of the dynamical time-step.

Front Page
Stephen Oxley 2002-01-19