SPAMCART: a code for smoothed particle Monte Carlo radiative transfer
We present a code for generating synthetic SEDs and intensity maps from Smoothed Particle Hydrodynamics simulation snapshots. The code is based on the Lucy (1999) Monte Carlo Radiative Transfer method, i.e. it follows discrete luminosity packets as they propagate through a density field, and then uses their trajectories to compute the radiative equilibrium temperature of the ambient dust. The sources can be extended and/or embedded, and discrete and/or diffuse. The density is not mapped onto a grid, and therefore the calculation is performed at exactly the same resolution as the hydrodynamics. We present two example calculations using this method. First, we demonstrate that the code strictly adheres to Kirchhoff’s law of radiation. Second, we present synthetic intensity maps and spectra of an embedded protostellar multiple system. The algorithm uses data structures that are already constructed for other purposes in modern particle codes. It is therefore relatively simple to implement.
keywords:radiative transfer – hydrodynamics – methods: numerical – ISM: dust
Monte Carlo Radiative Transfer (MCRT) is a stochastic method for simulating radiative transfer through a medium. Individual MCRT calculations are accurate but slow. Fortunately, they are trivial to parallelise and therefore well suited to modern multi-threaded CPUs. MCRT is often used to post-process simulation snapshots and can be adapted to solve a variety of radiative processes. For example, codes such as hyperion (Robitaille, 2011) and radmc-3d (Dullemond, 2012) calculate the dust and molecular line emissivities. These data can then be used to generate realistic synthetic observations. Other MCRT codes include mocassin (Ercolano et al., 2003), which models photoionisation fronts and emission line intensities from Hii regions, and torus (Harries, 2011), which models time-dependent radiative transfer.
Smoothed particle hydrodynamics (SPH) (Lucy, 1977; Gingold & Monaghan, 1977) is a mesh-free method of solving the equations of fluid dynamics. A discrete set of particles are used to model the density distribution. These particles are smoothed over one another using a kernel function with smoothing length . By letting vary with density, the particles can model a field spanning many orders of magnitude in density. This property makes SPH an attractive scheme for modelling astrophysical systems such as star forming clouds (e.g. Bate, 2009; Lomax et al., 2014, 2015), galactic discs (e.g. Dobbs, Bonnell & Pringle, 2006; Dobbs, Burkert & Pringle, 2011) and the cosmic web (e.g. Springel & Hernquist, 2003). Open-source SPH codes include gadget-2 (Springel, 2005) and gandalf (Hubber & Rosotti, 2013).
Performing MCRT calculations on density fields from SPH simulations usually involves mapping the particles onto an octree (e.g. Stamatellos & Whitworth, 2005; Robitaille, 2011; Rundle et al., 2010). This fundamentally changes the structure of the density field by adding noise and/or merging several particles into a single cell (we provide a brief analysis of this noise in Appendix A). Another option is to use the particle positions to construct a Voronoi tessellation (e.g. Hubber, Ercolano & Dale, 2016). This is less noisy than an octree, although the implementation is more complicated. Nevertheless, both methods alter the SPH density structure by replacing smoothed particles with uniform density cells.
In this paper we present the new code spamcart (short for Smoothed PArticle Monte CArlo Radiative Transfer). The code performs MCRT calculations using the properties of the particles’ smoothing kernels. Therefore, post-processing calculations are performed at exactly the same resolution as the SPH simulation. The algorithm is similar to other ray-tracing and Monte Carlo methods used in SPH, (e.g. Altay, Croft & Pelupessy, 2008; Forgan & Rice, 2010). However, this is the first time such an algorithm has been used to perform full MCRT calculations. The method utilises the pre-existing neighbour-finding/gravity tree found in modern SPH codes, and is therefore relatively simple to implement.
In §2 we describe the general Monte Carlo method we use for radiative transfer. In §3 we explain how this method can be applied to smoothed particles. In §4 we present some example calculations, including a physical benchmark. In §5 we list some future implementations and in §6 we present the summary.
2 Monte Carlo radiative transfer
The method presented here is an adaptation of the Lucy (1999) algorithm, modified to operate on smoothed particles rather than a grid of uniform density cells. Most of the operations involve drawing a random variate from the Uniform distribution in the interval .
2.1 Luminosity packet trajectories
A radiation source with luminosity emits luminosity packets. Each packet is emitted from an origin in a direction . If the radiator is an isotropic point source, is given by polar angles
If there is an external isotropic radiation field, is a random point on a closed convex surface and is given by polar angles
where is the angle between and the inward surface normal. Each packet has energy 111If we assume radiative equilibrium, then the value of is arbitrary and does not affect the calculation. and wavelength , drawn randomly from the spectral energy distribution of the source. For the surface of a star, approximated as a blackbody, the equation
is solved for , where is the temperature of the source.
An individual packet travels an optical depth
This corresponds to a distance , where
Here, is the density of the medium and is the mass extinction coefficient. Solving Eqn. 2.5 for is one of the more computationally demanding aspects of MCRT and is covered in the next section.
Once a packet has travelled distance , it is either absorbed or scattered at position
A packet is absorbed if
where is the albedo of the medium. Otherwise the packet is scattered. An absorbed packet is immediately re-emitted with a new wavelength drawn from the emissivity distribution of the medium. This is achieved by solving the equation
for , where is the monochromatic mass absorption coefficient and is the temperature of the medium at . If the packet is scattered, remains unchanged and is scattered away from by a random angle . Here, is drawn from the Henyey & Greenstein (1941) phase function:
where is the mean scattering cosine. Hence, relative to the pre-scattering direction,
The process of absorption then re-emission and/or scattering is repeated until the packet exits the system.
2.2 Emissivity distribution
For a medium in radiative equilibrium, the energy emission rate is equal to the energy absorption rate , i.e.
Here, is the Stefan-Boltzmann constant, is the local temperature and is the Planck mean absorption coefficient:
The average value of within a volume can be estimated by summing over the path lengths of luminosity packets which pass through ,
We can estimate , and therefore the mass emissivity distribution, by calculating and solving Eqn. 2.11 for . However, Eqn. 2.13 depends on the trajectories of the luminosity packets and the trajectories depend on the emissivity distribution. Therefore several iterations are required to reach an equilibrium solution.
3 SPAMCART algorithm
The main algorithm of the spamcart code consists of two tasks: (i) to calculate the trajectories of luminosity packets through an ensemble of smoothed particles and (ii) to estimate for each particle. Here, we detail how these tasks are performed.
For an ensemble of smoothed particles, each with position , mass and smoothing length , the density at any point in space is given by
where is the kernel function. Most kernel functions have compact support, i.e. they are only finite within smoothing lengths. We use the M4 cubic spline kernel (Monaghan & Lattanzio, 1985) which has (see Appendix B). The smoothing length of each particle is calculated so that
These equations are solved by iteration and we normally adopt (as suggested by Price & Monaghan, 2004).
3.2 Scatter Calculation
The value of some arbitrary quantity at an arbitrary position can be estimated by scattering the same quantity, from each particle to , via the kernel function:
The gradient may also be estimated from the gradient of the kernel function:
The gradient of the M4 kernel function is given in Appendix B .
3.3 Optical depth
The set of all points along a ray is defined by , where . Here, is the origin, is the direction unit vector and is the length of the ray. The optical depth at wavelength along the ray is given by
In order to calculate Eqn. 3.5, we must first identify all the particles which are intersected by the ray (see Fig. 2). The sum of all particle column density contributions is then used to calculate the total optical depth. Particle-ray intersections can be identified efficiently by walking a tree-structure (see Fig. 3) and opening cells which pass the slab test (e.g. Williams et al., 2005, see Appendix C). These trees are a standard element of SPH codes, used to optimise neighbour-finding and gravity calculations.
3.4 Propagating luminosity packets
We propagate a luminosity packet, starting at , travelling in direction . Initially, the ray has length . Here is a modest overestimate of the intended value of (see Eqn 2.5). A good choice of is provided by a second order Taylor expansion at , i.e.
Here, the inverse mean-free-path is estimated via a scatter calculation and is an overestimation factor which we set to .
We use a modified Newton-Raphson root-finding method to iteratively solve for (see Appendix D). We terminate the iterations when , usually within four or five iterations. For some values of , is less than and no solution exists. Here, either the luminosity packet has left the ensemble of SPH particles or is too short. In the first case, the calculation for that luminosity packet is complete. In the second case, a new ray is constructed with the same direction and new origin . The calculation is then repeated with .
3.5 Estimating the energy absorption rate
The energy absorption rate for an individual particle is estimated by summing the column densities along all packet trajectories:
Here is the column density along the path of luminosity packet through the particle. Estimating requires almost no computational expense as the values of have already been computed whilst propagating the luminosity packets.
The absorption rate from the previous iteration is used when randomly generating a new emission wavelength. The value of at an arbitrary position can be found via a scatter calculation. The radiative equilibrium calculation is complete when the change in (or temperature ) between iterations is less than a desired tolerance.
3.6 Generating intensity maps
We generate intensity maps of the dust emission via ray tracing. We construct a virtual rectangular screen facing the ensemble of particles. This screen is divided into pixels. We construct rays, each of which has infinite length and passes through the centre of a pixel with direction normal to the screen. A schematic of a single pixel is show in Fig. 4 . All intensity maps presented here have .
Each ray intersects particles. These are sorted into descending order of distance from the pixel. The intensity of the pixel is given by . This is calculated by iteration:
where . In the case where , is the background intensity.
The intensity from point sources (e.g. stars) may also be added to the map. This is achieved by (i) locating the pixel in which the source lies, (ii) calculating the optical depth between the source and the pixel and (iii) adding to the pixel intensity. Here is the monochromatic luminosity of the source and is the area of the pixel.
4 Example calculations
In the following calculations, we use the dust optical properties derived by (Li & Draine, 2001). Here, the dust is a mixture of carbonaceous and amorphous silicate grains. The size distribution follows Weingartner & Draine (2001) with and is normalised to give a dust-to-gas mass ratio of one percent. The code does not currently support imaging of scattered light (see §5), so we set the albedo to zero, i.e. .
We present two examples. The first is a benchmark which tests the ability of the algorithm to reach thermal equilibrium with a radiation field. In the second example, we use the outputs of an SPH simulation to generate spectra and intensity maps of an embedded multiple system.
4.1 Blackbody radiation field
4.1.1 Undiluted blackbody radiation field
In an undiluted blackbody radiation field, the intensity is equal to the Planck function . By construction, an absorbing/emitting object in this field has a surface intensity equal to . Following Kirchhoff’s law of radiation—a good absorber is an equally good emitter—the surface temperature must therefore be equal to . Furthermore, because the emission spectrum is identical to the absorption spectrum, the object is invisible with respect to the background intensity.
4.1.2 Set up
We set up a sphere of particles. The sphere has radius and the density profile of a critical Bonnor-Ebert sphere, truncated at dimensionless radius . In the first instance, we give the sphere mass (optically thin), in the second (optically thick). We place the sphere in an undiluted blackbody radiation field. This is simulated by constructing a virtual shell which encompasses all the particles and their smoothing kernels, i.e. . The shell is given luminosity . This luminosity is divided into luminosity packets, which are directed inwards towards the particles. The wavelength of each packet is drawn randomly from the Planck function . Each particle is given an initial temperature and spamcart is iterated until the mean change in temperature is less than .
Fig. 5 shows the spectral energy distribution (SED) though the centre of the sphere. In the optically thin case, we see that the SED converges on a blackbody after about three iterations. In the optically thick case, the SED converges after about five iterations.
Fig. 6 shows the particle temperatures after the first three iterations of the algorithm with the optically thin sphere. This includes calculations with both and . When , the particle temperatures settle on a narrow distribution of . When the temperatures distribution is . This demonstrates the Poisson-like uncertainties on Monte Carlo calculations, i.e. a factor of fewer luminosity packets increases the signal-to-noise ratio by a factor of .
4.1.4 Diluted blackbody radiation field
We repeat the optically thick simulation, this time with a blackbody, diluted by a factor of . In this instance, we expect the dust to reach local radiative equilibrium at some temperature . Here, when a luminosity packet is absorbed by the dust, it is usually re-emitted at a longer wavelength. The sphere should therefore appear to glow, relative to the background, at long wavelengths and cast a silhouette at short wavelengths.
Fig. 7 shows intensity maps of the optically thick sphere. The left column shows the results of the undiluted blackbody field, the right column shows the results for the diluted field. The sphere in the undiluted field is almost invisible relative to the background intensity at wavelengths between and . There is some visible noise at , of order ten percent. Here, the optical depth through the sphere is very high () and the intensity is very sensitive to temperature fluctuations across the the sphere’s outer surface. The sphere in the diluted field behaves as expected; it is brighter than the background at and darker than the background at .
4.2 Embedded sextuple system
We demonstrate a more realistic application of the algorithm by calculating the dust emission from a protostellar multiple system. High order multiple systems, i.e. , are observed amongst mature field stars (e.g. Tokovinin, 2008; Eggleton & Tokovinin, 2008). These systems are more common amongst pre-Main Sequence stars, embedded in star forming regions such as Taurus and Ophiuchus (e.g. Leinert et al., 1993; Ratzka, Köhler & Leinert, 2005; Kraus et al., 2011). Furthermore, similar systems form routinely in simulations of molecular clouds and cores (e.g. Delgado-Donate et al., 2004; Bate, 2009; Lomax et al., 2014, 2015).
4.2.1 Set up
We take a snapshot from one of the SPH simulations presented by Lomax et al. (2014)222These simulations use initial conditions drawn from distribution functions that reproduce the observed properties of the cores in Ophiuchus. The simulation used here is No. 52, with episodic radiative feedback.. Here, a core collapses and fragments into seven protostars (represented by sink particles). One of these is ejected from the core and the remaining six form a stable sextuple system. At (after the initial collapse), the protostars in the sextuple have masses between . The ejected protostar has mass . The remainder of the initial mass (, ) is distributed between discs and and a diffuse envelope.
We irradiate the system with an undiluted blackbody radiation field and a blackbody field, diluted by a factor of . This represents contributions from the Cosmic Microwave Background (CMB) and the galactic stellar population. The protostars are treated as blackbody sources with and ; these temperatures are estimated using the Stamatellos, Whitworth & Hubber (2011) episodic accretion model.
4.2.2 Intensity maps and spectra
Fig. 8 shows column density and intensity maps of the embedded system. The system is viewed both face-on and edge-on. In the face-on view, we see a balanced quadruple system (top right) in orbit with a binary system (bottom left), separated by .
In both face-on and edge-on views, the intensity at and roughly traces the column density. At , face-on, only the the inner regions of the discs are visible. When viewed edge-on, the discs are opaque and significant intensity is only seen above and below their midplanes.
Fig. 9 shows the emission spectrum of the system. When viewed face on, we may define different wavebands, dominated by different sources of radiation. Between , the dust is optically thin and most of the radiation is from the CMB. Between , most of the radiation is from dust emission. Starlight is largely unattenuated in these wavebands, but its fraction of the total emission is very low. At wavelengths less than , the majority of the emission is attenuated starlight. Similar dust emission is seen when the system is viewed edge on. However, the starlight is almost completely extinguished at . We note that the level of extinction at is probably exaggerated as we have not captured scattered light in these spectra. Nevertheless, in the face-on view, we still see the double peaked spectra typical of protostellar discs (e.g. Robitaille, 2011; Whitney et al., 2013).
5 Future developments
While the code is not yet heavily optimized, the algorithm is efficient enough to run on an Intel i5-4200U dual core CPU at in a reasonable time. A single iteration of the optically thick sphere in §4.1.1 requires per packet per CPU core. The optically thin sphere requires per packet per CPU core. A single iteration of the sextuple system in §4.2 requires approximately per packet per CPU core. We note that the optimal number of packets is not obvious. In this paper, we select sensible values through a process of trial and error.
A thorough analysis of the processing time required for this algorithm in different scenarios is beyond the scope of this paper. However, we can estimate the time scaling for simple systems such as a star embedded in a uniform density sphere. Here,
The flight path of a packet is proportional to the square of the average optical depth through the system . The number of particles along the path is proportional to the square of the kernel extent and the cube root of the total number of particles . The time taken to find each particles scales roughly with the log .
One may ensure that each particle is visited the roughly same number of times, regardless of , by setting . Ignoring , the time scaling is now . This is similar to that of an SPH simulation timestep. Therefore it may be feasible in the future to run this algorithm on-the-fly in an SPH simulation, especially in circumstances where .
At present, we are able to generate intensity maps of the dust emission, plus attenuated light from background and point sources. However, it is desirable to also capture light scattered by dust. This can be achieved by modifying Eqn. 3.9. Here, the term represents the amount of dust emission per unit solid angle, per unit wavelength, per unit mass. We may add an extra term which accounts for scattered light:
Here, is the mass scattering coefficient and is the angle between luminosity packet and the viewing angle. The sum is only over luminosity packets with wavelengths in the interval . The scattering term must be calculated during the Monte Carlo iteration, so the viewing angle and wavelength interval must be known beforehand.
We have developed a method for performing MCRT calculations directly on a distribution of SPH particles The algorithm operates differently from uniform-density cell methods, but the two schemes are mathematically equivalent. This allows an MCRT calculation to be performed on an SPH snapshot with (i) no loss in density resolution and (ii) no introduction of noise from mapping particles to cells.
We present a version of this algorithm that uses the Lucy (1999) method to compute (i) the propagation of luminosity packets through a medium and (ii) the radiative equilibrium temperature. The trajectories of the packets rely on the temperature of the medium so the calculation must be solved by iteration.
We provide two example calculations using the smoothed particle MCRT method. First, we show that a cloud is invisible when it is bathed in an undiluted blackbody radiation field. This holds for both optically thin and optically thick cases. Therefore the code obeys Kirchhoff’s law of thermal radiation. We note that this is a powerful test of any radiative transfer code and can be applied to any configuration. Conversely, a cloud which has reached radiative equilibrium with a diluted blackbody field glows at long wavelengths and casts a silhouette at short wavelengths. Second, we generate intensity maps and spectra of protostellar and dust emission from an embedded sextuple system. Here, a double peaked disc spectrum is seen when the system is observed face on. When viewed edge on, the opacity of the disc blocks nearly all of the starlight.
Future additions to the code include the addition of scattered light to synthetic observations and optimisations for optically thick regions of dust. The code, written in Fortran 2003/08 with OpenMP parallelisation, will be made publicly available in the near future.
OL and APW gratefully acknowledge the support of a consolidated grant (ST/K00926/1) from the UK STFC. We also thank the referee for their considerate and constructive comments.
- Altay, Croft & Pelupessy (2008) Altay G., Croft R. A. C., Pelupessy I., 2008, ACP08, 386, 1931
- Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
- Bate (2009) Bate M. R., 2009, MNRAS, 392, 590
- Delgado-Donate et al. (2004) Delgado-Donate E. J., Clarke C. J., Bate M. R., Hodgkin S. T., 2004, MNRAS, 351, 617
- Dobbs, Bonnell & Pringle (2006) Dobbs C. L., Bonnell I. A., Pringle J. E., 2006, MNRAS, 371, 1663
- Dobbs, Burkert & Pringle (2011) Dobbs C. L., Burkert A., Pringle J. E., 2011, MNRAS, 417, 1318
- Dullemond (2012) Dullemond C. P., 2012, RADMC-3D: A multi-purpose radiative transfer tool. Astrophysics Source Code Library, record ascl:1202.015
- Eggleton & Tokovinin (2008) Eggleton P. P., Tokovinin A. A., 2008, MNRAS, 389, 869
- Ercolano et al. (2003) Ercolano B., Barlow M. J., Storey P. J., Liu X.-W., 2003, MNRAS, 340, 1136
- Forgan & Rice (2010) Forgan D., Rice K., 2010, MNRAS, 406, 2549
- Gingold & Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
- Harries (2011) Harries T. J., 2011, MNRAS, 416, 1500
- Henyey & Greenstein (1941) Henyey L. G., Greenstein J. L., 1941, ApJ, 93, 70
- Hubber, Ercolano & Dale (2016) Hubber D. A., Ercolano B., Dale J., 2016, MNRAS, 456, 756
- Hubber & Rosotti (2013) Hubber D. A., Rosotti G. P., 2013, GANDALF (Graphical Astrophysics code for N-body Dynamics And Lagrangian Fluids). GitHub
- Kraus et al. (2011) Kraus A. L., Ireland M. J., Martinache F., Hillenbrand L. A., 2011, ApJ, 731, 8
- Leinert et al. (1993) Leinert C., Zinnecker H., Weitzel N., Christou J., Ridgway S. T., Jameson R., Haas M., Lenzen R., 1993, A&A, 278, 129
- Li & Draine (2001) Li A., Draine B. T., 2001, ApJ, 554, 778
- Lomax et al. (2014) Lomax O., Whitworth A. P., Hubber D. A., Stamatellos D., Walch S., 2014, MNRAS, 439, 3039
- Lomax et al. (2015) Lomax O., Whitworth A. P., Hubber D. A., Stamatellos D., Walch S., 2015, MNRAS, 447, 1550
- Lucy (1977) Lucy L. B., 1977, ApJ, 82, 1013
- Lucy (1999) Lucy L. B., 1999, A&A, 344, 282
- Min et al. (2009) Min M., Dullemond C. P., Dominik C., de Koter A., Hovenier J. W., 2009, A&A, 497, 155
- Monaghan & Lattanzio (1985) Monaghan J. J., Lattanzio J. C., 1985, A&A, 149, 135
- Price & Monaghan (2004) Price D. J., Monaghan J. J., 2004, MNRAS, 348, 139
- Ratzka, Köhler & Leinert (2005) Ratzka T., Köhler R., Leinert C., 2005, A&A, 437, 611
- Robitaille (2011) Robitaille T. P., 2011, A&A, 536, A79
- Rundle et al. (2010) Rundle D., Harries T. J., Acreman D. M., Bate M. R., 2010, MNRAS, 407, 986
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 312
- Stamatellos & Whitworth (2005) Stamatellos D., Whitworth A. P., 2005, A&A, 439, 153
- Stamatellos, Whitworth & Hubber (2011) Stamatellos D., Whitworth A. P., Hubber D. A., 2011, ApJ, 730, 32
- Tokovinin (2008) Tokovinin A., 2008, MNRAS, 389, 925
- Weingartner & Draine (2001) Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296
- Whitney et al. (2013) Whitney B. A., Robitaille T. P., Bjorkman J. E., Dong R., Wolff M. J., Wood K., Honor J., 2013, ApJS, 207, 30
- Williams et al. (2005) Williams A., Barrus S., Morley R. K., Shirley P., 2005, J. Graphics Tools, 10, 49
Appendix A Gridding Errors
Transposing an ensemble of SPH particles onto a grid introduces errors to the density field. Here, we perform a brief analysis on these errors. We generate an octree around the distribution of particles used to model the sextuple system in §4.2. First, we build the smallest cuboidal root cell that contains all of the particles. This cell is recursively subdivided into eight equal volume cells until each leaf cell contains particles or fewer. Fig. 10 shows the - projection of the particles and the octree with .
The density of each cell is computed using an SPH scatter calculation (see Eqn. 3.4) at the centre of the cell. We compare this with the particle density (scatter calculation) of each particle in the cell. Fig. 11 shows the relative difference between and for the ensemble of particles with different values of . Note that the particles are sorted into ascending order of density error. For comparison, we also include the relative difference between the particle scatter density and its gather density (see Eqn. 3.2)333The difference between the scatter and gather density calculations does not indicate that SPH is somehow wrong. It reflects the difference between the mass to volume ratio of an SPH particle (gather), and the superposition of multiple smoothing kernels at a given point in space (scatter).. In Table 1 we give some further statistics including the number of leaf cells, the maximum depth of the tree, the total mass of the grid (defined as the sum of the products of leaf cell volume and leaf cell density) and the 90th centile density error.
Gridding the density field with 1, 8 and 64 incurs a 90th centile density error of 24%, 45% and 76% respectively. Furthermore, the mass of the grid is roughly 80% to 90% that of the particle ensemble. We note that although the method presented in the paper is more accurate than gridding an ensemble of particles, MCRT calculations on octrees are mathematically and computationally simpler.
Appendix B Kernel function
The M4 kernel has support , where . The density at is given by
and the density gradient is given by
The column density from dimensionless impact parameter to is given by
In practice, Eqn. B.3 can either be computed on-the-fly or stored in a triangular lookup table in the the range and .
Appendix C k-d tree
We use a -d tree (in three dimensions) to find the particles intersected by a ray. We construct the tree by placing the entire ensemble of particles into a root cell. The root cell contains two branch cells. We select the dimension across which the particle position has the greatest variance. All particles with up to and including the median value are placed in the first cell and the remaining particles are placed in the second cell. This process is repeated recursively until each cell contains 8 or fewer particles, in which case they are leaf cells. For each cell, we calculate the Axis-Aligned Bounding Box (AABB) that encompasses all of the particle smoothing volumes within the cell.
We find the particles intersected by a ray by using the slab method (e.g. Williams et al., 2005). Consider a ray with origin , direction and length . Also, consider an AABB with lower limits and upper limits . Along each dimension , we define a slab (a slab is a space between two parallel planes) with lower limit and upper limit . The length of the ray segment, -, within the slab can be calculated:
The length of the ray through the AABB in all three dimensions, , is given by
The ray intersects the AABB if the following statements are true:
Starting at the root cell, we check for an AABB-ray intersection. If the result is true, we open the branch cells and recursively repeat the process until we encounter a leaf cell. If the leaf cell is intersected by the ray, the cell’s contents are added to a particle list. This method is not exclusive to -d trees and may be used with any axis aligned tree, such as an octree (e.g. Barnes & Hut, 1986).
Appendix D Root-finding method
Newton’s method can be used to find the length of a luminosity packet trajectory. The equation , where is the column density, can be solved rapidly by iterating
where . This performs well if is near the root of , but may fail to converge for poor initial choices of .
We modify this method so that convergence is guaranteed. First we bracket the root, i.e. find values and so that and are on opposite sides of . We assume that Eqn. 3.7 succeeds in overestimating and set and . We now define the parabolic curve which (i) passes through point with gradient and (ii) passes through point . This curve is described by the equation
and by construction must have a single root in the interval . Rearranging for the second derivative,
Solving the quadratic formula, the iteration function is now
If and are on opposite sides of , the new is set to the old . Otherwise remains the same.