Robust Velocity-induced Acoustic Oscillations at Cosmic Dawn

TL;DR

This paper introduces 21cmvFAST, a modified simulation tool that predicts robust velocity-induced acoustic oscillations in the 21-cm power spectrum during cosmic dawn, which upcoming observations could detect to confirm early universe acoustic phenomena.

Contribution

The paper demonstrates that velocity-induced acoustic oscillations are a distinct, analytically predictable feature in the 21-cm power spectrum, and introduces a simulation tool to study them.

Findings

VAOs are detectable at high significance by HERA.

VAOs have a unique shape distinguishable from density fluctuations.

The shape of VAOs can be analytically computed.

Abstract

The redshifted 21-cm line of hydrogen holds great potential for the study of cosmology, as it can probe otherwise unobservable cosmic epochs. In particular, measurements of the 21-cm power spectrum during cosmic dawn---the era when the first stars were formed---will provide us with a wealth of information about the astrophysics of stellar formation, as well as the origin of fluctuations in our Universe. In addition to their usually considered density fluctuations, dark matter and baryons possess large relative velocities with respect to each other, due to the baryon acoustic oscillations (BAOs) suffered by the latter, which suppress the formation of stars during cosmic dawn, leaving an imprint on 21-cm observables during this era. Here we present 21cmvFAST, a version of the publicly available code 21cmFAST modified to account for this effect. Previous work has shown that the inclusion of relative velocities produces an acoustic modulation on the large-scale 21-cm power spectrum during cosmic dawn. By comparing analytic calculations with simulations from 21cmvFAST, here we demonstrate that this modulation takes the form of robust velocity-induced acoustic oscillations (VAOs), during both the Lyman-$\alpha$ coupling era and the subsequent epoch of heating. The unique shape of these VAOs, which is determined by the acoustic physics that generated the relative velocities, can be analytically computed and is easily distinguishable from the usual density-sourced fluctuations. We find that, under a broad range of astrophysical assumptions, the VAOs are detectable at high significance by the upcoming HERA interferometer, which could therefore confirm the presence of acoustic oscillations at cosmic dawn.

Tables

Table 1. Table 1 : Projected signal-to-noise ratios for each combination of feedback and foreground assumptions, in the first table for the X-ray heating era ( z = 13 − 19 𝑧 13 19 z=13-19 ) and in the second table for the Lyman- α 𝛼 \alpha coupling era ( z = 23 − 27 𝑧 23 27 z=23-27 ). In all cases we assume three years of observation with HERA.

	Foregrounds
Feedback strength	Pessimistic	Moderate	Optimistic
High	2	4	14
Regular	10	15	62
Low	20	33	133

Table 2. Table 2 : Same as Table  1 , albeit with an X-ray luminosity of log 10 ⁡ ( L X / SFR ) = 40.8 subscript 10 subscript 𝐿 𝑋 SFR 40.8 \log_{10}(L_{X}/{\rm SFR})=40.8 , a factor of 6.3 larger than our fiducial case. This pushes all heating transitions Δ ​ z ∼ 2 similar-to Δ 𝑧 2 \Delta z\sim 2 earlier, so the EoH covers the range z = 15 − 21 𝑧 15 21 z=15-21 .

	Foregrounds
Feedback strength	Pessimistic	Moderate	Optimistic
High	2	4	23
Regular	7	11	47
Low	12	20	83

Equations

$${\delta }_c={\delta }_b=0$$

$${\mathbf{v}}_{\mathrm{cb}}(z)={\mathbf{v}}_{\mathrm{cb}}^i\frac{1+z}{1+z_{\mathrm{kin}}},$$

$${\mathbf{v}}_b$$

$${\mathbf{v}}_c$$

$${\delta }_c^{\prime }$$

$${\theta }_c^{\prime }$$

$${\delta }_b^{\prime }$$

$${\theta }_b^{\prime }$$

$${\delta }_T^{\prime }$$

$$c_s^2=\frac{k_B{T}_g}{{\mu }_b},$$

$$\frac{{\delta }_b}{{\delta }_m}=1+r_{\mathrm{LSS}}-k^2/k_F^2$$

$$M_F\equiv \frac{4\pi }{3}{\rho }_m^0{r}_F^3,$$

$$f_{\mathrm{gas}}(M,v_{\mathrm{cb}})=f_b^{(0)}{\left[1+(2^{a/3}-1){\left(\frac{M_F(v_{\mathrm{cb}})}{M}\right)}^a\right]}^{-3/a},$$

$$V_{\mathrm{cool}}(z,v_{\mathrm{cb}})={\left[{\left(V_{\mathrm{cool}}^{(0)}\right)}^2+{\alpha }^2{v}_{\mathrm{cb}}^2(z)\right]}^{1/2},$$

$$M_{\mathrm{cool}}(z,v_{\mathrm{cb}})=7\times 10^6{M}_{\odot }{\left(\frac{V_{\mathrm{cool}}(z,v_{\mathrm{cb}})}{10\mathrm{km}{\mathrm{s}}^{-1}}\right)}^3{\left(\frac{1+z}{20}\right)}^{-3/2},$$

$$M_{\mathrm{cool}}(z,v_{\mathrm{cb}},F_{\mathrm{LW}})=M_{\mathrm{cool}}(z,v_{\mathrm{cb}},0)\times \left[1+B(F_{\mathrm{LW}}{)}^{\beta }\right],$$

$$f_{\ast }(M)={\alpha }_{\ast }\log (M/M_{\mathrm{cool}})$$

$$P_m^{\mathrm{loc}}(k,z,v_{\mathrm{cb}})=P_{\zeta }(k){\int }_{-1}^1\frac{d\mu }{2}{\left|{\mathcal{T}}_m(k,z,v_{\mathrm{cb}},\mu )\right|}^2,$$

$${\sigma }_M^2(z,v_{\mathrm{cb}})=\int \frac{d^{3}k}{(2\pi {)}^3}|W(k{R}_M){|}^2{P}_m^{\mathrm{loc}}(k,z,v_{\mathrm{cb}}),$$

$$\nu (M,z,v_{\mathrm{cb}})={\left(\frac{{\delta }_{\mathrm{crit}}}{{\sigma }_M(z,v_{\mathrm{cb}})}\right)}^2,$$

$$\nu f(\nu )=\mathcal{A}(1+{\beta }_1{\nu }^{\prime -p}+{\beta }_G{G}_1)\sqrt{{\nu }^{\prime }}\exp (-{\nu }^{\prime }/2.0),$$

$$G_1=\exp \left[-\frac{(\log ({\nu }^{\prime })-d{)}^2}{2\times (0.6{)}^2}\right],$$

$$n(M,z,v_{\mathrm{cb}})=\frac{{\rho }_m^0}{M^2}\left(-2\frac{d\log {\sigma }_M}{d\log M}\right)\nu f(\nu ).$$

$$F_{\mathrm{coll}}(z,v_{\mathrm{cb}})={\int }_{M_{\mathrm{cool}}}^{\infty }dMMn(M)\frac{f_{\mathrm{gas}}}{{\rho }_b}{f}_{\ast },$$

$$F_{\mathrm{coll}}^{\mathrm{tot}}(\delta ,v_{\mathrm{cb}})=\frac{F_{\mathrm{coll}}(v_{\mathrm{cb}})}{\overline{F_{\mathrm{EPS}}}}{F}_{\mathrm{EPS}}(\delta ,v_{\mathrm{cb}}),$$

$$F_{\mathrm{EPS}}(\delta ,v_{\mathrm{cb}})=\mathrm{erfc}\left[\frac{{\delta }_{\mathrm{crit}}-\delta }{\sqrt{2[{\sigma }_{\mathrm{cool}}^2(v_{\mathrm{cb}})-{\sigma }_R^2]}}\right],$$

$$\frac{n_1}{n_0}=\frac{g_1}{g_0}{e}^{-T_{\ast }/T_s},$$

$$T_{21}=38\mathrm{mK}\left(1-\frac{T_{\gamma }}{T_s}\right){\left(\frac{1+z}{20}\right)}^{1/2}{x}_{\mathrm{HI}}(1+{\delta }_b)\frac{{\partial }_r{v}_r}{H(z)},$$

$${{\mathbf{v}}_{\mathrm{cb}}}^{(i)}(\mathbf{k})=\frac{i\mathbf{k}}{k}\sqrt{\frac{P_v(k)}{P_{\delta }(k,z_i)}}{\delta }^{(i)}(\mathbf{k}),$$

$$⟨\delta T_{21}(\mathbf{k})\delta T_{21}^{\ast }({\mathbf{k}}^{\prime })⟩=(2\pi {)}^3{\delta }_D(\mathbf{k}-{\mathbf{k}}^{\prime })P_{21}(\mathbf{k}),$$

$${\Delta }^2(k,z)=P_{21}(k,z)\frac{k^3}{2{\pi }^2},$$

$$⟨f(\mathbf{k})f^{\ast }({\mathbf{k}}^{\prime })⟩\approx (2\pi {)}^3{\delta }_D(\mathbf{k}-{\mathbf{k}}^{\prime })b_f^2{\Delta }_{v^2}^2\frac{k^3}{2{\pi }^2},$$

$${\delta }_{v^2}=\sqrt{\frac{3}{2}}\left(\frac{v_{\mathrm{cb}}^2}{v_{\mathrm{rms}}^2}-1\right),$$

Full text

I Introduction

The matter in our Universe clusters under its own gravitational pull, slowly growing the seeds of primordial fluctuations into the rich large-scale structure of our cosmos [1, 2, 3, 4]. Nonetheless, not all matter behaves equally. While the majority of matter is dark, cold, and collisionless, roughly a sixth of it is baryons, capable of interacting with photons through electromagnetic forces. These interactions give rise to the famous baryon acoustic oscillations (BAOs), which are imprinted onto the distribution of baryons thereafter and have been observed in the cosmic microwave background (CMB) [5], as well as at low-redshift surveys [6, 7, 8].

Interestingly, the lack of BAOs for the dark matter (DM) gives rise to a bulk relative velocity between it and the baryons [9]. This “streaming” velocity, with a root-mean-squared value of $v_{\rm rms}\approx 30$ km s*-1* at recombination, becomes supersonic thereafter due to the low baryonic temperature. The resulting supersonic motions between baryons and DM—which fluctuate on acoustic scales—have three suppressive effects on the formation of the first stars, as $(i)$ fewer haloes are formed [9, 10, 11], $(ii)$ the amount of gas in each halo available for cooling is smaller [12, 13, 14], and $(iii)$ it is more difficult for said gas to settle down into stars [15, 16, 17, 18, 19, 20]. Together, these effects modulate the number of star-forming galaxies in the early universe as a function of the DM-baryon relative velocity [21, 22, 23, 24, 25].

Even though we do not have direct access to the first star-forming galaxies at redshift $z\sim 20$, the hydrogen 21-cm line will allow us to indirectly map the distribution of the first stars, as their emission affects the 21-cm brightness temperature [26, 27, 28]. The presence of streaming velocities is often ignored in public 21-cm codes, such as 21cmFAST [29, 30], which is a good approximation if star formation mainly occurs in fairly massive haloes (with virial temperatures $T_{\rm cool}\gtrsim{\rm few}\times 10^{4}$ K), but will cause a severe overestimation if there is stellar formation in molecular-cooling haloes (with $T_{\rm cool}\lesssim 10^{4}$ K [31, 32, 33]). Here we present 21cmvFAST222https://github.com/JulianBMunoz/21cmvFAST , a modified version of 21cmFAST333https://github.com/andreimesinger/21cmFAST which properly accounts for the effects of the DM-baryon relative velocities.

Building upon the work of Refs. [9, 13, 24, 25], we include relative velocities in the initial conditions of 21cmvFAST, which then affects the amount of stellar formation at every cell in the simulation box. We additionally account for the averaged effect of Lyman-Werner (LW) feedback, which dissociates molecular hydrogen and raises the cooling temperature of haloes, reducing the effect of the streaming velocities at lower redshifts [34, 22, 35]. Thus, while for studies of the epoch of reionization, at $z\sim 6$, the streaming velocities are fairly unimportant (see, however, Ref. [36]), they have a large impact on the 21-cm signal during the cosmic-dawn era ($z=15-30$), as shown in Refs. [13, 24, 25, 24, 22, 23] and as we will further explore here.

We find that the DM-baryon relative velocities ($v_{\rm cb}$) produce a modest delay $\Delta z\sim 2$ on all the landmarks of cosmic evolution, as fewer stars form on average. More interestingly, however, the relative velocities show acoustic oscillations in their power spectrum, due to their BAO origin. These oscillations are imprinted in the distribution of the first stars, and therefore in the 21-cm power spectrum during cosmic dawn. This was first proposed in Ref. [13] during the Lyman-$\alpha$ coupling era (LCE), where the Lyman-$\alpha$ photons emitted by the first stars coupled the thermal and spin temperatures of neutral hydrogen [37, 38, 39]. A similar prediction was made in Refs. [25, 24] for the epoch of heating (EoH), where the X-rays emitted by the first galaxies heated the neutral hydrogen, thus affecting the 21-cm signal. Here we use 21cmvFAST to study the effect of the relative velocities across the entire cosmic-dawn era ($z=15-30$), comparing our simulations with analytic calculations.

Given the large astrophysical uncertainties during this epoch, we consider a broad range of possible feedback strengths. In all cases we find velocity-induced acoustic oscillations (VAOs) on the 21-cm power spectrum, with an amplitude growing in size during the LCE (peaking at $z\approx 24$ for our fiducial parameters), vanishing in the transition to the EoH, then peaking again during the heating era (at $z\approx 17$), and vanishing after the intergalactic medium (IGM) is fully heated ($z\lesssim 12$). This is in agreement with our analytic estimates. Furthermore, we confirm that the shape of the VAOs follows the well-understood fluctuations in $v_{\rm cb}^{2}$ [13, 40, 41], with damping at small scales due to the nonlocality of photon propagation. We parametrize this effect through a window function, as in Ref. [13], which dampens VAOs below $k\approx 0.05$ Mpc*-1* during the LCE, due to the large mean-free path of UV photons; and below $k\approx 0.2$ Mpc*-1* during the EoH, as X-ray photons typically get absorbed closer to their source [42]. Our analytically calculated VAOs and window functions very well match our numerical simulations during both the LCE and the EoH.

Finally, we study the observability of the VAOs with the upcoming hydrogen epoch of reionization array (HERA) [43]. We find that under most of our astrophysical scenarios, the VAOs are observable at large ($>5$) signal-to-noise ratios. This would confirm that the “small” molecular-cooling galaxies are driving stellar formation at cosmic dawn, and additionally yield the first measurement of baryon acoustic oscillations at a redshift intermediate between galaxy surveys and the CMB, which in Ref. [44] we propose as a novel standard ruler during cosmic dawn.

This paper is structured as follows. In Section II we review the effects of the streaming velocities on the formation of the first stars. We use these effects in Section III to find the changes in the 21-cm signal due to velocities, using 21cmvFAST to find 21-cm maps and power spectra. We quantify the effect of the velocities in detail in Section IV, and study its observability in Section V. We conclude in Section VI.

II The effect of velocities on the first stars

We begin by reviewing the effect of the streaming velocities on the first star-forming galaxies. This section heavily draws from Refs. [9, 12, 13, 18], so the reader might want to skip to Section III to see the impact on the observable 21-cm line. Throughout this work we use the best-fit cosmological parameters from Planck 2018 [5], and all lengths and wavenumbers will be comoving unless otherwise specified.

As in Ref. [9], we use a moving-mesh perturbation theory, which is valid for regions of radius $R\lesssim 3$ Mpc, where the velocity is coherent. Then, there exists a zero-order solution to the Boltzmann equations, with no overdensities,

[TABLE]

but a nonzero relative velocity between the DM and baryon fluids

[TABLE]

where $\mathbf{v}^{i}_{\rm cb}$ is the value of the DM-baryon relative velocity at kinematic decoupling, $z_{\rm kin}$ [45]. Additionally, we are free to choose a reference frame, which we set at the DM rest frame by setting

[TABLE]

which allows for easier computation at high wavenumbers $k$, where the baryon fluctuations oscillate rapidly [24]. Then, with this background solution we can rewrite the first-order Boltzmann equations from Ref. [9] as

[TABLE]

where $\theta_{i}$ is the velocity divergence of species $i$, $\mu\equiv\mathbf{k}\cdot\mathbf{v_{\rm cb}}/(kv_{\rm cb})$, $H$ is the Hubble expansion rate at redshift $z$, $v_{\rm cb}=|\mathbf{v_{\rm cb}}|$, and the prime denotes derivative with respect to redshift. Additionally, $\Gamma_{C}$ is the Compton interaction rate, $T_{\gamma}$ is the CMB temperature, and we include the effect of temperature fluctuations $\delta_{T}=\delta T_{g}/T_{g}$, where $T_{g}$ is the gas temperature. Here, the isothermal sound speed is given by [40]

[TABLE]

where $\mu_{b}$ is the average molecular mass of baryons and $k_{B}$ is the Boltzmann constant.

We solve Eqs. (4) for different values of $v_{\rm cb}^{i}$ and $\mu$, by setting the initial conditions to be $\delta_{b}=\mathcal{T}_{b}$, $\delta_{c}=\mathcal{T}_{c}$, and $\delta_{T}=0$, where $\mathcal{T}_{X}$ is the transfer function, which evolves the initial primordial fluctuations to over/underdensities in the fluid $X$, obtained with the Boltzmann solver CLASS [46]. With this we obtain the baryon and DM transfer functions $\mathcal{T}_{b}(k,z,v_{\rm cb}^{i})$ and $\mathcal{T}_{c}(k,z,v_{\rm cb}^{i})$, at arbitrary redshifts and relative velocities, which we will use to find the properties of the first star-forming haloes. Additionally, we define the matter transfer function $\mathcal{T}_{m}$ as usual, though due to their relative velocities the baryon and DM parts will have different phases.

Regions with large DM-baryon relative velocity $v_{\rm cb}$ form fewer stars through three mechanisms. First, less gas is able to collapse into the haloes, as the relative velocity allows gas to stream away. Second, the cooling of said gas in each halo is impeded by the relative velocity, raising the minimum halo mass necessary to form stars. Third, larger velocities result in lower matter fluctuations, due to the additional nonthermal “pressure”, and thus fewer haloes are formed. For a recent review of these effects see, for instance, Ref. [21]. We now compute each of these effects in turn.

II.1 Gas Fraction

As in Refs. [12, 13], we calculate the effect of the DM-baryon relative velocity on the fraction of gas accreted by DM haloes through the filter mass $M_{F}$. This quantity, first introduced in Ref. [47], allows us to compute the effect of the baryonic pressure, be it thermal or not, on the collapse of baryons to haloes, by using results from the linear regime. As in Refs. [48, 12] we fit the ratio

[TABLE]

to first order in $k^{2}$, where we use the notation of Ref. [48] in defining $r_{\rm LSS}$ as the low-$k$ large-scale structure (LSS) difference between baryons and all matter (which, however, we take at $k=1$ Mpc*-1* to avoid the BAO wiggles). We define the filter mass as [12]

[TABLE]

where $\rho_{m}^{0}$ is the matter density today, and the (comoving) filter scale is given by $r_{F}=\pi/k_{F}$.

While this fit is not perfect, it provides us with a good estimate of $M_{F}$ by using linear theory alone. Still, there are two important caveats to note here. The first is that the effect of the $v_{\rm cb}$-induced “pressure” in the baryon small-scale distribution has a different $k$-behavior from the usual thermal pressure, although we parametrize it in the same way. The second is that the filter mass is a real-space quantity, whereas the overdensities $\delta_{b}$ and $\delta_{m}$ that we computed are Fourier-space quantities. Therefore, we have to average over the (cosine of the) angle $\mu$ between $\mathbf{k}$ and $\mathbf{v_{\rm cb}}$ at some point, since real-space quantities cannot depend on it. We choose to average the ratio $\delta_{b}/\delta_{m}$ directly, instead of each density individually, although we have checked that this produces differences below 10% on the filter mass. Given the filter mass we can calculate the amount of matter in each halo that is composed of gas through the fit

[TABLE]

where $a=0.7$, and $f_{b}^{(0)}=f_{b}(1+3.2r_{\rm LSS})$ [12, 49].

Notice that, in addition to the effects calculated here, there can be resonant amplification at very high wavenumbers $k\sim 10^{3}$ Mpc*-1*, as proposed in Ref. [40], which are only relevant at scales below the haloes that we are interested in and thus can safely be ignored.

II.2 Minimum Halo Mass

In addition to lowering the fraction of baryonic gas in each halo, the relative velocity has a critical effect on the formation of stars in small-mass haloes through H2 molecular cooling. By softening the large gas overdensities that form in the core of these haloes, the relative velocities make it more difficult to form stars [16, 15, 17, 18, 20, 19]. We use the parametrization of Ref. [18], where this effect is cast as an increase in the minimum virial velocity that a halo requires to cool through molecular hydrogen as

[TABLE]

where we take the well-known ($v_{\rm cb}=0$) result of $V_{\rm cool}^{(0)}=4.0$ km s*-1* for molecular cooling, as it lies in between the results of the different simulations from Refs. [16, 15, 34]. The parameter $\alpha$ accounts for the effect of streaming velocities, and it is expected to lie in the range $\alpha=4-6$ [17]. We will study two different values within this range, with $\alpha=4$ as in Ref. [18] (fitted from the simulations of Refs. [16, 15]), or $\alpha=6$ as recent simulations seem to indicate [20, 19]. From $V_{\rm cool}$ we can define a minimum halo mass $M_{\rm cool}$ that can form stars as [49]

[TABLE]

during matter domination, and assuming our fiducial cosmological parameters, as only haloes with $M>M_{\rm cool}(z,v_{\rm cb})$ have virial temperatures large enough to produce stellar formation.

We additionally include the effect of Lyman-Werner (LW) feedback in the formation of stars. As gas collapses into stars, these will emit UV photons in the Lyman-Werner band (11.2$-$13.6 eV), which can easily photodissociate H2 [50, 51, 52, 53, 32]. This raises the minimum halo mass as the Universe evolves [34, 54, 55, 56], until it reaches the atomic-cooling threshold, which we take to be at $V^{\rm HI}_{\rm cool}=17$ km s*-1* [57, 58] (corresponding to $T_{\rm cool}=10^{4}$ K and $M_{\rm cool}\approx 3\times 10^{7}\,M_{\odot}$). The LW band overlaps with the Lyman-$\alpha$ band that produces Wouthuysen-Field coupling, which we discuss in Section III, albeit the relation between the two fluxes depends on astrophysical assumptions regarding the photon escape fraction, as well as self-shielding, among others [59, 60]. We will not attempt to model this process, and instead rely on the results of Refs. [34, 54, 35], where it was shown that the effect of LW feedback is to raise the mass of cooling haloes as

[TABLE]

where $M_{\rm cool}(z,v_{\rm cb},0)$ is the result without the LW feedback calculated above, and $F_{\rm LW}\equiv 4\pi J_{\rm LW}$ is the LW flux, with $J_{\rm LW}$ in the customary units of $10^{-21}$ ergs s*-1* cm*-2* Hz*-1* sr*-1*. Here, the two parameters $B$ and $\beta$ control the effect of the LW feedback on the cooling of gas in haloes.

Given the large astrophysical uncertainties in both the effect of velocities and the strength of the LW feedback, we will study four possible feedback strengths:

$\bullet$ No feedback, $B=0$. We do not expect this case to accurately reproduce our Universe, but it is a useful toy model where the effect of the velocities is most prominent in the 21-cm power spectrum.

$\bullet$ Low feedback, $B=4$ and $\beta=0.47$. This model follows the fits found in Refs. [34, 54], albeit with a coefficient $B$ roughly half of the best-fit value.

$\bullet$ Regular feedback, $B=7$ and $\beta=0.47$. This is the prescription in Refs. [22, 35], as a fit to the simulations of Refs. [34, 54], which results in a larger impact of $F_{\rm LW}$ on the minimum halo mass.

$\bullet$ High feedback, $B=7$ and $\beta=0.47$ again; however, here we assume that the effect of LW photons and the relative velocities is fully uncorrelated, as in Ref. [24], so that $V_{\rm cool}^{(0)}(F_{\rm LW})=V_{\rm cool}^{(0)}(0)\times\left[1+B\,(F_{\rm LW})^{\beta}\right]^{1/3}$ instead of the formulation of Eq. (11). This dampens the effect of velocities considerably, as we will see.

In each of these models we will take the average LW flux from the strong-feedback case in Ref. [22], where it is self-consistently computed through the amount of baryons collapsed into stars as a function of redshift, except for the low feedback strength, where we divide it by a factor of 2. Additionally, we take $\alpha=4$ in the regular- and high-feedback cases, and $\alpha=6$ in the no- and low-feedback cases, which enhances the $v_{\rm cb}$-induced fluctuations. We leave for future work to self-consistently compute the Lyman-Werner feedback in each cell of our 21cmvFAST simulation, similarly to how the Lyman-$\alpha$ flux is computed.

Note that we parametrize the unknown strength of the LW feedback through the poorly understood fitting function Eq. (11) ($B$ and $\beta$), where the coefficients are fit from a few (three) datapoints in Refs. [34, 54], as opposed to the redshift at which $F_{\rm LW}$ is evaluated, as in Ref. [22], since in Ref. [35] it was found that evaluating at the redshift of collapse is a good approximation unless $J_{21}$ grows extremely rapidly.

Correction to $f_{*}$

Following Refs. [34, 22], we assume that not all haloes above $M_{\rm cool}$ form stars with the same efficiency, as haloes near the threshold have a smaller fraction of their gas that can cool. We follow the parametrization in Ref. [22], where the fraction of baryons that form stars depends on the halo mass $M$ as

[TABLE]

for $M_{\rm atom}>M>M_{\rm cool}$, and $f_{*}=f_{*}^{0}$ for $M>M_{\rm atom}$, where $f_{*}^{0}$ is the usual input fraction of baryons in stars in 21cmvFAST. We set the scale at which the fraction plateaus at the mass $M_{\rm atom}$ of haloes that can form stars through atomic cooling, as in Ref. [22], which sets $\alpha_{*}=f_{*}^{0}/\log(M_{\rm atom}/M_{\rm cool})$ for $M_{\rm cool}<M_{\rm atom}$. We will take $f^{0}_{*}=0.1$ throughout, except in the unrealistic no-feedback case, where we do not include this log-correction, which we compensate by taking a lower value of $f^{0}_{*}=0.03$. For reference, $M_{\rm atom}\gg M_{\rm cool}$ at $z\geq 20$ for all feedback cases, whereas by $z=15$ the LW flux in the regular-feedback case is large enough to almost saturate, producing $M_{\rm cool}\approx M_{\rm atom}$.

II.3 Halo Mass Function

Through providing an additional source of pressure for baryons at small scales, the DM-baryon relative velocities can greatly affect the abundance—and bias—of the first star-forming haloes of our Universe, as first pointed out in Ref. [9], and confirmed with simulations in Ref. [10]. We follow their calculation here, albeit using the Reed halo mass function (HMF) from Ref. [61], which is similar to the Sheth-Tormen HMF [62], but better calibrated at the redshifts of interest ($z\sim 10-30$) to correctly account for a suppression of high-mass haloes.

We begin by calculating the local value of the matter power spectrum [9]

[TABLE]

by averaging over angles $\mu$, where $P_{\zeta}(k)$ is the primordial power spectrum. Then, we define the variance in matter overdensities that will form haloes of mass $M$ as

[TABLE]

where $R_{M}=[\Omega_{M}H_{0}^{2}/(2M)]^{-1/3}$ is the typical radius of these overdensities (with $H_{0}$ the Hubble expansion rate today and $\Omega_{M}$ the cosmic matter abundance), which we assume have a top-hat shape, so the window function is given by $W(x)=3x^{-2}[\cos(x)-\sin(x)x^{-1}]$. Note that we find $P_{m}$ by solving Eqs. (4) starting at the baryon drag era ($z_{d}=1060$), as suggested in Refs. [24, 63], as otherwise we would be underestimating the suppression due to velocities. Nonetheless, as 21cmFAST does not keep track of baryons and dark matter independently we cannot include any additional suppression due to differences in $\delta_{b}-\delta_{c}$ [64, 65], which could further enhance the effects of the streaming velocities.

We show the standard deviation of overdensities in star-forming haloes as a function of $v_{\rm cb}$, defined as $\sigma_{\rm cool}(z,v_{\rm cb})=\sigma_{M}(z,v_{\rm cb})$ for $M=M_{\rm cool}(z,v_{\rm cb})$, in Fig. 1, for redshifts $z=20$ and 30. We normalize all lines to their $v_{\rm cb}=0$ value, with values of $\sigma_{\rm cool}=\{0.52,0.42,0.41,0.41\}$ and $\{0.37,0.32,0.31,0.31\}$, for no feedback, low, moderate, and high feedback, at $z=20$ and 30, respectively. This figure shows that, in all cases, higher velocities require larger haloes to form stars, translating into lower values of $\sigma_{\rm cool}$. The size of this effect decreases with the strength of the LW feedback, as it raises the minimum cooling mass, and increases with redshift, as the velocities are more relevant at earlier times, although the size of the effect does not scale simply as $v_{\rm cb}(z)$, as it depends on its integrated history. Finally, we note that the “kinks” in some of the lines with feedback in Fig. 1 correspond to the transition to atomic-cooling haloes, where the effect of $v_{\rm cb}$ in $V_{\rm cool}$ (as in Eq. (9)) vanishes.

When computing the HMF we work with the dimensionless quantity [62, 66]

[TABLE]

where $\delta_{\rm crit}\equiv 1.69-\gamma\,(1+z)$ is the critical overdensity for collapse, with $\gamma=3\times 10^{-5}$ serving as a small correction to account for the effect of baryons at high redshifts [48]. The Reed HMF [61] has the following form:

[TABLE]

with $\beta_{1}=1.02$, $\beta_{G}=0.2$, $p=0.3$,

[TABLE]

$d=0.788$, and $\nu^{\prime}=a\nu$ with $a=0.764$. The amplitude $\mathcal{A}$ is found by normalizing the HMF at all masses, so that $\int d\nu f(\nu)=1$. Then, we find the HMF as a function of velocity to be

[TABLE]

We show this quantity (multiplied by $M^{2}$ to obtain the fraction of matter in haloes per unit log-mass) in Fig. 2, for three values of the relative velocity, where it is clear that regions with larger relative velocities form fewer small-mass haloes, whereas for atomic-cooling haloes with $M\gtrsim{\rm few}\times 10^{7}\,M_{\odot}$ the difference is less significant.

II.4 All Combined

Finally, we can combine all the aforementioned effects to calculate the fraction of baryons that collapse to star-forming haloes as [25]

[TABLE]

where $\rho_{b}$ is the comoving baryon abundance. This includes the velocity effect on the HMF (through $n(M)$); on the gas fraction $f_{\rm gas}$ of haloes; and on the minimum cooling mass $M_{\rm cool}$, all of which are functions of redshift and $v_{\rm cb}$. To illustrate the dependencies in this function, we show $F_{\rm coll}$ as a function of the minimum cooling mass $M$ in Fig. 2, for both the no-feedback and the regular-feedback cases. There we see how larger velocities decrease the number of haloes, the gas fraction in every halo, and increase the minimum cooling mass. These effects are less apparent in the case with feedback, although still present to a large degree.

We show the collapsed function of baryons that can form stars, from Eq. (19), in Fig. 3 as a function of redshift, for different values of $v_{\rm cb}$, as well as its average over velocities. From this figure we see that velocities have a larger impact at higher $z$, as their effects partially redshift away at later times, with the average $F_{\rm coll}$ being a factor of 2-5 smaller than the $v_{\rm cb}=0$ case. This illustrates the necessity of accounting for relative velocities if one is interested in molecular-cooling haloes, since otherwise one would form too many of them. This can be solved easily with 21cmvFAST.

As we will see, our observable—the 21-cm brightness temperature—depends on $F_{\rm coll}$ nonlocally, as photons emitted from a galaxy can travel far from its source. The approach of 21cmFAST to solve this problem is to use an excursion set approach [67, 68] to compute the amount of photons emitted at earlier times from distant sources. For that the code smooths the density field in spherical regions of different radii $R$ around each cell, and uses the averaged field to compute the collapse fraction, which tells us how many photons were emitted from this region. Computationally we use the extended Press-Schechter (EPS) formalism to account for over- and under-densities on top of velocity fluctuations [69]. Assuming that we take a patch of radius $R$, with averaged overdensity $\delta$ and relative velocity $v_{\rm cb}$, and $\sigma_{R}^{2}$ is the matter density variance in patches of radius $R$, we follow Refs. [69, 12] in writing the collapsed fraction of baryons in this patch as

[TABLE]

where $F_{\rm EPS}(\delta,v_{\rm cb})$ is the EPS collapsed fraction, given by

[TABLE]

where erfc is the complementary error function, and $\overline{F_{\rm EPS}}$ is the global average of $F_{\rm EPS}$. We follow the prescription in 21cmFAST by taking the nonlinear overdensities in $\delta$ (i.e., extrapolating the $z=0$ result backwards), and calculating $\overline{F_{\rm EPS}}$ numerically by averaging over our cells [29, 30].

Notice that in the limit $R\to\infty$, the collapsed fraction in Eq. (21) asymptotes to $F_{\rm EPS}(0,\left\langle v_{\rm cb}\right\rangle)$, which is not necessarily the same as $\left\langle F_{\rm EPS}(\delta,v_{\rm cb})\right\rangle$ (where $\left\langle\right\rangle$ denotes spatial average). Nonetheless, as we divide by the actual averaged value of $\overline{F_{\rm EPS}}$, any bias that could arise from this procedure is cancelled out. Additionally, as shown in Figs. 2 and 3, the averaged $F_{\rm coll}$ is almost identical to its value for $v_{\rm cb}=\left\langle v_{\rm cb}\right\rangle$. This is because the probability density function (PDF) for $v_{\rm cb}^{2}$ is highly peaked around $v_{\rm rms}^{2}$, where $v_{\rm rms}\approx 30$ km s*-1* is the root mean square of the Maxwell-Boltzmann distribution of $v_{\rm cb}$, so using the mode value gives very similar results to averaging over the PDF. Additionally, for computational efficiency we find the $\sigma_{\rm cool}(v_{\rm cb})$ (as plotted in Fig. 1) by calculating $F_{\rm coll}$ with a Press-Schechter HMF [70] including all velocity effects, as it takes the simple form in Eq. (21) for $\delta=\sigma_{R}=0$, and one can find $\sigma_{\rm cool}(v_{\rm cb})$ by inverting the erfc function.

III The 21-cm line

The distribution of the first stars during cosmic dawn cannot be directly probed. Fortunately, their presence can be indirectly inferred with measurements of the 21-cm line, which we now describe.

III.1 The Basics

We start by briefly reviewing the physics of the 21-cm hydrogen line (for a detailed review we encourage the reader to see Refs. [26, 27, 21]). Hydrogen, in addition to being the most abundant element in our Universe, possesses a rich hyperfine structure. Of particular interest is its spin-flip transition from the triplet to the singlet ground state, which produces a photon with a 21-cm wavelength. To determine whether cosmological hydrogen emits or absorbs 21-cm photons during a certain era, one has to calculate its spin temperature and compare with that of the CMB. Given a comoving number density of hydrogen atoms in the triplet ($n_{1}$) and singlet ($n_{0}$) states, we can define the local spin temperature $T_{s}$ through

[TABLE]

where $g_{1}=3$ and $g_{0}=1$ are the numbers of degrees of freedom of the triplet and singlet states, respectively, and $T_{*}\approx 0.068$ K is the temperature corresponding to the 21-cm hyperfine transition. If $T_{s}>T_{\gamma}$, the hydrogen will emit 21-cm photons, and if $T_{s}<T_{\gamma}$, it will absorb photons from the Raleigh-Jeans tail of the CMB. We can define the brightness temperature of the 21-cm line as [27]

[TABLE]

for our fiducial cosmological parameters, where $x_{\rm HI}$ is the neutral-hydrogen fraction (nearly unity for all eras of interest in this work), and $\partial_{r}v_{r}$ is the line-of-sight gradient of the velocity.

The standard history of the 21-cm line of hydrogen is as follows. During the dark ages collisions coupled the spin temperature to that of the gas, producing 21-cm absorption for $z\gtrsim 30$ [71]. This process, however, became inefficient as the Universe expanded and the gas cooled further, resulting in nearly no absorption until the cosmic-dawn era. In this epoch the first stars were formed, emitting abundant UV photons, which coupled $T_{s}$ to $T_{g}$ through the Wouthuysen-Field effect [37, 38, 39], by which Lyman-$\alpha$ photons resonantly scatter between hydrogen atoms, imprinting the kinetic temperature onto the hyperfine populations of hydrogen. This produces—once again, 21-cm absorption—and we call this period the Lyman-$\alpha$ coupling era (LCE). Later on, X-ray photons will be produced in a large-enough quantity to reheat the intergalactic medium to temperatures $T_{g}\gg T_{\gamma}$, switching the 21-cm signal to emission. We call this period the epoch of heating (EoH). Finally, during the epoch of reionization, at $z\lesssim 10$, there will be enough UV photons that all neutral hydrogen will be ionized, until no 21-cm signal is left.

We are interested in the 21-cm line during cosmic dawn ($z\sim 20$), when the spin temperature is set by the competition between couplings to the CMB and Lyman-$\alpha$ backgrounds, and the gas temperature is first dominated by adiabatic cooling of the gas, and then rises due to X-ray heating. Given the complexity of this signal, the most straightforward approach is to perform quasi-numerical simulations. We do so with 21cmvFAST.

III.2 Implementation on 21cmvFAST

We include the effects of relative velocities in the commonly used 21cmFAST code [29], by modifying the streamlined version called 21CMMC presented in Refs. [30, 72]. Here we outline the main changes to the code, and the results found, and we encourage the reader to visit the webpage with the modified code, where all changes are detailed.

Each realization of 21cmvFAST is composed of a box of side $L$ and cell size $R_{\rm cell}$, with a set of initial densities $\delta^{(i)}(\mathbf{x})$ given by the matter power spectrum, evolved using linear or second-order Lagrangian perturbation theory to lower redshifts [73]. This is an excellent approximation for the high redshifts ($z\gtrsim 10$) that we are interested in. We will always choose $R_{\rm cell}\leq 3$ Mpc, so that our moving-background perturbation theory in Eq. (4) holds. The first modification that we implement is to generate a box of initial DM-baryon relative velocities. We do so (keeping correlations with over/under-densities) by assigning each Fourier-space pixel with density $\delta^{(i)}(\mathbf{k})$ a velocity

[TABLE]

where $P_{v}$ and $P_{\delta}$ are the power spectra for $v_{\rm cb}$ and $\delta$, respectively, which we Fourier transform to real space, to find $\mathbf{v_{\rm cb}}(\mathbf{x})$ at each pixel (and thus $v_{\rm cb}(\mathbf{x})$ by taking its norm). We show a slice of the initial conditions in our simulation in Fig. 4, where it is clear that overdensities fluctuate on all scales, whereas the relative velocities are clearly coherent for scales below $\sim 10$ Mpc, and most of their power lies in BAO scales ($R\approx 150$ Mpc).

We use our calculations from Section II to find $F_{\rm coll}$ for each point in the box at $z=20$, assuming regular feedback strength and given our initial conditions, both including and neglecting the effects of the relative velocities, which we show in Fig. 4. This figure clearly showcases how regions with high velocity see a significant suppression in the amount of baryons collapsed to galaxies, and thus of stars formed. We have smoothed $F_{\rm coll}$ over different radii, to represent the characteristic distance traveled by different kinds of photons. X-ray photons can travel small distances before being absorbed, if their energies are low, so we have averaged over radii $R=3$ and 20 Mpc, corresponding to energies of 0.2 and 0.5 keV, respectively. On the other hand, during cosmic dawn the UV photons that produce Lyman-$\alpha$ coupling will typically travel as far as $R=100$ Mpc before redshifting into the Lyman-$\alpha$ transition, as they typically do not get absorbed. We see from Fig. 4 how smoothing over larger $R$ partially dampens the effect of velocities at smaller scales, as expected, although even for $R=100$ Mpc there is a visible signature of $v_{\rm cb}$ in the collapsed fraction, which will be inherited by the 21-cm power spectrum.

III.3 Global Signal

Given the implementation of 21cmvFAST outlined above, we can now calculate the 21-cm brightness temperature for any combination of cosmological and astrophysical parameters, including molecular-cooling haloes with a time-dependent LW feedback, and the effect of relative velocities.

In all of our simulations we take an X-ray luminosity of $\log_{10}(L_{X}/\rm SFR)=40$ (divided by the star-formation rate SFR), and a threshold in X-ray energy of $E_{0}=0.2$ keV [25], with a spectral index of $\alpha_{X}=1$, and do not modify the rest of the standard 21CMMC parameters. We will also explore the case of $E_{0}=0.5$ keV in App. A, as the first X-ray sources might have harder spectra than we assume [74, 75, 76]. Instead of a fixed cooling temperature $T_{\rm cool,Xray}=10^{4}$ K for X-ray haloes we use the prescription for $F_{\rm coll}$ and $\sigma_{\rm cool}$ outlined in Section II. Additionally, we have not modified the ionization part of 21cmFAST, assuming that the haloes we study here provide a negligible amount of UV photons that escape into the IGM. Therefore, for ionizations we employ the “Faint Galaxies” model of Ref. [77], where reionization is driven by galaxies with $T_{\rm cool,UV}=2\times 10^{4}$ K, with an ionizing efficiency of $\zeta=20$, and by the end of our simulations, at $z=12$, we always have $x_{e}\ll 1$.

We show in Fig. 5 the average 21-cm brightness temperature $\left\langle T_{21}\right\rangle$ in our simulation boxes as a function of redshift. The inclusion of relative velocities delays the formation of the first stars, shifting all landmarks from the “plain” molecular-cooling case (i.e., ignoring the effects of both $v_{\rm cb}$ and Lyman-Werner feedback) $\Delta z=2$ later, in agreement with the results in Ref. [22]. Accounting for LW feedback as well produces a further delay of $\Delta z=3$, due to the need for more-massive haloes to form stars at lower redshifts. As a consistency check, we also show in that plot the result for a simulation with fixed $v_{\rm cb}=v_{\rm avg}=\left\langle v_{\rm cb}\right\rangle=0.92v_{\rm rms}$, where no velocity-induced fluctuations are present. We see that in this case the average 21-cm temperature is very similar to the full case. So while $v_{\rm cb}$ affects the background evolution of the 21-cm line in a particular way, these effects can potentially be mimicked by time-dependent feedback on the star-formation rate, and are thus not unique. We also show, for comparison, the background evolution of $T_{21}$ in the unphysical case that LW feedback is saturated at all times, and only atomic-cooling galaxies are allowed to form stars. We elaborate further on this case in App. B.

III.4 Fluctuations and Power Spectrum

Fortunately, there is more information in the 21-cm line than the global signal carries. Including $v_{\rm cb}$ and LW feedback delays the formation of structure, which affects the 21-cm fluctuations during the cosmic-dawn era in two ways. First, larger velocities result in fewer Lyman-$\alpha$ photons, producing anisotropies during the LCE [13, 17, 23]. Secondly, even after Lyman-$\alpha$ coupling is saturated, larger velocities result in fewer X-ray photons being emitted, and therefore less IGM heating [25, 17]. Following these two effects we divide the cosmic dawn into two different eras, the EoH and the LCE, which have qualitatively different fluctuation behaviors. This is illustrated in Fig. 6, where we show the 21-cm temperature in slices 1 Mpc thick of three simulation boxes (each 300 Mpc in length), across the entire cosmic dawn, starting at $z=35$ and finishing at $z=12$, for the same three cases as Fig. 5 (i.e., ignoring both $v_{\rm cb}$ and LW feedback, ignoring only the feedback, and including the full set of effects). Here it is very clear how the inclusion of $v_{\rm cb}$ and LW feedback delays the onset of the LCE, as well as the transition to the EoH, as the first galaxies take longer to emit the necessary amount of UV and X-ray photons. We now study the 21-cm fluctuations during each of these two eras in detail, focusing on our regular-feedback case, with the full set of feedback effects, to explore the effect of relative velocities.

Throughout this section we will employ two types of simulations. We will show 21-cm maps with a box size of $L=300$ Mpc and a 1-Mpc resolution, as those are easiest to visually inspect.444We note that whenever we use small simulation boxes, with $L=300$ Mpc, we do so for illustration purposes only, as for small box sizes part of the $v_{\rm cb}$ power at large scales is lost. For instance, for $L=300$ Mpc almost 10% of the $v_{\rm cb}$ power is not properly included. These boxes share the density and velocity fields from Fig. 4. When computing power spectra, however, we will increase the box size to $L=900$ Mpc, with a 3-Mpc resolution, in order to compute the large-scale 21-cm fluctuations with smaller variance. In all cases where we compare simulations with and without $v_{\rm cb}$ we keep an average value $v_{\rm avg}$ for the latter, so that both cases have very similar global signals.

During the Lyman-$\alpha$ Coupling Era

We start with the LCE, which stretches from the formation of the first abundant sources of Lyman-$\alpha$ photons, to the transition to X-ray heating, where $\left\langle T_{21}\right\rangle$ is minimum. This corresponds to the redshift range $z=20-28$ in our fiducial model. During this era, higher velocities produce a shallower 21-cm signal by reducing the amount stars emitting UV photons able to produce Lyman-$\alpha$ coupling. This is, however, a fairly nonlocal effect, as these photons can travel distances as large as 100 Mpc before transitioning into the Lyman-$\alpha$ line, which will smear the VAOs [13].

We show 21-cm maps during the LCE in Fig. 7, both with and without velocity fluctuations. This Figure shows three representative redshifts, as marked by blue and green circles in Fig. 5, as landmarks of the LCE. The first redshift we show is $z=28.3$, at the onset of Lyman-$\alpha$ coupling. Here regions of large velocity produce shallower 21-cm absorption (shown in red), whereas regions with $v_{\rm cb}\approx 0$ produce deeper absorption (shown in blue). This effect is perhaps more apparent at $z=25.0$, halfway through the Lyman-coupling era (where $\left\langle T_{21}\right\rangle=0.5\times\left\langle T_{21}\right\rangle_{\rm min}$), where regions of large $v_{\rm cb}$ translate into red “islands” in Fig. 7. Finally, we show the 21-cm fluctuations at $z=20.3$, in the transition between Lyman-coupling and X-ray heating (where $\left\langle T_{21}\right\rangle$ is at its minimum). Here we find negligible large-scale VAOs. This is to be expected, since the effect of the velocities on the 21-cm brightness temperature during the EoH and LCE are opposite, and ought to cancel out at some intermediate redshift, as we will quantify further below.

The simulation boxes that we have shown so far clearly indicate that DM-baryon relative velocities are important for the evolution of the 21-cm signal. Nonetheless, to more precisely quantify their effect we will now study the power spectrum of the 21-cm signal. Given the fluctuations $\delta T_{21}$ in the 21-cm temperature, we define the 21-cm power spectrum through the two-point function as usual,

[TABLE]

where $\delta_{D}$ is the Dirac delta function. For simplicity we will use the amplitude of 21-cm fluctuations, defined as

[TABLE]

where we drop the subscript “21” and the vector symbol for $k$, as we ignore line-of-sight (LoS) anisotropies in this work. Hereafter we refer to $\Delta^{2}(k)$ as the 21-cm power spectrum for simplicity. In all cases the power spectrum of our simulations, and its errorbars, are obtained from a FFTW of our simulation boxes, as in the standard 21cmFAST [78, 29], albeit with a finer $k$-binning to allow us to better distinguish the large-scale acoustic oscillations.

We show in Fig. 8 the 21-cm power spectrum at the same three redshifts as in Fig. 7, where we see how the fluctuations behave as a smooth function of wavenumber $k$, but at large scales are imprinted with acoustic oscillations. These do not behave as the regular (density-induced) BAOs and are, in fact, velocity-induced acoustic oscillations (VAOs), sourced from the fluctuations on $v_{\rm cb}$. These VAOs are large at $z=25.0$, even dominating the power spectrum at large scales, but become largely irrelevant at lower $z$ as the Universe becomes fully Lyman-$\alpha$ coupled. The size of the VAOs also decreases at higher redshifts, although as long as there is some Lyman-$\alpha$ coupling the impact of the velocities will be important. The broadband (non-VAO) 21-cm power spectrum follows a similar evolution. We see in Figs. 7 and 8 how in the transition from the LCE to the EoH all the power resides at small scales, as the first X-ray emitting sources heat up the gas around them preferentially. We note that the regular density BAOs are included through the total-matter fluctuation $\delta^{(i)}(\mathbf{k})$, whose transfer function we obtain from CLASS [46], but their amplitude is too small (and smoothed by the complicated astrophysics of the cosmic dawn) to be observable in our power spectra.

During the Epoch of Heating

We now move to the EoH, which we define to begin at the minimum of $\left\langle T_{21}\right\rangle$, and to last until $T_{s}\gg T_{\gamma}$, where the entire box emits in 21 cm ($T_{21}>0$). For our fiducial parameters the EoH extends from $z=20$ to $z=15$, and throughout the entire EoH $T_{g}$ and $T_{s}$ are very efficiently coupled. During this era the 21-cm fluctuations are dominated by gas-temperature anisotropies due to the inhomogeneous X-ray heating of the IGM [79, 80]. The way the relative velocities will impact the 21-cm signal is, then, by changing the amount of X-ray heating in different regions. Patches with large velocities will form fewer stars, producing less heating, and therefore causing deeper 21-cm absorption.

As in the LCE case we will show results for three representative redshifts, now in Fig. 9. The first is $z=17.2$, halfway through the heating transition (where $\left\langle T_{21}\right\rangle=0.5\times\left\langle T_{21}\right\rangle_{\rm min}$). In this case there are marked 21-cm fluctuations at large scales due to the DM-baryon relative velocity, since regions with large relative velocities produce a more-negative $T_{21}$. The same is true at the transition to emission, defined as $\left\langle T_{21}\right\rangle=0$, which corresponds to $z_{0}=14.5$ in our model. The last case we show is at $z=12.0$, at the end of our simulations, at which point the IGM has been heated to $T_{s}\approx T_{g}\gg T_{\gamma}$, making the $T_{s}$ term of Eq. (23) negligible. In this case clearly the DM-baryon relative velocity does not have an important effect, as $T_{21}$ has saturated. The smooth (non-VAO) power follows roughly the same trajectory, as it is driven by the anisotropies in the distribution of stars, and thus of X-ray photons.

As in the previous section we plot the 21-cm power spectrum during the EoH in Fig. 10. For the highest redshift, $z=17.2$, we find very important VAOs at $k\sim 0.1$ Mpc*-1*, much larger than during the LCE and not severely damped. As time progresses, X-rays slowly heat up the IGM, making the power spectrum smaller (and the VAOs less pronounced) at $z=14.5$. We do not find a trace of VAOs in the $z=12.0$ case, as expected, since in our model $T_{21}$ has saturated by then.

Notice that here, and throughout this work, we ignore any explicit LoS anistropies in the signal. There are, however, some effects that can make the signal anisotropic, such as redshift-space distortions [81, 82] and the effect of the lightcone [83, 84] (as we only take co-evaluated boxes). These effects can be sizable on large scales [85], and potentially affect recovery of parameters from data [86, 87]. Nonetheless, the VAO signal is expected to be isotropic, so we leave detailed studies of the effect of LoS anisotropies for future work.

To summarize, at small scales ($k\gtrsim 0.3\,\rm Mpc^{-1}$) the velocities have little direct effect on the power spectrum (although of course they indirectly affect it by changing the background). Nonetheless, on large scales ($k=0.03-0.3\,\rm Mpc^{-1}$) the 21-cm signal presents distinct VAOs, produced by the effect of the DM-baryon relative velocities on the first galaxies.

IV Quantifying the Effect of $v_{\rm cb}$

So far we have established that the DM-baryon relative velocities affect the 21-cm signal during the entire cosmic-dawn era, even when accounting for Lyman-Werner feedback, producing striking VAOs on large scales. Let us, however, confirm that these acoustic oscillations are indeed produced by the relative velocities, and not over/under-densities. For that, we perform simulations in which we fix $v_{\rm cb}$ to its average value, $v_{\rm avg}$, and we compare in Fig. 11 the 21-cm power spectrum with our regular case at two redshifts, $z=17.2$ (halfway through the EoH) and $z=25.0$ (halfway through the LCE). This choice of a nonfluctuating $v_{\rm cb}=v_{\rm avg}$ keeps roughly the same 21-cm global signal, as well as small-scale power spectrum. Nonetheless, as is clear from Fig. 11, the addition of velocity fluctuations produces additional power on large scales, with acoustic wiggles. This is what we call the VAOs. We now move on to quantify the shape and amplitude of these VAOs, comparing them with analytic expectations.

IV.1 Shape of the Effect

The streaming velocities fluctuate with a power spectrum sourced by the BAOs [9], so that the power spectrum of any function $f(v_{\rm cb})$ (such as $F_{\rm coll}$ or $\delta T_{21}$) will have a component sourced by $v_{\rm cb}$ fluctuations. Two facts will greatly simplify our discussion about velocities moving forward. First, for several examples it was shown in Refs. [13, 40, 41] that

[TABLE]

where $b_{f}$ is a ($k$-independent) bias, and $\Delta^{2}_{v^{2}}$ is the power spectrum of the quantity

[TABLE]

which has no units, zero mean, and unit variance555Whenever we calculate $v_{\rm cb}$ or $\Delta^{2}_{v^{2}}$ it will always be at the drag era, unless otherwise specified.. This approximation is, in fact, accurate to the percent level for smooth functions $f(v_{\rm cb})$, as we confirm in App. C. Secondly, density and relative-velocity fluctuations are uncorrelated to first order, so we can simply add the 21-cm power spectrum induced by $\delta_{v^{2}}$ to the one due to density fluctuations to find the total signal. Given these two conditions, finding the impact of DM-baryon relative velocities on the 21-cm power spectrum amounts to calculating the shape and size of the term proportional to $\Delta^{2}_{v^{2}}$, which we have dubbed the VAOs. These two conditions additionally allow us to use the VAOs as a standard ruler during cosmic dawn, as we explore in Ref. [44].

Notice that in Fig. 11 the VAOs do not have the same shape during the EoH ($z=17.2$) and the LCE ($z=25.0$), as in the latter case they are damped on small scales ($k\gtrsim 0.1\,\rm Mpc^{-1}$). This is due to the propagation of Lyman-$\alpha$ photons, which typically travel much longer paths than their X-ray counterparts, erasing small-scale $v_{\rm cb}$ fluctuations, as predicted in Ref. [13]. We will detail how we include this effect for both X-ray and Lyman-$\alpha$ photons. We will follow an approach similar to Refs. [13, 24].

IV.1.1 Lyman-$\alpha$ Pumping

The number $n_{i}$ of X-ray or Lyman-$\alpha$ photons present at a point $\mathbf{x}$ at redshift $z$ depends on previous redshifts $z^{\prime}$ and positions $\mathbf{x}^{\prime}$ in the past lightcone. In particular, for Lyman-$\alpha$ photons this number is given by the integral over all previous redshifts

[TABLE]

assuming a flat spectrum of UV photons, where $\eta$ is conformal time, $N_{\alpha}$ is the number of Lyman-$\alpha$ photons emitted per collapsed baryon, and $c$ is the speed of light. Therefore, a fluctuation in $\dot{F}_{\rm coll}$ due to velocities will affect the number of photons up to a significant distance away. This nonlocality can be recast as a window function in Fourier space, where now the fluctuation in the number of Lyman-$\alpha$ photons can be written as

[TABLE]

where $\delta_{F}(\mathbf{k},z)$ is the fluctuation in $F_{\rm coll}$, and $W_{\alpha}(k,z)$ is the window function, which we calculate as [13]

[TABLE]

where $\mathrm{sinc}(x)=\sin(x)/x$, and we integrate up to the horizon redshift for Lyman-$\alpha$ absorptions, given by $z_{\rm hor}(z)=(32/27)(1+z)-1$, as photons from higher redshifts (with correspondingly higher frequencies) would have been absorbed as Lyman-$\beta$ [49]. Here, and throughout this discussion, we have taken the simplifying assumption that $z^{\prime}\sim z$, so that $\eta-\eta^{\prime}\approx(z^{\prime}-z)/H(z^{\prime})$, and we show the resulting window function in Fig. 12, which clearly damps fluctuations in Lyman-$\alpha$ photons for wavenumbers larger than $k\sim 0.01$ Mpc*-1*, due to the large mean-free path of these photons. Nonetheless, the window still shows some power, albeit highly suppressed, at $k=0.1$ Mpc*-1*, signaling that Lyman-$\alpha$ pumping can potentially produce 21-cm fluctuations on scales that can be reached by upcoming 21-cm observatories.

IV.1.2 X-ray Heating

For the case of X-rays the situation is slightly different, as their typical mean-free path is significantly shorter than that of Lyman-$\alpha$ photons. We can, nonetheless, define a window function for this scenario as well. We will follow the calculation in 21cmFAST, and assume that all photons with optical depths $\tau\geq 1$ are absorbed, whereas the rest travel freely [29]. This makes the mean-free path of an X-ray photon dependent on its energy. We can, then, define the minimum energy $E_{\rm min}(z,z^{\prime})$ that a photon can have to travel from $z^{\prime}$ to $z$, by solving for $\tau=1$, or equivalently, for $\lambda_{X}(E_{\rm min})=c(z^{\prime}-z)/H(z)$, where $\lambda_{X}$ is the comoving mean-free path (from Ref. [42]), to find

[TABLE]

capping it at $E_{\rm min}=E_{0}$ if lower. Here, and throughout, we assume that the spectrum of emission is log-flat in energy, (i.e., $L_{X}\propto E^{-\alpha_{X}}$ with $\alpha_{X}=1$), from a threshold energy $E_{0}=0.2$ keV, below which we assume all X-rays are self-absorbed near the sources, up to $E_{\rm max}=2$ keV, above which X-ray photons have mean-free paths large enough to barely heat up the IGM [75, 74, 72].

Under these assumptions, the gas heating is given by an integral of the form [29]

[TABLE]

where $E_{i}$ and $\sigma_{i}$ are the energy threshold and cross section for ionization. For our purposes it suffices to take the large-$E$ behavior of $\sigma_{i}\propto E^{-7/2}$, to find that the amount of energy deposited into the IGM by X-rays with a minimum energy $E_{\rm min}$ scales as $\epsilon_{X}\propto E_{\rm min}^{-7/2}$. We can, therefore, define the ratio $\mathcal{R}_{F}$ of the X-ray flux $\mathcal{F}$ deposited by photons above some energy $E_{\rm min}$, compared to the full spectrum down to $E_{0}$, as

[TABLE]

Then we can write the window function in analogy with the Lyman-$\alpha$ case simply as

[TABLE]

Using all these results we show the X-ray window function in Fig. 12, along with that of Lyman-$\alpha$ photons, each at their respective relevant redshifts. Interestingly, we see that X-ray photon anisotropies are not severely damped within the scales of interest ($k\lesssim 0.3$ Mpc*-1*), where velocity fluctuations are important. This is in agreement with the X-ray window function obtained for regular overdensities in Ref. [88]. We further illustrate this point in Fig. 13 by showing $\Delta^{2}_{v^{2}}$ (the power spectrum of $\delta_{v^{2}}$) multiplied by window functions (squared) due to either X-ray or Lyman-$\alpha$ propagation, at $z=17.2$ and $z=25.0$ respectively. We see how accounting for X-ray propagation dampens the small-scale power, while leaving the interesting region of $0.01\leq k\,{\rm Mpc}\leq 0.3$ roughly unaltered, whereas the propagation of Lyman-$\alpha$ photons dampens power at all relevant scales, making the power spectrum difficult to observe beyond $k\gtrsim 0.1$ Mpc*-1*.

IV.2 Size of the Effect

Now that we understand how photon propagation alters the shape of the $v_{\rm cb}$ fluctuations, let us fit for the amplitude of the VAOs in our simulations.

As mentioned at the beginning of this section, the VAOs are, to first order, statistically independent from density fluctuations [13]. We, therefore, choose to model the 21-cm power spectrum simply as

[TABLE]

with the VAOs parametrized as

[TABLE]

where $A_{\rm vel}(z)$ is a redshift-dependent amplitude of fluctuations, with units of mK2, which we will fit from the data, and we set the window function to be $W_{i}(k,z)=W_{X}(k,z)$ during the EoH, and $W_{\alpha}(k,z)$ during the LCE. The non-VAO component, $\mathcal{P}_{n}$, holds a wealth of information on astrophysical parameters [89, 90, 91, 43, 72, 92]. Nonetheless, as here we are only interested in the VAOs we will model $\mathcal{P}_{n}$ as a simple fourth-order polynomial,

[TABLE]

which captures the smooth behavior of the power spectrum without wiggles. This $\mathcal{P}_{n}$ is not designed to provide an excellent fit to the 21-cm power spectrum, but merely to allow us to marginalize over the smooth part of the signal.

We employed the fit from Eq. (36) in Figs. 8, 10, and 11, to compare with the simulation results. In these fits we only used modes with $k\geq 0.05$ Mpc*-1*, to avoid the sizable error bars of large-scale modes in the boxes. We see that in all cases our model provides a reasonable fit for the 21-cm fluctuations over the entire relevant range where velocities are important ($0.03\,\rm Mpc^{-1}\leq{\it k}\leq 0.3\,Mpc^{-1}$). For reference, at the redshifts that we showed in Figs. 8 and 10 the value of the VAO amplitude is $A_{\rm vel}=\{330,1100,36,270,27,0.02\}$ mK2, for decreasing redshift. This amplitude always increases the large-scale 21-cm power spectrum, allowing for an easier detection of this signal.

This $A_{\rm vel}$ is, however, not directly observable, as it depends on the assumptions about the shape of the VAOs (i.e., the window functions). In order to better quantify the size of the VAOs we define the variance in $T_{21}$ induced by the velocities alone as

[TABLE]

with $\Delta_{21,\rm vel}^{2}$ as defined in Eq. (37), which depends on $A_{\rm vel}$ and takes into account the suppression in fluctuations due to photon propagation. We show this VAO-sourced variance as a function of redshift in Fig. 14, where we can distinguish the two peaks due to the Lyman-$\alpha$ coupling and X-ray heating, the former peaking at $z\approx 24$ and the latter at $z\approx 18$. Additionally, we reach the expected result that at high redshifts there are not enough UV photons for efficient Lyman-$\alpha$ coupling (so the VAO amplitude all but vanishes above $z\gtrsim 30$), whereas at small redshifts ($z\lesssim 14$) the gas is heated up above $T_{\gamma}$, both of which render the effect of VAOs negligible. Interestingly, the VAO amplitude cancels in the transition from the LCE to the EoH. This is because the effect of velocities on the 21-cm line during the LCE and the EoH goes in opposite directions. During the LCE a large $v_{\rm cb}$ results in higher $T_{21}$, as fewer Lyman-$\alpha$ photons couple $T_{s}$ to $T_{g}$. On the other hand, during the EoH larger velocities produce a more-negative $T_{21}$, as fewer X-ray photons heat up the gas. This means that, by necessity, these two components have to cancel at some point, resulting in nearly zero VAO fluctuations, which in our simulations occurs at $z\approx 21$.

In addition to the VAO-sourced 21-cm variance $\delta T_{21,\rm vel}^{2}$, we show in Fig. 14 the VAO-only power spectrum $\Delta_{21,\rm vel}^{2}$ as a function of redshift at a reference scale $k_{\rm ref}=0.14\,\rm Mpc^{-1}$, computed with Eq. (37). Here the effect of the window functions is obvious, as the fluctuations on $T_{21}$ during the LCE are significantly smaller than during the EoH at this wavenumber, even though the effect of the velocities in $F_{\rm coll}$ is larger at higher redshifts. As a test, we also show in Fig. 14 results for a case in which we fix $v_{\rm cb}$ to its average value, obtained through the same fitting procedure as our main result. This curve can be taken as an estimate of the errors in both the fitting procedure and our simulations, and is significantly smaller than the amplitude in the case with velocities. We therefore conclude that our procedure is properly capturing both the shape and amplitude of the VAOs.

IV.3 Comparison with Analytic Calculations

We now estimate the expected size of the VAO signal from analytic arguments, extending on the work of Refs. [13, 24] by including X-ray and Lyman-$\alpha$ fluctuations jointly. During cosmic dawn we can safely neglect collisional hyperfine transitions, so the spin temperature will be given by [27]

[TABLE]

where $x_{\alpha}\propto n_{\alpha}$ accounts for the strength of Lyman-$\alpha$ pumping (see Eq. (29)). This allows us to recast Eq. (23) and write the fluctuations in the 21-cm brightness temperature as

[TABLE]

where we have ignored density fluctuations in this simplistic model, and $\kappa(z)$ is a correction factor that we will describe below.

As illustrated above, the DM-baryon relative velocities modulate the collapse fraction $F_{\rm coll}$, which in turn alters the 21-cm temperature through both the $x_{\alpha}$ and $T_{g}$ terms. In this section we will follow our simple model in which the amount of X-ray and Lyman-$\alpha$ photons are both proportional to $F_{\rm coll}$, so that

[TABLE]

where $T_{g}^{(\rm ad)}$ is the gas temperature considering only adiabatic cooling, and $\Delta T_{X}$ is the change in gas temperature due to X-ray heating. For practical purposes we use the box-averaged values of $x_{\alpha}$ and $\Delta T_{X}$ from our simulations, instead of computing them from first principles. Additionally, given that we are ignoring density fluctuations in this calculation, we will not obtain the correct background evolution. We thus include a correction factor $\kappa(z)=\left\langle T_{21}\right\rangle_{\rm sim}/\left\langle T_{21}\right\rangle_{\rm ana}$ in Eq. (41), in order to recover the correct average 21-cm signal.

With our toy model in Eq. (41), we can calculate the amplitude of the VAOs using the result of Refs. [13, 40, 41]—here recast as in Eq. (27)—and defining the 21-cm velocity bias as

[TABLE]

Then, the VAO component of the power spectrum will be

[TABLE]

which has the theoretically motivated shape, and an amplitude given by $b_{21}^{2}$. We compare the results for this analytic calculation with those of our simulations in Fig. 14, where we can see that the main features of the signal are well reproduced by this approximation. For instance, we see the rise in power from early times when the Lyman-$\alpha$ pumping turns on, peaking at $z\approx 25$, and the cancellation between the Lyman-$\alpha$ and X-ray terms at $z\approx 21$, followed by the increase during the X-ray heating phase, at $z\approx 18$.

V Total signal and observability

The inclusion of streaming velocities produces VAOs in the 21-cm power spectrum on large scales. We now explore how the signal changes when changing different assumptions, especially the poorly known Lyman-Werner feedback, as well as how observable these VAOs are by upcoming 21-cm interferometers.

V.1 Different Feedback Assumptions

We start by considering the VAOs under three additional feedback assumptions, different than the regular-feedback strength that we have focused on so far. The first is an unrealistic no-feedback case, which will show the largest 21-cm VAOs. The other two are realistic feedback strengths: low and high, slightly more optimistic and pessimistic than our fiducial case, as explained in Section II. We also explore the 21-cm power spectrum for the unphysical case that LW feedback is saturated at all times (so only atomic-cooling haloes can form stars) in App. B.

We first show, in Fig. 15, the evolution of the global 21-cm signal as a function of redshift for the four feedback cases. Here, barring the effect of the LW feedback, we have kept fixed all the fiducial parameters (except in the no-feedback case, where we have a lower $f_{*}=0.03$). Adding feedback always delays evolution, as argued in Section III, which makes the three realistic feedback cases delayed by $\Delta z\sim 3$ with respect to the no-feedback case. Within the three realistic cases, however, the low feedback is the fastest, due to the smaller effect of LW radiation, followed by the high-feedback case, which is affected by LW as much as the regular case, while the $v_{\rm cb}$ effect is smaller (due to the assumption of uncorrelated feedback), producing smaller halo suppression. The regular case is the most delayed.

Even though the global signal does not vary significantly between the different feedback cases, the large-scale 21-cm fluctuations will. We show the 21-cm power spectra for all feedback cases at the midpoint of the EoH in Fig. 16 (corresponding to $z=\{19.1,17.6,17.2,17.2\}$, in ascending order of feedback strength), and at the midpoint of the LCE in Fig. 17 ($z=\{27.7,25.5,24.5,25.5\}$)666We note that 21cmFAST has a redshift resolution of $\Delta z\approx 0.3$ during cosmic dawn, which determines our precision on all redshifts reported.. In all cases we see a similar small-scale $k\gtrsim 0.3\,\rm Mpc^{-1}$) power at a given redshift, while increasing the feedback lowers the large-scale power due to the VAOs. This is expected, since stronger feedback will prevent smaller-mass molecular-cooling haloes to form stars, thus rendering some of the suppression of the relative velocities ineffective. In particular, we see in Fig. 16 how the no-feedback case gives rise to large ($\Delta_{21}^{2}\gtrsim 400$ mK2) acoustic power at low-$k$, whereas that power is reduced to $\Delta_{21}^{2}\approx 200$ mK2 for the realistic feedback cases. The situation is similar during the LCE, albeit with smaller power overall. Interestingly, in all cases our curves provide a good fit (through Eq. (36)) to the simulation data at relevant scales ($k=0.03-0.3$ Mpc*-1*), showing that the effect of astrophysical feedback is to lower the amplitude of the VAO signal, but not to alter its shape, which is indeed that of $\Delta^{2}_{v^{2}}$.

Finally, we show the amplitude of the VAOs as a function of redshift in Fig. 18 for the same four feedback cases. All the curves follow a similar trend, with the VAO amplitude rising at $z\sim 25$ due to Lyman-$\alpha$ pumping, and again at $z\sim 18$ due to X-ray heating, with the transition in between having nearly vanishing velocity-induced power. Nonetheless, the size of the VAOs changes dramatically for different feedback assumptions. Focusing on the peak of each curve during the EoH, the power spectrum at a reference scale $k_{\rm ref}=0.14\,\rm Mpc^{-1}$ is reduced from 234 mK2 in the no-feedback case, to 138 mK2 for low feedback strength, and even further to 14 mK2 in the pessimistic case of high feedback. In the regular case studied in the previous section this amplitude is 64 mK2. This large range brackets the possible amplitude of the VAOs during cosmic dawn, depending on the thus far unknown physics of the first stellar formation. Nonetheless, in all feedback cases, even with pessimistic assumptions, the power spectrum shows nonzero $v_{\rm cb}$-induced 21-cm perturbations. Our results are in agreement with the overall trends in Ref. [23], albeit with some distinctions given the different choices of astrophysical and feedback parameters, which we will discuss below.

V.2 Observability

While the 21-cm signal due to the DM-baryon relative velocities is significant during both the EoH and the LCE, the large foregrounds make the task of observing—and characterizing—it challenging. We now forecast the sensitivity of large 21-cm observatories to the VAO signal, focusing on the upcoming Hydrogen Epoch of Reionization Array (HERA)777https://reionization.org  [43], for which we obtain the projected power-spectrum uncertainties with the publicly available 21cmSense [93, 94]888https://github.com/jpober/21cmSense . This code accounts for the $u$-$v$ sensitivities of each antenna in an array, and returns the errors in the 21-cm power spectrum, including cosmic variance.

In addition to the instrumental noise of the antennae, one of the main obstacles to any prospective 21-cm detection is the presence of large foregrounds, from both our Galaxy and the atmosphere, that contaminate the signal. These foregrounds swamp any primordial 21-cm signal in the majority of the Fourier plane, as wavenumbers with

[TABLE]

are part of the so-called “foreground wedge”, and not expected to be usable for cosmology [95, 96, 97, 98, 99]. Here $k_{||}$ and $k_{\perp}$ are the line-of-sight (LoS) and perpendicular wavenumbers, and $a$ and $b$ parametrize the severity of the foregrounds by accounting for how many wavenumbers are discarded. We follow the parametrization in 21cmSense, and define three foreground scenarios:

$\bullet$ Pessimistic, with $a=0.1\,h\,\rm Mpc^{-1}$, and $b$ given by the horizon limit, as usual.

$\bullet$ Moderate, with $a=0.05\,h\,\rm Mpc^{-1}$, and the same $b$.

$\bullet$ Optimistic, with $a=0$, and a smaller $b$ by a factor of $\sin(\theta_{\rm FWHM}/2)$, where $\theta_{\rm FWHM}$ is the full-width half-maximum of the beam.

Note that here we report wavenumbers multiplied by the reduced Hubble constant $h$, as those are the inputs of 21cmSense. In all cases we assume that the sensitivity of different baselines can be added coherently, and we take a system temperature of $T_{\rm sys}=100+120(\nu/150\,\rm MHz)^{-2.55}\,\rm K$, as in Ref. [43]. For reference, the wedge has been experimentally observed in Ref. [100], with $a$ and $b$ similar to the pessimistic case, although it is likely that a better understanding of the foregrounds and the instrument will improve these parameters.

We divide the observation range into four redshift bins, three of them in the EoH, centered at redshifts $z=14,16,$ and 18, each with a width $\Delta z=2$, and one in the LCE, centered at $z=25$ and with $\Delta z=4$. We use the HERA Split Core configuration, with 540 days of observation (corresponding to roughly three full years), and the standard settings of 21cmSense, albeit binning logarithmically in $k$ instead of linearly to better capture the VAO shape.

We calculate the signal-to-noise ratio (SNR) from each redshift slice summing simply over $k-$bins as

[TABLE]

where the binned errors $\delta(\Delta^{2}_{21})$ in the 21-cm power spectrum are the output of 21cmSense, and we consider the range $k=(0.05-0.5)\,h\,\rm Mpc^{-1}$. Note that we are only computing the SNR of the VAOs, although, of course, we include the full signal in $\delta(\Delta^{2}_{21})$ to properly account for cosmic variance. We show our results in Table 1 separately for the EoH and the LCE, considering each combination of foreground and feedback cases studied in this paper, where the SNRs of the EoH bins have been added in quadrature. We see that HERA has a high likelihood of observing the VAOs during the EoH, where our regular-feedback case yields SNRs of 10-62, and even the most pessimistic feedback assumptions predict a SNR larger than unity. The situation is significantly worse during the LCE, where only under highly optimistic assumptions will it be possible to detect the VAO signal at SNR larger than 5. The cases with the largest SNRs yield the most promising cosmological measurements from VAOs [44].

As a test, we have run simulations varying our fiducial X-ray luminosity by a factor of six upwards (so $\log_{10}(L_{X}/{\rm SFR})=40.8$), to determine how much this parameter affects the observability of the VAOs. We show our results for this case in Table 2. On the one hand, increasing the luminosity shifts all heating transitions earlier (for this case $\Delta z\approx 2$ away from our fiducial case, as for instance we find $z_{0}=17.7$), where feedback effects were less severe, thus increasing the (relative) size of the VAOs. On the other hand, the overall 21-cm absorption is shallower, as the gas is warmer by the time Lyman-$\alpha$ coupling is completed. We find that the velocity-induced power spectrum at its redshift peak $z=19.7$ (and scale $k_{\rm ref}=0.14\,\rm Mpc^{-1}$) is $\Delta_{\rm vel}^{2}(k_{\rm ref})=65$ mK2. This is very similar to our fiducial case, showing that the effect of weaker feedback at large redshifts is compensated by the shallower 21-cm absorption, which sets the overall scale of the signal. There is one more difference between these two cases, however. Since all heating transitions occur earlier for higher X-ray luminosity, an increased noise (due to the higher system temperature) will be present in these large-$z$ observations, which lowers the SNRs for this case in Table 2.

V.3 Comparison to Previous Work

We are not the first to compute the effect of the DM-baryon streaming velocities in the 21-cm power spectrum, although we are the first to do so with both a (public) simulation code and analytical results. This has allowed us to confirm that the well-understood shape of the VAOs is indeed present in 21-cm simulations. We will now compare our results with those of previous work on this topic.

On the theory side, the authors in Ref. [13] were the first to predict and calculate the velocity-induced 21-cm fluctuations due to Lyman-$\alpha$ pumping. In that work it was predicted that the acoustic shape of fluctuations would follow $\Delta^{2}_{v^{2}}$ with suppression on small scales, exactly as we found in our simulations, where the window function posited by Ref. [13] provides a good fit to our data. Additionally, that work showed a $v_{\rm cb}$-induced 21-cm amplitude of $\delta T_{21}\approx$ 3 mK, in line with (albeit smaller than) our result of $\delta T_{21}\approx$ 10 mK at our peak ($z=25$). Additionally, in Ref. [24] the VAOs during the EoH were analytically estimated, and were found to produce a 21-cm power spectrum of $\Delta^{2}_{21}\approx 100$ and 300 mK2 at $z=15$ and $20$, respectively. This is in line with the results of our simulations. Note, however, that in Refs. [13, 24] the non-VAO part of $\Delta_{21}^{2}$, which we parametrized as a smooth polynomial $\mathcal{P}_{n}$, was assumed to be proportional to the density power spectrum. This is an overly simplistic assumption, which does not resemble the real 21-cm signal, as is clear from the results of Refs. [29, 25, 22], or our own Fig. 11. The non-VAO part of the power spectrum behaves as a shallower function of wavenumber $k$, resulting in observable VAOs even at $k\sim 0.1$ Mpc*-1*.

On the simulation side, in Ref. [25] the effect of velocities in the 21-cm signal was first implemented, albeit only applied to the X-ray heating era. In that study it was found that the relative velocities can produce 21-cm fluctuations that dominate the signal on large scales, and have sharp BAO-like features, as we confirm here. Their results are comparable to our “no-feedback” case (purple line in Fig. 16), which, however, we argue is unlikely to represent reality. This code was improved upon in Ref. [22], where the effects of Lyman-Werner feedback were included. In this case we find broad agreement between their “strong-feedback” and our ”regular-feedback” models, as both have the same feedback parameters, although we find larger-amplitude oscillations, likely because of our homogeneous LW feedback. Finally, in Ref. [23] that work was extended to include fluctuations during the LCE, where the shape of the VAOs seems to be damped at $k\gtrsim 0.1$ Mpc*-1*, which we find here as well. We note two important differences between our work (based on 21cmvFAST) and that of Refs. [25, 22]. First, as in 21cmFAST we use the EPS formalism to compute the amount of photons emitted from galaxies some distance $R$ away by averaging densities and velocities (as described in Section II), which in Ref. [101] was found to agree well with full simulations), whereas in Refs. [25, 22] the collapsed fraction is calculated at each point and then averaged over $R$. Secondly, for computational simplicity we are only taking the average value of the Lyman-Werner flux across the whole simulation box, as opposed to computing it at every pixel, as done in Ref. [22].

We note that full hydrodynamical simulations of stellar formation, such as those of Refs. [102, 103, 104, 105, 106, 77, 107], are expected to be a closer match to reality than the quasi-numerical simulations done here (based on Refs. [29, 30]) or in Refs. [25, 22, 23, 108, 36, 109]. Nonetheless, the effects that we are studying chiefly affect haloes below the atomic-cooling threshold, with masses $M\lesssim 10^{7}\,M_{\odot}$, which are not resolved in current simulations [103]. A significant improvement in the mass resolution is required to resolve the haloes that are most suppressed by the streaming velocities, and we hope that this work stimulates the study of VAOs in those simulations.

VI Conclusions

In this work we have studied the effect of the dark matter (DM)-baryon relative velocities on the 21-cm line. These “streaming” velocities, sourced by the baryon acoustic oscillations (BAOs), become supersonic after recombination, and suppress stellar formation in molecular-cooling haloes (with masses $M\lesssim{\rm few}\times 10^{7}\,M_{\odot}$), which are expected to host the first stars during cosmic dawn. We performed quasi-numerical simulations with the 21cmvFAST code, which we make public. This allows us to properly include molecular-cooling haloes, as setting $T_{\rm cool}<10^{4}$ K in the regular 21cmFAST will result in a severe overestimation of the amount of star-forming haloes.

The fluctuations in relative velocities become imprinted onto the distribution of the first stars and, through their X-ray and UV emission, onto the 21-cm line. This generates velocity-induced acoustic oscillations (VAOs) in the 21-cm power spectrum on large scales, turning the cosmic dawn into a novel probe of acoustic physics, as the BAOs are imprinted onto the DM-baryon relative velocities on their genesis. We emphasize that the presence of these acoustic oscillations on the 21-cm line has previously been pointed out in Refs. [13, 24, 25, 22, 23]. Here we have compared, for the first time, analytic and simulation-based results, in order to isolate the shape and size of the VAOs. We have found that, to a good approximation, the VAOs follow a well-understood analytic shape, regardless of astrophysical effects such as photon propagation. We exploit this result in Ref. [44], where we propose using VAOs as a standard ruler at cosmic dawn. Interestingly, here we found that VAOs are observable with the upcoming HERA interferometer under a broad range of astrophysical assumptions. This would illuminate the physics of the cosmic-dawn era, showing that molecular-cooling haloes play an important role in stellar formation. In addition, the inclusion of streaming velocities can result in sizable non-Gaussianities, which we plan to explore in future work.

We note that the VAOs we study here are different from the “regular” (density) BAOs, as they do not originate from small wiggles in over/under-densities, but instead from $\mathcal{O}(1)$ fluctuations in the relative velocity. The density BAOs would also affect the 21-cm signal during different cosmic eras [110, 111, 112, 113, 114]. Nonetheless, density fluctuations affect the cosmic-dawn 21-cm line in a complicated fashion, hampering a precise reconstruction of the matter power spectrum, so we have not included them in the smooth $\mathcal{P}_{n}$ function in Eq. (36). Additionally, while we have included the effect of Lyman-Werner feedback, which increases the minimum mass required to form stars by dissociating molecular hydrogen, there may be other feedback processes, such as chemical [115] or other radiative feedback [116, 117, 118, 119, 120] that may increase the halo mass beyond even atomic-cooling haloes, rendering the effects of velocities negligible. Our simulations show, however, that if the VAO signal is present, it will always have the acoustic shape that we have modeled here.

To summarize, the cosmic-dawn era will provide us with a new arena to search for acoustic oscillations, in the form of VAOs. Their detection is within reach of upcoming 21-cm experiments, and will teach us a great deal about the haloes that hosted the first stars. Moreover, given the well-understood shape of the VAOs, they will allow us to separate cosmology from astrophysics, illuminating our knowledge about the end of the cosmic dark ages.

Acknowledgements.

We wish to thank Aaron Ewall-Wice, David Pinner, and Jonathan Pober for helpful discussions, as well as Yacine Ali-Haïmoud, Cora Dvorkin, Anastasia Fialkov, Bradley Greig, Marc Kamionkowski, Andrei Mesinger, and especially Ely Kovetz for comments on a previous version of this manuscript. We are also thankful to the creators of CAMB, CLASS, 21cmFAST and 21cmSense, as well as the libraries on which they depend, for making their codes publicly available. Some computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. This work was funded by a Department of Energy (DOE) Grant No. DE-SC0019018.

Appendix A Dependence on the X-ray Spectrum

In the main text we have assumed that the spectrum of emission of X-ray sources is given by a power-law with index $\alpha_{X}=1$ and a cutoff below energies $E_{0}=0.2$ keV [29, 30, 25]. Nonetheless, it is possible that the first abundant X-rays in our Universe were emitted by compact binary objects, which are expected to have a higher-energy cutoff in their spectrum, with $E_{0}\approx 0.5$ keV or even higher [121]. In that case, the mean-free path of X-ray photons will be elongated, erasing part of the small-scale power in the VAOs. We illustrated this in Fig. 4, where we considered the $F_{\rm coll}(v_{\rm cb})$ function averaged over distances of $R=20$ Mpc, typical of X-rays at 0.5 keV. While the small-scale power is erased, the majority of the velocity-induced fluctuations are on large scales, which are unaltered by the X-ray propagation.

To ensure that our results can be generalized to any X-ray spectrum, we have performed a simulation with and without velocities, considering a larger X-ray cutoff of $E_{0}=0.5$ keV, similar to that of Ref. [75]. As shown in Ref. [72], the X-ray spectral index $\alpha_{X}$ does not alter results significantly, so we do not change it, and neither do we alter the total soft-band ($E_{0}\leq 2$ keV) X-ray luminosity $L_{X}$, as we varied that parameter in the main text. In general, the change from a soft to a hard X-ray spectrum can delay the onset of the EoH significantly, and produce deeper 21-cm absorption [74]. This is not the case for $E_{0}\leq 0.5$ keV [72], since very low cutoffs result in very local X-ray heating, whereas for $E_{0}\gtrsim 1$ keV this delay can be important [74].

We recalculated the X-ray window function with the new value of $E_{0}$ using Eq. (35), confirming that the $v_{\rm cb}$ fluctuations are quickly erased for $k\gtrsim 0.2$ Mpc*-1*. We fit for the results of the new simulation (with regular feedback), which we show in Fig. 19 at the three critical redshifts during the Epoch of Heating, both with $v_{\rm cb}$ fluctuations and fixing $v_{\rm cb}=v_{\rm avg}$ (to keep the same background). First, we see that after heating is fully completed (at $z=12$) the velocities have a negligible effect on the 21-cm signal, as expected. Then, both at the transition redshift (where $\left\langle T_{21}\right\rangle=0$, corresponding to $z_{0}=14.8$) and the halfway point (defined as $\left\langle T_{21}\right\rangle=0.5\times\left\langle T_{21}\right\rangle_{\rm min}$ and occuring at $z=17.2$) we see that the effect of velocities is as large as in our fiducial case , whereas the VAOs are suppressed at $k\gtrsim 0.2$ Mpc*-1* (cf. Fig. 9). This would hamper the detection of 21-cm VAOs, as smaller scales are typically freer of foregrounds, and therefore easier to measure. Nonetheless, we see that even in this case there is a sizable effect at observable scales, which experiments like HERA could observe at a SNR of $\{3,7,50\}$ during the EoH under pessimistic, moderate, and optimistic foregrounds, respectively. This is slightly lower than in the case studied in the text, due to the shorter mean-free path of X-ray photons. As a test of the robustness of VAOs, we have attempted to model this case (with $E_{0}=0.5$ keV) with the window function of the case in the main text ($E_{0}=0.2$ keV). While this produces a slightly worse fit (indicating that $E_{0}$ can potentially be constrained by the shape of VAOs, were they to be detected), the acoustic peaks do not change locations, as expected.

Appendix B Atomic-cooling Haloes

Here we explore the effect of velocities in the 21-cm line assuming that Lyman-Werner feedback is fully saturated at all times, so that stars can only form through atomic cooling. This requires haloes with a circular velocity roughly above $V_{\rm cool}=17$ km s*-1* [57, 49, 58] (translating into a minimum virial temperature of $T_{\rm cool}=10,500$ K), which have masses of $M\gtrsim 3\times 10^{7}\,M_{\odot}$ at $z\sim 20$. For these haloes the effects of streaming velocities are not very significant, as they are formed of overdensities at scales larger than molecular-cooling haloes, and it is not expected that $V_{\rm cool}$ is directly modulated by $v_{\rm cb}$. Nonetheless, relative velocities can still leave a small imprint in the 21-cm line without molecular-cooling haloes, as we will explore here.

We perform a simulation under these assumptions, with the same fiducial parameters as in the main text. From this simulation, we showed in Fig. 5 the background evolution of $T_{21}$ when only atomic-cooling haloes are allowed, along with the cases studied in the main text. There we see that at high redshifts the atomic-only case severely underpredicts the amount of stars formed, delaying the beginning of the LCE by $\Delta z\approx 3$, and the transition to the EoH by $\Delta z\approx 1$. Nonetheless, as time evolves the LW feedback in our regular case becomes more important, so the atomic-only case transitions into emission only $\Delta z\approx 0.3$ later than our regular-feedback case. Moreover, by the end of our simulations ($z=12$) the regular-feedback case and the atomic-only are almost identical, as the LW feedback is nearly saturated, impeding star formation in molecular-cooling haloes. This showcases the need to properly include the redshift evolution of the LW feedback—plus relative velocities—to obtain the evolution of $T_{21}$ across the entire cosmic dawn.

Focusing on the atomic-only case for now, we show in Fig. 20 the (relative) difference between the 21-cm power spectrum with and without relative velocities at $z=16.5$ and $23.5$, corresponding to the mid points of the EoH and the LCE in this case, where we see that even for these heavy haloes there is a small VAO signal, at the few-percent level. In Fig. 20 we have performed simulations with the same initial conditions in the cases with and without relative velocities, so the “cosmic variance” (Poisson noise) mostly cancels when subtracting them [122], leaving a cleaner VAO signal. While the large-scale differences are due to the VAO, the small scale changes are due to the different backgrounds, and thus can potentially be mimicked by a redshift-dependent feedback.

We also compare our atomic-cooling case with the standard 21CMMC output in Fig. 21, as a cross check of our modifications to the code. Here we have set $v_{\rm cb}=0$, changed our transfer functions to match the fit from Ref. [123], and switched our halo mass function to that of Ref. [62], as those are the settings in 21CMMC [30, 29, 72, 86]. With these modifications we see that our prescription very well matches the standard 21CMMC output, as clear from Fig. 21, where both simulations share the same initial conditions. Nonetheless, in our prescription we compute $\delta_{m}(k,z)$ (and thus $F_{\rm coll}$) by evolving Eqs. (4), whereas in 21CMMC a scale-independent growth factor is assumed. This produces small background differences in $\left\langle T_{21}\right\rangle$ (of $\approx$ 1 mK), which induce a slight difference in the power spectra in Fig. 21 as well.

Appendix C Power Spectrum of Functions of $v_{\rm cb}$

In the main text we argued that the power spectrum of any function of the relative velocity is, to a good approximation, proportional to that of the quantity

[TABLE]

as first shown in Ref. [13]. Here we further strengthen their argument by directly calculating the power spectrum of a complete basis of functions, which can in principle reconstruct the power spectrum of any function of $v$.

We can write any function as

[TABLE]

where $a_{n}$ are some expansion coefficients, and we have defined

[TABLE]

so that $\left\langle\delta_{v}^{(n)}\right\rangle=0$ for all $n$. Analytic functions of $\mathbf{v_{\rm cb}}$ should only depend on even powers of $v_{\rm cb}$, although we will also consider odd values of $n$ here. The power spectrum of $f$ due to velocities will be given by

[TABLE]

so in principle it is sufficient to calculate all cross spectra

[TABLE]

to find the power spectrum of any arbitrary function $f$ of the velocities. It is this quantity $\Delta_{n,m}^{2}(k)$ that we will focus on here.

We compute, following the approach in Refs. [13, 41, 40], the cross spectrum $\Delta_{n,m}^{2}(k)$ of velocities, and show it in Fig. 22. From this figure it is apparent that all the curves are very similar, with acoustic peaks essentially at the same locations and small differences in tilts. Therefore, the power spectrum of any smoothly varying function of velocities can be well approximated by a linear combination of these very-similar curves. We have found that the size of the power spectrum of any function of $v$ is proportional to $\left.df(v)/dv\right|_{v=v_{\rm rms}}$, and only when setting this first derivative to zero we can find deviations from the $\Delta^{2}_{v^{2}}$ pattern that we have studied in this work. We conclude, therefore, that it is unlikely that the specific dependence of $T_{21}$ on $v_{\rm cb}$ can change the shape of the power spectrum of any function, whereas photon propagation throughout the Universe can have a larger impact.