# TRAPHIC - Radiative Transfer for Smoothed Particle Hydrodynamics Simulations

Mon. Not. R. Astron. Soc. 000, 1–31 (0000)

Printed 23 June 2008

A (MN L TEX style ?le v2.2)

TRAPHIC - Radiative Transfer for Smoothed Particle Hydrodynamics Simulations

arXiv:0802.1715v2 [astro-ph] 23 Jun 2008

Andreas H. Pawlik1? and Joop Schaye1? 1

Leiden Observatory, Leiden University, P.O. Box 9513, 2300RA Leiden, The Netherlands

ABSTRACT

We present traphic, a novel radiative transfer scheme for Smoothed Particle Hydrodynamics (SPH) simulations. traphic is designed for use in simulations exhibiting a wide dynamic range in physical length scales and containing a large number of light sources. It is adaptive both in space and in angle and can be employed for application on distributed memory machines. The commonly encountered computationally expensive scaling with the number of light sources in the simulation is avoided by introducing a source merging procedure. The (time-dependent) radiative transfer equation is solved by tracing individual photon packets in an explicitly photon-conserving manner directly on the unstructured grid traced out by the set of SPH particles. To accomplish directed transport of radiation despite the irregular spatial distribution of the SPH particles, photons are guided inside cones. We present and test a parallel numerical implementation of traphic in the SPH code gadget-2, speci?ed for the transport of mono-chromatic hydrogen-ionizing radiation. The results of the tests are in excellent agreement with both analytic solutions and results obtained with other state-of-the-art radiative transfer codes. Key words: methods: numerical – radiative transfer – hydrodynamics – cosmology: large-scale structure of the Universe – HII regions – di?use radiation

1

INTRODUCTION

Radiation is one of the fundamental constituents of our Universe. Its interaction with baryons may lead to an energy exchange that can both heat and cool the matter, initiating pressure forces that may strongly in?uence the subsequent hydrodynamical evolution. Radiation may also exert a direct pressure force upon the matter through the exchange of momentum. Radiative interactions are furthermore often the dominating process in governing the excitation and ionization state of atoms and molecules. The inclusion of the transport of radiation into hydrodynamical simulations may therefore provide the key for interpreting the outcome of physical experiments and observational campaigns. To perform hydrodynamical simulations, the Lagrangian technique Smoothed Particle Hydrodynamics (SPH; Gingold & Monaghan 1977; Lucy 1977) is often employed. In SPH, the continuum ?uid is discretized using a ?nite set of point particles, each carrying its own collection of variables. The deformation of the ?uid by internal and external processes manifests itself in a steady redistribution of the point particles in space. All it takes to determine the values of physical ?eld variables at a given point in the sim?

E-mail: pawlik@strw.leidenuniv.nl ? E-mail: schaye@strw.leidenuniv.nl

ulation box is to perform a weighted average over the values the relevant variables take on the particles in its surrounding. This elegant simplicity of the SPH technique is one of the many reasons for its success. Although the ?rst radiative transfer calculation was already included in SPH at the very birth of this numerical technique more than thirty years ago (Lucy 1977), the detailed treatment of radiation transport in SPH is still an enormous challenge. One of the main reasons for this is certainly the high dimensionality of the problem. In fact, the radiative transfer equation depends on no less than seven variables (three space coordinates, two angles, frequency and time). Moreover, existing numerical schemes to solve the radiative transfer equation have most often only been formulated for use with uniform grids. Despite these di?culties, there has been encouraging progress over the last years, with several interesting and rather di?erent approaches (see Section 3). In this paper we present a new method to solve the radiative transfer equation in SPH. We speci?cally designed our method for use in hydrodynamical simulations exhibiting a wide dynamic range in physical length scales and containing a large number of light sources. A prominent example for such simulations are cosmological simulations of large-scale structure formation. Performing cosmological simulations is an exceedingly demanding computational task. Di?culties in describing the growth

2

Andreas H. Pawlik & Joop Schaye

photon packets are traced in cones, we refer to our scheme as traphic - TRAnsport of PHotons In Cones. The introduction of cones is required in order to perform the transport of photon packets directly on the unstructured grid de?ned by the SPH particles. Although we have designed our radiative transfer scheme to be readily coupled to the hydrodynamic evolution of the matter in SPH simulations, here we limit ourselves to the description and testing of the radiative transfer scheme itself. We will report on fully coupled radiation-hydrodynamics SPH simulations employing traphic in future work. The outline for the rest of this paper is as follows. We start with a brief review of the principles of SPH that we consider crucial for the understanding of the present work. In Section 3 we then recall the radiative transfer equation and place our radiative transfer method in context by reviewing the approaches that have been used to solve the radiative transfer problem in SPH so far. In Section 4 we give a detailed description of the ideas behind our method. Throughout we will emphasize how we satisfy the requirements set by our primary aim of performing radiative transfer in simulations exhibiting a wide range in length scales and a large number of light sources. Thereafter, in Section 5, we apply our method to the transport of hydrogen-ionizing radiation, describe its numerical implementation in the parallel stateof-the-art SPH code gadget-2 and report its behaviour in test problems. Finally, we summarise our approach and conclude with an outlook.

of the initially tiny matter density perturbations produced before and during the event of recombination and their metamorphosis into the rich structure observable today do not only arise because of our ignorance of how to properly model the governing physical processes. Often, it is simply the lack of computational power which prevents us from faithfully representing the basic actions involved: Cosmological simulations are both time consuming and memory exhaustive. To overcome these computational challenges, one can make use of advanced techniques and resort to clever approximations, reducing the computational e?ort. Consider the wide range of scales encountered in cosmological hydrodynamical simulations. According to the hierarchical model of structure formation, the ?rst structures and building blocks of the evolving universe are expected to form at small scales. The non-linear evolution at these scales shapes the distribution of the matter at all times, thus necessitating simulations of high spatial resolution. On the other hand, we also require su?ciently large simulation boxes in order to properly account for the modulation of the small-scale nonlinear evolution by the large-scale structure formation processes (e.g. Barkana & Loeb 2004) and not to be deceived by the cosmic scatter. To accommodate these two antithetic demands while keeping the number of particles representing the matter low enough to be computationally manageable, spatially adaptive SPH simulations have been invoked. It is then immediately clear that when solving the transport of radiation along with the gravito-hydrodynamics of the matter, one requires the radiative transfer scheme to be adaptive, too. Even spatially adaptive cosmological SPH simulations, however, make use of hundreds of millions of particles and are therefore still extremely memory-consuming. It is then indispensable to distribute the computational load over a large number of machines. For this reason, we require a radiative transfer scheme that is parallel on distributed memory machines. When performing radiative transfer simulations, sources of light are assigned a special importance. As a result, the computation time of most of the available radiative transfer schemes scales linearly with the number of sources in the simulation box. However, in cosmological simulations, even at times as early as 1 billion years after the Big Bang, i.e. at a redshift of z ? 6, a non-negligible amount (? 1 per cent) of baryonic matter has undergone star formation. In addition to these stellar sources of light, the intergalactic gas emits photons too, producing a radiation component often referred to as di?use radiation. Hence, without breaking the linear scaling with the number of light sources, radiative transfer simulations at the resolution of state-of-the-art cosmological hydrodynamic simulations need to be dispatched to the realm of the future. In this paper we present a radiative transfer scheme for use in SPH simulations that is adaptive, parallel on distributed memory and that avoids the linear scaling of the computation time with the number of sources, making it ideal for application in large simulations covering a wide range of length scales and containing many sources. In our scheme we follow the propagation of individual photon packets. Hence, it is explicitly photon-conserving (Abel, Norman, & Madau 1999) and can be applied to solve the time-dependent radiative transfer equation. Because the

2

SMOOTHED PARTICLE HYDRODYNAMICS

Excellent reviews of SPH exist (see e.g. Monaghan 2005; Monaghan 1992; Price 2005). Here, we just brie?y outline the basic principles we consider critical for the understanding of our radiative transfer method. At the basis of SPH lies the representation of a ?eld A(r) by its integral interpolant AI (r), Z AI (r) = d3 r ′ A(r′ )W (r ? r′ , h), (1) where the smoothing length h determines the spatial resolution. The interpolation kernel W satis?es Z d3 r ′ W (r ? r′ , h) = 1 (2)

h→0

lim W (r ? r′ , h)

=

δD (r ? r′ ),

(3)

where δD is the Dirac delta function, such that AI (r) coincides with A(r) in the limit h → 0. The last two conditions do not ?x the functional form for W . An often adopted choice is the spherically symmetric compact spline 8 ` r ?3 ` r ?2 r 1 ? < 1` 6 h? + 6 h , 0 h 2 8 ′ r 1 r 3 W (r ? r , h) = (4) < h 1 2 1? h , 2 πh3 : r 0, > 1, h where r = |r ? r′ |. From here on, when referring to the interpolation kernel W , we assume that it is of this form. The discrete representation of a ?uid by SPH particles is achieved by approximating the integral interpolant (Eq. 1) by the summation interpolant, X mj A(rj )W (r ? rj , h), (5) AI (r) ≈ AS (r) = ρj j

TRAPHIC - Radiative Transfer for SPH

where the summation is over the particles of mass m and mass density ρ. For a self-contained numerical treatment, any ?eld needs to be discretized only at the positions of the particles i representing the ?uid, X mj A(rj )W (ri ? rj , h). (6) Ai ≡ AS (ri ) = ρj j Since the kernel W is compact, only local information needs to be accessed for the evaluation of the last sum, which can be carried out in a computationally e?cient manner in parallel on distributed memory machines. For the interpretation of Eq. 6, two main approaches can be taken, depending on the choice for h (Hernquist & Katz 1989). In the scattering approach each particle j is considered as a cloud of radius hj and contributes to the ?eld at the position of particle i with the weight W (ri ? rj , hj ). In the gathering approach, each particle i searches the sphere of radius hi (in the following referred to as the sphere of in?uence, or neighbourhood) for particles (its neighbours) and obtains an estimate of the ?eld value at its position by summing the ?eld values at the positions of the neighbours j, each weighted by W (ri ? rj , hi ). The two interpretations are identical only if the spatial resolution is ?xed throughout the ?uid, i.e. if hi = hj = h. However, the SPH formalism only unfolds its true strength by allowing the resolution to vary in space, according to the local density ?eld. This is usually achieved by either directly ?xing the number of neighbours Nngb for all particles, or by requiring that the kernel volume contains a constant mass (Springel & Hernquist 2002).

3

3

RADIATIVE TRANSFER IN SPH PREVIOUS WORK

Before describing traphic, our method to solve the radiative transfer equation in SPH, we brie?y recall the radiative transfer equation and give a short overview of the main numerical methods that have been employed so far to obtain its solution in SPH simulations. The classical equation of radiative transfer reads (see e.g. Mihalas & Weibel Mihalas 1984) 1 ?Iν + n · ?Iν = ?ν ? κν ρIν . (7) c ?t In Eq. 7, Iν ≡ Iν (r, n, t) is the monochromatic intensity (units ergs cm?2 s?1 Hz?1 sr?1 ) of frequency ν at position r, n is a unit vector along the direction of light propagation and c is the speed of light. Sources and sinks of radiation are described by the emissivity ?ν (units ergs cm?3 s?1 Hz?1 sr?1 ) and the mass absorption coe?cient κν (units cm2 g?1 ), respectively. In general, both ?ν and κν are functions of r, n and time t. If we de?ne the photon number density ψν such that ψν d?dν is the number of photons per unit volume with frequencies (ν, ν + dν) travelling with velocity c into a solid angle d? around n, the intensity is given by Iν = chp νψν , where hp is the Planck constant. Eq. 7 is a partial di?erential equation depending on seven variables, which in general is hard to solve. Di?erent approximations have been employed to obtain its solution, giving rise to di?erent numerical approaches. Below we discuss the main approaches that have been taken so far to accomplish the transport of radiation in SPH simulations.

In his study of optically thick proto-stars, Lucy (1977) modelled the transport of radiation as a di?usion process by including a corresponding term in the SPH formulation of the energy equation. Brookshaw (1994) pointed out drawbacks of the particular formulation used in Lucy (1977), and re-phrased the di?usion equation in SPH such as to reduce its sensitivity to particle disorder. This new formulation of the di?usion equation was employed by Viau, Bastien, & Cha (2006) to perform collapse simulations of centrally condensed clouds. Treating the transport of radiation in the di?usion limit is, however, only a good approximation in optically thick media. In the optically thin regime, in?nite signal propagation speeds result, since the di?usion equation is a partial di?erential equation of parabolic type. This defect can be remedied by supplementing the di?usion equation with a so-called ?ux-limiter. SPH simulations with ?ux-limited diffusion have been carried out by e.g. Herant et al. (1994), Whitehouse & Bate (2006), Fryer, Rockefeller, & Warren (2006) and Mayer et al. (2007), covering a wide range of physical applications, from supernovae explosions and protostellar collapse to the study of proto-planetary disks. The (?ux-limited) di?usion approach to solve the radiative transfer equation fails in complex geometries. In particular, opaque obstacles illuminated by a point source do not cast sharp shadows, because the shadow region is ?lled in by the di?usion. To solve the radiative transfer equation without being restricted to the optically thick regime or simple geometries, Monte Carlo techniques can be employed (Oxley & Woolfson 2003, Stamatellos & Whitworth 2005, Croft & Altay 2007, Semelin, Combes, & Baek 2007). In Monte Carlo radiative transfer, individual photon packets emitted by each source are directly followed as they pass through the matter, thus simulating the physical process of radiation transport as encountered in nature. While Monte Carlo simulations allow realistic shadows to be cast behind opaque obstacles, they are computationally very demanding, since the Poisson noise inherent to the statistical description of the radiation ?eld leads to a signal-to-noise ratio that grows only with the square-root of the number of photon packets emitted. Ray-tracing schemes keep the advantages of Monte Carlo radiative transfer simulations, whilst not being affected by the statistical noise. In short, in ray-tracing schemes rays are cast from each source throughout the simulation box. The optical depth to each point in space is calculated along theses rays, and the attenuation of the ?ux emitted by the sources is obtained. Variants of ray-tracing working with photons instead of ?uxes blur the di?erences with Monte-Carlo schemes. Ray-tracing schemes are most easily implemented on a regular grid superimposed on the SPH density ?eld. However, this implies that all information about substructure on scales smaller than the size of the grid cells is ignored. As a cure, adaptive grids have been invoked (e.g. Razoumov & Sommer-Larsen 2006, Alvarez, Bromm, & Shapiro 2006). In SPH, grid-less raytracing has been introduced by Kessel-Deynet & Burkert (2000) to investigate the e?ects of ionizing UV radiation by massive stars on the surrounding interstellar medium. The Str¨mgren volume method (Dale, Ercolano, & Clarke 2007, o

4

Andreas H. Pawlik & Joop Schaye

? them only locally, between a particle and its Nngb neigh?ngb to be di?erent from Nngb , the number bours. We allow N of neighbours used during the SPH calculations. In SPH, Nngb determines the adaptive spatial resolution of the hy? drodynamical simulation. Similarly, in our scheme, Nngb sets the adaptive spatial resolution of the radiation transport. It is the ?rst of two numerical parameters in our method (see also the left-hand panel of Fig. 1). Working directly on the set of SPH particles, our scheme exploits the full spatial resolution of the hydrodynamical simulation. A further advantage of our approach is that the radiation transport is automatically parallel on distributed memory, as long as the SPH simulation itself is parallel on distributed memory. This is in contrast to radiationhydrodynamical schemes that employ individual communication schemes for the transport of radiation and the SPH computations. Our radiative transfer scheme can hence be coupled to the SPH simulation without introducing any additional computational structures. The main challenge we face when performing radiation transport by propagating photons only from a particle to its neighbours is achieving a transport that is directed: Photons travel along straight rays, while the spatial distribution of the SPH particles is generally highly irregular. In the following we give an overview of the concepts employed in our scheme and describe how we overcome this and other challenges to accomplish the transport of radiation in SPH simulations. We distinguish two types of particles present in SPH simulations that are relevant for the transport of photons: Star particles and SPH, or gas, particles. Only gas particles can interact with photons, via absorptions and scatterings. We assume, however, that both star and gas particles can be sources of photons and will therefore refer to them as source particles, or simply sources. The transport of radiation at simulation time tr over a single radiative transfer time step ?tr starts with the emission of photons by the source particles. Subsequently, these photons are propagated downstream, radially away from the source, simultaneously for all sources. In our particle-toneighbour transport scheme photons are propagated from gas particle to gas particle. On their way, photons may interact with the matter ?eld discretized by the gas particles. Each gas particle can simultaneously receive photons, possibly emitted at di?erent times, from multiple (in principle all) sources in the simulation. Their further propagation would require a loop over all propagation directions. We explicitly avoid the resulting linear scaling of the computation e?ort with the total number of sources by introducing a source merging procedure. The distribution of photons amongst the neighbours of the source particles must respect the angular dependence of the source emissivity. We achieve this by introducing for each source a set of randomly oriented, tessellating emission cones1 (Section 4.2.1), as illustrated in the middle panel of Fig. 1. The fraction of the total number of photons emitted into each cone is proportional to its solid angle. The emission

compare also Johnson, Greif, & Bromm 2007) and the ionization front tracking technique employed by Yoshida et al. (2007) are related approaches. Another grid-less radiation tracing scheme that has been applied to SPH simulations is presented in Ritzerveld & Icke (2006). Inspired by the Kessel-Deynet & Burkert (2000) approach, Susa (2006) describes a radiation-hydrodynamics scheme for the transport of ionizing radiation in cosmological simulations of structure formation. One major drawback most ray-tracing schemes share with the Monte-Carlo technique is that the computational e?ort to solve the radiative transfer equation scales linearly with the number of sources in the simulation box (but see Section 4.2.3 later on). For this reason, SPH radiationhydrodynamical simulations employing the ray-tracing approach have mostly been restricted to the study of problems involving only a few sources. In the following sections we will describe traphic, a novel method to solve the radiative transfer equation in SPH. traphic employs a photon tracing technique that works directly on the discrete set of irregularly distributed SPH particles. It is designed to overcome the challenges set by our goal of carrying out large radiation-hydrodynamics SPH simulations exhibiting a wide dynamic range in length scales and containing a large number of light sources.

4

TRAPHIC - TRANSPORT OF PHOTONS IN CONES

In this section we give a detailed description of traphic, our method to solve the radiative transfer equation in SPH. We start the presentation of our radiative transfer scheme by sketching the main idea, in order to give an overview and to introduce the reader to the underlying concepts, which will be illustrated and explained in more detail in the subsequent sections. traphic can be applied to solve the radiative transfer equation in both two-dimensional and three-dimensional SPH simulations. When illustrating concepts with schematic ?gures we will, however, restrict ourselves to two dimensions for the sake of clarity. Although we ultimately aim at performing simulations in which the radiation transport is fully coupled to the hydrodynamics, here we assume that the SPH density ?eld is static and that the radiation does not exert any thermodynamic or hydrodynamic feedback on the matter. We will report on the coupling of radiation and hydrodynamics in future work. 4.1 Overview

We solve the radiative transfer equation (Eq. 7) in ?nite steps (tr , tr + ?tr ), where tr is the current simulation time, by tracing photon packets radially away from their location of emission until a certain stopping criterion is ful?lled. We do not introduce a grid on which the photons are propagated. Instead, we propagate the photons directly on the discrete set of particles in the simulation. For the transport of photons we employ the same particle-to-neighbour communication scheme that is already used in the SPH simulation to solve the equations of hydrodynamics. That is, we accomplish the global transport of photons by propagating

1 We use the word cone in a general sense. It does not necessarily describe, but includes in its meaning, a regular cone, which is de?ned as a pyramid with a circular cross section.

TRAPHIC - Radiative Transfer for SPH

direction (which we will equivalently refer to as the propagation direction) associated with each emitted photon packet is determined by the central axis of the emission cone. The number of cones used in the emission cone tessellation, Nc , is the second numerical parameter in our method. It sets the angular resolution. The random orientation given to the cones ensures that photons are not restricted to propagate along a ?xed number of directions and prevents artefacts introduced by the shape of the emission cones. Some of the emission cones may not contain any neighbours, as is the case for the bottom-right cone of the tessellation shown in the middle panel of Fig. 1. To nevertheless emit photons into the corresponding directions, we create additional neighbours to which the photons can then be transferred (Fig. 1, right-hand panel). We refer to these neighbours as virtual particles (ViPs), since they do not take part in the SPH simulation. The properties of the ViPs are obtained from those of the neighbouring SPH particles using SPH interpolation. ViPs compensate for the lack of neighbours in certain solid angles around the emitting source. Hence, by employing ViPs we introduce the freedom of choosing emission directions independently of the geometry of the SPH particle distribution. As we will argue in Section 4.2.1 and investigate in more detail in Appendix A, this freedom is a necessary requirement for achieving a desired angular dependence of the transport process, e.g. the isotropic emission of photons by star particles, in any particle-to-neighbour transport scheme. The photon packets distributed amongst the neighbours of the sources are further transferred downstream until they experience an interaction. In contrast to the emission of photons by source particles, gas particles distribute photons only amongst the subset of neighbours located in the downstream direction (see Fig. 2). We distinguish this directed particle-to-neighbour transport from the emission process by referring to it as transmission (Section 4.2.2). In more detail, each gas particle transmits photons downstream by distributing photon packets amongst the subset of neighbours located in certain regular cones only. These so-called transmission cones are centred on the emission (propagation) direction associated with each photon packet and have an opening angle determined by the angular resolution, i.e. by Nc . They con?ne the photon packets to the cone into which they were emitted by the corresponding source particle. As a result, the photon packets are propagated radially away from the sources, in a manner that is independent of the geometry of the SPH particle distribution. Like emission cones, transmission cones might not contain any neighbours. Again we employ ViPs to accomplish the photon transport into the corresponding directions. By keeping the opening angle of the cones ?xed, the solid angle subtended by a transmission cone as viewed from the original source decreases with the distance from the source. As a result, the photon transport becomes adaptive in angle, as in the ray-splitting technique employed in ray-tracing schemes (e.g. Abel & Wandelt 2002). We will demonstrate later on, when testing a numerical implementation of our scheme (Section 5.3.2), that the use of implicitly adaptive transmission cones con?ning the photon propagation yields sharp shadows behind opaque obstacles. Gas particles receive photon packets from other particles through the processes of emission and transmission

5

described above. A given gas particle can simultaneously receive multiple photon packets, which may each be associated with a di?erent emission direction. In fact, our assumption that all star and gas particles in the simulation box can be source particles would imply that the number of directions from which photon packets can be simultaneously received by a given gas particle would scale linearly with the total number of particles in the simulation. The computational effort to accomplish the subsequent photon transmission along the associated directions would then also scale linearly with the total number of particles. Since this consideration applies to the transmission of photons by any SPH particle in the simulation, the total computational e?ort to solve the radiative transfer equation would scale no less than quadratically with the total number of SPH particles. To avoid this computationally expensive scaling, we introduce a source merging procedure (cp. Fig. 4, which will be explained in more detail in Section 4.2.3). We make use of the fact that source particles that are seen as close (in angle) to each other from a given gas particle can be approximated by, or merged into, a single point source. We average the emission directions associated with the photon packets received from sources in solid angle bins de?ned by a so-called reception cone tessellation which, except for a random rotation, is identical to the cone tessellation employed for the emission process. Consequently, photon packets need to be transmitted into at most Nc directions. We distinguish two types of interactions with the matter ?eld sampled by gas particles and ViPs that photons may experience. Absorption interactions change the number of photons contained in each packet. Scattering interactions change the propagation direction (and possibly the frequency) of the photons in a packet. All interactions are taken into account in an explicitly photon-conserving manner. We give a detailed description of interactions in Section 4.3. 4.2 Transport of photons

We now expand on the description of our radiative transfer scheme given above. We comment in detail on the main concepts underlying traphic, i.e. the emission of photons by source particles, the transmission of photons by gas particles and the source merging procedure. The description of how these concepts are employed to advance the solution of the radiative transfer equation in time is deferred to Section 4.4. 4.2.1 Photon emission by source particles

In this section we describe the emission of photons by source particles. Star particles emit photons according to their intrinsic luminosity, which can for instance be obtained from population synthesis models. An example of the emission of photons by gas particles is the emission of recombination radiation by a recombining ion. Let us consider the emission of photons by source particle i, located at ri . We denote the number of photons emitted per unit time per unit frequency ν per unit solid angle ˙ ? around the unit direction vector n by Nν,i (n, r, t), or sim˙ ˙ ply Nν,n,i . With this notation, the number of photons Nν,i R ˙ ˙ emitted per unit time at frequency ν is Nν,i ≡ d? Nν,n,i , and the total number of photons emitted per unit time is R R ˙ ˙ Nγ,i ≡ dν d? Nν,n,i .

6

Andreas H. Pawlik & Joop Schaye

Figure 1. Emission of photon packets by source particles. Left-hand panel: A source particle (black star) and its neighbourhood (big grey ? ? disc). In this example the radius h of the neighbourhood is de?ned such that there are Nngb = 4 neighbouring gas particles (white discs). Note that SPH simulations typically use Nngb ≈ 48. Middle panel: The randomly oriented emission cone tessellation of the source particle (solid lines) is shown for the case Nc = 4. The dashed arrows indicate the central axes of the cones. Note that the bottom-right cone does not contain a neighbouring gas particle. Right-hand panel: The source particle has transferred its photon packets to its neighbours (black discs). The emission direction (small solid arrows) associated with each packet points in the direction of the axis of the cone in which the neighbour resides. For the empty cone, a virtual particle is created (red hatched disc), placed randomly along the central axis of the cone.

R t+?tr ˙ Source particle i emits t Nγ,i photons over the radiative transfer time step ?tr . In agreement with our particle-to-neighbour based transport approach, the emission process is modelled by distributing these photons ? amongst the Nngb nearest gas particles, residing in a sphere ? of radius hi centred on ri . This sphere is schematically depicted in the left-hand panel of Fig. 1. Source particle i distributes its photons amongst its neighbours with the help of a set of space-?lling emission cones, de?ned as follows (middle panel of Fig. 1). We tessellate the simulation box into Nc cones of (generally not identical) solid angles ?k , k = 1, 2, ..., Nc , with the apexes i located at the position ri of particle i. The number of neigh?k bours Nngb,i in each cone k is determined, taking values in ?k ? the range 0 Nngb,i Nngb . For what follows we consider the emission of photons for each frequency ν separately. To each neighbour j in cone k a photon packet of characteristic frequency ν is emitted. kj ?k The packet contains a fraction wi > 0 (j = 1...Nngb,i ) of ˙ ν,i to be emitted per unit time the total number of photons N at frequency ν, with ˙ d? Nν,n,i ?k wj = Ri × P ?k i , Nngb,i ˙ ν,n,i d? N wl R

l=1 i

kj wi

(8)

j where wi are weights attached to neighbour j in cone k. The ?rst factor on the right-hand side of Eq. 8 determines which fraction of the total number of photons is emitted into cone k, whereas the second factor controls the fraction of photons that is transferred to neighbour j, that is, it controls the distribution of photons amongst the particles within cone k. j Here we set wi = 1, i.e. equal weights for all neighbours in a given cone. The number of cones used in the tessellation, the parameter Nc , determines the angular resolution of the radiation transport. A cone tessellation may consist of cones

with di?erent solid angles2 . We therefore de?ne the angular resolution to be the size of the average solid angle, ? = 4π/Nc . For each emission cone k the central axis is de?ned, characterised by the unit vector nk pointing away from the source (Fig. 1, middle panel). We refer to this vector as the emission direction associated to photon packets emitted into cone k. When a photon packet is emitted to a neighbouring gas particle, the emission direction is transferred in addition to the number of photons it contains (Fig. 1, righthand panel). Since the orientation of the emission cone tessellation is randomised by applying a random rotation3 at each emission, there is no a priori limit on the values the emission directions can take. Note that while we transfer the emission direction, we do not transfer the position of the source. Photon packets are traced further downstream based on their emission direction only, as we will explain in Section 4.2.2. Each photon packet has also an associated clock t? . At emission the clock time is set to the time of the simulation, t? = tr . In Section 4.4 we will see that the clock can be employed to propagate the photon packets at the speed of light. Some cones may not contain any neighbouring gas particles. This is for instance the case for the bottom right cone in the middle panel of Fig. 1. For a ?xed number of neighbours these empty cones will occur more frequently if the angular resolution is higher (i.e. if the solid angles of the cones are smaller). On the other hand, for a ?xed angular resolution Nc , empty cones will occur more frequently if the ? number of neighbours Nngb is smaller. Spatial clustering of the neighbours also increases the probability of the occurrence of empty emission cones. In the absence of a neigh2 We note that in the numerical implementation of our radiative transfer scheme the tessellation is made of cones having equal solid angles, see Appendix B1. 3 For a numerical implementation of the random rotations we refer to Appendix B2.

TRAPHIC - Radiative Transfer for SPH

bouring gas particle photons cannot be transported along the emission direction of the empty cone. We therefore create a new neighbour, a so-called virtual particle (ViP), to which the photons are transferred. The ViP is placed along ? the cone axis at a random distance < hi from the source particle, such that the volume of the cone is uniformly sampled (see the right-hand panel of Fig. 1). The properties of ViPs (e.g. density) are determined from their neighbouring gas particles using SPH interpolation. ViPs are only employed for the transport of photons. We stress that ViPs are not used by the SPH simulation. We emphasize that the introduction of cones is essential for accomplishing the emission process in our particleto-neighbour transport scheme. To see this, consider the alternative of distributing the photons directly amongst the neighbours of gas particle i. This amounts to setting ?1 = 4π i and ?j = 0 for j > 1 in Eq. 8. Without any reference sysi tem, the emission directions associated to the photon packets would then be given by the unit vectors pointing from source i to neighbour j. In Appendix A we show that in this case the net emission direction in general would correlate with the direction towards the centre of mass of the neighbours. As a result, the emission process would depend on the clustering of the neighbours in space as set by the geometry of the SPH simulation. The use of emission cones combined with ViPs gives us the freedom to choose emission directions independently of the spatial clustering of the neighbours. This allows us to model any angular dependence of the source emissivity, within the bounds of the chosen angular resolution. In summary, source particles transfer photon packets to their neighbouring gas particles using a set of randomly oriented, tessellating cones. Each cone de?nes a direction of transport, i.e. the emission direction, associated with the packets. The random orientation applied to the cone tessellation ensures that photon packets are not transferred along a ?xed set of directions only. Virtual particles (ViPs) are placed in emission cones not containing any neighbours, to which the photon packets are then transferred to. The combination of emission cones and ViPs makes the radiation transport independent of the spatial distribution of the neighbours. 4.2.2 Photon transmission by gas particles and ViPs

7

Figure 3. The apex angle ω of the transmission cones as a function of the angular resolution Nc (see Eq. 9).

? ? set of neighbours Nt,ngb Nngb located downstream. These are the neighbours residing in a regular4 cone centred on the direction of propagation of the packet (see the middle panel of Fig. 2). We refer to this cone as transmission cone. The apex of the transmission cone is located at the position of gas particle i. The solid angle of the cone is set by the angular resolution, ?t = 4π/Nc ≡ ? , and determines the apex angle ω through the standard relation ? ? ?t 180 ω = 2 arccos 1 ? . (9) × 2π π We show this relation in Fig. 3. By de?nition, a neighbour j with position rj is interior to the transmission cone with apex at the position ri of the transmitting particle i if the inner angle between the transmission cone axis and the vector rj ? ri is less than ω/2. ? The received photon packet is split into Nt,ngb packets, each of which is transferred to one of the downstream neighbours j of particle i. The number of photons contained in j each packet is set by the weights wi , such that each neighj P k bour j receives a photon fraction wi / k wi of the parent photon packet, where the sum is over all downstream neighj ? bours k = 1, ..., Nt,ngb . Here we assume wi = 1, giving equal weight to each neighbour. If there are no downstream neighbours, i.e. if the transmission cone is empty, we simply create a neighbour as in the case of empty emission cones described in the previous section. That is, we place a virtual particle (ViP) along the emission direction of the photon packet, ? a random (in volume) radial distance < hi away from the transmitting gas particle. The properties of the ViP are determined by SPH interpolation from its neighbours, and the photon packet is transferred. Each photon packet inherits the emission direction from its parent packet (Fig. 2, right-hand panel). After all packets have been transferred to the downstream neighbours, the

4

In the last section we described how a source particle distributes photon packets amongst the gas particles in its neighbourhood. For each packet we de?ned an emission direction, independently of the spatial distribution of the neighbouring gas particles by employing a randomly oriented set of emission cones. In this section we describe how the packets are propagated through the simulation box along their emission directions, by employing directed particle-toneighbour transport, a process which we refer to as transmission. Consider a gas particle i, which receives a photon ? packet. Recall that the neighbours of particle i are the Nngb nearest gas particles located in the sphere with radius ? i , h centred on particle i (Fig. 2, left-hand panel). Analogous to the case of emission, particle i re-distributes the photons contained in the packet amongst its neighbours. In contrast to emission, photons are only distributed amongst the sub-

A regular cone is a pyramid with a circular cross-section.

8

Andreas H. Pawlik & Joop Schaye

Figure 2. Transmission of photon packets by gas and virtual particles. The particle positions shown are the same as in Fig. 1, as are the neighbourhood (dashed circle) and emission cone tessellation of the source particle (white star). The transmission of photons is illustrated for one of the neighbours that received radiation from the star in Fig. 1. Left-hand panel: The neighbourhood (big grey disc) ? of the transmitting gas particle (black disc) is de?ned. As in Fig. 1, Nngb = 4. The emission direction associated with the photon packet to be transmitted is shown as the short solid arrow. Middle panel: The transmission cone is shown, centred on the emission direction of the photon packet that is to be transmitted. In this cone one downstream neighbour is found. Right-hand panel: The photon packet is transmitted to the downstream neighbour, turning it into a transmitting particle. The cycle repeats with de?ning the neighbourhood for this particle. As a result, the photon packet is propagated downstream, radially away from the source particle.

transmission is ?nished. Subsequently, the transmission procedure can be applied again. The transmission process thus generally splits each photon packet into multiple packets, propagating with the same emission direction. With each subsequent transmission individual photon packets are con?ned to a smaller fraction of the solid angle traced out by the emission cone they were emitted into (see Fig. 2, middle panel). This is similar to the technique of ray splitting employed in ray-tracing codes. Hence, in traphic the photon transport takes place adaptively in angle. The transmission procedure speci?ed above for gas particles is also applied for ViPs. Again, the transmission cone of the ViP can either contain gas neighbours or be empty. If it is empty, the ViP creates another ViP, as in the case of transmission by gas particles. After the ViP has performed the transmission, it is deleted. ViPs are therefore temporary ? constructs. For Nc ? Nngb the total number of ViPs in the ? simulation is proportional to Nc , whereas for Nngb ? Nc the total number of ViPs approaches zero. The main purpose of the transmission cones is to con?ne the downstream propagation of photons to the solid angles into which they were emitted by the source particles, which is the main challenge for schemes using particle-toneighbour transport. Just as emission cones were introduced to make the emission of photons by source particles independent of the geometry of the SPH particle distribution, transmission cones are introduced to further propagate the photons downstream, independently of the geometry of the SPH particle distribution. As a consequence, shadows will be thrown behind opaque obstacles, and their sharpness will be in agreement with the chosen angular resolution. In fact, as we will see in Section 5.3.2, the shadows are much sharper than implied by the formal angular resolution. The cones making up the emission tessellation will generally be irregular - in three dimensions there exists no tessellation of space with regular cones. The use of regular cones for the downstream propagation of photon packets emitted into a certain irregular emission cone is therefore an approxi-

mation. In principle, each photon packet could be transmitted into a cone having the shape of the emission cone it originates from. However, tracing photons within irregular cones is computationally more demanding. Furthermore, in the next section we will see the necessity for combining multiple photon packets propagating along di?erent directions into a single photon packet that propagates along a correspondingly averaged direction. Since there is no unique way of de?ning the shape of a transmission cone associated with this average direction, we use a regular shape. 4.2.3 Photon merging by gas particles

In the previous section we described the transmission of a single photon packet arriving at a gas particle. In principle, each gas5 particle can simultaneously receive multiple photon packets, possibly emitted by di?erent sources at di?erent times, propagating along di?erent emission directions, a situation which is depicted schematically in the left-hand panel of Fig. 4. Accomplishing the transmission would then require executing a loop over all received packets for each transmitting gas particle. Clearly, for each gas particle this would imply a linear scaling of the computation e?ort with the number of packets that need to be transmitted, or in other words, with the total number of source particles present in the simulation. In cosmological simulations of structure formation, for instance, where a signi?cant number of star particles can usually be found already at high redshift, the simultaneous reception of multiple photon packets will be the rule rather than the exception. Even without taking into account the di?use radiation of the intergalactic gas resulting from radiative de-excitations, recombinations and scatterings, a linear scaling with the number of sources would quickly turn the

5

We do not consider ViPs here, since by construction each ViP can receive only a single photon packet.

TRAPHIC - Radiative Transfer for SPH

9

Figure 4. Merging of photon packets by gas particles. Left-hand panel: A gas particle (white disc) is receiving photon packets simultaneously from two transmitting gas particles (black discs). The emission direction associated with each photon packet is indicated by an arrow, whose length is given by the number of photons contained in the packet. The big grey discs indicate the neighbourhoods of the transmitting gas particles, the solid lines are the transmission cones (we assume Nc = 4). For clarity, we do not show particles and photon packets that are not directly involved in the merging event. Middle panel: A randomly oriented reception cone tessellation is de?ned for the receiving gas particle. We omitted the transmission cones and the neighbourhoods of the transmitting particles. Right-hand panel: The photon packets have been transferred to the receiving gas particle, turning it into a transmitter. Their emission vectors (grey arrows), after translation to the position of the receiver, are both found to fall into the upper-right reception cone. Note that the particles itself fall into di?erent reception cones. Their positions are irrelevant for the merging, since transmitted photon packets rather than transmitting particles are merged. The emission vectors are merged through a weighted average as described in Section 4.2.3, resulting in a single photon packet (black arrow).

radiation transport into a computational task that would be too expensive to be carried out even with the most advanced supercomputer available today. A few radiation tracing schemes have attempted to tackle this scaling problem by replacing sources that are close to each other with a single point source (e.g. Cen 2002; Razoumov & Sommer-Larsen 2006; Gnedin & Abel 2001). Here we introduce a photon packet merging procedure that strictly limits the number of packets that need to be transmitted per particle to at most Nc , without restricting the number of directions along which each individual packet can propagate. For each gas particle i receiving a photon packet we de?ne a so-called reception cone tessellation, as shown in the middle panel of Fig. 4. A reception cone tessellation is a tessellating set of Nc cones attached to a gas particle. It is identical to the set of Nc cones employed for the emission of photons by source particles, but generally has a di?erent orientation. As their name indicates, and in contrast to the emission cones, reception cones are used to collect photons. Whenever a photon packet p is transferred to gas particle i, the packet is binned into one of the reception cones by examining into which reception cone its emission direction np falls, after translation to the location of the gas particle. Photon packets whose emission directions fall in one and the same cone are merged, leaving only a single packet for transmission (Fig. 4, right-hand panel). The reception cone tessellation is given a random orientation to prevent artefacts. A merged photon packet created in cone k is assigned the new, merged emission direction nm,k . This new emission P direction is de?ned as the weighted average nm,k = P wp np / wp , where the sum is over all photon packets received in emission cone k (k = 1, 2, ..., Nc ), and the weights wp are the number of photons per packet p that are to be

transmitted into direction np to the downstream neighbours of particle i (Fig. 4, right-hand panel). Similarly, the merged photon packet is associated a merged clock t? by averagm,k ing over the individual clock times t? in reception cone k, p P P i.e. t? = wp t? / wp . p m,k As a result of the photon packet (resp. source) merging at most Nc packets have to be transmitted per gas particle, fully consistent with the angular resolution of our transport scheme. The computation e?ort required to perform radiative transfer simulations with traphic can therefore be readily controlled: In the most demanding situation, i.e. in the optically thin limit, when the whole box is ?lled with radiation, it merely scales with the product of spatial and ? angular resolution, N × Nngb ×Nc , where N is the total number of particles (SPH + stars). The merged photon packets are not associated with any existing source particle in the simulation. They should be thought of as emerging from an imaginary source, whose properties are implicitly de?ned by our merging prescription. The merged photon packets are traced further downstream, radially away from this imaginary source, into the merged emission direction, exactly as was described for non-merged emission directions in the previous section.

4.3

Photon interactions with gas particles

Photons are propagated radially away from each source until they interact with the matter represented by the SPH particles. Here we distinguish between absorptions and scattering interactions and describe how they are accounted for in traphic. We remind the reader that in this work we do not discuss the e?ect of interactions on the thermodynamical and hydrodynamical evolution of the SPH particles.

10

4.3.1

Andreas H. Pawlik & Joop Schaye

Absorption exp(?τij ) of the original number of photons. By construction, this procedure explicitly conserves photons, since 1 ? e?τij + e?τij = 1. (12)

Absorptions remove photons from the beam emitted by the source. This process is described in terms of the mass absorption coe?cient κν . From the formal solution of the radiative transfer equation it can easily be seen (e.g. Mihalas & Weibel Mihalas 1984) that the fraction of photons of frequency ν that is absorbed while travelling along the straight line connecting the positions r1 and r2 is given by 1 ? exp(?τ ), where Z r2 τ = dr κν (r)ρ(r), (10)

r1

Photons absorbed by ViPs are distributed amongst all neighbours of the ViP using (photon-conserving) SPH interpolation. This is necessary, because ViPs are a temporary construct employed for the transport of photons; permanent information is only stored at the particles in the SPH simulation. For further reference, we denote the total number of photons impinging on and the total number of photons absorbed by particle i over a single time step ?tr with ?Nin,i and ?Nabs,i , resp.

is the optical depth over the distance d = |r1 ? r2 |. In traphic, absorptions are accounted for by removing the photon fraction 1 ? exp(?τij ) from photon packets propagating between particle i and its neighbouring particles j. For the calculation of the optical depth τij between these two particles we need to numerically evaluate the integral Eq. 10. This is computationally expensive, since it involves the calculation of the density ρ (and the absorption coe?cient κν ) at a large number of points along the photon path using the SPH formulation. Similar to the approach of Kessel-Deynet & Burkert (2000), we therefore approximate the density ?eld near a gas or virtual particle i by the SPH density ρi evaluated at its position. Since photon packets are propagated between particle i, located at ri , and particle j, located at rj , along their emission direction, i.e. in the direction of the unit vector n, the distance d they cover is generally not equal to the distance between the particles. Instead, we set d = |n · (ri ? rj )|, i.e., we employ the projection of the particle distance along the emission direction. This is a valid approximation far away from the source from which the photon packets originate, but it generally fails close to the source. Therefore, we only employ the projected distance when propagating photons between a transmitting particle and its neighbours. On the other hand, when propagating photons from a source particle to its neighbours, we use the particle distance, i.e. we set d = |rj ?ri |, corresponding to radially outward propagation. The distance d is used to obtain the optical depth τij between particle i and particle j, τij = κν ρj d. This approximation neglects the existence of any substructure between particle i and its neighbouring particle j. However, this does ? not result in a loss of information if the radius hi of the ? neighbourhood of particle i is chosen such that hi hi ? (resp. Nngb Nngb ). The distance d is furthermore used to advance the clock t? of each photon packet according to the rule t? → t? + d/c. (11)

4.3.2

Scattering

By synchronising the clock t? with the simulation time tr , the clock can be employed to advance photon packets at the speed of light, as we discuss in Sec. 4.4. We note in passing that the clock t? can also be used to implement the cosmological redshifting of photons. The number of absorbed and transmitted photons is now calculated as follows. From the photon packet emitted or transmitted by particle i to particle j, a fraction 1?exp(?τij ) is absorbed by particle j. The photon packet is subsequently transmitted by particle j containing a fraction

In contrast to absorptions, scatterings change the direction and possibly the frequency of the interacting photons. A scattering event can be thought of as an absorption event followed by an immediate re-emission. Scatterings can hence be described by Eq. 7 after adapting the absorption coe?cient κν and the emissivity ?ν . The re-emission re-distributes the photons in angle (and possibly frequency). For this reason the scattered photons are sometimes referred to as diffuse. Adapting the emissivity for scatterings requires evaluating an integral over the intensity. Therefore, scatterings turn the radiative transfer equation (Eq. 7) into an integrodi?erential equation. In traphic, scatterings, that is the transport of di?use photons, are modelled by a combination of absorption and re-emission events. Photons to be scattered travelling between two particles are ?rst removed from the photon packet as described in the last section. For a true absorption event, the photon energy would be lost to the thermal bath provided by the matter. In contrast, in case of scattering, a gas particle that absorbed photons re-emits them, closely following the description for source particles in Section 4.2.1. In particular, we employ randomly oriented emission cones in combination with Eq. 8 to re-emit the absorbed photons. Since modelling the scattering process means following photons emitted earlier by a source, the clock of the photon packets that are to be re-emitted is not set to the simulation time tr , as was the case for the intrinsic emission of photons by source particles. Instead, it is determined by the clock of the photon packets that are scattered and the details of the scattering interactions (e.g. there could be a time delay between absorption and re-emission, or photons could have travelled a distance greater than the (projected) particle distance d when inferring the properties of the scattering process from a sub-resolution model). If the scattering details are taken into account by correspondingly adjusting d, Eq. 11 can be employed to ?nd the new clock time. Note that it is crucial for the modelling of scatterings to employ a radiative transfer scheme which does not scale with the number of sources, since every gas particle will soon re-emit photons.

4.4

Solving the radiative transfer equation

traphic solves the radiative transfer equation (Eq. 7) in discrete time steps ?tr , the size of which depends on the

TRAPHIC - Radiative Transfer for SPH

details of the problem under study and has to be decided on a case by case basis. Hence, we defer its discussion to the second part of this work, where we present an implementation of our method to solve speci?c problems. In this section we brie?y review the concepts described in the previous sections used to obtain the solution to the radiative transfer equation over a single time step before discussing how this solution is advanced in time. Starting at simulation time tr , the size of the radiative transfer time step ?tr is determined. Source particles emit ? photon packets to their Nngb neighbouring gas particles employing Nc cones. The total number of photons contained ˙ in the packets is determined by the number of photons Nν,i to be emitted per unit time in frequency bin ν and by the size of the time step ?tr . For each neighbouring gas particle the optical depth to the source particle is evaluated and the number of photons interacting with the matter on their way from the source particle is inferred. Some photons may be absorbed, some may be scattered; the remaining photons are left unimpeded for transmission. For the subsequent transmission of photons by the gas particles, photons associated to sources that are seen as close (in angle) to each other are merged, limiting the number of directions into which photons have to be transmitted to at most Nc . Next, for each transmitting and scattering ? gas particle the Nngb neighbours are found. Photons to be ? transmitted are distributed amongst the subset Nt,ngb of neighbours that are located downstream, as determined by the corresponding (merged) emission directions. Photons to ? be scattered are distributed amongst all Nngb neighbours following the prescription in Section 4.3. Virtual particles (ViPs), which are placed in cones that do not contain any neighbours, are deleted after having transmitted or scattered their photons. Photons that were absorbed by a ViP are distributed amongst the neighbours of the ViP using (photon-conserving) SPH interpolation. The cycle described in the previous paragraph repeats until a stopping criterion is satis?ed. The form of the stopping criterion depends on whether one solves the timeindependent or the time-dependent radiative transfer equation. When solving the time-independent radiative transfer equation, i.e. when neglecting the ?rst term on the left hand side of Eq. 7, one formally assumes c → ∞. Accordingly, the stopping criterion becomes independent of the speed of light. The cycle can then for example be repeated until all photons have been absorbed or have left the simulation box. In contrast, when solving the time-dependent radiative transfer equation, the stopping criterion is directly tied to the speed of light c: the cycle is repeated until all photons have propagated a distance c?tr . To decide whether or not a photon packet has travelled over the distance c?tr , we employ its clock t? . For further illustration it is convenient to explicitly follow the clock of a photon packet as it is propagated through the simulation box. Upon emission, the clock of the photon packet is initialized with the simulation time tr , as already mentioned in Sec. 4.2.1. After the emission, when the photon packet has been transferred from the source to one of the neighbouring SPH particles (or to a ViP), its clock is advanced according to Eq. 11. It is then checked whether the clock has reached or crossed the threshold value t? = tr + ?tr . th

11

(13)

If not, the photon packet is transmitted further downstream. As before, its clock is advanced according to Eq. 11. The photon packet is repeatedly transmitted and its clock is advanced until it has reached or crossed the threshold value de?ned in Eq. 13. At this point, the value displayed by its clock can in full generality be expressed as t? = t? + ?0 , th with ?0 0. The propagation error, ?0 , is typically larger than zero because the set of particles on which the photon packets are propagated is discrete, with the particle-to-neighbour distances generally not related to c?tr . Employing Eq. 13 to stop photon packets therefore results in the photon packets typically being propagated too far. Due to the Lagrangian nature of the SPH simulation, ?0 is individual for each photon packet. If d is the particle-to-neighbour distance corresponding to the emission or transmission event after which the photon packet was stopped, then ?0 < d/c. The propah gation error ?0 is therefore strictly limited, ?0 < maxi ? i /c, where the maximum is over all emitting and transmitting particles. Note that ?0 becomes smaller with higher spatial resolution. Once stopped, a photon packet is held (and hence its clock stands still) until the simulation time progresses to a value tr +?tr > t? , at which point the packet is transmitted and its clock is advanced again. The clocks of photon packets are thus kept in synchronization with the current simulation time, which e?ectively matches their propagation speed to become the speed of light6 . Photon packets may be held over several radiative transfer time steps, as ?0 may be (several times) larger than ?tr . While the photon packet is held, its propagation error with respect to the current simulation time steadily decreases with each passing radiative transfer time step. Since t? is ?xed while tr is increasing in steps of ?tr , the current propagation error decreases according to ?n = ?n?1 ? ?tr , where n indicates the number of radiative transfer time steps that have passed since the photon packet was stopped. Immediately before the photon packet is propagated again, its propagation error is strictly limited by the size of the time step, ??tr ?n < 0. The photon packet is then propagated such as to again bring the clock into synchronization with the simulation time, as described above. When multiple photon packets are merged into a single packet, the clock of the merged packet is determined by the averaging procedure described in Section 4.2.3. As a result, photons may be propagated over distances that di?er somewhat from the case without merging. We have seen above that, due to the synchronization of the photon packet propagation with the current simulation time tr , clocks display only values in the interval tr t? tr +?tr +?0 . The di?erence in the clocks of di?erent photon packets and hence the photon propagation error introduced by the merging procedure can therefore be controlled to become arbitrarily small by increasing the temporal (i.e. decreasing ?tr ) and spatial (leading to smaller ?0 ) resolution.

6

We note that the clocks can also be employed to propagate photon packets at speeds c di?erent from the speed of light (cp. the ? reduced speed of light approximation suggested in Gnedin & Abel 2001), by simply making the replacement c → c. ?

12

Andreas H. Pawlik & Joop Schaye

To employ the resampling, we think of particle i at position ri in the simulation box as being de-localised within its sphere of in?uence of radius hi centred on ri and assume that the probability of ?nding it in the volume d3 r around a particular point r in that sphere is given by the value of the interpolation kernel at that point (cp. the scattering approach of Section 2). For the radiative transfer on a static set of particles with positions ri we then periodically Monte Carlo re-distribute the particles within their sphere of in?uence according to this probability. We note that while we change the positions of the SPH particles, we keep all their other properties (e.g. density) ?xed in order to avoid introducing scatter in the physical variables. We describe the numerical implementation of the resampling in Appendix B3. In Sections 5.3.1 and 5.3.2 we will demonstrate that the application of this recipe leads to a signi?cant reduction of particle noise. We will, however, also see that for most applications the noise is small. We therefore do not employ the resampling technique in our simulations, unless explicitly stated. Since the resampling randomly o?sets the SPH particles from their positions provided by the hydrodynamical simulation, the apexes of the transmission cones attached to them are randomly o?set, too. The resampling procedure may therefore lead to a slight di?usion of photons out of the emission cone they were emitted into, e?ectively decreasing the angular resolution. Note that a similar di?usion of photons could occur when solving the radiative transfer equation coupled to the hydrodynamics (on which we will report in a future publication), because of the physical particle motion. We will study this di?usion in our numerical implementation of traphic in Section 5.3.2.

After the stopping criterion (time-independent or timedependent) has been ful?lled for all photon packets, the properties of gas particles are updated according to the radiative interactions that occurred. We emphasize that the properties of the gas particles (e.g. the ionization state) are therefore only updated at the end of each time step. Within each time step the order in which photon packets from different sources arrive at gas particles is therefore irrelevant. After updating the particle properties, the simulation time is advanced, tr → tr + ?tr . Finally, the size of the next time step is determined and the radiation is transported as described above. In summary, we have presented a radiative transfer scheme for use in SPH simulations that works directly on the unstructured grid formed by the discrete set of irregularly distributed SPH particles. traphic thus employs the full spatial resolution of the SPH simulation. It achieves directed transport of radiation by adaptively tracing photon packets in cones. traphic can be used to solve both the time-independent and the time-dependent radiative transfer equation in an explicitly photon-conserving way. Our scheme is by construction parallel on distributed memory machines if the SPH simulation itself is parallel on distributed memory machines. Furthermore, the computation time necessary to accomplish the radiation transport does not scale linearly with the number of sources. Instead, it merely scales with the product of spatial and angular resolution, making our scheme suitable for simulations containing a large number of sources as well as for taking into account a di?use radiation component. 4.5 Reduction of particle noise

The radiative transfer equation describes the propagation of photons within a continuum medium. In our scheme, however, photons are propagated on the discrete set of SPH particles, localised at irregular positions dictated by the SPH simulation. Consequently, numerical noise may arise from the discreteness and irregularity of the spatial distribution of SPH particles and this could in?uence the numerical solution of the radiative transfer equation obtained with traphic. To see how this particle noise can be reduced, it is helpful to employ a formal analogy (e.g. Monaghan 2005) between estimating the density in SPH simulations and estimating a probability density from sample points: We can consider the positions of the SPH particles as a random7 sample drawn from a probability density function proportional to the mass density. We can therefore Monte Carlo resample8 the density ?eld without sacri?cing its information content. Periodically assigning new positions to the SPH particles during the radiation transport by resampling the density ?eld should lead to a better representation of the continuum physics by the discrete set of SPH particles and hence lead to a reduction of particle noise.

We stress that SPH itself is not a Monte Carlo method, as ?rst noted in Gingold & Monaghan (1978). We only appeal to the Monte Carlo picture for use with the radiative transfer, not for the hydrodynamic evolution of the particles in the SPH simulation. 8 We explicitly note that the resampling does not a?ect the hydrodynamical simulation, for which the original positions are used.

7

5

APPLICATION - TRANSPORT OF IONIZING RADIATION

In this section we apply the radiative transfer scheme presented in Section 4 to the transport of ionizing radiation. Ionizing radiation is thought to play a key role in determining the ionization state and shaping the spatial distribution of the baryonic matter in the universe on both small and large scales. Examples include the triggering and quenching of star formation through radiative feedback from nearby ionizing stellar sources both in the early (e.g. Yoshida et al. 2007; Wise & Abel 2007; Johnson, Greif, & Bromm 2007; Alvarez, Bromm, & Shapiro 2006; Susa & Umemura 2006) and present-day universe (e.g. Gritschneder et al. 2007; Dale, Bonnell, & Whitworth 2007), the thin shell instability (for a recent simulation see Whalen & Norman 2007) and the growth and percolation of ionized regions generated by the ?rst stars and quasars during the so-called Epoch of Reionization (for recent simulations see e.g. Iliev et al. 2006a; Trac & Cen 2006; Kohler, Gnedin, & Hamilton 2007; Paschos et al. 2007). In the following we describe a numerical implementation of traphic, our radiative transfer scheme, speci?ed for the transport of mono-chromatic hydrogen-ionizing radiation, in the state-of-the-art SPH code gadget-2 (Springel 2005). We start in Section 5.1 with a short review of the photo-ionization process. In Section 5.2 we present the details of our numerical implementation, establishing the con-

TRAPHIC - Radiative Transfer for SPH

nection to the general description given in Section 4. Finally, in Section 5.3, we perform several radiative transfer simulations, comparing the performance of our numerical implementation of traphic in three di?erent test problems to both analytic solutions and numerical simulations carried out with other radiative transfer codes as reported in the literature. 5.1 Photo-ionization rate equation τeq ≡ τion τrec . τion + τrec

13

(20)

From Eq. 19 we see that the equilibrium ionized fraction is exponentially approached on the instantaneous ionization equilibrium time-scale τeq . We will employ this time-scale later on for the numerical integration of the rate equation.

5.2 Here we brie?y recall the principles of the photo-ionization and recombination process occurring for a hydrogen-only gas of mass density ρ exposed to hydrogen-ionizing radiation. We will employ the equations derived in this section in the description of the numerical implementation of traphic. Hydrogen-ionizing photons can be absorbed by neutral hydrogen. We approximate the frequency-dependence of the mass absorption coe?cient κν for hydrogen-ionizing radiation (e.g. Osterbrock 1989) by σν nHI κν ≡ (14) ρ ? ??3 ν Θ(ν ? ν0 ), (15) σν = σ0 ν0 with nHI = (1 ? χ)ρ/mH the neutral hydrogen number density, ν0 the Lyman-limit frequency of energy hp ν0 = 13.6 eV, σ0 = 6.3 × 10?18 cm2 the absorption cross-section for photons at the Lyman-limit, mH the mass of a hydrogen atom and Θ(x) the Heaviside step function; the ionized fraction is χ ≡ nHII /nH . The number of photo-ionizations per unit time per neutral hydrogen atom at a certain point in space is determined by the photo-ionization rate Γ, Z ∞ 4πJν (ν) Γ= dν σν , (16) hp ν 0 R where Jν ≡ d? Iν /(4π) is the mean ionizing intensity. The rate of change of the neutral fraction η ≡ 1 ? χ at this point is then d η χ ? . (17) η = α(T )ne χ ? Γη ≡ dt τrec τion In the last equation, α(T )ne is the number of recombinations occurring per unit time per ionized hydrogen atom, τrec is the recombination time-scale and τion is the photo-ionization time-scale.9 With the de?nition χ ≡ τrec /(τion +τrec ) we can rewrite ? Eq. 17 to read dχ χ?χ ? =? . dt τion χ ? (18)

Numerical implementation

We have adapted traphic for the transport of hydrogenionizing radiation according to the physics of photoionization presented in Section 5.1 and implemented it using a single frequency bin in the parallel N-body-Tree-SPH code gadget-2 (Springel 2005). The description of this implementation is the subject of this section.

5.2.1

Transport of ionizing photons

The transport of ionizing photons is performed in ?nite radiative transfer time steps of size ?tr , employing the absorption coe?cient κν given by Eqs. 14 and 15. Photon packets are transported in cones according to the description given in Section 4. At the end of each radiative transfer time step, i.e. at time tr + ?tr , where tr is the current simulation time, we know the number of ionizing photons ?Nin,i impinging on and the number of ionizing photons ?Nabs,i absorbed by particle i over the time interval ?tr (cp. Section 4.3). The photo-ionization rate Γi is obtained directly, that is without explicitly referring to the mean intensity Jν , using

t ηi r t ? ? mtr Xi r tr i Γi ?tr = ?Nin,i 1 ? exp(?τitr ) , mH

(21)

where Xi is the hydrogen mass fraction and ? ? ?Nabs,i τitr ≡ ? ln 1 ? ?Nin,i

(22)

is the a posteriori optical depth and we use superscripts to indicate the time at which quantities are evaluated. In the next section we describe how the photo-ionization rate is used to update the neutral fraction of particle i.

5.2.2

Solving the rate equation

In the simulation the neutral fraction associated with any gas particle i is assumed to be constant over the time step ?tr (see Section 4.4). The correspondingly discretized rate equation would read (cp. Eq. 17),

t t t ηi r +?tr ? ηi r = αtr ntr χtr ?tr ? Γtr ηi r ?tr i e,i i i

Setting dχ/dt = 0 yields the the equilibrium ionized fraction χeq = τrec,eq /(τion + τrec,eq ). Over time-scales that are short compared with τrec /|dτrec /dt| and ne /|dne /dt|, Eq. 18 constitutes a ?rst order linear homogeneous di?erential equation in χ? χ with constant coe?cients, whose solution reads ? χ(t) ? χeq

9

(23)

=

(χ(t0 ) ? χeq )e

? τ 0 eq

t?t

(19)

Note that in Eq. 17 collisional ionizations can easily be taken into account by replacing Γ with (Γ + C(T )ne ), where C(T )ne describes the number of collisional ionizations per unit time per neutral hydrogen atom. In this work, however, we assume that collisional ionizations are unimportant, setting C ≡ 0 throughout.

According to Eq. 19, however, the neutral fraction evolves on the time-scale τeq during the photo-ionization process. In order to accurately follow the evolution of the neutral fraction using Eq. 23 one would have to choose time steps ?tr ? τeq limited by the ionization equilibrium time-scale, over which our assumption of a constant neutral fraction is a good approximation. In our implementation we choose the radiative transfer time step ?tr independently of the time-scale τeq , which can be prohibitively small, by employing the following subcycling strategy at the end of each radiative transfer time step. We only assume that the ?ux dNin,i /dt of ionizing

14

Andreas H. Pawlik & Joop Schaye

because ionization equilibrium implies that the number of photo-ionizations dNin /dt(1 ? e?τ )?tr over the time interval ?tr exactly cancels the number of recombinations (1 ? η)2 NH nH α?tr over that same time interval. This balance, however, has a unique (and stable; see Eq. 19) solution for the neutral fraction. The equilibrium neutral fraction then also implies the correct photo-ionization rate, via Eq. 21. When one is interested in following the details of the evolution of the neutral fraction towards photo-ionization equilibrium, on the other hand, the dependence of the photoionization rate on the neutral fraction needs to be taken into account, as presented above.

photons impinging on particle i is constant over the time step ?tr , i.e. dNin,i /dt = ?Nin,i /?tr , in agreement with the temporal discretisation of the radiation transport in traphic. We then a posteriori follow the evolution of the neutral fraction in time ti ∈ (tr , tr + ?tr ) on sub-cycles δti ?tr ,

t t t ηi i +δti ? ηi i = αti nti χti δti ? Γti ηi i δti , i e,i i i

(24)

Our assumption of a constant ionizing ?ux implies that the photo-ionization rate Γ ∝ (1?e?τ )η ?1 (see Eq. 21). Hence, a change in the neutral fraction implies a change in the photoionization rate Γti , i ? – t 1 ? exp(?τiti ) ηi r (25) Γt i = Γ t r i i t , 1 ? exp(?τitr ) ηi i where Γtr and τitr are the photo-ionization rate and the i optical depth at the beginning of the sub-cycling, given by t t Eqs. 21 and 22, and τiti = τitr ηi i /ηi r . The number of ionizations ?Nsub,i occurring over ?tr P t t is then ?Nsub,i = (mtr Xi r /mH ) Γti ηi i δti , where the i i sum is over all sub-steps δti in (tr , tr + ?tr ). We set ti δti ≡ min(f τeq,i , tr + ?tr ? ti ), where f < 1 is a dimensionless factor. It can be demonstrated that for this choice of δti Eq. 24 is a su?ciently good approximation to Eq. 17 that the neutral fraction never violates the physical bound 0 ηi 1. Because the neutral fraction changes during the sub-cycling, the number of ionizations ?Nsub,i can be less than the number of photons ?Nabs,i that have been removed due to absorptions during the radiation transport over the time step ?tr based on the assumption of a constant neutral fraction. We then explicitly conserve photons by adding ?Nabs,i ? ?Nsub,i photons to the photon transport in the next radiative transfer step. When either the photo-ionization rate or recombination rate is high, τeq and hence δt will be very small (dropping the particle index i for simplicity). For δt ? ?tr the sub-cycling becomes computationally very expensive. We ?nd, however, that photo-ionization equilibrium is typically reached after only a few sub-cycles. Once photo-ionization equilibrium is reached, an explicit integration of the rate equation is no longer necessary. For a photon-conserving transport we still require knowledge about the number of photo-ionizations and recombinations occuring during the equilibrium phase. Both can, however, be obtained in a stroke, based on the number of photo-ionizations and recombinations occuring during the last sub-cycle step over which the rate equation was explicitly integrated. The importance of properly following the evolution of the photo-ionization rate in the presence of an evolving neutral fraction in the optically thick10 regime has been pointed out by Mellema et al. (2006). There, a time-averaged photoionization rate obtained from an iterative procedure was employed. The sub-cycling procedure (Eq. 25) presented here has the advantage that it is well-motivated also in the presence of recombinations. We note that whenever one is only interested in obtaining the equilibrium neutral fraction, the detailed handling of the photo-ionization rate is rather unimportant. This is

10

5.2.3

The time step ?tr

For τ ? 1 Eq. 25 implies that the photo-ionization rate is t constant, Γii = Γtr . i

Our main consideration when choosing the size of the radiative transfer time step for the simulations we are presenting in this work, is that we wish to accurately reproduce the analytical and numerical reference results we are comparing with. These results include the time-dependence of the size of ionized regions around ionizing sources. At early times, just after the sources start to emit ionizing photons, the ionized regions expand quickly into the neutral hydrogen ?eld surrounding the sources. To accurately reproduce this early phase of rapid growth, we necessarily have to employ time steps ?tr that are relatively small. The phase of rapid growth is, however, only of relatively short duration. The subsequent evolutionary stages of modest resp. slow growth, which account for most of the (simulation) time, are often more interesting. We show in Section 5.3.1 that whenever we are not interested in the very early phase of rapid growth we can, thanks to the photon-conserving nature of traphic, choose substantially larger time steps without a?ecting the outcome of our simulations. For all but one of the simulations we present in this work, we will be concerned with solving the timeindependent radiative transfer equation. In these simulations, we choose to propagate photon packets downstream from their location of emission over only a single interparticle distance per radiative transfer time step, unless stated otherwise. This approach is equivalent to solving the time-independent radiative transfer equation in the limit of small radiative transfer time steps, ?tr → 0. We have explicitly checked for all our simulations that the time step was chosen su?ciently small to be in agreement with this limit. Our treatment of the time-independent radiation transport reduces the computational e?ort for the simulation of problems for which the time-step has been ?xed to a small value, e.g. by considerations like those presented in the beginning of this section. In the limit that radiation completely ?lls the simulation box, the computational e?ort required to solve the time-independent radiative transfer equation over the time interval T by propagating photon packets only over a single inter-particle distance per radiative transfer time step ?tr is proportional to (cp. Section 4.2.3) ? NSP H × Nngb × Nc × T /?tr . This has to be compared to the computational e?ort required to follow all photon packets over each time step until they leave the simulation box, 1/3 ? which is proportional to NSP H × Nngb ×Nc ×NSP H ×T /?tr .

TRAPHIC - Radiative Transfer for SPH

This is larger by a factor of NSP H , which for typical simulations reaches values of the order of 100. In Sec. 5.3.3 we will present one simulation in which we solve the time-dependent radiative transfer equation. In this simulation we will employ the photon clock (Section 4.4) to accurately control the distance over which photon packets are propagated over each time step ?tr to match the light crossing distance c?tr . In this context it is interesting to note that in the case of small time steps, i.e. c?tr < Lbox , it is less expensive to solve the time-dependent than the timeindependent radiative transfer equation. This is because the computational e?ort to solve the time-dependent equation (again in the limit that radiation completely ?lls the simulation box and assuming that its boundaries are photon1/3 ? absorbing) scales as NSP H ×Nngb ×Nc ×c?tr /Lbox ×NSP H × T /?tr , which is smaller than the computational e?ort for obtaining the time-independent solution (assuming that over each radiative transfer time step photon packets are transported until they leave the box) by the factor11 c?tr /Lbox . 5.2.4 E?ective multi-frequency description

1/3

15

In our current implementation we use only a single frequency. Thus, we either assume that the ionizing radiation is mono-chromatic, or we assume ionizing radiation with a frequency spectrum Jν and provide an effective multi-frequency description using only a single frequency bin (for a textbook introduction to the numerical treatment of multi-frequency radiation, see e.g. Mihalas & Weibel Mihalas 1984). For the latter case, we de?ne a frequency-independent (grey) photo-ionization cross-section σ such that the total ? photo-ionization rate (Eq. 16) is conserved, Z ∞ Z ∞ 4πJν (ν) 4πJν (ν) Γ= dν σν ≡ σ ? , (26) dν hp ν hp ν 0 ν0 or, σ= ? Z

∞

settings, however, no analytic solution to the radiative transfer problem is known, mainly due to the complex geometries involved. We have therefore designed the test problems in Sections 5.3.1 and 5.3.3 to closely follow the description given in the Cosmological Radiative Transfer Codes Comparison Project (Iliev et al. 2006), which provides a set of very useful numerical reference solutions to compare with. Throughout we will assume that the density ?eld is static. To facilitate a direct comparison with Iliev et al. (2006), we present results after mapping physical quantities de?ned on the SPH particles to a regular grid, unless stated otherwise. This is done using a massconserving SPH interpolation similar to the one described in Alvarez, Bromm, & Shapiro (2006). We opted for the SPH interpolation since it is most consistent with the SPH simulation we are employing. For comparison, we repeated our analysis using TSC mass-conserving interpolation (Hockney & Eastwood 1988) but found no signi?cant di?erences. For all tests reported in this section we employ the onthe-spot approximation (e.g. Osterbrock 1989), in order to allow a direct comparison with Iliev et al. (2006). In the on-the-spot approximation di?use photons emitted during recombinations to the hydrogen ground energy level are assumed to be immediately re-absorbed by neutral hydrogen atoms close to the location of emission. The e?ect of diffuse recombination radiation can then be simply taken into account by setting the recombination coe?cient α to the socalled case B value αB , which for the temperature range of interest can be well approximated by ??0.7 ? T cm3 s?1 . (28) αB (T ) = 2.59 × 10?13 104 K We will report on the study of di?use radiation in which we will be explicitly following the scattering of di?use photons instead of employing the on-the-spot approximation in future work. In their simulations, Iliev et al. (2006) assumed that the speed of light is in?nite, i.e. they solved the timeindependent radiative transfer equation. For the comparison with these simulations we will therefore make the same assumption (recall Sec. 5.2.3 for the discussion of how we solve the time-independent radiative transfer equation). From now on, when referring to the radiative transfer equation, we therefore assume it to be of the time-independent form, unless stated otherwise. We will repeat one of the simulations presented in Sec. 5.3.3 to solve the time-dependent radiative transfer equation by employing the photon packet clocks as described in Sec. 4.4. 5.3.1 Test 1: Spherically-symmetric HII region expansion

dν

0

4πJν (ν) σν × hp ν

?Z

∞

dν

ν0

4πJν (ν) hp ν

–?1

.

(27)

5.3

Tests

In this section we report on the performance of our implementation in simple and well-de?ned test problems. The tests comprise the evolution of the ionized region around a single star in a homogeneous medium (Section 5.3.1), the casting of a shadow behind an opaque obstacle (Section 5.3.2) and the propagation of ionization fronts driven by the ionizing radiation of multiple stars in an inhomogeneous density ?eld obtained from a cosmological simulation (Section 5.3.3). We compare the results obtained with traphic with analytic solutions, where available. For most physical

11

Observe that for larger time steps, i.e. ?tr Lbox /c, and assuming photon-absorbing boundaries, the computational e?ort for solving the time-dependent radiative transfer equation equals the computational e?ort for solving the time-independent radiative transfer equation. It exceeds the computational e?ort for solving the time-independent radiative transfer equation with photon packets propagating only a single inter-particle distance per ra1/3 diative transfer time-step if c?tr > Lbox /NSP H .

In this section we consider the problem of a steady ion˙ izing source emitting Nγ mono-chromatic photons of frequency hp ν = 13.6 eV per unit time, which is turned on in a static, initially neutral, uniform medium of constant hydrogen number density nH . The analytical solution to this problem is well-known (for a textbook discussion see e.g. Dopita & Sutherland 2003). In equilibrium, the number of photons emitted by the source has to compensate the number of recombinations within the surrounding, stationary, ionized region, the

16

Andreas H. Pawlik & Joop Schaye

Figure 6. Test 1. Slice through the simulation box at z = Lbox /2 showing the neutral fraction at the end of the simulation (tr = 500 Myr), i.e. in photo-ionization equilibrium, for simulations with (bottom row) and without (top row) resampling of the density ?eld. The angular resolution increases from Nc = 8 in the left-most to Nc = 128 in the right-most column, as indicated in the panel titles. The spatial ? ? resolution is ?xed to NSP H = 643 , Nngb = 32 and is indicated by the cross of length 2 h in the upper left corner of each panel. Black contours show neutral fractions of η = 0.9, 0.5, log η = ?1, ?1.5, ?2, ?2.5, ?3, ?3.5, ?4, going from the outside in. The white dot-dashed circle indicates the Str¨mgren sphere, which should be, and is, just inside to the η = 0.5 contour. The colour scale is logarithmic and has o a lower cut-o? log η = ?5. Note that the resampling strongly suppresses the particle noise seen in the top-row panels.

Figure 7. Test 1. Spherically averaged neutral and ionized fractions within the Str¨mgren sphere at tr = 500 Myr for di?erent angular o resolutions, as indicated in the legend. The pro?les in the left-hand (right-hand) panel correspond to simulations without (with) resampling ? the density ?eld. The spatial resolution is ?xed (NSP H = 643 , Nngb = 32). The black solid curves corresponds to the exact pro?les, given by the solution of Eq. 33. The grey bands show the range of neutral and ionized fractions found by other codes as reported in Iliev et al. ? (2006). The vertical dotted line marks the radius r = 2 h , corresponding to the spatial resolution employed. In the right-hand panel, the additional (right-most) vertical dotted line indicates the radius corresponding to the spatial resolution 2 h of the SPH simulations, which is the scale on which particle positions are randomly displaced during the resampling.

TRAPHIC - Radiative Transfer for SPH

ηeq (r)nH 4πr 2 Z ˙ dν Nγ (ν)e?τν σν = χ2 (r)n2 αB , eq H

17

(33)

where the optical depth τν (r) is given by Z r τν (r) = nH σν dr ′ ηeq (r ′ ).

0

(34)

Figure 5. Test 1. The evolution of the ionization front for the angular resolutions Nc = 8, 16, 32, 64 and 128, as indicated in the ? legend. The spatial resolution is ?xed (NSP H = 643 , Nngb = 32). The top panel shows the position of the ionization front rI,num normalized by the Str¨mgren radius rs as a function of time. The o thick black solid curve shows a numerical reference solution obtained with a one-dimensional, grid-based radiative transfer code (see text). The grey curve shows the analytic reference solution, Eq. 32, which has been obtained by assuming χ ≡ 1 throughout the ionized region. The results from the numerical simulations employing traphic closely match the numerical reference solution. The bottom panel shows the position of the ionization fronts of the top panel divided by the analytic reference solution. Note that the analytic reference solution slightly di?ers from the numerical reference solution, due to the simplifying assumptions inherent to the analytic approach (see also the discussion of Eq. 33).

so-called Str¨mgren sphere. Assuming that the Str¨mgren o o sphere is fully ionized, i.e. χ ≡ 1, its radius rs is therefore given by the balance equation 4 3 ˙ Nγ = πrs αB (T )n2 . (29) H 3 The evolution of the spherical, fully ionized region towards the equilibrium Str¨mgren sphere is described by the o following equation for its radius rI , the ionization front, 4 3 drI ˙ = Nγ ? πrI αB (T )n2 . (30) H dt 3 Introducing the new variables ξ ≡ rI /rs and τ ≡ t/τs , where the Str¨mgren time-scale τs = 1/(αB nH ) is the recombinao tion time in a fully ionized gas, we arrive at the di?erential equation

2 4πrI nH

dξ 1 ? ξ3 = , dτ 3ξ 2 the solution of which reads rI (t) = rs (1 ? e?t/τs )1/3 .

(31)

(32)

Hence, the ionization front reaches the Str¨mgren sphere o after a few Str¨mgren times τs and stays static thereafter. o In reality the neutral fraction inside the ionized region is, however, not zero, but varies smoothly with the distance to the ionizing source. We therefore invoke the commonly employed de?nition of the ionization front as the radius at which the neutral fraction equals 50 per cent, i.e. η = 0.5. The equilibrium neutral fraction ηeq (r) = limt→∞ η(r) can be obtained by solving the equation (e.g. Osterbrock 1989)

We refer to the neutral fraction ηeq (r) obtained by direct numerical integration of Eq. 33 as the exact neutral fraction pro?le and to χeq (r) = 1 ? ηeq (r) as the exact ionized fraction pro?le. These pro?les are shown as black solid curves in Fig. 7. Note that for our choice of parameters, rI,eq = 1.05 rs . The ionization front obtained from the solution to Eq. 33 is thus at a slightly larger radius than the equilibrium ionization front obtained assuming the Str¨mgren o sphere is fully ionized (Eq. 32). The last observation indicates that the analytic solution given by Eq. 32 generally fails to provide an accurate reference solution that can be used to judge the validity of the numerical results obtained with a new radiative transfer scheme like traphic, due to its simpli?cation of the problem. We therefore employ a one-dimensional, explicitly photon-conserving, grid-based radiative transfer scheme that we have implemented ourselves to obtain an accurate numerical reference solution. We will make use of this numerical reference solution in our comparisons below. In the following we present a suite of radiative transfer simulations demonstrating that traphic is able to accurately follow the evolution of the ionization front around a single ionizing point source. For the setup of the numerical simulations we closely follow the description of Test 1 presented in Iliev et al. (2006), the only di?erence being the spatial resolution we employ. This allows us to directly compare our results to the results presented in the code comparison project. In particular, we choose a number density nH = 10?3 cm?3 and an ionizing luminosity ˙ of Nγ = 5 × 1048 photons s?1 . The gas is assumed to have a constant temperature T = 104 K. With these values, rs = 5.4 kpc and τs = 122.4 Myr. We set up the initial conditions in a simulation box of length Lbox = 13.2 kpc containing NSP H = 643 SPH particles12 . The ionizing source is located in the centre. The box boundary is photon-transmissive. We assign each SPH particle a mass m = nH mH L3 /NSP H . The positions of the Box SPH particles are chosen to be glass-like13 . Glass-like initial conditions yield a more regular distribution of the particles within the box as compared to Monte Carlo sampling of the density ?eld. The SPH smoothing kernel is computed and the SPH densities are found using the SPH formalism implemented in gadget-2, with Nngb = 48. We perform 5 simulations, increasing the angular resolution in factors of two from Nc = 8 to Nc = 128. The number of neighbours employed for the transport of radi? ation is ?xed to Nngb = 32. Hence all 5 simulations employ the same spatial resolution. The time step is set to

12 We note that Iliev et al. (2006) employed N 3 cell = 128 cells, with the ionizing source located in one of the corners of the box. 13 The glass-like distribution of particles is achieved by ?rst placing them randomly in the simulation box and thereafter letting them freely evolve under the in?uence of a reversed-sign (i.e. repelling) gravitational force until they settle down into an equilibrium con?guration, see White (1996).

18

Andreas H. Pawlik & Joop Schaye

Eq. 33, but to the pro?le that results after locally averaging the exact pro?le along the radial direction. To investigate the e?ect of particle noise on the neutral fraction pro?le we apply the resampling technique introduced in Section 4.5. We perform a series of 5 simulations that are identical to the simulations presented above, except that every 10th radiative transfer time step the particle positions are randomly perturbed within their SPH spheres of in?uence. The numerical implementation of the resampling is described in Appendix B3. The densities are not recalculated after the positions have been changed due to the resampling, because this would generate ?uctuations in the neutral hydrogen density which would increase the recombination rate and lead to a smaller Str¨mgren sphere. The reo sulting neutral fraction pro?les are shown in the right-hand panel of Fig. 7. All pro?les are now in close agreement with each other and with the exact result. The panels in the bottom row of Fig. 6 show the neutral fraction in a slice through the centre of the simulation box from the simulations with resampling. Clearly, resampling strongly suppresses the particle noise visible in the panels in the top row of Fig. 6, yielding nearly spherical neutral fraction contours. We now investigate the dependence of the equilibrium neutral and ionized fraction pro?le on the spatial resolution by varying NSP H , the number of particles employed in the simulation. Because traphic is explicitly photonconserving, we expect that the radiative transfer in a homogeneous medium is essentially independent of the spatial resolution (see e.g. the discussion in Mellema et al. 2006), except for a trivial averaging. For each of the simulations ? with angular resolution Nc = 8, 32 and 128 and Nngb = 32 presented above, we performed three additional simulations, at lower (NSP H = 163 , 323 ) and higher (NSP H = 1283 ) spatial resolutions, but otherwise identical to the ?ducial (NSP H = 643 ) case. We will focus on the Nc = 8 runs, but note that the Nc = 32 and Nc = 128 series shows the same behaviour. The equilibrium neutral and ionized fraction pro?les are shown in the left-hand panel of Fig. 8. For all spatial resolutions the ionization front is at nearly the same radius. It can furthermore be seen that when the spatial resolution is increased, the equilibrium neutral fraction follows the exact result more closely. The simulation employing the ?ducial spatial resolution (NSP H = 643 ) is almost converged. The di?erences in the neutral fraction pro?les between the simulations using di?erent numbers of particles are fully consistent with the corresponding spatial resolutions, as is demonstrated in the right-hand panel Fig. 8. There, we average the neutral fraction pro?les obtained in the simulations employing NSP H = 323 , 643 and 1283 particles over the spatial resolution element of the lowest resolution simulation (NSP H = 163 ), the size of which is indicated by the vertical line. The averaged pro?les match the neutral fraction pro?le obtained in the low spatial resolution simulation almost exactly. This shows that decreasing the spatial resolution does not introduce any artefacts. The solution obtained by traphic is the converged solution averaged over the adopted spatial resolution. Finally, we show how the size of the time step ?tr a?ects the outcome of our simulations. We again concentrate on the simulation with angular resolution Nc = 8 ? (and NSP H = 643 , Nngb = 32), noting that the simulations

?tr = 104 yr. In Fig. 5 we show the evolution of the ionization front radius, which we determined by taking the average over the positions of all particles that have a neutral fraction 0.4 < η < 0.6. For all 5 simulations, the position of the ionization front never deviates more than 5 per cent from the analytic solution, Eq. 32, comparable to what has been found with other codes as reported in the Cosmological Radiative Transfer Code Comparison Project (Iliev et al. 2006). Note that the deviations from the analytic solution can mainly be attributed to the fact that the analytic approach provides only an approximate expression for the radius of the ionization front, because it erroneously assumes ? χ ≡ 1. In fact, all simulations (except for the Nc = Nngb run, which shows a small deviation of less than two percent, the reason for which will become clear in our discussion below), nearly perfectly follow the numerical reference solution and approach the proper asymptotic limit rI,eq = 1.05 rs . The top row of Fig. 6 shows the neutral fraction in a slice through the centre of the simulation box at tr = 500 Myr, which marks the end of the simulation, for each of the 5 simulations. As we already noted, the ionization front is at the expected position. As is true for all other surfaces of constant neutral fraction shown, the ionization front clearly exhibits the expected spherically symmetric shape, although it is noisy in some of the simulations. The amount of noise depends on the ratio of the angular and spatial resolution employed. For Nc = 8 (left-most panel in the top row), the average number of neighbours per emis? sion or transmission cone is high, Nngb /Nc = 4 and, as a result, numerical noise arising from the representation of the continuous density ?eld with discrete SPH particles are suppressed. With increasing angular resolution the average number of neighbours per cone decreases, and the contours become more noisy. The noise level reaches a maximum at ? Nc = Nngb (middle panel in the top row). For higher angular resolutions, the probability of ?nding no neighbours inside an emission or transmission cone becomes high and a large number of ViPs are created. The ratio of the number of ViPs to the number of SPH particles enclosed by the ionization front for the simulation with angular resolution Nc = 8, 16, 32, 64 and 128 is ≈ 0, 0.003, 0.06, 0.5 and 0.9, resp. The ViPs placed in empty cones regularize the neutral fraction of the ionized density ?eld by distributing the pho? tons they absorb amongst their Nngb SPH neighbours using (photon-conserving) SPH interpolation. In the left-hand panel of Fig. 7 we plot the neutral and ionized fraction averaged in spherical shells as a function of distance to the star, for all 5 simulations. The grey bands indicate the neutral and ionized fraction pro?les obtained with other codes as reported in the cosmological radiative transfer code comparison project (Iliev et al. 2006). Except for the Nc = 32 run, for which the neutral contours were most noisy (see Fig. 6), all simulations agree very well with the results published in the comparison project. The deviations from the exact equilibrium neutral fraction pro?le, i.e. the result of the numerical integration of Eq. 33, can be explained by the particle noise seen in Fig. 6. Due to the fuzzy structure exhibited by the neutral fraction contours, a range of neutral fractions can simultaneously be found within each spherical shell. The pro?les obtained from the numerical simulation with traphic should therefore not be directly compared to the exact pro?le, i.e. the solution of

TRAPHIC - Radiative Transfer for SPH

19

Figure 8. Test 1. Spherically averaged neutral and ionized fractions within the Str¨mgren sphere at tr = 500 Myr. The simulations o ? all have the same angular resolution (Nc = 8) and employ the same number of neighbours (Nngb = 32), but use a di?erent number of SPH particles, increasing from NSP H = 163 to NSP H = 1283 in factors of 23 . The black solid curves corresponds to the exact pro?les, obtained from Eq. 33. The grey bands show the range of neutral and ionized fractions found by other codes as reported in Iliev et al. (2006). Left-hand panel: For each simulation the spatial resolution is marked by a vertical line (in colour and line-style identical to the ? corresponding pro?le) at radius 2 h . The higher the spatial resolution, the more closely the exact pro?le is approached. Right-hand panel: The pro?le corresponding to the lowest spatial resolution simulation (NSP H = 163 , blue dot-dashed line) is repeated from the ? left-hand panel. The pro?les of all other simulations have been averaged over the spatial resolution element 2 h of the lowest spatial resolution simulation, corresponding to the radius marked by the vertical line. The pro?les become nearly identical after smoothing them to the same resolution.

Figure 9. Test 1. E?ect of the size of the time step. We show the evolution of the ionization fronts for four simulations with Nc = 8, ? NSP H = 643 , Nngb = 32 and time steps ?tr = 0.01, 0.1, 1 and 10 Myr, resp, as indicated in the legend. After an initial phase, the evolution of the ionization fronts becomes independent of the size of the radiative transfer time step.

of higher angular resolution exhibit a similar behaviour. In Fig. 9 we show the evolution of the ionization front for four di?erent choices for the size of the radiative transfer time step, ?tr = 0.01, 0.1, 1 and 10 Myr. Note that the simulation with ?tr = 0.01 is identical to the simulation dis-

cussed in Fig. 5. In order to keep the angular sampling the same, at each radiative transfer step we split the emission of photons over 10, 100 and 1000 random orientations of the emission cone tessellation of the ionizing source for the simulations employing ?tr = 0.1, 1 and 10 Myr, resp. Photon packets that are emitted by the source in a certain orientation are transmitted further downstream and can propagate over multiple inter-particle distances. We follow each photon packet until it has either been absorbed or left the simulation box, to properly solve the time-independent radiative transfer equation for the large time steps under consideration. From Fig. 9 we see that the evolution of the ionization front for the simulations with time step ?tr = 0.1, 1 and 10 Myr is delayed with respect to the evolution of the ionization front for the simulation with time step ?tr = 0.01 Myr. This delay increases with the size of the time step, being barely visible for the simulation using ?tr = 0.1 Myr. The delay arises because the neutral fraction is only updated at the end of each radiative transfer time step. Photons that have been absorbed during the transport over a single time step but that have not been used to advance the neutral fraction during the subsequent sub-cycling of the rate equation are therefore only re-inserted in the transport process at the beginning of the next time step14 , and are thus de-

14

We note that we have also tried to re-insert these photons

20

Andreas H. Pawlik & Joop Schaye

an in?nitely thin opaque disc perpendicular to the x-axis at a distance of 0.08 Mpc along the x-axis from the star (thick blue vertical lines in Figs. 10-14). The y and z coordinates of the disc centre are identical to those coordinates of the star. Photons that cross the disc are removed. We performed a series of radiative transfer simulations (with time step ?tr = 104 yr), using di?erent choices for ? the parameters Nngb , which sets the spatial resolution if the total number of SPH particles is ?xed, and Nc , which sets the angular resolution. The results are shown in Fig. 10 - 14, displaying a slice through the simulation box at z = Lbox /2. In each panel, the black dash-dotted line shows the expected position of the ionization front (Eq. 32) at time tr = 80 Myr, which marks the end of the simulation. The black solid lines emerging from the star at the centre of the slice show the boundaries of the shadow expected to be thrown by the opaque disc (thick blue vertical line). In the top-left corner of each panel we indicate the spatial resolution by a cross ? of length 2 h corresponding to the average diameter of the ? sphere containing Nngb neighbours. The angular resolution is indicated by the angle ω enclosed by the black dashed line and the upper shadow boundary, where ω is determined using Eq. 9. It indicates the maximum angle photons can theoretically diverge from the shadow boundary into the shadow region, given the chosen angular resolution Nc . The position of the ionization front (the iso-surface for which the neutral fraction η = 0.5) at time tr = 80 Myr as obtained with traphic is shown by the red solid line. In all panels, that is for all spatial and angular resolutions, the ionization front is found at the proper position, in agreement with our ?ndings of Section 5.3.1. The shadow thrown by the opaque disc is always sharp. We now discuss the dependence of the results on the chosen spatial and angular resolutions. In Fig. 10 we show the ionization front obtained in sim? ulations employing a ?xed spatial resolution, Nngb = 32, but varying angular resolution, ranging from Nc = 8 in the leftmost to Nc = 128 in the right-most panel. The most prominent di?erence between the results of the di?erent simulations is the noisiness of the contour tracing out the ionization front. The angular resolution study presented here is very similar to the one in the last section. For the lowest angular resolution, Nc = 8, the ionization front is very smooth due to the large number of neighbours within each emission and transmission cone. The noise increases with the angular reso? lution until it reaches a maximum for Nc = Nngb . For higher angular resolutions, particle noise is e?ciently suppressed due to the large number of ViPs that are placed in empty cones and that distribute the photons they absorb amongst ? their Nngb SPH neighbours using (photon-conserving) SPH interpolation. From Fig. 10 it can furthermore be seen how the sharpness of the shadow thrown by the opaque disc depends on the angular resolution. For the lowest angular resolution, the shadow is somewhat blurred, though not nearly as much as the angular resolution would imply. Increasing the angular resolution yields slightly sharper shadows. However, if the ? angular resolution is increased beyond Nc = Nngb , the shadows become slightly less sharp. This is because the photons absorbed by ViPs are distributed amongst the neighbouring gas particles using SPH interpolation and the interpolation procedure does not know about the shadow region. The

layed. From Fig. 9 it can, however, be seen that after a few time steps (here ≈ 15), the ionization front catches up to agree with the ionization front obtained in the simulation using ?tr = 0.01 Myr. We have convinced ourselves that from then on the neutral fraction pro?le around the ionizing source is also nearly identical to the pro?le obtained in the simulation with ?tr = 0.01 Myr (Fig. 7). In summary, in this section we showed that traphic is able to reproduce the expected equilibrium neutral fraction around an ionizing source embedded in a homogeneous medium, as well as the dynamical evolution towards it. Because the radiative transfer is explicitly photon-conserving, the spatial resolution only determines the scale over which the converged solution is averaged. We have furthermore seen that the performance of traphic is stable with respect to variations in the size of the time step. Particle noise due to the discrete nature of SPH simulations is small except for ? the choice of parameters Nc = Nngb . The noise can be successfully suppressed by applying the resampling technique of Section 4.5, i.e. by periodically perturbing the positions of the SPH particles within their spheres of in?uence. In the next section we employ a test with broken symmetry to demonstrate the ability of traphic to correctly produce shadows behind opaque obstacles and we study fur? ther how the choice of the parameters Nc and Nngb a?ects the outcome of the numerical simulations. 5.3.2 Test 2: Breaking the spherical symmetry

In the absence of scattering interactions, photons propagate along straight lines into the direction set at the time of their emission. Consequently, opaque obstacles throw sharply de?ned shadows. In this section we are mainly interested in studying the properties of the shadow thrown by an opaque obstacle exposed to ionizing radiation from a single point source, as obtained with traphic. At the same time, we will extend the study of particle noise presented in Section 5.3.1 ? to include other choices for the parameter Nngb . The geometry of our test problem closely follows the description of the shadow test in Mellema et al. (2006). We ˙ consider a source emitting Nγ = 1054 photons s?1 , each of hydrogen-ionizing energy hp ν = 13.6 eV, residing in an initially neutral, static hydrogen-only ?eld of constant number density nH = 1.87×10?4 cm?3 and temperature T = 104 K. The Str¨mgren radius (Eq. 29) corresponding to this set of o parameters is rs = 0.965 Mpc and the Str¨mgren time is o τs = 654.3 Myr. For the numerical simulation a star particle is placed in the centre of a box of size Lbox = 1 Mpc. The boundaries of the box are photon-transmissive. The density ?eld is sampled using NSP H = 643 gas particles with mass m = nH mH L3 /NSP H at glass-like positions. The partiBox cle densities are found using the SPH interpolation implemented in gadget-2, with Nngb = 48. We furthermore place

immediately after they have been absorbed, by integrating the rate equation already at the end of each transport cycle (without updating the neutral fraction) to obtain the number of photons that have been erroneously counted for being absorbed (because of the assumption of a constant neutral fraction). The re-insertion of these photons in the transport process within the same radiative transfer time step over which they have been absorbed indeed reduced the delay of the ionization front observed in Fig. 9.

TRAPHIC - Radiative Transfer for SPH

21

Figure 10. Test 2. Slice through the simulation box at z = Lbox /2 showing the ionization front (red solid line) at time tr = 80 Myr around an ionizing source sitting in the centre of the simulation box. The black dot-dashed line shows the expected ionization front position. The thick blue vertical line indicates an obstacle opaque to ionizing photons and the black solid lines trace out the boundaries of the shadow this obstacle is expected to throw. The cross and the black dashed line indicate the spatial and angular resolution, ? respectively, as described in the text. The spatial resolution is ?xed to NSP H = 643 , Nngb = 32. The angular resolution increases from Nc = 8 in the left-most panel to Nc = 128 in the right-most panel, in factors of 2.

Figure 11. Test 2. Same as Fig. 10, but with the angular resolution ?xed to Nc = 32 and the spatial resolution decreasing, from ? ? NSP H = 643 , Nngb = 8 in the left-most panel to NSP H = 643 , Nngb = 128 in the right-most panel, increasing the number of neighbours ? Nngb in factors of 2.

slight di?usion of photons across the shadow boundary is in this case consistent with the spatial resolution. In Fig. 11 we show the results of the simulations where we ?xed the angular resolution to Nc = 32, but varied the ? spatial resolution by changing Nngb . The trends visible in Fig. 10 can also be observed here. The ionization front is ? most noisy and the shadow is sharpest for Nc = Nngb . For ?ngb > Nc , noise due to the discreteness of the particles N employed for the transport of photons is suppressed by the large number of neighbours per cone, but the shadow is slightly blurred. The shadow becomes sharper for a smaller number of neighbours, since generally not all of the solid angle will be covered by the neighbours, an e?ect that becomes more important for smaller numbers of neighbours. ? For Nngb < Nc the ionization front becomes smoother due to the regularizing e?ect of the SPH interpolation from ViPs, which also leads to a small di?usion of photons across the shadow boundary, consistent with the spatial resolution. ? In Figs. 12 and 13 we keep the ratio Nngb /Nc ?xed at ? ? Nngb /Nc = 2 and Nngb /Nc = 1/2, resp. In the ?rst case there are on average 2 neighbours per cone, whereas in the second case there is on average one neighbour in every second cone. From Fig. 12 it is clear that the shadow does get sharper when the angular resolution is increased, although the e?ect is small, since the shadow is always very sharp. ? Because we keep the ratio Nngb /Nc ?xed at 2, the number

of ViPs employed in the simulation stays low for all angular resolutions and the shadows are not visibly di?used by the SPH interpolation of absorbed photons from the ViPs. Furthermore, the noisiness of the ionization front remains constant throughout the parameter range. This is because ? ? the noise is primarily set by the ratio Nngb /Nc if Nngb > Nc . In Fig. 13, on the other hand, there is a substantial probability for creating a ViP per cone. Since the absolute number of ViPs present in the simulation increases with increasing angular resolution Nc , the noise decreases with Nc . In the last section (see bottom row of Fig. 6) we employed a resampling technique to decrease the noise exhibited by the neutral fraction contours. Recall from Section 4.2.2 that the apexes of the transmission cones are attached to the positions of the SPH particles. Hence, resampling results in slight shifts in the position of each transmission cone apex, on the scale of the spatial resolution employed in the SPH simulation. This shifting is expected to lead to a small di?usion of photons across the expected shadow boundary. As noted earlier in Section 4.5, such a di?usion due to particle resp. apex motion will also occur in radiationhydrodynamical simulations because of the movement of the SPH particles. It is therefore interesting to study the properties of the shadow thrown by an opaque obstacle in the case of resampling. In Fig. 14 we show the results of a simulation which em-

22

Andreas H. Pawlik & Joop Schaye

? Figure 12. Test 2. Same as Fig. 10, but ?xing the ratio between spatial and angular resolution to Nngb /Nc = 2 (NSP H = 643 ).

? Figure 13. Test 2. Same as Fig. 10, but ?xing the ratio between spatial and angular resolution to Nngb /Nc = 1/2 (NSP H = 643 ).

the resampling is performed without changing the hydrogen density, since this would lead to an enhanced recombination rate. Numerical noise is successfully suppressed by the random perturbations given to the positions of the SPH particles. Consequently, the ionization front appears signi?cantly smoother. The degree of photon di?usion into the shadow region is small and does not signi?cantly degrade the angular resolution of the radiative transfer. This is because the di?usion scale is set by the spatial resolution employed in the SPH simulation. Therefore, well-de?ned shadows will be thrown as long as the obstacle is spatially resolved. In summary, in this section we showed that traphic is able to produce a well-de?ned shadow behind an opaque obstacle, with the shadow sharpness in full agreement with the chosen spatial and angular resolutions. In fact, the shadows are much sharper than implied by the formal angular resolution, thanks to the angular adaptivity inherent to traphic. For a ?xed angular resolution, the shadows are sharpest for ? Nc = Nngb . They are slightly broadened by photon di?u? ? sion for both Nc < Nngb and Nc > Nngb , due to the increased coverage of the solid angles traced out by the transmission cones with SPH particles for an increasing number ? of neighbours Nngb and the SPH interpolation of the photons absorbed by ViPs, resp. We con?rmed our ?nding of ? the last section that unless Nngb = Nc , noise due to the discreteness of the particles on which the transport of photons

Figure 14. Test 2. Left-hand panel: Identical to the middle panel of Fig. 10 and 11. Right-hand panel: Same as left-hand panel, except for the fact that in the simulation we periodically resampled the density ?eld, resulting in a strong suppression of the particle noise seen in the left-hand panel. Note the (small amount of) diffusion of photons across the shadow boundary due to the motion of the transmission cone apexes.

ploys the same parameters as used for the simulation presented in the middle panels of Fig. 10 and 11, but with an additional resampling of the hydrogen density ?eld every 10th radiative transfer time step. As in the last section,

TRAPHIC - Radiative Transfer for SPH

takes place is small, since it is suppressed by either the large ? number of neighbours per cone (if Nc < Nngb ) or the large ? number of ViPs employed (if Nc > Nngb ). The resampling technique presented in Section 4.5 is very e?ective at suppressing particle noise. We have seen that resampling the density ?eld does not severely degrade the angular resolution, even though it leads to a small shift of the cone apexes. As long as the opaque obstacle is spatially resolved by the SPH simulation, a well-de?ned shadow will still be thrown. The ability to produce sharp shadows is one of the main requirements a radiative transfer code has to pass. The results of this section, together with the results of Test 1, which showed that traphic is able to reproduce the expected neutral fraction within a spherically symmetric HII region, indicate that traphic can be used to perform the transport of ionizing photons in arbitrary complex geometries. This is the subject of Test 3, which we will describe next.

23

5.3.3

Test 3: Cosmological density ?eld

Figure 15. Test 3: Slice through the density ?eld at z = Lbox /2. Contours show a neutral fraction of η = 0.5 (left-hand column) and η = 0.01 (right-hand column) at times t = 0.05 Myr (top row) and t = 0.2 Myr (bottom row). Red contours show the results of ? our ?ducial (Nc = 32, Nngb = 32) simulation. For comparison, we show the results of c2 -ray (green) and crash (blue), as reported in Iliev et al. (2006). The agreement is excellent.

In this test we study the propagation of ionization fronts around multiple sources in a static cosmological density ?eld. The parameters of this test are taken from Test 4 of the Cosmological Radiative Transfer Code Comparison Project (Iliev et al. 2006). For reference we repeat the test description here. The initial conditions are provided by a snapshot (at redshift z ≈ 8.85) from a cosmological N-body and gasdynamical simulation performed using the cosmological (uniform-mesh) PM+TVD code of Ryu et al. (1993). The simulation box is Lbox = 0.5 h?1 comoving Mpc on a side, uniformly divided into Ncell = 1283 cells. The initial temperature is ?xed at T = 100 K everywhere. The halos in the simulation box were found using a friends-of-friends halo ?nder with a linking length of 0.25. The ionizing sources are chosen to correspond to the 16 most massive halos in the box. We assume that these have a black-body spectrum Bν (ν, T ) with temperature T = 105 K. The ionizing photon production rate is assumed to be constant and assigned assuming that each source lives for ts = 3 Myr and emits fγ = 250 photons per atom during its lifetime. Hence, the number of ionizing photons emitted per unit time is M ?b ˙ Nγ = fγ , ?0 mH ts (35)

where M is the total halo mass, ?0 = 0.27, ?b = 0.043 and h = 0.7. For simplicity, all sources are assumed to switch on at the same time. The boundary conditions are photontransmissive. Outputs are produced at t = 0.05, 0.1, 0.2, 0.3 and 0.4 Myr. With respect to the original test setup described above, we require three changes. First, since our code does not yet solve for the temperature of the gas, we assume a constant temperature of T = 104 K for the ionized gas. Second, since our code currently treats only a single frequency (bin), we assume the grey photo-ionization cross-section Eq. 27, for which we ?nd σ = 1.49×10?18 cm?2 . The third change con? cerns the input density ?eld. Since our code works directly on the set of particles used in SPH simulations, we have to Monte Carlo sample the original input density ?eld in order to place particles in the box. We replace every grid cell i by

i NSP H = Mi /m SPH particles (randomly distributed within the volume of the grid cell), where Mi = ρi L3 /Ncell is the box mass of the cell and m is the mass of an SPH particle. If i NSP H is not an integer, we draw a random number from a uniform distribution on the interval (0,1) and place an additional particle if this number is smaller than the di?eri ence between NSP H and the nearest lower integer. We use NSP H = Ncell = 1283 . Since the Monte Carlo sampling only P i results in the approximate equality i NSP H ≈ NSP H , we adjust the particle P masses a posteriori to conserve mass, i.e. i m → m × NSP H / i NSP H . After the particles have been placed, we calculate their densities using the SPH formalism of gadget-2, with Nngb = 48. Note that Monte Carlo sampling the density ?eld with NSP H ? Ncell particles yields a smaller e?ective resolution than that of the grid input ?eld in low density regions (many grid cells will be left empty of particles), and to a spurious higher resolution in high density regions (cells are sampled with many particles, even though there is no substructure on the scale of a single cell in the input ?eld). Note also that because the initial conditions were speci?ed on a uniform grid, we do not bene?t from the intrinsic spatial adaptivity of traphic, e?ectively wasting computational resources. We performed a set of three radiative transfer simulations with angular resolutions Nc = 8, 32 and 128, which we refer to as the low angular resolution, ?ducial and high angular resolution simulations, resp. Every simulation used ? Nngb = 32 neighbours. The time step was set to ?tr = 100 yr. In Fig. 15 we show in red the neutral fraction contours η = 0.5 (left-hand column) and η = 0.01 (right-hand column) at times t = 0.05 Myr (top row) and t = 0.2 Myr (bot-

24

Andreas H. Pawlik & Joop Schaye

Figure 16. Test 3: E?ect of angular resolution. The same slice as shown in Fig. 15. Contours show a neutral fraction of η = 0.5 (left-hand column) and η = 0.01 (right-hand column) at time t = 0.05 Myr (top row) and t = 0.2 Myr (bottom row). Green, red and blue lines correspond to the low (Nc = 8), ?ducial (Nc = 32) and high (Nc = 128) angular resolution simulations, resp. The numbers of ViPs per SPH particle at the end of the simulations are ≈ 0.05, 1 and 2.8 for the low, ?ducial and high angular resolution simulations, resp. Note that the ?ducial simulation is already converged, even though its angular resolution Nc = 32 corresponds to a relatively large cone opening angle of ω ≈ 41 degrees.

Figure 17. Test 3: E?ect of resampling. The same slice as shown in Fig. 15. Contours show a neutral fraction of η = 0.5 (lefthand column) and η = 0.01 (right-hand column) at time t = 0.05 Myr (top row) and t = 0.2 Myr (bottom row). The red contours correspond to the ?ducial angular resolution simulation and are identical to the red contours in Figs. 15 and 16. The blue contours correspond to a simulation identical to the simulation employing the ?ducial angular resolution, except for the fact that in this simulation we periodically (every 10th radiative transfer time step) resampled the density ?eld to suppress the particle noise. Note that resampling does not visibly decrease the e?ective angular resolution.

tom row) for the ?ducial angular resolution. The ionization front (neutral fraction η = 0.5) is delayed by the dense ?laments, leading to the characteristic ”butter?y” shapes of the ionized regions. For comparison, we also show the results obtained with two other codes, the ray-tracing scheme c2 -ray (Mellema et al. 2006; green contours) and the Monte Carlo code crash (Maselli, Ferrara, & Ciardi 2003; Ciardi et al. 2001; blue contours), as published in the cosmological radiative transfer code comparison project (Iliev et al. 2006)15 . Both c2 -ray and crash are mesh codes, working directly on the uniform mesh input density ?eld provided by the PM+TVD code of Ryu et al. (1993). The agreement between the results of traphic and c2 ray resp. crash is very good, not only for the ionization front, but also for smaller neutral fractions. The contours from traphic are slightly noisier than those from c2 -ray, which is expected since in addition to the particle noise a?ecting the radiative transfer, the Monte-Carlo sampling noise imprinted on the density ?eld a?ects our simulations, particularly in the under-sampled low density regions, as

15

The performance of two more codes was reported in Iliev et al. (2006): ftte (Razoumov & Cardall 2005) and simplex (Ritzerveld & Icke 2006). For clarity and since they are very similar to the results obtained with c2 -ray and crash, we do not show the results obtained with ftte. We do not include the results of simplex in our comparison, since they di?er considerably from those obtained with all other codes.

already noted earlier. The noise level is, however, substantially lower than one would anticipate based on the tests presented in Sections 5.3.1 and 5.3.2. The most likely explanation for this welcome surprise is that the presence of multiple ionizing sources leads to a regularization in the distribution of the neutral fraction. Numerical noise arising from the representation of the continuous density ?eld by a discrete set of particles is therefore reduced. Di?erences between our results and those of c2 -ray resp. crash also arise through the di?erent treatment of the photon spectrum. Since the photo-ionization cross-section depends on frequency (Eq. 15), the thickness of ?nite ionization fronts (e.g. de?ned as 0.9 < η < 0.1) and hence the position of the particular contour η = 0.5 will in part be determined by the details of the numerical implementation of the multifrequency transport. In Fig. 16 we show the neutral fraction contours η = 0.5 (left-hand side) and η = 0.01 (right-hand side) at times t = 0.05 Myr (top row) and t = 0.2 Myr (bottom row) for the low (green contours) and the high (blue contours) angular resolution simulations. For comparison, the contours obtained in the ?ducial simulation are also shown (red). We note that the high angular resolution simulation yields neutral fraction contours that are almost identical to those obtained in our ?ducial simulation, indicating numerical convergence. The low angular resolution simulation, although still in good agreement with the high angular resolution simulation, fails to properly reproduce the expected neutral

TRAPHIC - Radiative Transfer for SPH

fraction contours when scrutinised in detail. In the low angular resolution simulation, neutral fraction contours are often slightly advanced instead of delayed by the dense ?laments. Although the e?ect is small, it becomes apparent when the contours are compared to the corresponding contours of the high angular resolution simulation. The last observation agrees with the discussion of anisotropies in particle-to-neighbour transport schemes sketched in Section 4.2.1 and presented in detail in Appendix A. In Appendix A we demonstrate that when photons are transported from a given particle to its neighbours, the net transport direction is generally strongly correlated with the direction towards the centre of mass of the neighbouring particles. As a result, the transport is partly governed by the spatial distribution of the SPH particles. For cosmological simulations this implies that photons are preferentially transported along dense ?laments. traphic propagates photons in cones to overcome this problem. The cones con?ne the photons to the solid angles they were emitted into, ensuring a correct transport of radiation on the scale of the chosen angular resolution. If the angular resolution is chosen too low to properly resolve the structures in the SPH density ?eld, the transport is no longer independent of the geometry of the SPH simulation and artefacts may occur, as we see in Fig. 16 for the low angular resolution simulation. Even with an angular resolution as low as Nc = 8, however, the artefacts are small. Nevertheless, it is clear that in order to properly solve the radiative transfer equation, the angular resolution must be chosen high enough to establish numerical convergence. Fig. 16 shows that an angular resolution of Nc = 32, which corresponds to a relatively large opening angle of ω ≈ 41 degrees (Eq. 9), is already converged. The reason why a relatively poor formal angular resolution suf?ces, is, as already noted in the discussion of the sharpness of shadows thrown by opaque obstacles in Section 5.3.2, that the photon transport with traphic is intrinsically adaptive in angle. In Fig. 17 we present results for a simulation that used the resampling technique presented in Section B3, but which was otherwise identical to the ?ducial simulation presented above. The neutral hydrogen densities of the SPH particles were not re-calculated according to the perturbed positions resulting from the resampling, but kept constant to avoid additional scatter in the density. The resampling leads to a reduction in the particle noise, which is most apparent for the η = 0.01 contours, since in case of the η = 0.5 the noise is already low, as noted above. Note that the resampling does not noticeable degrade the angular resolution. In Fig. 18 we show the evolution of the mean (over all particles i) ionized fraction, both volume-weighted, i.e. P 3 P 3 χv = ih P P i χi / i hi , and mass-weighted, i.e. χm = i m i χi / i mi . Again, the results obtained with traphic are in excellent agreement with the results obtained with c2 ray and crash. For the latter two we obtained the mean ionized fraction by averaging the ionized fraction reported in the cosmological radiative transfer code comparison project P P (Iliev et al.P 2006) over all grid cells i, i.e. χv = i χi / i 1 P and χm = i ρi χi / i ρi . The ratio of mass-weighted and volume-weighted mean ionized fractions is at early times slightly larger for the low angular resolution simulation than for the high angular resolution simulation, as can be seen in the bottom panel

25

Figure 18. Test 3: The volume- and mass-weighted mean ionized fractions, χV and χm , resp., averaged over the whole simulation box as a function of time, for the low, ?ducial (without and with resampling of the density) and high angular resolution simulation, as indicated in the legend. All results fall nearly on top of each other. Di?erences in χm /χV are only noticeable when χV is small. For comparison, we also show the results obtained with c2 -ray and crash as reported in Iliev et al. (2006).

of Fig. 18. This is another manifestation of the fact that particle-to-neighbour transport is generally not independent of the geometry of the SPH simulation, resulting in photons being preferentially transported into high (particle) density regions. It once more underlines the importance of the concept of emission and transmission cones (with su?ciently small solid angles) which traphic uses to accomplish the transport of radiation independently of the spatial distribution of the SPH particles. Finally, we show that for the present radiative transfer problem simulations solving the time-dependent radiative transfer equation give a signi?cantly di?erent result than the simulations solving the time-independent radiative transfer equation discussed above. We carried out a simulation using ? Nc = 32, Nngb = 32 and additionally employed the photon packet clocks to limit the distance over which photon packets can propagate during each time step. The size of the radiative transfer time step was set to ?tr = 10?3 Myr. This time step is a factor 10 times larger than the time step used for solving the time-independent radiative transfer equation in the simulations presented above, which required a smaller time step because of our particular treatment of the timeindependent radiative transfer equation (see Section 5.2.3). The locations of the ionization fronts (i.e. η = 0.5) obtained in this simulation are shown in Fig. 19 (blue curves), together with those obtained in the corresponding simulation solving the time-independent radiative transfer equation (red curves, cp. the left-hand panels of Fig. 15), at four di?erent simulation times t = 0.05, 0.1, 0.2 and 0.3 Myr. It is clear from Fig. 19 that the simulation solving the time-independent radiative transfer equation produces ionized spheres that are unphysically large. This is due to the well-known fact (see, e.g., the discussion in Abel, Norman, & Madau 1999) that ionization fronts may propagate at speeds larger than the speed of light, if the time-independent radiative transfer equation is solved. The di?erence between the two simulations is larger at early

26

Andreas H. Pawlik & Joop Schaye

Figure 19. Test 3: Slice through the density ?eld at z = Lbox /2. Contours show a neutral fraction of η = 0.5 at times t = 0.05, 0.1, 0.2 ? and 0.3 Myr (from left to right). Red contours show the results of our ?ducial (Nc = 32, Nngb = 32) simulation (at times 0.05 and 0.2 Myr these are identical to the red contours shown in the top-left resp. bottom-left panel of Fig. 15). Blue contours show the result of ? a simulation employing the same resolution (Nc = 32, Nngb = 32), but using the clock of the photon packets to solve the time-dependent radiative transfer equation. Note that in the simulation solving the time-independent radiative transfer equation the ionized regions are too large, since in this simulation the ionization fronts are initially propagating at speeds larger than the speed of light.

times than at late times, which is expected, since in equilibrium, i.e. tr → ∞, the results of both simulations must agree. In summary, in this section we studied the propagation of ionization fronts around multiple sources in a static cosmological density ?eld. We demonstrated the importance of the concept of cones which underlies the photon transport in traphic. Without the con?nement by transmission cones of su?ciently small solid angle, particle-to-neighbour transport is governed in part by the spatial distribution of the particles, resulting in the preferential transport of photons into high (particle) density regions. Thanks to the fact that traphic is adaptive in angle, a relatively modest formal angular resolution of Nc = 32 is already su?cient to obtain a converged solution. Since the setup of our simulations followed the description of the corresponding test in the cosmological radiative transfer code comparison project (Iliev et al. 2006), we were able to benchmark our radiative transfer scheme by comparing with the results obtained by the ray-tracing scheme c2 -ray (Mellema et al. 2006) and the Monte Carlo code crash (Maselli, Ferrara, & Ciardi 2003; Ciardi et al. 2001). We found excellent agreement in the positions of neutral fraction contours as well as the mass and volume-weighted mean ionized fractions. We have furthermore seen that for the test problem presented in this section, simulations solving the timeindependent radiative transfer equation lead to ionized regions that are unphysically large during their early evolution. This illustrates the importance of correctly accounting for the ?nite speed of light when performing radiative transfer simulations to study the morphology of ionized regions that are strongly out of equilibrium.

Due to the limits on the computational power available today, radiative transfer simulations at the resolution of stateof-the-art hydrodynamical simulations need to be adaptive both in space and in angle and parallel on distributed memory machines. traphic meets these requirements. The radiative transfer equation is solved by transporting photons in an explicitly photon-conserving manner directly on the discrete set of SPH particles. Photons are transported globally by propagating them locally, between a particle and its neighbours. For a ?xed number of SPH particles, the number of neigh? bours employed for the transport, the parameter Nngb , sets the adaptive spatial resolution. The conceptual similarity between the photon transport and the calculation of particle properties in SPH simulations allows a straightforward numerical implementation of traphic in SPH simulations for application on distributed memory machines. In traphic photon packets are transported inside cones. Because source particles emit into a set of cones that tessellates the simulation box, the angular dependence of the emission can be modelled independently of the spatial distribution of the neighbouring SPH particles. The number of cones, Nc , is the second parameter employed in traphic and determines the formal angular resolution of the radiative transfer simulation. After their emission, the photon packets are transported downstream from SPH particle to SPH particle. They are con?ned to the solid angle they were originally emitted into by regular cones of solid angle 4π/Nc . Mimicking the ray-splitting technique used in ray-tracing schemes, the photon transport is adaptive in angle. The e?ective angular resolution is therefore much higher than the formal one. Radiative transfer simulations employing traphic can thus be performed using only a relatively small number of cones without loss of accuracy. The propagation of photons inside cones that do not contain a neighbour requires the introduction of so-called virtual particles (ViPs). It is the concept of cones combined with ViPs that enables the directed transport of photons on the unstructured grid de?ned by the discrete set of irregularly distributed SPH particles. In addition to the number of ? cones Nc and the number of SPH neighbours Nngb , it is the ?ngb that directly in?uences the performance of ratio Nc /N

6

CONCLUSIONS

In this paper we have presented traphic (TRAnsport of PHotons In Cones), a novel radiative transfer scheme for SPH simulations. traphic has been designed for use in radiation-hydrodynamical simulations exhibiting a wide dynamic range of physical length scales and containing a large number of light sources, as for instance encountered in cosmological simulations. The requirements for this are high:

TRAPHIC - Radiative Transfer for SPH

traphic. It controls the amount of noise introduced by the representation of the underlying continuum physics with a discrete set of particles. This particle noise is small for both ? ? Nc < Nngb and Nc > Nngb due to the large number of neighbours per cone and the large number of ViPs, resp. ? For the choice of parameters Nc = Nngb the particle noise can be substantial. It can, however, be e?ciently suppressed by employing a density ?eld resampling strategy. A practical problem often encountered when solving the radiative transfer equation is the linear scaling of the computational e?ort with the number of sources. Such a scaling limits the application of most of the existing radiative transfer schemes to simulations containing only a few sources. In traphic this problem is circumvented by introducing a source merging procedure. Photon packets emitted by different sources are merged in accordance with the angular resolution employed, such that at any point in space at most Nc photon packets need to be propagated. This renders the photon transport e?ectively independent of the number of sources in the simulation and makes it feasible to include a di?use radiation component emitted by the SPH particles. In the most extreme case of radiation from an arbitrary number of sources completely ?lling the simulation box, the computational e?ort required by traphic merely scales as ? the product of spatial and angular resolution, N × Nngb ×Nc , where N is the total number of particles (SPH + stars). Applications of traphic include the solution of both the time-independent and time-dependent radiative transfer equation in large-scale simulations of cosmological reionization. Here we have speci?ed traphic for the transport of mono-chromatic hydrogen-ionizing radiation and described its implementation in the parallel SPH code gadget-2 (Springel 2005). We showed that traphic is able to accurately reproduce the expected growth of the ionized sphere around a single point source in a homogeneous medium and to cast sharp shadows behind opaque obstacles. Furthermore, we tested our scheme in a physical setting of complex geometry: The growth of ionized regions around multiple point sources in a cosmological density ?eld. Detailed comparisons with the results obtained by other codes (Iliev et al. 2006) for identical problems show excellent agreement. In this paper we have limited ourselves to the description of the radiative transfer scheme traphic and its implementation for use on static density ?elds. However, traphic can by design readily be coupled to the thermo- and hydrodynamic evolution of matter in SPH simulations. We will report on such a coupling and its extension to the transport of multi-frequency radiation in future work.

27

Groningen. This work was supported by Marie Curie Excellence Grant MEXT-CT-2004-014112.

REFERENCES Abel T., Norman M. L., Madau P., 1999, ApJ, 523, 66 Abel T., Wandelt B. D., 2002, MNRAS, 330, L53 Alvarez M. A., Bromm V., Shapiro P. R., 2006, ApJ, 639, 621 Barkana R., Loeb A., 2004, ApJ, 609, 474 Brookshaw L., 1994, MmSAI, 65, 1033 Cen R., 2002, ApJS, 141, 211 Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS, 324, 381 Croft R. A. C., Altay G., 2007, arXiv, 709, arXiv:0709.2362 Dale J. E., Ercolano B., Clarke C. J., 2007, MNRAS, 1056 Dale J. E., Bonnell I. A., Whitworth A. P., 2007, MNRAS, 375, 1291 Dopita M. A., Sutherland R. S., 2003, Astrophysics of the di?use universe, Springer Fryer C. L., Rockefeller G., Warren M. S., 2006, ApJ, 643, 292 Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375 Gingold R. A., Monaghan J. J., 1978, MNRAS, 184, 481 Gingold R. A., Monaghan J. J., 1982, J. Comput. Phys., 46, 429 Gnedin N. Y., Abel T., 2001, NewA, 6, 437 Goldstein H., 1980, Classical Mechanics, Addison-Wesley Publ. Co Gritschneder M., Naab T., Heitsch F., Burkert A., 2007, IAUS, 237, 246 Herant M., Benz W., Hix W. R., Fryer C. L., Colgate S. A., 1994, ApJ, 435, 339 Hernquist L., Katz N., 1989, ApJS, 70, 419 Hockney R. W., Eastwood J. W., 1988, Computer Simulations Using Particles, Taylor & Francis Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006, MNRAS, 369, 1625 Iliev I. T., et al., 2006, MNRAS, 371, 1057 Johnson J. L., Greif T. H., Bromm V., 2007, ApJ, 665, 85 Kessel-Deynet O., Burkert A., 2000, MNRAS, 315, 713 Kohler K., Gnedin N. Y., Hamilton A. J. S., 2007, ApJ, 657, 15 Lucy L. B., 1977, AJ, 82, 1013 Maselli A., Ferrara A., Ciardi B., 2003, MNRAS, 345, 379 Mayer L., Lufkin G., Quinn T., Wadsley J., 2007, ApJ, 661, L77 Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, NewA, 11, 374 Mihalas D., Weibel Mihalas B., 1984, Foundations of radiation hydrodynamics, Oxford University Press, New York Miles R. E., 1965, Biometrika, Vol. 52, No. 3/4, pp. 636 Monaghan J. J., 1992, ARA&A, 30, 543 Monaghan J. J., Rep. Prog. Phys., 68, 1703 (2005) Okabe A., Boots B., Sugihara K., Chiu S. N., 2000, Spatial Tessellations, Second edition, Wiley Osterbrock D. E., 1989, Astrophysics of gaseous nebulae and active galactic nuclei, Palgrave Macmillan Oxley S., Woolfson M. M., 2003, MNRAS, 343, 900 Paschos P., Norman M. L., Bordner J. O., Harkness R., 2007, arXiv, 711, arXiv:0711.1904

ACKNOWLEDGMENTS We thank Claudio Dalla Vecchia for stimulating discussions, enduring encouragement as well as technical help and Huub R¨ttgering for his valuable advice. We are grateful to Craig o Booth and Benedetta Ciardi for a thorough reading of the draft. A.H.P. thanks the University of Potsdam and the Max Planck Institut f¨r Astrophysik for their hospitality. Some u of the simulations presented here were run on the Cosmology Machine at the Institute for Computational Cosmology in Durham as part of the Virgo Consortium research programme and on Stella, the LOFAR BlueGene/L system in

28

Andreas H. Pawlik & Joop Schaye

both on the weights and on the spatial distribution of the neighbours. To study the angular dependence P the emission proof P cess, we introduce the vector sum s = i wi ui / i wi over all unit vectors ui ≡ ri /|ri |, where ri is the position of neighbour i. This vector sum can be interpreted as the net direction of the emission. Let us denote the angle enclosed by s and the vector to the centre of mass rcm of the neighP ? bours, rcm = i ri /Nngb , with α. We ask for the probability density function (pdf) p(cos α) of the cosine of α, where P P rj s · rcm i wi ui · Pj cos α = = P . (A1) |s| |rcm | | i wi ui || j rj | The meaning of the pdf can be understood by looking at its extremes. If p(cos α) is strongly peaked around cos α ≈ 1, the emission has a net direction biased towards the centre of mass direction, and photons will be preferentially transported into the high (particle) density regions. On the other hand, if p(cos α) is strongly peaked around cos α ≈ ?1, photons will be preferentially transported into the direction opposite to the centre of mass. Finally, a ?at pdf p(cos α) = 1/2 indicates that the net emission direction and the direction towards the centre of mass are uncorrelated. As an illustration, assume that all neighbours have identical weights wi = 1 and that their positions ri are drawn from a probability distribution that uniformly samples the ? surface of the sphere of radius h surrounding the emitter at O. To obtain p(cos α), we note that our assumptions imply that the centre of mass vector rcm and the vector sum s point in the same direction, so that p(cos α) = 2δD (cos α ? 1), where δD (x) is the Dirac δ -function. Hence, there is a perfect correlation between s and rcm . Another example is given by assuming that the particles of unit mass are distributed randomly within the sphere of ? radius h around the emitter. The resulting pdf p(cos α) as obtained through a Monte Carlo simulation is shown in its R cos α cumulative form F (cos α) = ?1 dx p(x) in Fig. A1 for ? two choices of the number of neighbours, Nngb = 16 (thin ? dotted curve) and Nngb = 64 (thick dotted curve; falling nearly on top of the thin dotted curve). Clearly, there is a very strong correlation between the net direction of emission and the direction towards the centre of mass of the neighbour distribution. For the important case of isotropic sources of photons the net emission direction should not be correlated with the direction towards the centre of mass, i.e. p(cos α) = 1/2, independent of cos α, and thus F (cos α) = (1 + cos α)/2. From here on, we will refer to this case as isotropic emission of radiation. We emphasize the statistical character of this statement. An individual particle may still emit anisotropically, which could be revealed by studying the amplitude |s|. For emission that is isotropic for individual particles, |s| = 0. Consider the neighbour distribution of the last example, for which we showed the emission pattern to be far from isotropic. Can the emission be isotropized by tuning the weights? For instance, consider solid angle weighting, de?ned as follows. We project all neighbours radially onto the unit sphere. To each neighbour we attach a weight proportional to its area of in?uence, which we de?ne by the collection of

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in C. The art of scienti?c computing, Cambridge University Press, 2nd ed. Price D., 2005, astro, arXiv:astro-ph/0507472 Razoumov A. O., Cardall C. Y., 2005, MNRAS, 362, 1413 Razoumov A. O., Sommer-Larsen J., 2006, ApJ, 651, L89 Ritzerveld J., Icke V., 2006, PhRvE, 74, 026704 Ryu D., Ostriker J. P., Kang H., Cen R., 1993, ApJ, 414, 1 Semelin B., Combes F., Baek S., 2007, arXiv, 707, arXiv:0707.2483 Soneira R. M., Peebles P. J. E., 1978, AJ, 83, 845 Springel V., Hernquist L., 2002, MNRAS, 333, 649 Springel V., 2005, MNRAS, 364, 1105 Stamatellos D., Whitworth A. P., 2005, A&A, 439, 153 Stewart G. W., 1980, SIAM Journal on Numerical Analysis, Vol. 17, No. 3, 403 Susa H., 2006, PASJ, 58, 445 Susa H., Umemura M., 2006, ApJ, 645, L93 Trac H., Cen R., 2006, astro, arXiv:astro-ph/0612406 Viau S., Bastien P., Cha S.-H., 2006, ApJ, 639, 559 Whalen D., Norman M. L., 2007, arXiv, 708, arXiv:0708.2444 Wise J. H., Abel T., 2007, arXiv, 710, arXiv:0710.3160 White S. D. M., 1996, Cosmology and large scale structure, Proceedings of the ”Les Houches Ecole d’Ete de Physique Theorique” (Les Houches Summer School), p. 349, Elsevier Scienti?c Publishing Company, Amsterdam Whitehouse S. C., Bate M. R., 2006, MNRAS, 367, 32 Yoshida N., Oh S. P., Kitayama T., Hernquist L., 2007, ApJ, 663, 687

APPENDIX A: THE ANISOTROPY OF PARTICLE-TO-NEIGHBOR TRANSPORT In Section 4.2.1 we introduced an emission cone tessellation to accomplish the emission of photons by a source particle to its neighboring gas particles. Then we also mentioned that the cones are required to achieve an emission in agreement with the angular dependence of the emissivity of the source, independently of the spatial distribution of the neighbors. In this appendix we explain why this is the case. In particular, assuming that particles have equal mass, we show that when transporting photons (or more generally, any extensive quantity) in particle simulations from a particle to its neighbors, the net transport direction is generally correlated with the direction towards the centre of mass of the neighbouring particles. As a result, the transport is partly governed by the geometry of the simulation, in addition to the intrinsic emissivity of the source. Consider a particle located at the origin O of a 3dimensional coordinate system (the emitter frame). The par? ticle emits photons to its Nngb nearest neighbours, residing ? centred on the emitting particle. in the sphere of radius h Throughout this section we assume that all particles have equal mass, m = 1. The emission is performed by transfer? PNngb ring the fraction wi / j=1 wj of the total number of pho? tons to neighbour i (i = 1, 2, ..., Nngb ), where the wi 0 are weights to be speci?ed. The emission properties will depend

TRAPHIC - Radiative Transfer for SPH

all points on the unit sphere which are closer to the projected neighbour for which we calculate the weight than to all other projected neighbours. Hereby, the distance of two points on the sphere is given by the arclength of the connecting geodesic (great circle). This de?nition of the area of in?uence results in a speci?c tessellation of the surface of the sphere, the so-called Voronoi tessellation (e.g. Okabe et al. 2000), and provides us with well-de?ned non-overlapping solid angles. The resulting cumulative distribution F (cos α) is shown in Fig. A1 (dashed curves). Although the weighted emission has weakened the correlation between s and rcm as compared to the case of equal weights, the net emission direction is still markedly peaked towards the direction to the centre of mass. We note that weighting by solid angle is the most natural choice and its inability to isotropize the transport might be surprising. Alternatively, one could introduce ad r hoc weights wi = wi , where r > 1 is some exponent and ? the wi are solid angle weights, to re-shape the distribution p(cos α) to become less and less peaked around cos α ? 1. However, the reason why the solid angle weights fail to isotropize the transport is not that the weights are inappropriate, but that the directions to the neighbours are ?xed and thus cannot be freely chosen along with the weights, as mentioned in Section 4.2.1. This strongly limits one’s ability to in?uence the shape of p(cos α). For illustration, consider an emitter having all its neighbours located in one hemisphere. The net emission direction s will then necessarily lie in that hemisphere, too, regardless of the chosen weights. For a statistical ensemble of emitters, s will then necessarily be correlated with the direction towards the centre of mass. This also explains the dependence of p(cos α) on the ? number of neighbours. The larger Nngb , the greater the probability of ?nding a neighbour in any chosen solid angle. Through this increase in directional sampling, p(cos α) can be expected to approach the uncorrelated case. This behaviour can indeed be observed in Fig. A1. We also compute the statistic p(cos α) for the case of neighbours clustered in space, which is more typical for cosmological simulations of structure formation. We did this both using a toy clustering model (Soneira & Peebles 1978) and for a cosmological density ?eld obtained from a ΛCDM hydrodynamical simulation (the same ?eld that is used in Section 5.3.3). For the weighting schemes we investigated the dependence of the shape of p(cos α) on the chosen weights is qualitatively very similar to that in the case of the random distribution of neighbours discussed earlier. We emphasize that our discussion is not limited to the transport of photons. Indeed, particle-to-neighbour transport is for example routinely employed in SPH simulations of galaxy formation, where a star particle distributes its metals and energy amongst its neighbours. Di?erent recipes for distributing metals and energy can be found in the literature. Most of these apply a weighting related to the SPH formalism. As an example, to each neighbour i we associate the weight wi = P

mi ? W (ri , h) ρi . mj ? j ρ W (rj , h)

j

29

Figure A1. Monte Carlo simulation of particle-to-neighbour transport. We show the cumulative probability distribution of cos α, where α is the angle between the net emission direction s and the direction towards the centre of mass rcm of the neighbours. The curves correspond to the following weighting schemes, with the number of neighbours indicated in brackets (from top to bottom): SPH (64), SPH (16), Solidangle (64), Solidangle (16), Equal (64), Equal (16). The last two curves are nearly on top of each other.

(A2)

? Here, h is the radius of the neighbourhood of the emitting particle (at O), ri = |ri | and W is the spline kernel Eq. 4.

Performing the transport on the same random neighbour distribution as before, we obtain the distribution p(cos α). As can be seen in Fig. A1 (solid curves), the direction of net transport and the direction to the centre of mass are found to be almost uncorrelated. We note that this result, the reason of which at the moment is not clear to us, is only a statistical statement. We calculated |s| and found that the emission of individual particles is still not close to isotropic. We arrive at qualitatively similar results for di?erent de?nitions of the SPH weights, both for random and clustered neighbour distributions. A more quantitative statement, however, must depend on both the speci?c properties of the spatial distribution of the neighbours, and the precise form of the adopted SPH weights. In summary, in this appendix we have shown that when transporting a quantity from a particle to its neighbours, the net transport direction is generally governed by the spatial distribution of the neighbours, in addition to the the intrinsic properties of the emitting particle. We emphasize that our observations apply to the particle-to-neighbour transport of arbitrary quantities in any particle simulation. We would like to remind the reader that in Section 4 we presented a speci?cally designed particle-based transport scheme, traphic, that overcomes the problems of particleto-neighbour based transport discussed in this appendix, not just statistically, but also on a particle by particle basis. We achieved this by employing cones to con?ne photons to the solid angle into which they are emitted, making the transport independent of the geometry of the SPH simulation, on the scale of the chosen angular resolution.

30

Andreas H. Pawlik & Joop Schaye

tessellation has a random orientation. The primary motivation for randomly rotating cones is to increase the angular sampling. Furthermore, randomly rotating the emission and reception cones leads to a reduction of artefacts arising from the particular de?nition we employ to construct the cone tessellation, as noted in the last section. Here we describe our numerical implementation of the random rotations. We can think of the set of cones that comprises a particular cone tessellation as a rigid body, to which we can attach a local Cartesian coordinate system with axes {e′ , e′ , e′ }. z x y The orientation of this coordinate system with respect to the canonical Cartesian coordinate system, e.g. the simulation box axes {ex , ey , ez }, is fully described by three variables, the Eulerian angles (e.g. Goldstein 1980). Eulerian angles are de?ned as the three successive angles of rotations that map the axes {ex , ey , ez } onto the axes {e′ , e′ , e′ }. There exist several conventions. In the x y z zxz convention we employ here, the initial system of axes {ex , ey , ez } is ?rst rotated counter-clockwise by an angle φ around the z-axis, with the resulting coordinate system labelled {eξ , eη , eζ }. Second, the coordinate system {eξ , eη , eζ } is rotated by an angle θ counter-clockwise about the ξ-axis, leaving the new coordinate system {e′ , e′ , e′ }. The third η ζ ξ and last rotation is carried out counter-clockwise by an angle ψ around the ζ ′ -axis, giving the desired {e′ , e′ , e′ } cox y z ordinate system. To obtain random Eulerian angles, we note that the invariant measure ? (the “volume element”) on SO(3), the group of proper rotations in R3 , in the zxz Eulerian angle parametrization reads (e.g. Miles 1965), ?(φ, θ, ψ)dφdθdψ = (B7) 1 sin θdφdθdψ. 8π (B10)

APPENDIX B: NUMERICAL IMPLEMENTATION B1 Cones

In this section we describe the numerical implementation of the cones employed for the emission and reception of photon packets in traphic. For the emission (Section 4.2.1), each source particle divides its neighbourhood using a set of tessellating emission cones. The same tessellation is also employed for the reception of photon packets by gas particles (Section 4.2.3). In the following, we employ spherical coordinates (r, φ, θ), which are related to the Cartesian components (rx , ry , rz ) of an arbitrary vector r through rx ry rz = = = r cos φ sin θ r sin φ sin θ r cos θ. (B1) (B2) (B3)

In our implementation, the emission (reception) cones sample the volume around each particle isotropically. Since the surface element of the unit sphere is given by ds = d(cos θ)dφ, this is achieved by distributing the cone boundaries16 uniformly (i.e. on a regular lattice with indices i, j) in (cos θ, φ). Thus, the boundaries of cone (i, j) are described by the four arcs of constant φij 1 φij 2

ij θ1 ij θ2

= = = =

i

2π , Nφ

0 2π , Nφ

i < Nφ , 0 i < Nφ ,

(B4) (B5) (B6)

(i + 1)

j ), 0 j < Nθ , Nθ j+1 arccos(1 ? 2 ), 0 j < Nθ . Nθ arccos(1 ? 2

Correspondingly, we de?ne the emission (reception) cone axes by φij φij + φij 1 2 , (B8) 2 ij ij θ1 + θ2 θij = , (B9) 2 Note that each of the Nc = Nφ × Nθ emission (reception) cones has the same solid angle ? = 4π/Nc , as can be seen from integrating over the surface element of the unit sphere within the boundaries (Eqs. B4-B7) of a single cone. We also implemented the tessellation used in Abel, Norman, & Madau (1999), which leads to cones that are more similar in shape. We could not ?nd any systematic di?erences in the test problems described in Sections 5.3.1 - 5.3.3 using either of the two types of tessellations. This is not surprising because any artefacts due to the shape of the cones will be suppressed by the random rotations of the emission (reception) cones we perform before each emission (reception), as we will describe in the following section. = B2 Random rotations

Random Eulerian angles are therefore obtained by drawing random variables u1 , u2 , u3 from a uniform distribution on the interval [0, 1] and applying the usual transformation (cp. Press et al. 1992), φ θ ψ = = = 2πu1 arccos(1 ? 2u2 ) 2πu3 . (B11) (B12) (B13)

We implement random rotations using rotation matrices, which are obtained from the random Eulerian angles. The matrix elements of a matrix representing a rotation associated with a given set of Eulerian angles can be readily calculated. For details, we refer the reader to Goldstein (1980). In principle, random rotation matrices can be obtained without drawing random Eulerian angles (e.g. Stewart 1980). However, we ?nd that it is faster to generate random Eulerian angles, and then calculate the corresponding rotation matrices. Moreover, storing the three Eulerian angles instead of the nine rotation matrix elements reduces the memory cost.

B3

Reduction of particle noise

Recall from Sections 4.2.2 and 4.2.3 that we apply a random rotation to each cone tessellation. Consequently, each cone

16

Strictly speaking, one should distribute the cone axes uniformly in (cos θ, φ), but this implies asymmetric cones.

For some of the simulations presented in Sections 5.3.1 5.3.3 we employ the resampling approach described in Section 4.5 to reduce numerical artefacts arising from the representation of the continuous density ?eld with a discrete set of particles. The resampling requires the sampling of the

TRAPHIC - Radiative Transfer for SPH

SPH kernel used in gadget-2, which is the spline given by Eq. 4. We approximate this kernel by a Gaussian, Wσ (r) = 1 exp(?r 2 /2σ 2 ). (2π)3/2 σ 3 (B14)

31

We ?nd that with a standard deviation of σ ? 0.29 × h, the Gaussian provides a reasonable ?t to the spline. A similar relation was employed in Alvarez, Bromm, & Shapiro (2006). When we resample the SPH density ?eld, all SPH particles are redistributed by randomly displacing them from the positions given by the hydrodynamical simulation, with the probability to ?nd a given particle displaced by the distance r given by Eq. B14. The (rare) displacements larger than h are discarded; in this case the original particle positions as determined by the hydrodynamical simulation are used.

A This paper has been typeset from a TEX/ L TEX ?le prepared by the author.