SPAMCART: a code for smoothed particle Monte Carlo radiative transfer

O. Lomax, & A. P. Whitworth
School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, UK

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.

radiative transfer – hydrodynamics – methods: numerical – ISM: dust
pagerange: SPAMCART: a code for smoothed particle Monte Carlo radiative transferDpubyear: 2015

1 Introduction

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.

3.1 Density

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

Figure 1: A diagram of line with origin , direction and length intersecting a particle at position with compact support . The other terms are defined in §3.3 .

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




and is the mass extinction coefficient of particle . A diagram of this system for a single particle is given in Fig. 1 . The analytical form of for the M4 kernel in given in Appendix B.

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.

Figure 2: A diagram of all particles which are intersected by a ray of length . The optical depth along the ray may be calculated using Eqn. 3.5 . Note that is greater than . This is because is calculated by iteration and the first iterate must be an overestimate. Also, because of the construction of Eqn. 3.5, particle column density beyond does not contribute to the total optical depth.

3.4 Propagating luminosity packets

Figure 3: A schematic of the cells of a particle tree. Each cell has an axis-aligned bounding box (AABB) which encompasses the smoothing volume of all the particles in the cell (see Appendix C for more details). If the ray intersects the AABB of a cell with sub-cells, the sub-cells are recursively checked. If the cell has no sub-cells (shaded pink), all of its particles are added to the ray.

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

Figure 4: Schematic of a single pixel (small square) of an intensity map (large square). The ray passes through the centre of the pixel at an angle normal to the surface of the map.

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

Figure 5: The spectral energy distribution (SED) through the centre of the Bonnor-Ebert spheres, with (top) and (bottom). In both cases, the sphere is illuminated by an undiluted blackbody. The left -axis gives the intensity calculated by successive iterations of the algorithm. The lowest red line shows the initial intensity, where each particle has a temperature of . The lines above show the SED converging towards a blackbody (black long-dashed line). The right -axis give the optical depth through the centre of the sphere (short-dashed black line).
Figure 6: Particle temperature versus radius for optically thin sphere. The top frame shows the temperature calculation with luminosity packets. The bottom frame shows the calculation repeated with luminosity packets. The particle temperatures from the first iteration are plotted as red squares, from the second iteration as green circles and from the third iteration as blue triangles.

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 .

4.1.3 Results

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 .

Figure 7: Intensity maps of a sphere of particles with and . The radial density profile follows that of a critical Bonnor-Ebert sphere. The left column shows the intensity when the sphere is illuminated by a blackbody radiation field. The right column shows the intensity when the sphere is illuminated by a blackbody radiation field, diluted by a factor of . The top row gives the intensity at , the middle row gives the intensity at and the top row gives the intensity at . The colour scale gives the intensity, normalised by the background intensity.

4.2 Embedded sextuple system

Figure 8: Intensity maps of dust emission from a simulated protostellar sextuple system. The left column shows a face-on view, the right column shows shows an edge-on view. The top row shows the column density of the system. The following three rows show the intensity at , and . Each frame has area . Note: the colour bar applies to the right column. Figures in the left column have been linearly scaled to fit in the same range.
Figure 9: Spectra of dust emission and starlight from the sextuple system. Values are integrated over the map area shown in Fig. 8 . The top frame shows the spectrum from a face-on view. The bottom frame shows the spectrum from an edge-on view. The solid red lines show the total emission. The dashed red line shows the emergent emission from the six protostars. The dashed black line shows the protostellar emission in the absence of dust extinction.

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

5.1 Performance

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  .

5.2 Features

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.

Other future developments include use of the partial diffusion approximation and modified random walk for regions of high optical depth (Min et al., 2009). The spamcart code will shortly be made open-source and uploaded to the GitHub repository.

6 Summary

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.

Figure 10: An - projection of the octree and particle ensemble with . The root cell is roughly cubic with an edge length of . The smallest leaf cells have an edge-length of .
Figure 11: The relative difference between the grid cell density and the SPH particle density against particle number. Particles are sorted in order of ascending error. From top to bottom, the solid coloured lines show , and . The dotted black line shows the relative difference between the gather and scatter density for each particle, i.e. .
1 154281 16 0.80 0.24
8 22632 13 0.87 0.45
64 3242 12 0.88 0.76
Gather 0.18
Table 1: Statistics of the particle octrees. The first column give the maximum number of particles per leaf. The second column gives the total number of leaf cells. The third column gives the number of cells in an equivalent regular grid. The forth column gives the ratio of the total grid mass to the total particle mass. The fifth column gives the 90th centile density error (see Fig. 11).

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


where and


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.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description