Enzo+Moray: Radiation Hydrodynamics Adaptive Mesh Refinement Simulations with Adaptive Ray Tracing
We describe a photon-conserving radiative transfer algorithm, using a spatially-adaptive ray tracing scheme, and its parallel implementation into the adaptive mesh refinement (AMR) cosmological hydrodynamics code, enzo. By coupling the solver with the energy equation and non-equilibrium chemistry network, our radiation hydrodynamics framework can be utilised to study a broad range of astrophysical problems, such as stellar and black hole (BH) feedback. Inaccuracies can arise from large timesteps and poor sampling, therefore we devised an adaptive time-stepping scheme and a fast approximation of the optically-thin radiation field with multiple sources. We test the method with several radiative transfer and radiation hydrodynamics tests that are given in Iliev et al. (2006); Iliev et al. (2009). We further test our method with more dynamical situations, for example, the propagation of an ionisation front through a Rayleigh-Taylor instability, time-varying luminosities, and collimated radiation. The test suite also includes an expanding H ii region in a magnetised medium, utilising the newly implemented magnetohydrodynamics module in enzo. This method linearly scales with the number of point sources and number of grid cells. Our implementation is scalable to 512 processors on distributed memory machines and can include radiation pressure and secondary ionisations from X-ray radiation. It is included in the newest public release of enzo.
keywords:cosmology — methods: numerical — hydrodynamics — radiative transfer.
Radiation from stars and black holes strongly affects their surroundings and plays a crucial role in topics such as stellar atmospheres, the interstellar medium (ISM), star formation, galaxy formation, supernovae (SNe) and cosmology. It is a well-studied problem (e.g. Mathews, 1965; Rybicki & Lightman, 1979; Mihalas & Mihalas, 1984; Yorke, 1986); however, its treatment in multi-dimensional calculations is difficult because of the dependence on seven variables – three spatial, two angular, frequency, and time. The non-local nature of the thermal and hydrodynamical response to radiation sources further adds to the difficulty. Depending on the problem of interest some simplifying assumptions may be made.
An important case was considered by Strömgren (1939) for an ultraviolet (UV) radiation source photo-ionising a static uniform neutral medium. When recombinations balance photo-ionisations, the radius of a so-called H ii region,
where is the ionising photon luminosity, is the recombination rate, and is the ambient hydrogen number density. Furthermore he found that the delineation between the neutral and ionised medium to be approximately the mean free path of the ionising radiation. His seminal work was expanded upon by Spitzer (1948, 1949, 1954) and Spitzer & Savedoff (1950), who showed that the ionising radiation heated the medium to K. If the density is equal on both sides of the ionisation front, then this over-pressurised region would expand and drive a shock outwards (e.g. Oort, 1954; Schatzman & Kahn, 1955). These early works provided the basis for the modern topic of radiation hydrodynamics of the ISM. A decade later, the first radiation hydrodynamical numerical models of H ii regions in spherical symmetry and plane-parallel ionisation fronts were developed (e.g. Mathews, 1965; Lasker, 1966; Hjellming, 1966). They described the expansion of the ionisation front and the evolution of its associated shock wave that carries most of the gas away from the source. At the same time, theoretical models of ionisation fronts matured and were classified by Kahn (1954) and Axford (1961) as either R-type (rare) or D-type (dense). In R-type fronts, the ionised gas density is higher than the neutral gas density, and in D-type fronts, the opposite is true. R-type fronts travel supersonically with respect to the neutral gas, whereas D-type fronts are subsonic. Furthermore “weak” and “strong” R-type fronts move supersonically and subsonically with respect to the ionised gas, respectively. The same terminology conversely applies to D-type fronts. “Critical” fronts are defined as moving exactly at the sound speed. These works established the evolutionary track of an expanding H ii illuminated by a massive star in a uniform medium:
Weak R-type: When the star (gradually) starts to shine, the ionisation front will move supersonically through the ambient medium. The gas is heated and ionised, but otherwise left undisturbed. This stage continues until .
Critical R-type: As the ionisation front moves outwards, it begins to slow because of the geometric dilution of the radiation. It becomes a critical R-type front, which is equivalent to an isothermal shock in the neutral gas.
Strong and weak D-type: The front continues to slow, becoming a strong D-type front, and then a critical D-type front. From this point forward, it is moving subsonically with respect to the ionised gas, i.e. a weak D-type front. Thus sound waves can travel across the ionisation front and form a shock. The ionisation front detaches from the shock, putting the shock ahead of the ionisation front.
Expansion phase: After the shock forms, the H ii region starts to expand, lowering the interior density and thus the recombination rates. This increases the number of photons available for ionising the gas. The sphere expands until it reaches pressure equilibrium with the ambient medium at .
In the 1970’s and 1980’s, algorithmic and computational advances allowed numerical models to be expanded to two dimensions, mainly using axi-symmetric to simplify the problem (e.g. Bodenheimer et al., 1979; Sandford et al., 1982; Yorke et al., 1983). One topic that was studied extensively were champagne flows. Here the source is embedded in an overdense region, and the H ii region escapes from this region in one direction. The interface between the ambient and dense medium was usually set up to be a constant pressure boundary. When the ionisation front passes this boundary, the dense, ionised gas is orders of magnitude out of pressure equilibrium as the temperatures on both sides of initial boundary are within a factor of a few. In response, the gas is accelerated outwards in this direction, creating a fan-shaped outflow.
Only in the past 15 years, computational resources have become large enough, along with further algorithmic advances, to cope with the requirements of three-dimensional calculations. There are two popular methods to solve the radiative transfer equation in three-dimensions:
Moment methods: The angular moments of the radiation field can describe its angular structure, which are related to the energy energy, flux, and radiation pressure (Auer & Mihalas, 1970; Norman et al., 1998). These have been implemented in conjunction with short characteristics (Stone et al., 1992, 2D), with long characteristics (Finlator et al., 2009), with a variable Eddington tensor in the optically thin limit (OTVET; Gnedin & Abel, 2001; Petkova & Springel, 2009), and with an M1 closure relation (González et al., 2007; Aubert & Teyssier, 2008). Moment methods have the advantage of being fast and independent of the number of radiation sources. However, they are diffusive and result in incorrect shadows in some situations.
Ray tracing: Radiation can be propagated along rays that extend through the computational grid (e.g. Razoumov & Scott, 1999; Abel et al., 1999; Ciardi et al., 2001; Sokasian et al., 2001; Whalen & Norman, 2006; Rijkhorst et al., 2006; Mellema et al., 2006; Alvarez et al., 2006a; Trac & Cen, 2007; Krumholz et al., 2007; Paardekooper et al., 2010) or particle set (e.g. Susa, 2006; Johnson et al., 2007; Pawlik & Schaye, 2008, 2010; Altay et al., 2008; Hasegawa et al., 2009). In general, these methods are very accurate but computationally expensive because the radiation field must be well sampled by the rays with respect to the spatial resolution of the grid or particles.
Until the mid-2000’s the vast majority of the three-dimensional calculations were performed with static density distributions. One example is calculating cosmological reionisation by post-processing of density fields from N-body simulations (Ciardi et al., 2001; Sokasian et al., 2001; McQuinn et al., 2007; Iliev et al., 2006, 2007). Any hydrodynamical response to the radiation field was thus ignored. Several radiative transfer codes were compared in four purely radiative transfer tests in Iliev et al. (2006, hereafter RT06). Only recently has the radiative transfer equation been coupled to the hydrodynamics in three-dimensions (e.g. Krumholz et al., 2007). In the second comparison paper (Iliev et al., 2009, hereafter RT09), results from these radiation hydrodynamics codes were compared. Even more rare are ones that couple it with magneto-hydrodynamics (e.g. Krumholz et al., 2007). The tests in RT06 and RT09 were kept relatively simple to ease the comparison.
In this paper, we present our implementation, enzo+moray, of adaptive ray tracing (Abel & Wandelt, 2002) in the cosmological hydrodynamics adaptive mesh refinement (AMR) code, enzo (Bryan & Norman, 1997; O’Shea et al., 2004). The radiation field is coupled to the hydrodynamics solver at small time-scales, enabling it to study radiation hydrodynamical problems. We have used this code to investigate the growth of an H ii region from a 100 Population III (Pop III) star (Abel et al., 2007), the early stages of reionisation from Pop III stars (Wise & Abel, 2008a), radiative feedback on the formation of high redshift dwarf galaxies (Wise & Abel, 2008b; Wise et al., 2010), ultraviolet radiation escape fractions from dwarf galaxies before reionisation (Wise & Cen, 2009), negative radiative feedback from accreting Pop III seed black holes (Alvarez et al., 2009), and radiative feedback in accreting supermassive black holes (Kim et al., 2011, in prep.). We have included enzo+moray in the latest public release of enzo111http://enzo.googlecode.com, and it is also coupled with the newly added MHD solver in enzo (Wang & Abel, 2009).
We have structured this paper as follows. In Section 2, we describe the mathematical connections between adaptive ray tracing and the radiative transfer equation. Furthermore, we detail how physics other than photo-ionisation and photo-heating are included. We then derive a geometric correction factor to any ray tracing method to improve accuracy. We end the section by describing a new computational technique to approximate an optically-thin radiation field with ray tracing and multiple sources. In Section 3, we cover the details of our radiation hydrodynamics implementation in enzo, specifically (1) the ray tracing algorithms, (2) coupling with the hydrodynamics solver, (3) several methods to calculate the radiative transfer timestep, and (4) our parallelisation strategy. We present our results from the RT06 radiative transfer tests in Section 4. Afterwards in Section 5, we show the results from the RT09 radiation hydrodynamics tests. In Section 6, we expand on these tests to include more dynamical and complex setups to demonstrate the flexibility and high fidelity of enzo+moray. Section 7 gives the results from spatial, angular, frequency, and temporal resolution tests. In Section 8, we illustrate the improvements from the geometric correction factor and our optically-thin approximation. We also show the effects of X-ray radiation and radiation pressure in this section. Finally in Section 9, we demonstrate the parallel scalability of enzo+moray. Last Section 10 summarises our method and results.
2 Treatment of Radiative Transfer
Radiation transport is a well-studied topic, and we begin by describing our approach in solving the radiative transfer equation, which in comoving coordinates (Gnedin & Ostriker, 1997) is
Here is the radiation specific intensity in units of energy per time per solid angle per unit area per frequency . is the Hubble constant, where is the scale factor. is the ratio of scale factors at the current time and time of emission. The second term represents the propagation of radiation, where the factor accounts for cosmic expansion. The third term describes both the cosmological redshift and dilution of radiation. On the right hand side, the first term considers the absorption coefficient , and the second term is the emission coefficient that includes any point sources of radiation or diffuse radiation. We neglect any terms in equation (2) that become important in the dynamic diffusion limit (), where is the characteristic size of the system. This occurs in relativistic flows or very optically thick systems, such as stellar interiors or radiation-dominated shocks (see Krumholz et al., 2007, for a rigorous derivation that includes terms to second-order).
Solving this equation is difficult because of its high dimensionality; however, we can make some appropriate approximations to reduce its complexity in order to include radiation transport in numerical calculations. Typically timesteps in dynamic calculations are small enough so that , therefore in any given timestep, reducing the second term to . To determine the importance of the third term, we evaluate the ratio of the third term to the second term. This is , where is the simulation box length. If this ratio is , we can ignore the third term. For example at , this ratio is 0.1 when = 53 proper Mpc. In large boxes where the light crossing time is comparable to the Hubble time, then it becomes important to consider cosmological redshifting and dilution of the radiation. Thus equation (2) reduces to the non-cosmological form in this local approximation,
We choose to represent the source term as point sources of radiation (e.g. stars, quasars) that emit radial rays that are propagated along the direction . Next we describe this discretisation and its contribution to the radiation field.
|Cell face area|
|Collisional ionization rate of hydrogen|
|CFL safety factor in timestep calculation|
|Distance from ray segment center to cell center|
|Distance from ray segment center to cell edge|
|Photon loss from absorption|
|Photon loss from Compton Scattering|
|Momentum change from radiation pressure|
|Ionization energy of absorber|
|Geometric correction factor|
|Shielding function for H|
|Photo-dissociation rate of H|
|Linear width of a HEALPix pixel|
|Column density of H|
|Number density of absorber|
|HEALPix pixels on level|
|Rays per cell|
|Normal direction of radiation|
|Number density of absorber|
|Distance from the radiation source to the|
|Distance from the radiation source to the|
|next cell boundary crossing|
|Ionization front velocity|
|Cell center coordinates|
|Next cell boundary crossing in the|
|Source position in the -th dimension|
|Secondary ionization factors|
|Case-B recombination rate|
|Angular resolution in units of rays per cell area|
|Mean free path|
|Solid angle associated with a ray|
|Cross-section of absorber|
|Angle associated with a ray|
2.1 Adaptive Ray Tracing
Ray tracing is an accurate method to propagate radiation from point sources on a computational grid, given that there are a sufficient number of rays passing through each cell. Along a ray, the radiative transfer equation reads
where is the photon number flux along the ray. To sample the radiation field at large radii, ray tracing requires at least rays to sample each cell with one ray, where is the radius from the source to the cell and is the cell width. If one were to trace rays out to , the cells at a smaller radius would be sampled by, on average, rays, which is computationally wasteful because only a few rays per cell, as we will show later, provides an accurate calculation of the radiation field.
We avoid this inefficiency by utilising adaptive ray tracing (Abel & Wandelt, 2002), which progressively splits rays when the sampling becomes too coarse and is based on Hierarchical Equal Area isoLatitude Pixelation (HEALPix; Górski et al., 2005). In this scheme, the rays are traced along normal directions of the centres of HEALPix pixels, which evenly divides a sphere into equal areas. The rays are initialised at each point source with the photon luminosity (ph s) equally spread across rays, where is the initial HEALPix level. We usually find = 0 or 1 is sufficient because these coarse rays will usually be split before traversing the first cell.
The rays are traced through the grid in a typical fashion (e.g. Abel et al., 1999), in which we calculate the next cell boundary crossing. The ray segment length crossing the cell is
where , , , and are the initial distance travelled by the ray, normal direction of the ray, the next cell boundary crossing in the -th dimension, and the position of the point source that emitted the ray, respectively. However before the ray travels across the cell, we evaluate the ratio of the face area of the current cell and the solid angle of the ray,
If is less than a pre-determined value (usually ), the ray is split into 4 child rays. We investigate the variations in solutions with in §7.2. The pixel numbers of the child rays are given by the “nested” scheme of HEALPix at the next level, i.e. , where is the original pixel number. The child rays (1) acquire the new normal vectors of the pixels, (2) retain the same radius of the parent ray, and (3) gets a quarter of the photon flux of the parent ray. Afterwards the parent ray is discontinued.
A ray propagates and splits until
the photon has travelled , where is the radiative transfer timestep,
its photon flux is almost fully absorbed () in a single cell, which significantly reduces the computational time if the radiation volume filling fraction is small,
the photon leaves the computational domain with isolated boundary conditions, or
the photon travels of the simulation box length with periodic boundary conditions.
In the first case, the photon is halted at that position and saved, where it will be considered in the solution of at the next timestep. In the next timestep, the photon will encounter a different hydrodynamical and ionisation state, hence , in its path. Furthermore any time variations of the luminosities will be retained in the radiation field. This is how this method retains the time derivative of the radiative transfer equation. The last restriction prevents our method from considering sources external to the computational domain, but a uniform radiation background can be used in conjunction with ray tracing in enzo+moray that adds the local radiation field to the background intensity.
2.2 Radiation Field
The radiation field is calculated by integrating equation (4) along each ray, which is done by considering the discretisation of the ray into segments. In the following section, we assume the rays are monochromatic. For convenience, we express the integration in terms of optical depth , and for a ray segment,
Here and are the cross section and number density of the absorbing medium, respectively. We use the cell-centred density in our calculations. Using trilinearly interpolated densities (see Mellema et al., 2006) did not produce improved results. In the static case, equation (4) has a simple exponential analytic solution, and the photon flux of a ray is reduced by
as it crosses a cell. We equate the photo-ionisation rate to the absorption rate, resulting in photon conservation (Abel et al., 1999; Mellema et al., 2006). Thus the photo-ionisation and photo-heating rates associated with a single ray are
where is the cell volume, is the photon energy, and is the ionisation energy of the absorbing material. In each cell, the photo-ionisation and photo-heating rates from each ray in the calculation are summed, and after the ray tracing is complete, these rates can be used to update the ionisation state and energy of the cells. Considering a system with only hydrogen photo-ionisations and radiative recombinations, these changes are very straightforward and is useful for illustrative purposes. The change in neutral hydrogen is
where is the recombination coefficient at 10 K in the Case B on-the-spot approximation in which all recombinations are locally reabsorbed, (Spitzer, 1978), and is the collisional ionisation rate. However, for more accurate solutions in calculations that consider several chemical species, the photo-ionisation rates are terms in the relevant chemical networks (e.g. Abel et al., 1997).
2.3 Geometric Corrections
For a ray tracing method to accurately, i.e. without non-spherical artifacts, compute the radiation field, the computational grid must be well-sampled by the rays. The main source of potential artifacts is the geometrical difference between the cell and the HEALPix pixel relevant in the angular integration of the intensity over the cell. In this section, we devise a correction scheme to account for these differences. Consider the solid angle and photon flux associated with a single ray, and assume the flux is constant across . There exists a discrepancy between the geometry cell face and HEALPix pixel when the pixel does not cover the entire cell face, which is illustrated in Fig. 1. This mismatch causes non-spherical artifacts and is most apparent in the optically thin case, where the area of the pixel is dominant over when calculating . One can avoid these artifacts by increasing the sampling to high values () but we have formulated a simple geometric correction to the calculation of the radiation field. This correction is not unique to the HEALPix formalism but can be applied to any type of pixelisation.
The contribution to and must be corrected by the covering factor of the pixel with respect to the cell. When the pixel is fully contained within the cell face, . Because the geometry of the pixel can be complex with curved edges, we approximate by assuming the pixel is square. The covering factor is thus related to the width of a pixel, , and the distance from the ray segment midpoint to the closest cell boundary , which is depicted in Fig. 1. To estimate , we first find the distance from the midpoint of the ray segment to the cell centre in orthogonal directions,
where is the distance travelled by the ray in each orthogonal direction. The distance to the closest cell boundary is . Thus the covering factor is related to the square of the ratio between the and by
One half of the pixel is always contained within the cell, which results in the factor of 1/2. Finally we multiply and by but leave the absorbed radiation untouched because this would underestimate the attenuation of the incoming radiation. Using calculated like above, the method is no longer photon conserving. In our implementation, we felt that the spherical symmetry obtained outweighed the loss of photon conservation. However we show that there are no perceptible deviations from photon conservation in §4.1 and §7.1.
We briefly next describe how to retain photon conservation with a geometric correction. Notice that we compute by only considering the distances in orthogonal directions. A better estimate would consider the distance between the cell boundary and ray segment midpoint in the direction between the midpoint and cell centre of . We find that the method outlined here provides a sufficient correction factor to avoid any non-spherical artifacts and deviations from photon conservation. Furthermore in principle, the ray should also contribute to any neighboring cells that overlap with , which is the key to be photon conservative with such a geometric correction.
2.4 Optically Thin Approximation
In an optically thin medium, radiation is only attenuated by geometric dilution in the local approximation to Equation (2), i.e. the inverse square law. With such a simple solution, the tracing of rays is wasteful, however these rays must be propagated because the the medium farther away can be optically thick. Here we describe a method that minimizes the computational work of ray tracing in the optically thin regime by exploiting this fact. Each ray tracks the total column density and the equivalent total optical depth traversed by the photon. If after the ray exits the cell, we calculate the photo-ionisation and photo-heating rates directly from the incoming ray instead of the luminosity of the source.
Note that the photon number in the ray has already been geometrically diluted by ray splitting. Here and are the radii from the radiation source to the cell centre and where the ray exits the cell. Thus the last factor corrects the flux to a value appropriate for the cell centre. The photo-heating ray is calculated in the same manner as the general case, . This should only be evaluated once per cell per radiation source. No photons are removed from the ray. With this method, we only require one ray travel through each cell where the gas is optically thin, thus reducing the computational expense.
We must be careful not the overestimate the radiation when multiple rays enter a single cell. In the case of a single radiation source, the solution is simple – only assign the cell the photo-ionisation and photo-heating rates when . However in the case with multiple sources, this is no longer valid, and we must sum the flux from all optically thin sources. Only one ray per source must contribute to a single cell in this framework. We create a flagging field that marks whether a cell has already been touched by an optically thin photon from a particular radiation source. Naively, we would be restricted to tracing rays from a single source at a time if we use a boolean flagging field. However we can trace rays for 32 sources at a time by using bitwise operations on a 32-bit integer field. For example in C, we would check if an optically thin ray from the n-th source has propagated through cell i by evaluating (MarkerField[i] n & 1). If false, then we can add the optically thin approximation [equation (14)] to the cell and set MarkerField[i] |= (1 n); to mark the cell.
2.5 Additional Physics
Other radiative processes can also be important in some situations, such as attenuation of radiation in the Lyman-Werner bands, secondary ionisations from X-ray radiation, Compton heating of from scattered photons, and radiation pressure. We describe our implementation of these physics next.
2.5.1 Absorption of Lyman-Werner Radiation
Molecular hydrogen can absorb photons in the Lyman-Werner bands through the two-step Solomon process, which for the lowest ro-vibrational states already consists of 76 absorption lines ranging from 11.1 to 13.6 eV (Stecher & Williams, 1967; Dalgarno & Stephens, 1970; Haiman et al., 2000). Each of these spectral lines can be modelled with a typical exponential attenuation equation (Ricotti et al., 2001), but Draine & Bertoldi (1996) showed that this self-shielding is well modeled with the following relation to total H column density
where is in units of cm. To incorporate this shielding function into the ray tracer, we store the total H column density and calculate the H dissociation rate by summing the contribution of all rays,
where is the effective cross-section of H (Abel et al., 1997). To account for absorption, we attenuate the photon number flux by
where is the H column density in the current cell.
2.5.2 Secondary Ionisations from X-rays
At the other end of the spectrum, a high-energy ( eV) photon can ionise multiple neutral hydrogen and helium atoms, and this should be considered in such radiation fields. Shull & van Steenberg (1985) studied this effect with Monte Carlo calculations over varying electron fractions and photon energies up to 3 keV. They find that the excitation of hydrogen and helium and the ionisation of He ii is negligible. The number of secondary ionisations of H and He is reduced from the ratio of the photon and ionisation energies () by a factor of
where is the electron fraction. The remainder of the photon energy is deposited into thermal energy that is approximated by
and approaches one as . Thus in gas with low electron fractions, most of the energy results in ionisations of hydrogen and helium, and in nearly ionised gas, the energy goes into photo-heating.
2.5.3 Compton Heating from Photon Scattering
High energy photons can also cause Compton heating by scattering off free electrons. During a scattering, a photon loses of energy, where is the electron temperature. For the case of monochromatic energy groups, we model this process by considering that the photons are absorbed by a factor of
which is the equivalent of the photon energy decreasing. Here is the optical depth to Compton scattering, and is the non-relativistic Klein-Nishina cross section (Rybicki & Lightman, 1979). The Compton heating rate is thus
which has been used in Kim et al. (2011).
2.5.4 Radiation Pressure
Another relevant process is radiation pressure, where the absorption of radiation transfers momentum from photons to the absorbing medium. This is easily computed by considering the momentum
of the absorbed radiation from the incoming ray, where is the normal direction of the ray. We do not include radiation pressure on dust currently. The resulting acceleration of the gas because of radiation pressure is
where is the gas density inside the cell. This acceleration is then added to the other forces, e.g. gravity, in the calculation in an operator split fashion.
3 Numerical Implementation in Enzo
In this section, we describe our parallel implementation of the adaptive ray tracing method into enzo. enzo is a parallel block-structured AMR (Berger & Colella, 1989) code that is publicly available (Bryan & Norman, 1997; O’Shea et al., 2004). First we explain the programming design of handling the “photon packages” that are traced along the adaptive rays. We use the terms photon packages and rays interchangeably. Next we focus on the details of the radiation hydrodynamics and then the importance of correct time-stepping. Last we give our parallelisation strategy of tracing rays through an AMR hierarchy. This implementation is included in the v2.0 public version of enzo.
3.1 Programming Design
Each photon package is stored in the AMR grid with the finest resolution that contains its current position. The photon packages keep track of their (1) photon flux, (2) photon type, (3) photon energy, (4) the time interval of its emission, (5) emission time, (6) current time, (7) radius, (8) total column density, (9) HEALPix pixel number, (10) HEALPix level, and (11) position of the originating source, totaling 60 (88) bytes for single (double) precision. When enzo uses double precision for grid and particle positions and time, items 4-7 and 11 are double precision.
We only treat point sources of radiation in our implementation; therefore all base level photon packages originate from them. As they travel away from the source, they generally pass through many AMR grids, especially if the simulation has a high dynamic range. This is a challenging programming task as rays are constantly entering and exiting grids. Before any computation, the number of rays in a particular grid is highly unpredictable because the intervening medium is unknown. Furthermore, the splitting of parent rays into child rays and a dynamic AMR hierarchy add to the complexity. Because of this, we store the photon packages as a doubly linked list (Abel & Wandelt, 2002). Thus we can freely add and remove them from grids without the concern of allocating enough memory before the tracing commences.
We illustrate the underlying algorithm of the ray tracing module in enzo in Fig. 2 and the ray tracing algorithm is shown in Fig. 3. The module is only called when advancing the finest AMR level. We describe its steps below.
Step 1.– Create new photon packages on the initial HEALPix level from point sources. Place the new rays in the highest resolution AMR grid that contains the source.
Step 2.– Initialise all radiation fields to zero.
Step 3.– Loop through all AMR grids, tracing any rays that exist in it. For each ray, the following substeps are taken.
Step 3a.– Compute the ray normal based on the HEALPix level and pixel number of the photon package with the HEALpix routine pix2vec_nest. One strategy to accelerate the computation is to store ray segment paths in memory (Abel & Wandelt, 2002; Krumholz et al., 2007); however this must be recomputed if the grid structure or point source position changes. We do not restrict these two aspects and cannot employ this acceleration method.
Step 3b.– Compute the position of the ray (), the current cell coordinates in floating point and its corresponding integer indices. Here is the position of the point source, is the distance travelled by the ray, and is the ray normal.
Step 3c.– Check if a subgrid exists under the current cell. If so, move the ray to a linked list that contains all rays that should be moved to other grids. We call this variable PhotonMoveList. Store the destination grid number and level. Continue to the next ray in the grid (step 3a). We determine whether a subgrid exists by creating a temporary 3D field of pointers that either equals the pointer of the current grid if no subgrid exists under the current cell or the child grid pointer that exists under the current cell. This provides a significant speedup when compared to a simple comparison of a photon package position and all of the child grid boundaries. Note that this is the same algorithm used in enzo when moving collisionless particles to child grids.
Step 3d.– Compute the next cell crossing of the ray and the ray segment length across the cell (Equation 5).
Step 3e.– Compare the solid angle associated with the ray at radius with a user-defined splitting criterion (Equation 6). If the solid angle is larger than the desired minimum sampling, split the ray into 4 child rays (§2.1). These rays are inserted into the linked list after the parent ray, which is subsequently deleted. Continue to the next ray (step 3a), which will be the first child ray.
Step 3f.– Calculate the geometric correction (Equation 13), the optical depth of the current cell (Equation 7), photo-ionisation and photo-heating rates (Equations 9 and 10), and add the column density of the cell to the total column density of the ray.
Step 3g.– Add the effects of any optional physics modules (§2.5) – secondary ionisations from X-rays, Compton heating from scattering, and radiation pressure.
Step 3h.– Update the current time (), photon flux (, Equation 8), and radius of the ray ().
Step 3i.– If the photon flux is zero or the total optical depth is large (), delete the ray.
Step 3j.– Check if the ray has travelled a total distance of in the last timestep. If we are keeping the time-derivative of the radiative transfer equation, halt the photon. If not (i.e. infinite speed of light), delete the photon.
Step 3k.– Check if the ray has exited the current grid. If false, continue to the next cell (step 3b). If true, move the ray to the linked list PhotonMoveList, similar to step 3c. If the ray exits the simulation domain, delete it if the boundary conditions are isolated; otherwise, we change the source position of the ray by a distance -sign(n[i]) * DomainWidth[i], where n is the ray normal, and i is the dimension of the outer boundary it has crossed. The radius is kept unchanged. In essence, this creates a “virtual source” outside the box because the ray will be moved to the opposite side of the domain, appearing that it originated from this virtual source.
Step 4.– If any rays exist in the linked list PhotonMoveList, move them to their destination grids and return to step 3. This requires MPI communication if the destination grid exists on another processor.
Step 5.– If all rays have not been halted (keeping the time-derivative of the radiative transfer equation), absorbed, or exited the domain, return to step 3.
Step 6.– With the radiation fields updated, call the chemistry and energy solver and update only the cells with radiation, which is discussed further in §3.3.
Step 7.– Advance the time associated with the photons by the global timestep (for its calculation, see §3.4). If does not exceed the time on the finest AMR level, return to step 1.
3.2 Energy groups
In our implementation, photon packages are mono-chromatic, i.e. energy groups (Mihalas & Mihalas, 1984, Ch. 6), and are assigned a photon type that corresponds whether it is a photon that (1) ionises hydrogen, (2) singly ionises helium, (3) doubly ionises helium, (4) has an X-ray energy, or (5) dissociates molecular hydrogen (Lyman-Werner radiation). One disadvantage of mono-chromatic rays is the number of rays increase with the number of frequency bins. However this allows for early termination of rays that are fully absorbed, which are likely to have high absorption cross-sections (e.g. H i ionisations near 13.6 eV) or a low initial intensity (e.g. He ii ionising photons in typical stellar populations). The other approach used by some groups (e.g. Trac & Cen, 2007) is to store all energy groups in a single ray. This reduces the number of the rays generated and the computation associated with the ray tracing. Unless the ray dynamically adjusts its memory allocation for the energy groups as they become depleted, this method is also memory intensive in the situation where most of the energy groups are completely absorbed but a few groups still have significant flux.
In practice, we have found that one energy group per photon type is sufficient to match expected analytical tests (§7.3). For example when modeling Population III stellar radiation (e.g. Abel et al., 2007; Wise & Abel, 2008b, for hydrogen ionising radiation only), we have 3 energy groups – H i, He i, He ii – each with an energy that equals the average photon energy above the ionisation threshold.
3.3 Coupling with Hydrodynamics
Solving the radiative transfer equation is already an intensive task, but coupling the effects of radiation to the gas dynamics is even more difficult because the radiation fields must be updated on a time-scale such that it can react to the radiative heating, i.e. sound-crossing time. The frequency of its evaluation will be discussed in the next section.
enzo solves the physical equations in an operator-split fashion over a loop of AMR grids. On the finest AMR level, we call our radiation transport solver before this main grid loop in the following sequence:
Solve for the radiation field with the adaptive ray tracer
Update species fractions and energies for cells with radiation with a non-equilibrium chemistry solver on subcycles (Equation 25).
For each grid:
Solve for the gravitational potential with the particle mesh method
Update species fractions and energies for cells without radiation with a non-equilibrium chemistry solver on subcycles (Equation 25).
Update particle positions
Star particle formation
All grids: Update solution from children grids
Since the solver must be called many times, the efficiency of the radiation solver is paramount. After every radiation timestep, we call the non-equilibrium chemistry and energy solver in enzo. This solves both the energy equation and the network of stiff chemical equations on small timesteps, i.e. subcycles (Anninos et al., 1997). The timestep is
where is the electron number density, is the specific energy, and is the hydrodynamic timestep. This limits the subcycle timestep to a 10% change in either electron density, neutral hydrogen density, or specific energy. In simulations without radiation, enzo calls this solver in a operation-split manner after the hydrodynamics module for grids only on the AMR level that is being solved. In simulations with radiative transfer, the radiation field can change on much faster time-scales than the normal hydrodynamical timesteps.
For example, a grid on level might have no radiation in its initial evaluation, but the ionisation front exists just outside its boundary. Then radiation permeates the grid in the time between , and the energy and chemical state of the cells must be updated with each radiation update to advance the ionisation front accurately. If one does not update these cells, it will appear that the ionisation front does not enter the grid until the next hydrodynamical timestep! Visually this appears as discontinuities in the temperature and electron fraction on grid boundaries. One may avert this problem by solving the chemistry and energy equations for every cell on every radiative transfer timestep, but this is very time consuming and unnecessary, especially if the radiation filling factor is small.
We choose to dynamically split the problem by cells with and without radiation. In every radiation timestep, the chemo-thermal state of only the cells with radiation are updated. For the solver subcycling, we replace with in Equation 25 in this case. Once the radiative transfer solver is finished with its timesteps, the hydrodynamic module is called, and then the chemo-thermal state of the cells without radiation are updated on a subcycle timestep stated in Equation 25.
For cells that transition from zero to non-zero photo-ionisation rates, the initial state that enters into the chemistry and energy solver does not correspond to the current time of the radiation transport solver , but either time if the grid level is the finest level because its chemo-thermal state has not been updated or time on all other levels. In principle, one could first revert the cell back to time and then update to with the chemistry and energy solver if the cell is on the finest level. However in practice, the time-scales in gas without radiation are small compared to the ionisation and heating time-scales when radiation is introduced. Therefore, we do not perform this correction and find that this does not introduce any inaccuracies in both test problems (see §4) and real world applications.
3.4 Temporal evolution
There have been several methods of choosing a maximal timestep to solve radiation transfer equation while retaining stability and accuracy. We describe several methods to calculate the radiative transfer timestep in this section. With a small enough timestep, the solution is accurate (ignoring any systematic ones), but the solver is slow. Furthermore for very small timesteps, the photon packages only advance a short distance, and they will exist in every cells with radiation and are stored between timesteps, excessively consuming memory. On the shortest time-scale, one can safely set the timestep to the light-crossing time of a cell (Abel et al., 1999; Trac & Cen, 2007) but encounters the problems stated above.
If the timestep is too large, the solution will become inaccurate; specifically, ionisation fronts will advance too slowly, as radiation intensity exponentially drops with a scale length of the mean free path
past the ionisation front. For example in our implementation, the chemo-thermal state of the system remains constant as the rays are traced through the cells. In the case of a single H ii region, the speed of the ionisation front is limited to approximately .
3.4.1 Minimizing neutral fraction change
Another strategy is restricting the neutral fraction to change a small amount, i.e. for a single cell,
where is the maximum fraction change in neutral fraction. Shapiro et al. (2004) found that this limited the speed of the ionisation front. We can investigate this further by evaluating the ionisation front velocity in a growing Strömgren sphere without recombinations,
Using and , we can make substitutions respectively on both sides of Equation (28) to arrive at the ionisation front velocity .
We have implemented this method but we only consider cells within the ionisation front (by experiment we choose ) because we are interested in evolving ionisation fronts at the correct speed. In the ionised region, the absolute changes in neutral fraction are small and will not significantly affect the ionisation front evolution. In other words, may be large but , thus we can safely ignore these cells when determining the timestep without sacrificing accuracy.
We search for the cell with the smallest based on Equation 27. In principle, one could use this value without modifications as the timestep, but there is considerable noise both spatially and temporally. In order to stabilise this technique, we first spatially smooth the values of with a Gaussian filter over a cube. Because we only consider the cells within the ionisation front, we set to the hydrodynamical timestep outside the front during the smoothing. After we have smoothed , we select the minimum value as . Significant noise in can exist between time steps as well. Because the solution can become inaccurate if the timestep is allowed to be too large, we restrict the next timestep to be less than twice the previous timestep,
Fig. 4 shows the smooth evolution of in a growing Strömgren sphere when compared to the values of .
3.4.2 Time averaged quantities within a timestep
Mellema et al. (2006) devised an iterative scheme that allows for large timesteps while retaining accuracy by considering the time-averaged values (, , , ) during the timestep. Starting with the cells closest to the source, they first calculate the column density to the cell. Then they compute the time-averaged neutral density for the cell and its associated optical depth, which is added to the total time-averaged optical depth. With these quantities, one can compute a photo-ionisation rate and update the electron density. This process is repeated until convergence is found in the neutral number density. In a test with a Strömgren sphere, they found analytical agreement with less timesteps than a method without time-averaging. Another advantage of this method is the use of pre-calculated tables for the photo-ionisation rates as a function of optical depth, based on a given spectrum. This minimizes the energy groups needed to accurately sample a spectrum. We are currently implementing this method into enzo+moray.
3.4.3 Physically motivated
A constant timestep is necessary when solving the time-dependent radiative transfer equation in enzo+moray. It should be small enough to evolve ionisation fronts accurately, as discussed earlier. The timestep can be based on physical arguments, for example, the sound-crossing time of an ionised region at K. To be conservative, one may choose the sound-crossing time of a cell (e.g. Abel et al., 2007; Wise & Abel, 2008b). Alternatively, the diameter of the smallest relevant system (e.g., an accretion radius, transition radius to a D-type ionisation front, etc.) in the simulation may be chosen to calculate the sound-crossing time.
If the timestep is too large, the ionisation front will propagate too slowly, but it eventually approaches the correct radius at late times (see §7.4). This does not prevent one from using a large timestep, particularly if the system is not critically affected by a slower I-front velocity. One example is an expanding H ii region in a power-law density gradient. After a brief, initial R-type phase, the I-front becomes D-type phase, where the ionisation and shock front progress jointly at the sound speed of the ionised region. A moderately large (0.1 Myr) timestep can accurately follow its evolution. However after the I-front passes a critical radius (Franco et al., 1990), the I-front detaches from the shock front, accelerates, and transitions back to an R-type front. This can also occur in champagne flows when the ionisation front passes a density discontinuity. The I-front velocities in these two stages differ up to a factor of . Although the solution is accurate with a large timestep in the D-type phase, the I-front may lag behind because of the constant timestep. After a few recombination times, the numerical solution eventually approaches the analytical solution. If such a simulation focuses on the density core expansion and any small scale structures, such as cometary structures and photo-dissociation regions, one can cautiously sacrifice temporal accuracy at large scales for computational savings.
3.4.4 Change of incident radiation
Ionisation front velocities can approach significant fractions of the speed of light in steep density gradients and in the early expansion of the H ii region. If the ionisation front position is critical to the calculation, the radiation transport timestep can be derived from a non-relativistic estimate of the ionisation front velocity
based on the incident radiation field at a particular position222See Shapiro et al. (2006) for the exact calculation of a relativistic ionisation front. Neglecting relativistic terms do not affect the solution because the front velocity is only considered in the timestep calculation.. Alternatively, the propagation of the ionisation front can be restricted by limiting the change in specific intensity to a safety factor , resulting in a timestep of
We consider the specific intensity after the ray travels through the cell, so , where is the optical depth through the cell. The time derivative of is
which can be expressed in terms of local optical depth and neutral fraction,
In practice, we have found that a ceiling of 3 can be placed on the optical depth, so optically thick cells do not create a very small timestep. We still find excellent agreement with analytical solutions with this approximation. We show the accuracy using this timestep method in §7.4.
3.5 Parallelisation Strategy
Parallelisation of the ray tracing code is essential when exploring problems that require high resolution and thus large memory requirements. Furthermore, enzo is already parallelised and scalable to processors in AMR simulations, and in unigrid calculations. enzo stores the AMR grid structure on every processor, but only one processor contains the actual grid and particle data and photon packages. All other processors contain an empty grid container. As discussed in Step 4 in §3.1, we store the photon packages that need to be transferred to other grids in the linked list PhotonMoveList. In a single processor (serial) run, moving the rays is trivial by inserting these photons into the linked list of the destination grid. For multi-processor runs, we must send these photons through MPI communication to the processors that host the data of the destination grids. We describe our strategy below.
The easiest case is when the destination grid exists on the same processor as the source grid, where we move the ray as in the serial case. For all other rays, we organise the rays by destination processors and send them in groups. We also send the destination grid level and ID number along with the ray information that is listed at the beginning of §3.1.
For maximum overlap of communication and computation, which enables scaling to large numbers of processors, we must employ “non-blocking” MPI communication, where each processor does not wait for synchronisation with other processors. We use this technique for the sending and receiving of rays. Here we desire to minimize the idle time of each processor when it is waiting to receive data. In the loop shown in Fig. 2 with the conditional that checks whether we have traced all of the rays, we aggressively transport rays that are local on the processor, and process any MPI receive calls as they arrive, not waiting for their completion in order to continue to the next iteration. We describe the steps in this algorithm next.
Step 1.– Before any communication occurs, we count the number of rays that will be sent to each processor. The MPI receive calls (MPI_Irecv) must have a data buffer that is greater than or equal to the size of the message. We choose to send a maximum of (= 10 in enzo v2.0) rays per MPI message. Therefore, we allocate a buffer of this size for each MPI_Irecv call. We then determine the number of MPI messages and send this number in a non-blocking fashion, i.e. MPI_Isend.
Step 2.– Pack the photon packages into a contiguous array for MPI communication while the messages from Step 1 completes.
Step 3.– Process the number of photon messages that we are expecting from each processor, sent in Step 1. Then post this number of MPI_Irecv calls for the photon data. Because we strive to make the ray tracing routine to be totally non-blocking, the processors will most likely not be synchronised on the same loop (Steps 3–5 in §3.1). Therefore, there might be additional MPI messages waiting to be processed. We check for these messages and aggressively drain the message stack to determine the total number of photon messages that we are expecting and post their associated MPI_Irecv calls for the photon data.
Step 4.– Send the grouped photon data with MPI_Isend with a maximum size of photons.
Step 5.– Place any received photon data into the destination grids. We monitor whether the processor has any rays that were moved to grids on the same processor. If so, this processor has rays to transport, and we do not necessarily have to wait for any MPI receive messages and thus use MPI_Testsome to receive any messages that have already arrived. If not, we call MPI_Waitsome to wait for any MPI receive messages.
Step 6.– If all processors have exhausted their workload, then all rays have been either absorbed, exited the domain, or halted after travelling a distance . We check this in a similar non-blocking manner as the calls in Step 1.
Lastly we have experimented with a hybrid OpenMP/MPI version of enzo, where workload is partitioned over grids on each MPI process. We found that parallelisation over grids for the photon transport does not scale well, and threading over the rays in each grid is a better approach. Because the rays are stored in a linked list in each grid, we must manually split the list into separate lists and let each thread work on each list.
4 Radiative Transfer Tests
Tests plays an important role in creating and maintaining computational tools. In this section, we present tests drawn from the Cosmological Radiative Transfer Codes Comparison Project (Iliev et al., 2006), where results from 11 different radiative transfer codes compared results in four test problems. The codes use various methods for radiation transport: ray tracing with short, long, and hybrid characteristics, Monte Carlo casting; ionisation front tracking (Alvarez et al., 2006b); variable Eddington Tensor formalism (Gnedin & Abel, 2001). They conducted tests that investigated (1) the growth of a single Strömgren sphere enforcing isothermal conditions, (2) the same test with an evolving temperature field, (3) shadowing created by a dense, optically thick clump, and (4) multiple H ii regions in a cosmological density field. In all of the tests presented here, we use the method of restricted neutral fraction changes (§3.4.1) for choosing a radiative transfer timestep. We cast 48 rays (HEALPix level 1) from the point source and require a sampling of at least rays per cell.
4.1 Test 1. Pure hydrogen isothermal H ii region expansion
The expansion of an ionising region with a central source in a uniform medium is a classic problem first studied by Strömgren (1939). This simple but useful test can uncover any asymmetries or artifacts that may arise from deficiencies in the method or newly introduced bugs in the development process. In this problem, the ionised region grows until recombinations balance photo-ionisations [equation (1)]. The evolution of the radius and velocity of the ionisation front has an exact solution of
where is the recombination time.
We adopt the problem parameters used in RT06. The ionising source emits ph s of monochromatic radiation at 13.6 eV and is located at the origin in a simulation box of 6.6 kpc. The ambient medium is initially set at K, , , resulting in kpc and Myr. The problem is run for 500 Myr. In the original tests, the temperature is fixed at K; however, our solver is inherently tied to the chemistry and energy solver. To mimic an isothermal behavior, we set the adiabatic index , which ensures an isothermal state but not a fixed ionisation fraction outside of the Strömgren sphere.
In Fig. 5, we show (a) the evolution of the neutral and ionisation fraction as a function of radius at 10, 30, 100, and 500 Myr, and (b) the growth of the ionisation front radius as a function of time and its ratio with the analytical Strömgren radius [equation (35)]. The ionisation front has a width of kpc, which is in agreement with the inherent thickness of kpc, given a 13.6 eV mono-chromatic spectrum. There are small kinks in the neutral fraction at 1.5 and 3 kpc that corresponds to artifacts created by the photon package splitting at these radii. However these do not affect the overall solution. One difference between our results and the codes presented in RT06 is the increasing neutral fraction outside of the H ii region. This occurs because the initial ionised fraction and temperature is set to and K, which are not the equilibrium values. Over the 500 Myr in the calculation, the neutral fraction increases to 0.2, which is close to its equilibrium value. In the right panel of Fig. 5, the ionisation front radius exceeds by a few percent for most of the calculation. This difference happens because the analytical solution [equation (35)] assumes the H ii region has a constant ionised fraction. The evolution of the ionised fraction as a function of radius can be analytically calculated (e.g. Osterbrock, 1989; Petkova & Springel, 2009), causing the ionisation front radius to be slightly larger, increasing from 0 to 3% in the interval 80–350 Myr. Our results are in excellent agreement with this more accurate analytical solution. We show the distribution of neutral fraction in the domain for = 10, 100, and 500 Myr in Fig. 6 that agrees well with the results in RT06. In Fig. 7, we show a slice of the neutral fraction through the origin. Other than the ray splitting artifacts that generate the plateaus at 1.5 and 3 kpc, one sees spherical symmetry without any noise in our solution.
4.2 Test 2. H ii region expansion: temperature evolution
This test is similar to Test 1, but the temperature is allowed to evolve. The radiation source now has a blackbody spectrum with a K. The initial temperature is set at 100 K. The higher energy photons have a longer mean free path than the photons at the ionisation threshold in Test 1. Thus the ionisation front is thicker as the photons can penetrate deeper into the neutral medium. Here we use 4 energy groups with the following mean energies and relative luminosities: , .
In Fig. 8, we show the radially averaged neutral and ionised fraction at 10, 30, 100, and 500 Myr, and the total neutral fraction of the domain in Fig. 9. Compared with Test 1, the ionisation front is thicker, as expected with the harder spectrum. Artifacts originating from ray splitting, similar to Test 1, appear at and 3 kpc as kinks in the neutral fraction. The total neutral fraction decreases to 0.67 over 4 Myr, which is in agreement with the analytical expectation and other codes in RT06. Fig. 10 shows the distribution of ionized fraction and temperature in Test 2. The overall trends agree with the codes presented in RT06 with the exception of the ray splitting artifacts that appear as slight rises in the distribution at and . In Fig. 11, we show slices of neutral fraction and temperature through the origin at 10 and 100 Myr. Here one sees the spherically symmetric H ii regions and a smooth temperature transition to the neutral ambient medium. In Fig. 12, we show the ratio of the ionisation front radius in our simulation and . Before , lags behind , initially by 10% and then increases to ; however afterwards, this ratio asymptotes to a solution that is 4% greater than . This behavior is approximately the median result in RT06, where this ratio varies between 1 and 1.1, and the early evolution of is under-predicted by almost all of the codes. If we use one energy group with the mean energy (29.6 eV) of a K blackbody, we find that = 1.08, which is representative of the codes in the upper range of RT06.
4.3 Test 3. I-front trapping in a dense clump and the formation of a shadow
The diffusivity and angular resolution of a radiative transport method can be tested with the trapping of an ionisation front by a dense, neutral clump. In this situation, the ionisation front will uniformly propagate until it reaches the clump surface. Then the radiation in the line of sight of the clump will be absorbed more than the ambient medium. If the clump is optically thick, a shadow will form behind the clump. The sharpness of the ionisation front at the shadow surface can be used to determine the diffusivity of the method. Furthermore the shadow surface should be aligned with the outermost neutral regions of the clump, which can visually assess the angular resolution of the method.
The problem for this test is contained in a 6.6 kpc box with an ambient medium of and K. The clump is in pressure equilibrium with and K. It has a radius of kpc and is centred at kpc. The ionized fraction of the entire domain is zero. In RT06, the test considered a plane parallel radiation field with a flux ph s cm originating from the plane. Our code can only consider point sources, so we use a single radiation source located in the centre of the boundary. The luminosity of ph s corresponds to the same flux at 5 kpc. The location where the ionisation front trapping can be calculated analytically (Shapiro et al., 2004), and with these parameters, it should halt at approximately the centre of the clump. We use the same four energy groups as in Test 2.
In Fig. 13, we show neutral fraction and temperature of a one-dimensional cut through the centre of the dense clump at at Myr. The ionisation front is halted at a little more than halfway through the clump, which is consistent with the analytical expectation. The hardness of the K blackbody spectrum allows the gas outside of the ionisation front to be photo-heated. Where the gas is ionised, the temperature is between 10,000 and 20,000 K, but decreases sharply with the ionised fraction. Fig. 14 depicts the average ionised fraction and temperature inside the dense clump, which both gradually increase as the ionisation front propagates through the overdensity. Our results are consistent with RT06. The distribution of ionized fraction and temperature within the clump are shown in Fig. 15 and agree well with the other codes in RT06. Finally we show slices of neutral fraction and temperature in the plane in Fig. 16. The neutral fraction slices prominently show the sharp shadows created by the clump and demonstrates the non-diffusivity behavior of ray tracing. The discretisation of the sphere creates one neutral cell on the left side of the sphere. This inherent artifact to the initial setup carries through the calculation. We did not smooth the clump surface like in some of the RT06 codes, in order to remove this artifact. It is seen in the neutral fraction and temperature states at all times and is not a caused by our ray tracing algorithm.
4.4 Test 4. Multiple sources in a cosmological density field
The last test in RT06 involves a static cosmological density field at . The simulation comoving box size is 0.5 Mpc and has a resolution of 128. There are 16 point sources that are centred in the 16 most massive haloes. They emit ionising photons per baryon in a blackbody spectrum with an effective temperature K, and they live for Myr. Thus the luminosity of each source is
where is the halo mass, , and . The radiation boundaries are isolated so that the radiation leaves the box instead being shifted periodically. The simulation is evolved for 0.4 Myr.
Fig. 17 depicts the distribution of the neutral fraction and temperature of the entire domain and shows good agreement with the codes presented in RT06. We show the growth of the H ii regions by computing the mass-averaged and volume-averaged ionised fraction in Fig. 18. Initially is larger than , and at kyr, the becomes larger. This is indicative of inside-out reionisation (e.g. Gnedin, 2000; Miralda-Escudé et al., 2000; Sokasian et al., 2003), where the dense regions around haloes are ionised first, then the voids are ionised last. At the end of the simulation, 65% of the simulation is ionised. However by visual inspection in the slices of electron fraction (Fig. 19), there appears to be very good agreement with C-ray and FTTE. By first glance, our result appears to be different than the RT06 because of the color-mapping. Our results are also in good agreement with the multi-frequency version of TRAPHIC (Pawlik & Schaye, 2010, see also for better representations of the electron fraction slices). In the slices of electron fraction and temperature, Fig. 19, the photo-heated regions are larger than the ionised regions by a factor of 2–3 because of the hardness of the K blackbody spectrum.
5 Radiation Hydrodynamics Tests
We next show results from radiation hydrodynamics test problems presented in Iliev et al. (2009, hereafter RT09). They involve (1) the expansion of an H ii region in a uniform medium, similar to Test 2, (2) an H ii region in an isothermal sphere, and (3) the photo-evaporation of a dense, cold clump, similar to Test 3. We turn off self-gravity and AMR in accordance with RT09.
5.1 Test 5. Classical H ii region expansion
Here we consider the expansion of an H ii region into a uniform neutral medium including the hydrodynamical response to the heated gas. The ionised region has a greater pressure than the ambient medium, causing it to expand. This is a well-studied problem (Spitzer, 1978) with an analytical solution, where the ionisation front moves as
where is the sound speed of the ionised gas and is the in Equation 35. The bubble eventually reaches pressure equilibrium with the ambient medium at a radius
where and are the ionised and ambient temperatures, respectively. These solutions only describe the evolution at late times, and not the fast transition from R-type to D-type at early times.
The simulation setup is similar to Test 2 with the exception of the domain size kpc. Here pressure equilibrium occurs at kpc, which is not captured by this test. However more interestingly, the transition from R-type to D-type is captured and occurs around kpc. The test is run for 500 Myr.
The growth of the ionisation front radius is shown in Fig. 20, using both K and as ionisation front definitions, compared to the analytical solution [equation (39)]. We define this alternative measure because the ionisation front becomes broad as the D-type front creates a shock. Densities in this shock, as seen in Fig. 21, are high enough for the gas to recombine but not radiatively cool. Before Myr, the temperature cutoff overestimates by over 10%; however at later times, it provides a good match to the growth at late times. With the criterion for the ionisation front, the radius is always underestimated by . This behavior was also seen in RT09.
Fig. 21 shows the progression of the ionisation front at times = 10, 200, and 500 Myr in radial profiles of density, temperature, pressure, and ionised fraction. The initial H ii region is over-pressurised and creates an forward shock wave. The high-energy photons can penetrate through the shock and partially ionises and heats the exterior gas, clearly seen in the profiles. As noted in RT09, this heated exterior gas creates an photo-evaporative flow that flows inward. This interacts with the primary shock and creates the double-peaked features in the density profiles at 200 and 500 Myr. Fig. 22 shows slices through the origin of the same quantities, including neutral fraction. These depict the very good spherical symmetry of our method. The only apparent artifact is a very slight diagonal line, which is caused by the HEALPix pixelisation differences between the polar and equatorial regions. This artifact diminishes as the ray-to-cell sampling is increased.
5.2 Test 6. H ii region expansion in an isothermal sphere
A more physically motivated scenario is an isothermal sphere with a constant density core, which is applicable to collapsing molecular clouds and cosmological haloes. The radial density profile is described by
where is the radius of the core. If the Strömgren radius is smaller than the core radius, then the resulting H ii region never escapes into the steep density slope. When the ionisation front propagates out of the core, it accelerates as it travels down the density gradient. There exists no analytical solution for this problem with full gas dynamics but was extensively studied by Franco et al. (1990). After the gas is ionised and photo-heated, the density gradient provides the pressure imbalance to drive the gas outwards.
This test is constructed to study the transition from R-type to D-type in the core and back to R-type in the density gradient. Thus the simulation focuses on small-scale, not long term, behavior of the ionisation front. The simulation box has a side length kpc with core density , core radius pc (15 cells), and temperature K throughout the box. The ionisation fraction is initially zero, and the point source is located at the origin with a luminosity of ph s cm. The simulation is run for 75 Myr.
Because this problem does not have an analytical solution, we compare our calculated ionisation front radius and velocity, shown in Fig. 23, to the RT09 results. Their evolution are in agreement within 5% of RT09. As in Test 5, we use an extra definition of K for the ionisation front. We compute the ionisation front velocity from the radii at 50 outputs, which causes the noise seen Fig. 23.
For the first Myr, the radiation source creates a weak R-type front where the medium is heated and ionised but does not expand because . When becomes subsonic, the medium can react to the passing ionisation front and creates a shock, leaving behind a heated rarefied medium. This behavior is clearly seen in the radial profiles of density, temperature, ionised fraction, and pressure in Fig. 24. The inner density decreases over two order of magnitude after 25 Myr. To illustrate any deviations in spherical symmetry, we show in Fig. 25 slices of density, temperature, neutral fraction, and ionised fraction at the final time. The only artifact apparent to us is the slight broadening of the shock near the and planes. This causes the ionisation front radius to be slightly smaller in those directions. In the diagonal direction, the neutral column density through the shock is slightly smaller, allowing the high-energy photons to photo-ionise and photo-heat the gas to and K out to pc from the shock. The reflecting boundaries are responsible for this artifact because this is not seen when the problem is centred in the domain, removing any boundary effects.
5.3 Test 7. Photo-evaporation of a dense clump
The photo-evaporation of a dense clump in a uniform medium proceeds very differently when radiation hydrodynamics is considered instead of a static density field. The ionisation front first proceeds as a very fast R-type front, then it slows to a D-type front when it encounters the dense clump. As the clump is gradually photo-ionised and heated, it expands into the ambient medium. The test presented here is exactly like Test 3 but with gas dynamics. In this setup, the ionisation front overtakes the entire clump, which is then completely photo-evaporated.
Fig. 26 shows cuts of density, temperature, neutral fraction, and pressure in a line connecting the source and the clump centre at , 10, and 50 Myr. At 1 Myr, the ionisation front has propagated through the left-most 500 pc of the clump. This heated gas is now over-pressurised, as seen in the pressure plot in Fig. 26, and then expands into the ambient medium. This expansion creates a photo-evaporative flow, seen in many star forming regions (e.g. M16; Hester et al., 1996) as stars irradiate nearby cold, dense overdensities. These flows become evident in the density at 10 Myr, seen both in the line cuts and slices (Fig. 27). They have temperatures up to 50,000 K. At this time, the front has progressed about halfway through the clump, if one inspects the neutral fraction. However the high energy photons have heated all but the rear surface of the clump. At the end of the test ( Myr), only the core and its associated shadow is neutral, seen in Fig. 28. The core has been compressed by the surrounding warm medium, thus causing the higher densities seen at Myr. The non-spherical artifacts on the inner boundary of the warm outermost shell are caused by the initial discretisation of the sphere, as discussed in §4.3. Next Fig. 29 shows the distributions of temperature and flow Mach number in the entire domain at 10 Myr and 50 Myr, showing similar behavior as the codes in RT09. Lastly in Fig. 30, we show slices of the flow Mach number at 10 Myr and 50 Myr, showing the supersonic photo-evaporative flows that originate from the clump.
6 Radiation Hydrodynamics Applications
We have completed presenting results from the RT06 and RT09 test suites. We expand on these test suites to include more complex situations, such as a Rayleigh-Taylor problem illuminated by a radiation source, champagne flows, an irradiated blast wave, collimated radiation, and an H ii with a variable source that further demonstrate its capabilities and accuracy. Last we use the new MHD implementation in enzo v2.0 in the problem of a growing H ii region in a magnetic field.
6.1 Application 1. Champagne flow from a dense clump
Radiation-driven outflows from overdensities, known as champagne flows, is a long studied problem (e.g. Yorke, 1986, §3.3). To study this, we use the same setup as Bisbas et al. (2009) – a spherical tophat with an overdensity of 10 and radius of 1 pc in a simulation box of 8 pc. The ambient medium has cm and K. The radiation source is offset from the overdensity centre by 0.4 pc. It has a luminosity of ph s and a K blackbody spectrum. The resulting Strömgren radius is 0.33 pc, just inside of the overdense clump. These parameters are the same used in Bisbas et al. (2009). The entire domain initially has an ionised fraction of . We do not consider self-gravity. The simulation has a resolution of 128 on base grid, and we refine the grid up to 4 times if a cell has an overdensity of , where is the AMR level. The simulation is run for 150 kyr.
We show slices in the x-y and x-z planes of density in Fig. 31 at kyr. In the direction of the clump centre, the ionisation front shape transitions from spherical to parabolic after it escapes from the clump in the opposite direction. At kyr, the surface of the H ii region is just contained within the overdensity. In the x-z plane, there are density perturbations only above a latitude of 45 degrees. We believe that these are caused by the mismatch between HEALPix pixels and the Cartesian grid, even with our geometric correction. After the ionisation front escapes from the clump in the negative x-direction, these perturbations grow from Rayleigh-Taylor instabilities as the gas is accelerated when it exits the clump. As the shock propagates through the ambient medium, it is no longer accelerated and has a nearly constant velocity, as seen in Test 6. Thus these perturbations are not as vulnerable to Rayleigh-Taylor instabilities at this point. The ambient medium and shock are always optically thick, even in the directions of the bubbles. Bisbas et al. found that the shock fragmented and formed globules; however we find the density shell is stable against such fragmentation. To investigate this scenario further, our next tests involve radiation driven Rayleigh-Taylor instabilities.
6.2 Application 2. Irradiated Rayleigh-Taylor instability
Here we combine the classic case of a Rayleigh-Taylor instability and an expanding H ii region. The Rayleigh-Taylor instability occurs when a dense fluid is being supported by a lighter fluid, initially in hydrostatic equilibrium, in the presence of a constant acceleration field. This classic test alone evaluates how subsonic perturbations evolve. We consider the case of a single-mode perturbation. The system evolves without any radiation until the perturbation grows considerably and then turn on the radiation source. These tests demonstrate that enzo+moray can follow a highly dynamic system and resolve fine density structures.
We run two cases, an optically-thick and optically-thin case. In the former, we take the parameter choices from past literature (e.g. Liska & Wendroff, 2003; Stone et al., 2008) by setting the top and bottom halves of the domain to a density and , respectively. The velocity perturbation is set in the -direction by
We set the acceleration field and the adiabatic index . We use a domain size of with a resolution of (64, 64, 192). For hydrostatic equilibrium, we set with . In order to consider a radiation source with a ionising photon luminosity of ph s, we scale the domain to a physical size of (0.5, 0.5, 1.5) pc; time is in units of Myr; density is in units of , resulting in an initial temperature of K. The radiation source is placed at the centre of lower -boundary face and starts to shine at Myr.
The optically-thin case is set up similarly but with three changes: (1) a density contrast of 10, (2) a luminosity of ph s, and (3) the source is born at 6.5 Myr. The time units are decreased to 200 kyr so that K. Note that in code units, pressure is unchanged. We adjust the physical unit scaling because we desire an optically thin bottom medium with K and . Furthermore, the ionisation front remains R-type before interacting with the instability. A possible physical analogue could be a radiation source heating and rarefying the medium below.
The and -boundaries are periodic, and the -boundaries are reflecting. These will cause artificial features, in particular, because of the top reflecting boundary; nevertheless, these tests provide a stress test on a radiation hydrodynamics solver. We show the evolution of the density, temperature, and ionised fraction of the optically thick and optically thin cases in Figs. 32 and 33. The initial state of the Rayleigh-Taylor instability is shown in the left panels.
In the optically thick case, a D-type front is created, which is clearly illustrated by the spherical density enhancement at 0.02 Myr. The shock then passes through the instability at 0.25 Myr and reflects off the upper -boundary. This and complex shock reflections create a Richtmyer-Meshkov instability (see Brouillette, 2002, for a review), driving a chaotic jet-like structure downwards. The radiation source photo-evaporates the outer parts of this structure. The interaction between the dense cool “jet” and the hot medium further drives instabilities along the surface, which can be seen when comparing Myr and Myr slices. At the latter time, the jet cannot reach the bottom of the domain before being photo-evaporated. Eventually this structure is completely destroyed, leaving behind a turbulent medium between the hot and cold regions.
The optically thin problem is less violent than the optically thick case because the R-type front does not interact with the initial instability as strongly. The radiation source provides further buoyancy in the already K gas. The gas first to be ionised and photo-evaporated are the outer regions of the instability. The enhanced heating also drives the upper regions of the instability, making the top interface turbulent. It then reflects off the upper -boundary and creates a warm K, partially ionised (), turbulent medium, seen in the slices Myr. The slices of electron fraction also show that the dense gas is optically thick.
6.3 Application 3. Photo-evaporation of a blastwave
A supernova blastwave being irradiated by a nearby star is a likely occurrence in massive-star forming regions. In this test, we set up an idealised test that mimics this scenario. The ambient medium has a density and temperature K. The domain size is 1 kpc. We use 2 levels of AMR with a base grid of 64 that is refined if the density or total energy slope is greater than 0.4. The blastwave is initialised at the beginning of the Sedov-Taylor phase when the mass of the swept-up material equals the ejected material. It has a radius of 21.5 pc, a total energy of erg, and total mass of , corresponding to eV per particle or K. The radiation source is located at the centre of the left -boundary and has a luminosity of erg s. We use a K blackbody spectrum with 2 energy groups (16.0 and 22.8 eV). The source turns on at 2.5 Myr at which point the blastwave has a radius of 200 pc. The simulation is run for 7.5 Myr.
Fig. 34 shows the ionisation front overtaking and disrupting the blastwave. We show the blastwave before the source is born at 2.5 Myr. The interior is rarefied () and is heated to K by the reverse shock. At Myr, the ionisation front is still R-type, and it ionises the rear side of the dense shell. Because the interior is ionised and diffuse, the ionisation front rapidly propagates through it until it reaches the opposite shell surface. Shortly afterward, the ionisation front transitions from R-type to D-type at a radius of 0.5 kpc, seen in the formation of a shock in the 5 Myr density panel. This transition occurs by the construction of the problem not by the interaction with the blastwave. The surfaces of the blastwave that are perpendicular to the ionisation front have the highest column density and thus are last to be fully ionised. The pressure forces from warm ambient medium and blastwave interior compress these surfaces, photo-evaporating them in the process, similar to Test 7. They survive until the final time Myr. As the R-type ionisation front interacts with the blastwave interior, the density perturbations create ionisation front instabilities (Whalen & Norman, 2008) that are seen on the H ii region surface at the coordinate . Behind the ionisation front, the dense shell of the blastwave is photo-evaporated, and a smooth overdensity is left in the initial blastwave centre.
6.4 Application 4. Collimated radiation from a dense clump
Some astrophysical systems produce collimated radiation either intrinsically by relativistic beaming or by an optically-thick torus absorbing radiation in the equatorial plane. The latter case would be applicable in a subgrid model of active galactic nuclei (AGN) or protostars, for example. Simulating collimated radiation with ray tracing is trivially accomplished by only initialising rays that are within some opening angle .
We use a domain that is 2 kpc wide and has an ambient medium with cm, K, . We place a dense clump with , K, , and pc, at the centre of the box. Radiation is emitted in two polar cones with with 768 (HEALPix level 3) initial rays, a total luminosity of erg s, and a 17.6 eV mono-chromatic spectrum. This results in Myr and pc, just outside of the sphere. The base grid has a resolution of 64, and it is refined with the same overdensity criterion as Application 1. We run this test for 25 Myr.
We illustrate the expansion of the H ii region created by the beamed radiation in Fig. 35. Before Myr, the H ii region is conical and contained within the dense clump, depicted in the Myr snapshot of the system. At this time, the ionisation front is transitioning from R-type to D-type in the transverse direction of the cone. This can be seen in the minute overdensities on the H ii transverse surface. When it breaks out of the overdensity, a champagne flow develops, where the ionisation front transitions back to a weak R-type front. The cloud surface is a constant-pressure contact discontinuity (CD) with a density jump of 100. After the front heats the gas at the CD, there exists a pressure difference of . In response, the high density gas accelerates into the ambient medium and heats it to K. Additionally a rarefaction wave travels towards the clump centre. At later times, the transverse D-type front continues through the clump, eventually forming a disc-like structure at the final time. The polar champagne flows proceed to flow outwards and produces a dense shell with a diffuse ( cm) and warm (5000 K) medium in its wake.
6.5 Application 5. Time variations of the source luminosity
Our implementation retains the time derivative of the radiative transfer equation [equation (2)] if we choose a constant ray tracing timestep, which saves the photon packages between timesteps if . This effect only becomes apparent when the variation time-scale of the point source is smaller than the light crossing time of the simulation. Furthermore, the timestep should resolve the variation time-scale by at least a few times. This property might be important in large box simulations with variable sources, e.g. AGN radiative feedback. To test this, we can use an exponentially varying source with some duty cycle . In a functional form, this can be described as
where , and controls the width of the radiation pulse. To illustrate the effects of source variability, we remove any dependence on the medium by considering an optically-thin uniform density . We take Mpc, which has a light crossing time of 3.3 Myr. A source is placed at the origin with ph s and Myr. We use a radiative transfer timestep of 50 kyr to resolve the duty cycle by 10 timesteps. The simulation is run for 6 Myr, so the radiation propagates throughout the box.
The variability of the source is clearly illustrated in the photo-ionisation rates, shown in Fig. 36. The shells of relative maximum corresponds to radiation that was emitted when the source was at its peak luminosity. They are separated by Mpc and are geometrically diluted with increasing radius. Averaged over shells of the same width, photo-ionisation rates decrease as .
Next we test the hydrodynamical response to a varying source by repeating Test 5. We set the peak luminosity erg s that is a factor of 4 more luminous than Test 5, so the average luminosity is erg s. The spectrum is mono-chromatic with an energy of 29.6 eV. We set the variation time-scale kyr and use a constant radiative transfer timestep kyr. The simulation is run for 200 Myr. We show the radial profiles of density, temperature, ionised fraction, and pressure in Fig. 37. The variable source has little effect on the overall growth of the H ii region. It has the approximately the same radius as Test 5 at Myr when run with a mono-chromatic spectrum (see §7.3). At early times, the variable source creates density perturbations with an average size of 500 pc inside the ionised region, seen in the Myr profiles. They do not create any instabilities and are smoothed out over its sound crossing time of Myr.
6.6 Application 6. H ii region with MHD
Another prevalent physical component in astrophysics is a magnetic field. We utilise the new magnetohydrodynamics (MHD) framework (Wang & Abel, 2009) in enzo v2.0 that uses an unsplit conservative hydrodynamics solver and the hyperbolic cleaning method of Dedner et al. (2002). We show a test problem with an expanding H ii region in an initially uniform density field and constant magnetic field. We use the same problem setup as Krumholz et al. (2007) – , K, pc with a resolution of 256. This ambient medium is threaded by a magnetic field . The Alfvén speed is 2.6 km s. The radiation source is located in the centre of the box with a luminosity ph s with a 17.6 eV mono-chromatic spectrum, resulting in a Strömgren radius pc. The simulation is run for 1.58 Myr. The hydrodynamics solver uses an HLL Riemann solver (Harten et al., 1983) and piecewise linear method (PLM) reconstruction (van Leer, 1977) for the left and right states in this problem.
As the H ii region grows in the magnetised medium, shown in Figs. 38 and 39, it transforms from spherical to oblate as it is magnetically confined in directions perpendicular to the magnetic field. This occurs at Myr because the magnetic pressure exceeds the thermal pressure, and the gas can only flow along field lines. Krumholz et al. observed some carbuncle artifacts along the ionisation front; whereas we see smooth density gradients, which is most likely caused by both the geometric correction to the ray tracing (§2.3) and the diffusivity of the HLL Riemann solver when compared to Roe’s Riemann solver used in Krumholz et al. (2007), who also use PLM as a reconstruction method. The evolution of the magnetic field lines evolve in a similar manner as their results.
7 Resolution Tests
Resolution tests are important in validating the accuracy of the code in most circumstances, especially in production simulations where the initial environments surrounding radiation sources are unpredictable. In this section, we show how our adaptive ray-tracing implementation behaves when varying spatial, angular, frequency, and temporal resolutions.
7.1 Spatial resolution
Here we use Test 1 (§4.1) as a testbed to investigate how the evolution of the Strömgren radius changes with resolution. We keep all aspects of the test the same, but use resolutions of 16, 32, 64, and 128. In Fig. 40, we show the ratio , similar to Fig. 5, using these different resolutions. The radii in the and runs evolve almost identically. Compared to these resolutions, the lower and resolution runs only lag behind by 1% until 300 Myr, and afterwards it is larger by 0.5% than the higher resolution cases. This shows that our method gives accurate results, even in marginally resolved cases, which is expected with a photon conserving method. Furthermore this demonstrates that the geometric correction does not significantly affect photon conservation.
7.2 Angular resolution
The Cartesian grid must been sampled with sufficient rays in order to calculate a smooth radiation field. To determine the dependence on angular resolution, we consider the propagation of radiation through an optically thin, uniform medium. The radiation field should follow a profile. As the grid is less sampled by rays, the deviation from should increase. This test is similar to Test 1, but the medium has , K, and . The simulation is only run for one timestep because the radiation field should be static in this optically-thin test.
We consider minimum ray-to-cell ratios . Slices of the photo-ionisation rates through the origin are shown in Fig. 41 for these values of . In this figure, we limit the colormap range to a factor of 3 to show the nature of the artifacts in more contrast. Unscaled, the rates in the figures would span 4 orders of magnitude. When , the cell-to-cell variations are apparent because there are not enough rays to sufficiently sample the radiation field, even with the geometric correction factor , whose improvements are shown later in §8.1. At , these artifacts disappear, leaving behind a shell artifact where the radiation fields do not smoothly decrease as 1/. At higher values of , this shell artifact vanishes as well.
One measure of accuracy is the deviation from an 1/ field because this problem is optically-thin. To depict the increase in accuracy with ray sampling, we take the difference between the calculated photo-ionisation rate and a 1/ field, and then plot the standard deviation of this difference field versus angular resolution in Fig. 42. We plot this relation for resolutions of , , and and find no dependence on spatial resolution, which is expected because we control the angular resolution in terms of cell widths, not in absolute solid angles. We find that the deviation from an inverse square law decreases as .
7.3 Frequency resolution
The ionisation front radius is within 5–10% of analytical solutions in Tests 1, 2, and 5 with only one energy group; however a multi-frequency spectrum can create differences in the reactive flows. We use Test 5 (§5.1; an expanding H ii region with hydrodynamics) to probe any differences in the solution when varying the resolution of the spectrum. In RT09, ZEUS-MP was used to demonstrate the effect of a multi-frequency spectrum on the dynamics of the ionisation front in this test. Instead of a single shock seen in the mono-chromatic spectrum, the shock obtains a double-peaked structure in density and radial velocity. We rerun Test 5 with a K blackbody spectrum sampled by