The reduced speed of light approximation has been employed to speed up radiative transfer simulations of reionization by a factor of gtrsim 510. However it has been shown to cause significant errors in the HIionizing background near reionizations end in simulations of representative cosmological volumes. This can bias inferences on the galaxy ionizing emissivity required to match observables such as the Lyalpha forest. In this work we show that using a reduced speed of light is to a good approximation equivalent to rescaling the global ionizing emissivity in a redshiftdependent way. We derive this rescaling and show that it can be used to correct the emissivity in reduced speed of light simulations. This approach of rescaling the emissivity after the simulation has been run is useful in contexts where the emissivity is a free parameter. We test our method by running full speed of light simulations using these rescaled emissivities and comparing them with their reduced speed of light counterparts. We find that for reduced speeds of light tildec geq 0.2 the 21 cm power spectrum at 0.1 leq k hrm Mpc1 leq 0.2 and key Lyalpha forest observables agree to within 20 throughout reionization and often better than 10. Positiondependent timedelay effects cause inaccuracies in reionizations morphology on large scales that produce errors up to a factor of 2 for tildec leq 0.1. Our method enables a factor of 5 speedup of radiative transfer simulations of reionization in situations where the emissivity can be treated as a free parameter.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
Accurate simulations of reionization using the reduced speed of light approximation
{'Christopher Cain a': 'a School of Earth and Space exploration, Arizona State University, Tempe, AZ 85281, USA  \nE-mail: clcain3@asu.edu', 'Abstract.': "The reduced speed of light approximation has been employed to speed up radiative transfer simulations of reionization by a factor of ≳ 5 -10 . However, it has been shown to cause significant errors in the HI-ionizing background near reionization's end in simulations of representative cosmological volumes. This can bias inferences on the galaxy ionizing emissivity required to match observables, such as the Ly α forest. In this work, we show that using a reduced speed of light is, to a good approximation, equivalent to re-scaling the global ionizing emissivity in a redshift-dependent way. We derive this re-scaling and show that it can be used to 'correct' the emissivity in reduced speed of light simulations. This approach of re-scaling the emissivity after the simulation has been run is useful in contexts where the emissivity is a free parameter. We test our method by running full speed of light simulations using these re-scaled emissivities and comparing them with their reduced speed of light counterparts. We find that for reduced speeds of light ˜ c ≥ 0 . 2 , the 21 cm power spectrum at 0 . 1 ≤ k/ [ h Mpc -1 ] ≤ 0 . 2 and key Ly α forest observables agree to within 20% throughout reionization, and often better than 10% . Position-dependent time-delay effects cause inaccuracies in reionization's morphology on large scales that produce errors up to a factor of 2 for ˜ c ≤ 0 . 1 . Our method enables a factor of 5 speedup of radiative transfer simulations of reionization in situations where the emissivity can be treated as a free parameter.", '1 Introduction': "Radiative transfer (RT) simulations are a crucial tool for studying the Epoch of Reionization (EoR). During this period, the first generation of hydrogen ionizing sources drove ionization fronts (I-fronts) into the neutral intergalactic medium (IGM), which later overlapped and produced a fully ionized IGM by z = 5 -6 [1-3]. RT simulations can self-consistently track the growth and overlap of ionized regions during reionization, and model in detail the physical conditions in the intergalactic medium (IGM) during the process [4-7]. They also allow the physics of galaxies to be directly linked to IGM conditions, enabling observables that probe reionization to constrain the physics of high-redshift galaxies.  \nDespite being the most physically accurate method available to study reionization, RT simulations are also the most computationally expensive. Formally, the RT equation is a 7 -dimensional problem - three spatial, two angular, one frequency, and one time. Moreover, representative cosmological volumes ( ≳ [100Mpc] 3 ) typically contain millions of ionizing sources, further increasing the cost of RT. Several methods, including moment-based RT [8-10] and adaptive ray tracing [11, 12] have been developed to speed up the angular component of the problem and eliminate the direct scaling of computation time with the number of sources. For many applications, only one to a few frequency bins are required to accurately model IGM conditions [e.g. 6, 13, 14]. Sub-grid treatments of the clumping of the IGM on ∼ kpc scales [e.g. 15, 16] can allow for simulations to be run at relatively coarse spatial resolution, which helps to counter the computational scaling with number of spatial dimensions.  \nAn common trick used to speed up the time-dependence of the RT equation is the reduced speed of light approximation (RSLA). The time step required to resolve the transfer  \nof radiation on an Eulerian grid is  \nt rad ≤ ∆ x cell c (1.1)  \nwhere ∆ x cell is the cell width and c is the speed of light. Often, t rad is much smaller than any other timescale (such as hydro-dynamical or chemical timescales), meaning that it could be much larger without affecting the treatment of these other processes. Indeed, if the RT is done in post-processing (as in many applications), the only relevant timescales are those associated with the RT equation itself. In addition to t rad , these are (1) the cell-crossing time of ionization fronts (I-fronts) and (2) the timescale over which the source population and underlying density field evolve. The typical speed of I-fronts ( v IF ) is 1 -2 orders of magnitude lower than c [17], and the galaxy population evolves on timescales of 10 s -100 s of Myr, much longer than typical values of t rad ( ≲ 1 Myr).  \nThe RSLA introduces a smaller speed of light, ˜ c < c , chosen as small as possible so that t rad remains the smallest timescale in the problem. Because the number of time-steps is inversely proportional to t rad , this gives a direct computational speedup by a factor of c/ ˜ c . The RSLA has been used in a number of contexts to dramatically increase computational efficiency of reionization RT simulations at various spatial scales [18, 19].  \nHowever, several recent works [20-23] have noted that the RSLA fails in certain important regimes of reionization modeling. Ref. [20] argued that the RSLA should be used only to model local fluctuations in the ionizing background during reionization, but never the mean ionizing background itself. Supporting this point, Ref. [22] found that using the RSLA leads to an under-estimate of the ionizing background close to and after the end of reionization by as much as a factor of ˜ c/c . Most recently, Ref. [23] found that reionization models calibrated to reproduce recent Ly α forest observations at 5 < z < 6 (from Ref. [24]) using the RSLA yielded biased results for the galaxy ionizing emissivity required to match these measurements. In this work, we seek to understand the reasons for these issues and find a way to account for them.  \nThe findings of Ref. [23] hint at a possible pathway to a solution. They found that in simulations with different values of ˜ c , the global ionizing emissivity ˙ N em ( z ) could be calibrated to reproduce measurements of the Ly α forest transmission at z ≤ 6 (their Fig. 10). This suggests that perhaps the RSLA is, at least approximately, equivalent to a redshift-dependent re-scaling of ˙ N em . In this work, we will show that this is indeed the case, and we will provide a formula for this re-scaling that can be applied to ˙ N em for any simulation that uses the RSLA. In situations where ˙ N em is being set by hand (as in Ref. [23]), or is a function of some set of free parameters, this re-scaling can be applied to ˙ N em in 'post-processing'. Then, the simulation can be treated as if it had used the full speed of light and the re-scaled ˙ N em .  \nThis work is organized as follows: in §2, we use a simple analytic model to understand why the RSLA fails in certain situations, and derive the aforementioned re-scaling of ˙ N em . In §3, we briefly describe our numerical methods for running RT simulations of reionization. In §4-5 we use simulations to test how well our approach works, and we conclude in §6. Throughout, we assume the following cosmological parameters: Ω m = 0 . 305 , Ω Λ = 1 -Ω m , Ω b = 0 . 048 , h = 0 . 68 , n s = 0 . 9667 and σ 8 = 0 . 82 , consistent with Ref. [25] results. All distances are quoted in co-moving units unless otherwise specified. Throughout, we will work in units in which the speed of light is unity ( c = 1 ).", '2.1 When and why does the RSLA fail?': "Wewill first build intuition for when and why the RSLA fails using a simple analytic argument. In the early stages of reionization, the universe contains a collection of ionized 'bubbles' surrounding isolated clusters of ionizing sources, and is otherwise fully neutral. The mean photo-ionization rate, Γ HI , inside one of these bubbles is  \nΓ HI = N γ cσ HI (2.1)  \nwhere N γ is the mean number density of ionizing photons inside the bubble, c is the speed of light, and σ HI is the HI-ionizing cross-section averaged over the ionizing spectrum. At some time t after reionization starts, N γ is given by  \nN γ ( t ) = ∫ t 0 dt ' [ ˙ N em ( t ' ) -˙ N abs ( t ' )] (2.2)  \nwhere ˙ N em is the emission rate of photons inside the bubble by sources and ˙ N abs is the absorption rate. The latter can be written as two terms,  \n˙ N abs = ˙ N ion abs + ˙ N rec abs (2.3)  \nwhere ˙ N ion abs is the absorption rate by neutral gas at the edges of the bubble, and ˙ N rec abs is due to recombinations within the bubble.  \nWe will first consider the limit where the absorption rate is dominated by ionizations of neutral H atoms, such that ˙ N ion abs >> ˙ N rec abs . This limit is a good approximation early in reionization when the ionized bubbles are much smaller than their recombination-limited Stromgren volumes and most ionized regions have not yet overlapped. The absorption rate can then be written as  \n˙ N abs ( t ) = ˙ N ion abs ( t ) ≈ ˙ N em ( t -R/c ) (2.4)  \nwhere we have approximated the ionized region as a spherical bubble with radius R . The quantity t -R/c is the so-called 'retarded time' at the boundary of the ionized region (the ionization front or I-front) relative to the sources near the bubble's center. That is, photons emitted at time t -R/c will reach the bubble's edge and be absorbed at time t .  \nWe next assume that R/c is sufficiently small that both ˙ N em and R can be approximated as constant across this length of time. In this limit, Eq. 2.2 evaluates to  \nN γ = ˙ N em × ( t -[ t -R/c ]) = ˙ N em R c (2.5)  \nPutting this into Eq. 2.1, we get  \nΓ HI = ˙ N em R c cσ HI = ˙ N em Rσ HI (2.6)  \nIn this limit, Γ HI is independent of the speed of light, and we can get the right Γ HI inside the bubble using the RSLA. Indeed, replacing R in Eq. 2.6 with the mean free path to ionizing photons, λ , reveals that we have reproduced the so-called 'local source approximation' with its well-known relationship between Γ HI , λ , and ˙ N em [e.g. Ref. 26].  \nNext, we consider the opposite limit, in which ˙ N ion abs = 0 . In this case, either the bubble has reached its limiting Stromgren volume (in which case ˙ N rec abs ≈ ˙ N em ), or the universe is entirely ionized. Either way, the right-hand side of Eq. 2.2 will be independent of the speed of light, since neither ˙ N em nor ˙ N rec abs depend on it 1 . Thus, N γ is independent of c , and we have  \nΓ HI ∝ c (2.7)  \nThus, using the RSLA in this limit would result in Γ HI being incorrect by the same factor by which c is reduced.  \nThe first approximation - that ˙ N ion abs dominates the absorption rate - is valid near the beginning of reionization, as we mentioned earlier. The reason Γ HI is independent of c in this limit is that the time required for photons to reach neutral gas is proportional to 1 /c , such that N γ inside the bubble(s) is also proportional to this factor. This cancels the factor of c in Eq. 2.1, which gives the correct Γ HI if a reduced c is used, even though N γ is incorrect. In the other limit, the fact that N γ is independent of c means that it will still be correct when a reduced c is used, causing Γ HI to be incorrect. The first limit is expected to apply near the beginning of reionization, and the second near and after its end. The limiting behavior described here is exactly what was found in Ref. [22] (see also Ref. [21]).", '2.2 Fixing the problem': "In this section, we will show that, to a good approximation, the use of the RSLA equivalent to a redshift-dependent re-scaling of ˙ N em . Consider two simulations, one using the full speed of light ˜ c = c = 1 and the other a reduced speed of light ˜ c < 1 . Let these simulations have global ionizing emissivities ˙ N c em and ˙ N ˜ c em , respectively. Our goal is to find a relationship between these four quantities such that these two simulations have identical physical properties including reionization history, net absorption rate, Γ HI , etc. In particular, if they have the same absorption rate (that is, ˙ N c abs = ˙ N ˜ c abs ), then from Eq. 2.2, we have  \nN c γ ( t ) -∫ t 0 dt ' ˙ N c em ( t ' ) = N ˜ c γ ( t ) -∫ t 0 dt ' ˙ N ˜ c em ( t ' ) (2.8)  \nwhere we have used the same speed-of-light notation for N γ that we did for ˙ N em . If we also assume that both simulations have the same mean Γ HI , then from Eq. 2.1 we have  \nN c γ cσ HI = N ˜ c γ ˜ cσ HI → N c γ = N ˜ c γ ˜ c c (2.9)  \nProvided Eq. 2.9 holds at all times, we can differentiate both sides with respect to time to get  \n˙ N c γ = ˙ N ˜ c γ ˜ c c (2.10)  \nDifferentiating both sides of Eq. 2.8 with respect to time and substituting for N c γ using Eq. 2.10 yields the main result of this paper,  \n˙ N c em ( t ) = ˙ N ˜ c em ( t ) -˙ N ˜ c γ ( t ) ( 1 -˜ c c ) (2.11)  \nWe have explicitly written the time dependence to emphasize that Eq. 2.11 holds at all times.  \nEq. 2.11 expresses the emissivity in the simulation with the full speed of light in terms of c , ˜ c , N ˜ c em , and ˙ N ˜ c γ . The first two of these are just numbers, and the other two are known after the simulation with the reduced speed of light has been run. This allows us to calculate ˙ N c em without running any simulations with the full speed of light . Again, ˙ N c em is the emissivity would produce (approximately) the same ionization history, Γ HI , etc. as the reduced speed of light run, had the full speed of light been used. Of course, since we had to assume that this relationship exists in order to derive Eq. 2.11, we have not yet proven that this statement is true. We will do so using RT simulations in §4-5.  \nThe usefulness of Eq. 2.11 is found in situations where ˙ N em is being treated as a 'free parameter', which can be 'corrected' after the simulation has been run. Indeed, in works such as Refs. [3, 6, 16], ˙ N em is a 'free function' of redshift that is set by hand to match some set of observables - in those works, the mean transmission of Ly α forest at z ≤ 6 . In many semi-numerical frameworks [e.g. 21cmFAST, Refs. 29, 30], ˙ N em is uniquely determined by a set of free parameters that are being constrained using a set of observations. In this case, a re-scaling of ˙ N em corresponds to changing the mapping between these parameters and ˙ N em in a way that can be quantified. We emphasize the Eq. 2.11 does not constitute a correction to the RSLA simulation results themselves, but rather simply a re-scaling of the global emissivity after the simulation has been run. As such, it is unfortunately not directly applicable in situations where ˙ N em is a self-consistent prediction of e.g. an underlying galaxy model [5, 7, 19]. In this situation, one would require a method to correct the simulation quantities themselves, which is beyond the scope of this work.", '3 Numerical Methods': "Our simulations are run with the adaptive ray-tracing RT code FlexRT, first described in [16, 31] and tested in Ref. [32]. Here, we will briefly discuss the relevant features of the code, referring the reader to these works for details. FlexRT solves the RT equation in post-processing on a time-series of cosmological density fields. We model the ionizing sources in a simple manner by extracting halos from an N-body simulation with the same large-scale structure as the density field and assigning them ionizing emissivities (see below). The RT equation itself is solved using an adaptive ray tracing method similar to the one described in Ref. [11] and employed in the code of Ref. [12]. Sub-resolved ionization fronts are tracked using the 'moving-screen' approximation, and the temperature behind them ( T reion ) is estimated using the flux-based method described in Ref. [17] 2 . FlexRT solves for the opacity of the ionized IGM to ionizing photons using a novel sub-grid model based on high-resolution hydro/RT simulations of IGM gas dynamics similar to those described in Refs. [18, 27]. Recent improvements to this model, and details about forward-modeled observables are detailed in Ref. [23].  \nIn §4, we will use the 'Reference' setup in Ref. [23], which assumes a fully dynamically relaxed IGM and mono-chromatic RT. The physical properties of this model are summarized in their Fig. 2 and associated text. Throughout in §4, the only parameters we will change relative to their Reference model are ˙ N em and ˜ c . We will generalize and test Eq. 2.11 in the context of multi-frequency RT simulations in §5. Our simulations have a box size of L = 200 h -1 Mpc with N RT = 200 3 RT cells. Ly α forest properties are calculated (following §3.7 of  \nTable 1 . Summary of different kinds of simulations referred to in this work, and the names used for each (left-most column). The second column from left gives the emissivity used in that simulation, and the third column the speed of light. The right-most column briefly describes each type of simulation.  \nRef. [23]) on a high-resolution ( N = 2048 3 ) hydrodynamics simulation with the same largescale initial conditions used for the N-body simulation from which we get the halos. In this work, we will consider reduced speeds of light of ˜ c = 0 . 05 , 0 . 1 , 0 . 2 , 0 . 3 , and 0 . 5 , alongside simulations that use the full speed of light ( ˜ c = c = 1 ).", '4.1 Strategy & Terminology': 'The claim of Eq. 2.11 can be summarized as follows: Given an RT simulation run with a reduced speed of light ˜ c and emissivity ˙ N ˜ c em ( z ) that has a mean photon number density N ˜ c γ ( z ) , a simulation using the full speed of light and emissivity ˙ N c em ( z ) should have the same reionization history, ionization morphology, IGM transmission properties, etc. as the original simulation. We can verify whether this claim is true by comparing two types of simulations. We first run simulations using the RSLA, with ˜ c < 1 . We refer to these simply as reducedc runs, with emissivity ˙ N ˜ c em . Next, we apply the re-scaling in Eq. 2.11 to ˙ N ˜ c em to obtain ˙ N c em , and run another simulation with the full speed of light using ˙ N c em . We refer to these as fullc , re-scaled ˙ N em runs. Eq. 2.11 is accurate to the extent that the properties of these two types of simulations agree with each other. Lastly, for purposes of comparison, we will also include a third type of simulation, in which we use the full speed of light but do not apply any re-scaling to ˙ N ˜ c em . We refer to these as fullc , un-scaled ˙ N em runs. We summarize the properties of these three types of simulations in Table 1.', '4.2 Ionization history and morphology': "We begin by studying how well the reionization history and morphology (shapes and sizes) of ionized regions match in the reducedc and fullc , re-scaled ˙ N em runs. For this section, we scale up ˙ N em from the reference model of Ref. [23] until reionization ends at z ≈ 6 when c = 1 . This is our fullc , un-scaled ˙ N em run. Next, we use the same ˙ N em to run a series of reducedc runs with ˜ c = 0 . 05 , 0 . 1 , 0 . 2 , 0 . 3 , and 0 . 5 . Lastly, we evaluate Eq. 2.11 for  \nFigure 1 . Summary of fullc , un-scaled ˙ N em and reducedc simulations used to study how well Eq. 2.11 works with respect to the reionization history and morphology. Left: the reionization history for each value of ˜ c , using a scaled-up version of ˙ N em from the reference model from Ref. [23], which ends reionization at z = 6 for ˜ c = c = 1 (black solid curve, the fullc , un-scaled ˙ N em run). Reducing ˜ c pushes the end of reionization later, such that it has not finished by z = 4 . 8 for ˜ c = 0 . 05 . Middle: number density of ionizing photons, N γ . The thick lines denote the actual values, while the faded lines show ˜ cN γ . Consistent with our analytic argument in §2.1, at the beginning of reionization ˙ N γ is a factor of 1 / ˜ c higher in the reducedc runs than in the fullc , un-scaled ˙ N em case (Eq. 2.5). Near reionization's end, they approach similar values. Right: Ionizing emissivity used in all these simulations ( ˙ N ˜ c em ).  \n<!-- image -->  \neach of these and then run a set of fullc , re-scaled ˙ N em simulations. Again, the degree to which the reducedc and corresponding fullc , re-scaled ˙ N em runs agree with each other determines how accurate is Eq. 2.11.", '4.2.1 Ionization history': "In Figure 1, we show the volume-averaged ionized fraction x V ion (left), the ionizing photon number density ˙ N γ (middle), and the ionizing emissivity ˙ N em (right) for our fullc , unscaled ˙ N em run (black solid curve) and our reducedc runs (colored curves, see legend). In the left panel, we see that the reionization histories are initially similar, but start to diverge as reionization enters its later stages. The end of reionization is delayed by ∆ z ≈ 1 . 5 in the ˜ c = 0 . 05 simulation relative to the ˜ c = 1 case. The bold curves in the middle panel show ˙ N γ for each simulation, while the faded curves show ˜ cN γ . Early in reionization, ˙ N γ scales linearly with 1 / ˜ c , consistent with the early reionization limit described by Eq. 2.5. However, as reionization ends, ˙ N γ for all simulations approaches the same value, reflecting the behavior expected near reionization's end in §2.1. The right-most panel shows ˙ N ˜ c em , which is the same for all the simulations shown here.  \nIn the left panels Figure 2, we compare each of our reducedc runs to their fullc , re-scaled ˙ N em counterparts. We show the former as bold curves and the latter as faded curves, with the color and line style identifying the value of ˜ c . The upper left panel shows x V ion , zoomed in on the second half of the reionization history for clarity. We find excellent agreement between the two even for our lowest value of ˜ c = 0 . 05 , such that the different sets of curves are difficult to see on the plot. The bottom left panel shows the linear differences of x V ion in the reducedc and fullc , re-scaled ˙ N em runs over the entire reionization history. For all but the ˜ c = 0 . 05 case, the difference is always ≲ 0 . 01 , and for ˜ c ≥ 0 . 3 it is < 0 . 005 .  \nIn the upper right panel, we show ˙ N c em calculated using Eq. 2.11, and the bottom right panel shows the ratio with ˙ N ˜ c em (the black solid curve). This ratio differs from unity the most at the end of reionization, dropping to ≈ 0 . 5 ( 0 . 8 ) by z = 5 for ˜ c = 0 . 05 ( 0 . 5 ). The fact  \nFigure 2 . Comparison of the reionization history and emissivity in our reducedc and fullc , rescaled ˙ N em simulations. Upper Left: the same reionization histories (again), but zoomed in on the last half of reionization and including the fullc , re-scaled ˙ N em results for each value of ˜ c (faded curves). The bold and faded curves agree very well for all values of ˜ c , such that they are difficult to tell apart on the plot. Lower Left: the linear differences between x V ion for the reducedc and corresponding fullc , re-scaled ˙ N em runs. This is less than 0 . 01 for all except ˜ c = 0 . 05 . Upper Right: emissivities for each of the fullc , re-scaled ˙ N em runs ( ˙ N c em from Eq. 2.11, colored curves) compared to ˙ N ˜ c em (black solid curve). Bottom Right: ratio of ˙ N c em and ˙ N ˜ c em for each ˜ c . The ratio is smallest at the end of reionization, dropping to 0 . 5 ( 0 . 8 ) for ˜ c = 0 . 05 ( 0 . 5 ).  \n<!-- image -->  \nthat the reducedc and fullc , re-scaled ˙ N em runs have such different ˙ N em but similar reionization histories validates Eq. 2.11. As we will see, the small differences in reionization history arise because Eq. 2.11 applies to globally averaged quantities, and does not account for spatial fluctuations in the distribution of photons in the IGM.  \nThese results show that, at least with respect to the global ionization history, the claim of Eq. 2.11 is reasonably accurate. Namely, that using the RSLA is equivalent to a redshiftdependent re-scaling of ˙ N em . We emphasize again that Eq. 2.11 does not prescribe a way to directly correct the reducedc simulation results for the effects of the RSLA. Put another way, we do not prescribe a way to convert the reducedc results into those of the fullc , un-scaled ˙ N em results. As such, our method is not directly applicable when ˙ N em is a prediction of e.g. some underlying galaxy model.", '4.2.2 Ionization morphology': 'Next, we apply the same analysis to the morphology of ionized and neutral regions, and statistics that depend on it. Figure 3 shows maps of the ionization field when the universe is 80% ionized for reducedc (top row) and fullc , un-scaled ˙ N em runs (middle row). Dark regions are neutral and white ones are ionized. The bottom row shows the difference maps for the top and middle rows. Red (blue) denotes regions that are more neutral (ionized) in the reducedc simulations.  \nFrom right to left (decreasing ˜ c ) in the top row, the neutral regions (islands) grow slightly more extended and porous. Equivalently, the smallest ionized bubbles grow larger, and the largest ones smaller, when ˜ c decreases. We see no such trend, however, in the middle row. From this, we deduce that the morphological differences in the top row are not due to the reducedc runs having different reionization histories, since such differences would appear  \nFigure 3 . Visualization of the ionization field at 80% ionized in reducedc (top row) and fullc , re-scaled ˙ N em runs (middle row) for ˜ c = 0 . 05 , 0 . 1 , 0 . 2 , and 0 . 3 (left to right). We show the difference maps for each value of ˜ c in the bottom row. We see in the top row that as ˜ c decreases (right to left), the neutral regions (islands) become more extended and porous at fixed ionized fraction. We do not see this in the middle row, indicating that it is not due to differences in the reionization history between reducedc runs. Instead, it is because the RSLA induces a longer delay between emission and absorption for photons emitted in large vs. small ionized bubbles. The re-scaling in Eq. 2.11 corrects for the delay only on average, and does not account for position dependence. The result is that the morphologies in the reducedc and fullc , re-scaled ˙ N em runs diverge as ˜ c decreases.  \n<!-- image -->  \nfor the fullc , re-scaled ˙ N em also. They are instead a position-dependent artifact of the RSLA. Photons emitted in large ionized regions experience a longer delay between emission and absorption than do those emitted in small bubbles. Eq. 2.11 accounts for this time-delay only on average, and thus misses the differences between photons emitted in different regions. As a result, the largest (smallest) bubbles are too small (large) in the reducedc runs. We can most clearly see this effect in the difference maps, where it is clearly visible even for ˜ c = 0 . 3 . In this work, we do not attempt to correct for this position-dependent effect, since doing so would be complicated by the non-locality of the ionizing radiation field. As we will see, this effect will set the minimum value of ˜ c for which Eq. 2.11 can be used.  \nFigure 4 . Power spectrum of the ionization field and 21 cm signal for simulations with ˜ c = 0 . 05 , 0 . 1 , 0 . 2 , and 0 . 3 (columns from left to right). Top Row: power spectrum of the ionization field, ∆ ion , vs. wavenumber k at an ionized fraction of 80% . The bold colored curves show results for the reducedc simulations, and the faded colored curves the fullc re-scaled runs. The faded gray curves shows (at fixed redshift) the fullc , un-scaled ˙ N em runs. 2nd Row: the same, but for the 21 cm power spectrum, ∆ 21 . The fullc , un-scaled ˙ N em results differ dramatically from the other two, reflecting its lower neutral fraction at fixed redshift. However, the reducedc and fullc , re-scaled ˙ N em runs differ in their shape: the former has significantly less power on large scales (small k ) and slightly more on small-scales (large k ). 3rd Row: ∆ 21 vs redshift at k = 0 . 2 and 0 . 1 h Mpc -1 (top and bottom set of curves, respectively). The bold colored curves are noticeably below the faded ones, reflecting the deficit in large-scale power in the reducedc runs. 4th & 5th Rows: ratios of the reducedc and fullc re-scaled results in the 2nd and 3rd rows, respectively. The shaded regions denote ± 20% from unity. In the bottom row, we show the ratio for k = 0 . 1 h Mpc -1 . Decreasing ˜ c increases differences at low k , so that for ˜ c = 0 . 1 the difference is ≈ 20% at all redshifts at k = 0 . 1 h Mpc -1 . At ˜ c ≥ 0 . 2 , we find ≲ 15% differences at k = 0 . 1 h Mpc -1 , which get smaller at larger k .  \n<!-- image -->  \nFigure 4 quantitatively compares the large-scale fluctuations in the ionization field in all three types of simulations. In the top row, we show the dimensionless ionization power spectrum ∆ ion ( k ) at 80% ionized for ˜ c = 0 . 05 , 0 . 1 , 0 . 2 , and 0 . 3 (left to right). The bold colored lines, faded colored lines, and faded gray lines indicate results for reducedc , fullc , re-scaled ˙ N em , and fullc , un-scaled ˙ N em runs, respectively. The second row is the same thing, but instead showing the 21 cm power spectrum 3 , calculated using the tools21cm  \npackage of Ref. [33].  \nThe power in the fullc , un-scaled ˙ N em run is much lower than the others - increasingly so at lower ˜ c . This reflects the fact that the reducedc runs are more delayed relative to the fullc , un-scaled ˙ N em case the smaller ˜ c is. So, at a fixed neutral fraction in the former, the neutral fraction in the latter decreases for smaller ˜ c . The fullc , re-scaled ˙ N em sims are much closer to the reducedc results, but differ noticeably in their shape. The latter have significantly less power at large scales (small k ), and slightly more power at small scales (large k ). This is a direct result of the morphological errors shown in Figure 3, with the reduced sizes of the largest bubbles driving the decrease in smallk power. In the third row, we show ∆ 21 vs. redshift at fixed k = 0 . 2 and 0 . 1 h Mpc -1 (top and bottom set of curves, respectively), key targets for experiments such as HERA [34-36]. The suppression in power is consistent across the entire reionization history, indicating that even early in reionization the RSLA can cause significant errors in ionization morphology.  \nThe fourth row shows the ratio of the bold and faded colored curves in the second row, with the shaded region denoting ± 20% . For ˜ c = 0 . 3 , the differences between the reducedc and fullc re-scaled runs never exceeds 20% even at the largest scales at 80% ionized. For lower ˜ c , the drop in power at small k increases, reaching a factor of ≈ 2 at the box scale for ˜ c = 0 . 05 . For ˜ c = 0 . 2 the deviation is still only ∼ 10% at k = 0 . 1 h Mpc -1 , but at ˜ c = 0 . 1 this increases to 20% , and to 30% for ˜ c = 0 . 05 . The bottom row shows the ratio of the k = 0 . 1 h Mpc -1 curves in the third row. Here we see that for ˜ c = 0 . 2 , the differences are never larger than 15% at any redshift. For ˜ c = 0 . 1 , the difference exceeds 20% at z ∼ 9 , and for ˜ c = 0 . 05 , it reaches nearly a factor of 2 .  \nWe see that the utility of Eq. 2.11 is limited by the fact that it does not account for the spatial fluctuations in the delay between emission and ionization of ionizing photons in the reducedc simulations. To summarize, ∆ 21 at k ≤ 0 . 1 h Mpc -1 can deviate by 20% between the reducedc and fullc , re-scaled ˙ N em runs for ˜ c = 0 . 1 , and by up to 40% for ˜ c = 0 . 05 . These differences drop to ≤ 15% for ˜ c = 0 . 2 and ≤ 10% for ˜ c = 0 . 3 . As such, we caution against applying our method for ˜ c < 0 . 2 in contexts where the large-scale reionization morphology is important.', '4.3 QSO observations at z ≤ 6': 'In this section, we apply the same kind of test to observables derived from high-redshift quasar absorption spectra at z ≤ 6 . These include the Ly α forest, which has recently yielded detailed insights into the tail end of reionization at 5 < z < 6 [3, 6, 24, 30, 37], and the ionizing photon mean free path [38-41]. These observables are sensitive not only to ionization morphology, but also to the large-scale spatial fluctuations in the ionizing background Γ HI and the IGM temperature, which are affected by the RSLA.', '4.3.1 Global Observables': "The evolution of the mean transmission of the Ly α forest is a tight boundary condition on the end of reionization. Several studies [3, 14, 23, 42] have 'calibrated' their reionization simulations to agree by construction with this observable, then asked whether other observables predicted by the simulation agree well with observations. Here, we will follow this procedure and assess how applicable Eq. 2.11 is in the context of forest-anchored studies. We first run a set of reducedc simulations with ˙ N ˜ c em ( z ) calibrated so that the mean transmission of the Ly α forest matches the recent observations of Ref. [24]. Note that unlike in the previous section, this procedure yields a different ˙ N ˜ c em for each value of ˜ c , since we are requiring all  \nFigure 5 . Tests with simulations calibrated to match the Ly α forest at z = 6 . Upper Left: ⟨ F Ly α ⟩ vs. redshift for all three types of simulations described in Table 1 for ˜ c = 0 . 1 , 0 . 2 , and 0 . 3 . The formatting of the curves is the same used in Figure 4 (see text and legend). By construction, all the reducedc results agree with each other and the measurements from Ref. [24]. The fullc , re-scaled ˙ N em , ˜ c = 0 . 3 run is indistinguishable from these. The ˜ c = 0 . 2 run is slightly above the others, but still in reasonable agreement with the measurements. However, the ˜ c = 0 . 1 run diverges significantly from the measurements. All the fullc , un-scaled ˙ N em have too much transmission, since they end reionization earlier than the reducedc runs. Bottom Left: ratio of ⟨ F Ly α ⟩ in the fullc , re-scaled ˙ N em and reducedc runs, with the shaded region denoting ± 20% . We find ≲ 5% agreement for ˜ c = 0 . 3 and ≲ 20% agreement for ˜ c = 0 . 2 , but for ˜ c = 0 . 1 , the simulations disagree by a factor of 1 . 3 -2 . Upper Right: ˙ N em for the reducedc (bold curves) and fullc , re-scaled ˙ N em (faded curves) simulations. As previously mentioned, the reducedc runs have different ˙ N em since they are all separately calibrated to match the forest. The fullc , re-scaled ˙ N em runs have similar ˙ N em , which is expected since they use the full speed of light and are supposed to match the same observable, ⟨ F Ly α ⟩ . Bottom Right: ratio of ˙ N em in the fullc , re-scaled ˙ N em and reducedc runs. They are similar early in reionization but diverge by up to a factor of 3 by its end.  \n<!-- image -->  \nreducedc runs to match the same observable. For each ˜ c , we apply Eq. 2.11 to ˙ N ˜ c em and run fullc , re-scaled ˙ N em simulations using c = 1 . In some plots, we will also include results from fullc , un-scaled ˙ N em runs. We will judge the accuracy of Eq. 2.11 in the same way as in the previous section - by how well the reducedc and fullc , re-scaled ˙ N em runs agree with each other.  \nWe first compare the mean transmission of the Ly α forest, ⟨ F Ly α ⟩ , with measurements from Ref. [24] in the top left panel of Figure 5. The bold colored lines show reducedc runs for ˜ c = 0 . 1 , 0 . 2 , and 0 . 3 . Faded colored lines (with matching line styles) show the fullc , re-scaled ˙ N em results, and faded gray lines (again, with matching line styles) show the fullc , un-scaled ˙ N em results. The bottom panel shows the ratios of the fullc , rescaled ˙ N em and reducedc results (with the shaded line again denoting ± 20% from unity). All the reducedc results agree well with the measurements and each other by construction. The fullc , re-scaled ˙ N em run with ˜ c = 0 . 3 is indistinguishable from its reducedc counterpart. For ˜ c = 0 . 2 , ⟨ F Ly α ⟩ for the fullc , re-scaled ˙ N em run slightly above the measurements, but is still within at most 25% of its reducedc counterpart. However, for  \nFigure 6 . The same as Figure 5, but showing λ 912 mfp (left) and T 0 (right). To avoid cluttering the plot, we omit the fullc , un-scaled ˙ N em results. For both ˜ c = 0 . 3 and 0 . 2 , we see agreement to within a few percent in both quantities between the reducedc and fullc , re-scaled ˙ N em runs. However, for ˜ c = 0 . 1 , two diverge by up to ≈ 35% for λ 912 mfp and ≈ 6% for T 0 . These results further confirm our earlier observations that our method should only be applied for ˜ c ≥ 0 . 2 .  \n<!-- image -->  \n˜ c = 0 . 1 , the fullc , re-scaled ˙ N em run diverges by as much as a factor of 2 . All three of the fullc , un-scaled ˙ N em lie well above the measurements, since their reionization histories end too early.  \nIn the upper right panel, we show the emissivities for the reducedc and fullc , rescaled ˙ N em runs ( ˙ N ˜ c em and ˙ N c em , respectively). The bottom panel shows their ratio. We see that ˙ N ˜ c em is quite different for different values of ˜ c , since all runs were calibrated to agree with the same observable. On the other hand, ˙ N c em is quite similar between the different values of ˜ c . This is indeed to be expected if Eq. 2.11 is accurate, since each ˙ N c em is supposed to match the same observable, ⟨ F Ly α ⟩ , with the full speed of light. In the lower right panel, we see that the ˙ N c em and ˙ N ˜ c em are similar early in reionization, but differ by a factor of 2 -3 by its end.  \nIn Figure 6, we make the same comparison for the lyman limit mean free path 4 ( λ mfp 912 , MFP, left) and the IGM temperature at mean density T 0 (right). To avoid cluttering the plot, we omit results from the fullc , un-scaled ˙ N em runs. We also include several sets of observations for comparison (see references in caption). For ˜ c = 0 . 2 and 0 . 3 , the reducedc and fullc , re-scaled ˙ N em agree to within a few percent in both λ 912 mfp and T 0 , as the bottom panels show. For ˜ c = 0 . 1 , the disagreement is much larger - up to 35% for λ 912 mfp and ≈ 6% for T 0 . Thus, the trend here is qualitatively similar to that observed in Figures 4 and 5 namely, that Eq. 2.11 gives a reasonable match between reducedc and fullc , re-scaled ˙ N em simulations for ˜ c = 0 . 2 , but not for ˜ c = 0 . 1 . This further confirms that our approach can only be reliably used down to ˜ c = 0 . 2 .", '4.3.2 Large-scale fluctuations': 'We saw in §4.2.2 and §4.3.1 that the applicability of Eq. 2.11 is likely limited to ˜ c ≳ 0 . 2 by spatially-varying time-delay effects caused by the RSLA on large scales. In this section,  \nFigure 7 . Distribution of Ly α forest effective optical depths over 50 h -1 Mpc segments of the forest, P ( < τ 50 eff ) , at z = 5 . 4 , 5 . 6 , 5 . 8 , and 6 . 0 . The line styles and colors have the same meanings as in Figure 6. At z = 5 . 8 and 6 , when the neutral fraction is > 10% , the reducedc results differ from each other at the few percent level. This is unsurprising since, as we showed earlier, they differ in large-scale ionization morphology. The ˜ c = 0 . 2 and 0 . 3 fullc , re-scaled ˙ N em agree with their reducedc counterparts to within a few percent, but this deviation grows noticeably larger for ˜ c = 0 . 1 , especially at z = 5 . 4 and 5 . 6 . This suggests that the effect of the RSLA on fluctuations in Γ HI and T may be affecting these differences when the neutral fraction is < 10% .  \n<!-- image -->  \nwe study this effect in more detail in the context of the Ly α forest. At the end of reionization, large-scale fluctuations in forest properties are set by three quantities that the RSLA affects: the ionization morphology [3], large-scale fluctuations in Γ HI [44], and fluctuations in temperature [45]. We will focus on the latter two in this section.  \nIn Figure 7, we show the cumulative distribution function (CDF) of τ 50 eff , P ( < τ 50 eff ) , the effective optical depth over 50 h -1 Mpc segments of the forest, at z = 5 . 4 , 5 . 6 , and 5 . 8 , and 6 . 0 . This statistic has been used in a number of works to argue that reionization must have ended later than z = 6 [3, 6, 30, 37, 42, 46]. The format of the curves is the same as that in Figure 6. All the reducedc runs, and the fullc , re-scaled ˙ N em runs with ˜ c = 0 . 2 and 0 . 3 (faded blue-dotted and green dot-dashed curves) mutually agree to within a few percent at all redshifts. We do see few-percent level differences between different reducedc runs at z = 6 and 5 . 8 , which get smaller at z = 5 . 6 and 5 . 4 . This is unsurprising, since we have already shown that the reducedc runs differ from each other in large-scale ionization morphology (Figures 3-4), which also affects the τ 50 eff CDF when the IGM is still partially neutral. We  \n<!-- image -->  \n<!-- image -->  \nFigure 8 . Visualization of the effect of the RSLA on large-scale fluctuations in Γ HI (in units of 10 -12 s -1 ). We show maps of the reducedc , fullc , re-scaled ˙ N em , and fullc , un-scaled ˙ N em runs in the upper left, upper right, and lower left panels, respectively, for ˜ c = 0 . 2 . The fullc , un-scaled ˙ N em run has a higher mean Γ HI , fewer neutral islands (white regions), and weaker spatial fluctuations in Γ HI , since reionization ends earlier in that run than in the others. The reducedc and fullc , re-scaled ˙ N em runs are in reasonably good visual agreement. However, the former has noticeably larger Γ HI fluctuations, with the bright (faint) regions being brighter (fainter) than in the fullc , re-scaled ˙ N em run. In the lower right, we plot log-space PDFs of Γ HI normalized by its mean. We show results for all three values of ˜ c , with the same line styles and colors used in Figure 5. We have offset the ˜ c = 0 . 1 ( 0 . 3 ) PDFs by 1 . 2 dex to the left (right) for clarity. In all three cases, the PDFs are much narrower in the fullc , un-scaled ˙ N em runs than in the others. As ˜ c decreases, the PDFs in the reducedc become grow wider than their fullc , re-scaled ˙ N em counterparts, owing to time-delay effects (see text).  \n<!-- image -->  \nalso see that the fullc , re-scaled ˙ N em ˜ c = 0 . 1 run differs significantly from the others, especially at z = 5 . 6 and 5 . 4 , when the neutral fraction is < 10% . This suggests that the effect of the RSLA on large-scale fluctuations in Γ HI and/or T may also play a significant role in driving these differences.  \nWe visualize the effect of the RSLA on large-scale fluctuations in Γ HI in Figure 8. We show maps of Γ HI (in units of 10 -12 s -1 ) at z = 5 . 65 , when the universe is ≈ 10% neutral in the reducedc and fullc , re-scaled ˙ N em runs. The top left, top right, and lower left  \n<!-- image -->  \npanels show maps for the reducedc , fullc , re-scaled ˙ N em , and fullc , un-scaled ˙ N em runs, respectively, for ˜ c = 0 . 2 . We see that the fullc , un-scaled ˙ N em run has a much higher mean Γ HI and weaker Γ HI fluctuations than the other two, owing to its earlier end to reionization. Compared to this run, the reducedc and fullc , re-scaled ˙ N em runs are in relatively good agreement. However, the spatial fluctuations in Γ HI are noticeably stronger in the reducedc run. Once again, this is due to time-delay effects caused by the RSLA. Because photons are delayed longest in reaching the edges of the largest ionized regions, the differences in N γ between them and the smaller bubbles are amplified. In addition, photons far from sources were emitted at a retarded time t ret = R/ ˜ c , where R is the ionized bubble size. Again, these effects are accounted for only on average by the re-scaling of Eq. 2.11, which results in larger spatial Γ HI fluctuations in the reducedc runs than in their fullc , re-scaled ˙ N em counterparts.  \nIn the lower right panel, we show the PDFs of Γ HI , normalized by its mean (in log space). Here, we show results for all three values of ˜ c , using the same color and line style convention adopted in Figure 5. To make the plot readable, we have offset the ˜ c = 0 . 1 ( 0 . 3 ) PDFs to the left (right) by 1 . 2 dex. In all cases, the gray curves are much narrower than the others, reflecting the reduced Γ HI fluctuations in the fullc , un-scaled ˙ N em runs. As ˜ c decreases, the PDFs of the reducedc runs grow noticeably wider than those of their fullc , re-scaled ˙ N em counterparts, with the difference becoming very significant for ˜ c = 0 . 1 . These differences contribute to those in P ( < τ 50 eff ) seen in Figure 7.  \nFigure 9 shows the same thing as Figure 8, but for IGM temperature. In the lower right panel, we offset the ˜ c = 0 . 1 and 0 . 3 results by ± 30000 K for clarity. The voids, which have ionized most recently, are hottest in the fullc , un-scaled ˙ N em run, which ends reionization the earliest and most rapidly. The reducedc run displays the smallest such T enhancements. The fastest I-fronts around the largest ionized regions, which produce the largest T reion [17, 47], grow move more slowly than they should in the reducedc runs for aforementioned reasons. We see in the lower right panel that the difference between the reducedc and fullc , re-scaled ˙ N em in the highT tail of the PDF becomes significant for ˜ c = 0 . 1 , explaining the corresponding differences in T 0 in Figure 6.  \nWe have repeated the tests in this section, in whole or in part, for several variations of the properties of the ionized sources and the IGM assumed in this work. We have first varied the clustering of ionizing sources assumed in our fiducial scenario, which assumes that the ionizing emissivity of halos scales with their UV luminosity. We have tested models wherein the bulk of the ionizing photons are emitted by the only the least massive (least clustered) and most massive (most clustered) halos. We have also applied our method to multi-frequency RT simulations, for which we show results in the next section. Lastly, we tested models that assume more sub-grid recombinations than the Reference model of Ref. [23], and we have varied the reionization history. We find that our results remain essentially the same for all of these variations of our fiducial scenario - namely, that Eq. 2.11 is accurate (to 20% or better) for ˜ c = 0 . 2 , but that using ˜ c = 0 . 1 leads to much larger errors.  \nThese results further corroborate our conclusions from §4.2. We find that, for ˜ c ≥ 0 . 2 , reducedc simulations calibrated to match the Ly α forest match well the properties of their fullc , re-scaled ˙ N em counterparts. This confirms that Eq. 2.11 can be applied to simulations using the RSLA with ˜ c ≥ 0 . 2 . However, for ˜ c = 0 . 1 , we see significant differences between the reducedc and fullc , re-scaled ˙ N em runs. These can largely be attributed to position-dependent time-delay effects that arise from using a reduced speed of light, for which Eq. 2.11 does not account. As such, for applications in which the morphology of  \nFigure 9 . The same as Figure 8, but for the IGM gas temperature. The main difference between the three maps is the strength of the temperature enhancements near neutral islands, where gas has been recently reionized and is hottest. The fullc , un-scaled ˙ N em run has the hottest temperatures in the voids, since reionization ends earliest and most rapidly. The reducedc run, by contrast, has the weakest T fluctuations. This is because the slowed growth of the largest ionized regions (see §4.2.2) reduces the speed of the fastest I-fronts, which lowers T reion around the largest bubbles. In the lower right panel, we see a significant difference between reducedc and fullc , re-scaled ˙ N em runs for ˜ c = 0 . 1 at the highT end of the PDF, which explains the differences in T 0 seen in Figure 6.  \n<!-- image -->  \nreionization and/or quasar-based observables at z ≤ 6 are important, we caution against applying our re-scaling method for ˜ c smaller than 0 . 2 .', '5 Generalization to multi-frequency simulations': "In the previous section, we tested Eq. 2.11 in mono-chromatic simulations. Here, we will generalize Eq. 2.11 to apply to multi-frequency simulations and demonstrate that the method works just as well in this case, but with one key additional caveat. For our purposes, the main effect of including multi-frequency RT is to change the temperature and ionizing background in the ionized IGM [23], but not the morphology of ionized regions. As such, we will only show our Ly α forest-focused tests (§4.3) in this section.  \nFigure 10 . Same as Figure 5, but showing results for our tests with multi-frequency RT, using Eq. 5.1 to re-scale the emissivity. We find very similar results in all panels to those in Figure 5, showing that our method works nearly as well for multi-frequency RT simulations as for single-frequency ones. See text for details.  \n<!-- image -->  \nA straightforward multi-frequency generalization of Eq. 2.11 is to simply apply the arguments of §2.2 in each frequency bin separately. Then we have  \n˙ N c em ( t, ν ) = ˙ N ˜ c em ( t, ν ) -˙ N ˜ c γ ( t, ν ) ( 1 -˜ c c ) (5.1)  \nA subtlety of Eq. 5.1 is that in general, ˙ N ˜ c em ( t, ν ) and ˙ N ˜ c γ ( t, ν ) do not share the same shape in frequency space. Indeed, this is the case during reionization because the absorption crosssection of HI is frequency-dependent, which results in hardening of the radiation field by the IGM. As such ˙ N c em ( t, ν ) will have a softer ionizing spectrum than ˙ N ˜ c em ( t, ν ) , since ˙ N ˜ c γ ( t, ν ) has a harder spectrum. This complicates somewhat the physical interpretation of the re-scaling.  \nIn this section, we use the same setup used for multi-frequency RT simulations described in Ref. [23], assuming that the intrinsic spectrum of ionizing sources is a power law of the form J ν ∝ ν -α with α = 1 . 5 . We use 5 frequency bins, with central frequencies chosen such that (initially) all bins contain an equal fraction of the emitted ionizing photons. The condition for choosing the frequency bins is that the average HI-ionizing cross-section, σ HI , should be the same as that of an α = 1 . 5 power law. As in the previous section, we ran reducedc simulations with ˜ c = 0 . 1 , 0 . 2 , and 0 . 3 , and we apply Eq. 5.1 to run their fullc , re-scaled ˙ N em counterparts. In Figure 10, we show ⟨ F Ly α ⟩ in the same format as Figure 5 for our multifrequency tests, except that we omit the fullc , un-scaled ˙ N em results. We find results very similar to those in Figure 5. The mean flux in the reducedc and fullc , re-scaled ˙ N em runs remains within 10% for ˜ c = 0 . 3 , and within 30% or better for ˜ c = 0 . 2 . For ˜ c = 0 . 1 , we again find differences as large as a factor of 2 in ⟨ F Ly α ⟩ . In the right column, we see that ˙ N ˜ c em and ˙ N c em , and the relationship between them, are very similar in the multi-frequency runs as in the mono-chromatic case. This shows that our method works just as well when applied to multi-frequency simulations with respect to the mean Ly α forest transmission.  \nFigure 11 . Same as Figure 7, but showing P ( < τ 50 eff ) for our multi-frequency RT tests. We find results very similar to those of the mono-chromatic runs, showing that our re-scaling procedure works well in the multi-frequency case. See text for details.  \n<!-- image -->  \nIn Figure 11, we show P ( < τ 50 eff ) , in the same format as Figure 7, for our multi-frequency tests. Again, we see results very similar to those in Figure 7. The differences between the bold and faded red curves are slightly larger than they are in Figure 7, suggesting that spatial fluctuations in the spectrum of the ionizing background may be slightly worsening the timedelay effects discussed in §4.3.2. However, this effect is not significant enough to change the level of accuracy obtained with ˜ c = 0 . 2 and 0 . 3 , further validating that for these values of ˜ c , our method can safely be applied to simulations with multi-frequency RT.  \nAs mentioned, applying Eq. 5.1 in each frequency bin results in ˙ N c em ( ν ) having a different spectral shape than ˙ N ˜ c em ( ν ) , in a way that depends on redshift. We illustrate this in Figure 12, where we quantify the spectral shape of ˙ N em before and after the re-scaling. The black solid curve in the left panel shows the fraction of the total emissivity in each bin in our multifrequency reducedc run with ˜ c = 0 . 2 (recall that all bins share the same fraction of the photons in this case). In the right panel, the black dotted curve shows 1 N freq ˙ N c em , the fraction that would be in each bin if ˙ N c em and ˙ N ˜ c em shared the same frequency dependence. The colored dot-dashed lines show ˙ N c em ( ν ) for each bin, with redder (bluer) colors denoting lower (higher) photon energies (see annotations). The re-scaling removes significantly more photons at higher energies, causing ˙ N c em ( ν ) to get softer with decreasing redshift.  \nIn the right panel, we calculate the 'effective' power law index, α eff , for each case vs.  \nFigure 12 . Illustration of how our re-scaling procedure changes the intrinsic spectrum of the sources when applied to multi-frequency RT simulations. Left : Ionizing emissivity in each bin in the reducedc run (black solid) compared to the spectrum in the fullc , re-scaled ˙ N em run. The black dotted line shows the average fraction of photons in each bin in the latter, while the colored dot-dashed lines show the actual fraction in each of the 5 frequency bins after re-scaling. Redder (bluer) lines indicate lower (higher) photon energies, as the annotations indicate. The re-scaling removes more high-energy photons from ˙ N c em due to IGM filtering. Right : effective spectral index α eff (see text) for the multi-frequency reducedc and fullc , re-scaled ˙ N em simulations. The latter has a softer spectrum (larger α eff ) which gets softer with time, reaching 2 . 5 by z = 5 .  \n<!-- image -->  \nredshift. We define α eff such that a power law spectrum with that index would have the same σ HI as the source spectrum in the simulation. For the reducedc run, this is trivially 1 . 5 , but for the fullc , re-scaled ˙ N em run, α eff increases with time, reaching ≈ 2 . 5 by z = 5 . Note that this effect arises from the fact that the IGM absorbs (filters) low-energy photons more readily than high-energy ones, such that ˙ N γ has a harder spectrum than ˙ N em . This effect will thus be less significant in simulations with less filtering. We showed in Appendix B of Ref. [23] that our simulations are likely predict levels of filtering on the high end of expectations, so the effect may be somewhat exaggerated here. Still, it represents a complication to applying Eq. 5.1 to multi-frequency simulations that should be carefully considered and quantified in any studies that use this approach.", '6 Conclusions': "In this work, we have studied the effect of the reduced speed of light approximation (RSLA) on radiative transfer simulations of reionization. We used a simple analytic model to show (1) when and why the RSLA produces inaccurate results, especially near reionization's end, and (2) that using the RSLA is, to a good approximation, equivalent to a redshift-dependent re-scaling of the global ionizing emissivity of reionization's sources (Eq. 2.11). We have run simulations with the reduced speed of light ( reducedc runs) and with the full speed of light after applying the re-scaling prescribed by Eq. 2.11 ( fullc , re-scaled ˙ N em runs). We have assessed how accurate Eq. 2.11 is by comparing these two sets of simulations in a number of different physical properties and observables. Our main findings are summarized below:  \n- · For ˜ c as small as 0 . 05 , our reducedc and fullc , re-scaled ˙ N em simulations agree in the global reionization history within a linear difference of ≲ 0 . 01 at all redshifts.  \nThis confirms that that the re-scaling of ˙ N em prescribed by Eq. 2.11 is equivalent to using the RSLA with respect to the global reionization history.  \n- · We have compared the morphology of ionized regions and the 21 cm power spectrum in our reducedc and fullc , re-scaled ˙ N em runs. We find that the sizes of the largest ionized regions are suppressed in the reducedc runs relative to their fullc , re-scaled ˙ N em counterparts. This is caused by position-dependent time-delay effects caused by the RSLA for which Eq. 2.11 does not account. The net effect is to suppress ionization power in the reducedc simulations. We find that differences in the 21 cm power spectrum at k = 0 . 1 h Mpc -1 is at most 15% ( 10% ) for ˜ c = 0 . 2 ( 0 . 3 ), but can exceed 20% for ˜ c = 0 . 1 and 30% for ˜ c = 0 . 05 . As such, we recommend applying Eq. 2.11 only for ˜ c ≥ 0 . 2 when ionization morphology is important.\n- · We have made the same comparison with respect to observables derived from highredshift quasar spectra, with a focus on the Ly α forest at z ≤ 6 . We find that when ˙ N em is calibrated to reproduce the mean transmission of the Ly α forest at z ≤ 6 , the resulting ˙ N em histories are very different for different ˜ c . For reducedc runs calibrated in this way with ˜ c = 0 . 2 , and 0 . 3 , we find that their fullc , re-scaled ˙ N em agree in the mean transmission to 20% or better. However, for ˜ c = 0 . 1 the difference can be as large as a factor of 2 . We find qualitatively similar results for the ionizing photon mean free path and the IGM temperature at mean density. Both agree to within ≤ 5% for ˜ c = 0 . 2 and 0 . 3 , but differ by up to 35% and 6% , respectively, for ˜ c = 0 . 1 .\n- · We looked at the large-scale fluctuations in the Ly α forest opacity in our two sets of simulations. For ˜ c ≥ 0 . 2 , the distribution of effective optical depths over 50 h -1 Mpc segments of the forest agrees to within a few percent between reducedc and fullc , re-scaled ˙ N em simulations. Once again, we find that these differences are considerably larger for ˜ c = 0 . 1 . We explored this result in more depth by looking at the large-scale fluctuations in Γ HI and T during the end stages of reionization. reducedc simulations over-produce large-scale fluctuations in Γ HI relative to their fullc , re-scaled ˙ N em counterparts thanks to the same position-dependent time-delay effects responsible for affecting the ionization morphology. In addition, reducedc simulations slightly underestimate the temperatures of recently ionized, hot voids near reionization's end, owing to their lower mean I-front speeds.\n- · Lastly, we generalized our method to apply to multi-frequency RT simulations. We found that with respect to the Ly α forest, the method works just as well with multifrequency simulations as with single-frequency ones. However, because of IGM filtering effects, ˙ N c em ends up having a softer spectrum than ˙ N ˜ c em , complicating the physical interpretation of the re-scaling procedure. We caution that when applying our method to multi-frequency simulations, this effect should be quantified and its implications carefully considered.  \nThe approach described in this work is potentially useful in the following situations: either (1) when ˙ N em is entirely a free function of redshift that is calibrated to match some observable (as in Refs. [14, 23]) or (2) when ˙ N em is uniquely determined by some set of free parameters that are being marginalized over, as in Ref. [30]. In the first case, ˙ N em can be calibrated using the RSLA at a reduced computational cost, and Eq. 2.11 can be applied to the end result to recover the 'true' ˙ N em . In the second case, the re-scaling prescribed by  \nEq. 2.11 would be equivalent to changing the mapping between ˙ N em and the free parameters that determine it in a quantifiable way. Unfortunately, our method is not directly applicable for situations where ˙ N em cannot simply be re-scaled after the simulation has been run. This is the case if ˙ N em is a self-consistent prediction a model that has no or very few free parameters, such as the THESAN [7, 19], CoDa [5, 13], or SPHINX [48, 49] simulations. Our results further suggest that for the first two applications, the RSLA should only be used (in conjunction with Eq. 2.11/5.1) for ˜ c ≥ 0 . 2 , to ensure that the effects of the RSLA on large-scale fluctuations in the ionization and radiation fields are minimized. This corresponds to a speed-up factor of up to 5 , which will help enable considerably faster searches of the reionization parameter space using RT simulations.", 'Acknowledgments': "The author acknowledges support from the Beus Center for Cosmic Foundations while this work was ongoing. He also acknowledges helpful conversations with Anson D'Aloisio, Garett Lopez, and Shikhar Asthana.", 'References': "- [1] X. Fan, M.A. Strauss, R.H. Becker, R.L. White, J.E. Gunn, G.R. Knapp et al., Constraining the Evolution of the Ionizing Background and the Epoch of Reionization with z ~ 6 Quasars. II. A Sample of 19 Quasars , The Astronomical Journal 132 (2006) 117 [ astro-ph/0512082 ].\n- [2] B.E. Robertson, R.S. Ellis, S.R. Furlanetto and J.S. Dunlop, Cosmic Reionization and Early Star-forming Galaxies: A Joint Analysis of New Constraints from Planck and the Hubble Space Telescope , The Astrophysical Journal Letters 802 (2015) L19 [ 1502.02024 ].\n- [3] G. Kulkarni, L.C. Keating, M.G. Haehnelt, S.E.I. Bosman, E. Puchwein, J. Chardin et al., Large Ly α opacity fluctuations and low CMB τ in models of late reionization with large islands of neutral hydrogen extending to z < 5.5 , Monthly Notices of the Royal Astronomical Society 485 (2019) L24 [ 1809.06374 ].\n- [4] N.Y. Gnedin, Cosmic Reionization on Computers. I. Design and Calibration of Simulations , The Astrophysical Journal 793 (2014) 29 [ 1403.4245 ].\n- [5] P. Ocvirk, D. Aubert, J.G. Sorce, P.R. Shapiro, N. Deparis, T. Dawoodbhoy et al., Cosmic Dawn II (CoDa II): a new radiation-hydrodynamics simulation of the self-consistent coupling of galaxy formation and reionization , Monthly Notices of the Royal Astronomical Society 496 (2020) 4087 [ 1811.11192 ].\n- [6] L.C. Keating, L.H. Weinberger, G. Kulkarni, M.G. Haehnelt, J. Chardin and D. Aubert, Long troughs in the Lymanα forest below redshift 6 due to islands of neutral hydrogen , Monthly Notices of the Royal Astronomical Society 491 (2020) 1736 [ 1905.12640 ].\n- [7] E. Garaldi, R. Kannan, A. Smith, V. Springel, R. Pakmor, M. Vogelsberger et al., The THESAN project: properties of the intergalactic medium and its connection to reionization-era galaxies , Monthly Notices of the Royal Astronomical Society (2022) [ 2110.01628 ].\n- [8] N.Y. Gnedin and T. Abel, Multi-dimensional cosmological radiative transfer with a Variable Eddington Tensor formalism , New Astronomy 6 (2001) 437 [ astro-ph/0106278 ].\n- [9] D. Aubert and R. Teyssier, A radiative transfer scheme for cosmological reionization based on a local Eddington tensor , Monthly Notices of the Royal Astronomical Society 387 (2008) 295 [ 0709.1544 ].  \n- [42] L.C. Keating, G. Kulkarni, M.G. Haehnelt, J. Chardin and D. Aubert, Constraining the second half of reionization with the Ly β forest , Monthly Notices of the Royal Astronomical Society 497 (2020) 906 [ 1912.05582 ].\n- [43] J. Chardin, M.G. Haehnelt, D. Aubert and E. Puchwein, Calibrating cosmological radiative transfer simulations with Ly α forest data: evidence for large spatial UV background fluctuations at z ∼ 5.6-5.8 due to rare bright sources , Monthly Notices of the Royal Astronomical Society 453 (2015) 2943 [ 1505.01853 ].\n- [44] F.B. Davies and S.R. Furlanetto, Large fluctuations in the hydrogen-ionizing background and mean free path following the epoch of reionization , Monthly Notices of the Royal Astronomical Society 460 (2016) 1328 [ 1509.07131 ].\n- [45] A. D'Aloisio, M. McQuinn and H. Trac, LARGE OPACITY VARIATIONS IN THE HIGH-REDSHIFT LY α FOREST: THE SIGNATURE OF RELIC TEMPERATURE FLUCTUATIONS FROM PATCHY REIONIZATION , The Astrophysical Journal 813 (2015) L38.\n- [46] T.R. Choudhury, A. Paranjape and S.E.I. Bosman, Studying the Lyman α optical depth fluctuations at z ∼ 5.5 using fast semi-numerical methods , Monthly Notices of the Royal Astronomical Society 501 (2021) 5782 [ 2003.08958 ].\n- [47] C. Zeng and C.M. Hirata, Nonequilibrium Temperature Evolution of Ionization Fronts during the Epoch of Reionization , The Astrophysical Journal 906 (2021) 124 [ 2007.02940 ].\n- [48] J. Rosdahl, H. Katz, J. Blaizot, T. Kimm, L. Michel-Dansac, T. Garel et al., The SPHINX cosmological simulations of the first billion years: the impact of binary stars on reionization , Monthly Notices of the Royal Astronomical Society 479 (2018) 994 [ 1801.07259 ].\n- [49] H. Katz, S. Martin-Alvarez, J. Rosdahl, T. Kimm, J. Blaizot, M.G. Haehnelt et al., Introducing SPHINX-MHD: the impact of primordial magnetic fields on the first galaxies, reionization, and the global 21-cm signal , Monthly Notices of the Royal Astronomical Society 507 (2021) 1254 [ 2101.11624 ]."}