Impact of the reduced speed of light approximation on the post-overlap neutral hydrogen fraction in numerical simulations of the epoch of reionization

Context. The reduced speed of light approximation is used in a variety of simulations of the epoch of reionization and galaxy formation. Its popularity stems from its ability to drastically reduce the computing cost of a simulation by allowing the use of larger and therefore fewer timesteps to reach a solution. This approximation is physically motivated by the fact that ionization fronts rarely propagate faster than some fraction of the speed of light. However, no global proof of the physical validity of this approach is available and possible artefacts resulting from this approximation therefore need to be identified and characterized to allow its proper use. Aims. In this paper we investigate the impact of the reduced speed of light approximation on the predicted properties of the intergalactic medium. Methods. To this end we used fully coupled radiation-hydrodynamics RAMSES-CUDATON simulations of the epoch of reionization. Results. We find that reducing the speed of light by a factor 5 (20, 100) leads to overestimating the post-reionization average volume-weighted neutral hydrogen fraction by a similar factor ∼5 (20, 100) with respect to full speed of light simulations. We show that the error is driven by the hydrogen photon chemistry by considering the analytical solution for a strongly ionized hydrogen gas in photoionization equilibrium. In this regime, reducing the speed of light has the same effect as artificially reducing the photon density or the hydrogen photoionization cross section and leads to an underestimated ionizing intensity. We confirm this interpretation by running additional simulations using a reduced speed of light in the photon propagation module, but this time we keep the true speed of light in the chemistry module. With this set-up, the post-reionization neutral hydrogen fractions converge to the full speed of light value, which validates our explanation. Increasing spatial resolution beyond a cell size of 1 kpc physical, so as to better resolve Lyman-limit systems, does not significantly affect our conclusions.


Introduction
The epoch of reionization (hereafter EoR), is the period during which the ultra-violet light of the first stars and galaxies ionized the intergalactic medium (IGM), in expanding HII regions of increasing size, finally merging together to produce the almost uniformly ionized IGM we see today. This process ended about 1 billion years after the Big Bang, and is tightly coupled to the formation and growth of galaxies during this period. This period is the focus of a broad effort throughout the astrophysical community, aiming to understand in detail the timeline and identify the main sources contributing to the reionization of the Universe. In order to interpret the current data, and to prepare for the wave of new data expected from current and future redshifted 21cm experiments (e.g. LO-FAR, MWA, HERA, SKA), and upcoming space observatories such as JWST, accurate models of galaxy formation during the epoch of reionization are required. Such models also have implications for near-field cosmology, as it has been shown that the reionization history of the local group may affect the properties of its population of lowmass galaxies (Ocvirk & Aubert 2011;Ocvirk et al. 2013Ocvirk et al. , 2014Gillet et al. 2015).
Modelling the epoch of reionization and its sources is a notoriously difficult problem. One one hand, large boxes in the 100 Mpc or more are required in order to properly sample the variety of environments and rare large voids or early massive galaxies or galaxy clusters, such as in the cosmic dawn simulation suites Ocvirk et al. (2016); Aubert et al. (2016). On the other hand, very high spatial resolution, down to the parsec scale, is necessary to describe the interstellar medium (hereafter ISM) of galaxies, and the complex processes it drives and reacts to, such as star formation, supernova feedback, and the ionizing escape fraction (Geen et al. 2016;Trebitsch et al. 2017;Gavagnin et al. 2017;Butler et al. 2017;Costa et al. 2018). The combination of large volume and high spatial resolution makes numerical simulations of galaxy formation during the EoR extremely challenging in principle, and computationally costly with the rich set of physics required.
In particular, ionizing radiation is a key process in these simulations, both for its impact in reionizing the Universe, and for the additional feedback it provides, and is especially costly to model. Two main classes of solutions have been adopted for treating hydrogen-ionizing radiation. Ray-casting algorithms involve sampling the radiation field produced by each source by a Monte-Carlo technique and follow each photon packet through the computational domain and its interaction with the hydrogen gas (Semelin et al. 2007;Baek et al. 2009;Pawlik & Schaye 2008). Another popular solution consists in considering the photons as a fluid and using the M1 closure relation (Aubert & Teyssier 2008, 2010. This approach has led to the development of fully coupled radiation-hydrodynamics (hereafter RHD) galaxy formation codes, such as RAMSES-CUDATON , EMMA (Aubert et al. 2015), and RAMSES-RT (Rosdahl et al. 2013). In these codes, the radiative timestep is set by the Courant condition, and is typically 100-1000 times shorter than the hydrodynamical timestep of the simulation. This simply reflects the fact that light propagation is much faster than any other process in this framework. This alone makes full RHD simulations very expensive. While some authors have alleviated this problem by accelerating the radiative transfer (hereafter RT) module using GPUs 1 (Aubert & Teyssier 2010;Aubert et al. 2015;Ocvirk et al. 2016), others have resorted to the so-called "reduced speed of light" approximation (hereafter RSL), in order to bring the radiative timestep closer to the hydrodynamic timestep and reduce the overhead due to RT (Kimm & Cen 2013;Rosdahl et al. 2013;Gnedin 2014;Gnedin & Kaurov 2014;Katz et al. 2017;Aubert et al. 2016;Katz et al. 2018). This can drastically reduce the RT computing time. As a consequence, RSL is extremely popular and has allowed groups to perform simulations which would be impossible using the true speed of light. Using a reduced speed of light may appear justified in dense regions (in particular in the ISM), where the ionization front (I-front) is much slower than the real speed of light. However, beyond this, no clear proof of the global physical validity of this approach has ever been provided. For instance, in optically thinner regions of average cosmic density or voids, ionization fronts may be faster than the reduced speed of light used, and the timing and geometry of reionization could be affected. This velocity mismatch has been investigated by Bauer et al. (2015); Gnedin (2016) and Deparis et al. (2018, in prep.), and it appears that using c > 10% of its true value is enough to describe the fronts' propagation without significant lag. However, the validity of RSL with respect to predicting the neutral hydrogen fraction, has not been much studied before in RSL cosmological simulations, and this is what we aim to investigate in the present paper. In Sec. 3, we present simulations of the epoch of reionization performed with RAMSES-CUDATON using 100%, 20%, 5% and 1% of the speed of light, and show that RSL simulations systematically overpredict the postreionization neutral hydrogen fraction. We then interpret this result using analytical arguments and further simulations. We finally close the paper with our conclusions.

Simulations
We used RAMSES-CUDATON ), a coupling between RAMSES (Teyssier 2002) and ATON/CUDATON (Aubert & Teyssier 2008, 2010. Thanks to the CUDA optimization of ATON, RAMSES-CUDATON is able to take advantage of GPU acceleration and routinely performs RHD simulations with the full speed of light. It is therefore an adequate testbed for investigating artefacts in RSL simulations, which are less computationally expensive. The setup used here is 1 Graphics Processing Units a 4 h −1 Mpc box discretized on a fixed 256 3 grid with parameters similar to the Cosmic Dawn I simulation performed on Titan at Oak Ridge National Laboratory ), but with several differences which we detail here: (i) the cosmology adopted is compatible with Planck 2015 (Planck Collaboration et al. 2015). We have adopted Ω Λ = 0.693, Ω m = 0.307, Ω b = 0.048, h = 0.677, and the power spectrum is described by σ 8 = 0.8288, with an index n = 0.963, (ii) the star formation efficiency is ǫ ⋆ = 0.033, promoting earlier reionization, (iii) the stellar particle ionizing radiation escape fraction is f esc = 0.1, and (iv) the temperature criterion for star formation has been removed, allowing star formation in all cells with a gas density higher than 50 times the average baryon density δ ⋆ =50 ρ b . In CoDaII, similar parameters allowed us to reproduce a number of observable constraints of the EoR. Therefore we consider them to be adequate for the purpose of this letter. For our 4 h −1 Mpc box, the resulting reionization history is shown as the thick black line in Fig.  1. Although the box successfully reionizes around z=6.5 for our full speed of light run, the timing of reionization is not of great importance for this experiment, since we are more interested in investigating the post-reionization ionized fraction when varying the speed of light.

Reduced speed of light, dual speed of light
The ATON RT module consists of 2 main sub-modules: an advection (i.e. photon propagation) module, and a thermochemistry module, which is purely local and models the interaction between the photons and the hydrogen gas. We will therefore consider the 2 following speeds of light: c pr : the speed of light used in the propagation module, -c ch : the speed of light used in the thermochemistry module.
With these definitions, the usual RSL approximation is obtained by using the same reduced speed of light c rsl =c pr =c ch in both modules. This is the approach of Rosdahl et al. (2013); Trebitsch et al. (2017); Katz et al. (2017) for instance. As we will see, this has a strong impact on the neutral ionized fraction after reionization. We also introduce the dual speed of light approximation (hereafter DSL), where the speed of light c pr is reduced in the propagation module, but maintained to its true value c ch =c in the chemistry module. We performed 1 reference "true" simulation with c pr =c ch =c , 3 RSL simulations with c pr =c ch = 0.2c, 0.05c and 0.01c, respectively. And finally we performed 3 DSL simulations with the true speed of light c ch =1 in the chemistry module but varying the propagation speed c pr =0.2c, 0.05c and 0.01c. The latter will help us understand the shortcomings of the RSL runs. The results are presented in the next section.

Results
In this section we examine the impact of the RSL approximation on the post-reionization neutral fraction for different values of the speed of light.

RSL: reduced speed of light
The reionization histories obtained in the RSL approximation with c rsl =c pr =c ch =c, 0.2c, 0.05c and 0.01c are shown   in the left panel of Fig. 1. The end of reionization, marked by a sudden drop in neutral fraction, is clearly delayed with respect to the true solution (thick black line), although for c rsl =0.1c the lag is rather small. These temporal aspects will be studied in more detail in Deparis et al. (2018, in prep.). Here, instead, we focus on the neutral hydrogen fraction. First of all, our run with the full speed of light yields a neutral fraction in disagreement with the observations of Fan et al. (2006). This is actually commonplace among simulations using the full speed of light (Aubert & Teyssier 2010; Petkova & Springel 2011;So et al. 2014;Bauer et al. 2015;Ocvirk et al. 2016). It is commonly accepted that resolving Lyman-limit systems is necessary to describe properly the population of absorbers in the IGM and predict the correct average neutral hydrogen fraction (Miralda-Escudé et al. 2000). In particular, Rahmati & Schaye (2017) give a requirement of a spatial proper resolution of 1 − 10 kpc. This is compatible with the resolution of our simulations and the Cosmic Dawn simulations in general , where the cell size is 3.38 physical kpc at z=5.5 and smaller at higher redshifts. The lower than observed hydrogen neutral fraction obtained with a full speed of light is therefore not likely to be be due to a lack of spatial resolution. Instead, Mutch et al. (2016); Madau (2017) suggests that both accounting for Lyman-limit systems and some evolution in source emissivity is required to predict average neutral hydrogen fractions in agreement with Fan et al. (2006) after overlap. We find that the best agreement with Fan et al. (2006) is found for the 0.05c run. More importantly, we find that the postreionization average neutral fractions are spread in a wide sequence parametered by the adopted speed of light. In this framework, the true solution is naturally the full speed of light run (thick black line of Fig. 1).
In all RSL runs, the neutral fraction after reionization is too high compared to the true solution, by a factor ∼ 5 (∼ 20, ∼ 100) for c rsl =0.2c (c rsl =0.05c, c rsl =0.01c respectively). We can understand this by considering the analytical expression for the neutral fraction of a hydrogen gas in photo-ionized equilibrium, in the strongly photo-ionized limit (i.e. x HI < 0.01 for instance). In this case, the neutral fraction behaves as: where α B is the case B recombination coefficient, ρ H is the total hydrogen density, ρ i is the ionizing photon density, σ is the ionizing cross-section of hydrogen, and c is the speed of light. This expression shows that, in this regime, the neutral hydrogen fraction follows an inverse linear dependence with the speed of light c: dividing c by a factor 5 results in a neutral fraction increased by the same factor 5 with respect to the true solution using the full speed of light, and this is indeed what we find in the RSL simulations shown in the left panel of Fig. 1. The other terms in this expression are subject to only small variations between our simulations: -The case B recombination coefficient α B (T ) is a function of temperature, which in our post-reionization IGM, is the temperature of post-reionization photo-ionized hydrogen, i.e. 5 000 -20 000 K.
-We find an average post-reionization (z = 5) ionizing photon density ρ i ∼ 0.5 photon per comoving m 3 for all RSL runs, as shown in the left panel of Fig. 2.
These 2 parameters are therefore not likely to cause the order of magnitude offsets seen in the RSL runs. This trend of RSL simulations yielding higher neutral fractions with respect to full speed of light simulations is also seen in the literature, although it has not been clearly identified, explained and reported until now. For instance, this result is similar to the z < 6 region of Fig. 13 of Bauer et al. (2015), hereafter B15, although the authors only comment on the lagging reionization history of their RSL runs. We note that our explanation of the B15 result is different from that given by Gnedin (2016). However, our results are quantitatively in better agreement with B15, and this favors our interpretation over Gnedin (2016).

DSL: dual speed of light
It appears from Eq. 1 that the large error in postreionization equilibrium fraction obtained with RSL may be entirely due to the use of a reduced speed of light in the thermo-chemistry solver, which is similar to artificially reducing the reaction cross-section, or the photon density. Here we confirm this explanation with 3 additional simulations, using a reduced speed of light in the propagation module c pr =0.2c, 0.05c, and 0.01c, but keeping this time the true speed of light c ch =c in the chemistry module. This setup assumes simultaneously 2 different speeds of light, hence the term "dual speed of light" approximation or DSL.
The resulting reionization histories are shown in the right panel of Fig. 1. In stark contrast to the RSL runs, the average neutral fractions after reionization converge this time to the true solution, clearly demonstrating that the reduced speed of light in the chemistry solver is indeed the origin of the large overestimate found with RSL.
Although the DSL approximation produces the correct post-reionization neutral fraction, it may well have its own caveats and artifacts. For instance, a slower propagation speed artificially increases the photon density around sources compared to the true speed of light case. In RSL, this is exactly compensated by a reduced speed of light in the photoionization rate, and this is why RSL yields correct neutral fractions for Stromgren spheres, as shown by Rosdahl et al. (2013). Conversely, in DSL, where the photoionization rate uses the true speed of light, the increased photon density will likely produce over-ionized Stromgren spheres. Moreover, Fig. 1 shows that the lag in reionization redshift is even worse in DSL than in RSL. All DSL runs reionize significantly later than their RSL counterpart. Therefore we can not at this stage advocate for the use of DSL as a replacement of RSL. The latter will overestimate the average neutral hydrogen fraction in the post-reionization photon bath, while DSL will overestimate the ionization of hydrogen before overlap, while Stromgren spheres are still growing under the influence of local sources. Our simulations will be further analyzed to understand and illustrate the merits and shortcomings of the two approaches in both phases in a future paper.

Conclusions
We have investigated the impact of the RSL approximation on the properties of the intergalactic medium using fully coupled radiation-hydrodynamics RAMSES-CUDATON simulations. We find that reducing the speed of light by a factor 5 (20, 100) leads to overestimating the postreionization neutral hydrogen fraction by the same factor ∼5 (20, 100, respectively) with respect to our reference simulation using the full speed of light. We consider the simple analytical expression for a hydrogen gas in photo-ionization equilibrium to show that the error is driven by the hydrogen -photon chemistry: reducing the speed of light has the same effect as artificially reducing the reaction cross-section or reducing the photon density, and results in underestimating the ionizing flux. We confirm this interpretation by running additional simulations using a reduced speed of light in the propagation module, but keeping this time the true speed of light in the chemistry module. With this setup, dubbed "dual speed of light" because of the simultaneous use of 2 different speeds of light, the post-reionization neutral hydrogen fractions converge to the true value, which validates our explanation.
Thus, we have exposed an overlooked, fundamental issue of the reduced speed of light approximation, which should be kept in mind when interpreting results obtained in this framework.
Indeed, the neutral hydrogen fraction is an extremely important parameter governing a whole chain of processes, such as gas cooling (as long as metals do not become the main coolant), and therefore star formation, but also possibly the ionizing escape fraction of galaxies, and hence the geometry and the timing of reionization. Therefore, errors in computing the neutral fraction such as reported here for the reduced speed of light approximation could partially offset this chain of processes and yield inexact predictions. Later on, after reionization, the neutral fraction governs the properties of the Lyman-α forest and may impact the properties of simulated Lyman-α emitters as well. However, in moderately ionized regions (x HI > 0.1) and for a given photon density, the dependency between x HI and the adopted speed of light c becomes shallower, and in quasi-neutral regions, the equilibrium neutral fraction may become insensitive to the adopted speed of light. Therefore, depending on the exact physical properties of the main absorbers responsible for the Lyman-α forest, they may or may not be affected by the RSL approximation. Further simulations are needed to investigate this issue in detail, which will be the topic of a future paper. As a major goal for the community, it is of prime importance that RHD codes achieve consistent predictions of the average neutral fraction, during and after the end of reionization. Also, since the neutral hydrogen fraction is such an important parameter, we urge all colleagues and authors to always show this quantity rather than, or along with, the average ionized fractions. Moreover, authors should be very careful to properly spell out the type of reduced speed of light approximation used, and its implementation in the photon propagation steps as well as in the thermochemistry module. This information is essential to understand and gauge the significance of the flurry of studies published these recent years and ongoing.