Native Synthetic Imaging of Smoothed Particle Hydrodynamics density fields using gridless Monte Carlo Radiative Transfer
aa r X i v : . [ a s t r o - ph . E P ] A p r Mon. Not. R. Astron. Soc. , 1–10 (0000) Printed 30 October 2018 (MN L A TEX style file v2.2)
Native Synthetic Imaging of Smoothed Particle Hydrodynamicsdensity fields using gridless Monte Carlo Radiative Transfer
Duncan Forgan ⋆ and Ken Rice Scottish Universities Physics Alliance (SUPA), Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, Scotland, UK
Accepted 2010 April 13. Received 2010 April 10; in original form 2010 January 29
ABSTRACT
An algorithm for creating synthetic telescope images of Smoothed Particle Hydrodynamics(SPH) density fields is presented, which utilises the adaptive nature of the SPH formalismin full. The imaging process uses Monte Carlo Radiative Transfer (MCRT) methods to modelthe scattering and absorption of photon packets in the density field, which then exit the systemand are captured on a pixelated image plane, creating a 2D image (or a 3D datacube, if thephotons are also binned by their wavelength). The algorithm is implemented on the densityfield directly: no gridding of the field is required, allowing the density field to be described toan identical level of accuracy as the simulations that generated it.Some applications of the method to star and planet formation simulations are presentedto illustrate the advantages of this new technique, and suggestions as to how this frameworkcould support a Radiative Equilibrium algorithm are also given as an indication for futuredevelopment.
Key words: methods: numerical, observational, scattering, stars: imaging, radiative transfer
The physics of dust is of paramount importance when modellingstar and planet formation: dust signatures are seen in both molec-ular clouds and circumstellar discs, and are a crucial componentin their observed properties over a wide range of wavelengths. In-deed, the existence of dust is essential if planets are to be formedinside these discs. The presence of dust has important effects on theradiative signatures that can be detected by astronomers: dust willscatter and polarise at short wavelengths, as well as reprocessingthis radiation to longer wavelengths. The nature of the scatteringand polarisation will depend on the geometry of both the systemand the dust in the system, as well as the physical properties of thedust itself (composition, morphology, grain size, etc).In recent years, telescopes/networks such as HST, SCUBA,MERLIN and Spitzer have efficiently probed stellar systems fromthe IR to the radio, with future ground and space-based mis-sions such as Herschel, ALMA, SCUBA II, JWST, SPICA and e-MERLIN improving the quality of this data. Astronomers will soonbe faced with a wealth of new, high-fidelity astronomical data onstellar objects across a wide spectral range, allowing detailed stud-ies of the evolution of dusty disc systems, in particular the interplaybetween disc dust and disc gas. For numerical simulations to informthese observations, the simulation of radiative transfer (RT) in these ⋆ E-mail: [email protected] systems must be pursued, so that imaging of numerical simulationscan provide theoretical insights to observations.Monte Carlo Radiative Transfer (MCRT) has become apopular tool in simulating and imaging astrophysical systems(Whitney & Hartmann 1992; Wood et al. 1996; Pinte et al. 2006).The technique involves the tracking of photon packets in themedium, allowing them to scatter and absorb in the medium (aswell as attaining a non-zero polarisation). The photons are trackeduntil they escape the medium, and can then be captured on an im-age plane. Traditionally, the method requires the density field to bedefined to as high an accuracy as possible - this is usually achievedby gridding the density field in 3 dimensions.Smoothed Particle Hydrodynamics (SPH) is a Lagrangianmethod which represents a fluid by a particle distribution (Lucy1977; Gingold & Monaghan 1977). Each particle is assigned amass, position, internal energy and velocity. From this, state vari-ables such as density can then be calculated at any position inthe system by interpolation - see reviews by Monaghan (1992,2005). It has been successful in modelling astrophysical sys-tems of various scales and geometries, from protostellar sys-tems to galactic systems and beyond. Its key advantage is theability to follow the change of density through many orders ofmagnitude adaptively. Gravity can be included in SPH calcula-tions, and optimised using traditional hierarchical tree methods(Hernquist & Katz 1989). The formalism can also be modified to c (cid:13) Duncan Forgan and Ken Rice simulate magneto-hydrodynamics (MHD) (Hosking & Whitworth2004; Price & Monaghan 2004, 2005; Price & Bate 2007).Radiative physics is crucial to the simulation of any astro-physical fluid. Radiative transfer in SPH has a long history, withefforts ranging from simple parametrisations for the radiative cool-ing time (e.g. Rice et al. 2003), optical-depth dependent radiativecooling (e.g. Stamatellos et al. 2007), through to flux-limited dif-fusion models (e.g. Whitehouse & Bate 2004; Bastien et al. 2004;Viau et al. 2006; Whitehouse & Bate 2006; Mayer et al. 2007;Forgan et al. 2009). While these algorithms are extremely well-suited to calculating gas temperatures during runtime, they lack theprecision and insight offered by radiative transfer techniques whichdo not average over frequency, such as MCRT.MCRT techniques can be applied to the SPH densityfield, and have been on several occasions. Early attempts be-gan by binning the particle distribution onto a grid - how-ever, the choice of geometry strongly influences the finalgridded field, and often adaptive meshes are required tocorrectly represent the matter distribution (Oxley & Woolfson2003; Kurosawa et al. 2004; Stamatellos & Whitworth 2005;Bisbas et al. 2009). Later efforts have utilised ray-tracing directlyin SPH fields (Kessel-Deynet & Burkert 2000; Altay et al. 2008;Pawlik & Schaye 2008), allowing the full power of the SPH for-malism to be applied to the calculation of optical depths. However,these methods currently only calculate optical depths along the rayfor the purpose of photoionisation, etc, and do not account for de-tailed scattering and polarisation, which are important in imagingsmall-scale systems such as protoplanetary discs.This paper introduces an algorithm for imaging SPH simu-lations of star and planet formation directly using MCRT, with-out requiring gridding, and including scattering and polarisation. Itutilises the same techniques that a traditional gridded MCRT codeuses, with the advantage that it can trace rays in the density fieldwith the same adaptive capability as SPH, such that it can modelradiative effects with at least the same resolving power. The paperis organised as follows: section 2 will outline the algorithm, sec-tions 3 and 4 will describe some applications of the technique toimaging numerical simulations, and section 5 will summarise theresults of this work.
For completeness, a summary of MCRT is given here.
Photon pack-ets are emitted from sources (these can be point sources or diffuseemission from the density field itself). These packets (hereafter re-ferred to simply as photons) then travel through the medium, andinteract with the medium stochastically, reproducing the statisticalscattering properties of the medium. The exact formalism is sum-marised below.
Photons can be emitted either from a central source (e.g. a pro-tostar), or as diffuse emission from the density field itself, whichis only done if the temperature structure of the diffuse componentis known (e.g. as output from a simulation).
Radiative equilibrium methods in MCRT (e.g. Bjorkman & Wood 2001) can calculate thetemperature structure directly from any luminosity source (e.g. stel-lar, accretion, external radiation fields). However, as the SPH sim- ulations used in this work produce temperatures for each SPH par-ticle self-consistently, the system is assumed to already be in tem-perature equilibrium. It is also assumed that the dust is thermallycoupled to the gas, and that T dust = T radiation . This is suitablefor most purposes, except where significant stellar irradiation dom-inates the radiation field (not modelled in the simulations describedin this paper). However, as has been already said, radiative equilib-rium techniques can be used to sidestep this issue: an implemen-tation of radiative equilibrium in SPH fields is discussed in a latersection.It is reasonable to (initially) assume all objects emit accordingto a blackbody spectrum, which then gives: L star = 4 π R s Z ν max ν min B ν ( T s ) dν (1) L gas = M gas Z (cid:20)Z ν max ν min ǫ ν ( r ) dν (cid:21) d r , (2)for source and diffuse emission respectively (given the frequencyrange [ ν min , ν max ] of interest), where ǫ ν ( r ) is the emissivity in-terpolated from the SPH particle field, i.e. ǫ ν ( r ) = 4 π X j κ ν B ν ( T j ) m j ρ j W ( r − r j , h ) , (3)where the sum over j indicates a sum over nearest neighbours, κ ν isthe absorptive opacity, T j is the temperature of each SPH particle, m j and ρ j are the masses and densities respectively, and W isthe smoothing kernel (see section 2.2.1 for more detail on the SPHformalism).There is an extra subtlety regarding the frequency distributionof photons emitted from the gas. The frequency distribution willdepend on the local emissivity, which as shown above is an inter-polative sum. While in theory the full form of equation (2) shouldbe used to calculate gas luminosity (and the frequency distributionof the photons it emits), this paper approximates the sum using thecontribution from one particle only, i.e. L gas = X i L i = X i πM i Z ν max ν min κ ν B ν ( T i ) dν. (4)This saves computational expense, and is sufficiently accurate forthe examples in this paper, bearing in mind a) that the approxima-tion breaks down only in the inner regions of the disc, which arealready under-resolved for other reasons (see section 4.2) and b)the effective resolution of the images is too low for this effect to besignificant.The luminosity of each object (whether pointmass or SPH par-ticle) defines how many photon packets are emitted from that objectusing N γ,object = N γ,tot (cid:18) L object L tot (cid:19) (5) In order to correctly reproduce the scattering of photons in themedium, the cumulative distribution function (CDF) of interaction F ( τ ) = 1 − e − τ (6)must be reproduced. In essence, this implies sampling an opticaldepth using τ scatter = − − log(1 − ζ ) (7) c (cid:13) , 1–10 ative Synthetic Imaging of SPH density fields using MCRT where ζ is a random number between 0 and 1. The photon willtravel a distance L along a constant vector path (hereafter ray ), suchthat τ scatter = Z L ρ χ ν dℓ, (8)(where χ ν is the total opacity), and then be either absorbed or scat-tered. The majority of computation required in any MCRT code isthe solution of equation (8). As the density field of an astrophysicalsystem is in general non-trivial, analytic solutions cannot be ap-plied, and numerical procedures must be used. The algorithm mustbe able to calculate the optical depth along the ray for arbitrarydistances - this requirement is what demands the most CPU time.To simplify matters, most MCRT codes use gridding as a means ofdefining the density field simply in space.Once the photon has reached the scattering location, its direc-tion after scattering must be calculated. The angle of scattering isdefined by the scattering matrix M which acts on the Stokes vec-tor S = ( I, Q, U, V ) , where I is the intensity, ( Q, U ) are the linearpolarisations at 45 ◦ to each other, and V is the circular polarisation: S ′ = R ( π − i ) M R ( − i ) S. (9)The R matrices are Mueller matrices, which describe rotations toand from the observer’s frame. They are defined as: R ( ψ ) = ψ sin 2 ψ − sin 2 ψ cos 2 ψ
00 0 0 1 (10)The scattering matrix M is dependent on the dominant source ofscattering in the medium. It can be expressed as a function of sev-eral scattering parameters, M i : M (Θ) = a M M M M M − M M M (11)The scattering angle Θ and azimuthal angle φ are sampled ran-domly from the scattering matrix, using the cumulative distributionfunctions: F (Θ) = R Θ0 M sin Θ ′ d Θ ′ R π M sin Θ ′ d Θ ′ (12) F Θ ( φ ) = 12 π (cid:18) φ − (cid:18) M − M M + M (cid:19) P φ (cid:19) (13)Where P = p Q + U I (14)In general, whether the photon is absorbed or scattered, the totalnumber of photons is conserved by forcing absorbed photons tobe immediately re-emitted. As MCRT deals with photon packets ,energy can be conserved by reducing the energy of each photonpacket after a scattering/absorption event by a factor equal to thelocal albedo. This equates with the concept of a fraction of the pho-tons in the packet being absorbed by the medium, and the comple-mentary fraction being scattered. Other methods can also be used,e.g. “killing” a photon if the local albedo is less than a randomlysampled value, which saves computing emission from low-intensitypackets. However, this work uses the former method rather than thelatter. In the following work, the absorbed photons do not affect the Figure 1.
Defining an image plane . local temperature structure - it is assumed that the equilibrium tem-perature is known (having been produced by the SPH simulation).
When photons exit the medium, they are captured on an imageplane, oriented at user-specified angles θ v , φ v to the system, at afixed distance d (see Figure 1). They are then binned by their ( x, y ) coordinates on this plane to provide a pixelated image, averagedover solid angle (analogous to imaging in a CCD). If spectra areof interest to the user, then these can also be obtained by binningin λ (or indeed, an entire datacube can be obtained by binning inall three). “Classic” MCRT methods do not specify a single viewingangle, and instead bin the photons over several lines of sight, whichallow the construction of a series of image planes from one simu-lation. However, such “multi-plane” simulations will have to runmany more photons than a “one-plane” simulation, to maintain acomparably low level of random error associated with each pixel .For this reason the code used in this work adopts the “one-plane”imaging scheme.This numerical apparatus requires the code to be run once forevery desired orientation of the image plane. It is possible to avoidthis issue by defining a series of image planes around the system,but the numerical resolution of the Monte Carlo method is reducedas the pixel to photon ratio is increased, so this is in general avoidedto prevent significant computational expense.The image plane size is defined in physical units, e.g. an im-age that captures features a distance of r au from the centre ofthe system. This prevents photons from features outside r (or pho-tons which escape with a vector which does not intersect the plane)from being recorded in the final image. Having specified r and thedistance d , this sets the angular size of the image. Coupled withinformation regarding the angular resolution of the telescope be-ing simulated, this defines how well resolved the image is, and thecalculated flux. In lieu of more sophisticated telescope simulation,the image is simply smoothed pixel by pixel with a Gaussian withwidth equal to the angular resolution. While obviously far short ofa true synthetic image, it remains a useful first-order approximationfor this work. This is also a factor when comparing one-plane simulations with differentpixel resolutions.c (cid:13) , 1–10
Duncan Forgan and Ken Rice
As has been said previously, the key component of any MCRT code(and the source of the greatest computational burden) is the calcu-lation of optical depths. This calculation requires a specification ofthe density field at all locations, which can be given by the SPH for-malism. The definition of a scalar field A ( r ) using an SPH particledistribution is A ( r ) = X j A j m j ρ j W ( r − r j , h ) (15)where A j is the value of A at the location of particle j , m j is themass of particle j , ρ j is the density of particle j , and W is thesmoothing kernel (or interpolating kernel). The summation is typi-cally over N neigh nearest neighbour particles to the location r . The smoothing length of particle j , h j , is defined such that a sphere ofradius h j will contain the N neigh nearest neighbour particles to j . For example, to calculate density, substitute A = ρ : ρ ( r ) = X j m j W ( r − r j , h ) . (16)The sphere that contains the N neigh nearest neighbours (i.e. asphere of radius h j ) is referred to as the smoothing volume .Thereis a subtlety to equation (16) that should be noted, relating to whichvalue of h to use. There are now two means by which to esti-mate density: the first is the so-called “gather” method, where thesmoothing length h = h i is defined for the location r i , and ρ ( r i ) = X j m j W ( r i − r j , h i ) (17)Where the index j refers to all particles which are contained withina radius h i of the location r i .The second method (which is used in this work and in SPHRAY
Altay et al. 2008) is the “scatter” method. The smoothing length h = h j is used - in this formalism, the density at any one loca-tion is calculated by adding the contributions from particles whosesmoothing volume intersects the location: ρ ( r i ) = X j m j W ( r i − r j , h j ) . (18)In the context of ray tracing, the density along the ray is affectedonly by particles with smoothing volumes that intersect it (seeFigure 2). By determining which particles intersect the ray, the restof the particle distribution can be ignored for the purposes of cal-culating optical depth, reducing computational expense (whereaswith the gather method, the ensemble of particles contributing tothe calculation changes significantly with position, and requiresthe inclusion of a larger subset of SPH particles to perform thecalculation).Using the scatter method, the column density Σ along the ray is Σ = Z L ρ ( r ) dℓ = Z L N X j =1 [ m j W ( r i − r j , h j )] dℓ (19)Which can be rearranged to give Σ = N X j =1 (cid:20)Z L m j W ( r i − r j , h j ) dℓ (cid:21) (20) Figure 2.
Illustrating the “scatter” method. Particles only contribute to thedensity along the ray if their smoothing volume intersects it.
Figure 3.
Optical Depths through a single smoothing volume .
The integral is now decomposed into N integrals, where N is thenumber of particles intersected by the ray. Each integral is definedby the impact parameter b (see Figure 2). The calculation itselfcan be performed for a smoothing volume of h = 1 , and scaledupwards (this is due to the construction of the smoothing kernel).The entire optical depth calculation has been decomposed into therepetition of a single algorithm for calculating the optical depththrough a single smoothing volume. This calculation will now beexpounded. Consider Figure 3. The ray (with direction vector n ) intersects thesphere with impact parameter b . If the ray penetrates a distance x into the sphere (out of a total possible distance s ), then the inte-gral can be defined analytically, given the functional form of W .defining ˜ r = r/h (21) ˜ b = b/h (22) ˜ x = x/h (23) r f = p ( s/ − x ) + b (24) c (cid:13) , 1–10 ative Synthetic Imaging of SPH density fields using MCRT Figure 4.
Schematic of an Axis Aligned Bounding Box (AABB). For agiven set of occupants of a tree cell (solid line), the AABB (dashed line)is set to the minimum dimensions required to completely contain theirsmoothing volumes. (and ˜ r f = r f /h ), the integral in equation (20) becomes: I = ( R r f W (˜ r ′ ) d ˜ r ′ x < s/ R b W (˜ r ′ ) d ˜ r ′ + R ˜ r f ˜ b W (˜ r ′ ) d ˜ r ′ x > s/ (25)The column density through a smoothing volume (for any impactparameter and any distance into the sphere) can now be calculated;multiplying by an opacity then provides an optical depth. Typically, W is constructed using cubic splines (Monaghan 1992). The kernelused in this work is W (˜ r ) = − ˜ r + ˜ r ˜ r < (2 − ˜ r ) < ˜ r <
20 ˜ r > (26)This provides compact support (i.e. it reduces to zero outside thesmoothing volume), and is simple to integrate. With a prescription for calculating optical depth for a single spherein place, a scheme for calculating ray/sphere intersections must beconstructed. To this end, the code creates a data object called a raylist , which stores (in order of intersection) all particles that theray (given its origin and direction vector) will intersect. Once thelist is created, the optical depth can be calculated quickly usingequation (18).The construction of the raylist must be computationally effi-cient for the code to be effective. The procedure is similar to thatimplemented by Altay et al. (2008) in
SPHRAY ; the code constructsan octree to spatially index the particles efficiently (as there maybe density changes over several orders of magnitude). The cellseither contain child cells, or particles (the leaf cells ). The tree isconstrained to have a maximum number of particles in each leaf.All cells have an associated Axis Aligned Bounding Box (AABB),which is the minimum box size, aligned to the three cartesian axes,to contain all the smoothing volumes of the particles in the cell (seeFigure 4). These AABBs are necessary as tree nodes may contain aparticle, but not its entire smoothing volumeThis allows the determination of intersections between the rayand the cells (or more correctly, their AABBs). Starting with theroot cell, each child cell is tested for intersection, constituting awalk through the tree. If a leaf cell is intersected by the ray, then the particles in the leaf are tested for intersection (by calculatingtheir impact parameters). This ensures that only a minimum frac-tion of the particles in the system need testing for intersection. Thisillustrates the necessity of AABBs; tree nodes may contain a par-ticle, but not its entire smoothing volume. Thus, calculating inter-sections between a ray and tree nodes may miss contributions tothe density field from smoothing volumes that cross node intersec-tions. Tests for intersections between rays and AABBs are carriedout using the ray slopes algorithm (Eisemann et al. 2007), whichhas been shown to be faster than other commonly used methods,such as using Pl¨ucker coordinates (Mahovsky & Wyvill 2004).
An important facet of an MCRT code is the determination of thescattering location of the photon. In general, the scattering locationwill occur inside a smoothing volume, and possibly at a locationwhere the density depends on the contributions from several parti-cles. Therefore, when attempting to determine the scattering loca-tion, it is important to define four classes of particle:(i) Particles that do not intersect the ray ( unlisted )(ii) Particles that intersect, but do not contain the location ofemission ( distant-listed )(iii) Particles that intersect, and contain the location of emissionin front of them ( front-listed )(iv) Particles that intersect, and contain the location of emssionbehind them ( back-listed )The classes are illustrated in Figure 5. Particles of class (i) obvi-ously do not affect the calculation - particles of class (ii) are ac-counted for simply. Particles of classes (iii) and (iv) will have dif-fering effects on the optical depth calculation, and will require sep-arate treatments.The scattering location is determined by iteration: firstly, theoptical depth is calculated particle by particle using the raylist un-til the optical depth exceeds the randomly selected optical depth τ scatter at particle k . Then, the optical depth is calculated from thebeginning of the sphere for particle ( k − (ensuring that all poten-tial contributors before and after this location are accounted for), it-erating over distance until the answer converges on τ scatter . As theoptical depth always increases with distance, convergence can beachieved with simple algorithms and relatively little computation.This code uses a recursive bisector algorithm to perform the itera-tion. Starting from the path length between the beginning of sphere ( k − to the end of sphere k , this value is halved recursively untilthe correct optical depth is obtained (to within some tolerance) oruntil the path length reaches a minimum value (defined as a fractionof the smallest smoothing length in the simulation). To confirm the raytracing component of the code was working cor-rectly, a simple test case was devised (Figure 6). Consider a uni-form density sphere, with radius R , density ρ and total opacity χ . A point source is located a distance D from the centre of thesphere. It emits rays at an angle θ from the vector connecting thesource and the sphere’s centre. The optical depth τ ( θ ) therefore hasan analytic solution: c (cid:13) , 1–10 Duncan Forgan and Ken Rice
Figure 5.
The four classes of particle in the system: unlisted (top left), distant-listed (top right), front-listed (bottom left) and back-listed (bottom right). The”x” denotes the emission location of the photon.
Figure 6.
Schematic of the raytracing experiment. The optical depth of theray can be calculated analytically to compare with the code’s output. τ ( θ ) = 2 ρ χ q R − D sin ( θ ) (27)Three SPH snapshots were generated, each containing a uniformdensity sphere (mass M ⊙ , R = 2133 au). Three were generatedto check convergence: snapshot 1 used SPH particles; snap-shot 2 used × ; snapshot 3 used . A point source wasplaced at a distance D = 4000 au, and the optical depth alongthe ray (assuming an opacity of unity) was calculated as a func-tion of θ . The numerical results are compared with the analyticalresult (scaled such that the maximum optical depth is 1) in Figure7. The column density is subject to the underlying random noise(at a level of around 5%) associated with generating an SPH snap-shot (which has not undergone any settling). However, the resultsvary little with increasing particle number, showing that the columndensity has converged even for the relatively low particle numberof . Figure 7.
Results of the raytracing experiment. The black solid line indi-cates the analytical result; the coloured dashed lines show simulation resultsfor increasing particle number.
Greaves et al. (2008) imaged HL Tau using the Very Large Array(VLA) with a resolution of 0.08” at a wavelength of 1.3cm (corre-sponding to a spatial resolution of 10 au at HL Tau’s distance of ∼
140 pc). They discovered excess emission at ∼ au, which theyidentified as a candidate protoplanet in the earliest stages of forma-tion, a possible example of protoplanetary disc fragmentation form-ing a bound object (Boss 1997). To lend weight to their hypothesis,they conducted SPH simulations (containing 250,000 particles) ofan unstable star-disc system with similar parameters to HL Tau, inwhich a clump forms at ∼ au, with a similar mass as that de-duced for the candidate (see Figure 8). The simulation uses an SPHcode based on the work of Bate et al. (1995). It employs individual c (cid:13) , 1–10 ative Synthetic Imaging of SPH density fields using MCRT particle timesteps, and models hydrodynamics and gravity, with ra-diative cooling as prescribed by Stamatellos et al. (2007). What cana telescope like ALMA be expected to see in HL Tau? Will the can-didate protoplanet be resolvable at millimetre wavelengths?The SPH simulation from Greaves et al. (2008) was there-fore imaged, using the same dust scattering parameters as usedin Wood et al. (2002) . The wavelength and resolution parameterswere set to be representative of ALMA : the wavelength rangeis [0 . , . cm, with an effective resolution of 0.01”. The star is . M ⊙ , and the disc is . M ⊙ , with an initial surface density pro-file of Σ ∝ r − A dust to gas ratio of 0.005 (corresponding to 50%solar metallicity) was assumed throughout. Dust sublimation is pre-scribed by enforcing any SPH particle with temperature greaterthan 1600 K to have zero dust mass. Temperatures were calculatedusing the equation of state outlined in Forgan et al. (2009) (moreinformation can be found in Stamatellos et al. 2007; Boley et al.2007). The star’s emission was modelled by a blackbody sourcewith radius R ⊙ and temperature 2000K.Figure 8 shows that the clump, which reaches temperatures ofaround 1500 K at its centre is detectable by its emission, althoughsomewhat fainter than the main star-disc system. In fact, it shouldbe expected that the clump is fainter still, if dust sublimation ismore appropriately modelled as opposed to a one temperature cut-off. The m=2 spiral mode attached to the clump is also apparent;however with a flux of around 0.001 mJy/beam, it is unlikely thatALMA will be able to image it, highlighting the need for resolutionof order ∼ and sufficient sensitivity in the detectionof disc spiral arms. The peak emission is around 0.5 Jy, which iswithin a factor of two of the observed SCUBA measurements of HLTau . The remaining discrepancy is most likely due to uncertaintyin the selected disc and star parameters (as well as the dust data).To summarise, the clump is indeed detectable with ALMA, but atelescope with better sensitivity is required to detect the spiral armscorrectly. It is expected that in an average stellar cluster, at least one star-discsystem will undergo a close encounter with another star in its life-time (Clarke & Pringle 1991; Forgan & Rice 2009b). The effect ofan encounter perturbs the disc, modifying its surface density pro-file and triggering enhanced accretion. The effect of the enhancedaccretion produces a stellar outburst. This feature was linked to theFU Orionis phenomenon (Bonnell & Bastien 1992; Pfalzner 2008),but it has been shown that such encounters are too infrequent to re-produce the correct observational signatures of FU Ori outbursts(Forgan & Rice 2009b).Forgan & Rice (2009a) carried out a series of SPH simula-tions to study the effects of a close encounter on a protostellardisc’s dynamics and structure. As a second application of the imag-ing techniques discussed here, the reference simulation describedin Forgan & Rice (2009a) was imaged at two points in its evolu-tion: the initial pre-encounter phase, where the disc remains in amarginally stable, self-gravitating state (Lodato & Rice 2004); andduring periastron, where the disc is heated strongly by the sec-ondary’s motion through it.The primary is . M ⊙ (modelled by a blackbody source with Note that the examples shown here approximate χ ν = κ ν radius R ⊙ and temperature 2000K), and the secondary is . M ⊙ (modelled by a blackbody source with radius . R ⊙ and temper-ature 1000K). The disc is . M ⊙ , with an initial surface densityprofile of Σ ∝ r − (and dust to gas ratio of 0.01, i.e metallicityequal to solar). The imaging parameters are the same as those ofthe HL Tau example above ( [0 . , . cm, with an effective reso-lution of 0.01”).Disc asymmetry plays an important role in the imaging of thissystem (Figures 9 and 10). In Figure 9, the disc displays a non-zeroellipticity due to its spiral structure, with shadowing indicating theincreased surface density of the disc corresponding to the spiralsthemselves. Again, however, these shadows are only apparent if thesensitivity of the telescope is 1 µ Jy/beam, well beyond the ALMA’scurrent continuum sensitivity estimates. This ellipticity is enhancedduring the encounter (Figure 10), with the semi major axis of theellipse aligned with the orbital vector of the primary and secondary.The erasure of the strong spiral structure during the encounter re-sults in a detectable tidal arm with flux around 1 mJy/beam. Thestellar emission is dwarfed by the disc at these wavelengths. Theinner regions of the disc are not hot enough for dust sublimation toopen a resolvable gap (the intrinsic inner gap is very much belowthe resolution limit). Again, the resolution of spiral structure in aself-gravitating disc will not be possible with ALMA unless highersensitivities are achieved, but it will be possible to detect enhancedemission from a tidal spiral wave generated as the result of a stellarencounter.
As the SPH systems being imaged by this code will be in gen-eral disordered and impossible to render analytically, a true run-time scaling is not feasible. However, an example scaling assumingsimple geometry can be calculated.In general, the runtime T goes as T ∼ N γ N steps (28)Where N γ is the number of photons emitted by the code, and N steps represents the computational expense required to track onephoton from emission to capture. This can hence be written T ∼ N γ < N ray > N scatt (29)Where < N ray > is the mean number of particles intersectedby any ray, and N scatt indicates how many times a photon willscatter before it exits. Generally, N scatt ∼ < τ > , and < N ray > ∼ N p P intersect (30)Where N p is the number of SPH particles, and P intersect isthe probability that any one particle is intersected by the ray. Thiscan be estimated for simple geometries: assuming the SPH particledistribution is a sphere, then P intersect ∼ − e − Rnσ s (31)Where R is the sphere’s radius, n = N p /V ∼ N p /R is thenumber density of the sphere, and σ s ∼ π < h > is the averagecross-section of the SPH particles. In a homogeneous sphere, < h > ∼ (cid:18) ρ (cid:19) ∼ RN p (32)Combining these results gives c (cid:13) , 1–10 Duncan Forgan and Ken Rice
Figure 8.
Left: Surface density plot of the HL Tau simulation to be imaged. Right: image of the HL Tau simulation, integrated over the wavelength range [0 . , . cm. The green cross indicates the pixel with maximum intensity. Figure 9.
Left: Surface density plot of the disc pre-outburst. The secondary is out of frame. Right: image of the simulation, integrated over the wavelengthrange [0 . , . cm. The green cross indicates the pixel with maximum intensity Figure 10.
Left: Surface density plot of the disc at periastron. Right: image of the simulation, integrated over the wavelength range [0 . , . cm. The greencross indicates the pixel with maximum intensity c (cid:13) , 1–10 ative Synthetic Imaging of SPH density fields using MCRT T ∼ N γ < τ > N p (cid:18) − e − N p (cid:19) (33)This runtime scaling shares common features with grid-basedMCRT codes, in particular the dependence on optical depth andnumber of particles/cells. For illustration, the simulations used inthis work typically took approximately 100 CPU hours to createan image, with N p = 500 , , N gamma = 10 in systems withsignificant optically thick components. SPH simulations are typically subject to several resolution con-ditions, principal among them the Jeans Resolution criterion ofBate & Burkert (1997) to correctly resolve fragmentation. Imag-ing these simulations is no different - it must be demonstrable thatthe local density field is sufficiently resolved to allow a satisfactorycalculation of the optical depth (more specifically the optical depthto the photosphere, as this is where most of the received emissionwill originate). In terms of the smoothing length, h , and the meanfree path λ mfp , this condition is h < λ mfp = 1 ρκ . (34)As the value of h is inversely proportional to the local number den-sity of SPH particles, the resolution of an MCRT image of an SPHdensity field is closely linked to the total number of particles inthe simulation, an intuitive result. This condition is satisfied in theouter disc easily, as the mean free path is of order the system scaleheight, (which is resolved by several smoothing lengths). However,within the inner 5-10 au of the disc, two separate issues arise:(i) The smoothing length exceeds the local mean free path. Thisis partially compensated by SPH’s smoothing prescription - thedensity field is defined continously across all space, and thereforeaffords a kind of sub- h resolution. However, this smoothing erasesinformation about fluctuations in the density field below this lengthscale, and this will lead to an overestimation of the escaping flux.(ii) The disc is not well resolved several scale heights above themidplane, where the photosphere is located. This under-resolvingwill result in a decrease in scattering events at the photosphere,affecting the scattered light component of the flux.These two issues show that inside ∼ au, the flux es-caping from the disc will be subject to resolution-based error.Users concerned with the inner regions of these discs couldadopt a particle splitting prescription to boost the resolution (e.g.Kitsionas & Whitworth 2002). However, as this work is concernedmainly with millimetre emission from the cooler, well-resolvedouter disc, this emission will be adequately resolved, and the is-sues described above will not play a significant role. The work described here has directly utilised the temperature struc-tures produced in the SPH simulations, but such outputs are not arequirement for imaging them using MCRT. The temperature foreach fluid element can be calculated using MCRT with so-calledradiative equilibrium methods (e.g. Bjorkman & Wood 2001). Inessence, these methods involve absorbed photons incrementing thelocal temperature field at the absorption point T ( r j ) by an amount ∆ T ( r j ) , and then being re-emitted. As many photons are passed through the system and are absorbed and re-emitted at various lo-cations, the temperature structure will relax to the correct value. Ingrid-based methods, the procedure of increasing the local tempera-ture is simple - the increment ∆ T is added to the grid cell the pho-ton is absorbed in. In SPH fields, the situation is somewhat morecomplex. A scalar field X at a location r j is defined as: X ( r j ) = X i X i m i ρ i W ( r i − r j , h i ) (35)(Where we have used the “scatter” interpretation). We can substi-tute for the temperature increase ∆ T : ∆ T ( r j ) = X i ∆ T ij m i ρ i W ( r i − r j , h i ) (36)The left hand side of equation (36) can be calculated using the tra-ditional radiative equilibrium methods, and is known. The desiredquantities are the ∆ T ij , that is the individual increases in temper-ature each SPH particle must be assigned. Equation (36) thereforeacts as a constraint on the ∆ T ij , but is insufficient to close the sys-tem. A method of closure is presented here as proof that the radia-tive equilibrium framework can indeed be incorporated - its appli-cation is left for future work. Closure can be achieved by specifyingthe following ansatz : ∆ T ij = A ij ∆ T g ( r ij ) (37)Where A ij is some real number, r ij is the separation of particle i from the location j and g ( r ij ) is some function that satisfies: g (0) = 1 , (38) g ( r ij > h i ) = 0 (39)i.e. particles with smoothing volumes that do not intersect the loca-tion j do not contribute. This will satisfy equation (36) provided X i m i ρ i A ij g ( r ij ) W ( r ij , h i ) = 1 (40)This has recast the problem from solving for ∆ T ij to the prefac-tors A ij and function g ( r ij ) . This may appear to be a sidewaysstep: however, giving a specific example is illuminating. To obtaina suitable function g ( r ij ) , consider the decay of flux F along apath of optical depth τ : F = F e − τ (41)Using the Stefan-Boltzmann law, F = σT yields σ (∆ T ij ) = σ (∆ T ) e − τ ij (42)Inspecting the above equation suggests the following form for g : g ( r ij ) = e − τ ij / (43)Assuming that A ij = A for all ( i, j ), the normalisation conditionbecomes A = 1 P i m i ρ i e − τ ij / W ( r ij , h i ) (44)The system has now been closed with an ansatz that is physicallymotivated, and reproduces expected behaviour (particles with a re-duced optical depth to the location will receive a greater tempera-ture boost than others in more optically thick regions).Other radiative equilibrium methods can also be recast in SPHform, such as that of Lucy (1999), which computes grid cell tem-peratures by counting the path lengths of any photons that passthrough it during one iteration, retrieving the temperature structure c (cid:13) , 1–10 Duncan Forgan and Ken Rice rapidly within a few iterations. By replacing grid cells with parti-cle smoothing volumes, it can be seen that this method is indeedamenable to SPH fields also.
This paper has outlined a means for applying Monte Carlo Ra-diative Transfer (MCRT) techniques directly to Smoothed ParticleHydrodynamics (SPH) density fields of star and planet formation.In doing so, it gives a means for synthetic telescope images to bemade, allowing the flexibility of SPH simulations to be retained incalculating optical depths, scattering and polarisation. This workhas shown two applications of the method to SPH simulations ofprotostellar discs: non-trivial features in the imaging of these envi-ronments are well-traced, and the resulting images give new oppor-tunities for theoretical input to observed discs. It can be shown whatresolutions and sensitivities are required to observe these features.ALMA appears to be able to resolve and image perturbed structureslike tidal arms in discs undergoing close encounters. But, while ithas sufficient resolution, it cannot image unperturbed spiral struc-ture (or the shadows it casts) without an improvement in sensitivityby at least a factor of 100. Its inherently graphical nature allowsthe code to be greatly optimised, whether by standard parallelisa-tion methods or by the use of Graphical Processing Units (GPUs)to handle the ray intersections.While used for imaging in this work, the method is not re-stricted to this alone. As the algorithm only modifies how the opti-cal depth is calculated (in order to work in SPH fields), it can be ap-plied to any circumstances traditional gridded MCRT has been ap-plied to. This includes radiative equilibrium simulations where thetemperature structure is calculated from stellar emission, or movingbeyond continuum emission to perform line radiative transfer cal-culations. In effect this places SPH on a par with grids or meshesfor MCRT techniques.
ACKNOWLEDGEMENTS
The authors would like to thank Matthew Bate and Barbara Whit-ney for useful comments which greatly strengthened this work.Plots of the SPH simulations were made using
SPLASH (Price2007). All simulations presented in this work were carried out us-ing high performance computing funded by the Scottish Universi-ties Physics Alliance (SUPA).
REFERENCES
Altay G., Croft R. A. C., Pelupessy I., 2008, MNRAS, 386, 1931Bastien P., Cha S., Viau S., 2004, in G. Garcia-Segura, G. Tenorio-Tagle, J. Franco, & H. W. Yorke ed., Revista Mexicana deAstronomia y Astrofisica Conference Series Vol. 22 of RevistaMexicana de Astronomia y Astrofisica Conference Series, SPHwith radiative transfer: method and applications. pp 144–147Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362Bate M. R., Burkert A., 1997, MNRAS, 288, 1060Bisbas T. G., W¨unsch R., Whitworth A. P., Hubber D. A., 2009,A&A, 497, 649Bjorkman J. E., Wood K., 2001, ApJ, 554, 615Boley A. C., Hartquist T. W., Durisen R. H., Michael S., 2007,ApJL, 656, L89Bonnell I., Bastien P., 1992, ApJL, 401, L31 Boss A. P., 1997, Science, 276, 1836Clarke C. J., Pringle J. E., 1991, MNRAS, 249, 584Eisemann M., Grosch T., Mller S., Magnor M., 2007, J. GraphicsTools, 12, 35Forgan D., Rice K., 2009a, MNRAS, 400, 2022Forgan D., Rice K., 2009b, MNRAS, in pressForgan D., Rice K., Stamatellos D., Whitworth A., 2009, MN-RAS, 394, 882Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375Greaves J. S., Richards A. M. S., Rice W. K. M., Muxlow T. W. B.,2008, MNRAS, 391, L74Hernquist L., Katz N., 1989, ApJS, 70, 419Hosking J. G., Whitworth A. P., 2004, MNRAS, 347, 994Kessel-Deynet O., Burkert A., 2000, MNRAS, 315, 713Kitsionas S., Whitworth A. P., 2002, MNRAS, 330, 129Kurosawa R., Harries T. J., Bate M. R., Symington N. H., 2004,MNRAS, 351, 1134Lodato G., Rice W. K. M., 2004, MNRAS, 351, 630Lucy L. B., 1977, AJ, 82, 1013Lucy L. B., 1999, A&A, 344, 282Mahovsky J., Wyvill B., 2004, J. Graphics Tools, 9, 35Mayer L., Lufkin G., Quinn T., Wadsley J., 2007, ApJL, 661, L77Monaghan J. J., 1992, ARA&A, 30, 543Monaghan J. J., 2005, Rep. Prog. Phys, 68, 1703Oxley S., Woolfson M. M., 2003, MNRAS, 343, 900Pawlik A. H., Schaye J., 2008, MNRAS, 389, 651Pfalzner S., 2008, A&A, 492, 735Pinte C., M´enard F., Duchˆene G., Bastien P., 2006, A&A, 459,797Price D. J., 2007, Publications of the Astronomical Society ofAustralia, 24, 159Price D. J., Bate M. R., 2007, MNRAS, 377, 77Price D. J., Monaghan J. J., 2004, MNRAS, 348, 139Price D. J., Monaghan J. J., 2005, MNRAS, 364, 384Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003,MNRAS, 339, 1025Stamatellos D., Whitworth A. P., 2005, A&A, 439, 153Stamatellos D., Whitworth A. P., Bisbas T., Goodwin S., 2007,A&A, 475, 37Viau S., Bastien P., Cha S., 2006, ApJ, 639, 559Whitehouse S. C., Bate M. R., 2004, MNRAS, 353, 1078Whitehouse S. C., Bate M. R., 2006, MNRAS, 367, 32Whitney B. A., Hartmann L., 1992, ApJ, 395, 529Wood K., Bjorkman J. E., Whitney B. A., Code A. D., 1996, ApJ,461, 828Wood K., Wolff M. J., Bjorkman J. E., Whitney B., 2002, ApJ,564, 887 c (cid:13)000