A Full Treatment of Peculiar Velocities on the Reionization Light Cone
Emma Chapman Mario G. Santos


Surface-acoustic-wave-controlled optomechanically induced transparency in a hybrid piezo-optomechanical planar distributed Bragg-reflector-cavity system
    Zhong-Yang Wang
    
========================================================================================================================================================




Accurate simulations of the 21-cm signal from the Epoch of Reionization require the generation of maps at specific frequencies given the values of astrophysical and cosmological parameters. The peculiar velocities of the hydrogen atoms producing the 21-cm radiation result in a shift in the observed frequency of the 21-cm radiation and changes the amplitude of the signal itself. 
This is not an effect we can remove but instead needs to be accurately modelled to ensure we infer the correct physical parameters from an observation. 
We describe the full calculation of the distortion of the 21-cm signal, and propose a new code that integrates the 21-cm intensity along the line of sight for each individual light cone pixel to fully capture the intensity contributions from multiple redshifts. This algorithm naturally deals with the typical divergences found in standard approaches, allowing for large optical depths and 21-cm absorption events at high redshifts. 
We find the new method results in up to a 60% decrease in power on the largest scales on the sky, and an increase of over 80% on the smallest scales on the sky. We find that the new implementation of the light cone results in a longer tail of bright temperatures in the brightness temperature distribution, as a result of the successful circumventing of a previous cap that had to be implemented to avoid a divergence in the brightness temperature. We conclude that this full treatment of the evolution of the light cone pixel can be an important effect.




cosmology: theory – dark ages, reionization, first stars – diffuse radiation – methods: statistical.




§ INTRODUCTION


380,000 years after the Big Bang electrons and protons in our
Universe combined into neutral hydrogen and the Universe became
transparent to photons. As the Universe continued to expand and cool it
settled into its `Dark Ages', as the underlying dark matter structure
formed and the neutral hydrogen was gravitationally attracted to the
densest clumps of dark matter. Eventually, the density of this hydrogen was
such that fusion could begin and the first stars were
born. The Epoch of Reionization is the time following the Dark Ages
when the first ionizing sources in our Universe, whether mini-quasars,
stars or galaxies, began to emit radiation and ionize the
surrounding hydrogen. Eventually these patches of ionized hydrogen
grew such that they overlapped and the majority of the Universe was
ionized, by around one billion years after the Big Bang. 

Up until recently, the study of this field remained mostly theoretical,
apart from a few indirect constraints.  The Thomson optical depth as
measured by Planck can give us an integral constraint relating to the
transparency of the Universe to photons. The latest constraint is
τ=0.054 ± 0.007 <cit.>, which is rather lower than
previous estimates from WMAP, relieving the tension from previous higher
values and setting an upper constraint on the extent of reionization to Δ z < 2.8. Quasars are another key probe of reionization. Quasar spectra are sensitive to intervening neutral hydrogen with the Lyman-alpha (Lyα) resonance in particular causing absorption lines in the broadband spectra of the quasar. These lines combine over redshift to form the Ly-α forest or, at higher redshifts when the forest saturates, a Gunn-Peterson trough, giving us an indication of the neutral fraction of the Universe <cit.>. Quasar absorption lines can be used to study the global ionizing emissivity budget, metal buildup and timing of reionization together implying that reionization comes to an end around 5.7 ≤ z ≤ 7.7 <cit.>. More specifically, quasars observed at very high redshifts can provide their own bounds for the neutral fraction at those redshifts, albeit along that line of sight only (e.g. <cit.>, <cit.>).

The 21-cm spin-flip transition of neutral hydrogen, either
in the form of an absorption or an emission relative
to Cosmic Microwave Background (CMB) blackbody
spectrum, holds the promise of revealing the detailed
astrophysical processes occurring during the Epoch of
Reionization (EoR) and the Cosmic Dawn (CD). 
For a review see e.g. <cit.>.
Several low frequency radio telescopes plan to observe this signal in order to detect the EoR (e.g. LOFAR[http://www.lofar.org/], MWA[http://www.mwatelescope.org], GMRT[http://gmrt.ncra.tifr.res.in], PAPER[http://eor.berkeley.edu], HERA[http://http://reionization.org] and the future and SKA1-LOW[http://skatelescope.org]). This neutral hydrogen signal extends back into the higher redshifts, when the first ionizing objects are still forming, in the Cosmic Dawn. Recently, a claim has been made of a possible detection of the absorption feature in the 21-cm global signal from the Cosmic Dawn by the EDGES telescope[https://www.haystack.mit.edu/ast/arrays/Edges/] <cit.>. This detection is still open to interpretation and an independent confirmation is required to ensure the absorption signal is real and not, for example, a conspiracy between the frequency dependence of the instrument and foreground contamination. Nevertheless, the wealth of information that we soon expect to access from all these experiments makes it the all more important to make sure we have an accurate calculation of the observed signal.


The 21-cm signal depends on several ingredients that can be calculated through simulations, in particular the fraction of neutral hydrogen and the spin temperature (more relevant at higher redshifts). 
Several codes have been developed for this, either full radiative transfer (e.g. ) or fast numerical simulations (e.g. ). However, turning these values into the observed 21-cm signal requires in principle a full radiative transfer calculation of the background radiation along the line of sight which depends on subtle and interconnecting effects. Firstly,
the peculiar velocities of the hydrogen atoms result not
only in the red-/blue-shifting of the 21-cm radiation but a boost/suppression of
the signal by a significant magnitude <cit.>. 
The resulting anisotropy in the power spectrum has the potential to reveal
rich cosmological and astrophysical
information (e.g. , , ,
, ) and so needs to be
carefully understood.
Another observational effect is the cosmological
redshift of the 21-cm photons due to the expansion of the Universe,
λ_obs = λ_21 (1+z), which results in a set of observations
in frequency corresponding to not only different positions in space but also to different times in the evolution of the Universe. This is the so called light
cone effect. In order to construct a light cone one must take only the information from photons which have traveled the same light travel time. Analyzing a given observed volume without taking this evolution over time into account will produce biased results.

Both the light cone effect <cit.> and
the effect of the peculiar velocities on the 21-cm signal <cit.> have been well studied. However, the meeting of the two is non standard. The usual approach is for brightness temperature cubes to be created at a given redshift with the peculiar velocity contribution included in the calculation and then, in a separate step, a light cone is assembled. This approach leads to problems when the peculiar velocity gradient is large, even leading to infinities that have to be dealt with in an ad hoc way. Moreover it doesn't take into account the fact that there could be contributions from multiple real space simulation box slices to the same map since the light cone map and redshift space distortion effects are coupled. <cit.> describes an algorithm to deal with the peculiar velocity effects and the above "infinities" but although the coupling to the light cone is mentioned in Sec. 6.2.2, it wasn't implemented. One recent exception is <cit.> who divided the simulations into thin shells, re-assigning 21-cm events to an observed frequency including the effects of peculiar velocities. However, their approach does not consider that a 21-cm event may only partially contribute to a frequency map, having redshifted out of the observed frequency range as the Universe expands, and they also only consider the signal in emission at low redshifts, when the spin temperature is much larger than the CMB.

In this paper, we revisit the effect of the peculiar velocities on the light cone and present a new algorithm that naturally deals with the possible large gradients in the peculiar velocities, avoiding the issue of infinities in the brightness temperature altogether. The code allows for multiple scattering events along the line of sight (e.g. large optical depths) and evolution of the physical quantities along the line width when such width translates to a large time interval in the evolution of the Universe. Moreover, it includes the calculation of these effects during absorption, at high z, and is appropriate for fast semi-numerical simulations.

In Section <ref> we review the basic physics of the 21-cm radiation, including the standard 21-cm brightness temperature calculation, introduce the light cone effect and discuss
the common way of including this effect, uncoupled to peculiar velocities. In Section <ref> we rederive the brightness temperature equation considering how the 21-cm intensity adjusted for redshift space distortions is observed due to the light cone effect, thus taking into account the inherently coupled nature of the two effects. The new open source algorithm resulting from this new derivation is compared with other methods in the literature in Section <ref>. We summarise the semi-analytic
code which we are adapting in Section <ref>. In Section <ref> we compare images and power spectra of the new algorithm with the traditional uncoupled treatment of the light cone and peculiar velocity effects.





§ THE 21-CM SIGNAL - STANDARD CASE




We start by re-deriving the fundamentals of the 21-cm signal for the
standard case. Although this has been discussed previously in the literature (see e.g. ), we believe it will be helpful to review the steps again in detail, in order to make it clear what approximations are usually made and what is new in the algorithm we are introducing.

The 21-cm signal corresponds to the intensity, I_21(ν_0), observed at frequency ν_0, due to the combined effects of emission and absorption of radiation along the line of sight by neutral hydrogen from its 21-cm line. This line is particularly useful to probe the high redshift Universe, because of the abundance of neutral hydrogen as well as the weak energy associated with the line, which can be easily excited without the need for strong (and therefore rare) background sources. The observed intensity is usually considered against the background radiation, I_ bg(ν_0) (the radiation one would observe if there was no intervening HI), e.g. I_21(ν_0)-I_ bg(ν_0), so that it can be positive or negative depending on the amount of background radiation absorbed versus the radiation emitted due to the line transition. Although some lines of sight might have a strong background source, this should be very rare on cosmological scales and the majority of the signal will have the cosmic microwave background as the background. It is common to define the brightness temperature:

    δ T_b(ν_0)≡ I_21(ν_0)λ_0^2/2 k_B - T_ CMB(z=0)

where λ_0 is the wavelength of the observed radiation, k_B the Boltzmann constant and T_ CMB is the CMB temperature today (2.725 K). Note that the calculation neglects angular fluctuations in the CMB temperature itself as the fluctuations are smaller than the brightness temperature fluctuations. From now on we will concentrate on the observed intensity I_21(ν_0) noting that equation <ref> should be considered when referring to the 21-cm signal brightness temperature. 





 §.§ Changes in intensity


The basic quantity we need to calculate is the variation in intensity dI as the "light ray" crosses a patch of neutral hydrogen (note that the basic radiative transfer relations we use can be found in standard textbooks, for example, ). This change can be due to three effects: i) absorption/scattering of photons in the 21-cm resonance along the patch; ii) emission of 21-cm radiation from the patch and iii) Doppler effects due to the expansion of the Universe and peculiar velocities. Throughout the paper, we will call a 21-cm emission/absorption along the path of the light ray towards the observer, a "21-cm event". Ignoring for the moment the last contribution, the change in intensity at a frequency (ν) across an infinitesimal positive line element ds, can be written as:


    dI(ν) = -k_21(ν) I(ν) ds + j_21 (ν) ds,
 

where j_21 is the emissivity coefficient and k_21 is the absorption coefficient. The first term on the right accounts for absorption and the last term for emission (ν is the frequency during the event, which
needs to be within the 21-cm line). We will always consider the changes along a given line of sight as the radiation propagates towards the observer. Therefore, s will be the line of sight "proper" length, but starting at some far away distance and increasing towards the observer (e.g. positive ds as it goes down the line of sight towards the observer). This is opposite from the standard definition of line of sight distance.

Following , the emissivity j_21 and absorption coefficient k_21, can be written to a good approximation as:


    j_21 (ν) ≈ [1.6 × 10^-40 J  s^-1 sr^-1] n_ HIϕ_21 (ν )


and


    k_21(ν) ≈ (2.6 × 10^-19  K Hz m^2)  n_ HIϕ_21(ν)/T_S,


where ϕ_21(ν) is the 21-cm line profile (with the integral over frequency normalized to one), T_s the spin temperature and n_ HI is the number density of HI atoms in total (in proper units unless otherwise stated). 
Note that both k_21 and j_21 are functions of frequency, time and spatial position. They are only non-zero when there is neutral hydrogen at the position we are considering and when the frequency is within the 21-cm line.

In an expanding Universe, the frequency of the light ray will redshift/blueshift due to the uniform expansion and the peculiar velocity. Again for an infinitesimal positive line element ds (in proper units), going down the line of sight towards the observer, the frequency of the photon is adjusted by dν:


    dν =- νds/c[H(z) + dv/ds],


where H(z) is the Hubble parameter at the redshift considered and dv/ds is the "proper" gradient of the peculiar velocity along the line of sight (both dv and ds are in proper units). Note that the gradient will hold the same value whether we take ds increasing or decreasing away from the observer (since the axis will be reversed). For |dv/ds| < H, the term in brackets will be positive and ν will decrease when s increases towards the observer, as expected (the usual redshift).

The intensity of the "light ray" will also change from an initial (I_i) to a final value (I_f = I_i + dI) according to this shift in frequency from ν_i to ν_f = ν_i + dν: I_f(ν_f) =  (ν_f/ν_i)^3 I_i(ν_i). Using the equation above for an infinitesimal ds, we obtain:

    dI(ν) = - 3 I(ν) ds/c[H + dv/ds].

To first order in ds, we can take the frequency at the beginning or the end of the patch ds. Neglecting the peculiar velocity term, this translates into the standard evolution: I∝ a^-3 (where a is the scale factor and we note that according to our current definition, we have ds=c dt). In principle, this shift in intensity along ds should be added to equation <ref>, though we will consider if this is necessary in the next subsection. 






 §.§ Approximations


Using equations <ref>, <ref> and <ref> above, one can in principle integrate the effects on the light ray along each line of sight in order to calculate the observed 21-cm signal. There are however several approximations usually made in the literature which we can now incorporate:


  * Take the 21-cm line profile ϕ_21 to be a top hat function centred at the ν_21 frequency, e.g. ϕ_21(ν) = 1/Δν_21 across the 21-cm line, where Δν_21 is the line width.

  * Assume there is only one absorption or emission 21 cm event for a given observed frequency ν_0. This would happen at redshift z=ν_21/ν_0 -1.

  * Assume |dv/ds|/H ≪ 1, Δν_21/ν_21≪ 1 and any physical quantities (such as j_21 and k_21) are constant across the line width. 


From equation <ref>, the proper distance traveled by the light ray across the event can now be written as 


    Δ s=Δν_21/ν_21cH^-1/1+1/H dv/ds.
 

where we are already taking the frequency along the event as the 21-cm line frequency. The natural width of the 21-cm line is quite small ∼ 10^-16 Hz, but thermal broadening will increase this by a large factor, Δν_21∼ν_21/c√(2 k_B T_K/m_H) (with T_K, the temperature of the gas and m_H the hydrogen atom mass). During reionization, the gas will be heated to temperatures of ∼ 10^4 K, which will generate line widths of up to ∼ 60 kHz. Still, the quantity Δν_21/ν_21 will be quite small (< 10^-4). The same cannot be said of |dv/ds|/H as we will see later.

Going back to equation <ref>, we can rewrite the change in intensity across the 21-cm event as


    Δ I(ν) = j_21 (ν) Δ s [1 - I/S_21],


where we are now using Δ to emphasize that we are considering a non-infinitesimal patch since we are considering the event across the full line width as explained in the assumptions above. Here, S_21≡j_21/k_21 is the specific intensity of the 21-cm transition, which is related to the spin temperature through T_S=S_21λ_21^2/(2 k_B). Putting the values together, we have


    Δ I  = [1.127× 10^-49J s^-1Hz^-1sr^-1] c H^-1 n_ HI/1+1/H dv/ds[1 - I/S_21],
 


One last approximation is related to the shift in intensity in an expanding Universe, according to equation <ref>. This is obviously the main effect when there is no 21-cm event but it needs to be checked if we can neglect it during the event itself. During the emission/absorption, the largest frequency shift will be the length of the 21-cm line. We can then write Δ I/I ∼ 3 Δν_21/ν_21∼ 10^-4 (at most). For the situation where the initial radiation is the CMB, the contribution to emission from the 21-cm event will be Δ I/I ∼ 5× 10^-2, which is much larger. If the I/S_21 term dominates, the contribution would be even larger. Therefore, it is safe to ignore the shift in intensity due to the shift in frequency across the event and we will do so throughout the rest of the paper. The bottom line for this is that, after we compute the emission/absorption across each event, we just need to rescale the intensity by the factor (ν_0/ν_21)^3 to take into account the expansion effect (and take the initial frequency to be just the 21-cm one).

As described, the standard approach only assumes one event across the line of sight. Therefore the initial intensity is just the CMB one, e.g. Δ I = I_f - I_i = I_f - I_ CMB. Using the Rayleigh Jeans approximation to transform the intensity into a temperature, the observed brightness temperature is finally:

    δ T_b(ν_0)    ≈    (1.82 × 10^-28  K m^2) n_ HI/1+z×
          cH^-1/1+1/H dv/ds[1 - T_ CMB/T_S],


where the relevant quantities are calculated at the redshift z=ν_21/ν_0-1. Using typical values, this will give an observed signal of about ∼ 14 mK from emission at z=10 with the Universe 50% ionized. 

The equation above represents the standard calculation used in practically all 21-cm papers and simulations. Using for instance the simulation box outputs at a given z (n_ HI, dv/ds and T_S) we can generate an observed δ T_b box for this same z. We can then perform any analysis on it, such as the power spectrum statistics. However, as discussed below, even in the framework of these approximations, such an approach is only valid for a small interval in redshift as otherwise there will be cosmic evolution (for instance, an 8 MHz interval corresponds to about dz∼ 0.5, where evolution can already be non-negligible).







 §.§ The standard light cone





The so called light cone effect results from the fact that positions along the
line-of-sight are observed in the frequency direction. In this way, the observer cannot build a volume at a single moment in time. While simulations will output a 3D realization of the Universe at a chosen redshift,
observations will build up a volume where the light from positions farther away, corresponding to a lower frequency, will be from an earlier stage of reionization. In this way, statistical quantities constructed from the simulation outputs, such as the power spectrum, will only be an accurate representation of the observations for narrow slices where cosmic evolution can be ignored. A more accurate representation is to directly construct the light cone from simulations. The basic method employed to construct such a simulated observation cube from a series of simulation boxes is as follows (also refer to Fig. <ref>):



  * Consider a series of simulation cubes representing the Universe
  between redshifts z_ min and z_ max, all of comoving side length L
  and made up of N cells along one side. We label the cells along the
  line-of-sight as k_i where 0 ≤ i < N and k_0 is nearest to the observer. 
  The box comoving resolution is Δ r = L/N. 

  * The simulation boxes all occupy the same region of comoving space in the Universe and we
  place the nearest edge (i.e. k_0) of all the boxes at a comoving
  distance x_min, which is the distance to the smallest output redshift, z_ min, with x defined through:

    x(z) = c/H_0∫^z_0 dz'/√(Ω_m(1+z')^3+Ω_r(1+z')^2 + Ω_Λ).

In the absence of peculiar velocities, each frequency of observation, ν,
corresponds to a simulation box output at a certain
redshift through the standard relation 1+z =
ν_21/ν. If there is no simulation box at that precise
redshift, one can be obtained through interpolation.


  * From the box z_ min, we take the nearest slice along the
  line-of-sight, k_0, of the simulation
  box, and this becomes the slice at the highest frequency, ν_ max, in our
  light cone.

  * For observed frequency ν_i, the 21-cm signal corresponds to 
  redshift z_i:

    1+z_i =
    ν_21/ν_i
 and at a comoving distance of x_i. We then
load box z_i and take the slice at comoving distance x_i. This
will have a cell address of:

k_i = (x_i - x_min)/Δr .

The slice located at line-of-sight address k_i becomes the observed
frequency map at frequency ν_i. 

  * In order to build up a light cone relating to observed
  frequencies between ν_ min and ν_ max, a simulation box size L
  is required such that L > x(z(ν_ min)) - x(z(ν_max)), though this
  can of course be satisfied by using boundary conditions.


We refer to the light cone
built up using this prescription from the real space brightness temperature boxes, as the standard light cone. 




§ THE 21-CM SIGNAL - FULL APPROACH




As we saw in the previous section, standard methods to calculate the 21-cm signal follow two separate processes:


  * Calculate brightness temperature boxes using equation <ref>.

  * Either ignore cosmic evolution and just proceed the analysis using the box at a given average z or slice different boxes to construct a light cone as described in section <ref>.


This approach ignores that formally the brightness temperature is only well defined from the point of view of the observer, e.g. the signal does not exist in isolation at any redshift. Therefore, a proper algorithm should generate the line of sight cone and brightness temperature at the same time.
Moreover, as we can see in equation <ref>, the result is not well defined when 1/H dv/ds = -1 which can occur in high density places. In fact, even for other values close to -1 the situation is ill defined as we will discuss later, so that a more in-depth algorithm should be used. 







 §.§ Integration along a line of sight



To calculate I_21(ν_0) from a given direction n̂, we should consider the initial background radiation, I_ bg(ν) at some high redshift z and then integrate the effects of absorption and emission along the line of sight n̂ from redshift z to 0. The chosen starting redshift should be high enough to guarantee that there is no 21-cm event before it, so that every event is counted. The initial intensity will be at a comoving distance from the observer given by x(z) as above.

As we integrate along the line of sight from z to 0 (or from x(z) to x=0), we need to calculate the small changes in intensity dI across the small steps dx. This is usually referenced instead in terms of the proper physical distance ds=dx/(1+z). The change dI will depend on the position s (or redshift) and incident intensity as well as the frequency. The question is then how to generate the signal for a given angle, n̂ and observed frequency, ν_0. Our approach goes through the following steps:


  * For a given observed frequency, ν_0, consider a high enough initial redshift (1+z_i > ν_21/ν_0) so that we are certain that there is no 21-cm event above that redshift contributing to this frequency (it is far away from the 21-cm line). The actual value was obtained by comparing to simulations and will be discussed later.

  * The initial redshift will correspond to a comoving distance x(z_i). We then go down along the line of sight towards the observer, looping through x_j+1=x_j-dx for a fixed, positive dx. This corresponds to a proper interval ds=dx/(1+z). Redshifts will change accordingly: z_j+1=z_j-Hds/c (1+z_j).

  * The initial starting frequency is given by ν_i = ν_0 (1+z_i) [1+(v_i-v_0)/c], 
where v_i is the peculiar velocity at the initial point and v_0 is the observer velocity which is just a dipole and can be ignored.

  * The frequency is then updated through equation <ref>: ν_j+1 = ν_j [1 - (H(z) + dv/ds) ds/c].

  * We then go through the line of sight until the frequency intersects with the 21-cm line and we have an event.

  * The initial intensity is just: I_i(ν_i) = I_CMB(ν_i). However, if there is no 21-cm event, the intensity will just change due to the normal expansion of the universe. Once there is a 21-cm event, the initial intensity (before any event) can be assumed to be the CMB value at the 21-cm frequency, I_CMB(ν_21). It will then be updated according to equation <ref> We can safely neglect the redshift effect on the intensity during the event as discussed in the previous section. This update is done for all the events along the line of sight, using the intensity value from the previous step.

  * At low enough redshifts, 1+z_f < ν_21/ν_0 (once we are sure there will be no more events) we stop the integration along the line of sight. The final intensity is then updated through I_21(ν_0) = (ν_0/ν_21)^3 I_f(ν_21), where I_f can be the result of several 21-cm events (integrated using equation <ref>). Note that at any event, the frequency for I will always be ν_21 so that no frequency update is needed between events, except for the final update when moving to the observed frequency ν_0.

  * The final brightness temperature is then calculated using equation <ref>.


For a given observing frequency, in principle we would need to do the calculation along the whole line of sight as peculiar velocities could move the frequency into the 21-cm line. This would be too slow and, as already indicated in the steps above, we only integrate between a given redshift interval. The target redshift is, as usual, 1+z=ν_21/ν_0 and factoring in the line width and peculiar velocities, we only expect 21-cm events for an interval of:

    Δ z/1+z≈Δν_21/ν_21+2v/c,

where v/c should correspond to the largest expected values of the peculiar velocity. This term dominates over Δν_21/ν_21. In terms of the observed frequency resolution this would correspond to Δν∼ν_21/(1+z)(2v/c). From the simulation, we get values at most of about v/c∼ 1.5× 10^-3 (see Fig. <ref>). At z=8 this would correspond to potentially filling maps in a Δν∼ 0.5MHz range. This will decrease with increasing redshift. To be conservative, we choose Δν = 0.5MHz throughout the simulation. We do not expect any 21-cm event outside this range. For a frequency map at ν_0 = 150 MHz which would usually correspond to only a redshift of z=8.47, we will integrate events along the line of sight at points within the range z=8.44 - 8.50.



One final issue is the actual calculation of the change in intensity using equation <ref> since the step in frequency might be bigger or smaller than the 21-cm line width itself (since our step size is not infinitesimally small). A fixed comoving step size dx, can correspond to different frequency steps and dν can be even zero if the velocity gradient dv/ds = -H. Again, the range of frequencies for which the 21-cm event will happen is set by the line profile, which we take as a top hat centered around ν_21, with ϕ_21(ν)=1/Δν_21. Considering the intersection of the frequency patch with the line width, the maximum and minimum frequencies for the event  will depend on the frequency evolution. If ν_j+1 < ν_j as it is usually the case, since we are going down in redshift, we have:

    ν_max   =   min{ν_j, ν_21+Δν_21/2}
    ν_min   =   max{ν_j+1, ν_21-Δν_21/2},

If ν_j+1 > ν_j, we have instead:

    ν_max   =   min{ν_j+1,ν_21+Δν_21/2}
    ν_min   =   max{ν_j, ν_21-Δν_21/2}.

These conditions are represented graphically for ν_j+1<ν_j in Fig. <ref>.




The intensity during a 21-cm event is updated through:

    I_j+1 = I_j+j_21(ν_21) ds' [1-I_j/S_21(ν_21)],

where the possibilities for ds' are


  * if there is no event, e.g. neither ν_j and ν_j+1 are within the line width, then ds'=0,

  * if dν = 0, e.g. dv/ds∼ -H (and ν_j is within the line), then ds'=ds,

  * otherwise:

    ds'=ν_max-ν_min/ν_21cH^-1/|1+1/H dv/ds|.



The situation above avoids the case when the patch goes to infinity as in the standard approach (equation <ref>). Note that ds' can never be larger than the chosen ds. Moreover, by using a small ds in the algorithm, we can implicitly avoid the situation where the simulation values can actually change across the patch. Finally note that when the line width is completely "inside" the patch, we have: ν_max-ν_min = Δν_21 and we recover equation <ref>. 





 §.§ Required functions


In order to implement equation <ref> we need to have several ingredients: the hydrogen number density, n_ HI, the line profile, ϕ_21, the specific intensity S_21, the peculiar velocity, v, as well as the peculiar velocity gradient dv/ds. As explained above, the line profile is assumed to be a top hat centred on the ν_21 frequency and with width given by thermal broadening. Therefore, if we want to be more accurate in the use of the line width, we also need the temperature of the gas, T_K.

The specific intensity can be expressed in terms of the spin temperature, S_21(z) = T_S(z) 2 k_ B/λ_21^2 with

    1-T_ CMB(z)/T_ S(z)=x_ tot/1+x_ tot(1-T_ CMB(z)/T_ K(z)),

where x_ tot=x_α+x_ c is the sum of the radiative (x_α≥ 0) and collisional (x_ c≥ 0) coupling parameters (for further details, see e.g. Section 2.4 in ). At "lower" redshifts (z≲ 15), most implementations take T_S(z)≈ T_K(z) and assume T_K(z)∼ 10^4 K. 

We see that in order to calculate the signal, we will need the redshift and space dependent functions: n_ HI, v, dv/ds, T_ K, x_α and x_ c. These will have to be provided by a simulation, which will generate boxes of each component for several redshifts. In our case, we used the output from the code simfast21 as described later, but any other output can in principle be plugged in to our line of sight signal generator.





 §.§ Choosing the resolution


The algorithm described above can be used to calculate the signal for any given observed frequency and angle. However, the observational maps are usually specified in a grid with a fixed angular and frequency resolution. The first question that then needs to be addressed is what resolution to choose for the map we want to construct as well as the optimal integration step ds along the line of sight.

In terms of angular resolution, the optimal choice is a pixel with an angular size equal or slightly smaller than the angular size of the comoving simulation cells at the highest redshift we want to consider. This will guarantee that the signal is the same across the angular pixel (a smaller angular size could be used to take into account any sky curvature effects). In terms of frequency, again the pixel size should be related to the cell resolution of the simulation. A particular issue that might raise some confusion is whether such frequency resolution should take into account the line width, Δν_21. If we consider a single patch, ds, this will indeed generate an observed intensity that is only nonzero for a frequency interval of about Δν_21/(1+z) (not factoring in peculiar velocities). This could imply that assuming the same intensity over the pixel for map pixels of larger sizes in frequency would be wrong. However, observed frequencies separated by more than Δν_21/(1+z) will have contributions from patches along the line of sight which are very similar as long as the separation is smaller than the simulation cell resolution, meaning that indeed, it is the simulation cell resolution that needs to be used when choosing the map frequency resolution. In particular, we should take a pixel frequency size corresponding to the simulation cell size at its highest redshift.

Another issue in the implementation is whether the choice of the 21-cm line width makes a large impact in the final result since our algorithm indeed takes into account this value which can be subject to some uncertainty. In reality, the line width showing in j_21 from Eq. <ref> will usually cancel with the frequency difference in Eq. <ref>. First, as already seen, if the frequency change across ds is larger and surround the line width, then ν_ max-ν_ min=Δν_21. On the other hand, when it is smaller than the line width, further steps down the line of sight will integrate Eq. <ref> up to Δν_21 (assuming all other quantities to be constant) so that the width would again cancel out. This situation also explains when the line width becomes important: when it is large enough that the other quantities in the calculation are no longer constant. That is, when the HI patch contributing to a given frequency becomes larger than the simulation cell size. This is particularly relevant for regions with high optical depth.

As the relevant point is to make sure quantities are constant across the patch, for this implementation we choose the comoving integration step, dx to be a fifth of the simulation comoving cell size. A final check is whether such size is consistent with the assumption of small optical depth across our small ds=dx/(1+z) (which allows us to use equation <ref>). For that, we need k_21 ds ≪ 1. Again, using typical values at z∼ 10 we get k_21 ds∼ 1  K/T_ S (with ds = 0.05 Mpc and a line width Δν_21=60 KHz). So, we should be always in the limit of low optical depth for practical values of the spin temperature. Note that this is true at each ds, not necessarily over the full integration over the line width. 



§ COMPARISON WITH PREVIOUS WORK




The algorithm presented above takes a different perspective on the calculation of the 21-cm signal and it is therefore important to compare to what has been done before. First, due to the intrinsic radiative transfer process of the calculation, discussion of the 21-cm signal should go hand in hand with the light cone implementation. Our approach concentrates on the generation of the observed map. How to then resample this map in order to make it convenient for data analysis is a second stage process. The map itself can, for instance, have much finer resolution in the line of sight direction, or even have a non equally spaced grid. Other approaches tend to work directly on the simulation boxes and usually rely on cubic cells.

Besides our direct light cone/peculiar velocity coupling and map making process, our approach has two other main differences: First, we use the intensity from the previous iteration when integrating the radiative transfer equation, so multiple events can occur along the line of sight. This is specially relevant at high redshifts during absorption, which again, as far as we are aware, as not been dealt with before in terms of the proper RSD calculation. Second, we take into account the width of the 21-cm line, which as we will see below, can be relevant in certain situations.













 §.§ The standard approach


Before we discuss some of the other methods that have been presented to address the peculiar velocities, we would like to first compare to what has been the standard way to deal with this. The basic assumption is the use of equation <ref>. There are then 3 situations where our approach will give different results:


  * When dv/ds≲ -H. In this case the standard approach is to set a minimum threshold to dv/ds to avoid the singularity (or negative temperatures), since this should only occur for high density, non-linear regions, which are expected to be scarce. However, as seen in Fig. <ref>, the number of cells can be quite substantial. Capping the values in such a way, will generate high temperature pixels in the maps and have an impact on the power spectrum.

  * When the assumption of thin optical depth is no longer valid, which is implicit in the standard equation. Our simulation naturally takes this into account by always using a small integration step ds and using the output intensity from one event as input to the next event.

  * When the 21-cm event "spreads" over a large enough region that the quantities are no longer constant (which would usually also mean a large optical depth). According to our calculation, we should only have 21-cm events across a redshift width of Δ z/(1+z) ∼ 0.003. This means that we should only expect differences from this effect if the simulation values change over this width. It also requires that the experiment has a frequency resolution better than this, which is usually the case.


Our method, takes into account all these effects.
Another approach to the standard method would be to realize that, after calculating the values for each cell in the box with a given size Δ r, translating to a frequency map would generate a pixel with frequency size:

    Δν_0/ν_0=Δ r 1+1/H dv/ds/cH^-1.
 
This would generate a map with pixels with different widths in the frequency direction and a rebinning to a common grid with appropriate weights would have to be done. This is the motivation for the method we discuss next.






 §.§ Cell shuffling


Other methods have been proposed to deal with the peculiar velocity effects and avoid the singularity described above. Several schemes for constructing redshift space boxes from simulation data can be found in <cit.>. Although it is referred there that the method can be applied to the arbitrary T_S case, the algorithms described assume T_S≫ T_ CMB. Moreover, the input intensity is always assumed to be the CMB one, so that effectively, multiple scattering events are neglected. In this situation, we can imagine each HI atom as an emitter and the 21-cm signal is basically proportional to the HI number density in redshift space. The scheme that is more readily comparable is the Mesh-to-Mesh Real-to-Redshift-Space-Mapping scheme, e.g. MM-RRM as described in <cit.>. Note again that these methods take the approach of going through the simulation box and calculating the signal instead of moving directly to the map making process as we do.

Similar versions of the MM-RRM scheme have been used in <cit.>, <cit.> and <cit.>. We focus here on the method presented in <cit.>. We assume that the code has produced a series of real space brightness temperature boxes, each representing the same comoving volume of the Universe at different cosmological redshifts. These boxes are assumed to consist of N cells along the line of sight. Their algorithm is as follows:


  * Divide each cell along the line of sight into n sub-cells, with each sub-cell assigned a brightness temperature 1/n'th of the original cell value.

  * Interpolate velocity and density fields onto this subgrid.

  * Shift the brightness temperature sub-cells according to the relation between real and redshift space: s = r+1+z_obs/H(z_obs)v_||(t,r)r̂, where v_|| is the peculiar velocity component along the line of sight.

  * Regrid the brightness temperature box to original resolution.


This algorithm is then carried out for each output box separately. <cit.> takes a slightly different approach in that the boundaries of each subgrid cell are shifted around instead of just the center. In this way, it allows for stretching as we consider in our case. Moreover, when regridding onto the regular, redshift-space space grid, only the overlapping part of the stretched cell is added. The approach in <cit.> is then equivalent to assuming that the sub-grid cells are so small that stretching can be neglected and that they fall completely in one or the other redshift-space cells. A somewhat non-trivial issue that seems to be ignored in these schemes is the intrinsic width of the 21-cm line, which implies that a give HI patch can contribute to different frequencies. 

In order to construct a light cone, the standard algorithm as described in Section <ref> is then used. However, contrary to our map making approach, in this case all points from the same redshift box are assumed to contribute to the same frequency. In reality, there may be events which occur at different cosmological redshifts but are observed at the same frequency due to the peculiar velocity effects. The two stage process (peculiar velocity+light cone) is used in <cit.> where it is pointed out that the light cone correction needs to be made before the peculiar velocity correction.

<cit.> also consider the building of a light cone including the effects of peculiar velocities. Their algorithm go particle by particle within the simulation, assigning an observed frequency including the effects of peculiar velocities. Their approach splits the simulation into thin shells, which can be an issue for large peculiar velocities as it will not be assigned as actually being observed as if produced in another shell. In addition, they only consider the signal in emission at low redshifts when the spin temperature is much larger than the CMB.

In summary, we avoid the traditional two-step process, instead correctly treating the two corrections as inter-dependent, naturally  allowing for large optical depths and absorptions at higher z. 

Our algorithm presents a number of improvements compared to previous implementations: we allow for any value of T_S, we consider multiple scattering events (where the input intensity is the result of a previous event, which can be relevant during absorption) and we implicitly deal with events that go over a "wide" comoving length (or time interval), by doing the integration along the line of sight in fixed ds steps. Moreover we calculate the line width from thermal broadening in each step and factor that in the 21-cm signal calculation.






§ SIMULATION









Though the new light cone code can
be applied to any simulation boxes, we apply it to simfast21[https://github.com/mariogrs/simfast21]. In response to the need for large field-of-view EoR simulations to
compare with observations, the semi-numeric scheme simfast21
was introduced by <cit.>. In the original simulation, while
peculiar velocity gradients were included in the calculation of the brightness
temperature, the effect on the observed light cone was not considered
and the simulation output was a real-space simulation box at each
redshift of interest. Our algorithm introduces a modification to this
code to output a light cone where the peculiar velocities are fully included. This modification replaces the original brightness temperature calculation, leaving all other parts of the code intact. For a full description of the original simfast21 code, <cit.> should be consulted.

In the original simfast21 implementation, after generating all the relevant boxes, the brightness
temperature would also be calculated at each redshift box, taking into account the peculiar velocity gradient:


    δ T_b (z)    =    [23 mK] x_HI(1+δ(z)) ( 1 -
      T_CMB(z)/T_s) 
          ×(h/0.7)^-1(Ω_b h^2/0.02)
    [(0.15/Ω_m
          h^2)(1+z/10)]^0.5
          ×1/1+ 1/H dv/dr
 

As explained before, we see that the equation above breaks down when dv/dr = -aH. This can occur when either the
halo is virialized or at the turn around point of halo collapse just
before virialization.
To avoid singularities, simfast21 defined a velocity gradient limit of dv/dr > -0.7H. Though this only affects a small number of cells, there has been evidence that this prescription could affect
more cells than necessary. <cit.> found that significantly more
cells had  dv/dr < -0.7H than were actually optically thick.
In the new brightness temperature code, we avoid the need for such
stark cut-offs completely by doing the integration along the line of sight for a fixed comoving dx.

We generated simulations to compare both approaches. We start by producing the quasi-linear density, velocity and ionization fields using the standard simfast21 code. We output boxes every 0.5 in redshift between z = 6 and 25. The initial
conditions are created on a 1200^3 grid while a smoothed grid of
400^3 is used for the final boxes. All boxes are 200 Mpc per side and we choose an output light cone field of view of 1 square degree. If a field of view is equivalent to a distance larger than 200 Mpc at a certain redshift, the field of view is filled by assuming the boxes to be tiled on the sky. Since we expect the difference to show up on small scales and along the frequency direction, such field of view is adequate.

The evolution of the neutral fraction is shown in Fig. <ref>. 
We also show histograms of the peculiar velocities and peculiar
velocity gradients at a redshift of 8 in Fig. <ref>. These
do not vary significantly over the redshift range. 


We then create a light cone of frequency resolution 0.1 MHz and 1 degree field of view with 9 arcsec resolution, for both the standard method and our new method. We would like to underline that in this paper we are always comparing our new results to the results using the standard light cone. This guarantees that we are comparing like with like since for instance the power spectrum without the light cone effect can be quite different.
Moreover, at each position along the line of sight x(z), corresponding to a certain redshift z, the values we use are interpolated between the two redshift simulation boxes that bracket the redshift z (taken at the same position in the box). We include this effect in the standard light cone method by interpolating the brightness temperature output from simfast21 as we go down the line of sight. This way we guarantee that any difference in the two methods is only due to the new brightness temperature calculation algorithm. We describe below the differences that can be observed.








§ RESULTS





We show the difference between the standard and new light cones within both the cosmic dawn (90 MHz) and EoR (150 MHz) in Fig. <ref>. Note that since we are only comparing 2D maps here, the standard light cone result is similar to just taking a slice of the coeval simulation box at a given redshift. The only difference will be if we have chosen a redshift between two coeval simulation boxes, in which case the slice would just be a linear interpolation between the two coeval boxes. When comparing the slices we see that the standard light cone maps appear smoothed in comparison to the new light cone maps. While the large scale structure is easily recognizable as similar between maps, the small scale structure is noticeably different, with evidence of significant differences in pixel brightness temperature values. This trend occurs across the frequency range tested.

To look at the difference in pixel values more closely, we plot the histogram of pixel values in the light cone maps for each frequency in Figure <ref>. It is immediately clear that the new light cone method results in longer tails into bright temperatures, especially in emission. This is as a result of being able to consistently include high peculiar velocity gradient pixels in the new implementation. When dv/ds approaches a value of -H, the brightness temperature diverges in the standard derivation of the brightness temperature. This was circumvented in previous implementations of simfast21 by using a cut off preventing dv/dr=-1. The new light cone method naturally accounts for cells approaching this peculiar velocity gradient limit by considering an integration step along the line of sight as in Equation <ref>. When dv/ds=-H then ν_max-ν_min=0 and there is no divergence. In this case, ds'=ds and, considering just the signal in emission, the intensity would be given by j_21Δ s (as in equation <ref>). Assuming this condition to be valid throughout the cell size (0.5 comoving Mpc in our simulation), then the maximum value in emission would be about 140 mK at z∼ 10 (well above the imposed cutoff in the old algorithm). This is consistent with what we see in the histograms, with a few cells in the box with larger values probably due to dark matter fluctuations (or the condition holding for more than 0.5 Mpc).

At 150 MHz and above we see reionization has begun by the lack of pixels showing 21-cm absorption and also the large peak at 0 mK representing ionized regions. The more compact nature of the brightness temperature pixel distribution, as a result of the cut off of these high peculiar velocity gradient pixels in the standard light cone implementation, results in the appearance of smoother light cone maps.

We calculate the 3D power spectrum of the new and standard light cones in order to compare how the simulations change with scale. The 3D power spectrum is defined as Δ_3D^2 = k^3 P(k)/2π^2 V where V is the volume of the sub-band of the light cone being processed, P(k)=<δ(k)δ^*(k)> and δ(k) is calculated using the standard Fourier transform definition. The average is calculated over spherical shells in k. We choose sub-bands of 10 MHz in frequency range for the 3D power spectra and all frequencies referred to in the context of a 3D power spectrum refer to the central frequency of that sub-band.

In Figure <ref> we show the power spectra at 90 and 150 MHz and in Figure <ref> we plot the percentage difference between the new and old light cone ((P_new-P_standard)/P_standard× 100) for a larger range of frequencies. Though there has been some discussion on how the 3D power spectrum risks omitting important information on the EoR <cit.> we include it here for completeness. 
The overall trend is that the difference increases towards larger k (smaller scales). This is consistent with having a larger distribution of pixel values for the new algorithm as seen in the histograms. The extra pixels with large values are also expected to increase shot noise and therefore the power on smaller scales. From this plot it is clear that at various scales we see power losses of up to 45% and power gains of up to 70% when implementing the new light cone approach. This is a significant change in the power spectrum. 

The evolution of this percentage difference with frequency is plotted in Fig. <ref>, highlighting the change in evolution which occurs at a peak at 110 MHz. Comparing to the 110 MHz histogram, we also saw a large difference in the pixel distribution similar to that shown for the 150 MHz histogram shown in Figure <ref>, with a clear maximum value cutoff for the old method. This difference is probably due to the fact that we are at the beginning of reionization, when the the non-linearities become more important but are still in the neutral phase (e.g. not set to zero), so the effect of our algorithm when dv/ds∼ -H is more obvious, allowing for a lot more higher pixel values overall at all scales. We emphasize again that these power spectra are calculated from the standard lightcone or the new light cone output, not simulation boxes at a fixed z.





We plot light cone line-of-sight slices along frequency in Figure <ref>. It is
evident that the general way in which reionization progresses is largely
unaffected by the new light cone. This is not surprising since the 
21-cm events at an observed frequency are not expected to come from
widely different redshift boxes. Moreover, the differences are harder to pick up here due to the x-axis resolution and the dynamic range of the color scale. In any case, we do see stronger differences at the onset of reionization as previously detected. These differences are more easily seen with standard statistics such as the cylindrical power spectrum which we analyse next.





A 2D cylindrical power spectrum in (k_⊥, k_∥) allows us to more easily
separate out the line-of-sight effects of the peculiar velocities. Such separation will be crucial to isolate the different astrophysical components from the cosmology as discussed in <cit.>.
The power spectrum is calculated by 3D Fourier transforming the 10 MHz light cone volume
and binning the pixels according to both the perpendicular k_⊥
and parallel k_∥ Fourier scales.
 In Fig. <ref> we show the percentage difference between the cylindrical power spectrum of the standard and new light cones. While the patterns of signal inflation and signal loss across the scales is consistent between the spherical and cylindrical power spectra, it is additionally apparent that these gains and losses are not symmetrical in k-space.
 


Anisotropies in the power spectrum are often measured in terms of
(k,μ) where μ is the angle with the line-of-sight, μ =
k_∥/|k|. For an isotropic power spectrum we would expect P_μ(k,μ)
to be a flat line for a given k.

The standard and new light cone power spectra, Δ^2_3D(k,μ), are shown in Fig. <ref>. We see power loss at almost all scales, with the exception of 110 MHz. This is consistent with what we have seen before since here we only consider scales k<1.5 h/Mpc.
The shape of the power spectrum remains the same except at 70 MHz and 150 MHz where we see larger differences in directions approaching μ=0. To investigate this further we plot the μ=0 bin only of the power spectrum for the standard and new light cones and the percentage difference between them in Figure <ref>. We see that for most frequencies there is a drop in the power in the μ=0 bin across all scales. For 70 MHz and 110 MHz we see a power increase on the largest k scales. 

It is clear that any fitting of a μ-decomposed polynomial
such as described in <cit.> would result in differing
values for the cosmological and astrophysical power spectrum terms
motivating this full inclusion of peculiar velocities into the light cone.







§ CONCLUSIONS

We have presented a new open source code (<https://github.com/chapmanemma/RSD_LC_github>) which calculates the redshift space light cone from the standard real space simulation output of any reionization
code, while properly taking into account the coupling with peculiar velocities. We set out a physically motivated treatment of peculiar velocities, both in the 21-cm intensity equation and in the calculation of the observed
frequency for each 21-cm event. This requires the individual assessment
of each location along a line-of-sight according to whether a 21-cm event would be observed at
the frequency of observation and how the intensity must be adjusted as
a result of the interaction of cosmological redshift and the peculiar velocity. We have shown that this full treatment produces significant statistical differences when compared to the light cone assembled by taking slices from coeval boxes equivalent to the redshift-frequency conversion. 

Our algorithm naturally deals with high velocity gradients, large optical depths/multi-scattering events and is applied at low frequencies/high redshifts during absorption. It also factors in the line width in the calculation. We find that the new method increases the number of bright pixels present in the light cone due to the natural dispensing of the divergence in the brightness temperature equation, when dv/ds = -H. Previously, simfast21 implemented a cut off on cells approaching these values, resulting in a loss of bright temperature cells in the light cone. With the new implementation we no longer lose this information. 

We find a general trend of power
loss with the application of the new light cone on the larger scales on the sky, while at the smallest scales we see a power increase. We see an overall power increase at 110 MHz (though with the same trend) which we believe is connected to the early stage of reionization.

We conclude that to be
confident of any analysis made from simulated light cone data, one must
be sure to fully incorporate the peculiar velocities as shown here. Imposing a hard cutoff on the peculiar velocities will dramatically change the pixel distribution as seen in the histograms.
In future work we plan to use a higher resolution simulation and multiple reionisation models in order to further assess when the difference in power spectrum, which results from using this full treatment of the light cone, is most significant. 







§ ACKNOWLEDGMENTS

EC acknowledges the support of a Royal Astronomical Society Research Fellowship and studentship funding from the Science and Technology Facilities Council in the preliminary stages of this project. In later stages EC acknowledges the Royal Society Dorothy Hodgkin Fellowship.
MGS acknowledges support from the South African Square Kilometre Array Project and National Research Foundation (Grant No. 84156).

The authors would like to thank Garrelt Mellema, Hannes Jensen and Kanan Datta for useful discussions.

mnras  

