© The Authors 2023

1. Introduction

The Cosmic Dawn (CD) is the period at the beginning of the Epoch of Reionisation (EoR) when the first stars formed. A more quantitative definition, born from the study of the 21-cm signal emitted by the intergalactic medium (IGM) during the EoR, is to say that the CD corresponds to the period when the fluctuations of the 21-cm signal were dominated not by fluctuations of the ionisation or density fields, but rather by fluctuations of the spin temperature of hydrogen, which were in turn regulated by fluctuations of the gas kinetic temperature and of the strength of the Wouthuysen-Field coupling by Lyman-α photons (Wouthuysen 1952; Field 1958). In practice, for the more standard models, the CD corresponds to an averaged ionisation fraction of the IGM of less than a few percent. The corresponding redshift range is model-dependent: the CD occurs earlier if low-mass halos are able to form stars efficiently and later if not. If, as can be expected, atomic cooling halos (i.e. halos able to cool below their virial temperature through atomic cooling only, in other words halos with a mass ≳10⁸ M_⊙) do form stars efficiently, the CD occurred at z ≳ 15, meaning that the signal must be observed at frequencies ≲90 MHz. In this regime, the signal is seen in absorption against the Cosmic Microwave Background (CMB) and exhibits brightness temperature fluctuations on large scales with an amplitude up to several tens of millikelvins.

So-called global experiments attempt to detect the signal averaged over the whole sky through its frequency dependence only. Separation from foregrounds and possible instrumental effects is then an especially challenging task considering the limited leverage offered by the available information. Consequently, at this stage, such experiments report incompatible results: the EDGES experiment claims a detection around 78 MHz (Bowman et al. 2018) while the SARAS 3 experiment excludes the same signal with a 95% confidence level (Singh et al. 2022). Interferometers, on the other hand, have the ability to measure angular fluctuations of the signal and thus, in principle, produce a full tomography of the signal. This, however, requires very high sensitivity and not even the Square Kilometer Array (SKA) will be able to produce tomographic images during the CD (see, e.g. Mellema et al. 2013; Koopmans et al. 2015). The 3D isotropic power spectrum benefits from a better signal-to-noise ratio while retaining more information than the global signal. Nevertheless, published upper limits on current instruments able to probe the CD (LWA, MWA, LOFAR, and AARTFAAC) are orders of magnitude above the expected level of the signal (Ewall-Wice et al. 2016; Gehlot et al. 2019, 2020; Eastwood et al. 2019; Yoshiura et al. 2021). In this work, we compare our modelled signals to the expected sensitivity of the New Extension in Nançay Upgrading LOFAR (NenuFAR), (Mertens et al. 2021) and of SKA.

From the upper limits or a detection of the signal, the parameters of astrophysical models can be inferred using either a classical Markov chain Monte Carlo (MCMC) approach (e.g. Greig & Mesinger 2015, 2017) or with methods involving some aspects of machine learning (Shimabukuro & Semelin 2017; Gillet et al. 2019; Schmit & Pritchard 2018; Jennings et al. 2019; Doussot et al. 2019; Cohen et al. 2020; Hortúa et al. 2020; Bevins et al. 2022; Zhao et al. 2022; Bye et al. 2022; Abdurashidova et al. 2022). In all cases, the modelling of the signal is a fundamental step of the inference process: either at each step of the MCMC approach or for building a learning sample in methods based on supervised learning. The maximum likelihood values and posterior distribution of the astrophysical parameters is affected by the approximations made in the modelling step. The modelling approaches fall into two groups: fast semi-numerical methods (Thomas et al. 2009; Santos et al. 2010; Mesinger et al. 2011; Visbal et al. 2012; Fialkov et al. 2014) and slower full radiative transfer simulations (e.g Mellema et al. 2006; Baek et al. 2010; Zahn et al. 2011; Semelin et al. 2017). The semi-numerical methods make a number of approximations while the simulations are limited by their resolution given the available computing power. During the CD, the two processes that shape the 21-cm signal are the Wouthuysen-Field coupling and the heating of the IGM by X-rays (see Furlanetto et al. 2006, for a review). In this work we focus on the modelling of the Wouthuysen-Field coupling, also called Lyman-α coupling.

The 21-cm differential brightness temperature is related to the hydrogen spin temperature by the following:

$$\begin{array}{cc}\hfill \delta T_b& =27x_{{\mathrm{H}}_{\mathrm{I}}}(1+\delta )\left(\frac{T_{\mathrm{s}}-T_{\mathrm{cmb}}(z)}{T_{\mathrm{s}}}\right){(1+\frac{1}{H(z)}\frac{\mathrm{d}{v}_{||}}{\mathrm{d}{r}_{||}})}^{-1}\hfill \\ \hfill & \times {\left(\frac{1+z}{10}\right)}^{\frac{1}{2}}\left(\frac{{\mathrm{\Omega }}_b}{0.044}\frac{h}{0.7}\right){\left(\frac{{\mathrm{\Omega }}_m}{0.27}\right)}^{\frac{1}{2}}\mathrm{mK},\hfill \end{array}$$(1)

where x_HI is the local neutral fraction of hydrogen, δ is the overdensity of the gas, T_s is the local spin temperature of hydrogen, T_cmb is the CMB temperature, H(z) is the Hubble parameter, $\frac{\mathrm{d}{\mathit{v}}_{||}}{\mathrm{d}{r}_{||}}$ is the velocity gradient along the line of sight, z is the redshift, and the usual notation for cosmological parameters is used. The local value of the spin temperature is the result of three competing processes: thermalisation with the CMB; collisions with other particles that drive it to the local kinetic temperature of the gas; and pumping by Lyman-α photons (see Furlanetto et al. 2006, for details) that drives it to the colour temperature of the radiation around the Lyman-α wavelength. The large (∼10⁶) Gun-Peterson optical depth of the IGM during the CD allows the radiation spectrum around the Lyman-α line centre to reach thermodynamical equilibrium with the gas through the many scatterings, and thus the colour temperature of the radiation spectrum is almost identical to the gas kinetic temperature. As a result the spin temperature can be written as follows:

$$\begin{array}{c}\hfill T_S^{-1}=\frac{T_{\mathrm{cmb}}^{-1}+x_{\alpha }{T}_K^{-1}+x_{\mathrm{c}}{T}_K^{-1}}{1+x_{\alpha }+x_{\mathrm{c}}},\end{array}$$(2)

where x_c is the collisional coupling coefficient (negligible at z < 25) and x_a is the Lyman-α coupling coefficient defined by

$$\begin{array}{c}\hfill x_{\alpha }=\frac{4P_{\alpha }{T}_{21}}{27A_{10}{T}_{\mathrm{cmb}}}\end{array}$$(3)

and

$$\begin{array}{c}\hfill P_{\alpha }=4\pi {\displaystyle \int }{J}_{\nu }(\nu )\sigma (\nu )\mathrm{d}\nu ,\end{array}$$(4)

where T₂₁ is the 21-cm hyperfine transition excitation temperature, A₁₀ is the corresponding spontaneous emission coefficient, J_ν(ν) is the local angle-averaged specific intensity, and σ(ν) is the Lyman-α line cross-section. As we can see, J_ν(ν) is the main quantity, along with T_K, that can induce fluctuations in T_S. Thus modelling the Lyman-α coupling means modelling J_ν.

Handling the full radiative-transfer equation in an expanding non-homogeneous universe when line scattering and velocity gradients are involved is a daunting task. Even under simplifying assumptions (e.g. Loeb & Rybicki 1999), the equation is complex. It is easier to get a picture of the involved processes by calculating at the propagation of a single photon. A detailed discussion can be found in Semelin et al. (2007); here, we only summarise the main aspects. The physics of the transfer is encapsulated in a single quantity, the optical depth:

$$\begin{array}{c}\hfill \tau ={\displaystyle {\int }_0^L{\int }_{-\infty }^{\infty }{n}_{\mathrm{HI}}(l)P(u_{\Vert },T_{\mathrm{gas}}(l))\sigma (\nu (l)\times [1-\frac{v_{\Vert }(l)+u_{\Vert }}{c}])\mathrm{d}{u}_{\Vert }\mathrm{d}l,}\end{array}$$(5)

where n_HI is the local neutral hydrogen number density, σ is the Lyman-α scattering cross-section in the atom rest frame, ν(l) is the redshifting photon frequency in the global rest frame, v_∥(l) is the component of the local gas velocity in the global rest frame parallel to the direction of propagation, u_∥ is the parallel component of the atom’s velocity in the gas local rest frame, and P is the probability for the scattering atom to have a velocity u_∥ as dictated by the local thermal velocity distribution. There are two typical values that shape the Lyman-α transfer in the high-redshift universe. A photon redshifting from far in the blue wing of the line accumulates a τ ∼ 1 while still 10 cMpc away from the location where it would redshift into the core of the line (at z ∼ 10, in a homogeneous universe). Thus, scatterings in the wing of the line introduce a ∼10 cMpc diffusion scale compared to a free-streaming case (e.g. Chuzhoy & Zheng 2007; Semelin et al. 2007). Conversely, when the photon reaches the core of the line, it has a mean free path of less than 1 ckpc and the optical depth to redshift out of the line is of the order of 10⁶. Thus, wing scatterings determine where a photon reaches the core, and core scatterings dominate the overall budget of P_α but they are mainly local. This picture applies to a homogeneous expanding universe.

However, we can see how the local value of the gas density, ionisation state, temperature, and velocity are all involved in the computation of τ. Using local values instead of cosmic averages obviously modifies the computed τ value. In this work we quantify the impact on the computed P_α.

In the semi-numerical approach (Mesinger et al. 2011; Fialkov et al. 2014), J_ν(ν_α) at redshift z is computed using a series of fast Fourier transforms (FFTs) that implement the convolution of a propagation kernel with the emissivity fields at redshifts z_emit > z. In the original method, the kernel is built from a free-streaming approximation in a homogeneous medium: photons travel in a straight line until they cosmologically redshift into the line core. The homogeneity assumption makes the kernel spherically symmetric (with a Dirac radial dependence peaked at a radius whose value is determined by z, z_emit and the cosmology). A recent improvement (Reis et al. 2021) modifies the shape of the kernel to include the effect of scatterings in wings of the Lyman-α line, although still assuming a homogeneous medium. The SPINTER code presented in this work follows the same approach.

The Lyman-α coupling can also be computed through Monte Carlo radiative transfer simulations (see e.g. Semelin et al. 2007; Baek et al. 2009; Vonlanthen et al. 2011, using the LICORICE code). In this case P_α, the number of scatterings per atom per second is directly evaluated. However, computing an average of 10⁶ scatterings per photon and a sufficient number of photons (to control the sampling noise) reaching the core of the line in each resolution element at each desired output redshift is not computationally feasible yet. As a consequence photons are propagated (through an inhomogeneous medium) until they redshift into the core of the line and then a prescribed number of scatterings is tallied locally, considering that the spatial diffusion in the core of the line is negligible. Baek et al. (2009) give a prescription for this number, based on Monte Carlo simulations in controlled environments.

In this work we show that the prescription by Baek et al. (2009) for the number of scatterings in the core was incomplete and should be modified in the presence of gas velocities. More generally, we re-examine the impact of the different approximations made in the semi-numerical approach on the resulting 21-cm signal, using the (corrected) radiative transfer simulation as a proxy for the ground truth. We focus especially on the effect of gas velocities that seem to have the largest impact. In Sect. 2, using a Monte Carlo modelling including all core scatterings, we quantify the impact of gas velocities in an environment which was designed on purpose. In Sect. 3, we present the semi-numerical SPINTER code and remind the reader of the main features of the radiative transfer LICORICE code. In Sect. 4, we evaluate the impact of the various approximations in SPINTER in a realistic CD setup. Section 5 presents our conclusions.

2. The impact of gas bulk velocity on the Lyman-α coupling

2.1. The treatment of gas bulk velocities in existing methods

The radiative transfer in the Lyman-α line has been studied analytically in a cosmological context in several works (e.g. Rybicki & dell’Antonio 1994; Chen & Miralda-Escudé 2004; Chuzhoy & Shapiro 2006; Furlanetto & Pritchard 2006; Hirata 2006; Meiksin 2006). In most cases, the authors assume a homogeneous and isotropic universe, thus neglecting both the effect of Doppler shifts from the gas bulk velocity and the effect of fluctuations of the density and temperature. Many of these studies focus on the back-reaction from atomic recoil and spin exchanges on the Lyman-α spectrum near the centre of the line, using a Fokker-Planck formalism. These results should still apply if gas bulk velocities are accounted for by modifying the J_∞ (the angle-averaged specific intensity by photon number ‘far’ from the line centre, or alternatively, in the absence of a back-reaction) and using an effective local Hubble flow that includes the divergence of the peculiar velocity field. Loeb & Rybicki (1999) studied the transfer around a point source in a uniform Hubble flow, thus removing the homogeneity assumption on the emissivity field but not on the medium of propagation. They note that corrections from peculiar velocities may be necessary.

The semi-numerical approach (Furlanetto 2006; Santos et al. 2008; Mesinger et al. 2011; Fialkov et al. 2014) considers the actual, non-homogeneous emissivity field from cosmological sources, but it still estimates the effects of propagation through a homogeneous and isotropic universe. As such, they computed a local J_∞ that does not account for gas bulk velocities. Full Monte-Carlo radiative transfer simulations (Semelin et al. 2007; Baek et al. 2009 and subsequent works) do include Doppler shifts from gas bulk velocities. However, when computing the Lyman-α intensity in a cosmological box, to limit the computational cost, the radiative transfer is halted when photons reach the core of the line and a prescribed number of scatterings is tallied locally. Thus the effect of gas velocities are fully accounted for in the wings of the line (thus affecting where the photon subsequently reaches the core), but they should also enter through the prescription of the number of scatterings in the core of the line. This last contribution had not actually been implemented in LICORICE until now.

2.2. Expected magnitude of the effect

To clarify the issue, we can distinguish two regimes where gas velocities have an impact on the radiative transfer in the Lyman-α line: in the blue wing of the line and in the core of the line.

2.2.1. Effect in the wing

In the blue wing, where the medium is relatively transparent, Doppler shifts arising from velocity gradients along the path of the photons add their own contribution to the cosmological redshifting, modifying the time and location where the photon eventually reaches the core of the line in the local rest frame of the gas. The ratio of the Doppler to cosmological frequency shifts is, at least locally, equal to

$$\begin{array}{c}\hfill r=\frac{\mathrm{d}{\mathit{v}}_{||}}{\mathrm{d}l}\times \frac{1}{H(t)},\end{array}$$(6)

where H(t) is the Hubble parameter. A constant ratio r along the radial path of a photon emitted from a point source would cause the radius – where the photon reaches the core of the line – to shrink (or expand if r is negative) by a factor 1 − r (for small r values). Indeed, at a distance L from the source, the accumulated velocity difference would be H(t)rL, corresponding to a Doppler shift $\mathrm{\Delta }\nu =-\nu \frac{H(t)rL}{c}$ to be added to the cosmological shift $\mathrm{\Delta }\nu =-\nu \frac{H(t)L}{c}$. Since the frequency at the centre of the core is constant, the (1 + r) factor generated by the gas velocity should be compensated for by a (1 − r) factor applied to the distance from the source L. Photons that would reach the core of the line when crossing a surface element dS at a distance L from the source do so at a distance L(1 − r). The corresponding surface element is modified by a factor (1 − r)², changing the local flux by ∼1 + 2r (if r is small). In the same way, the local photon number density is modified by ∼1 + 3r.

2.2.2. Effect in the core

In the core of the line, where the medium is extremely opaque, the average number of scatterings before the photon is finally shifted to the red wing of the line is determined by the frequency shift between two scatterings (and by the gas density). In a uniform expanding universe, the number of scatterings is, on average, equal to the well-known Gunn-Peterson optical depth ${\tau }_{\mathit{GP}}=\frac{3\mathrm{\Lambda }{\lambda }_{\alpha }^3{n}_{\mathrm{HI}}}{8\pi H(t)}$, where n_HI is the number density of neutral hydrogen, λ_α is the wavelength at the centre of the line, and Λ is the natural line width (Gunn & Peterson 1965). As we can see, the number of scatterings scales as H⁻¹, which has been confirmed by Monte Carlo simulations (Baek et al. 2009). This scaling is not directly apparent in the analytical solutions to the Fokker-Planck equation given, for example, in Furlanetto & Pritchard (2006): in Sect. 2.3 we show that the reason is that the H⁻¹ scaling is encapsulated in the value of the angle-averaged intensity far from the line centre J_∞. This is corroborated by the analytical expressions in Rybicki & dell’Antonio (1994) and Chugai (1980) that exhibit the H⁻¹ scaling. The correction from the back-reaction of the gas on the Lyman-α spectrum (that also depends on H(t)) comes on top of this primary scaling and can be computed independently. In a non-uniformly expanding medium, one can expect any expansion or contraction due to gas velocities to act in the same way as a local Hubble expansion, and thus the number of scatterings per photon to be determined by the effective local value of the Gunn-Peterson optical depth:

$$\begin{array}{c}\hfill {\tau }_{\mathit{GP}}^{\mathrm{loc}}={\tau }_{\mathit{GP}}\frac{H}{H+\mathrm{div}(\mathit{v})/3}.\end{array}$$(7)

This ansatz is tested in Sect. 2.3. We note that, while ${\tau }_{\mathit{GP}}^{\mathrm{loc}}$ is the total number of scatterings, more than 99.9% of these occur within a few thermal widths from the core in frequency, where the mean free path is very short, and thus tallying the total number at the location where the core is reached is a very accurate approximation.

2.2.3. Expected magnitude of the velocity gradients

We now estimate the expected relative amplitude of the velocity gradient and Hubble parameter during the CD. From Newtonian perturbation theory, we know that, in comoving coordinates, the overdensity δ is related to the comoving velocity by the continuity equation: $\frac{\partial \delta }{\partial t}+\mathrm{\nabla }.\mathit{v}=0$. During the matter-dominated era, δ is proportional to the expansion factor a. Thus $H(t)=\frac{\partial \delta }{\partial t}\frac{1}{\delta }$. We can estimate the typical amplitude of density fluctuations at redshift 10 at the scale of the simulation resolution L ∼ 1 cMpc from structure formation theory as σ(L)∼0.35 (where σ(L) is the variance of the density field smoothed on scale L). Then $\frac{\mathrm{\nabla }.\mathrm{v}}{H(t)}\sim 0.35$. We checked that this is indeed the typical value that we find in the LICORICE simulations used in this work. We note that the relative effect on the number of scatterings per photon in the core of the line is one-third of this value, and that, in the wing, 0.35 can only be an upper limit for the ratio r, in particular configurations.

Moreover, we should mention that (i) the velocity effect in the wing of the line typically occurs on 10 − 100 cMpc scales, where σ(L) is a few percent at most; (ii) the effect in the core of the line, which does indeed occur on small scales, is then be smoothed on the scale of the instrument resolution, typically 10 cMpc for the SKA, where σ(L)∼0.05 − 0.10 at those redshifts; and (iii) σ(L) decreases as redshift increases so the effect would be smaller at redshift 20. Nevertheless, this rough estimate hints at a non-negligible contribution.

2.3. Validating the $\frac{1}{H+\mathrm{div}(\mathit{v})/3}$ scaling ansatz

As no complete analytical approach exists to describe the radiative transfer in the Lyman-α line in a non-homogeneously expanding medium, we turn to a Monte Carlo simulation to validate our ansatz. We use a simplified version of LICORICE, where the spatial dependence of physical fields is prescribed with analytical formulas, thus removing the need for grids and the difficulty of defining a spatial resolution. The Monte Carlo code computes the optical depth along the path of the photons taking both the Hubble expansion and the Doppler shifts from the gas velocity field into account. Thermal motion of atoms are included, scatterings are assumed isotropic in the rest frame of the atom, and atomic recoil is computed (but not the back-reaction from spin exchange). More details are given in Semelin et al. (2007) and Baek et al. (2009). To evaluate the effect of the gas velocity gradients on the Wouthuysen-Field coupling, we define a geometrically simple situation that nullifies other potential sources of fluctuations.

We consider a cosmological volume with uniform gas density and temperature (computed from a pure adiabatic cosmological evolution). The photons are emitted with a position (x_ini, 0, 0) with x_ini uniformly sampled between −130 and −70 cMpc, and an initial direction of propagation along the x axis. The initial frequency is chosen such that the photons redshift into the Lyman-α line at z = 10 after having travelled 100 cMpc. The redshift at emission is sampled in a (narrow) range, such that the photon frequency at the desired output redshift (z = 10 in our case) falls in a range a few hundreds of thermal widths around ν_α. To reduce the computing time, we ran the test at one-tenth of the actual gas density (thus reducing the number of scatterings for each photon to ∼0.8 × 10⁵). We have checked that in the absence of velocity fields and atomic recoil, this setup produces a local photon number density at z = 10 (and thus angle-averaged specific intensity) for any x coordinate between −20 and 20 cMpc that is spectrally flat in a range of 100 thermal widths centred on ν_α.

Then we consider several possible velocity fields:

• Case 1: no peculiar velocities

• Case 2:

  1. if x ∈ [ − 3, 3]       v = 3yH  u_y

  2. if x ∈ [ − ∞, −3] v = 0

  3. if x ∈ [3, ∞]        v = 0

• Case 3:

  1. if x ∈ [ − 3, 3]       v = yH  u_y + zH  u_z

  2. if x ∈ [ − ∞, −3] v = 0

  3. if x ∈ [3, ∞]        v = 0

• Case 4:

  1. if x ∈ [ − 3, 3]       v = (x + 3)H  u_x

  2. if x ∈ [ − ∞, −3] v = 0

  3. if x ∈ [3, ∞]        v = 6H  u_x.

In the above formulas, H is the Hubble parameter at redshift z = 10, and the coordinates are in units of comoving megaparsec. In cases 2 (resp. 3), the divergence of the velocity field in the slab x ∈ [ − 3.,3] equals (resp. equals two-thirds of) the contribution from the Hubble flow (that is 3H), and is 0 elsewhere. In these cases, where the velocities are perpendicular to the initial direction of propagation, there should be little effect from the free-streaming regime, but a strong effect from the core scatterings. In case 4, effects from the free-streaming regime and from the core scattering should combine. We note that the non-zero uniform velocity at x > 3 was chosen only to ensure continuity at x = 3. Since the gradient is zero in that region, we do not expect any net effect on J_ν. The velocities considered here are obviously larger than expected in a typical cosmological situation (especially being coherent on such large scales), but the goal here is simply to make the effect more visible to validate the ansatz.

In Fig. 1 we show the normalised, angle-averaged specific intensity J_ν(x, ν) at z = 10 in cases 1 (no velocities) and 3 (a slab with velocities perpendicular to the initial direction of propagation). The spectral number density of photons was actually computed, but it differs from J_ν(x, ν) only by a factor $\frac{4\pi }{c}$. In the left panel (no velocities), the effect of the back-reaction from atomic recoil is clearly visible as a depletion of the spectrum around ν_α (it is important to keep in mind that we are operating at one-tenth of the cosmological gas density, so the feature is narrower than at the nominal density). On the right, the slab between −3 and 3 cMpc is expanding perpendicularly to the x direction, creating an effective expansion rate 1.66 times as large as in the rest of the universe. We can check that, far from the line centre, J_ν(x, ν) is depleted by a factor of approximately 0.6 (see below for a quantitative view), confirming the effect of velocities on the J_∞ of the Fokker-Planck theory. The magnitude of the depletion is consistent with the scaling ansatz, as in case 3, in the centre slab, $\frac{1}{H+\mathrm{div}(\mathit{v})/3}=\frac{3}{5H}$ compared to $\frac{1}{H}$ in case 1. The spectrum is further depleted near the centre of the line by the gas back-reaction.

Fig. 1. Angle-averaged specific intensity near the Lyman-α line centre as a function of the position and frequency in two idealised setups described in the main text as case 1 (left panel) and case 3 (right panel). The angle-averaged specific intensity is normalised to 1 far from the line centre and where there are no effects from velocities.

Figure 2 shows cuts of the previous maps along the frequency direction, that is normalised spectra. Spectra outside and inside the slab with non-zero velocities are plotted. The scaling by a factor of 0.6 resulting from the × 1.66 effective expansion rate is confirmed. Figure 3 shows cuts along the spatial direction, both in the wings of the line and in the core. In the wing cuts, cases 2 and 3 show a depletion by a factor of 0.5 and 0.6 in the non-zero velocity slab, in line with the ansatz and their effective expansion rate. The effect seems to be, as expected, sensitive to div(v) only and not to the specific topology of the velocity field. The depletion is similar in the core, showing that the additional impact of velocities on the amplitude of the back-reaction trough the Gunn-Peterson optical depth is small. The not-so-sharp transitions at the boundaries of the velocity slab are due to the scatterings in the wings of the line that create a diffusion in the location where photons reach the core of the line. Case 4 is slightly more complex to interpret. The velocity field, which is oriented along the free-streaming direction of propagation, induces an effect both in the free-streaming regime and on the core scatterings. The velocity gradient along the (main) direction of propagation in the non-zero velocity slab equals the Hubble parameter. Thus photons that reach the gas-rest-frame core in the [−3,3] slab would otherwise reach the core in a [ − 3, 9] slab. The same is true for reaching any narrow rest-frame range of frequencies. This is just from redshifting and Doppler effects and this alone would boost J_ν by a factor of 2. However, the diffusion regime effect applies another $\frac{3H}{3H+\mathrm{div}(\mathit{v})}=3/4$ correction, leading to the ×1.5 observed correction. This agreement is consistent with both a $1+\frac{\mathrm{d}{\mathit{v}}_{||}}{\mathrm{d}l}/H$ correction factor in the wing and the (H + div(v)/3)⁻¹ scaling for the number of scatterings in the core. We now need to evaluate the impact of the correction due to the velocity field in a realistic cosmological case, where the amplitude of the velocities are typically smaller than in this setup.

Fig. 2. Normalised spectra around the Lyman-α line in different regions of case 1 and 3 (see main text) with either zero velocities (full lines) or an expanding peculiar velocity field (dashed line). These correspond to vertical cuts in Fig. 1

Fig. 3. Spatial fluctuations of the angle-averaged specific intensity at the centre of the Lyman-α line and in the wings for all four cases described in the main text, characterised by different velocity fields in the shaded slab. These correspond to horizontal cuts in Fig. 1.

3. Numerical methods

To run the cosmological test, we used two different codes: LICORICE, a full radiative transfer code, and SPINTER, a fast code using FFTs that applies a kernel to the emissivity field that takes scatterings into account but assumes homogeneity and isotropy for the gas.

3.1. The full radiative transfer with LICORICE

The LICORICE code performs the full Monte Carlo 3D radiative transfer in the Lyman-α line. It is described in detail in Semelin et al. (2007), Baek et al. (2009) and Vonlanthen et al. (2011), and has been used in a number of subsequent papers to compute the 21-cm signal during the CD (e.g Semelin et al. 2017). Here, we give a few relevant features of the code and refer the reader to the above references for a complete description. The Lyman-α part of the code performs Monte Carlo ray-tracing on a uniform grid. Directions, frequencies, and target optical depths of photon packets are sampled from the relevant distributions. Source luminosities are implemented in such a way that they are not affected by the Monte Carlo sampling noise (i.e. photons are assigned to sources in a deterministic way, in proportion to their luminosity). Optical depths are computed along the path of the photons, taking the local density, ionisation state, and temperature of the gas into account. The effect of the local proper velocity field and of cosmological redshifting on the frequencies of the photon (in the gas rest frame), and thus on the value of the cross-section for scattering, are taken into account. Indeed, the full Lyman-α line profile, including the wings, is used. Cascades from higher Lyman lines are also included (see Vonlanthen et al. 2011). Although it is optional in the code, in typical 21-cm simulations, propagation is stopped at the location where photons enter the core of the line and a calibrated number of scatterings is assigned to the corresponding cell. This number is the average number of scatterings required to redshift through the line core at the local density (see Baek et al. 2009). Indeed, actually propagating all photons until they redshift out of the line would typically be a thousand times more expensive (in proportion to the number of computed scatterings), reaching the 10⁷ CPU hour range for a grid with moderate resolution. The back-reaction from the gas on the local Lyman-α spectrum, due to the atomic recoil and spin exchange, is computed using the method described in Hirata (2006).

As the LICORICE code results serve in this work as a proxy for the ground truth, it is important to mention the limitations of the code. The most obvious one is the Monte Carlo sampling noise. This is even more the case than in situations where the photon packets can deposit a fraction of their content in each cell along their path, such as for radiative transfer of ionising photons or X-rays. For Lyman-α transfer, the scatterings occur essentially in the core of the line where the diffusion length is of the order of 1 ckpc, and the contribution of each photon is concentrated in one cell (while wing scatterings are important to determine in which location a redshifting photon reaches the line core, their contribution to the total scattering budget is negligible). Consequently, the relative level of the noise scales as ${\left(\frac{N_{\mathrm{phot}}}{N_{\mathrm{cell}}}\right)}^{-\frac{1}{2}}$, where N_phot is the number of photon packets that reach the core of the line in the time interval over which we want to estimate the average of the Lyman-α coupling, and N_cell is the number of resolution elements. We see that reaching an average noise level of 1% on a 256³ grid already requires ∼1.6 × 10¹¹ photons. In 21-cm simulations of the CD, we usually accept larger levels of noise and average the coupling estimation over a time interval of the order of 20 Myr. The SPINTER code, on the other hand, has no Monte Carlo noise and yields an estimate of the instantaneous coupling. Thus, for this work, for the cosmological comparison tests, we used a ∼1.5 Myr averaging time and we ensured a 1–2% typical noise level¹. Reaching that level of noise for a single target redshift on a 256³ grid in a cosmological box (200 h⁻¹ Mpc size) requires around 3000 single-core hours².

3.2. The semi-numerical approach with the SPINTER code

The analytical formula for estimating the average Lyman-α intensity in a homogeneous and isotropic universe and neglecting wing scatterings is described in Furlanetto (2006). Santos et al. (2010) and Mesinger et al. (2011) implemented a generalised version that accounts for a non-homogeneous emissivity field. The method to obtain the intensity field at a target redshift comes down to using a series of FFTs to convolve the emissivity field at higher redshifts with a Dirac-like spherical kernel whose radius depends on the emission and target redshifts. Reis et al. (2021) improved on that by using a spherically symmetric kernel that accounts for wing scatterings in a homogeneous medium. SPINTER implements a similar approach. In the following sections, we describe a formal framework for using a kernel that includes wing scattering that, to our knowledge, has not been formulated analytically before and we provide some details on the SPINTER implementation.

3.2.1. A theoretical framework

In the following, all quantities are comoving unless stated otherwise. We introduce ϵ(r, t, ν) dV dν dt, the number of photons emitted isotropically in a volume dV around position r, between times t and t + dt and between frequencies ν and ν + dν (that is ϵ is the emissivity by number). Furthermore, n is the number density of photons such that n(r, t, ν) dV dν is the number of photons in a volume dV around position r with a frequency between ν and ν + dν, at time t. We can write the relation between the two quantities resulting from radiative transfer in a very general way as follows:

$$\begin{array}{c}\hfill \mathcal{D}n=ϵ.\end{array}$$(8)

Under fairly general conditions, 𝒟 is a linear operator (e.g. if two-photon processes are negligible). If a procedure to calculate the Green function of operator 𝒟 can be devised, then n can easily be computed for any field ϵ. We now do so, first in the cosmological free-streaming case and then introducing resonant scattering.

3.2.2. The free-streaming case

We first assume that no absorption or scattering occurs. Instead of the time variable t, we use the redshift z. They are related by $\mathrm{d}t=-\frac{1}{H(z)}\frac{\mathrm{d}z}{1+z}$. We consider a pulse-like source term δ(r)δ(z − z₀)δ(ν − ν₀) and the corresponding Green function G for operator 𝒟:

$$\begin{array}{c}\hfill \mathcal{D}G(\mathit{r},z,\nu )=\delta (\mathit{r})\delta (z-z_0)\delta (\nu -{\nu }_0).\end{array}$$(9)

Then from the general theory of Green functions, we know that for any source field ϵ, n can be computed as follows:

$$\begin{array}{c}\hfill n(\mathit{r},z,\nu )={\displaystyle \int }ϵ({\mathit{r}}_0,z_0,{\nu }_0)G(\mathit{r},{\mathit{r}}_0,z,z_0,\nu ,{\nu }_0)d^3{r}_0\mathrm{d}{\nu }_0\mathrm{d}{z}_0.\\ \hfill \end{array}$$(10)

From the physics of free-streaming radiative transfer in a uniformly expanding universe, we know that photons emitted at redshift z₀ and frequency ν₀ will, at time t and redshift z, have redshifted to frequency $\nu =\frac{1+z}{1+z_0}{\nu }_0$, and will be located on a sphere of radius $r_1(z,z_0)={\int }_z^{z_0}\frac{c}{H(z^{\prime })}\mathrm{d}{z}^{\prime }$ centred on the emission point. Thus the Green function is necessarily of the form

$$\begin{array}{c}\hfill G=A\delta \left(\Vert \mathit{r}-{\mathit{r}}_0\Vert -r_1(z,z_0)\right)\delta (\nu -\frac{1+z}{1+z_0}{\nu }_0).\end{array}$$(11)

To establish the expression of A, we can use the conservation of the number of photons. More precisely, the number of photons emitted before time t is equal to the total number of photons present in the universe at time t (because we do not have any absorption terms). That is,

$$\begin{array}{c}\hfill {\displaystyle {\int }_V{\int }_{\nu }{\int }_{-\infty }^tϵ(\mathit{r},t^{\prime },\nu )d^{3}r\mathrm{d}{t}^{\prime }\mathrm{d}\nu ={\int }_V{\int }_{\nu }n(\mathit{r},t,\nu )d^{3}r\mathrm{d}\nu .}\end{array}$$(12)

Injecting Eqs. (10) and (11), we obtain the expression for A and we can write the full expression of the Green function as follows:

$$\begin{array}{c}\hfill G=\frac{1}{H(z_0)(1+z_0)}\frac{1}{4\pi r_1^2}\delta (\Vert \mathit{r}-{\mathit{r}}_0\Vert -r_1)\delta (\nu -\frac{1+z}{1+z_0}{\nu }_0).\end{array}$$(13)

Then, we can use the Green function to write the number density of photons generated by a source function ϵ. Defining r′=r₀ − r and n as a vector with norm 1, and performing the integration on ν₀ and in the radial dimension of r′,

$$\begin{array}{c}\hfill n(\mathit{r},z,\nu )=-{\displaystyle \int }\frac{1}{H(z_0)(1+z)}\frac{1}{4\pi }ϵ(\mathit{r}+r_1\mathit{n},z_0,\frac{1+z_0}{1+z}\nu ))\mathrm{d}\mathrm{\Omega }\mathrm{d}{z}_0.\end{array}$$(14)

The physical, angle-averaged specific intensity J is related to the photon comoving number density by $J={(1+z)}^3\frac{c}{4\pi }n$. Then,

$$\begin{array}{c}\hfill J={(1+z)}^2\frac{c}{{(4\pi )}^2}{\displaystyle {\int }_z^{+\infty }{\int }_{\mathrm{\Omega }}ϵ(\mathit{r}+r_1\mathit{n},z_0,\frac{1+z_0}{1+z}\nu )\frac{1}{H(z_0)}\mathrm{d}\mathrm{\Omega }\mathrm{d}{z}_0.}\end{array}$$(15)

If the emissivity ϵ is considered homogeneous at a fixed redshift, we can perform the angular integration and recover the expression often used in analytical models:

$$\begin{array}{c}\hfill J(z,\nu )={(1+z)}^2\frac{c}{4\pi }{\displaystyle {\int }_z^{+\infty }ϵ(z_0,\frac{1+z_0}{1+z}\nu )\frac{1}{H(z_0)}\mathrm{d}{z}_0.}\end{array}$$(16)

3.2.3. Introducing resonant scattering

The goal here is to perform an approximate computation of the angle-averaged specific intensity at the centre of the Lyman-α line, J_α, taking scatterings in the wings of the line into account (but not the back-reaction from the interaction with hydrogen atoms in the core of the line, this is handled in post-treatment). We make twwo simplifying assumptions. First we consider that the frequency of photons is unchanged during scatterings in the cosmological frame. This means we ignore the effects of thermal velocities of atoms and proper gas bulk velocities. The validity of this assumption is evaluated in Sect. 4. As a result, the frequency part of the Green function remains unchanged, as frequency shifts are only due to cosmological redshifting. Our second assumption is to consider that the spatial part of the Green function is isotropic. This is true only if the medium around the pulse source is homogeneous. By extension, it assumes that the effect of the gas density fluctuations on J_α is small. This assumption is also subsequently checked.

Then, introducing the function F to isolate the frequency dependence, the Green function takes the following form:

$$\begin{array}{c}\hfill G=F(\Vert \mathit{r}-{\mathit{r}}_0\Vert ,z,z_0)\delta ({\nu }_{\alpha }-\frac{1+z}{1+z_0}{\nu }_0)\end{array}$$(17)

with the normalising condition

$$\begin{array}{c}\hfill {\displaystyle \int }F(r,z,z_0)4\pi r^2\mathrm{d}r={\left[H(z_0)(1+z_0)\right]}^{-1}.\end{array}$$(18)

The main innovation in Reis et al. (2021) and in SPINTER is to compute the function F(r, z, z₀) numerically and approximately with a Monte Carlo simulation of wing scatterings for a single isotropic source emitting N photons at redshift z₀ in a homogeneous medium. We note that F depends on the gas density in principle, in practice we take it to be the average density of the universe.

Then the photon number density at redshift z can be written as follows:

$$\begin{array}{c}\hfill n(\mathit{r},z,v_{\alpha })={\displaystyle \int }F(\Vert \mathit{r}-{\mathit{r}}_{\mathbf{0}}\Vert ,z,z_0)ϵ({\mathit{r}}_0,z_0,{\nu }_0)(-\frac{1+z_0}{1+z})d^3{r}_0\mathrm{d}{z}_0,\end{array}$$(19)

with ${\nu }_0=\frac{1+z_0}{1+z}{\nu }_{\alpha }$.

3.2.4. Implementation

The first step in SPINTER is to compute the Green functions. In theory, a different Green function should be computed for each Lyman line since upper lines also contribute to J_α through cascades. In practice, the cross-section of the lines above Lyman-α is smaller and wing scatterings do not occur often. As a consequence, we use the free-streaming approximation for the upper lines. For the Lyman-α line, we use the resonant scattering formalism. For a given output redshift z, the two Green functions were tabulated as functions of the radius r and the emission redshift z₀ using a simple Monte Carlo ray-tracing method, subject to the assumption described in Sect. 3.2.3. The evaluation is fast (i.e. does not require too many Monte Carlo photons) because we assume that the Green functions have a radial dependence only. The thickness of the radial bins was set by the spatial resolution of the desired outputs. The z₀ bins were matched to the radial bins assuming a straight-line propagation.

Subsequently the photon number density n was computed as a sum over all the redshift bins of convolutions of the emissivity field with a kernel computed using the Green functions. The kernel is the sum of the Green function contributions from all included Lyman lines and also accounts for the periodic boundary conditions (by summing the contributions of several replica of the sources, and thus several evaluations of the same Green function, if required by the box size and values of horizons of the Lyman lines). The convolutions were computed in Fourier space. From n, J_α and x_α were computed. The back-reaction can also be computed (Hirata 2006).

SPINTER is written in Fortran. Parts of SPINTER are parallelised for shared memory architectures using OpenMP: the initial Monte Carlo computation of the Green functions and the computation of the kernel to use in the convolutions. For the FFTs, we used the Intel MKL library. Computing the target field(s) at a single redshift on a 256³ grid requires 0.5 CPU hours (on 2015 Intel Xeon E7-8857 CPUs), and a few minutes using several cores. The difference with the 3000 hours required for a LICORICE run is striking, but it is important to bear in mind that the ∼1% noise level used in LICORICE for this comparison would not be called for in many applications.

To study the effect of different approximations, both SPINTER and LICORICE are used with different setups. These are sumarized in Table 1.

4. Evaluating the impact of approximations in SPINTER

4.1. Cosmological test setup

The numerical approaches employed in SPINTER and LICORICE are so different that there are some subtleties involved in comparing their results. Guided by our final goal to be able to robustly model the contribution of the Wouthuysen-Field coupling to the 21-cm signal during the CD, we chose to compare the 3D field of the coupling coefficient x_α at fixed redshift, before the back-reaction (which can be evaluated in post-treatment at a negligible CPU cost) in a cosmological simulation box. The various fields that determine x_a (source emissivity, neutral gas density, temperature, and velocity of baryonic matter) are provided by a high-resolution radiative hydrodynamics simulation, HIRRAH-21, that resolves halos down to ∼4 × 10⁹ M_⊙ in a 200 h⁻¹ cMpc simulation box (Doussot & Semelin 2022). The native resolution of those fields, 2048³, were down sampled to 256³ (reducing the Monte Carlo noise at 2048³ resolution would be prohibitively costly). A difficulty is that SPINTER makes an evaluation of the instantaneous x_α based only on the past history of the emissivity field; whereas, LICORICE, due to the nature of Monte Carlo radiative transfer, evaluates an average of x_α over a redshift interval and takes the past evolution of all the other fields into account as they impacted the propagation of the photons from their emission point to the location where they redshift into the core of the Lyman-α line. To simplify the interpretation of the results, we used fields from HIRRAH-21 at an initial redshift z_ini and froze them until the redshift where we wanted to estimate x_α, z_final. The emissivity is assumed to be zero before z_ini. We typically chose z_final such that photons emitted just below Lyman-β at z_ini had enough time to redshift down to Lyman-α. Moreover, in the case of LICORICE, we only selected for actual-propagation photons that are expected to reach a Lyman line in the narrow δz = ±0.02 range around z_final. Thus we only needed to sample a narrow frequency range in the spectrum of the sources, whose boundary changes with the emitting redshift. In practice, photons emitted with a frequency around ${\nu }_{\mathrm{ini}}=\frac{1+z_{\mathrm{ini}}}{1+z_{\mathrm{final}}}{\nu }_{\alpha }$ in a range $\delta \nu =\frac{\delta z}{1+z_{\mathrm{final}}}{\nu }_{\mathrm{ini}}$ do not necessarily reach the local rest frame ν_α in the target redshift range due to Doppler shifts from the local gas velocities. We still propagated them until they reached ν_α and counted them, considering that they replace photons that would have reached the local ν_α in the target redshift range by having been emitted with a frequency slightly outside of the initial frequency range (we used a flat source spectrum). We were then able to bring the Monte Carlo noise level to just a few percent, even though the output was averaged only over a δz = ±0.02 interval.

We use z_ini = 12.2, which is a rather low value, to study the CD. Indeed, at this redshift the volume-averaged x_α is ∼1.4 in the HIRRAH-21 simulation. The reason for this late coupling is the rather high value of the minimum resolved-halo mass: 4.×10⁹ M_⊙. HIRRAH-21 is a full radiative-hydrodynamics simulation and reaching this mass resolution in a 200 h⁻¹cMpc box is already a challenge. At the same time, x_α ∼ 1 is the regime where fluctuations in the coupling are most likely to dominate the brightness temperature fluctuations. We believe that the effects we exhibit would be at least qualitatively similar for models where the studied regime occurs at higher redshift.

4.2. Comparing SPINTER and LICORICE Lyman-α coupling maps

We shall assume that the signal is seen in absorption at a redshift where the kinetic temperature of the neutral gas, T_K, is substantially smaller than the CMB temperature, x_α is of the order of ∼1, and the collisional coupling is negligible. Then we can simplify the computation of the spin temperature to $T_S^{-1}\sim \frac{x_{\alpha }{T}_K^{-1}}{1+x_{\alpha }}$. Since we then also have $\delta T_b\propto T_S^{-1}$, we can find that $\delta T_b\propto \frac{x_{\alpha }}{1+x_{\alpha }}$. Thus studying $\frac{x_{\alpha }}{1+x_{\alpha }}$ gives us a good first idea of the impact of the different ways of modelling x_α on the 21-cm brightness temperature. In Fig. 4 we show $\frac{x_{\alpha }}{(1+x_{\alpha })}$ maps (slices with single cell thickness) for z_ini = 12.2, for various modelling choices.

Fig. 4. Maps of the coupling coefficient xα (no back-reaction) at z = 12.2 in our test setup (see main text) for six different associations of code and modelling methods. The bottom left panel tick labels are in comoving megaparsec (the box size is 200 h−1cMpc).

The first (expected) result is that the map produced with SPINTER, including the effects of wing scatterings, is more contrasted than the map for SPINTER when the wing scatterings are ignored. Indeed, it has been shown before that the wing back-scatterings create a steeper radial profile of the Lyman-α flux around a point source (Chuzhoy & Zheng 2007; Semelin et al. 2007). This larger contrast is also found in Reis et al. (2021). The second result is that there is very little difference, apart from Monte Carlo sampling noise, between the map produced with SPINTER with wing scatterings and the one produced with LICORICE when the fluctuations of the HI number density and temperature, but not of the velocity, are included. This seems to indicate that using homogeneous density and temperature fields in SPINTER is an acceptable approximation. We verify this using the power spectrum.

However, we do see substantial changes in the maps when the effect of the gas velocities on the propagation of Lyman-α photons is included. To better understand those changes, we have separated two contributions for the effect of velocities: the effect in the wings that changes the location where the photons reach the core (which depends on the gradient of the velocity along the direction of propagation) and the effect in the core that changes the number of scatterings for each photon before redshifting out of the line (which depends on div(v)). The main effect in the wings is to weaken the coupling close to the sources: there, a large fraction of the local photons comes from the neighbouring source, and since the gas velocity field is converging towards the source the Doppler effect creates a blueshift that counterbalances the cosmological redshifting. Thus photons have to travel farther from the source to reach the gas rest-frame Lyman-α frequency. The second visible effect of velocity in the wings is to create rather small-scale fluctuations in the voids, where the coupling is the weakest. In the core of the line and near the sources, velocities have the opposite effect as the one they have in the wings: the negative velocity divergence increases the number of scatterings and thus the coupling intensity. The effect in the core also creates fluctuations in the voids. The visible consequences of the total velocity effect (wings and core) is i) an overall small decrease of the average coupling (a few percent), ii) a weakened coupling close to the sources, and iii) small-scale fluctuations in the voids. To gain a better grasp on the fluctuations in the voids, we now analyse the correlation to the density field.

4.3. Correlation of the velocity-induced fluctuations of x_α to the density field

We show in Sect. 2 that we expect the velocity gradients (along the line of propagation or acting through the divergence) to have an impact on the Lyman-α coupling. In the linear regime, the growth of the density field is proportional to the divergence of the velocity. So, we can expect the fluctuations in x_α caused by gas velocities to correlate with the density field, at least in some regions. We explore this possible correlation in Fig. 5 where we plotted 2D histograms of the number of grid cells in a given bin of x_α and overdensity δ, with different contributions from the velocities. The histograms are shown separately for regions more than 6 cMpc away from any source and for regions closer than 1.2 cMpc to a source. This split is more relevant than a split based on the density: as we see in the plots, overdense and underdense regions are found both near and far from the sources.

Fig. 5. 2D histograms of the number of grid cells as a function of overdensity and log10(xα). The upper row was computed including only cells that are located less than 1.2 cMpc from a source. The low row was computed including only cells that are located more than 6 cMpc away from any source. The different column corresponds to different modelling methods (the same as in Fig. 4).

The modification to x_α caused by including the velocity effect on core scatterings is entirely determined by the local value of the divergence of the velocity (see Eq. (7)). Overdense (underdense) regions typically have negative (positive) velocity divergence that should result in increased (decreased) x_α. This is exactly what we observe in the third column of Fig. 5: an increased positive correlation between x_α and δ is observed compared to the case without velocities (first column). This positive correlation exists both near and far from the sources.

In the case where the effect of velocities is included only for the propagation in the wings, the effect depends on the gradient of the velocity along the line of sight. Thus the local effect is direction dependent. Near the sources, we can expect a dominant contribution from photons travelling radially from the closest source and thus, as stated in Sect. 2, a net decrease of x_α, compared to a case without velocities. Typical spherical collapse predicts increasing velocity gradients towards the higher-density centre, so in our case a stronger decrease of x_α. We do observe this anti-correlation between x_α and δ in Fig. 5. What is less anticipated is that the anti-correlation persists far from the sources, where the radiation field is more isotropic; in the case of a perfectly isotropic radiation field, we would expect a zero net effect as photons travelling in opposite directions would experience opposite effects from the velocity field. This seems to indicate that a correlation between the velocity field (and thus the density structures, i.e. pancakes and filaments) and the main direction of propagation of photon is still effective far from the sources. The last column of Fig. 5 indicates that the correlation and anti-correlation induced by the effect of velocities on the transfer in the core and in the wings do cancel out to a large extent when the two effects are combined.

4.4. Impact on the power spectrum of the Lyman-α coupling

Figure 6 shows the 3D isotropic power spectrum of $\frac{x_{\alpha }}{1+x_{\alpha }}$ in the same six cases as Fig. 4. The same effects that we identified on the maps are present in the power spectra. Quantitatively, including wing scatterings in SPINTER boosts the power on all scales by a factor of 2 or more. LICORICE without gas velocities gives results nearly identical to SPINTER with wing scatterings, except for a small boost at small scales either due to residual Monte Carlo noise or some limited self-shielding effects in high-density regions (Semelin et al. 2007). Including the effect of gas velocities only for the number of scatterings in the core boosts the power on all scales (from a factor 1.5 on large scales up to a factor 3 on small scales), while only including the effect of gas velocities on the propagation in the wings decreases the power on large scales (by a factor ∼3) and increases it on small scales (by a similar factor). The total effect of velocities, combining the two contributions, is to decrease the power on large scales and increase it on scales corresponding to wavenumbers larger than 1 h cMpc⁻¹. We note that the magnitude of the variations between the different modellings may depend on the history and morphology of the Lyman-band emissivity field (as we check by looking at the same quantities at different redshifts). What matters is that these variations exist and cannot be ignored, at least in some specific regimes.

Fig. 6. 3D isotropic power spectra of $\frac{x_{\alpha }}{1+x_{\alpha }}$ in a realistic cosmological setting at z = 12.2 for the different code – modelling method combinations.

4.5. Impact on the 21-cm brightness temperature power spectrum

While the previous study of the impact of gas velocities on x_α is interesting because this quantity is close enough to the physics of radiative transfer such that we can readily interpret the results, it is not sufficient to estimate whether the full modelling of gas velocities would change how we infer astrophysical knowledge from 21-cm power spectrum observations during the CD. Indeed, in addition to x_α fluctuations, neutral hydrogen density fluctuations and gas kinetic temperature fluctuations determine the brightness temperature power spectrum. Moreover, these fluctuations are clearly not uncorrelated and have varying relative contributions depending on the redshift and on the astrophysical model, as parameterised for example by f_X which quantifies the intensity of heating by X-rays (see e.g. Semelin et al. 2017, for a full definition).

In Fig. 7 we present the 3D isotropic power spectrum of δT_b computed from the HIRRAH-21 fields at z_ini = 12.2 for f_X = 1 and f_X = 0. The X-ray contribution is a mix of sources with a soft spectrum (such as active galactic nuclei) and a hard spectrum (such as X-ray binaries) in equal contribution (a parameter r_H/S = 0.5 as defined in Semelin et al. 2017). The HIRRAH-21 simulation was run with f_X = 1. The brightness temperature for f_X = 0 can be evaluated in post-treatment by recomputing the local kinetic temperature of the neutral gas assuming an adiabatic evolution from a homogeneous universe at z ∼ 130, the redshift of thermal decoupling between the gas and the CMB. On large scales we observe an effect similar to what we found for the $\frac{x_{\alpha }}{1+x_{\alpha }}$ power spectrum: a boost when including wing scatterings in SPINTER and a depletion when including velocities in LICORICE. What happens on smaller scale seems to depend on the presence of X-ray heating. When no X-ray heating is present (f_X = 0), the δT_b power spectrum shows similar reactions to the different approximations as the Lyman-α coupling power spectrum. When X-ray heating is present (f_X = 1), the δT_b power spectrum shows little sensitivity on small scales to the x_α modelling.

Fig. 7. 3D isotropic 21-cm signal power spectrum for a realistic cosmological setting at z = 12.2 (see main text) plotted for different modelling methods and codes. In the left panel, the IGM has been heated with a fiducial amount of X-rays (fX = 1); in the right panel, the IGM has been subject to adiabatic cooling and heating only (fX = 0). The expected SKA thermal noise level for 1000 h of observation (see the main text for further details) has also been plotted to quantify how the Lyman-α modelling could bias parameter inference performed on SKA observations.

In the simulation with f_X = 1, the neutral IGM has been heated to an average temperature of ∼7 K, a 2 K increase from the ∼5 K in the f_X = 0 case. Due to the presence of soft X-ray with a limited mean free path in the neutral IGM, we can expect the heating close to the sources to be even larger, locally initiating the heating transition towards a signal in emission. Then, the power spectrum of the 21-cm signal may be dominated on small scales by the contribution of these heated bubbles, and thus insensitive, on these scales, to the Lyman-α fluctuations. It is likely that the scale below which the heating fluctuations become dominant depends on the relative contribution of hard and soft X-rays as, typically, a harder spectrum results in heating on larger scales.

In Fig. 7, we also plotted, following the methodology of McQuinn et al. (2006), the expected thermal noise level for 1000h of observation with SKA at z = 12 (10 MHz bandwidth, width of k bins equal to the k value at the centre of the bin, T_sys = 100 + 300( ν / 150 MHz )^−2.55  K, and 30-λ low k cutoff). Sample variance is not included in the plot. For reference, in this case, it equals 10% of the signal power spectrum at k = 0.03 h cMpc⁻¹ and already decreases to 1.6% at k = 0.1 h cMpc⁻¹. As we can see, SKA should be able to easily discriminate between the various modelling methods. This means that the inference of astrophysical parameters from SKA observations are biased if an approximate Lyman-α coupling modelling is used.

4.6. Observability with NenuFAR and SKA

The HIRRAH-21 simulation, despite its high resolution for a fully coupled radiative transfer simulation, does not resolve halos with masses less than ∼4 × 10⁹ M_⊙, whereas we know that halos with masses down to 10⁸ M_⊙ should efficiently form stars from hydrogen atomic cooling. This results in a rather late onset of the CD and of reionisation. A complementary assessment of the potential impact of Lyman-α coupling modelling should focus on higher redshift scenarios to determine how Lyman-α modelling affects the interpretation of observations. We present such an assessment here making the assumption that the amplitude of the effects of including wing scatterings versus including velocities remains of the same order at higher redshift. Indeed, we show only the effect of wing scatterings because running a full radiative transfer of the Lyman-α with an averaging of the coupling over a sufficiently narrow redshift interval is very difficult while considering a full cosmological evolution (it is important to remember that in the previous cosmological test, the fields were frozen at z = 12.2). If we were to use the redshift interval that we typically use in LICORICE simulations for averaging the coupling, we would not be able to easily distinguish between the effect of velocities and the effect of averaging over a large interval.

The models plotted in Fig. 8 were computed with a modified version of LICORICE that takes into account halos below the simulation resolution using the conditional mass function formalism in a way similar to 21cmFAST (see Gillet et al. 2021, for an alternative approach to including unresolved star formation). Resolved halos were treated as before. An unresolved collapsed fraction (including the effect of Poisson noise on the conditional mass function) was computed for each particle and it contributes to the star formation. Full details about the implementation will be given in Meriot & Semelin, in prep. The plotted models use a 200 h⁻¹cMpc box and 256³ particles. While they only resolve halos with masses larger than 2.×10¹² M_⊙, the conditional mass function treatment allows for star formation in unresolved halos with masses down to either 4 × 10⁷ M_⊙ or 10⁸ M_⊙. A weak X-ray contribution was used, f_X = 0.1, to maximise the absorption signal. It is important to note, however, that using f_X = 1 would not decrease the signal much as it is still insufficient to start the heating transition at the redshift of interest (z = 16.5). The ’exotic’ model includes an additional homogeneous radio background following Fialkov & Barkana (2019), with an intensity parameter A_r = 9.6 and spectral index β = −2.6 (resulting in an additional background 10 times stronger than the CMB at the redshift of interest), such that the corresponding global signal has an amplitude corresponding to the claimed EDGES detection by Bowman et al. (2018; however, see also Singh et al. 2022). All models use an f_α = 2 (see Semelin et al. 2017, for a definition), a phenomenological parameter that allows one to vary the average intensity of the Lyman-α coupling. The power spectra are plotted at z = 16.5.

Fig. 8. 3D isotropic power spectra of the 21-cm signal at z = 16.5 for three different models (see the main text for further details) plotted for two different ways of modelling the Lyman-α coupling (both without the contribution of gas velocities) and compared to the expected thermal noise of observations on the NenuFAR and SKA radio interferometers. The dotted black line shows the contribution from sample variance in the case of the exotic signal.

The redshift was selected as the lowest possible value for an observation with the NenuFAR radio interferometer while averaging over a 10 MHz bandwidth. Indeed, NenuFAR³ (Zarka et al. 2012, 2020) has a 10-85 MHz operating bandwidth. In 2019, the Cosmic Dawn Key Science Program⁴ (ES01, P.I. L. Koopmans, B. Semelin, F. Mertens), hereafter CD KSP, started on NenuFAR. It accumulated > 1300 h of single field observation by the end of 2022 (Mertens et al. 2021). The expected thermal noise of NenuFAR is plotted together with the models. With NenuFAR still being under deployment, the CD KSP observations have been acquired with varying (expanding) configurations: 475 h with a 56-station core configuration and 460 h with an 80-station core configuration; we also expect ∼140 more hours with the 80-station core configuration and 150 h with the full 96-station core configuration by the end of 2023. The CD KSP noise level plotted in the figure reflects this mix. It uses the same Δk = k bin width and 10 MHz bandwidth as for the SKA, but a 10-λ low k cutoff (due to the different station configuration). It does not include contributions from the k_∥ = 0 modes, following Mertens et al. (2020). The system temperature, T_sys = T_inst + T_sky, uses the same T_sky as for SKA but T_inst = 874 K, the value given at 80 MHz by the Nenupy software (Loh & Girard 2020). The NenuFAR – 5000 h noise level assumes 5000 h of observation on the full 96-station core configuration, and the SKA 1000 h uses the same parameters as for Fig. 7 but at redshift z = 16.5.

The exotic signal is detectable by NenuFAR CD KSP at wavenumbers k < 0.1 h cMpc⁻¹, at a level likely sufficient to discriminate between the two Lyman-α coupling modellings. To quantify the discrimination power, one would need to compare the power spectrum of the difference of the models to the thermal noise and sample variance. This would, however, not be enough to yield an estimate on the induced bias on the model parameters, which is the final issue. We do not, at this time, have a full framework to do this computation. The two more standard signals are only marginally detectable with NenuFAR in 5000 h at the largest scales and the robustness of this detection would be sensitive to the Lyman-α modelling. We note that the plotted ‘non-exotic’ models are still rather well suited for a detection (high f_α and low minimal halo mass for efficient star formation). They are, however, incompatible with the Bowman et al. (2018) claimed detection, since they have a sky-averaged absorption signal of only ∼40 mK at this redshift. Thus, NenuFAR can in principle put constraints on the more optimistic models in terms of signal strength, but an accurate determination of which models can be excluded and at what level would have to rely on an accurate Lyman-α coupling modelling.

Figure 8 also shows that a 1000-hour observation with SKA should be able to not only detect the signal for the plotted models on a large range of wavenumbers, but it should be sensitive enough to discriminate between the different Lyman-α coupling modellings. We note that a 100 h survey (the medium survey in Koopmans et al. 2015), with a ×10 higher thermal noise for the power spectrum, would already detect the signal in these favourable cases.

5. Conclusions

The modelling of the Lyman-α coupling of the hydrogen spin temperature to the gas kinetic temperature has a long history. While the theory has been established more than fifty years ago (Wouthuysen 1952; Field 1958), the necessity of computing the local spectrum around the Lyman-α line has led, as least initially, to drastic simplifications. Initially, the coupling was simply assumed to saturate very quickly and thus predictions during the CD, when the spatial fluctuations of this coupling are a dominant contribution to the 21-signal, were unreliable. Then, both semi-numerical methods and full radiative transfer codes were developed to compute the coupling and produce more accurate predictions for the 21-cm signal during the CD. Recently, Reis et al. (2021) have improved the modelling of semi-numerical codes to include the effect of wing scatterings.

In this work, we have performed a careful comparison of the results from the semi-numerical method and a full radiative transfer code. We find that, ignoring the role of gas peculiar velocity, the two methods agree to a good level, even though the semi-numerical method assumes a propagation in a homogeneous universe. However, we find that when the Doppler effect from gas velocities is included (as it is by default) in the full radiative transfer code, the resulting 21-cm signal power spectrum is modified by up to a factor of ∼2 at some scales and redshifts that depend on other model parameters. We have presented some theoretical estimates of the expected amplitude of the effect of velocities on the local Lyman-α flux. We have shown that the effect of velocities can be analysed by distinguishing two contributions. The first occurs while the photons are still in the wing of the Lyman-α line; it mainly changes the location where the photons redshift into the core. The second contribution occurs during the propagation in the core of the line; it changes the number of scatterings a photon undergoes before redshifting out of the core. Both effects tend to counter each other, but they do not necessarily balance out. Finally we have shown that various 21-cm signals resulting from different levels of Lyman-α modelling, with everything else being identical, can be distinguished with high significance by SKA; while the same is true for the NenuFAR CD KSP, but only for exotic models including an additional radio background.

At this stage, it is unclear how the effect of velocities could be included in a fast semi-numerical modelling. We have obtained very good results by training a neural work to produce a correction to be applied to the output of SPINTER to reproduce the output of LICORICE. However, this network was trained with data from a specific model at a specific redshift (using half of the data cube for training and half for testing). To be confident with using such a network, it would need to be trained with a variety of models and different redshifts. That would require running many LICORICE simulations to build the training sample, with stringent constraints on redshift interval averaging. That would be a challenge in terms of computing time. Furthermore, if such a training sample were available, a network could probably be trained, not to compute a correction to SPINTER, but probably as an emulator to LICORICE. It is not impossible that using a modified kernel in SPINTER, which would implement an average effective radial velocity profile (probably redshift dependent), could lead to an improved agreement. A theoretical framework remains to be formulated for such an approach to be more than just phenomenological.

Acknowledgments

Some of the numerical simulations in this work were performed using HPC resources from GENCI-CINES (grant 2018-A0050410557). LVEK acknowledges the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 884760, ‘CoDEX’). RB acknowledges the Israel Science Foundation (grant No. 2359/20), the Vera Rubin Presidential Chair in Astronomy and the Packard Foundation.

References

All Tables

All Figures

	Fig. 1. Angle-averaged specific intensity near the Lyman-α line centre as a function of the position and frequency in two idealised setups described in the main text as case 1 (left panel) and case 3 (right panel). The angle-averaged specific intensity is normalised to 1 far from the line centre and where there are no effects from velocities.
In the text

	Fig. 2. Normalised spectra around the Lyman-α line in different regions of case 1 and 3 (see main text) with either zero velocities (full lines) or an expanding peculiar velocity field (dashed line). These correspond to vertical cuts in Fig. 1
In the text

	Fig. 3. Spatial fluctuations of the angle-averaged specific intensity at the centre of the Lyman-α line and in the wings for all four cases described in the main text, characterised by different velocity fields in the shaded slab. These correspond to horizontal cuts in Fig. 1.
In the text

	Fig. 4. Maps of the coupling coefficient xα (no back-reaction) at z = 12.2 in our test setup (see main text) for six different associations of code and modelling methods. The bottom left panel tick labels are in comoving megaparsec (the box size is 200 h−1cMpc).
In the text

	Fig. 5. 2D histograms of the number of grid cells as a function of overdensity and log10(xα). The upper row was computed including only cells that are located less than 1.2 cMpc from a source. The low row was computed including only cells that are located more than 6 cMpc away from any source. The different column corresponds to different modelling methods (the same as in Fig. 4).
In the text

	Fig. 6. 3D isotropic power spectra of $\frac{x_{\alpha }}{1+x_{\alpha }}$ in a realistic cosmological setting at z = 12.2 for the different code – modelling method combinations.
In the text

Fig. 7. 3D isotropic 21-cm signal power spectrum for a realistic cosmological setting at z = 12.2 (see main text) plotted for different modelling methods and codes. In the left panel, the IGM has been heated with a fiducial amount of X-rays (fX = 1); in the right panel, the IGM has been subject to adiabatic cooling and heating only (fX = 0). The expected SKA thermal noise level for 1000 h of observation (see the main text for further details) has also been plotted to quantify how the Lyman-α modelling could bias parameter inference performed on SKA observations.

In the text

Fig. 8. 3D isotropic power spectra of the 21-cm signal at z = 16.5 for three different models (see the main text for further details) plotted for two different ways of modelling the Lyman-α coupling (both without the contribution of gas velocities) and compared to the expected thermal noise of observations on the NenuFAR and SKA radio interferometers. The dotted black line shows the contribution from sample variance in the case of the exotic signal.

In the text