Impact of the reduced speed of light approximation on ionization front velocities in cosmological simulations of the epoch of reionization

Coupled radiative-hydrodynamics simulations of the epoch of reionization aim to reproduce the propagation of ionization fronts during the transition before the overlap of HII regions. Many of these simulations use moment-based methods to track radiative transfer processes using explicit solvers and are therefore subject to strict stability conditions regarding the speed of light, which implies a great computational cost. It can be reduced by assuming a reduced speed of light, and this approximation is now widely used to produce large-scale simulations of reionization. We introduce a new method for estimating and comparing the ionization front speeds based on maps of the reionization redshifts. We applied it to a set of cosmological simulations of the reionization using a set of reduced speeds of light, and measured the evolution of the ionization front speeds during the reionization process. We find that ionization fronts progress via a two-stage process, the first stage at low velocity as the fronts emerge from high density regions and a second later stage just before the overlap, during which front speeds increase close to the speed of light. Using a set of small 8Mpc/h^3 simulations, we find that a minimal velocity of 0.3c is able to model these two stages in this specific context without significant impact. Values as low as 0.05c can model the first low velocity stage, but limit the acceleration at later times. Lower values modify the distribution of front speeds at all times. Using larger 64Mpc/h^3 volumes that better account for distant sources, we find that reduced speed of light has a greater impact on reionization times and front speeds in underdense regions that are reionized at late times. The same quantities measured in dense regions with slow fronts are less sensitive to c values.


Introduction
The reionization of the Universe is driven by the propagation of ionization fronts created by the first astrophysical light sources. This process is not homogeneous as it reflects the distribution of absorbers and sources created by the rise of large-scale structures: the typical associated picture consists of a network of HII regions that eventually overlap at the end of reionization, by z ∼ 6.
As a consequence, great efforts are currently being made to model this complex process using cosmological simulations by including radiative transfer physics. Many of these works rely on a moment-based description of radiative transfer (see, e.g., Bauer et al. 2015;Gnedin 2014;Ocvirk et al. 2016;Rosdahl et al. 2018;Aubert et al. 2018), which is described as a set of conservative equations on radiative quantities (e.g., radiative density, radiative flux, radiative pressure). In many aspects, this technique is akin to modeling a "photon fluid".
A practical advantage of this description is the possibility to use ready-made generic methods or existing modules developed for hydrodynamics. However, moment-based methods have to assume a finite speed of light as they model the actual propa-gation of radiation density. Since this fluid description is often handled through an explicit approach, this choice of methodology is not without consequences (even though implicit solvers do exist; see, e.g., González et al. 2007). The main constraint for explicit moment-based solvers is related to the Courant-Friedrichs-Lewy (CFL) condition. It imposes that the higher the velocity of a process is, the higher the temporal resolution has to be to guarantee numerical stability (see, e.g., Aubert & Teyssier 2008;Rosdahl et al. 2013). The time step ∆t must satisfy where C CFL < 1 (in 1D), ∆x is the spatial resolution, and v is the maximum speed of the considered process. For high velocities, and in the most extreme case for light speed velocities, such solvers are conditionally stable only if the time resolution is consequent, leading to a large number of time steps and a great computational cost. Several methodologies have emerged to reduce this cost. The first consists in taking advantage of the recent evolution in hardware; for example, using graphical processing units (GPUs) can divide the computational time by almost two orders of magnitude due to their highly parallel capability (Aubert & Teyssier 2010). The second method is the subject of this study: the reduced speed of light approximation (RSLA). The RSLA considers that light propagates at a fraction of its actual speed v =c · c, where c is the physical speed of light (299 792 458 m s −1 ) and c ≤ 1 is the RSLA factor. The main idea of this technique relies on the fact that the typical velocities of dynamical processes (between 100 and 10 4 km s −1 ) can be several orders of magnitude lower than the physical speed of light. As the speed of light is significantly greater than any other speed, lowering its value in a simulation should not have a significant impact, but should allow a comfortable gain on the computational cost. For example, if the reionization properties remain similar while usingc = 1 orc = 0.1, a factor of ten in the computation of radiative processes can be gained just by reducing the number of radiative time steps. This gain can be appreciable as radiation represents approximately 80% of the computational cost in this kind of simulation whenc = 1. This is the approach used by Katz et al. (2017), among others, where different values ofc are used in different regimes of density within a given experiment, thus capturing typical front speeds in the relevant regimes while remaining computationally reasonable (see also Rosdahl et al. 2018). However, the consequences of such an approximation and their dependencies (e.g., on implementation, resolution, or box size) are still debated (Gnedin & Abel 2001;Gnedin 2016;Bauer et al. 2015;Rosdahl et al. 2018, Ocvirk et al., in prep.).
In this study, we investigate the validity of this approximation and we explore how reducing the speed of light influences the evolution of ionization fronts during the reionization epoch. First, we present a set of radiative hydrodynamics (RHD) simulations with different speeds of light, and introduce a method for computing the ionization front (hereafter I-front) speeds using reionization redshift maps. Then, using different diagnostics (average ionization history, ionization maps, front speeds), we quantify the impact of RSLA on the global reionization.

Methodology
Simulations were produced with EMMA (Electromagnetisme et Mécanique sur Maille Adaptative) an adaptive mesh refinement (AMR) cosmological code with fully coupled RHD . The light was considered a fluid; its propagation was tracked using the moment-based M1 approximation (Levermore 1984;González et al. 2007;Aubert & Teyssier 2008). We ran a set of small simulations at a similar resolution of production runs such as CODA (Ocvirk et al. 2016) or CROC (Gnedin 2014). We considered 8 h −1 cMpc ∼ 12 cMpc 3 volumes, simulated from redshift z = 150 to the end of the reionization (5 < z < 6). Dark matter was resolved with 256 3 particles with a mass of 3.4 × 10 6 M . Hydrodynamics and radiation were simulated on a grid of 256 3 resolution elements for a coarsest spatial resolution of 46 ckpc. The grid was refined according to a semi-Lagrangian scheme when a cell contained a mass greater than eight dark matter particles; the refinement was not allowed if the spatial resolution of the newly formed cells was lower than 500 pc. Our stellar mass resolution was set to be equal to 7.2 × 10 4 M . Our emissivity model follows Starburst99, with a population of 10 6 M having a Top-Heavy initial mass function and a Z = 0.001 metallicity. To limit the number of unknowns, we did not consider supernovae feedback and the corresponding module was turned off. All simulations were run with the exact same parameters (including star formation) and initial conditions except forc, the ratio of the simulated speed of light to the real speed of light. The boundary conditions were assumed to be periodic for all processes, including radiative transfer. Radiation from source outside the volume was therefore modeled implicitly by exact replicates of the inner distribution of sources. We assess the impact of a less trivial contribution of external sources in Sect. 4. We did not consider the influence of an additional homogeneous background.
We ran a set of six simulations withc = 1,c = 0.3,c = 0.1, c = 0.05,c = 0.02, andc = 0.01. Thec value was used to solve the set of conservative equations of the M1 approximation (on radiative densities and fluxes) in all regimes and in the calculations of photoionization rates; we followed the same methodology as Rosdahl et al. (2013) and Bauer et al. (2015). We note that this methodology differs from that of Gnedin (2016), who distinguishes between the radiation background and the fluctuations superimposed. In this casec is only used to propagate the fluctuations. Figure 1 presents the cosmic star formation histories (SFH) for differentc values. This measure shows that the RSLA does not change the cosmic SFH. As ionizing photons are emitted by newly formed stars, the photons budget is not impacted by the RSLA and all the simulations have to propagate, at least on average, the same quantity of radiative energy. Figure 1 also presents the average volume weighted ionization state as a function of time X V (t)

Cosmic star formation and ionization histories
with the local ionization fraction We observe a direct link between the ionization history and the RSLA. The slower the light is, the later the reionization occurs. We define the reionization redshift as the redshift when the volume weighted hydrogen neutral fraction decreases below X V = 10 −4 . At these times the percolation process of HII regions (also known as the overlap) is complete. Figure 1 presents this reionization redshift as a function of the reduced value of the speed of light. Reionization redshifts converge asc is close to unity and the slope flattening indicates that anyc > 0.1 leads to similar instants of reionization. The other way round, reionization redshifts rapidly decline with decreasingc values. For instance,c = 0.1 presents a delay of ≈60 Myr and less than 0.5 in redshift, but these values are extended to ≈425 Myr and almost 2.2 in redshift for thec = 0.01 run. The different reionization histories obtained here suggest that I-fronts may at some point become as fast as the assumed speed of light, and the RSLA may therefore have an impact on simulation predictions. We assess this point directly in the next sections.

Ionization maps
Ionization redshift maps offer an interesting tool for analyzing the whole reionization process using a single data field. Redshift maps are obtained by keeping in memory, for each cell, the time when its ionization fraction crosses a given threshold. The maps were computed on the fly at each simulation time step to get the highest temporal resolution possible; this is in contrast to similar studies where such maps are computed in post-processing from snapshots with a sparser time sampling (see, e.g., Ocvirk et al. 2013). In this study, the cell reionization time is considered to be the first time when the volume weighted ionization fraction goes above 50%. Figure 2 shows three maps of reionization redshifts for threec values, and the reionization redshift probability den-sity function (PDF) can be found in Fig. 3 for all runs. Maps are made from slices with a thickness corresponding to one coarse cell (46 ckpc), taken at the same coordinate in the z-axis for the three simulations. The slice is chosen to contain the first cell to cross the ionization threshold in thec = 0.1 simulation.
Qualitatively, the reionization maps forc = 0.01, 0.1, 1 present common features. The sources are located at the same places, consistent with the fact that the ionizing feedback does not significantly change global star formation processes in our simulations ( Fig. 1). Radiation escapes high density regions in comparable butterfly shapes when it propagates freely in directions that offer less resistance perpendicular to filamentary regions. Thec = 0.1 run is quite similar to thec = 1. Dark green isocontours (z ≈ 6) are located at the same positions and have the same shape and extension. Greater differences can be noted at smaller redshifts when radiation reaches the underdense regions (in blue); these voids reionize at later times forc = 0.1 compared toc = 1, as expected from lower velocities. When the same map forc = 0.01 is considered the differences are striking: the underdense regions, distant from the main sites of photon production, are reionized at much later times compared to the two previous experiments with higher velocities.
The quantitative difference is particularly clear in the redshift PDFs (see Fig. 3): at high redshifts all PDFs are similar, but the differences become more notable at later times (starting from z ∼ 11), and the cumulative effect leads to a significant delay in the global reionization for lowc. This is in agreement with the measurements made on the volume weighted ionization history ( Fig. 1): runs with high speed of light reionize earlier because radiation is able to ionize voids more rapidly.

Ionization front speeds
The reionization maps represent a time at each point of the simulation in each cell. Using the spatial gradient of this map, we can estimate the time required to reionize each cell per unit length. We use a discrete gradient where i, j, k are cartesian cell indexes and ∆x is the 1D spatial extent of a reionization map cell. This gradient represents the time needed to ionize a given distance (e.g., in [yr pc −1 ]), the inverse is therefore equivalent to a velocity (e.g., in [pc yr −1 ]) : Here V reio can be interpreted as an estimator of the speed of the I-front passing through a cell. Figure 4 shows maps of I-front speeds obtained with this method for three runs with different reduced speeds of light. In the maps presented in Fig. 4, the concentric "rings" with high and low speeds fronts are induced by the consecutive events of star formation. I-fronts can only expand if there are internal radiative sources, but star formation is not continuous and happens in a bursty manner. Fronts slow down during the period between two consecutive star formation events. We investigate this effect further in the appendix.
It should be noted that the velocity obtained using Eqs. (5) and (6) is only an approximate estimator; it produces an adequate measure in the case of a perfect plane parallel front, although it can be less accurate in the more realistic case of multiple, multi-directional fronts. In particular, limitations appear ciated reionization speeds become artificially high and even infinite in the case of neighboring cells with identical reionization times. Example of such situations are as follows: -In overdense regions where adjacent cells have simultaneous star forming events. -In underdense regions where two (or more) I-fronts start to overlap. As a consequence, measured front speeds can be up to several times greater than the speed of light in the simulation. In the PDFs of I-front speeds shown in Fig. 5, we observe that limited speeds of light do not prevent the presence of I-front speeds greater than their values. For instance, there are still front speeds of 10 −1 whenc = 10 −2 , but these values represent a small integrated fraction of the PDFs (∼8 × 10 −5 at most in all our experiments), thanks in particular to the important time resolution of the reionization maps. Therefore, the probability of finding velocities greater thanc is low and the impact of the reduced speed of light on front propagations can be measured using this estimator. Furthermore, our focus is on the relative differences of different simulations obtained via identical methodology rather than measures of absolute values for the velocities.

I-front speed as a function of redshift
The results from the two previous sections can be combined: at each map position we have access to a reionization time/redshift and a front speed. Therefore, front speeds can be binned according to the associated reionization redshift, providing the distributions of front speeds as a function of redshift. Figure 6 presents this distribution and Fig. 7 the evolution of the average front speeds for differentc.
We first focus on the 2D histograms that show the redshift evolution of the distribution of front speeds (see Fig. 6). In thẽ c = 1 run, we observe that for z > 8, the distribution of speeds remains between ≈10 −1 c and ≈10 −4 c. This phase is then followed by a quick acceleration "plume" with velocities close to c. We therefore observe a two-stage process: first, light escapes overdense regions at an average quasi-constant moderate speed and second, the percolation of ionized bubbles allows radiation A142, page 4 of 10 to reach underdense regions. In these regions during these later stages, the front speeds are not significantly limited by the inter- action with matter. The light is free to fill the voids and front speeds are able to reach high values.
In thec = 0.1 run, the first phase is not impacted by the RSLA, which is expected since the front speeds are already lower than ≈ 10 −1 c for thec = 1 run. The main differences appear in the second phase: the z < 8 peak is reduced by the RSLA, whereas thec = 1 experiment presents fronts that are able to reach speeds greater than 0.1c. By settingc = 0.1, the velocities of these fronts are limited to ≈0.1c in an artificial manner: the end of the reionization is slightly delayed. At low speed of light values, in thec = 0.01 run, the differences are even more drastic. Not only are front speeds severely limited during the late accelerating phase, but a difference can also be seen in the first phase (z > 8) as extreme front speed values are artificially decreased by an extreme RSLA. It leads to an important cumulative delay in the reionization redshift.
Considering now the average front speeds as a function of redshift (see Fig. 7), we clearly distinguish the two stages: an almost constant value at early times, followed by an acceleration for redshift z < 8. In the z > 8 quasi-constant phase, the average front speed goes from 4.1 × 10 −3 for thec = 1 run to 2.2 × 10 −3 for thec = 0.01 run. For reduced speed of light values lower thañ c ∼ 0.05, the impact seems strong enough to shift down the average front speed value compared to other experiments for z < 10. The late-stage acceleration for redshifts z < 8 is significantly impacted by the RSLA, and as expected a high speed of light value leads to a more sudden "flash" of the volume at the end of the reionization than with lower values. This final acceleration corresponds to the light reaching underdense regions, and in all cases the average front speed is limited by a reduced speed of light, even thoughc = 0.3 is quite consistent withc = 1. To sum up, reduced speed of lightsc > 0.05 can be considered as satisfying to reproduce the first "slow" stage at quasi-constant speed for front propagations, but will likely induce a delay for the overlap ifc < 0.3.
These findings also provide additional insights into the influence (or lack of influence) of RSLA as a function of local matter distribution by relating these two stages to two different regimes of matter density. The first phase is linked to the initial stages of the production of radiation by the first sources in dense regions.  In this regime the RSLA can be a good approximation down toc = 0.05. By extension, we can assume that speeds of light as low as this value can be used to study overdense objects at any time, for example to study the escape fraction of ionizing radiation of galaxies. The second phase with high I-front speeds is related to the overlap of HII regions as radiation propagates in the more diffuse intergalactic medium (IGM) or in voids: in this regime of densities, when radiation propagates freely, I-front speeds can be as high as c and a reduced speed of light is always a limiting factor. For instance, the remote influence of one galaxy on another distant one, via the propagation of radiation through diffuse matter, is likely to be improperly timed at any time.
Numerical parameters such as the box size or the spatial resolution are likely to impact our quantitative predictions, due to this local density influence. Increasing the resolution should lead to denser regions and potentially slower I-fronts in the first phase. Increasing the box size provides greater underdense regions or voids and could induce faster I-fronts in the second accelerating phase. Overall, further investigations are needed to assess the density or environmental dependence of I-front speeds and how the RSLA impacts the prediction of simulations: we initiate these investigations in the next section.

Probing different regimes of volumes and densities
In order to probe different regimes of volumes or densities, we produced an additional set of reionization simulations with a larger volume (64 Mpc h −1 ) 3 and a lower mass resolution (512 3 dark matter particles with a 2×10 8 M mass and 512 3 cells in the initial grid + 5 refinement levels). Reionization maps and front velocities were computed in an identical fashion, thus providing front characteristics in regimes of greater volumes and lower resolution. Light speedsc = 1 andc = 0.1 were used, bracketing the regimes where convergence was achieved in the previous smaller simulations. The stellar emission was set to reionize the full volume by z = 6 for thec = 0.1 model; all the other parameters were kept identical to the main set of models discussed previously. With this setup, thec = 1 simulation is reionized at z = 7; both simulations (64 Mpc h −1 ) 3 reionize later than their small box equivalents using the samec values. The redshift evolutions of the front speed distributions are shown in Fig. 8. Overall, this different regime of volume and resolution leads to a global behavior similar to that encountered previously: velocities evolve according to a two-stage process with an initial slow-front regime followed by an increase as the overlap of HII regions take place. In both cases the average front velocities result in values close to the speed of light, and in both cases values higher thanc can be measured, due to the same limitations A142, page 6 of 10 of the methodology discussed previously (see Sect. 3.3). Also, thẽ c = 1 model experienced an earlier reionization compared to thẽ c = 0.1 model, as already observed for the smaller simulations (see Fig. 1). One notable difference compared to the small volume models can be noted in the initial slow front stage: velocities are higher in the current, less resolved simulations. Notably, local velocities as low as 10 −4 c cannot be measured in the two simulations, whereas such values were found in the fiducial set of smaller simulations with higher mass resolution. Likewise, the typical average front velocity is ∼10 −2 , higher than previously found. This regime is driven by the densest regions around the first sources and the lack of high density contrast in these large volume simulations allows the fronts to experience less "resistance" as they escape into the more diffuse IGM.
At later stages close to the overlap, front velocities tend to be higher for both values of the speed of light compared to the simulations made in smaller volume. In fact, spurious values greater thanc are more frequent in both cases, even though they still represent a small fraction of the total number of cells (<0.001% in both cases). High values of front velocities, close toc, are found at earlier times compared to the smaller boxes, even though reionization takes place later; the late stage evolution is shallower and less sudden in the 64 Mpc h −1 simulations than models produced in 8 Mpc h −1 boxes. Large underdense regions allow fast propagation at larger redshifts through a more diffuse IGM.
Thanks to the 64 Mpc h −1 volume probed here, we can also assess the influence of external sources on the previous smallbox experiments. For bothc values, we split the (64 Mpc h −1 ) 3 volumes into 8 × 8 × 8 = 512 smaller (8 Mpc h −1 ) 3 cubes. Each of these smaller subvolumes would experience a reionization driven by the properties of its own sources, but also influenced by external fronts created by external sources.
Since bothc = 1 andc = 0.1 simulations share the same set of initial conditions, such subvolumes can be compared from one simulation to another. The comparisons are shown in Fig. 9. First, we compare the median reionization time of each subvol- ume for both values of the speed of light. These times present a strong correlation: the volume that reionizes the earliest in one simulation also reionizes at early times in the other. The progression of the reionization is driven by the buildup of structures and sources, a process only weakly modified by the speed of light. In addition, it should be noted that the densest regions (the density of each subvolume being color-coded) are reionized first, as expected if sources tend to appear first in these regions, whereas the most underdense are reionized last. Despite the tight correlation, an offset from the 1:1 correspondence can be seen between the subvolume reionization times of thec = 1 andc = 0.1 simulations. This offset is not systematic, as can be seen when comparing the two simulations on a cell-by-cell basis (see Appendix A).
It results instead from the presence of regions that experience a delayed reionization in all the (8 Mpc h −1 ) subvolumes.
In addition, a departure from a linear correlation can be seen for the volumes that reionize last, corresponding to the volumes that are the most underdense. In these cases, the reionization delay of the reduced speed of light simulation is more significant than that observed in denser subvolumes. In the bottom panel of Fig. 9 we show the comparison of median front velocities, taken from the distribution of velocities at all times, in the two simulations. Volumes that reionize the earliest (and therefore the densest ones) present median front speeds that differ only slightly when comparing thec = 1 andc = 0.1 experiments. On the other hand, volumes that reionize the latest present a larger departure from the 1:1 correspondence; underdense regions promote high velocities through a diffuse IGM, and are therefore more sensitive to the choice ofc.
This large sample of subvolumes gives us some insight into the impact of cosmic variance, finite-volume effects, and external sources. Dense regions are likely to host sources at early times that create slow fronts, and such volumes are self-reionized. Another possibility could be that these regions are swept by slow fronts created in even denser, nearby regions. In this regime of density, we find that ac = 0.1 speed of light has only a moderate impact on front velocities and reionization times. On the other hand, underdense regions are reionized at later times and are therefore likely to be externally reionized by sweeping fronts; in a reduced speed of light regime, such fronts would arrive at later times and be slower. In such cases, a high speed of light would be necessary. The transition between these two regimes takes place approximately at the average cosmic density that was enforced by construction in our previous set of 8 Mpc h −1 simulations; a conservative stance would therefore push for a highc if external fronts have to be taken into account in this specific case. In any case, and as demonstrated here, effects of cosmic variance, finite volume, or density regimes can be quite important in general.

Conclusions
We described a method used to estimate ionization front speeds from reionization maps produced by cosmological simulations. With this method we find that reducing the speed of light modifies the propagation of fronts in our simulations.
We find that the propagation of radiation in our models of reionization is a two-stage process: a first phase with quasiconstant, moderate average I-front speeds, and a second stage with an acceleration phase during which the average I-front speed increases greatly, close to the speed of light. Using a set of small (8 Mpc h −1 ) 3 experiments we found the following: -During the early constant speed phase, the average value does not depend onc forc > 0.05. -The second acceleration phase is always impacted by a reduction in the speed of light, but convergence can be achieved forc ∼ 0.3 to reasonably reproduce the redshift evolution of average I-front speeds. These quantitative prescriptions are likely to depend on the type of simulation (resolution, volume, external sources). In order to get some insight on the robustness of these results, we performed the same analysis on a set of similar (8 Mpc h −1 ) 3 volumes extracted from a larger (64 Mpc h −1 ) 3 simulation with a lower mass resolution. We find that in the slow phase there are greater front speeds, but still well below 0.1c, and that large variations in median front velocities and reionization times can be found in underdense regions that experience late reionization when comparing thec = 1 andc = 0.1 experiments. The influ-ence of fronts created by distant sources is therefore likely to be impacted by the choice of speed of light.
Overall, our finding shows that quantitative prescriptions on the speed of light can be quite sensitive to the experimental setup. In addition the inclusion of bright sources creating large-scale fluctuations in the radiation field, such as quasars (see, e.g., Chardin et al. 2015) or X-rays with large mean free paths (see, e.g., Baek et al. 2010), may change the quantitative results obtained here on front velocities. We also did not consider the impact of an additional homogeneous radiation background, which could affect propagation in the late-reionized, low density, regions, i.e., the regimes where the greatest discrepancies are being measured for different values ofc. Even though this lack of background is consistent with many previous studies made on such scales where the reionization is driven by the population of simulated sources, such a background could dramatically lower the impact of reduced speed of light by controlling the ionization state of vast regions. Our conclusions should therefore be considered in the strict, but widely accepted context of reionization models driven by the "simulated sources".
The two-stage evolution of front velocities is likely to be a robust finding, as demonstrated for example by recent similar findings by D' Aloisio et al. (2018), and could potentially be used to define domains of validity of reduced speed of light experiments, depending on the sizes, densities, or redshifts probed. However, beyond its impact on the propagation of fronts, nontrivial effects of the reduced speed of light could exist, for example on the photo-suppression of distant star formation, the relative delays between different kind of feedbacks (see, e.g., Trebitsch et al. 2017), or the gas thermochemistry (see, e.g., Ocvirk et al. 2018). Systematic studies have yet to be generalized.

Appendix A: Impact of the speed of light on the front speed spatial distribution
As discussed in the main text, ripple-like features can be seen in the spatial distribution of front speeds, which we attribute to imprints of the episodic nature of photon production by the sources. In this appendix, we briefly discuss this effect and how it may be affected by reduced values of the speed of light. In Fig. A.1 we present four different experiments made using the same setup as described in Sect. 4, consisting of reionization simulations in a large volume (64 Mpc h −1 , 512 3 base resolution), withc ∼ 1 andc ∼ 0.1. Each pair of simulations was run using two different recipes of photon production from the stellar particles. The first recipe follows the Starburst99 model described in Sect. 4, leading to a steady production rate of ionizing photons for ∼3 Myr before a sharp decline. The second recipe spreads the initial steady production rate of photons over 10 Myr, while keeping the same total amount of produced photons. Therefore, we end up with a set of four simulations consisting of two values for the speed of light and two types of photon production: "short and bright" and "long and faint". Figure A.1 presents the front velocity spatial distribution around the same source in the four different setups. Figure A.2 shows front velocity profiles extracted from the same distributions.
In the four configurations the typical ripples can be observed in the velocity field. At large distances from the source, fronts can achieve greater speeds in thec = 1 simulations, presented as lighter regions on the maps, as expected. If we first consider the reduced speed of light experiment,c = 0.1, the episodic nature of the photon production leaves a different imprint, depending on its parameters. The short and the long source scenarios present the same overall spatial distribution of front velocities, however the ripples are smeared out in the case of a "long and faint" production of photon. Conversely, the "short and bright" model presents a greater number of ripples with more contrast; the bursty episodes of photon production can be detected in the structure of the velocity field. The difference is particularly striking in the innermost part of the ionized region. Meanwhile the same innermost parts in thec = 1 experiments do not seem to be impacted by the difference in photon production history, and present very similar profiles for both types of sources (see Fig. A.2), which suggests that for this peculiar source the delay between two trains of photons is greater than 10 Myr, because otherwise the ripples would be smeared out.
However, it can be seen that velocities are quite similar in the four experiments at small distances (see Fig. A.2), with a measured velocity close to 0.005c <c, and is therefore independent of the choice ofc. Ripples in the velocity field are fluctuations around low values rather than strong peaks. This suggests that the visibility of the ripples, which in all cases are seen in the slow progression regime, is not related to differences in long-distance propagations.
Overall, in the current experiments, the differences induced by the source models remain marginal and do not have a strong impact on reionization times distributions. We do not exclude that there might be regimes where the cumulative effects of these local differences in front velocities could lead to strong divergence in the front propagation on large scales, for different models of photon production, and when reduced values ofc are used, which would require specific experiments to be demonstrated. However, since most of the observed differences take place in the slow-front regime, we think this lack of convergence is unlikely for reasonable stellar models.

Appendix B: Cell-by-cell comparison of reionization times and front velocities with different speeds of light
In this section we briefly describe the raw data that led to Fig. 9. Instead of showing the reionization times and front velocities averaged in (8 Mpc h −1 ) 3 boxes extracted from large (64 Mpc h −1 ) 3 simulations, Fig. B.1 compares the same quantities on a cell-by-cell basis for thec = 0.1 andc = 1.0 large simulations. In particular, it should be noted that no systematic offset from the 1:1 correspondence is observed; the tendency of quantities averaged over a (8 Mpc h −1 ) 3 to drift away from this correspondence in Fig. 9 is simply due to the systematic presence of underdense regions within that experience delayed reionizations whenc = 0.1 is used.