RAMSES-RT: Radiation hydrodynamics in the cosmological context

We present a new implementation of radiation hydrodynamics (RHD) in the adaptive mesh refinement (AMR) code RAMSES. The multi-group radiative transfer (RT) is performed on the AMR grid with a first-order Godunov method using the M1 closure for the Eddington tensor, and is coupled to the hydrodynamics via non-equilibrium thermochemistry of hydrogen and helium. This moment-based approach has the large advantage that the computational cost is independent of the number of radiative sources - it can even deal with continuous regions of emission such as bound-free emission from gas. As it is built directly into RAMSES, the RT takes natural advantage of the refinement and parallelization strategies already in place. Since we use an explicit advection solver for the radiative transport, the time step is restricted by the speed of light - a severe limitation that can be alleviated using the so--called"reduced speed of light"approximation. We propose a rigorous framework to assess the validity of this approximation in various conditions encountered in cosmology and galaxy formation. We finally perform with our newly developed code a complete suite of RHD tests, comparing our results to other RHD codes. The tests demonstrate that our code performs very well and is ideally suited for exploring the effect of radiation on current scenarios of structure and galaxy formation.


INTRODUCTION
With the surging interest in reionization and the first sources of light in the Universe, and also thanks to a steadily increasing computational power, cosmological simulation codes have begun to include ionizing radiative transfer (RT) in the last decade or so. This is generally seen as a secondorder component in most astrophysical processes, but important nonetheless, and is obviously very important in the context of simulating reionization. Due to the challenges involved, most implementations have started out with the post-processing of ionizing radiation on simulations including only dark matter, but a few have begun doing coupled radiation hydrodynamics (RHD), which model the interplay of radiation and gas. E-mail: joki@strw.leidenuniv.nl It is highly desirable to follow self-consistently, with RHD simulations, the time-evolution and morphology of large-scale intergalactic medium (IGM) reionization and at the same time the smaller scale formation of the presumed sources of reionization; how galaxy formation is regulated by the ionizing radiation being released, how much of the radiation escapes from the galaxies to ionize the IGM, how first generation stars are formed in a metal-free environment and how radiative and supernovae feedback from those stars affect the inter-galactic medium. The galaxies and the IGM are indeed inter-connected via the ionizing radiation: the photons released from the galaxies affect the state of the surrounding gas via ionization and heating and may even prevent it from falling in or condensing into external gravitational potentials, especially small ones (e.g. Wise & Abel 2008;Ocvirk & Aubert 2011), which can then in turn significantly alter the ionization history.
The importance of RT and RHD is of course not limited to the epoch of reionization. Stars keep emitting ionizing radiation after this epoch and their radiative feedback likely has an effect on the post-reionization regulation of starformation (e.g. Pawlik & Schaye 2009;Hopkins, Quataert & Murray 2011), the mass distribution of stellar populations (Krumholz, Klein & McKee 2012) and even gas outflows (Hopkins et al. 2012).
Radiation hydrodynamics are complex and costly in simulations. The inclusion of coupled radiative transfer in hydrodynamical codes in general is challenging mainly because of the high dimensionality of radiative transfer (space-, angular-, and frequency dimensions) and the inherent difference between the typical timescales of radiative transfer and non-relativistic hydrodynamics. Simulating the interaction between small and large scales (so relevant to the epoch of cosmic reionization) makes things even worse: one wants to simulate, in a statistically significant region of the Universe (i.e. of the order of 100 comoving Mpc across) the condensation of matter in galaxy groups on Mpc scales, down to individual galaxies on kpc scales, followed by the formations of stellar nurseries in those galaxies on pc scales, and ultimately the formation of stars on sub-pc scales and then the effect of radiation from those stars back to the large scale IGM. This cycle involves size differences of something like 9 to 10 orders of magnitude -which is too much for the most advanced codes and computers today, actually even so without the inclusion of radiative transfer.
Due to these challenges, simulations typically focus on only a subset of these scales; either they consider reionization on large scales and apply sub-resolution recipes to determine stellar luminosities and UV escape fractions, or they ignore the cosmological context and focus on star formation and escape fractions in isolated galaxies or even isolated stellar nurseries.
A number of large scale 3D radiative transfer simulations of reionization have been carried out in recent years (e.g. Gnedin & Ostriker 1997;Miralda-Escudé, Haehnelt & Rees 2000;Gnedin 2000;Ciardi, Stoehr & White 2003;Sokasian et al. 2004;Iliev et al. 2006b;Zahn et al. 2007;Baek et al. 2010;Aubert & Teyssier 2010;Petkova & Springel 2011b), though they must all to some degree use subgrid recipes for star formation rates, stellar luminosities and UV escape fractions, none of which are well constrained. The ionization history in these simulations thus largely depends on these input parameters and resolution -some in fact use the observational constraints of the ionization history to derive constraints on these free parameters (e.g. Sokasian et al. 2004;Baek et al. 2010;Aubert & Teyssier 2010;Petkova & Springel 2011b). Furthermore, most of these works have used a postprocessing RT strategy instead of RHD, which neglects the effect the ionizing radiation has on the formation of luminous sources.
The primary driver behind this work is the desire to understand the birth of galaxies and stars during the dark ages, and how they link with their large scale environment. We have thus implemented a RHD version of the widely used cosmological code RAMSES (Teyssier 2002), that we call RAMSES-RT, with the goal of running cosmological RHD simulations, optimized for galactic scale radiation hydrodynamics. RAMSES is an adaptive mesh refinement (AMR) code, which greatly cuts costs by adaptively allowing the resolution to follow the formation of structures. The RHD implementation takes full advantage of the AMR strategy, allowing for high resolution simulations that can self consistently model the interplay of the reionizing Universe and the formation of the first galaxies.
Some of the goals we will be able to tackle with this implementation are: • Study radiative feedback effects in primordial galaxies. These galaxies are by definition young and small, and the first stars are thought to be gigantic and very bright due to the lack of metals. The ionizing radiation from these first stars is likely to have a dramatic effect on the galaxy evolution. This is closely associated with the formation of molecules, needed to form the first stars, which is sensitive to the radiation field. Radiative feedback effects also appear to be relevant in lower-redshift galaxies, and likely have a considerable impact on the initial mass function of stellar populations (Krumholz, Klein & McKee 2012).
• Investigate the escape of ionizing photons from early galaxies, how it affects the ionization history and external structure formation, e.g. the formation of satellite galaxies.
• Study the emission and absorption properties of galaxies and extended structures. Observable properties of gas are highly dependent on its ionization state, which in turn depends on the local radiation field (e.g. Oppenheimer & Schaye 2013). To predict it correctly, and to make correct interpretations of existing observations, one thus needs to model the ionization state consistently, for which RHD simulations are needed.
• Improve sub-resolution recipes: of course we have not implemented a miracle code, and we are still nowhere near simulating simultaneously the 9 to 10 orders of magnitude in scale needed for fully self-consistent simulations of reionization. Sub-resolution strategies are still needed, and part of the objective is to improve those via small-scale simulations of stellar feedback (SNe, radiation, stellar winds).
It is useful here to make clear the distinction between continuum and line radiative transfer: our goal is to study the interplay of ionizing radiation, e.g. from stellar populations and AGN, and the interstellar/intergalactic gas. We consider continuum radiation, because the spectra of stars (and AGN) are smooth enough that emission and absorption processes are not sensitive to subtle rest-frame frequency shifts, be they due to local gas velocities or cosmological expansion.
On the other side is line transfer, i.e. the propagation of radiation over a narrow frequency range, usually corresponding to a central frequency that resonates with the gas particles. An important example is the propagation of Lyα photons. Here, one is interested in the complex frequency and direction shifts that take place via scattering on the gas particles, and gas velocities and subtle frequency shifts are vital components. Line transfer is mostly done to interpret observational spectra, e.g. from Lyα emitting/absorbing galaxies (e.g. Verhamme, Schaerer & Maselli 2006), and is usually run in post-processing under the assumption that the line radiation has a negligible effect on the gas dynamics (through this assumption is not neccessarily true; see Dijkstra & Loeb 2009).
There is a bit of a grey line between those two regimes of continuum and line radiation -some codes are even able to do both (e.g. Baek et al. 2009;Pierleoni, Maselli & Ciardi 2009;Yajima et al. 2012). Our implementation deals strictly with continuum radiation though, as do most RHD implementations, for the sake of speed and memory limitations. We do approximate multi-frequency, but only quite coarsely, such that simulated photons represent an average of photons over a relatively wide frequency range, and any subtle frequency shifts and velocity effects are ignored.

Radiative transfer schemes and existing implementations
Cosmological hydrodynamics codes have traditionally been divided into two categories: Smoothed Particle Hydrodynamics (SPH) and AMR. The drawbacks and advantages of each method have been thoroughly explored (e.g. Agertz et al. 2007;Wadsley, Veeravalli & Couchman 2008;Tasker et al. 2008) and we now believe that both code types agree more or less on the final result if they are used carefully with recently developed fixes and improvements, and if applied in their regimes of validity. On the radiation side, it is quite remarkable that we have the same dichotomy between raytracing codes and moment-based codes. Comparative evaluations of both methods have been performed in several papers (Iliev et al. 2006a & Abel 2011), and here again, each method has its own specific advantage over the other one. Comparing both methods in the coupled case (RHD) within the more challenging context of galaxy formation, such as in the recent Aquila comparison project, (Scannapieco et al. 2012), remains to be done.

Ray-based schemes
Here the approximation is made that the radiation field is dominated by a limited number of sources. This allows one to approximate the local intensity of radiation, Iν , as a function of the optical depth τ along rays from each source. The simplest solution is to cast rays, or long characteristics from each source to each cell (or volume element) and sum up the optical depth at each endpoint. With the optical depths in hand, Iν is known everywhere and the rates of photoionization, heating and cooling can be calculated. While this strategy has the advantage of being simple and easy to parallelize (each source calculation is independent from the other), there is a lot of redundancy, since any cell which is close to a radiative source is traversed by many rays cast to further-lying cells, and is thus queried many times for its contribution to the optical depth. The parallelization is also not really so advantageous in the case of multiprocessor codes, since rays that travel over large lengths likely need to access cell states over many CPU nodes, calling for a lot of inter-node communication. Furthermore, the method is expensive: the computational cost scales linearly with the number of radiative sources, and each RT timestep has order O(Nsources N cells ) operations, where Nsources is the number of radiative sources and N cells is the number of volume el-ements. Implementation examples include Abel, Norman & Madau (1999), Cen (2002) and Susa (2006).
Short characteristics schemes overcome the redundancy problem by not casting separate rays for each destination cell. Instead, the calculation of optical depths in cells is propagated outwards from the source, and is in each cell based on the entering optical depths from the inner-lying cells. Calculation of the optical depth in a cell thus requires some sort of interpolation from the inner ones. There is no redundancy, as only a single ray segment is cast through each cell in one time-step. However, there is still a large number of operations and the problem has been made inherently serial, since the optical depths must be calculated in a sequence which follows the radiation ripple away from the source. Some examples are Nakamoto, Umemura & Susa (2001), Mellema et al. (2006), Whalen & Norman (2006) and Alvarez, Bromm & Shapiro (2006).
Adaptive ray tracing (e.g. Abel & Wandelt 2002;Razoumov & Cardall 2005;Wise & Abel 2011) is a variant on long characteristics, where rays of photons are integrated outwards from the source, updating the ray at every step of the way via absorption. To minimize redundancy, only a handful of rays are cast from the source, but they are split into sub-rays to ensure that all cells are covered by them, and they can be merged again if need be.
Cones are a variant on short characteristics, used in conjunction with SPH (Pawlik & Schaye 2008 and the moving-mesh AREPO code (Petkova & Springel 2011a). The angular dimension of the RT equation is discretized into tesselating cones that can collect radiation from multiple sources and thus ease the computational load and even allow for the inclusion of continuous sources, e.g. gas collisional recombination.
A hybrid method proposed by Rijkhorst et al. (2006) combines the long and short characteristics on patch-based grids (like AMR), to get rid of most of the redundancy while keeping the parallel nature. Long characteristics are used inside patches, while short characteristics are used for the inter-patches calculations.
Monte-Carlo schemes do without splitting or merging of rays, but instead reduce the computational cost by sampling the radiation field, typically both in the angular and frequency dimensions, into photon packets that are emitted and traced away from the source. The cost can thus be adjusted with the number of packets emitted, but generally this number must be high in order to minimize the noise inherent to such a statistical method. Examples include Ciardi et al. (2001), Maselli, Ferrara & Ciardi (2003), Altay, Croft & Pelupessy (2008), Baek et al. (2009), and Cantalupo & Porciani (2011). An advantage of the Monte-Carlo approach of tracking individual photon packets is that it naturally allows for keeping track of the scattering of photons. For line radiation transfer, where doppler/redshift effects in resonant photon scattering are important, Monte-Carlo schemes are the only feasible way to go -though in these cases, post-processing RT is usually sufficient (e.g. Cantalupo et al. 2005;Verhamme, Schaerer & Maselli 2006;Laursen & Sommer-Larsen 2007;Pierleoni, Maselli & Ciardi 2009).
Ray-based schemes in general assume infinite light speed, i.e. rays are cast from source to destination instantaneously. Many authors note that this only affects the initial speed of ionization fronts (I-fronts) around points sources (being faster than the light speed), but it may also result in an over-estimated I-front speed in underdense regions (see §6.5), and may thus give incorrect results in reionization experiments where voids are re-ionized too quickly. Some ray schemes (e.g. Pawlik & Schaye 2008;Petkova & Springel 2011a;Wise & Abel 2011) allow for finite light speed, but this adds to the complexity, memory requirement and computational load. With the exception of the cone-based methods (and to some degree the Wise & Abel (2011) implementation), which can combine radiation from many sources into single rays, ray-based schemes share the disadvantage that the computational load increases linearly with the number of radiative sources. Moment methods can naturally tackle this problem, though other limitations appear instead.

Moment-based radiative transfer
An alternative to ray-tracing schemes is to reduce the angular dimensions by taking angular moments of the RT equation (Eq. 2). Intuitively this can be thought of as switching from a beam description to that of a field or a fluid, where the individual beams are replaced with a "bulk" direction that represents an average of all the photons crossing a given volume element in space. This infers useful simplifications: two angular dimensions are eliminated from the problem, and the equations take a form of conservation laws, like the Euler equations of hydrodynamics. They are thus rather easily coupled to these equations, and can be solved with numerical methods designed for hydrodynamics. Since radiation is not tracked individually from each source, but rather just added to the radiation field, the computation load is naturally independent of the number of sources.
The main advantage is also the main drawback: the directionality is largely lost in the moment approximation and the radiation becomes somewhat diffusive, which is generally a good description in the optically thick limit, where the radiation scatters a lot, but not in the optically thin regime where the radiation is free-streaming. Radiation has a tendency to creep around corners with moment methods. Shadows are usually only coarsely approximated, if at all, though we will see e.g. in section §6.4 that sharp shadows can be maintained with idealized setups and a specific solver.
The large value of the speed of light is also an issue. Moment methods based on an explicit time marching scheme have to follow a Courant stability condition that basically limits the radiation from crossing more than one volume element in one time-step. This requires to perform many timesteps to simulate a light crossing time in the free-streaming limit, or, as we will see later, to reduce artificially the light speed. Implicit solvers can somewhat alleviate this limitation, at the price of inverting large sparse matrices which are usually ill-conditioned and require expensive, poorly parallelized, relaxation methods.
The frequency dimension is also reduced, via integration over frequency bins: in the grey (single group) approximation the integral is performed over the whole relevant frequency range, typically from the hydrogen ionization frequency and upwards. In the multigroup approximation, the frequency range is split into a handful of bins, or photon groups, (rarely more than a few tens due to memory and computational limitations) and the equations of radiative transfer can be solved separately for each group. Ray-tracing schemes also often discretize into some number of frequency bins, and they are usually more flexible in this regard than moment-based schemes: while the spectrum of each source can be discretized individually in ray-tracing, the discretization is fixed in space in moment-based schemes, i.e. the frequency intervals and resulting averaged photon properties must be the same everywhere, due to the field approximation.
In the simplest form of moment-based RT implementations, so-called flux limited diffusion (FLD), only the zero-th order moment of the radiative transfer equation is used, resulting in an elliptic set of conservation laws. A closure is provided in the form of a local diffusion relation, which lets the radiation flow in the direction of decreasing gas internal energy (i.e. in the direction opposite of the energy gradient). This is realistic only if the medium is optically thick, and shadows cannot be modelled. The FLD method has been used by e.g. Krumholz et al. (2007), Reynolds et al. (2009) and Commerçon et al. (2011), mainly for the purpose of studying the momentum feedback of infrared radiation onto dusty and optically thick gas, rather than photoionization of hydrogen and helium. Gnedin & Abel (2001) and Petkova & Springel (2009) used the optically thin variable Eddington tensor formalism (OTVET), in which the direction of the radiative field is composed on-the-fly in every point in space from all the radiative sources in the simulation, assuming that the medium between source and destination is transparent (hence optically thin). This calculation is pretty fast, given the number of relevant radiative sources is not overburdening, and one can neglect these in-between gas cells. Finlator,Özel & Davé (2009) take this further and include in the calculation the optical thickness between source and destination with a long characteristics method, which makes for an accurate but slow implementation. A clear disadvantage here is that in using the radiation sources to close the moment equations and compute the flux direction, the scaling of the computational load with the number of sources is re-introduced, hence negating one of the main advantages of moment-based RT. González, Audit & Huynh (2007), AT08 and Vaytet, Audit & Dubroca (2010) -and now us -use a different closure formalism, the so-called M1 closure, which can establish and retain bulk directionality of photon flows, and can to some degree model shadows behind opaque obstacles. The M1 closure is very advantageous in the sense that it is purely local, i.e. it requires no information which lies outside the cell, which is not the case for the OTVET approximation.
As shown by Dubroca & Feugeas (1999), the M1 closure has the further advantage that it makes the system of RT equations take locally the form of a hyperbolic system of conservations laws, where the characteristic wave speeds can be calculated explicitly and are usually close, but always smaller than the speed of light c. Hyperbolic systems of conservation laws are mathematically well understood and thoroughly investigated, and a plethora of numerical methods exist to deal with them (e.g. Toro 1999). In fact, the Euler equations are also a hyperbolic system of conservation laws, which implies we have the RT equations in a form which is well suited to lie alongside existing hydrodynamical solvers, e.g. in RAMSES.

From ATON to RAMSES-RT
The ATON code (AT08, Aubert & Teyssier 2010) uses graphical processing units, or GPUs, to post-process the transfer of monochromatic photons and their interaction with hydrogen gas. GPUs are very fast, and therefore offer the possibility to use the correct (very large) value for the speed of light and perform hundreds to thousands of radiation sub-cycles at a reasonable cost, but only if the data is optimally structured in memory, such that volume elements that are close in space are also close in memory. It is ideal for post-processing RT on simulation outputs that are projected onto a Cartesian grid, but hard to couple directly with an AMR grid in order to play an active part in any complex galaxy formation simulation. Even so, we have in the newest version of the ATON code included the possibility to perform fully coupled RHD simulations using a Cartesian grid only (this usually corresponds to our coarser grid level in the AMR hierarchy), where RT is performed using the ATON module on GPUs.
In our RAMSES-RT implementation we use the same RT method as ATON does -the moment method with the M1 Eddington tensor closure. The biggest difference is that RAMSES-RT is built directly into the RAMSES cosmological hydrodynamics code, allowing us to perform RHD simulations directly on the AMR grid, without any transfer of data between different grid structures. Furthermore, we have expanded the implementation to include multigroup photons to approximate multifrequency, and we have added the interactions between photons and helium. We explicitly store and advect the ionization states of hydrogen and helium, and we have built into RAMSES-RT a new non-equilibrium thermochemistry model that evolves these states along with the temperature and the radiation field through chemical processes, photon absorption and emission. Finally, for realistic radiative feedback from stellar populations, we have enabled RAMSES-RT to read external spectral energy distribution (SED) models and derive from them luminosities and UV "colors" of simulated stellar sources.
We have already listed a number of RT implementations, two of which even function already in the RAMSES code (AT08, Commerçon et al. 2011), and one might ask whether another one is really needed. To first answer for the ATON implementation, it is optimized for a different regime than RAMSES-RT. As discussed, ATON prefers to work with structured grids, but it cannot deal well with adaptive refinement. This, plus the speed of ATON, makes it very good for studying large scale cosmological reionization, but not good for AMR simulations of individual halos/galaxies, e.g. cosmological zoom simulations, where the subject of interest is the effect of radiative feedback on the formation of structures and galaxy evolution, and escape fractions of ionizing radiation. The Commerçon et al. (2011) implementation is on the opposite side of the spectrum. Being based on the FLD method, it is optimized for RHD simulations of optically thick protostellar gas. It is a monogroup code that doesn't track the ionization state of the gas. Furthermore, it uses a rather costly implicit solver, which makes it hard to adapt to multiple adaptive time stepping usually used in galaxy formation problems.
Thus there aren't so many cosmological RHD implementations out there, and there should be room for more. The main advantage of our implementation is that our method allows for an unlimited number of radiative sources and can even easily handle continuous sources, and is thus ideal for modelling e.g. the effects of radiative feedback in highly resolved simulations of galaxy formation, UV escape fractions, and the effects of self-shielding on the emission properties of gas and structure formation, e.g. in the context of galaxy formation in weak gravitational potentials.
The structure of the paper is as follows: in §2 we present the moment based RT method we use. In §3 we explain how we inject and transport radiation on a grid of gas cells, and how we calculate the thermochemistry in each cell, that incorporates the absorption and emission of radiation. In §4 we present two tricks we use to speed up the RHD code, namely to reduce the speed of light, and to "smooth" out the effect of operator splitting. In §5 we describe how the radiative transfer calculation is placed in the numerical scheme of RAMSES, and demonstrate that the radiation is accurately transported across an AMR grid. In §6, we present our test suite, demonstrating that our code performs very well in coupled radiation hydrodynamics problems and finally, §7 summarizes this work and points towards features that may be added in the future. Details of the thermochemistry and additional code tests are described in the appendix.

MOMENT-BASED RADIATIVE TRANSFER WITH THE M1 CLOSURE
Let Iν (x, n, t) denote the radiation specific intensity at location x and time t, such that is the energy of photons with frequency over the range dν around ν propagating through the area dA in a solid angle dΩ around the direction n. The equation of radiative transfer (e.g. Mihalas & Mihalas 1984) describes the local change in Iν as a function of propagation, absorption and emission, where c is the speed of light, κν (x, n, t) is an absorption coefficient and ην (x, n, t) a source function. By taking the zeroth and first angular moments of (2), we can derive the moment-based RT equations that describe the time-evolution of photon number density Nν and flux Fν (see e.g. AT08): where Pν is the radiative pressure tensor that remains to be determined to close the set of equations. Here we have split the absorption coefficient into constituent terms, njσνj, where nj is number density of the photo-absorbing species j (=Hi, Hei, Heii), and σνj is the ionization cross section between ν-frequency photons and species j. Furthermore we have split the source function into (e.g. stellar, quasar) injection sources,Ṅ ν , and recombination radiation from gas, N rec ν . Here we only consider the photo-absorption of hydrogen and helium, which is obviously most relevant in the regime of UV photons. However, other absorbers can straightforwardly be added to the system. Eqs.
(3)-(4) are continuous in ν, and they must be discretized to be usable in a numerical code. AT08 collected all relevant frequencies into one bin, so the equations could be solved for one group of photons whose attributes represent averages over the frequency range. For a rough approximation of multifrequency, we split the relevant frequency range into a number of photon groups, defined by where (νi0, νi1) is the frequency interval for group i. In the limit of one photon group, the frequency range is (νi0, νi1) = (νHi, ∞); with M > 1 groups, the frequency intervals should typically be mutually exclusive and set up to cover the whole H-ionizing range: [ν00, ν01 : ν10, ν11 : ... : νM0, νM1] = [νHi, ∞[.
Integrating the RT equations (3) and (4) over each frequency bin corresponding to the group definitions yields M sets of four equations: where σ N ij represent average cross sections between each group i and species j, defined by 1 We simplify things however by defining the group cross sections as global quantities, assuming a frequency distribution of energy J(ν) for the radiative sources (e.g. a black-body or some sophisticated model). The cross sections are thus in practice evaluated by where hν is photon energy (with h the Planck constant). Likewise, average photon energies within each group are evaluated by¯ and furthermore, for the calculation of photoionization heating 2 , energy weighted cross sections are stored for each group -absorbing species couple: In RAMSES-RT, σ N ij , σ E ij and¯ i can be either set by hand or evaluated on-the-fly from spectral energy distribution tables as luminosity weighted averages from in-simulation stellar populations, using the expressions from Verner et al. (1996) for σν,Hi, σν,Hei and σν,Heii.
For each photon group, the corresponding set of equations (6)-(7) must be closed with an expression for the pressure tensor P. This tensor is usually described as the product of the photon number density and the so-called Eddington tensor D (see Eq. 12), for which some meaningful and physical expression is desired. Some formalisms have been suggested for Dν . Gnedin & Abel (2001), Finlator,Özel & Davé (2009), and Petkova & Springel (2009 have used the so called optically thin Eddington tensor formalism (OTVET), in which P is composed on-the-fly from all the radiation sources, the main drawback being the computational cost associated with collecting the positions of every radiative source relative to every volume element. Instead, like AT08 (and González, Audit & Huynh 2007 before them), we use the M1 closure relation (Levermore 1984), which has the great advantages that it is purely local, i.e. evaluating it in a piece of space only requires local quantities, and that it can retain a directionality along the flow of the radiative field. In our frequency-discretized form, the pressure tensor is given in each volume element for each photon group by where the Eddington tensor is and are the unit vector pointing in the flux direction, the Eddington factor and the reduced flux, respectively. The reduced flux describes the directionality of the group i radiation in each point, and must always have 0 fi 1. A low value means the radiation is predominantly isotropic, and a high value means it is predominantly flowing in one direction.
Photons injected into a point (via an increase in photon density only) initially have zero reduced flux and thus are isotropic. Away from the source, the moment equations and M1 closure develop a preferred outwards direction, i.e. the reduced flux tends towards one. Beams can be injected by imposing unity reduced flux on the injected photons. In this case, the M1 closure correctly maintains unity reduced flux (and χ = 1) along the beam (see demonstrations in Fig. 1 and Sections 5.3, 6.4, and 6.8). For the arguments leading to these expressions and a general discussion we point the reader towards Levermore (1984) and González, Audit & Huynh (2007) and AT08.

THE RADIATIVE TRANSFER IMPLEMENTATION
We will now describe how pure radiative transfer is solved on a grid -without yet taking into consideration the hydrodynamical coupling. The details here are not very specific to RAMSES-RT and are much like those of AT08.
In addition to the usual hydrodynamical variables stored in every grid cell in RAMSES (gas density ρ, momentum density ρu, energy density E, metallicity Z), RAMSES-RT has the following variables: First, we have the 4 × M variables describing photon densities Ni and fluxes Fi for the M photon groups. Second, in order to consistently treat the interactions of photons and gas, we track the non-equilibrium evolution of hydrogen and helium ionization in every cell, stored in the form of passive scalars which are advected with the gas, namely For each photon group, we solve the set of equations (6)-(7) with an operator splitting strategy, which involves decomposing the equations into three steps that are executed in sequence over the same time-step ∆t, which has some pre-determined length. The steps are: (i) Photon injection step, where radiation from stellar and other radiative sources (other than gas recombinations) is injected into the grid. This corresponds to theṄ i term in (6).
(ii) Photon transport step, where photons are propagated in space. This corresponds to solving (6)-(7) with the RHS being equal to zero.
(iii) Thermochemistry step, where the rest of the RHS of (6)-(7) is solved. This is where the photons and the gas couple, so here we evolve not only the photon densities and fluxes, but also the ionization state and temperature of the gas.

The injection step
The equations to solve in this step are very simple, whereṄ i is a rate of photon injection into photon group i, in the given cell. Normally, the injected photons come from stellar sources, but they could also include other point sources such as AGN, and also pre-defined point sources or even continuous "volume" sources 3 . Given the time t and time-step length ∆t, the discrete update in each cell done for each photon group is the following sum over all stellar particles situated in the cell: where n denotes the time index (n = t and n + 1 = t + ∆t), fesc is an escape fraction, V is the cell volume, m , τ and Z are mass, age and metallicity of the stellar particles, respectively, and Πi is some model for the accumulated number of group i photons emitted per solar mass over the lifetime (so far) of a stellar particle. The escape fraction, fesc is just a parameter that can be used to express the suppression (or even boosting) of radiation from processes that are unresolved inside the gas cell. RAMSES-RT can read stellar energy distribution (SED) model tables to do on-the-fly evaluation of the stellar particle luminosities, Πi. Photon cross sections and energies can also be determined on-the-fly from the same tables, to represent luminosity-weighted averages of the stellar populations in a simulation. Details are given in Appendix B.

The transport step
The equations describing free-flowing photons are i.e. (6)-(7) with the RHS = 0. Note that we have removed the photon group subscript, since this set of equations is solved independently for each group over the time-step. We can write the above equations in vector form where U = [N, F] and F(U) = [F, c 2 P]. To solve (20) over time-step ∆t, we use an explicit conservative formulation, expressed here in 1D for simplicity, where n again denotes time index and l denotes cell index along the x-axis. F l+1/2 and F l−1/2 = F (l−1)+1/2 are intercell fluxes evaluated at the cell interfaces. Simple algebra gives us the updated cell state, the form of (22) (see e.g. Toro 1999), as long as the Courant time-step condition is respected (see §4.1). Following AT08 and González, Audit & Huynh (2007) we implement two flux functions which can be used in RAMSES-RT.
One is the Harten-Lax-van Leer (HLL) flux function (Harten, Lax & Leer 1983), are maximum and minimum eigenvalues of the Jacobian ∂F/∂U. These eigenvalues mathematically correspond to wave speeds, which in the case of 3D radiative transfer depend only on the magnitude of the reduced flux f (14) and the angle of incidence of the flux vector to the cell interface. This dependence has been calculated and tabulated by González, Audit & Huynh (2007), and we use their table to extract the eigenvalues. The other flux function we have implemented is the simpler Global Lax Friedrich (GLF) function, which corresponds to setting the HLL eigenvalues to the speed of light, i.e. λ − = −c and λ + = c, and has the effect of making the radiative transport more diffusive. Beams and shadows are therefore better modelled with the HLL flux function than with the GLF one, whereas the inherent directionality in the HLL function results in radiation around isotropic sources (e.g. stars) which is noticeably asymmetric, due to the preference of the axis directions. Fig. 1 illustrates the difference between the two flux functions in some idealized 2D RAMSES-RT tests, where we shoot off beams and turn on isotropic sources. Here the photon-gas interaction is turned off by setting all photoionization cross sections to zero (σ N j = σ E j = 0 for any species j). It can be seen that the HLL flux function fails to give isotropic radiation (far left) and that the GLF function gives more diffusive beams (second from left). Note also how the diffusivity of beams with the HLL flux function is directiondependent. A horizontal or vertical beam is perfectly retained while a diagonal one "leaks" to the sides almost as much as with the GLF function, which has the advantage of being fairly consistent on whether the beam is along-axis or diagonal. The right frames of the figure give an idea of how the radiative transport behaves in the case of multiple sources, i.e. with opposing beams and neighboring isotropic sources. The two opposing beams example is a typical configuration where the M1 closure relation obviously fails, creating a spurious source of radiation, perpendicular to the beam direction: Since opposing fluxes cannot cross each other in a single point in the moment approximation, the radiation is "squeezed" into those perpendicular directions. It is unclear to us how much of a problem this presents in astrophysical contexts. Beams, which clearly represent the worst case scenario, are not very relevant, but multiple nearby sources are. We generally prefer to use the GLF flux function, since we mostly deal with isotropic sources in our cosmological/galac-tic simulations, but the choice of function really depends on the problem. There is no noticeable difference in the computational load, so if shadows are important, one should go for HLL. AT08 have compared the two flux functions in some of the benchmark RT tests of Iliev et al. (2006a) and found that they give very similar results. We do likewise for the test we describe in §6.7, and come to the same conclusion.

The thermochemical step
Here we solve for the interaction between photons and gas. This is done by solving (6) and (7) with zero divergence and stellar injection terms.
Photon absorption and emission have the effect of heating and cooling the gas, so in order to self-consistently implement these interactions, we evolve along with them the thermal energy density ε of the gas and the abundances of the species that interact with the photons, here Hi, Hei and Heii via photoionizations and Hii, Heii (again) and Heiii via recombinations. We follow these abundances in the form of the three ionization fractions xHII, xHeII and xHeIII, that we presented in Eqs. (15). The set of non-equilibrium thermochemistry equations solved in RAMSES-RT consists of: In the photon density and flux equations, (25) and (26), we have replaced the photon emission rateṄ rec i with the full expression for recombinative emissions from gas. Here, α A j (T ) and α B j (T ) represent case A and B recombination rates for electrons combining with species j (= Hii, Heii, Heiii). The b rec ji factor is a boolean (1/0) that states which photon group j-species recombinations emit into, and ne is electron number density (a direct function of the H and He ionization states, neglecting the contribution from metals).
The computational approach we use to solving Equations (25)-(30) takes inspiration from Anninos et al. (1997). The basic premise is to solve the equations over a sub-step in a specific order (the order we have given), explicitly for those variables that remain to be solved (including the current one), but implicitly for those that have already be solved over the sub-step. Eqs. (25) and (26) are thus solved purely explicitly, using the backwards-in-time (BW) values for all variables on the RHS. Eq. (27) is partly implicit in the sense that it uses forward-in-time values for N and F, but BW values for the other variables. And so on, ending with Eq. (30), which is then implicit in every variable except the one solved for (xHeIII). We give details of the discretization of these equations in Appendix A.

The 10% thermochemistry rule
For accuracy, each thermochemistry step is restricted by a local cooling time which prohibits any of the thermochemical quantities to change by a substantial fraction in one time-step. We therefore sub-cycle the thermochemistry step to fill in the global RT time-step (see next section), using what can be called the 10% rule: In each cell, the thermochemistry step is initially executed with the full RT timestep length, and then the fractional update is considered. If any of the evolved quantities (Ni, Fi, ε, ionization fractions) have changed by more than 10%, we backtrack and do the same calculation with half the time-step length. Conversely, if the greatest fractional change in a sub-step is < 5%, the timestep length is doubled for the next sub-step (without the backtracking).
Together, the quasi-implicit approach used in solving the thermochemistry, and the 10% rule, infer that photons are in principle conserved only at the 10% level 4 . This is because the thermochemistry solver is explicit in the photon density updates (i.e. uses before-timestep values of ionization fractions), but the following ionization fraction updates are implicit in the photon densities (i.e. they use aftertimestep values for the photon densities). Thus, in the situation of a cell in the process of being photoionized, the ionization fractions are underestimated at the photon density updates and the photon densities are underestimated at the ionization fraction updates. Conversely, if the cell gas is recombining, the recombination-photon emission is slightly overestimated, since before-timestep values for the ionization fractions are used for the emissivity. However, judging from the performance in RT tests ( §6) and thermochemistry tests ( §C), this does not appear to be cause for concern.
3.3.2 The on-the-spot approximation The photon-emitting recombinative term, the second RHS sum in (25), is optionally included. Excluding it is usually referred to as the on-the-spot approximation (OTSA), meaning that any recombination-emitted photons are absorbed "on the spot" by a near-lying atom (in the same cell), and hence these photon emissions cancel out by local photon absorptions. If the OTSA is assumed, the gas is thus not photoemitting, and the case A recombination rates are replaced with case B recombination rates in (25)-(30), i.e. photonemitting recombinations straight down to the ground level are not counted. The OTSA is in general a valid approximation in the optically thick regime but not so when the photon mean free path becomes longer than the cell width.
It is a great advantage of our RT implementation that it is not restricted to a limited number of point sources. The computational load does not scale at all with the number of sources, and photon emission from gas (non-OTSA) comes at no added cost, whereas it may become prohibitively expensive in ray-tracing implementations.

TIMESTEPPING ISSUES
RT is computationally expensive, and we use two basic tricks to speed up the calculation. One is to reduce the speed of light, the other is to modify slightly the traditional operator splitting approach, by increasing the coupling between photon injection and advection on one hand and thermochemistry and photo-heating on the other hand.

The RT timestep and the reduced speed of light
In each iteration before the three RT steps of photon injection, advection and thermochemistry are executed, the length of the time-step, ∆tRT, must be determined. We use an explicit solver for the radiative transport (21), so the advection timestep, and thus the global RT timestep, is constrained by the Courant condition (here in 3D), where ∆x is the cell width. This time-step constraint is severe: it results in an integration step which is typically 300 times shorter than in non-relativistic hydrodynamical simulations, where the speed of light is replaced by a maximum gas velocity (∼ 1000 km/s) in Eq. 45. In a coupled (RHD) simulation, this would imply a comparable increase in CPU time, either because of a global timestep reduction (as we chose to implement, see Sec. 5), or because of many radiative sub-steps (as is implemented e.g. in ATON 5 ). In the case of radiative transfer with the moment equations, there are two well-known solutions to this problem. The first solution is to use an implicit method instead of an explicit one to solve the transport equation, which means using forward-in-time intercell fluxes in (21), i.e replacing F n ≡ F t with F n+1 ≡ F t+∆t . This seemingly simple change ensures that the computation is always stable, no matter how big the time-step, and we can get rid of the Courant condition. However: (i) It doesn't mean that the computation is accurate, and in fact we still need some time-stepping condition to retain the accuracy, e.g. to restrain any quantity to be changed by more than say 10% in a single time-step. Furthermore, such a condition usually must be checked by trial-and-error, i.e. one guesses a timestep and performs a global transport step (over the grid) and then checks whether the accuracy constraint was broken anywhere. Such trial-and-error time-stepping can be very expensive since it is a global process. (ii) Replacing F t with F t+∆t is actually not simple at all. Eqs. 21 become a system of coupled algebraic equations that must be solved via matrix manipulation in an iterative process, which is complicated, computationally expensive, and of limited scope (i.e. can't be easily applied to any problem). Due to these two reasons we have opted out of the implicit approach. It is absolutely a valid approach however, and used by many (e.g. Petkova & Springel 2009;Commerçon et al. 2011).
The second solution, which we have chosen instead, is to keep our solver explicit, and relax the Courant condition by changing the speed of light to a reduced light speed cr c, the payoff being that the time-step (45) becomes longer. This is generally referred to as the reduced speed of light approximation (RSLA), and was introduced by Gnedin & Abel (2001). The idea of the RSLA is that in many applications of interest, the propagation of light is in fact limited by the much slower speed of ionizing fronts. In such situations, reducing the speed of light, while keeping it higher than the fastest I-front, will yield the correct solution at a much reduced CPU cost. In the following section, we provide a framework to help judge how accurate the RSLA may be in various astrophysical contexts.

A framework for setting the reduced light speed value
In the extremely complex framework of galaxy formation simulations, the accuracy of the results obtained using the RSLA can really only be assessed by convergence tests. It is nonetheless useful to consider a simple idealized setup in order to derive a physical intuition of where, when, and by how much one may reduce the speed of light. In this section, we thus discuss the expansion of an ionized region around a central source embedded in a uniform neutral medium.
We consider a source turning on and emitting ionizing photons at a rateṄ into a homogeneous hydrogen-only medium of number density nH. An expanding sphere of ionized gas forms around the source and halts at the Strömgren radius rS within which the rate of recombinations equals the source luminosity: where α B ∼ 2.6 × 10 −13 cm 3 s −1 is the case-B recombination rate at T ∼ 10 4 K, and where we have assumed that the plasma within rS is fully ionized. The relativistic expansion of the I-front to its final radius rS is derived in Shapiro et al. (2006), and may be expressed as: where w = t/trec is time in units of the recombination time trec = (nHα B ) −1 , y = rI/rS is the position of the I-front in units of rS, and the factor q ≡ tcross/trec ≡ rS/(ctrec) describes the light crossing time tcross across the Strömgren radius in units of the recombination time, and basically encompasses all the free parameters in the setup (source luminosity, gas density, and temperature). Writing q ∝Ṅ 1/3 n 1/3 H , we see that in many astrophysical contexts, q stays in the range ∼ 10 −3 − 10 −2 (see Table 1), simply because we are generally either interested in the effect of bright sources (e.g. a whole galaxy) on relatively low-density gas (e.g. the IGM) or of fainter sources (e.g. an O-star) on high-density gas (at e.g. molecular-cloud densities).
Let us now discuss briefly the evolution of an I-front given by Eq. 33 for illustrative values of q: • q = 0 (blue curve of Fig. 2): this is the limiting nonrelativistic case, which assumes an infinite speed of light (tcross = 0). In this case, the I-front expands roughly as y ∝ w 1/3 (its speed decreases as w −2/3 ) almost all the way to rS, which it reaches after about a recombination time.
• q = 1 (red curve of Fig. 2): here (and for all q > 1), the I-front basically expands at the speed of light all the way to rS, which it thus reaches after a crossing time (which is equal to a recombination time in this case).
• q = 10 −3 (green curve of Fig. 2): in this typical case, the I-front starts expanding at the speed of light, until w ∼ (q/3) 3/2 . It then slows down and quickly reaches the limiting q = 0 behavior after a crossing time (at w ∼ q). The I-front then reaches rS after a recombination time (at w ∼ 1).
An important feature appearing in the two latter cases is that for any physical setup q > 0, the I-front is always well described by the q = 0 limit after a crossing time (i.e. w q). We can use this feature to understand the impact of reducing the speed of light in our code. Say we have a physical setup described by a value q0. Reducing the speed of light by a factor fc < 1 (cr = fcc) implies an increase by a factor 1/fc of the effective crossing time, and the effective q in our experiment becomes q0/fc. The solution we obtain with cr will be accurate only after an effective crossing time, i.e. after w = q0/fc. Before that time, the reduced-light-speed solution will lag behind the real one.
How much one may reduce the speed of light in a given numerical experiment then depends on the boundary conditions of the problem and their associated timescales. Call Figure 2. I-front expansion in a Stromgren sphere for a set of values of the dimensionless crossing time q. The blue curve shows the infinite-light-speed limit (q = 0). The green curve shows a typical case with q = 10 −3 , and the red curve shows the q = 1 case, as discussed in the text. The thin grey curves show other values of q, spanning the range 10 −4 − 10 in steps of one dex. The grey lines in the top-left corner of the plot show slopes corresponding to an expansion at the speed of light (dot-dashed line) or as (t/trec) 1/3 (dashed line). For any q > 0, the I-front radius is accurately described by the q = 0-limit after a crossing time.
τsim the shortest relevant timescale of a simulation. For example, if one is interested in the effect of radiative feedback from massive stars onto the ISM, τsim can be set to the lifetime of these stars. If one is running a very short experiment (see Sec. 6.5), the duration of the simulation may determine τsim. Given this timescale contraint τsim, one may reduce the speed of light by a factor such that the I-fronts will be correctly described after a timelapse well shorter than τsim, i.e. tcross/fc τsim. In other words, one may typically use fc = min(1; ∼ 10 × tcross/τsim). We now turn to a couple of concrete examples.

Example speed of light calculations
In Table 1 we take some concrete (and of course very approximate) examples to see generally what values of fc are feasible. We consider cosmological applications from intergalactic to inter-stellar scales and setups from some of the RT code tests described in §6.

Reionization of the inter-galactic medium
Here we are concerned with the expansion of ionization fronts away from galaxies and into the IGM, as for example in the fourth test of Iliev et al. (2006a) (hereafter Il06). In this test, the IGM gas density is typically nH = 10 −4 cm −3 , and the sources haveṄ = 7 10 52 s −1 . In such a configuration, the Strömgren radius is rS ∼ 600 kpc, corresponding to a crossing time tcross ∼ 2 Myr. Because of the low density of the gas, the recombination time is very long ( 1 Gyr), and we are thus close to the q = 10 −3 case discussed above (the green curve in Fig. 2).
Test 4 of Il06 is analyzed at output times τsim,1 = 0.05 Myr and τsim,2 = 0.4 Myr (see Fig. 19). In both cases, τsim < MW ISM 10 −1 2 10 50 0.9 3 10 −3 1.2 2 10 −3 1 1 3 10 −2 MW cloud 10 2 2 10 48 2 10 −3 6 10 −6 1 10 −3 5 10 −3 0.1 80 6 10 −4 Iliev tests 1,2,5 10 −3 5 10 48 5.4 2 10 −2 122.3 1.4 10 −4 10 8 10 −2 2 10 −2 Iliev test 4 10 −4 7 10 52 600 2 1200 2 10 −3 0.05 4 10 −5 1 tcross, and we cannot reduce the speed of light to get an accurate result at these times, because the expanding front has not yet reached the q = 0 limit. Interestingly, we cannot increase the speed of light either, as is done in Il06 with C 2 -Ray which assumes an infinite light-speed. From Fig. 2, it is clear that this approximation (the blue curve) will overpredict the radius of the front. We can use the analysis above to note that had the results been compared at a later output time τsim > 2 Myr, the infinite-light-speed approximation would have provided accurate results. It is only ten times later, however, that reducing the speed of light by a factor ten would have provided accurate results. We conclude that propagating an I-front in the IGM at the proper speed requires to use a value of the speed of light close to the correct value. This is especially true in Test 4 of the Il06 paper (last row of Table 1). This confirms that for cosmic reionization related studies, using the correct value for the speed of light is very important.

Inter-stellar medium
There is admittedly a lot of variety here, but as a rough estimate, we can take typical densities to be nH ∼ 10 −1 cm −3 in the large-scale ISM and nH ∼ 10 2 cm −3 in star-forming clouds. In the stellar nurseries we consider single OB stars, releasingṄOB ∼ 2 × 10 48 photons per second, and in the large-scale ISM we consider groups of (∼ 100) OB stars. The constraining timescale is on the order of the stellar cycle of OB stars (τsim ∼ 10 Myr), and less for the stellar nurseries. In these two cases, which are representative of the dense ISM inside galactic disks, we see in Table 1 that the allowed reduction factor for the speed of light is much larger (fc 10 −4 to 10 −3 ). This is due to two effects acting together: the gas density is higher, but the sources are fainter, since we are now resolving individual stellar clusters, and not an entire galaxy. Tests 1, 2 of Il06 and test 5 of its' RHD sequel (Iliev et al. 2009) are also representative of such a favorable regime to use the reduced speed of light approximation (second to last row in Table 1). This rigorous analysis of the problem at hand confirms that propagating I-front in galaxy formation simulation can be done reliably using our current approach, while cosmic reionization problems are better handled with GPU acceleration and the correct speed of light.

Smoothed RT
A problem we had to face, while performing RAMSES-RT galaxy formation runs, as well as the various test cases presented here, is that there is often a small number of cells, usually along I-fronts, or close to strong radiation sources, that execute a huge number of thermochemistry subcycles in a single RT time-step. This is in part fault of the operatorsplitting approach used, where the RT equations have been partly decoupled. Specifically, the photon density updates happen in three steps in this approach (see Fig. 3, top). The photon injection step always increases the number of photons, usually by a relatively large amount, and the transport step does the same when it feeds photons into cells along these I-fronts. The thermochemistry step in the I-front cells has the exact opposite effect: the photon density decreases again via absorptions. If the photon-depletion time is shorter than the Courant time, we have a curious situation where the cell goes through an inefficient cycle during the thermochemistry subcycles: it starts neutral with a large abundance of photons (that have come in via the transport and/or photon injection steps). It first requires a number of subcycles to evolve to a (partly) ionized state, during which the photon density is gradually decreased. It can then reach a turnaround when the photons are depleted. If the RT timestep is not yet finished, the cell then goes into a reverse process, where it becomes neutral again. This whole cycle may take a large number of thermochemical steps, yet the cell gas ends up being in much the same state as it started.
In reality, the ionization state and photon density would not cycle like this but would rather settle into a semiequilibrium where the rate of ionizations equals that of recombinations.
For the purpose of saving up on computing time and reducing the number of thermochemistry subcycles, we have implemented an optional strategy we call smoothed RT that roughly corrects this non-equilibrium effect of operator splitting (see Fig. 3, bottom). In it, the result of (N i , F i ) from the transport and injection steps in each cell is used to infer a rate for the thermochemistry step, rather than being set as an initial condition. We use the pre-transport, preinjection values of N t i and F t i as initial conditions for the thermochemistry, but instead update the thermochemistry equations (25) and (26) to where the new terms at the far right represent the rates at which the photon densities and fluxes changed in the transport and injection steps, i.e.
where N t i and F t i (N i and F i ) denote a cell state before (after) cell injection (Eq. 16) and transport (Eqs. 18 and 19) Figure 3. Sketch plots showing a photon density evolution over a global RT time-step with normal RT (top) and smoothed RT (bottom). In normal RT the photon density is updated to N during photon transport (a) and injection (b). This is then used as an initial state for thermochemistry (c). It is often the case that the photons are depleted over the global time-step ∆t, in a process which takes many thermochemistry subcycles. In smoothed RT, the photon density state is not updated by the transport and injection steps, but rather the difference is used to infer a photon injection rate for the cell, which is gradually added during each thermochemistry substep. This can dramatically reduce the needed number of chemistry substeps.
have been solved over ∆t. The injection and transport steps are unchanged from the normal operator splitting method, except for the fact that the cell states are not immediately updated to reflect the end results of those steps. The results of the injection step go only as initial conditions into the transport step, and the end results of the transport step are only used to calculate the photon density and flux rates of change via Eqs. (36) and (37). Only after the thermochemistry step does a cell get a valid state that is the result of all three steps. The idea is that when the photons are introduced like this into the thermochemistry step, they will be introduced gradually in line with the subcycling, and the photon density vs. ionization fraction cycle will disappear as a result and be replaced with a semi-equilibrium, which should reduce the number of subcycles and the computational load. The total photon injection (or depletion) will still equal N i − N t i , so in the limit that there are no photoionizations or photonemitting recombinations, the end result is exactly the same photon density (and flux) as would be left at the end of the transport and injection steps without smoothed RT.
The advantage of the smoothing approach is perhaps best explained with an example: consider a cell with a strong source of radiation and gas dense and neutral enough that the timescales of cooling, ionization and/or recombination are much shorter than the global timestep length, ∆t. This could either be a source containing a stellar particle or a cell along an ionization front. Without smoothing, the photoionization rate in the cell can change dramatically as a result of photon injection/transport. The thermochemistry step thus starts with a high rate of photoionzations which gradually goes down in the thermochemistry sub-cycling as the gas becomes more ionized and the photons are absorbed. With smoothing, this dramatic change in the photoionization rate never happens, thus requiring fewer thermochemistry sub-cycles to react. A situation also exists where the smoothing approach slows down the thermochemistry: if a cell contains a strong source of radiation, but diffuse gas (i.e. long timescales compared to ∆t for cooling, ionization and/or recombination), the non-smoothed approach would result in little or no thermochemistry sub-cycling, whereas the smoothed approach would take many sub-cycles just to update the radiation field and effectively reach the final result of the injection and transport steps. The gain in computational speed is thus quite dependent on the problem at hand, and also on the reduced light speed, which determines the size of the RT time-step, ∆t. We've made a comparison on the computational speed between using the smoothed and non-smoothed RT in a cosmological zoom simulation from the NUT simulations suite (e.g. Powell, Slyz & Devriendt 2011) that includes the transfer of UV photons from stellar sources. Here, smoothed RT reduces the average number of thermochemistry subcycles by a factor of 6 and the computing time by a factor 3.5. So a lot may indeed be gained by using smoothed RT.
One could argue that the ionization states in I-fronts are better modelled with smoothed RT, since the cycle of photon density and ionization fraction is a purely numerical effect of operator splitting. We have intentionally drawn a slightly higher end value of Ni in the smoothed RT than non-smoothed in Fig. 3: whereas non-smoothed RT can completely deplete the photons in a cell, smoothed RT usually leaves a small reservoir after the thermochemistry, that more accurately represents the "semi-equilibrium value".
Of course an alternative to smoothed RT, and a more correct solution, is to attack the root of the problem and reduce the global time-step length, i.e. also limit the transport and injection steps to the 10% rule. Reducing the global time-step length is highly impractical though; the main reason for using operator splitting in the first place is that it enables us to separate the timescales for the different steps.
The same method of smoothing out discreteness that comes with operator splitting (in the case of pure hydrodynamics) has previously been described by Sun (1996), where it is referred to as "pseudo-non-time-splitting".

RADIATION HYDRODYNAMICS IN RAMSES
RAMSES (Teyssier 2002) is a cosmological adaptive mesh refinement (AMR) code that can simulate the evolution and interaction of dark matter, stellar populations and baryonic gas via gravity, hydrodynamics and radiative cooling. It can run on parallel computers using the message passing interface (MPI) standard, and is optimized to run very large numerical experiments. It is used for cosmological simulations in the framework of the expanding Universe, and also smaller scale simulations of more isolated phenomena, such as the formation and evolution of galaxies, clusters, and stars. Dark matter and stars are modelled as collisionless particles that move around the simulation box and interact via gravity. We will focus here on the hydrodynamics of RAMSES though, which is where the RT couples to everything else. RAMSES employs a second-order Godunov solver on the Euler equations of gravito hydrodynamics in their conservative form, where t is time, ρ the gas density, u the bulk velocity, φ the gravitational potential, E the gas total energy density, P the pressure, and Λ represents radiative cooling and heating via thermochemistry terms (resp. negative and positive), which are functions of the gas density, temperature and ionization state. In RAMSES, collisional ionization equilibrium (CIE) is traditionally assumed, which allows the ionization states to be calculated as surjective functions of the temperature and density and thus they don't need to be explicitly tracked in the code. E is divided into kinetic and thermal energy density (ε) components: The system of Euler equations is closed with an equation of state which relates the pressure and thermal energy, where γ is the ratio of specific heats. The Euler equations are adapted to super comoving coordinates, to account for cosmological expansion, by a simple transformation of variables (see §5.4).
The Euler equations are solved across an AMR grid structure. Operator splitting is employed for the thermochemistry source terms, i.e. Λ is separated from the rest of the Euler equations in the numerical implementationwhich makes it trivial to modify the thermochemistry solver, i.e. change it from equilibrium to non-equilibrium.
The basic grid element in RAMSES is an oct (Fig. 4), which is a grid composed of eight cubical cells. A conservative state vector U = (ρ, ρu, E, ρZ) is associated with each cell storing its hydrodynamical properties of gas density ρ, momentum density ρu, total energy density E and metal mass density ρZ. (One can also use the primitive state vector, defined as W = (ρ, u, P, Z).) Each cell in the oct can be recursively refined to contain sub-octs, up to a maximum level of refinement. The whole RAMSES simulation box is one oct at = 1, which is homogeneously and recursively refined to a minimum refinement level min, such that the coarse (minimum) box resolution is 2 min cells on each side. Octs at or above level min are then adaptively refined during the simulation run, to follow the formation and evolution of structures, up to a maximum refinement level max, giving the box a maximum effective resolution of 2 max cell widths per box width. The cell refinement is gradual : the resolution must never change by more than one level across cell boundaries.

RAMSES multi-stepping approach
With AMR multi-stepping, the resolution is not only adaptive in terms of volume, but also in time, with different timestep sizes on different refinement levels. A coarse timestep, over the whole AMR grid, is initiated at the coarse level, min, as we show schematically in Fig. 5. First, the coarse time-step length ∆t min is estimated via (the minimum of) Courant conditions in all min cells. Before the coarse step is executed, the next finer level, min + 1, is made to execute the same time-step, in two substeps since the finer level Courant condition should approximately halve the time-step length. This process is recursive: the next finer level makes its own time-step estimate (Courant condition, but also ∆t ∆t −1 ) and has its next finer level to execute two substeps. This recursive call up the level hierarchy continues to the highest available level max, which contains only leaf cells and no sub-octs. Here the first two substeps are finally executed, with step lengths ∆t max ∆t min /2 max− min . When the two max substeps are done, the max − 1 time-step is re-evaluated to be no longer than the sum of the two substeps just executed at max, and then one max − 1 step is executed. Then back to level max to execute two steps, and so on. The substepping continues in this fashion across the level hierarchy, ending with one timestep for the coarsest level cells (with a modified time-step length ∆t min ).
At the heart of RAMSES lies a recursive routine called amr step( ) which describes a single time-step at level , and is initially called from the coarsest level ( min). To facilitate our descriptions on how the RT implementation is placed into RAMSES, we illustrate the routine in pseudocode format in Listing 1, where we have excluded details and bits not directly relevant to RHD (e.g. MPI syncing and loadbalancing, adaptive refinement and de-refinement, particle propagation, gravity solver, star formation, and stellar feedback).
First, the recursion is made twice, solving the hydrodynamics over two sub-steps at all finer levels. Then the Euler equations are solved over the current coarse time-step, for all cells belonging to the current level. It is important to note here that the hydrodynamical quantities are fully updated at the current level in the hydro solver, but there are also intermediate hydro updates in all neighboring cells at the next coarser level. The coarser level update is only partial though, because it only reflects the intercell fluxes across inter-level boundaries, and fluxes across other boundaries (same level or next coarser level) will only be accounted for when the coarser level time-step is fully advanced. Until Having put down the basics of AMR hydrodynamics, we are now in a position to add radiative transfer.

RAMSES-RT
In RAMSES-RT, each cell stores some additional state variables. Here U = (ρ, ρu, E, ρZ, ρxHII, ρxHeII, ρxHeIII, Ni, Fi), where xHII, xHeII and xHeIII are the hydrogen and helium ionization fractions, which are advected with the gas as passive scalars (in the hydro solver), and Ni, Fi represent the 4M variables of photon density and flux for each of the M photon groups. Note that this represents a hefty increase in the memory requirement compared to the hydrodynamics only of RAMSES: the memory requirement for storing U (which is the bulk of the total memory in most simulations) is increased by a factor of 1.5(1+4/9M ), where the 1.5 represents the ionization fractions and the parenthesis term represents the photon fluxes and densities. Thus, with three photon groups, the memory requirement is increased by roughly a factor 3.5 compared to a traditional RAMSES simulation. t + ∆t rt min t Figure 6. Diagram of the amr step in RAMSES-RT. This is much like the normal amr step in RAMSES, except that the time-step length has the extra constraint of the light speed Courant condition, and each level step also performs photon injection, RT transport and thermochemistry over the same time-step and level.
Listing 2: The AMR step in RAMSES-RT. Given the time-scale difference between hydrodynamics and radiative transfer, the obvious approach to performing RHD is to sub-cycle the three radiative transfer steps (injection, advection, thermochemistry) within the hydrodynamical step. There is, however, a major drawback to this approach, which is that it is incompatible with AMR multistepping: the RT sub-cycling must be done before/after each hydrodynamical AMR step at the finest refinement level only, and since light can in principle cross the whole box within the fine level hydrodynamical timestep, the RT subcycling must be done over the whole grid, over all levels. However, the partial hydrodynamical flux between cells at level boundaries always leaves some cells between the fine level steps with an intermediate (i.e. partially updated) gas state. This makes the thermochemistry ill-defined in those cells, since it needs to update the gas temperature in every cell, and for this to work the temperature must have a well defined and unique value everywhere. There are three ways around this: First is to perform the RT subcycling only after a coarse hydrodynamical step, but here potentially thousands of finescale hydro steps would be executed without taking into account the thermochemistry.
Second, is to prohibit AMR multi-stepping, which makes the whole grid well defined after each step and thus allows for RT sub-cycling over the whole box. Multi-stepping is however one of the main advantages of AMR, and essen-tially allows us to refine in time as well as space, so this isn't really an option.
We thus default to the third strategy, which we use in RAMSES-RT. Here we drop the subcycling of RT within the hydro step and perform the two on the same timestep length, which is the minimum of the RT and hydro timestep. Thus, with each hydro step, at any level, the RT steps are performed over the same level only. The basic scheme is illustrated in Fig. 6, and the pseudocode for the updated amr step is shown in Listing 2. Obviously, the main drawback here is the timescale difference, which can be something like a factor of 100 − 1000, meaning the number of hydrodynamical steps is increased the same factor and the run-time accordingly (plus numerical diffusion likely becomes a problem with such small hydrodynamical time-steps). However, if we also apply a reduced speed of light, we can shrink this factor arbitrarily, down to the limit where the hydrotimestep is the limiting factor and the only increase in computational load is the added advection of photons (which is considerably cheaper for one photon group than the hydrodynamical solver) and the non-equilibrium thermochemistry (which typically has a computational cost comparable to the equilibrium solver of RAMSES, provided we use RT smoothing). The question, which we have tried to answer in §4.2, is then how far we are allowed to go in reducing the light speed.
Parallelization is naturally acquired in RAMSES-RT by simply taking advantage of the MPI strategies already in place in RAMSES.

Radiation transport on an AMR grid
In RAMSES-RT, the radiation variables are fully incorporated into the AMR structure of RAMSES. The ionization fractions and photon densities and fluxes are refined and de-refined along with the usual hydro quantities, with a choice of interpolation schemes for newly refined cells (straight injection or linear interpolation). The radiative transfer, i.e. injection, transport and thermochemistry, is multi-stepped across the level hierarchy, thus giving AMR refinement both in space and time. Inter-level radiation transport is tackled in the same way as the hydrodynamical advection, i.e. transport on level includes partial updates of neighbouring cells on level − 1. Update of the finer level cell RT variables over level boundaries involves the RT variables in a coarser cell, which are evaluated, again with the same choice of interpolation schemes. RAMSES-RT includes optional refinement criteria on photon densities, ion abundances and gradients in those, in addition to the usual refinement criteria that can be used in RAMSES (on mass and gradients in the hydrodynamical quantities).
Of the seven standard RT and RHD tests described in Section 6, five include active or inactive grid refinement, demonstrating that the radiation hydrodynamics perform robustly in conjunction with (on-the-fly) cell refinements/de-refinements. In addition, we demonstrate in Fig. 7 how radiation flux is well retained across changes in grid refinement. The upper left map of the figure shows a beam of radiation in a 2D RAMSES-RT experiment, where we use the HLL flux function and deactivate radiation-gas interactions (with zero photoionization cross sections). The beam is injected into two cells in the bottom left corner by impos-ing a unity reduced photon flux of 3×10 10 photons s −1 cm −1 , corresponding to a photon density of 1 cm −2 , at an angle of 26.5 • from the horizontal. The beam traverses a circular region of 2 successive levels of increasing refinement, going from refinement level 6 to 8, i.e. effective resolutions of 64 2 to 256 2 cells. We use here straight injection (i.e. no interpolation) for inter-level cell fluxes, but linear interpolation gives identical results. The snapshot is taken at t = 3.04 × 10 −11 s, just before the beam has had time to cross to the right edge of the 1 cm wide box. To the right of the map we plot photon flux profiles, cN , across the coloured vertical lines in the map. The beam experiences diffusion, as can be seen by the widening of the flux profiles, but this is exclusively due to the intercell flux function and independent of the refinement changes. The far left plot shows the integrals across each flux profile, i.e. the total photon flux across each line. The values are consistent until around x = 0.6, and then reduce to zero towards the edge of the beam. We've verified that if the test is let to run for double the time, i.e. about 6 × 10 −11 s, the total flux is consistent throughout the whole box width to about 1 in 10 4 , so photons are very well conserved across the changes in refinement.
To further demonstrate flux conservation, the lower panel in the same figure shows an identical experiment except that the beam is horizontal, such that it can be perfectly maintained with the HLL flux function. To stay just under a light crossing time, we consider a shapshot at 2.6 × 10 −11 s. Here again, the flux is well preserved towards the edge of the beam, and we have verified that in two crossing times, the total flux is retained perfectly to the number precision, which here is 7 decimals.
We also consider another beam with the same setup, shown in Fig. 8, where instead of a static refinement region, the grid is actively refined on inter-cell gradients in photon density N . According to the criterion, two adjacent cells at positions i and i + 1 are refined if Straight injection (no interpolation) is used here for interlevel fluxes and cell refinements, but the results are identical when linear interpolation is used for inter-level fluxes and cell refinements. The snapshot here is taken at 3.3 × 10 −11 s (∼ a crossing time). The plot on the far right shows the flux conservation across different x-coordinates. (Note the total flux is slightly different from that in Fig. 7 because of the different geometry of the beam injection.) The total flux is again well maintained towards the beam edge. We verified that in two light crossing times, the discrepancy of the beam flux at different x−coordinates levels out to within 0.03%. These simple beam experiments demonstrate that the code accurately transports radiation across (even dynamically) changing refinement levels. The main errors are the artificial diffusion of radiation on the grid, which is not caused by refinement, but rather by the inter-cell flux function, and the dipole approximation inherent to the M1 closure, which does not allow opposing streams of radiation to pass through one another. Note though that while the diffusion is artificial, the total flux is well maintained, i.e. energy is conserved.

Speed of light
The AMR transport tests also demonstrate that radiation in RAMSES-RT propagates at the correct speed, i.e. at the speed of light. In each beam map (Figs. 7-8), a black line has been plotted over the beam, starting at the beam injection and ending at the light-crossing distance, i.e. t × c, where t is the snapshot time. Qualitatively it can be seen that the beam ends roughly at the same position as the black line, and in the flux plots on the far right side of each beam map it can be seen that the beam has roughly half the original flux at this end position. The far end of the beam is smooth over a few cell widths rather than discontinuous, because of numerical diffusion.

Cosmological settings
RAMSES uses super-comoving variables to allow for the impact of the cosmological expansion on the Poisson equation, the equations of hydrodynamics (38-40) and particle propagation (Martel & Shapiro 1998;Teyssier 2002): a change is made from the physical variables to super-comoving ones with u,ε = a 5 ΩmρcH 4 0 L 2 ε, where H0 is the Hubble constant, Ωm the matter density parameter, L the comoving width of the simulation box (physical width at a = 1), and ρc the critical density of the Universe. When these variables are used instead of the physical ones, the cosmological expansion is accounted for, while all relevant equations remain unchanged, Euler equations included.
For consistency, and to partly account for the effect of cosmological expansion on the radiative transfer, the additional change is made in RAMSES-RT to super-comoving RT variables for the photon transport: c.
The dilution (∝ a −3 ) of photon number density is thus accounted for, while it can easily be verified that Eqs. (6)- (7) remain unchanged with the new variables -including the reduced flux (14) used in the M1 tensor (12). Note that when reduced light speed is used, the photons will be over-diluted in cosmological simulations, since the time taken for them to get from source to destination will be overestimated. Note also that wavelength stretching with redshift, which in reality adds a fourth power of a to the dilution of Nγ, is not accounted for here. This is actually nontrivial to do: one could add one power of a to the definitions ofÑ andF, but it would be a very crude approximation of the wavelength dilution, as the wavelength shift that should feed photons from one group to the next is neglected. In any case, this effect is likely to be important only in the context of reionization, where the photons have a chance of travelling cosmological distances before they are absorbed. While cosmological diffusion and redshifting is difficult to account for in ray-tracing methods, where the radiation is typically traced as far as it can get in one moment in time, momentbased approach are more straightforwardly able to model these effects (e.g. Ricotti, Gnedin & Shull 2002;Petkova & Springel 2009;Finlator, Dave & Ozel 2011).

RADIATIVE TRANSFER TESTS
The tests described in this section come from two papers that were born out of a series of workshops on radiative transfer. Tests with simple analytic results to compare to are hard to engineer in radiative transfer, so the solution was to instead make simple tests where the correct result is not necessarily well known but the results of many different codes can instead be compared. Thus it is likeliest that the correct results are usually where most of the codes agree, and if a code stands out from all or most of the others in some way, this would most likely be a problem with that particular code. These tests have become sort of benchmark tests for RT codes, and most publications that present new implementations use some or all of these tests for validation.
The first paper is Iliev et al. (2006a), hereafter known as Il06 -it describes four RT post-processing tests, i.e. with the hydrodynamic advection turned off, and shows the results for 11 RT codes. The second paper is Iliev et al. (2009), hereafter known as Il09 -it describes three additional tests, and results for 9 codes, where the RT is coupled to the hydrodynamics.
The tests results from Il06 and Il09 are normally downloadable on the web, but at the time of this writing the links have been down for some time. However, Ilian Iliev has been kind enough to provide all test results for one of the codes, the grid based short characteristics code C 2 -Ray, which is described in detail in Mellema et al. (2006). We thus present here RAMSES-RT results with comparisons to those of C 2 -Ray. The inclusion of the C 2 -Ray results in the plots shown here should be useful to guide the eye if one then wants to compare with the other codes in Il06 and Il09.
As prescribed by the test papers, all tests use hydrogen only gas. We use smooth RT in the RAMSES-RT runs for all tests, but remark that turning off the smoothing has no discernible effect on the results (only calculation speed). Unless noted otherwise in the following tests, the GLF intercell flux function is used ( §3.2), and the the on-the-spot approximation is applied ( §3.3.2). In all except test 1, where the radiation is monochromatic, the radiation energy distribution is assumed to be a T eff = 10 5 K blackbody, which is approximated with three photon groups bordered by the hydrogen and helium ionization energies: A reduced speed of light fraction of fc = 1/100 is used unless otherwise noted. AT08 contain an analysis of the effect of different light speeds in the first three tests from Il06, and find the results start diverging non-negligibly somewhere between fc = 10 −2 and 10 −3 , which matches well with our analysis in §4.2. The prescribed resolution in the tests is 128 3 cells, but in most tests we use adaptive refimenent for demonstrative purposes, with a coarse resolution of 64 3 cells, and an effective resolution of 128 3 cells. We use a Courant factor of 0.8, so the RT timestep is set by where ∆x is the cell width and cr the reduced light speed.
Taking as an example the test 1 setup, which has a box width of 6.6 kpc, a simulation time of 500 Myr, and a reduced light speed fraction fc = 10 −2 , this translates into a (fine level) timestep length of ∼ 4, 500 yr, so ∼ 10 5 fine-level steps need to be computed to run the test.

Il06 test 0: The basic thermochemistry physics
This is essentially a one-cell test of the non-equilibrium thermochemistry and not radiative transfer per se, so it doesn't really count with the rest of the comparison project tests (hence test zero). It is important nontheless since thermochemistry is a major new component in RAMSES-RT.
We start with completely neutral hydrogen gas with density nH = 1 cm −3 and temperature T = 100 K at t = 0. A photo-ionizing flux of F = 10 12 s −1 cm −2 with a 10 5 K blackbody spectrum is applied to the gas and maintained until t = 0.5 Myr at which point it is switched off. The run is continued for a further 5 Myr, allowing the gas to cool down and recombine. The run-time is separated into 500 logarithmically equally spaced timesteps, and the thermochemistry solver sub-cycles these timesteps adaptively (see §3.3.1). The photon flux is not evolved, i.e. it is kept fixed (until 0.5 Myr) thoughout the integration. The resulting evolution of the neutral fraction and temperature of the gas is shown in Fig. 9. The evolution closely follows that of the codes described in Il06, with the exception of SimpleX and FFTE which stand out somewhat, and we don't see any sign of the stiffness-induced oscillations that can be seen in the Crash code test.

HI fraction
10 −3 10 −2 10 −1 10 0 500 Myr Figure 10. Il06 test 1. Map of the neutral fraction in a box slice at z = 0, at 500 Myr. Overplotted is the AMR grid, which is refined on the fly during the experiment from 64 3 to 128 3 cells effective resolution. Maximum refinement stays on the corner source throughout the run, and it adaptively follows the expansion of the I-front.

Il06 test 1: Pure hydrogen isothermal HII region expansion
A steady monochromatic (hν = 13.6 eV) source of radiation is turned on in a homogeneous neutral gas medium, and we follow the resulting expansion of a so-called Strömgren sphere of ionized gas. Heating and cooling is turned off and the temperature is set to stay fixed at T = 10 4 K. The box is a cube of width L box = 6.6 kpc. The gas density is nH = 10 −3 cm −3 and the initial ionization fraction is xHI = 1.2 × 10 −3 , corresponding to collisional ionization equilibrium. The radiative source is in the corner of the box and the emission rate isṄγ = 5 × 10 48 photons s −1 . The simulation time is tsim = 500 Myr. To demonstrate on-thefly AMR at work (and speed up the runtime), we use a base resolution of 64 3 cells, but allow for one level of further refinement, i.e. to the effective prescribed resolution of 128 3 cells. Typically, AMR refinement is applied on mass-related criteria, since massive structures are usually the objects of interest in simulations. However, since the density field is homogeneous in this test, we apply refinement on gradients in xHI and xHII: two adjacent cells at positions i and i + 1 are refined if 2 x where x is either xHI or xHII. The Strömgren radius, rS, is the radius of the ionization front (I-front) from the center when steady state has been reached, and in the case of fixed density and temperature it has the simple analytical result shown in Eq. 32. In this result the I-front evolves in time according to where trec = (nHα B HII ) −1 is the recombination time. For the parameters of this experiment, trec = 122.4 Myr and rS = 5.4 kpc. Fig. 10 shows maps at 500 Myr of the neutral fraction, with the grid refinement overplotted, in a box slice at z = 0. The Strömgren sphere is nicely symmetric and qualitatively it can be seen to agree well with results from the RT codes described in Il06 (their Fig. 6). Fig. 11a shows the evolution of the I-front position and velocity with RAMSES-RT (solid blue), compared with the analytic expression (green dot-dashed) and the result for the C 2 -Ray code (red dashed), which is typical for the RT code results presented in Il06 and does not stand out particularly in this test. Our result can be seen to match the C 2 -Ray one, though we have an initial lag due to the reduced speed of light that can best be seen in the top plot showing the fraction of the numerical result's I-front radius versus rS. The analytic rI is typically ahead of rS by 5%, which is simply because the analytic result is step-like with complete ionization within rS and none outside, whereas the real result has a gradually evolving ionization profile with radius. Indeed, Pawlik & Schaye (2008) computed the exact analytic result to this problem, accounting for an equilibrium neutral fraction inside the Strömgren sphere, and found an equilibrium I-front radius which is exactly 1.05 rS. Fig. 11b shows spherically averaged radial profiles of the gas ionization state at 30 and 500 Myr. Again we see a good match with the C 2 -Ray result. There is still a little lag in the I-front position at 30 Myr due to the RSLA and xHI is somewhat lower inside the Strömgren sphere in RAMSES-RT. However, the C 2 -Ray result stands out a little in this test in Il06 as being most effective at ionizing the gas within the Strömgren sphere (i.e. has the lowest values of xHI), and the RAMSES-RT result is typical of the Il06 codes' results in this plot.
A further comparison is made in Fig. 11c, here comparing ionization fraction histograms at three simulation times. Again the RAMSES-RT result closely matches the C 2 -Ray one, whose histograms fall into a group with the codes IFT, Flash-HC and FFTE that stand out a little in Il06 (Fig. 9) as having less frequent intermediate neutral fractions than the other codes.
Finally for this test, Fig. 12 shows a comparison with C 2 -Ray of the globally averaged neutral fraction as a function of time. It is a close match, and the C 2 -Ray result is here typical for the Il06 codes. All in all, there is nothing out of the ordinary in the RAMSES-RT result for Il06 test 1, except for a slight initial delay of the I-front which is to be expected due to the RSLA.
We note that performing this test with the full prescribed 128 3 resolution, rather than using AMR like we've done here, has no discernible effect on the results. In the AMR run, the number of fine level cells is maximally (at the end of the run) 15% of the number of fine level cells in the non-AMR run, and the computation time is 30% of that in the non-AMR run. The cost of the experiment (with AMR) is on the order of 50 cpu hours 6 , which is a lot for a simple test in which little actually happens: for much of the run, the I-front is moving towards a stand-still at speeds which are much slower than our reduced speed of light (fc = 0.01), so barring the RT Courant condition, the timesteps taken could have dramatically increasing length towards the end of the test. Implicit transport solvers can take advantage of this (almost) static situation by on-the-fly adapting the timestep length (which is in the case of implicit solvers not constrained by the Courant condition), so presumably an implicit solver can run this test (and most of the tests described in this work) with considerably less computation than we do. However, in more realistic cosmological scenarios, such steady regimes simply do not happen over times longer than the typical age of stellar populations, which is on the order of 10 Myr (50 times shorter than the run-time for this test). Furthermore, stellar populations typically are turning on and off on even shorter timescales than that thoughout the simulation volume, which limits the dynamical time of ionization fronts even further. This presumably constrains the main advantage (possible long time steps) of implicit solvers severely, since even though they are not constrained by Courant-like conditions, they still need to resolve dynamical timescales.

Il06 test 2: HII region expansion and the temperature state
The setup here is the same as in Il06 test 1, except for the following points: • We allow for cooling and photo-heating of the gas, i.e. the temperature is no longer constant, and the analytic result, Eq. 32 no longer applies (because of the non-constant recombination rate).
• The initial temperature is 100 K.
• The initial ionization fraction of the gas is xHII = 10 −6 . It should be fully neutral according to the test recipe in Il06, but this is (the default) minimum value for xHII in RAMSES-RT, that exists in order to keep bounds on the subcycling of the thermochemistry. In any case, the specific value is not critical to the test results, as long as it is low.
• The radiation source is a T = 10 5 K blackbody, modeled with the three photon groups defined by (44). The emission rate is the same as before,Ṅγ = 5 × 10 48 photons s −1 .
• We don't use grid refinement in this test. The grid is homogeneous and the resolution is 128 3 grid cells, as prescribed in Il06.
Slice maps at z = 0 of the neutral fraction and temperature are shown in Fig. 13. Both the ionization and heating fronts are smooth and symmetric, and the maps agree qualitatively with other codes in Il06 (Figs. 11-14). In comparison with the same test with ATON (AT08, Fig. 3), both fronts are clearly much thicker here, which is due to our multifrequency implementation (whereas ATON used one photon group). More detailed comparison with the Il06 codes can be made through the ionization state and temperature plots in Fig. 14a, where we include the C 2 -Ray result. The ionization state profile develops very similarly to that of C 2 -Ray, though we have less ionization on both sides of the front, especially on the outer side where the difference in xHII is as high as a factor of ten. Presumably this is due to the different implementations of multi-frequency photo-heating and cooling. The thermal profiles are also similar to C 2 -Ray, though we have considerably lower (up to a factor of two) temperatures on the inside of the I-front, and conversely higher temperatures on the outside. As can be seen in Fig.  17 in Il06, C 2 -Ray has the strongest heating of any code on the inside of the I-front in this test and most codes have stronger heating on the outside, so our thermal profiles (as the ionization state profiles) are fairly typical of the ones presented in Il06 for this test. Fig. 14b shows the evolution with time of the ionization front, compared with C 2 -Ray and the analytic result from test 1. The front moves more slowly here than in test 1 due to the lower initial temperature, so we no longer lag behind in the initial front propagation. Our front propagates slightly further than in C 2 -Ray, and ends at almost exactly the same radius as the FFTE code, which has the furthest expanding I-front of any code in this test in Il06. Still the difference between the codes is small, with the ratio between the numerical and analytic results (rnum/r analyt ) ranging between 1.01 and 1.11. Fig. 14c shows histograms of the ionized fraction and temperature at different times in the test for RAMSES-RT and C 2 -Ray. The ionized fraction histograms are quite similar, the biggest difference being a higher fraction of almost completely neutral gas xHII 10 −2 in RAMSES-RT, which we already saw in Fig. 14a (top) beyond the I-front. The temperature histogram for RAMSES-RT differs a bit from C 2 -Ray in having less extreme temperatures (C 2 -Ray has both hotter gas and colder gas) but are very similar to those for the codes ART, RSPH and Crash in Il06.
Finally, Fig. 15 shows the time evolution of the volume averaged neutral fraction in RAMSES-RT and C 2 -Ray, and here we see a close match. There is quite a lot of discrepancy between the different codes in the analogue plot in Il06 (Fig.  20), with 3 groups of results, and our result closely follows those of C 2 -Ray, Crash and RSPH.  As with test 1, there is nothing out of the ordinary in the RAMSES-RT result for Il06 test 2, except perhaps for an ever so slightly further advanced I-front than most codes in Il06 have.

Il06 test 3: I-front trapping in a dense clump and the formation of a shadow
This test considers self-shielding within a dense gas cloud bombarded on one side by UV radiation, and the shadow trailing on the 'dark' side -something which may find place with clouds close to sites of star-formation. The setup is as follows: the simulation box has width L box = 6.6 kpc. We place a spherical cloud of gas in the center of the (y, z)-plane, with radius r cloud = 0.8 kpc, and it's center at (xc, yc, zc) = (5, 3.3, 3.3), as seen in Fig. 16, top left, showing an (x, y)−slice of gas density through the middle of the box. Outside the gas cloud we have n out H = 2 × 10 −4 cm −3 , T out = 8000 K and x out HII = 0, and inside we have n cloud H = 200 n out H = 4 × 10 −2 cm −3 , T cloud = 40 K and x cloud HII = 10 −6 . We apply a constant ionizing photon flux F = 10 6 s −1 cm −2 from the x = 0 boundary of the box (left in the Fig. 16 maps), and run for 15 Myr. We use a light speed fraction of fc = 10 −1 . This is ten times higher than the "norm" in the RT tests, but it is needed for the light to have reached the cloud in the first snapshot under consideration, at 1 Myr. In order to best capture the formation of a shadow behind the cloud, we apply the HLL flux function in this test rather than the usual GLF function, and we use the OTSA. We have run identical tests though, one with the GLF flux function, and one where we use the HLL flux function but don't assume the OTSA, and we show maps of those experiments for a qualitative comparison. As usual, the resolution prescribed by Il06 is 128 3 cells, but here we apply static AMR refinement such that the coarse resolution is 64 3 cells, but a rectangular region that encompasses the gas cloud and the shadow behind it has one level of additional refinement, making the effective resolution in the cloud and its shadow 128 3 cells. The refinement region is shown in the top panel of Fig. 16, plotted over a density map that shows the spherical gas cloud. The fraction of volume at the fine resolution is 4%, and the computation time for the test is roughly a quarter of a an analogous uniform grid run (about 32/130 cpu hours for the AMR/non-AMR runs). Fig. 16 shows slices at z = 0.5 L box of the neutral fraction and temperature at 1 and 15 Myr. From second top to bottom row are shown RAMSES-RT+HLL, RAMSES-RT+HLL without the OTSA, RAMSES-RT+GLF (with the OTSA) and C 2 -Ray. The I-front travels fast through the diffuse medium outside the cloud, but moves much more slowly inside it, and a shadow is cast behind it. As the UV radiation eats its way into the cloud, ionizing and heating it, the shadow also very slowly diminishes in width because some photons manage to cross through the edges of the cloud. The RAMSES-RT+HLL maps compare very well with C 2 -Ray, though the shadow is slightly thinner at 15 Myr and there is stronger heating inside the shadow; this could be due to differences in the multifrequency approach and/or photoheating. Without the OTSA, the shadow is diminished from the sides due to photons being cast from the surrounding gas. Using the GLF flux function has much the same effect as not assuming the OTSA, though the shadow is considerably more diminished here. The result with HLL but without the OTSA is the most physical of the RAMSES-RT results, as one should expect recombination photons to be cast into the shadow. Fig. 17a shows the evolution of the position and speed of the I-front through the center of the (y, z)−plane. In solid blue we plot the RAMSES-RT result and in dashed red is the C 2 -Ray result for comparison. Horizontal dotted lines mark the edges of the cloud. There is a large initial delay in the I-front compared to C 2 -Ray, which is because in the diffuse gas outside the cloud, the I-front speed is limited by the reduced speed of light. After the I-front gets into the cloud (lower dotted line) it quickly catches up and then evolves in a similar fashion in the two codes. If compared to the rest of the codes in Il06, it turns out that the evolution of the I-front in C 2 -Ray slightly stands out from the rest of the codes (e.g. a small upwards 'bump' in the front position at log(t/trec) ∼ 0, and a slightly shorter distance of the Ifront from the origin at the end of the simulations), and most of the others in fact evolve very similarly to that of RAMSES-RT. The comparison appears best with RSPH, which has the furthest extended I-front at the end-time of 15 Myr. The same can be said for the speed of the front. If we look away from the initial ∼ 0.2 Myr, when our I-front has to catch up, the speed compares reasonably to C 2 -Ray, and quite well to the other codes in Il06. Fig. 17b shows the evolution of the mean ionized fraction and temperature inside the cloud, compared between RAMSES-RT and C 2 -Ray. The evolution is similar between the two codes in both cases. Compared with the other codes in Il06, the evolution of the ionized fraction is most similar to RSPH, IFT and Coral, while the temperature in RAMSES-RT is consistently a little higher than in most codes (all except Coral and Flash which stand out quite a lot in mean temperature). Fig. 17c shows profiles of the ionization state and temperature along the x-axis at the center of the (y, z)−plane at 1, 3 and 15 Myr. The ionization state profile in RAMSES-RT is similar in most respects to that of C 2 -Ray, though it extends a bit further at the end of the run-time. There is initially less ionization on the far side of the front in RAMSES-RT, but at the end of the run this is reversed and we have slightly more ionization on the far side in RAMSES-RT. This 'shift' can be explained by the temperature profiles: at early times the cloud is efficiently shielding the far side from even the high-energy photons in both codes, but at the end of the RAMSES-RT run the shielding buffer in the cloud is thin enough that the high-energy photons can get through, hence efficiently heating the gas inside the buffer as well as in the shadow, and the gas in the shadow becomes slightly ionized as a consequence. The analogue ionization state profiles for the other codes in Il06 are mostly similar to ours. Most of them are actually closer to the RAMSES-RT than the C 2 -Ray profile, with the exception of Crash which has a much more underdeveloped I-front and less ionization, and FFTE and IFT which have an almost step-wise xHII-profile on the far side of the I-front. The temperature profiles differ pretty widely between the codes. RAMSES-RT doesn't particularly stand out, though, and is most similar to that of Coral at 15 Myr. The temperature profile for RAMSES-RT also differs notably from that of ATON, where the shielded region inside the cloud is thicker and more step-like both in the ionized fraction and temperature, due to the monochromatic radiation.
Finally, Fig. 18 shows histograms of the neutral fraction and temperature at 1, 3 and 15 Myr for RAMSES-RT and C 2 -Ray. The comparison (also with the other codes in Il06) is qualitatively similar, though there is quite a difference between the individual codes in these plots.
As with the previous tests, RAMSES-RT performs well here and we don't really have anything out of the ordinary in our results. One should keep in note though that here we've used the non-diffusive HLL flux function, whereas in most cosmological simulations it would be more natural to use the more diffusive GLF function to have better spherical symmetry around radiative stellar sources, which comes with the price of less pronounced and shorter lived shadows than HLL. The survival of shadows in more realistic scenarios remains an open question, but considering the effects of recombination radiation, and the likelihood of any transparent region to have ionizing sources shining from different directions, it seems unlikely to us that shadowing is an efficient way of shielding gas from ionizing radiation.

Il06 test 4: Multiple sources in a cosmological density field
This test involves the propagation of ionization fronts in a static hydrogen-only density field taken from a cosmological simulation snapshot at redshift 9. The density cube is 128 3 cells and its width is 500 h −1 co-moving kpc (corresponding to 50 h −1 physical kpc). The Hubble factor is h = 0.7. The initial temperature is fixed at 100 K everywhere. 16 radiative sources are picked out corresponding to the most massive halos in the box and these are set to radiate continuously for 0.4 Myrs. The mass-dependent radiation intensity for each halo is given in a downloadable table (from the RT comparison project website). Unlike in Il06, we don't apply the OTSA in this test, i.e. we include the radiative transfer of recombination radiation, but we've verified that this has no discernible effect on the results. Our analysis from §4.2 indicates that a reduced light speed gives incorrect results in this test. Thus we use a full light speed here (i.e. fc = 1), and for comparison with the codes from Il06, which implicitly assume infinite light speed, we make an analogue run with a hundred-fold light speed (fc = 100). Fig. 19 shows box slices, at z = 0.5 L box , of the neutral fraction and temperature at times 0.05 and 0.4 Myr. Shown are our two runs with different light speed fractions (top and bottom row), and for comparison we show the result for the C 2 -Ray code, from Il06 7 : the I-fronts and photo-heating in our fc = 1 run clearly lag behind the C 2 -Ray result, and there is also less heating of the ionized gas. This is in accordance with the ATON results described in AT08, where a similar delay was found. They prescribed this delay to the fact that ATON is monochromatic, but since our multi-frequency approximation (three photon groups) gives results that are still much more similar to the ATON results than those of C 2 -Ray, especially in terms of the neutral fraction maps, we are inclined to blame the delay on another factor, which is the speed of light. Our results with the speed of light set to one-hundred times the physical value are shown in the bottom row of Fig. 19 and here the results are considerably closer to those of C 2 -Ray in terms of the propagation of heating-and I-fronts, although the maximum temperature in the ionized gas is still colder in comparison. All four codes considered in the Il06 4 test use an infinite effective speed of light and this may give premature fronts in the immediate vicinity of the sources and also further away in under-dense regions. Thus we are perhaps not really dealing with a delay in RAMSES-RT, but rather premature fronts in the Il06 codes. As AT08 note, we are far from reaching a static state in the fronts in this experiment in the run-time of 0.4 Myr and we should expect the different light speed runs to converge to similar results when static state is reached. This is further corroborated by our I-front light crossing time analysis from Section 4.2.
The smaller degree of photo-heating in the ionized gas compared to the C 2 -Ray results is in line with the temperature profiles from the previous tests (e.g. Fig. 14a), and presumably is a consequence of the different ways multifrequency is approximated. Another notable difference in the maps in Fig. 19 is that our fronts are smoother and less jagged than those in C 2 -Ray. This is an effect of the photon diffusion inherent in the GLF flux function used here. Like AT08 we find that using HLL instead gives more jagged fronts. Fig. 20a shows the evolution of the mass-and volumeweighted ionized fractions, compared for the different runs. The RAMSES-RT run with the physical light speed gives ionized fractions which are close (both mass-and volumeweighted) to the ATON ones, whereas increasing the light speed by a factor of hundred from the physical value gives results closer to C 2 -Ray (as well as the three other codes that ran this test in Il06). Presumably we would converge further towards C 2 -Ray in the limit of infinite light speed, but computational time constraints do not allow to pursue that investigation. This is a further hint that the correct speed of light is important in the non-steady regime of ionization fronts.
Finally, Fig. 20b shows neutral fraction and temperature histograms at three times in the test. Again there is a strong discrepancy between the RAMSES-RT run with fc = 1 and C 2 -Ray, especially at early times, and the gap all but closes when fc = 100 is used instead with RAMSES-RT. There remains some difference though in the minimum/maximum temperature, being smaller/larger for C 2 -Ray than for our fc = 100 run, presumably because of our rather crude multifrequency approximation.
To summarize, there is notable discrepancy between the RAMSES-RT results and those presented in Il06, in that the RAMSES-RT ionization front lags behind, which appears to be due to a finite speed of light. This is corroborated to some degree by other papers in the literature: Wise  cally do a comparison between finite and infinite light speed, with the finite one resulting in a delay which is substantial, though maybe a bit less than ours, and they do comment on the ionization bubbles in this test being unphysically large with infinite light speed methods. Other sources appear to be in conflict with our conclusions, though: Petkova & Springel (2011a) use a finite light speed and get results which seem to compare well with those of C 2 -Ray.

Il09 test 5: Classical HII region expansion
We now come to the tests described in second radiative transfer codes comparison paper by Iliev et al. (2009), which we denote as Il09. This paper provides 3 code comparison tests to add to those in Il06, but with the important difference that whereas the Il06 tests are pure radiative transfer post-processing tests with fixed density fields, the tests in Il09 are RHD tests, i.e. with the radiative transfer directly coupled to the gas-dynamics. Thus we now switch from the context of post-processing RT to hydro-coupled RHD. Here, the pressure buildup in photo-heated gas causes it to expand. Typically, the I-front is initially R-type, where it expands much faster than the gas response to it, which means RT-postprocessing is a fairly good approximation. The I-front then begins to slow down when it approaches the Strömgren radius, but gets moving again when the gas catches up to it, and then the front is D-type, i.e. moves along with the expanding gas. As before we compare our RAMSES-RT tests results with those of the grid-based short characteristics ray-tracing code C 2 -Ray , here coupled to the Capreole code, which employs a Riemann solver for the hydrodynamics. As the Capreole+C 2 -Ray combination is sensitive to numerical instabilities appearing in Il09 test 6, we compare also in that particular test to C 2 -Ray coupled to the Eule-rian TVD solver of Trac & Pen (2004) (that combination was not used in any other tests). The test numbers continue from the Il06 paper, thus we now come to Il09 test 5, which concerns the expansion of an ionization front due to a point source in an initially uniform-density medium. The initial setup, much like that of Il06 test 2, is as follows.
The box cube is L box = 15 kpc in width. The gas is hydrogen only as usual, initially homogeneous with density nH = 10 −3 cm −3 , temperature 100 K, and ionization fraction xHI = 10 −6 (Il09 prescribes xHI = 0). The radiative source is in the corner of the box and the emission rate iṡ Nγ = 5 × 10 48 photons s −1 . We don't apply the OTSA in this test, i.e photons are emitted from gas recombinations. The simulation time is 500 Myr. The base resolution of the box is 64 3 cells and we apply on-the-fly refinement on nH and xHII gradients (see Eq. 46), so that the ionization front has the prescribed effective resolution of 128 3 cells.
We first compare volume dissections at z = 0 in the simulation cubes at 100 and 500 Myr, for the RAMSES-RT and C 2 -Ray results, shown in Fig. 21. The maps show, from left to right, the neutral fraction, pressure, temperature, density and mach number, M ≡ v/cS, where cs = 1.4 P/ρ is the sound speed. (Unfortunately the M output is missing from the C 2 -Ray results we've downloaded.) In these maps, the RAMSES-RT results look very similar to those of C 2 -Ray. The xHI-maps show stronger ionization immediately around the corner source in the C 2 -Ray result, and correspondingly the temperature and density maps show this corner gas is also hotter and more diffuse in the C 2 -Ray result than in RAMSES-RT. Conversely, the photo-heating region is somewhat further-reaching in the RAMSES-RT result than in C 2 -Ray, as can be seen in the pressure and temperature maps. These small differences are likely due to the different approaches in approximating multi-frequency. Notably, the C 2 -Ray maps stand out in a very similar way when com- pared to most of the corresponding maps from other codes in Il09, i.e. a stronger effect close to the radiative source but shorter-reaching photo-heating.
To paint a more quantitative picture, Fig. 22 compares radial profiles of the same quantities (xHI, P , T , nH and M ) for RAMSES-RT and C 2 -Ray at 10, 200 and 500 Myr. The ionization state profiles (top left) indeed show C 2 -Ray to ionize the gas more strongly close to the radiative source, but RAMSES-RT to ionize more strongly beyond the I-front. The I-front itself is however at very similar positions at all times. The pressure and temperature plots show the same thing, but apart from these minor differences at the extreme ends the shapes are very similar. The density plots show that C 2 -Ray has more has more diffuse gas close to the source as a result of the stronger photoheating, and also it appears to have a more pronounced backflow peak around 200 Myr (this double peak is a temporary effect of photo-heating by high-energy photons beyond the I-front). The smaller backflow peak in RAMSES-RT is perhaps in part a relic of on-the-fly refinement, though most of the codes in Il09 actually have backflow peaks similarly smaller than that of C 2 -Ray. Unfortunately we can't compare the Mach profiles directly, but the RAMSES-RT profiles do look very similar in shape to those presented in Il09 (see their Fig. 15).
Finally, Fig. 23 shows how the position and velocity of the I-front (defined as where the radial average of xHII is equal to 0.5), for RAMSES-RT and C 2 -Ray. The plots for the two codes are virtually identical, the only noticeable difference being a slight initial lag in the front speed. One might attribute this to the reduced speed of light in the RAMSES-RT run, but actually most other codes described in Il09 have a very similar lag in the initial front speed compared to C 2 -Ray.
The fraction of the volume refined to the effective resolution of 128 3 cells is 28% at the end of the run, and the computational time is roughly half that of an analogous uniform grid run. The runs clock in at about double the cpu hours of test 1, even though test 1 had roughly twice the number of timesteps to perform, due to a smaller box width. This gives a qualitative idea of the added cost of adding two more photon groups (test 1 had one group) and coupling with the hydrodynamics, which totals to about four times the computational load. All in all, the RAMSES-RT results for this test compare very well with most of the codes presented in Il09.
The RAMSES-RT result differs slightly from that of C 2 -Ray in some aspects, most notably in the form of weaker photoheating and ionization close to the radiative source and wider I-fronts. However, these are precisely the aspects where C 2 -Ray stands out from the other codes presented in Il09.
6.7 Il09 test 6: HII region expansion in a r −2 density profile This test mimics a radiative source going off in a dense cloud, e.g. a stellar nursery. The setup is much like that of the preceding test 5, the main difference being that the gas is here inhomogeneous, the box is much smaller, L box = 0.8 kpc in width, and the radiative corner source is a hundred times more luminous, i.e. it radiates atṄγ = 5 × 10 50 photons s −1 .
As in the previous test we don't apply the OTSA. The base resolution is 64 3 cells, but on-the-fly refinement on nH and xHII gradients ensures the prescribed effective resolution of 128 3 cells at ionization and shock fronts. The initial temperature is 100 K everywhere and the running time is 75 Myr. The dense cloud is centered on the corner source and is set up with a spherically symmetric, steeply decreasing power-law density profile with a small flat central core of gas number density n0 = 3.2 cm −3 and radius r0 = 91.5 pc: The Strömgren radius for the core density, given by Eq. 47, is rS ≈ 70 pc, which lies within the flat core. Thus, the I-front makes an initial transition from R-type to D-type within the core, and then may accelerate back to R-type as it expands into decreasingly dense gas outside the core.
We first compare the evolution of the position and speed of the I-front, which is plotted in Fig. 24 for RAMSES-RT and the Capreole+C 2 -Ray combination. The I-front moves very quickly (R-type) to ≈ 70 pc within the first fraction of a Myr, stops for while and then starts to expand again with the flow of the gas. Both the speed and position compare well with C 2 -Ray. The initial speed in C 2 -Ray has an apparent lag which is due to under-sampling in the front positions at early times, as noted by Il09. Other code results which are better sampled in Il09 show initial speeds that are virtually identical to the RAMSES-RT plot, especially those of the RH1D code. The final front position in RAMSES-RT is slightly further out than that of C 2 -Ray, though very similar to at least three of the codes in Il09 (Flash-HC, Licorice and RSPH). It also appears that the C 2 -Ray front is starting to accelerate slightly at the end, whereas the RAMSES-RT front is about to approach constant speed; RAMSES-RT also agrees with most other Il09 codes on this point. Fig. 25 shows the overall structure of ionization and the gas at 25 Myr, here with a comparison between RAMSES-RT (upper two rows) and TVD+C 2 -Ray (bottom row). The Capreole+C 2 -Ray version of this test is sensitive to so-called 'carbuncle' numberical instabilities (see Sec. 4.2 in Il09), so we compare here to the more stable and symmetric combination of C 2 -Ray coupled to the Eulerian TVD solver of Trac & Pen (2004) (used only in this test). In addition to the default RAMSES-RT run with on-the fly AMR refinement, we show here in the middle row results from an identical RAMSES-RT run with the base resolution set to 128 3 cells and AMR refinement turned off. There are slight spherical asymmetries appearing in the top row maps, in particular the xHII, T and Mach maps, and the middle row maps are presented here to show that (the first) two of these are purely artifacts of on-the-fly AMR refinement. The slightly square shape of the inner region in the Mach map however does not seem to be due to refinement and is likely rather a grid artifact which is amplified by the radially decreasing density. It should also be noted that the other plots produced for this test (I-front, Fig. 24 and radial profiles, Fig. 26) are absolutely identical regardless of whether on-the-fly refinement is used or the full resolution applied everywhere, suggesting that AMR refinement produces very robust results. The difference in runtime between the AMR and non-AMR runs is actually not much: the AMR run completes in about 2/3 of the ∼ 640 cpu hours taken for the non-AMR run. This lack of speedup is due to a combination of a large portion of the grid being refined (∼ 60% by volume when most), a shallow refinement hierarchy (one level of refinement) and overhead in refinement-related computations.
As usual the I-front is considerably wider in RAMSES-RT than in the C 2 -Ray results, though we don't find the same discrepancy as in the previous test between the photoheating intensity close to the source (also, there is no such discrepancy here between C 2 -Ray and the other codes in Il09). The two maps furthest to the right, of density and Mach number, show the expanding shell of dense gas due to photoheating. Here the shell appears considerably thinner in RAMSES-RT than in TVD+C 2 -Ray, and indeed TVD+C 2 -Ray appears to have the thickest density shell of any of the codes in Il09 (Capreole+C 2 -Ray included, but here there are also severe asymmetries). The RAMSES-RT maps compare well with the C 2 -Ray ones, and to most of the maps in Il09, and don't show any I-front instabilities that seem to have a tendency to come up in this test (and Il09 do show that these are numerical and not physical instabilities). Fig. 26 shows a comparison between RAMSES-RT and TVD+C 2 -Ray for radially averaged profiles at 3, 10 and 25 Myr of the ionization state, pressure, temperature, density and Mach number. The comparison is generally very good. The I-front (and corresponding density shock) lag a little behind in C 2 -Ray, but it actually lags a little behind all but one code in this test in Il09, and RAMSES-RT is spot-on compared with those others in every respect.
All in all, RAMSES-RT thus performs well on this test, and no problems appear that are worth mentioning.
6.8 Il09 test 7: Photo-evaporation of a dense clump The setup of this test is identical to test 3 in Il06, where UV radiation is cast on a gas cloud, creating a shadow behind it and a slowly-moving I-front inside it. Here however, since the hydrodynamics are turned on, photo-heating causes the cloud to expand outwards and simultaneously contract at the center. We recap the setup: The box is L box = 6.6 kpc in width. A spherical cloud of gas with radius r cloud = 0.8 kpc is placed at (xc, yc, zc) = (5, 3.3, 3.3) kpc from the box corner. The density and temperature are n out H = 2 × 10 −4 cm −3 and T out = 8000 K outside the cloud and n cloud H = 200 n out H = 4 × 10 −2 cm −3 and T cloud = 40 K inside it. From the x = 0 boundary a constant ionizing flux of F = 10 6 photons s −1 cm −2 is emitted towards the cloud. The simulation time is 50 Myr, considerably longer than the 15 Myr in the corresponding pure RT test. The base resolution is 64 3 cells, but on-the-fly refinement on nH, xHI and xHII gradients ensures the prescribed effective resolution of 128 3 cells at ionization and shock fronts. In order to best capture the formation of a shadow behind the cloud, we focus on a RAMSES-RT run with the HLL solver, but we also show some results with the usual GLF solver. Fig. 27 shows slices in the xy-plane through the middle of the box of various quantities at 10 and 50 Myr, for the RAMSES-RT result and C 2 -Ray for comparison 8 . As in the corresponding pure RT test, it can be seen from the xHI maps that the shadow behind the cloud is less conserved with RAMSES-RT than with C 2 -Ray, though the HLL solver does a much better job though than GLF. However, the diffusion of photons doesn't have a large impact on the resulting dynamics, or even the propagation of the I-front along the axis of symmetry. The shadow becomes thinner towards the end of the run with all codes in Il09, though it is thinner than most in RAMSES-RT+HLL, and it pretty much disappears in RAMSES-RT+GLF. The shadow thickness in RAMSES-RT+HLL is still comparable at 50 Myrs to the results of RSPH, Zeus-MP and Licorice in Il09. The pressure maps of RAMSES-RT+HLL, C 2 -Ray and other codes in Il09 are very similar both at 10 and 50 Myrs, though C 2 -Ray, and also to some extent Flash-HC and Licorice have a fork-like shape inside what remains of the shadow at 50 Myr. The other codes have the same shape as RAMSES-RT+HLL in this region. The temperature maps are similar as well, though the backwards-expanding cloud shell seems to be slightly less shock-heated in RAMSES-RT than most other codes. The shell expands in a very similar way for the two codes, as can be seen in the density and Mach slices. The expansion goes a bit further, though, in RAMSES-RT. Also, the expanding cloud seems to develop a slightly hexagonal shape in RAMSES-RT, an effect which is not apparent in any of the codes in this test in Il09 (though there is a hint of it in the Flash-HC result). It can only be speculated that this is a grid artifact. To be sure it doesn't have to do with the on-the-fly refinement we ran an identical experiment with a base resolution of 128 3 cells and no refinement in RAMSES-RT+HLL. The RAMSES-RT+HLL maps and plots presented here are virtually identical to this non-refinement run, except of course for graininess in the slice maps. None of these discussed effects (hexagons and a slightly over-extended I-front compared to other codes) are thus due to on-the-fly refinement. As in the previous test, the speed-wise gain in using AMR is not a lot: the AMR run completes in about half of the ∼ 64 cpu hours taken for the non-AMR run. Again the relatively modest speedup is due to a combination of a large portion of the grid being refined (∼ 30% by volume when most), a shallow refinement hierarchy, and refinement-related overhead. With deeper refinement hierarchies in cosmological simulations, the speedup can be much greater, but a quantitative demonstration is beyond the scope of this paper.
Next we turn our attention to the evolution of the position and speed of the I-front along the x-axis of symmetry through the box. This is presented for the RAMSES-RT (HLL and GLF) and C 2 -Ray runs in Fig. 28a. The I-front propagation is considerably different between RAMSES-RT and C 2 -Ray, but actually C 2 -Ray considerably stands out here from other codes in Il09. For the first 7 Myrs or so, the RAMSES-RT front lags behind that of C 2 -Ray and in fact all the codes in Il09. This is due to the reduced speed of light: before hitting the cloud, the photons have to travel from the left edge of the box through a very diffuse medium -so diffuse that here the I-front speed apparently is approaching the speed of light, or is at least considerably faster than the one-one-hundredth of the light speed which is used in the RAMSES-RT run. However, once the I-front in the RAMSES-RT run has caught up, the reduced light speed should have a negligible effect on the results. After roughly 7 Myr, the RAMSES-RT I-front overtakes C 2 -Ray front, and stays ahead of it for the remainder of the run. This however is also the case for most of the codes in Il09; their I-front is ahead of the C 2 -Ray front, and four out of six codes end up with the I-front at ∼ 5.6 kpc. The RAMSES-RT+HLL front ends up at ∼ 5.7 kpc, so slightly ahead of what is typically found in Il09. Using the GLF solver instead of HLL has the effect that the I-front disappears soon after 40 Myrs, which is due to diffusive photons eating into the shadow from it's edges, but up to that point the I-front evolution is much the same. RAMSES-RT also reproduces the retreat of the I-front between roughly 30 and 40 Myrs, which is seen in all runs in Il09. This momentary negative speed is due to the expansion of the cloud and the D-type movement of the I-front with the gas. Fig. 28b shows histograms of the gas temperature and Mach number at 10 and 50 Myr in the RAMSES-RT+HLL and C 2 -Ray runs. The shapes of the histograms are very similar between the two codes (and are also very similar to RAMSES-RT+GLF, which is not shown).
Finally, Fig. 29 shows a comparison between RAMSES-RT and C 2 -Ray of profiles along the x-axis of symmetry of the various quantities at 1, 10 and 50 Myrs. The profiles compare badly at 1 Myr, but as already discussed this is simply due to the I-front having not caught up at this early time when using the reduced speed of light. At later times the profiles generally compare well, though we see these effects which have already been discussed, of a further expanding densityfront out of the original cloud, and a further progressed Ifront. The RAMSES-RT profile plots show a staircase effect which is most obvious in the 50 Myrs plot at the radial interval 0.45 r/L box 0.75: this is simply due to the grid being unrefined at this x-interval along the axis of symmetry, i.e. at the effective base resolution of 64 3 cells per box width. The run with the full resolution and no AMR refinement shows no staircases, but otherwise the results are identical to those shown here.
We have made an alternative run with RAMSES-RT+HLL with the speed of light fraction set to fc = 1/10 rather than the default 1/100, and here the initial evolution of the I-front position and radial profiles at 1 Myr are almost identical to those of C 2 -Ray. At later times the results are very much in line with those where fc = 1/100, except the I-front position is slightly more advanced at 50 Myr, or at 5.78 kpc rather than at 5.71 kpc.
In summary, RAMSES-RT performs well on this test with no apparent problems. The reduced light speed (fc = 1/100) has very little effect on the results and on-the-fly refinement gives results which are identical to the fully refined simulation with a homogeneous 128 3 cells grid. Even using the diffusive GLF solver retains much of the results (I-front development, cloud expansion), except that the I-front disappears a bit prematurely.

RT test conclusions
RAMSES-RT performs very well on all the tests from Il06 and Il09, with no discrepancies to speak of from expected results or those from other codes.
The most notable discrepancies clearly result from the reduced speed of light approximation, which leads to I-fronts that are initially too slow compared to full speed of light runs -or infinite speed, as is the case for many of the codes compared against from the RT comparison project. In test 4, the high-z cosmological field, we actually demostrated the reverse, where the codes we compared to had considerably premature I-fronts as a consequence of their infinite light speed approximations. Our shadows are considerably shorter lived with the GLF intercell flux function than those of the other codes (most of which use ray-tracing schemes). This can be fixed for problems involving shadows and idealized geometries by using the HLL flux function instead, but as we showed in §3.2 the sacrifice is that isotropic sources become anisotropic. Many codes in the RT comparison project show various instabilities and asymmetries in ionization fronts; no such features are manifested in the RAMSES-RT results.

DISCUSSION
In this paper, we have presented a new implementation of radiation hydrodynamics in the RAMSES code. It is based on a moment representation of the radiation field, where we have used the M1 closure relation to define a purely local variable Eddington tensor. Because the resulting system is a set of hyperbolic conservation laws, we have exploited the Godunov methodology to design a time-explicit, strictly photon conserving radiation transport scheme. The resulting algorithm is first order accurate in space and time, and uses various Riemann solvers (GLF and HLL) to compute radiation fluxes. The main novelty compared to our previous implementation (AT08) is the coupling between gas and radiation, resulting in a fully consistent radiation hydrodynamics solver, and the introduction of adaptive mesh techniques in the radiation transport step, making use of both the AMR and parallel computing capabilities of RAMSES. Overall, the code was quite easy to implement, owing to the explicit nature of the time integration scheme. The price to pay is the need to resolve the propagation of hyperbolic waves traveling at or close to the speed of light. Among many different options available to overcome this constraint, we have chosen to use the "reduced speed of light" approximation. This approximation is valid when the propagation speed of I-fronts is still slower than the (reduced) light speed. We have developed a recipe to assess the validity of this approximation, based on the light-crossing time of Strömgren spheres. We have verified that this framework indeed allows us to estimate in advance the speed-of-light reduction factor reliably. We have shown, for example, that in cosmological problems, such as cosmic reionization, using the correct value for the speed of light is crucial, and using either a reduced or an infinite speed of light (like in some ray-tracing codes) might result in large inaccuracies.
This new algorithm has already been used in galaxy formation studies, exploiting the coupling between radiation and hydrodynamics offered by RAMSES-RT. In Rosdahl & Blaizot (2012), we studied the impact of ionizing radiation in determining the thermal state of cold filaments streaming into high redshift galaxies, allowing us to make accurate observational predictions and demonstrating a possible link between cold streams and Lyman-alpha blobs. More recently, we have also explored the role of ionizing radiation in the overall efficiency of stellar feedback (Geen et al. 2013, Powell et al. 2013. Beyond ionizing radiation, possible extensions of RAMSES-RT are the inclusion of photo-dissociating radiation and the thermochemistry of molecules, as well as the effect of dust as an additional source of opacity and thermal regulation inside star forming galaxies. This would require introducing additional photon groups (such as a far UV and IR photons) and the associated microphysics, but the overall methodology would remain very similar.
In order to improve the current algorithm, we have many possibilities ahead of us. One obvious development is to develop a second-order sequel of our current first order Godunov solver. Second-order Godunov schemes, both in time and space, are used routinely in hydrodynamics codes (such as the MUSCL scheme in RAMSES). This might reduce significantly the rather large diffusivity of our current implementation. However, since photo-ionization and photodissociation problems are governed to a large extent by the thermo-chemistry, it is not clear how much the accuracy of the results would depend on the advection scheme. A second route we would like to explore in the future is the optional introduction of radiation sub-cycles during each adaptive hydro step. This is quite challenging since it would in principle require decoupling in time of the various AMR levels, resulting in the loss of strict photon conservation. In some cases, however, it is advantageous to sacrifice the exact number conservation of photons in favour of modelling the correct speed of light with many radiation sub-cycles. In any case, this would offer us a new tool with greater flexibility. Along the same lines, because of the fundamentally different propagation properties of I-fronts in the IGM on one hand and deep inside galaxies on the other, we could couple RAMSES-RT to ATON: use ATON to transport radiation on the coarse grid with GPUs at the full speed of light, and use RAMSES-RT on the fine AMR levels at a reduced light speed. This would require us to define two photon group populations that mirror each other: a large scale, low density photon population that propagates at the correct speed of light and makes use of GPU acceleration (if available), and a small scale, high density photon population that makes use of the "reduced speed of light" approximation. Coupling properly the two photon group populations will of course be quite challeng-ing and at the heart of this new avenue of research. A last development we have in mind is the introduction of radiation pressure as a new channel of coupling radiation with hydrodynamics. This is highly relevant for studies focusing on radiation pressure on dust, from both young star clusters and supermassive black holes. density, momentum density, energy density, hydrogen and helium ion abundances, photon densities and fluxes 9 ) and evolves numerically over a time-step ∆t the thermochemical state UT = ε, xHII, xHeII, xHeIII, N1, .. NM , F1, .. FM (A1) (where ε = E − 1/2ρu 2 is the thermal energy density), i.e. solves the set of 4 + 2M coupled equations where S ≡UT . Due to the stiffness of the thermochemistry equations, it is feasible to solve them implicitly, i.e. using S(U t+∆t T ) on the right hand side (RHS), which guarantees stability and convergence of the solver. However a fully implicit solver is complicated in implementation, computationally expensive and not easily adaptable to changes, e.g. a varying number of photon groups or additional ion/chemical abundances. Instead we take an approach inspired by Anninos et al. (1997). The idea is to solve one equation at a time in a specific order, and on the RHS use forward-in-time (FW) values, i.e. evaluated at t + ∆t, wherever available, but otherwise backwardsin-time (BW) values, evaluated at t. So for the first variable we choose to advance in time, there are no FW variables available. For the next one, we can use the FW state of the first variable, and so on. In that sense the method can be thought of as being partially implicit.
The cell-thermochemistry is called once every RT timestep of length ∆tRT , but in each cell it is split into local sub-steps of length ∆t that adhere to the 10% rule, where ∆UT is the change in UT during the sub-step. The RT step thus contains a loop for each cell, which calls the thermo step(UT, ∆t) routine once or more often: first with ∆t = ∆tRT , then possibly again a number of times to fill in ∆tRT if the first guess at ∆t proves too long to meet the condition set by (A3).
The thermo step(UT, ∆t) routine performs the following tasks: (i) N and F update (ii) E update (iii) xHI update (iv) xHeII and xHeIII update (v) Check if we are safe to use a bigger time-step Tasks (ii) to (iv) are in the same order as in Anninos et al. (1997), but they don't include radiative transfer in their code, so there is no photon update. The argument we have for putting it first rather than anywhere else is that the photon densities appear to be the most dynamic variables and so are also most likely to break the time-step condition (A3). This we want to catch early on in the thermochemistry step so we avoid doing calculations of tasks (ii) to (v) that turn out to be useless because of the too-long time-step.
We now describe the individual tasks. Temperature dependent interaction rates frequently appear in the taskstheir expressions are given in Appendix E. The temperature can at any point be extracted from the energy density and ionization state of the gas via where γ is the ratio of specific heats (usually given the value of 5/3 in RAMSES, corresponding to monatomic gas), mH the proton mass, kB the Boltzmann constant and µ is the average mass per particle in the gas, in units of mH .

(i) Photon density and flux update
The photon number densities and fluxes, Ni and Fi, are updated one photon group i at a time. For the photon density the equations to solve are whereṄi represents the time derivative of Ni given by the RT transport solver (which is nonzero only if the smoothed RT option is used), Ci represents photon-creating re-combinations, and Di represents photon-destroying absorptions. The creation term is non-existent if the OTSA is used (emitted photons are assumed to be immediately reabsorbed), but is otherwise given by

Ci =
Hii,Heii,Heiii where the b rec ji factor is a boolean (1/0) that states which photon group j-species recombinations emit into and α A j and α B j are the temperature dependent case A and B recombination rates for the recombining species. The photon destruction factor is given by where cr is the (reduced) light speed and σ N ij is the crosssection between species j and photons in group i.
Photon emission from recombination is assumed to be spherically symmetric, i.e. to go in all directions. It is therefore purely a diffusive term, and the photon flux equation only includes the photo-absorbtions: whereḞi is the time derivative used only in smoothed RT and the destruction factor remains as in (A7). Equations (A5) and (A8) are solved numerically using a partly semi-implicit Euler (SIE) formulation, in the sense that they are semi-implicit in the photon density and flux but otherwise explicit (in temperature and the ion abundances). A tiny bit of algebra gives: where all the variables at the RHS are evaluated at the beginning of the time-step, i.e. at t. For each photon group update, the 10% rule is checked: if the cool step routine returns with an un-updated state but instead a recommendation for a new time-step length ∆tnew = 0.5 ∆t, so the routine can be called again with a better chance of completing.
(ii) Thermal update Due to the dependency of µ on the ionization fractions it is easiest to evolve the quantity where µ can be extracted via with X and Y = 1 − X the hydrogen and helium mass fractions, respectively. Here we ignore the metal contribution to µ, which in most astrophysical contexts is negligible. The temperature is updated by solving where Λ ≡ε = H + L, H is the photoheating rate and L the cooling rate. These rates are calculated as follows: The photoheating rate H is a sum of the heating contributions from all photoionization events: where ν is photon frequency, F (ν) local photon flux and j photoionization energies. With the discretization into M photon groups, (A15) becomes H = Hi,Hei,Heii where¯ i, σ N ij , σ E ij and are the photon average energies, average cross sections and energy weighted cross sections, respectively, for ionization events between group i and species j (see Eqs. 9-11).
The primordial cooling rate L is given by Heiii (T ) ne nHeIII + θ(T ) ne (nHII + nHeII + 4nHeIII) where the various cooling processes are collisional ionizations ζ, collisional excitations ψ, recombinations η, dielectronic recombinations ω, bremsstrahlung θ and Compton cooling , all analytic (fitted) functions of temperature taken from various sources. The complete expressions are listed (with references) in Appendix E. If the OTSA is used, the η A coefficients are replaced with η B .
The temperature update (A14) is solved numerically using semi-implicit formulation in Tµ, using FW values of photon densities and BW values of H and He species abundances. The temperature is updated to where K ≡ (γ−1)m H ρk B . The temperature-derivative, Λ ≡ ∂L ∂Tµ , is found by algebraically differentiating each of the primordial cooling rate expressions in the case of L (and using ∂L ∂Tµ = µ ∂L ∂T ). The temperature derivative of the heating rate is zero.
With T t+∆t µ in hand, the time-stepping condition is cool step is re-started with half the time-step length. In tests we've found that the usual time-step constraint given here is not enough to ensure stability, as the temperature in some cases oscillates, even in a divergent way. Λ and Λ are both evaluated backwards in time, i.e. at t, and the large difference that can exist in these values from t to t + ∆t appears to cause these instabilities. To fix that we include also a first-order time-step constraint on the temperature, i.e. if the time-step length is halved. With this fix, we have not seen further temperature oscillations, but there is no guarantee that numerical instabilities are eliminated.
where βHi(T ) is the rate of collisional ionizations by electrons and α A HII (T ) the case A hydrogen recombination rate, which is replaced here by α B HII if the OTSA is used. In terms of ionization fraction, (A21) becomes where we have in the second line separated the rates into Hii creation C and destruction D, and in the third line collected multiples of xHII.
To prevent stiffness-induced instabilities, we have gone for an approach which is semi-implicit in xHII: the 10% rule to be obeyed, successive calls to cool step will quickly fix that by halving ∆tTC until the rule is no longer broken, and only then will cool step start to return updated values of UT .

APPENDIX B: STELLAR UV EMISSION AND DERIVED PHOTON ATTRIBUTES
In the photon injection step in RAMSES-RT ( §3.1), the task is to inject photons into each grid cell corresponding to the luminosities of stellar particles that reside in it. Here we describe how we derive these luminosities from stellar energy distribution (SED) models, along with the photoionization cross sections and energies for each photon group.

B1 Stellar luminosities
Stellar particles in RAMSES represent stellar populations, so it makes sense to use SED models to infer their luminosities. RAMSES-RT can read SED tables at startup and derive from them stellar luminosities for photon injection, as well as photon group attributes that can be updated to reflect the average emission from the stellar particles populating the simulation. We have hitherto used the SED model of Bruzual & Charlot (2003) (BC03), but it can be replaced with any other model, e.g. Starburst99 (Leitherer et al. 1999), as long as the file format is adjusted to match. The model should come in the form of spectra, J λ (τ, Z), giving emitted energy in Solar units per Solar mass per wavelength per second, binned by stellar population age and metallicity. In Fig. B1 we show J λ from BC03 and Starburst99 at solar metallicity for various population ages.
Age and metallicity dependent population luminosity L, given in number of photons emitted per second into photon group i, is calculated from the SED model by where Jν = c/ν 2 J λ(ν) . The cumulative population luminosity is then Since both the photon injection and the calculation of photon group attributes are done on the fly, Li(τ, Z) and Πi(τ, Z) must be evaluated as quickly as possible for given stellar particle ages and metallicities. Values of Li and Πi are therefore only calculated from the SED spectra via (B1) and (B2) at simulation startup, and tabulated with equallyspaced logarithmic bins of age and metallicity, so that they can be retrieved with minimum computational effort via linear interpolation, e.g. when injecting photons into cells via Eq. 17.

B2 Photon group attributes
There are three sets of global attributes for each photon group. These are average photon energies¯ i, average photoionization cross sections σ N ij and energy weighted cross sections σ E ij , that are defined in §2 (Eqs. 9-11). For an age and metallicity dependent reference spectrum Jν (τ, Z), these arē Jν dν Jν dν .
Since there are three ionizeable species in the current implementation of RAMSES-RT, each photon group has three values of σ N and three of σ E . These attributes can be set as run parameters to reflect some typical stellar spectra, e.g. a blackbody or a SED. It can also be left to RAMSES-RT to set them on the fly to reflect the in-simulation stellar populations, using the expressions (B3)-(B5), with the loaded SED spectra representing Jν and the expressions from Verner et al. (1996) for σνj (see Appendix E4). Due to the averaged nature of the photon groups, we must however suffice to set the group attributes to reflect the average stellar emission in the simulation, weighted by the stellar luminosities 10 . If this option is used, the photon group attributes are updated every n coarse time-steps (where n is an adjustable parameter) by polling all the stellar particles in the simulation and setting for each group i and species j, The values of each stellar particle's Li(τ , Z ),¯ ij (τ , Z ), σ N ij (τ , Z ) and σ E ij (τ , Z ) are interpolated from tables that are generated at startup via (B1) and (B3)-(B5).
Although one is free to use many photon groups to resolve frequencies, it is practical to only use a handful, due to limitations in memory and computation. We typically use three photon groups in our simulations, representing Hi,

0
Hei and Heii ionizing photons, as indicated by vertical lines in the plots of Fig. B1. The stellar luminosities, instantaneous and accumulated, average cross sections and energies for those groups are plotted in Fig. B2 for BC03 (top) and Starburst99 (bottom), as calculated via (B1)-(B5). From the luminosity plots (top rows), it can be seen that the stellar populations emit predominantly for the first ∼ 3 − 6 Myrs and the luminosity drastically goes down as the most massive stars in the population begin to expire.

APPENDIX C: NON-EQUILIBRIUM THERMOCHEMISTRY TESTS
To validate the non-equilibrium thermochemistry in RAMSES-RT we ran one-cell thermochemistry tests, that start at some initial state (temperature, ionization state, photon flux) and evolve over roughly a Hubble time. We are interested here in verifying that our implementation is correct and error free and also in comparing equilibrium vs. nonequilibrium cooling -e.g. Cen & Fang (2006) report that the methods can produce significantly different results. We compare against the equilibrium thermochemistry of RAMSES which has been modified to use the exact same heating, cooling and interaction rates as RAMSES-RT.
We test to see (i) whether the thermochemistry of RAMSES-RT is stable, i.e. if the stiffness of the equations results in any sudden divergence or 'wiggles' in the evolution of the gas, (ii) whether RAMSES-RT evolves the ionization fractions towards the correct states predicted by the equilibrium solver of RAMSES, and (iii) whether the RAMSES and RAMSES-RT evolve the temperature towards the same final value.
There are four tests: first we disable cooling and evolve only the ionization states of hydrogen and helium at different constant temperatures in a zero UV radiation field, and see if we reach equilibrium ionization states (predicted by RAMSES). Then we turn on a constant UV radiation field and again see if we reach equilibrium states. Then we turn on cooling, and for two sets (zero, nonzero radiation field) see if the temperature evolution is comparable to RAMSES equilibrium cooling from the same initial conditions.

C1 Ionization convergence at constant temperature and zero ionizing photon flux
In the first test, cooling is turned off and we check for a range of densities, temperatures and initial ionization states whether we get a convergence of the ionized fractions towards their equilibrium states, as predicted by RAMSES, assuming zero flux of ionizing photons. Fig. C1 shows the results. Each panel of 3 × 6 plots in the figure represents an evolution given the constant temperature written to the right of the panel, and shows how the ionized fractions, xHII, xHeII and xHeIII, evolve from different (color-coded) starting states xi = xHII = xHeIII (the HeII fraction always starts at zero). A black dashed line in each plot shows the equilibrium ionization fraction for the given temperature and species (which is gas density independent in the case of zero ionizing flux). Each column of plots represents a (non-evolving) hydrogen number density.
The non-equilibrium ionization fractions always evolve towards the equilibrium ones, at a rate which depends on gas density, as expected. It can even take longer than the age of the Universe to reach equilibrium for the most diffuse gas (nH 10 −6 cm −3 ), which indeed is a significant difference from the equilibrium assumption. If we zoom in around the equilibrium states we find a difference between the calculated equilibrium state and the evolved one which is typically around one in ten-thousand -this simply corresponds to the allowed error in the iterative equilibrium calculation, and can be decreased at will by reducing this error margin.  Figure B2. Heii, Hei, and Hi ionizing luminosities and photon group attributes derived from the BC03 (top panel) and Starburst99 SED models (bottom panel), as functions of age (x-axis) and metallicity (colors). The plot columns represent the three photon groups. Top rows show stellar luminosity, in the number of photons that goes into each group per second per solar mass. Second rows show accumulated number of photons emitted. Third rows show the average photon energies per interaction. Bottom rows show average cross sections per interaction.
C2 Ionization convergence at constant temperature and nonzero ionizing photon flux This is the same as the previous test, except now we apply a constant flux of 10 5 ionizing photons s −1 cm −2 through the cell, assuming the spectrum of a blackbody at 10 5 K. Fig. C2 shows the results. The black dashed lines in each plot show the equilibrium state which now is density dependent -the denser the gas the harder it is for the radiation field to battle against recombinations. Again the nonequilibrium ionized state always evolves towards the equilbrium one at a gas density dependent rate, though note that  here it takes a maximum of ∼ 10 Myr, which is much shorter than it can take in the zero photon flux case.
C3 Temperature convergence with zero ionizing photon flux Now cooling is turned on, and we compare the RAMSES-RT non-equilibrium temperature evolution with that of equilibrium RAMSES (though keep in mind it has been adjusted to contain the exact same cooling rates as RAMSES-RT). Each of the 5 rows in fig. C3 shows cooling for a range of decreasing initial temperatures, from top to bottom. The color-codings (initial ionization states) and columns (hydrogen number densities) are the same as before. The solid colored lines show non-equilibrium cooling in RAMSES-RT and the black dashed lines represent equilibrium cooling in RAMSES starting from the same temperature. Clearly the temperature evolution is quite similar between equilibrium/non-equilibrium cooling, especially if the initial ionization fraction is 'correct', i.e. if it matches the equilibrium one at the initial temperature.
The final temperature reached in the non-equilibrium case is usually a bit lower than in the equilibrium case. This is independent of gas density and initial temperature (as long as the initial temperature allows for cooling to happen). The reason for this is that the non-equilibrium ionization evolution lags behind the instantaneous equilibrium one, so there is always a somewhat larger reservoir of electrons in the non-equilibrium case. Electrons are the primary cooling agents, and complete electron depletion completely stops cooling, so it makes sense that if the electrons deplete more slowly, cooling is more effective and can bring the gas to a lower final temperature.

C4 Temperature convergence with nonzero ionizing photon flux
This is the same as the previous test, except now we apply a constant flux of 10 5 ionizing photons s −1 cm −2 , assuming the spectrum of a blackbody at 10 5 K. The results are shown in Fig. C4. Things are much the same as before, except that the non-equilibrium temperature seems to converge to a value which is much closer to the the equilibrium one -because of the ionizing flux there is always a reservoir of electrons both in the equilibrium and non-equilibrium evolution, which makes for a much closer match in the final temperature.
Although the final temperature reached is identical between the two methods, the evolution towards that final temperature can be quite different, depending on the initial ionization states.
A zoom-in on one of the plots is shown in Fig. C5, and reveals that there is very little difference between the final temperatures reached. The little difference there is results from interpolation from cooling-rate tables in RAMSES equilibrium cooling and it can be decreased further by increasing the size of these tables. x i =0.0 x i =0.2 x i =0.5 x i =0.8 x i =1.0 Figure C5. Close-up of temperature convergence, for the UV inclusive test with initial temperature T ≈ 10 5 K and n H = 10 −2 cm −3

C5 Thermochemistry tests conclusions
The main conclusions of the one-cell thermochemistry tests are: • We always eventually reach the equilibrium ionization state with the non-equilibrium method...
• ...but this can take a very long time to happen for diffuse gas, even more than a Hubble time.
• Non-equilibrium temperature evolution of the gas is quite dependent on the initial ionization fraction of the gas at intermediate temperatures and low densities...
• ...but in the end we reach the same or at least a very similar temperature as in the equilibrium case.
• The convergence of the non-equilibrium solver towards the results of the equilibrium solver of RAMSES, given the same cooling rate expressions, suggests that our thermochemistry solver is robust and correct.

APPENDIX D: ON MULTI-STEPPING IN THE AMR LEVEL HIERARCHY
As discussed in §5.2, solving hydrodynamics over an AMR grid with a multi-stepping approach always leaves ill-defined states at inter-level boundaries between the start and finish of the coarse level timestep. This imposes severe constraints on how the RT can be coupled to the hydrodynamics, and essentially means that RT cannot be sub-cycled within multistepping hydrodynamics. Here we will clarify this point in detail.
Hydrodynamic advection across the boundaries of a cell is performed in an operator-split fashion, such that the advection is solved separately across a discretized timestep for each boundary. In order for the solver to be consistent, i.e. for the result at the end of the timestep to be independent of the order in which the boundaries are accounted for, the solver must work from the same initial cell state U for all the inter-cell updates. Thus, a copy is first made of the original cell states involved, i.e.

U → U,
where we can term U the source state and U the destination state. Using U as source terms for the intercell fluxes, the advection can be solved with some computational method (e.g.
lmax-1 lmax 3 5 4 Figure D1. Level gas state updates via intercell fluxes also perform partial gas updates in neighbouring cells at level − 1. The example shown corresponds to the hierarchy from Fig. 5.
Steps 3 and 4 at the finest level also include partial updates of neighbouring max −1 cells, but these neighbour cell states are not fully updated until all the intercell fluxes are taken into account, which is in step 5 from Fig. 5.
Godunov solver for the hydrodynamics in RAMSES and an HLL/GLF flux function for the RT advection in RAMSES-RT), which performs the update on U. To take a concrete example, each RT advection update, Eq. 22, uses U for the update (the LHS term and the first RHS term), but the intercell fluxes are derived from U, i.e. F(U). Once all the updates (6 per cell) have been collected, the cell update is made final by: In the amr step( ) hierarchy in RAMSES, such copies are made of all cells before the AMR recursion, and the update is made final after the recursion has returned and the hydro solver has been called at the current level, i.e. advection has been performed over the timestep over the current level and all finer levels.
This allows cell states to be updated not only at the current level, but also (twice) in all neighbouring cells at the next coarser level. The coarser level update is only partial though, because it only reflects the intercell fluxes across inter-level boundaries, and fluxes across other boundaries (same level or next coarser level) will only be accounted for when the coarser level time-step is advanced. Until then, these coarser level neighbour cells have two gas states, U and U. This is shown schematically in Fig. D1.
If RT subcycling is to be done at each AMR fine-level step, over the whole grid, the question is, which cell state do we use for the thermochemistry, i.e. the interaction between photons and gas, in those inter-level boundary cells?
Choosing one but not the other leads to an obvious and severe inconsistency between the source and destination states. If the thermochemistry does the update on U, then a gas element which is transported from one cell to a neighbour during the following hydro transport is not thermochemically evolved over the time-step, because it originates from U. If instead the update is done on U, a gas element which stays still in any cell over the following hydro transport step is not thermochemically evolved over the time-step. One might then just update both states via thermochemistry, i.e. apply it on each cell twice. This does not really make sense for these inter-level intercell boundary cells that have U = U, as U doesn't represent a true state but is rather an intermediate and temporary quantity that exists between well-defined times. Also, it would be really nontrivial to implement: applying thermochemistry on each of the states also implies transporting the photons through two different states in each cell, which creates alternative timelines for the radiative transfer! Thus, subcycling RT within multi-stepping hydrodynamics in a conservative way is not possible (or at least nontrivial), which has led us to disallow RT subcycling within the hydro timestep in our implementation.

APPENDIX E: INTERACTION RATE COEFFICIENTS ADOPTED IN RAMSESRT
Here we collect the rate coefficients used in RAMSES-RT for hydrogen and helium interactions, which are fitted functions taken from various sources. These are, in order of appearance, collisional ionization rates, recombination rates, cooling rates (collisional ionization, recombination, collisional excitation, bremsstrahlung, Compton and dielectric recombination), and photoionization cross sections.

E1 Collisional ionization rate coefficients
Those are in units of [cm 3 s −1 ] and are taken from Cen (1992), with temperature everywhere assumed in Kelvin:

E2 Recombination rate coefficients
These are all taken from Hui & Gnedin (1997). For readability, we use the following unitless functions: