Velocity Acoustic Oscillations on Cosmic Dawn 21 cm Power Spectrum as a Probe of Small-scale Density Fluctuations

We investigate the feasibility of using the velocity acoustic oscillations (VAO) features on the Cosmic Dawn 21 cm power spectrum to probe small-scale density fluctuations. In the standard cold dark matter (CDM) model, Population III stars form in minihalos and affect the 21 cm signal through Lyα and X-ray radiation. Such a process is modulated by the relative motion between dark matter and baryons, generating the VAO wiggles on the 21 cm power spectrum. In the fuzzy or warm dark matter models for which the number of minihalos is reduced, the VAO wiggles are weaker or even fully invisible. We investigate the wiggle features in the CDM with different astrophysical models and in different dark matter models. We find that (1) in the CDM model the relative streaming velocities can generate the VAO wiggles for broad ranges of parameters f *, ζ X , and f esc,LW ζ LW, though for different parameters the wiggles would appear at different redshifts and have different amplitudes. (2) For the axion model with m a ≲ 10−19 eV, the VAO wiggles are negligible. In the mixed model, the VAO signal is sensitive to the axion fraction. For example, the wiggles almost disappear when f a ≳ 10% for m a = 10−21 eV. Therefore, the VAO signal can be an effective indicator for small-scale density fluctuations and a useful probe of the nature of dark matter. The Square Kilometre Array-low with ∼2000 hr observation time has the ability to detect the VAO signal and constrain dark matter models.


INTRODUCTION
The hyperfine transition line of the neutral Hydrogen, named 21 cm line, is one of the most powerful tools that can efficiently map the 3D large-scale structures of the Universe from the recombination era to present day.This signal depends on the density of neutral Hydrogen, the spin temperature and the CMB temperature.The brightness temperature of the signal writes (e.g.Mesinger et al. 2011; Barkana & Loeb 2001;Furlanetto et al. 2006;Pritchard & Loeb 2012) 15 (1) where  HI is the neutral fraction of the hydrogen;  is the overdensity of the gas;  CMB is the CMB temperature at redshift ; Corresponding author: Bin Yue; Yan Gong; Xuelei Chen yuebin@nao.cas.cn;gongyan@nao.cas.cn;xuelei@cosmology.bao.ac.cnΩ m , Ω b and ℎ are the matter relative density and Hubble constant respectively.The spin temperature of Hydrogen atoms is determined jointly by the scattering with CMB photons and Ly photons, and the collisions between baryonic particles, where   and  c are Ly coupling (the Wouthuysen-Field effect, Wouthuysen 1952;Field 1958) and collision coupling coefficients respectively;  k is the kinetic temperature of the inter-galactic medium (IGM).The 21 cm signal depends on the Ly, the X-ray, and the ionizing radiation fields.At the Cosmic Dawn, these fields are built up mainly by the formation of Pop III stars in minihalos.Therefore the 21 cm signal reflects not only the large-scale distribution of the neutral hydrogen, but also the mechanisms that suppress or promote the formation of minihalos and the star formation therein.
After the decoupling of dark matter, and before the recombination era, the baryons are coupled tightly with the photons, while the dark matter evolves almost independently.As a result, there are large-scale relative motions between the dark matter and baryonic matter.The root mean square (rms) of the relative streaming velocities ≈ 30 km s −1 at the recombination era (Tseliakhovich & Hirata 2010).After recombination, the rms of the relative streaming velocities decays as  rms =  2 db 1/2 ≈ 30(1 + )/(1 +  rec ) km s −1 as the Universe expands, where  rec is the redshift of recombination.Because of the baryon acoustic oscillations before recombination, the relative streaming velocity was also coherently oscillating, and it imprints oscillatory wiggles on the power spectrum of the large-scale streaming velocity field after recombination (e.g.Tseliakhovich et al. 2011).
The relative streaming velocity varies in space, and it can prevent the gas from collapsing into dark matter halos with circular velocities smaller than or comparable to the local streaming velocity.Many simulations show that Pop III stars form first in minihalos with mass ∼ 10 5  ⊙ through H 2 cooling (e.g.Haiman et al. 1996;Tegmark et al. 1997;Yoshida et al. 2008).Such H 2 -cooling minihalos have typical circular velocities of ∼ 3 km s −1 , comparable with the relative velocity which has  rms ≈ 0.6 km s −1 at  = 20.The formation of Pop III stars in those minihalos located in regions with streaming velocities ≳ 1 rms would be significantly suppressed (Greif et al. 2011;Fialkov et al. 2012).The spatial distribution of the star-forming minihalos will reflect the large-scale features of the relative streaming velocity field (Tseliakhovich & Hirata 2010;Tseliakhovich et al. 2011).
The radiation fields of the Pop III stars (i.e.Ly, X-ray, and ionizing photons) should show large-scale features corresponding to the structures of the streaming velocity field, and finally this would in turn affect the large-scale distribution of the 21 cm signal (Dalal et al. 2010;Visbal et al. 2012;Fialkov et al. 2013;Ali-Haïmoud et al. 2014;Muñoz 2019a;Muñoz et al. 2020Muñoz et al. , 2022;;Schauer et al. 2023;Conaboy et al. 2022;Hegde & Furlanetto 2023).Furthermore, the dark matter/baryon relative motion also suppresses the growth of small-scale perturbations, and influences the 21 cm forest (Shimabukuro et al. 2023).
Basically, the overall shape and amplitude of the 21 cm power spectrum will be changed by the relative streaming velocities.More importantly, the large-scale coherently oscillating structures would leave imprints on the 21 cm signal -generate some wiggles on the 21 cm power spectrum (Mc-Quinn & O'Leary 2012;Ali-Haïmoud et al. 2014).Such velocity acoustic oscillations (VAO) features are useful tools for constraining the cosmology (e.g.Muñoz 2019b;Sarkar & Kovetz 2023) and astrophysical characteristics of Pop III stars.
Although the cold dark matter (CDM) model has been taken as the standard model of cosmology, the dark matter could alternatively be warm or fuzzy (e.g.Newton et al. 2021;Dekker et al. 2022;Dayal & Giri 2023;Iršič et al. 2017;Nadler et al. 2021;Marsh & Niemeyer 2019).A distinguishing character of these models is their small-scale density fluctuations.In the warm dark matter (WDM) or fuzzy dark matter (FDM) models, small-scale fluctuations are suppressed and the abundance of minihalos are significantly reduced.As a result, the build-up of Ly, X-ray, and ionizing radiation fields is delayed, and the global evolution and spatial fluctuations of the 21 cm signal change.Therefore the 21 cm signal can be a useful tool to distinguish the different dark matter models (e.g.Shao et al. 2023;Hibbard et al. 2022).In the WDM or FDM model, the number of small halos is reduced, their contribution to the ionizing photons budget is also reduced.To be consistent with the constraints on reionization history, massive halos may hold a larger amount of stellar mass.This can explain the high- stellar mass excess (Gong et al. 2023;Lin et al. 2023).Note however, since the clumpiness of the IGM is also reduced in the WDM model, the reionization process is not necessarily delayed (Yue & Chen 2012).
In addition to the relative streaming velocities, the Lyman-Werner (LW) radiation (e.g.Omukai 2001;Machacek et al. 2001;Wise & Abel 2007;O'Shea & Norman 2008) and Xray radiation (e.g.Ricotti 2016;Park et al. 2021aPark et al. ,b, 2023) ) from previously-formed Pop III stars can also suppress the Pop III stars formation in minihalos, boosting the critical mass of minihalos that can host Pop III stars.Since the influence of relative streaming velocities on the 21 cm signal is more sensitive to the smaller minihalos, the LW and Xray feedback effects may significantly reduce the influence of relative motion.Sarkar et al. (2022) investigated the 21 cm global spectrum and power spectrum in the FDM model, in the presence of dark matter/baryon relative motion and LW feedback, and the CMB and Ly heating.They predicted that HERA should be able to constrain the FDM masses up to ∼ 10 −19 − 10 −18 eV.Flitter & Kovetz (2022) further proposed that HERA can detect the FDM down to a fraction of ∼ 1% in the mixed dark matter model.In Vanzan et al. (2023) they investigated the feasibility of using the 21 cm angular power spectrum, influenced by both small-scale suppression and the relative motion, to constrain the ultra-light axions for various array configurations.
The large-scale VAO wiggles on the 21 cm power spectrum are modulated by the formation of Pop III stars in minihalos, whose abundance is rather sensitive to the small-scale fluctuations of the density field.Such wiggles, regardless of the overall shape and amplitude of the 21 cm power spectrum, can be a good indicator for small-scale density fluctuations.Hotinli et al. (2022) proposed that the large-scale VAO wig-gles on the 21 cm power spectrum can be used to measure the mass of ultra-light axions, and predicted that the HERA experiment can detect such signal when the axion mass  a ≳ 10 −18 eV, much more sensitive than current constraints.
The LW feedback is the main negative mechanism that reduces the amplitude of the VAO wiggles.Sometimes it even makes the VAO signal less efficient in distinguishing the CDM model with strong feedback, and the WDM/FDM model.The effect depends on the properties of Pop III stars.Currently, there are only indirect constraints on the properties of Pop III stars (e.g.Xing et al. 2023;Maiolino et al. 2023;Mebane et al. 2020;Chatterjee et al. 2020;Bevins et al. 2022).It is therefore worthwhile to consider the evolution and amplitude of the VAO signal for wide ranges of Pop III star parameters, and investigate the feasibility of distinguishing CDM and non-CDM dark matter models (e.g. the FDM or WDM), or constraining the fraction of the non-CDM component in mixed dark matter models, by the upcoming large radio interferometer arrays such as the SKA-low telescope.This is the basic motivation of this work.
The outline of this paper is: in Sec. 2 we introduce the methods to include the streaming velocivites and LW feedback effects in calculating the 21 cm signal from the Cosmic Dawn.In Sec. 3 we present our results of the VAO signal on the 21 cm power spectrum for various astrophysical parameters and axion masses.In Sec. 4 we give the summary and conclusion.

The collapse fraction field modulated by relative streaming velocities
The power spectrum of dark matter (d) or baryon (b) where  pri =      is the power spectrum of the primordial overdensity  pri ;  d/b (, ) is the transfer function of dark matter or baryon at redshift .Throughout this paper, we compute the transfer functions using the AxionCAMB program (Hlozek et al. 2015), which is a modified version of the cosmology code CAMB (Lewis et al. 2000).
According to the linearized continuity equation the velocity of dark matter (v d ) or baryon (v b ) in Fourier space is where  = 1/(1 + ) is the scale factor.To make predictions of the 21cm signal before reionization, we modify the original version of 21cmFAST (Mesinger et al. 2011) to include the effect of relative motion between dark matter and baryon.
The velocity fields of both the dark matter and the baryon is obtained according to Eq. ( 5), then the speed of the relative motion v db (x, ) = v d (x, ) − v b (x, ).We generate random density and velocity fields on 120 3 cells in a box with comoving side length 2000 Mpc.By comparing the results of simulations with the same box size but different resolutions, we find that such a choice is enough for capturing the converged large-scale VAO wiggles, see the Sec.4.2.2.
The streaming velocities between dark matter and baryon could prevent the IGM gas from being accreted into smaller minihalos, so that Pop III stars can never form in such halos (e.g.Fialkov et al. 2012).Moreover, the LW radiation also suppresses the star formation in smaller minihalos and reduces the VAO wiggles on the 21 cm signal (e.g.Fialkov et al. 2013).Schauer et al. (2021) and Kulkarni et al. (2021) first performed high-resolution numerical simulations to find the critical minihalo mass for Pop III stars formation in the presence of LW feedback and relative streaming velocities simultaneously.Since Kulkarni et al. (2021) explored the LW specific intensity range wider than Schauer et al. (2021), we adopt their critical halo mass as the fiducial model.We will discuss the results using Schauer et al. (2021)  , (6) where  LW,21 is the specific intensity of LW radiation in units of 10 −21 erg s −1 cm −2 Hz −1 sr −1 ,  rec db is the dark matter-baryon streaming velocity at recombination era. =20 and  are fitted well by and with parameters: ( =20 ) 0 = 1.96 × 10 5  ⊙ ,  1 = 0.80,  2 = 1.83,  3 = −0.06,and  0 = 1.64,  1 = 0.36,  2 = −0.62, 3 = 0.13.The halo collapse fraction for a cell with overdensity  cell , mass  cell and streaming velocity  db is where ρm is the mean density of the total matter in the Universe.The conditional halo mass function / in the cell is calculated from the Extended Press-Schechter formalism (Lacey & Cole 1993;Mo & White 1996;Sheth & Tormen 2002), where  lin () is the linear growth factor and  is the variance of density fluctuations at mass scale .In the standard CDM model the collapse barrier is simply mass-independent,   () = 1.686/ lin ().In the non-CDM model, however, the collapse barrier depends on mass.We follow the Marsh & Silk (2014) to define a mass-dependent critical overdensity for the non-CDM (axion or axion-CDM mixed) model, with the growth factor where   = /[3/(4)] 1/3 ,  0 = 0.002ℎ Mpc −1 and  ℎ = 300.The transfer function  non−CDM is calculated by axionCAMB.
Different from the original steps of 21cmFAST where the density field is first smoothed and then the collapse fraction is calculated by the smoothed density, here we first calculate the collapse fraction for each cell of our random field by the EPS formalism (Lacey & Cole 1993;Mo & White 1996;Sheth & Tormen 2002), then we smooth the collapse fraction field on different scales directly by Fast Fourier Transform (FFT), to calculate the flux of X-ray and Ly photons from Pop III stars.
The specific intensity of the LW radiation at redshift   LW () = where  esc,LW is the escape fraction of the LW radiation,  up = 13.6/11.2(1+ ) − 1,  () is the Hubble parameter and  is the light speed.The LW emissivity is given by where   is the number density of baryon atoms in the present Universe,  LW is the total released LW energy per stellar atom throughout the stellar lifetime, divided by the LW bandwidth.21cmFAST adopted the Pop III star spectrum in Barkana & Loeb (2005), the corresponding  BL05 LW ∼ 1.1 × 10 −21 erg/Hz.But we will explore the VAO signal for a wide range of  LW values.Finally, we calculate the LW intensity of each cell using (15) Regarding the  esc,LW , it depends on both stellar mass and host halo mass (Schauer et al. 2015(Schauer et al. , 2017)).The result of the moderate model in Schauer et al. (2015) is  esc,LW ∼ 0.5.In this paper, we treat  esc,LW  LW as a joint parameter.
In addition to the LW radiation, the X-ray radiation also influences the formation of Pop III stars.Different from the LW radiation that almost purely plays a negative role, the effect of X-ray radiation is more complicated.Strong X-ray radiation will, of course, suppress the formation of Pop III stars, since it heats and ionizes the IGM.However, moderate X-ray radiation may promote the formation because it generates lots of free electrons that help H 2 to form (e.g.Haiman et al. 2000;Ricotti 2016;Park et al. 2021a,b, 2023 andreferences therein).For this reason, we will not model Xray feedback on Pop III star formation in our paper.Instead, we limit our investigation in the redshift range where the Xray feedback is not the dominant effect, say  k ≲ 10 3 K.In all our cases, when the VAO wiggles reach maximum amplitude,  k is still smaller than or comparable to the CMB temperature.We therefore expect the X-ray feedback will not change our conclusion on the VAO signal.
In Fig. 1 we plot slices of the velocity field |v db (x, )| and density field Δ(x, ) at  = 20, for the CDM model.The results of axion models are quite similar, as they have quite similar power spectra on large scales.We also plot the collapse fraction  coll (>  crit , x, ) fields of the same slice for the CDM model and axion model with  a = 10 −20 eV, ignoring the LW feedback.The collapse fraction field is modulated by both the density field and the relative velocity field, so in principle, it should be tightly correlated with both of them.However, the latter mechanism relies on the contribution of halos with mass ∼  crit in  coll .For the axion model, small halos are suppressed, therefore the VAO feature in the collapse fraction field is much less obvious than in the CDM model.In the CDM model, the correlation coefficient between |v db (x, )| field and  coll (x, ) field is −0.44, and that between Δ(x, ) field and  coll (x, ) field is 0.84.In the axion model with  a = 10 −20 eV, the correlation coefficients are -0.03 and 0.82 respectively.Indeed, for such an axion model, the collapse fractions are almost independent of the relative streaming velocities.The rms of noise flux per beam for an interferometer array with  antennas is

The radio telescope sensitivity
where  B is the Boltzmann constant,  sys is the system temperature,  e is the effective area of either a single dish or the full array,  is the frequency channel width and  s is the on-source integration time.Using the total observation time where Ω survey is the survey area and Ω FoV is the field of view of the array, and the survey speed figure of merit (SSFOM) (Braun et al. 2019), the noise is The noise power spectrum where  beam = / max is the beam width,  max is the maximum baseline length, and is the comoving volume of the real space voxel. () is the comoving distance up to .Throughout this paper, we adopt  max = 2000 m, as most of the SKA-low stations are located within such distance1.
The variance of the measured 21 cm power spectrum is where   () is the number of -modes in the -bin; the window function (Battye et al. 2013) denotes the rapid decline of the measured power spectrum below the resolution.Suppose a survey has area Ω survey and bandwidth Δ, then the survey volume and where  ln  is the relative width of the selected -bin and we take  ln  = 0.12.

The modulations of spin temperature by Ly𝛼 scattering and X-ray heating
Since the VAO wiggles reflect the inhomogeneity of the Ly radiation and X-ray heating, according to Eq. ( 1) and the formation history of Pop III stars, in principle, we expect to see the VAO signal in the following scenarios: 1.If the Ly coupling effect starts to work while the IGM is still neutral and unheated, say   ∼ 1,  k ∼  ad k and  HI ∼ 1, where  ad k is the kinetic temperature of the IGM following adiabatic evolution.In this case, the 21 cm power spectrum reflects the structures of the Ly radiation field.This is the pure Ly-modulated mode; 2. If the Ly coupling has already been saturated, the IGM is obviously (but not heavily) heated but still not yet reionized, say   ≫ 1,  ad k <  k ∼ O ( CMB ) and  HI ∼ 1.In this case, the 21 cm power spectrum reflects the structures of the X-ray radiation field.This is the pure X-ray-modulated mode; 3. If both the Ly coupling and IGM heating have already been saturated, and the IGM has already been obviously (but not fully) reionized, say   ≫ 1,  k ≫  CMB and  HI < 1.In this case, the 21 cm power spectrum reflects the structures of the ionizing radiation field.This is the pure reionization-modulated mode.
If these modulations work at distinct and well separated epochs, we would observe the following sequence: at the beginning when both the Ly and X-ray radiation are weak, there is no VAO signal.Then the VAO signal appears when Ly modulation starts to work.This signal disappears later on when Ly coupling saturates, i.e.   ≫ 1.The VAO signal will appear again when the X-ray modulation works.And disappears again when  k ≫  CMB .Next, the VAO signal may appear again when the reionization modulation works, and finally disappears at the end of the reionization era.In practice, however, these three modulation mechanisms may work together, showing joint effects, depending on the choice of parameters of star formation efficiency  * , X-ray production efficiency   that refers to the cumulative number of X-ray photons produced per Solar mass, and the ionization efficiency  ion .
The 21 cm signal is modulated by star formation in minihalos.When the IGM is heated to  k ≳ 10 3 K, the Jeans mass (Ripamonti et al. 2007) exceeds the atomic cooling mass, and star formation in minihalos is fully suppressed and the VAO signal disappears soon.In this paper, we limit our investigation in the stage with  k ≲ 10 3 K.In Fig. 2 we plot the 21 cm global spectrum in different dark matter models.Here we set  * = 0.005 since it is generally believed that the star formation is less efficient in minihalos (e.g.Trenti & Stiavelli 2009).We adopt  X = 5 × 10 55  −1 ⊙ , roughly corresponding to  X = 0.05 in Furlanetto et al. (2006).It is seen that, for such star formation efficiency and X-ray efficiency, the 21 cm global signal achieves the maximum absorption ∼ −100 mK at  ∼ 20 in the CDM model.In the next paragraph, we will see that the VAO wiggles reach the maximum amplitude at  ∼ 17.We note that, however, the value of   for Pop III stars is still quite uncertain.For example, the SARAS-2 result (Singh et al. 2017) disfavors models with very low Xray efficiency.In practice, generally, theoretical models with X-ray efficiency spanning a wide range are allowed (Cohen et al. 2017).Moreover, we set  esc,LW  LW =  BL05 LW .For the axion model with  a = 10 −18 eV, the 21 cm global spectrum is almost identical to the CDM.However, for  a ≲ 10 −19 eV, the global evolution of the 21 cm signal is delayed.
In the top panel of Fig. 3 we plot the 21 cm power spectrum at redshifts when   ∼ 1.0, for the CDM model and axion models with different masses.It is clearly seen that in the CDM model and models with massive axion particles (i.e. a ≳ 10 −18 eV), the streaming velocities produce wiggles on the 21 cm power spectrum, through both Ly scattering We also check the detectability by using the uncertainties estimated from Eq. ( 20).We ignore the foreground here but present discussions of its influence in Sec.4.2.3.The VAO signal in the CDM model is marginally detectable for SKA1-low with Ω survey = 20 deg 2 and  obs = 2000 hour.The total signal-to-noise ratio is ≈ 6. See the hatched regions in the bottom panel of Fig. 3. Since the VAO signal appears at large scale, the cosmic variance dominates the uncertainties when the survey is as small as 20 deg 2 .However, for SKA1low this cannot be reduced by simply increasing the survey area, because given the total observation time  obs , if the survey area is larger, then the on-source integration time is smaller.So larger area survey may have smaller uncertainties at large scale where the uncertainties are dominated by sample variance, but will have larger uncertainties at small scale where the uncertainties are dominated by the instrumental noise.For this reason, we consider the SKA2-low that has higher sensitivity, and Ω survey = 200 deg 2 .See the filled 2 One can also fit the wiggles by the analytical formula proposed by Muñoz (2019a).regions in the bottom panel of Fig. 3.For such a survey, the total signal-to-noise increases to ≈ 22.
In both CDM and axion models, there are relative motions between the dark matter and baryon.However, in the axion model, or any other dark matter lacking small-scale density fluctuations, there is a lack of minihalos and Pop III stars, therefore the relative motions have no chance to modulate the Ly radiation field and X-ray field.Therefore the VAO features are not encoded in the 21 cm signal.Indeed, the VAO signal is a good indicator of the small-scale density fluctuations.
3.2.The choice of  * ,   and  esc,LW  LW For the CDM model, in Fig. 4, 5, and 6 we show the extracted VAO wiggles vs.  and  for varying  * ,   and  esc,LW  LW respectively.Parameters are given in each panel.
We see that the VAO signal increases with increasing  * , and decreases with increasing LW feedback.Interestingly, for   the behaviors are more complicated.When   ≲ 1 × 10 58  −1 ⊙ , the VAO signal decreases with increasing X-ray radiation.This is because as the IGM is heated, the spin temperature is more close to the CMB temperature, so the 21 cm signal becomes weaker.However, for higher   the VAO signal starts to increase again.This is because the X-ray is also a strong Ly source.Roughly one-third of the X-ray energy is finally deposited into Ly photons.Therefore X-ray can generate not only inhomogeneous heating but also an inhomogeneous Ly background that promotes the VAO signal.Moreover, for higher   such a mechanism works earlier.This also promotes the VAO signal since at higher redshifts the relative streaming velocities are larger.Also for this reason, increasing the   will shift the VAO signal to higher redshifts, but cannot wipe out it.Note that here we show the results for   ≳ 1 × 10 58  −1 ⊙ just for the purpose of exploring the parameter space that allows the existence of the VAO wiggles.It is not clear if Pop III stars can really produce such amount of X-ray photons.
Even though a very strong LW radiation field is artificially set,  esc,LW  LW is 10 times  BL05 LW , there are still VAO features on the 21 cm power spectrum (the bottom right panel in Fig. 6).This is because: first, the LW radiation field is not just simply proportional to the  esc,LW  LW , since stronger LW radiation will suppress the star formation in smaller halos, reducing the production amount of the LW radiation in return.Second, in the Kulkarni et al. (2021) model, the critical dark matter halo mass for H 2 cooling is less sensitive to the LW radiation at high redshifts.For example, at  ∼ 25, the critical mass is still just ∼ 10 6  ⊙ even though  LW,21 is as high as ∼ 30.
Nevertheless, for the CDM model, there is VAO signal for all  * ,  X and  esc,LW  LW parameters shown there.But if the signal appears too early, e.g. * = 0.005 and  X ≳ 1 × 10 59  −1 ⊙ , all the frequencies are below the lower limit of the SKA-low band.Therefore, to investigate the detectability of the VAO signal, it is important to understand the full evolution of the VAO signal over different redshifts.
For larger  * , the Ly coupling and X-ray heating occur earlier, therefore the VAO signal is more obvious, as the streaming velocities decline with decreasing redshift.The redshift at which the VAO signal first peaked can be a useful indicator of early star formation efficiency.In Fig. 7 we plot the redshift at which the first peak at  = 0.05 Mpc −1 reaches maximum amplitude, as a function of  * , see the left -axis.We also plot the amplitude of the first peak, see the right -axis.
On the other hand, reionization also sets limits on the  * , since  ion =  *  esc,ion  ion 1 1+ nrec .For Pop III stars the total released ionizing photons per stellar atom  ion = 4.4 × 10 4 (Barkana & Loeb 2005), the escape fraction of ionizing photons  esc,ion ∼ 20 − 40% (Kimm et al. 2017), and the mean recombination times per ionized atom nrec ∼ 5 − 10 ( Mao et al. 2020).Currently, the contribution to reionization from Pop III stars is not directly constrained, but it is generally believed that this is a small fraction, with the bulk part of the reionization produced by photons from the second generation stars (Pop II).In our model, for all  * and   values presented in this paper, the ionized fraction ≲ 10% when the VAO signal reaches the maximum.
In Fig. 8 we plot the extracted VAO wiggles for various axion masses.We see that when  a = 10 −18 eV, the VAO wiggles are still as obvious as that of the CDM model.However, for  a ≲ 10 −19 eV, the VAO signal is negligible at any redshifts, no matter what  * value we choose.The light "wiggles" in the redshift spectra plot for  a ≲ 10 −19 eV are produced by pure density fields, i.e. the Baryon Acoustic Oscillations (BAO) effect.The  locations of peaks and troughs of the BAO wiggles are similar (but not the same) to the VAO wiggles, however, the amplitudes are much smaller.The relative height of the BAO peak at  =0.05 Mpc −1 is ≲ 5%, while for the VAO signal the peak at  =0.05 Mpc −1 can be as high as ∼ 20% in the CDM model.
In summary, in the CDM model, we always see the VAO features for wide ranges of astrophysical parameters.For some parameters, for example  * ∼ 0.005 and   ∼ 1 × 10 59  −1 ⊙ , the frequencies of the VAO wiggles are below the lower bound of the SKA-low band.For stronger LW feedback, for example  esc,LW  LW ≫ 10 BL05 LW , Pop III star formation in minihalos with ≲ 10 6  ⊙ could be fully suppressed and the VAO signal will be reduced significantly.In such cases, finding reasonable constraints on these parameters from other observations is important.For example, the 21 cm global spectrum can put constraints on the thermal history of the Universe that is tightly related to the Pop III stars (Mirocha et al. 2018).Moreover, the unresolved Cosmic X-ray background and Cosmic Near-Infrared background can also put upper limits on Pop III stars properties (Xu et al. 2016), although in practice it is challenging to distinguish the contribution of Pop III stars and the contributions of first Pop II galaxies and low- faint foreground galaxies (Helgason et al. 2014).Nevertheless, for most literatures, the adopted parameters  * ≪ 0.1,   ≪ 1 × 10 59  −1 ⊙ and  esc,LW  LW ≲ 10 BL05 LW (e.g.Gilmore 2012; Gessey-Jones et al. 2022;Chantavat et al. 2023;Bera et al. 2023;Sun et al. 2021;Qin et al. 2021;Incatasciato et al. 2023;Ventura et al. 2023;Mittal & Kulkarni 2022;Hegde & Furlanetto 2023).That is to say, for parameters favored by most people, it is quite likely that the VAO signal does exist and its frequencies are in the SKA-low band.
On the other hand, for the dark matter model that lacks small-scale density fluctuations and minihalos, e.g.axion models with  a ≲ 10 −19 eV, the VAO signal is negligible  in all the frequency range.So indeed the VAO features are useful tools as a probe of small-scale density fluctuations.

The critical minihalo mass for Pop III stars formation
So far we use the Kulkarni et al. (2021) critical mass (Eq.( 6) ) as the threshold for Pop III star formation in minihalos in the presence of both LW radiation and relative streaming velocities.Here we check if the VAO signal still exists when we replace Eq. ( 6) with other definitions for the critical mass.
They only performed simulations with  LW,21 ≤ 0.1.Since our LW radiation field is inhomogeneous, we will use the above equations anyway.To avoid the risk, we only investigate the models with  esc,LW  LW ≤  BL05 LW here.In Fig. 9 we show the extracted VAO wiggles for  esc,LW  LW = 0.01 BL05 LW , 0.1 BL05 LW and  BL05 LW respectively, using Eq. ( 25) as the critical minihalo mass for Pop III star formation.Obviously, the amplitude of the wiggles decreases with increasing LW feedback strength.For  esc,LW  LW = 0.01 BL05 LW , the maximum peak is ∼ 10 mK.However, when  esc,LW  LW =  BL05 LW , the maximum peak is ∼ 2 mK, comparable to the SKA2-low uncertainties level with Ω survey = 200 deg 2 .So for the Schauer et al. (2021) critical mass, when  esc,LW  LW ≳  BL05 LW , the VAO features on the 21 cm power spectrum are not detectable, unless increasing the observation time.Fialkov et al. (2012) proposed that in the presence of relative streaming velocities however the LW radiation is negligible, the critical mass  crit ( db , | LW,21 = 0) is the virial mass corresponding to the circular velocity where  circ ( db = 0) = 3.714 km s −1 , and   = 4.015 for the optimal fit.When the LW feedback works, analogous to Machacek et al. (2001), the above critical mass is boosted by a factor 1 + 6.96(4 LW,21 ) 0.47 and the final critical mass is (Fialkov et al. 2013) Using the above critical mass, we show the extracted VAO wiggles for  esc,LW  LW =  BL05 LW , 10 BL05 LW , and 100 BL05 LW respectively in Fig. 10.The VAO features still exist even though the LW radiation field is strong, say  esc,LW  LW ∼ 10 BL05 LW .
This is partially because the Ly radiation field is not just influenced by the instantaneous LW feedback.Actually, even for extremely strong LW radiation, say  LW,21 ∼ O (10), which suppresses the formation of Pop III stars in minihalos, there is still a large fraction of the Ly photons from Pop III stars formed earlier when the LW feedback is still weak.However, when  esc,LW  LW ≳ 10 BL05 LW , the VAO signal is reduced significantly and not detectable for the SKA2-low.
Since Eq. ( 29) just simply re-scales the critical mass for zero-LW radiation, the relation crit ( db = 0,  LW,21 ≠ 0) (30) always holds, that the dependence of the critical mass on the relative streaming velocities is not influenced by the LW radiation.It may underestimate the role played by the LW feedback.
We can also define an alternative critical mass that is more sensitive to the LW feedback.From the simulations in Machacek et al. (2001), in the presence of LW radiation (31) The mass is then translated into circular velocity  circ,0 ( LW,21 | db = 0).Then, analogous to Eq. ( 28), the circular velocity for the critical mass in the presence of both LW radiation and relative streaming velocities is Then the viral mass corresponding to the above circular velocity is the new critical mass.
This critical mass is more sensitive to the LW radiation compared with Eq. ( 29).The extracted VAO wiggles for  esc,LW  LW = 0.1 BL05 LW ,  BL05 LW and 10 BL05 LW are shown in Fig. 11.We see that for this new critical mass, the VAO wiggles are reduced significantly when  esc,LW  LW ≳  BL05 LW .As a summary, in Fig. 12 we show the four critical masses as a function of  LW,21 for various relative streaming velocities.Their behaviors are quite different.For example, the Kulkarni et al. (2021) and Fialkov et al. (2013)  masses are sensitive to the relative streaming velocities even the  LW,21 is as high as ∼ 1.So the models using those critical masses will predict stronger VAO wiggles.The Schauer et al. ( 2021) critical mass is sensitive to the relative streaming velocities when the  LW,21 ≲ 0.1.So for this critical mass, the VAO wiggles are strong for weak LW radiation.For the critical mass corresponding to Eq. ( 32), it is sensitive to the relative streaming velocities when  LW,21 ≲ 0.01.For higher LW specific intensity, the dependence on relative streaming velocities decreases gradually (curves for different velocities are approaching each other).As a result, the VAO wiggles are reduced by the LW feedback.
To correctly calculate the VAO wiggles, it is crucial to model the critical mass appropriately.Currently, the effects of LW feedback are still not yet fully understood, particularly the self-shielding effect.Nevertheless, the Kulkarni et al. (2021) and Schauer et al. (2021) critical masses are the first results from simulations including both LW feedback and relative streaming velocities.Based on our results using these critical masses, it is likely that for broad ranges of astrophys-ical parameters, there are VAO features on the 21 cm power spectrum.

The VAO signal in various dark matter models
The locations of the wiggles are determined by the largescale distribution of the relative velocity field, so in the axion model their amplitudes are smaller but the locations are the same as that of the CDM model.In Fig. 13 we plot the VAO wiggles in different dark matter models, with the SKA1-low uncertainties for survey area 20 deg 2 and SKA2-low uncertainties for survey area 200 deg 2 respectively.
In principle, to make the forecast for the parameter constraints, we should set all astrophysical and cosmological parameters as free parameters, and involve the 21cmFAST in the MCMC procedure (Greig & Mesinger 2015).It will require quite a lot of computational cost.This is not the focus of this paper.In this paper, we focus on inspecting the existence of VAO wiggles in the CDM model with different astrophysical parameters, and the axion models with different masses.However, if we only focus on two free parameters For SKA1-low, the main constraints are from 0.01 Mpc −1 ≲  ≲ 0.1 Mpc −1 .At smaller scale the instrumental noise is large, at larger scale the cosmic variance is large.In the bottom panel of Fig. 3 and in Fig. 13, seemingly at smaller  the uncertainties are smaller, however, the relative uncertainties increase rapidly with decreasing , due to the cosmic variance.In such constraints, we do not involve the influence of the foreground, so this is an optimistic forecast.The foreground may bias the measured large-scale power spectrum, however, it can be somewhat corrected for, see the discussions in Sec.4.2.3.There is a degeneracy between the  a and  * in the SKA1-low constraints, since both lower  * and lower  a can result in weaker VAO wiggles.Nevertheless, assuming input parameter  a = 10 −18 eV and  * = 0.005, the SKA1-low can help to rule out the axion models with  a ≲ 10 −19 eV and  * ≲ 5 × 10 −3 , or with  a ≳ 10 −19 eV and  * ≳ 10 −2 .The SKA2-low can break the degeneracy and tightly constrain the  * and lower limit of  a .We remind again that this is a demonstration because we freeze the parameters   ,  esc,LW  LW and the cosmology, however it is still useful.
For mixed dark matter models, in Fig. 15 we show the VAO wiggles in dark matter models with CDM and different fractions of various axions.For axion models with  a = 10 −21 , 10 −20 and 10 −19 eV, the critical   above which VAO signal disappears are ∼ 10%, ∼ 20%, ∼ 60% respectively.For  a = 10 −18 eV there is always VAO signal even when  a reaches up 100%.Similar to Fig. 14, if we use fixed values for  * ,   and  esc,LW  LW , the predicted constraints on  a −  a is shown in Fig. 16.If in different dark matter models the VAO wiggles have the similar shape, only the normalizations are different, then the amplitude of the strongest peak at  = 0.05 Mpc −1 is a very useful indicator for dark matter properties.In Fig. 17 we show this amplitude vs. dark matter models.It shows more straightforwardly that the VAO signal almost disappears when  a ≲ 10 −19 eV.In Fig. 18 we show the amplitude of the strongest peak vs.  a , for different axion masses.To avoid a crowded panel, we use two panels.
For   = 5×10 55  −1 ⊙ and  esc,LW  LW =  BL05 LW , we provide an approximation for the amplitude of the strongest VAO peak,  2 peak = Eq. ( 33) can be used to quickly estimate the VAO signal.

The conclusion
The relative streaming velocities between dark matter particles and baryon particles have large-scale coherent structures.Such velocity can suppress the formation of Pop III stars in minihalos, reducing their production of Ly, X-ray, and ionizing photons, therefore modifying the 21 cm brightness temperature field.Therefore the structures of the streaming velocities are encoded in the 21 cm signal.Intuitively, the streaming velocities generate the VAO wiggles on the 21 cm power spectrum of the Cosmic Dawn.We modified the 21cmFAST code to add the streaming velocities and the LW feedback, then used it to investigate the VAO signal in CDM and axion dark matter models, and the axion-CDM mixed model.We found: • In the CDM model, there are always VAO wiggles on the 21 cm power spectrum for wide ranges of Pop III star parameters.The wiggles are detectable for SKA2low with ∼ 2000 hour integration time; • In the axion model, when  a ≲ 10 −19 eV, the VAO wiggles become negligible; • In the mixed models with  a = 10 −21 , 10 −20 and 10 −19 eV, when  a ≳ 10%, 20% and 60%, the VAO wiggles are neglibilbe.For  a = 10 −18 eV, there are always obvious VAO wiggles even when  a is 100%; • Using the SKA2-low, the VAO signal can be used as a good tool to distinguish dark matter models with different particle properties.
Considering the uncertainties of astrophysical parameters, to detect the small-scale density fluctuations and distinguish the CDM and the axion models reliably, for the SKA-low it is necessary to survey the full frequency range 50 MHz ≲  0 ≲ 130 MHz.

The influence of X-ray
We did not model the X-ray feedback in our paper, as this effect cannot be simply described by the change of critical minihalo mass.We checked that, however, in all our cases, when then VAO wiggles reach the maximum amplitude, the kinetic temperature of the IGM is still smaller than or comparable to the CMB.Furthermore, at this time the Jeans mass is boosted by ≲ 10, and the filtering mass (Gnedin 2000) is boosted by ≲ 1.3 by X-ray heating.Since both the Jeans mass and the filtering mass are just the minimal mass for minihalos that can hold gas, they are much smaller than the critical mass for Pop III star formation (e.g.Barkana & Loeb 2001;Glover 2013).If they are just boosted by modest fac- tors, star formation in minihalos above the critical mass is unlikely influenced.We believe that the X-ray feedback will not change our conclusion.However, to fully model such an effect, future hydrodynamic simulations including both the relative streaming velocities, the LW feedback, and the X-ray feedback simultaneously will be required.

The convergence of simulations with different resolutions
In this paper, all the fiducial simulations have the cell size 16.7 Mpc.This is much larger than the coherence scale of the relative streaming velocity field, which is ∼3 Mpc (e.g.Tseliakhovich & Hirata 2010).Since the initial conditions of the density field and velocity field are generated in Fourier space, we can reproduce their structures above the cell size, but lose the structures within the cell.When calculating the collapse fraction, it implies that for each cell we replace ⟨  coll ( db )⟩ cell with the approximation  coll (⟨ db ⟩ cell ), where ⟨⟩ cell denotes the average inside the cell.Although such a kind of approximation is quite popular in practice, it may result in some errors.
In this section, we check for the convergence of our simulations for various resolutions.In Fig. 19, we show the VAO wiggles at the same redshift ≈ 17 in simulations with cell size 16.7 Mpc (our fiducial simulation), 5 Mpc, and 2.9 Mpc, respectively.All simulations have the same parameters except the resolutions.From this figure, we see that in the overlapped  range and  ≳ 0.03 Mpc −1 , all simulations, including our fiducial simulation, have results consistent with the simulation with cell size down to the coherence scale of the relative streaming velocity field, ≈ 3 Mpc.At smaller , however, the discrepances become more obvious.

The impact of the foreground
To evaluate the impact of the foreground on the VAO signal, we generate mock 21 cm data cubes with survey area Ω survey = Figure 15.The VAO wiggles in dark matter models where CDM mixed with different fractions of axions with various masses.The hatched regions refer to the uncertainties for SKA1-low observation with survey area 20 deg 2 , while the color-filled regions refer to SKA2-low with 200 deg 2 .100 deg 2 , beam size  beam = 390 ′′ and frequency channel width  = 0.2 MHz.The 21 cm maps are centered at  = 17, when the VAO wiggles are strongest in our fiducial model.We generate two samples with bandwidth Δ = 20 MHz and 40 MHz respectively.To focus on the performance of foreground removal, we ignore the signal evolution within the bandwidth.We also ignore the beam effects.We take the foreground from the GSM sky model3 (de Oliveira-Costa et al. 2008;Zheng et al. 2017) near the Locakman hole where the Galactic foreground is lowest, and add it to the 21 cm cubes.We also add the thermal noise of SKA2-low for  obs = 2000 hour.
Since the foreground overwhelmingly dominates over the 21 cm signal and is spectrally smooth, we use the SVD algorithm (e.g.Yue et al. 2015) to find the principal components of the spectrum and remove the first 3 strongest principal components from each line of sight.The original 21 cm power spectrum and the power spectrum of the foreground-removed data cubes, and the corresponding VAO wiggles, are plotted in Fig. 20.We find that, the SVD foreground removal algorithm generates some bias (the first peak is systematically underestimated) on the recovered VAO wiggles, and the performance depends on bandwidth.For larger bandwidth, the bias is smaller.However, we realize that the foreground removal algorithms are complicated and still in development.For example, although the PCA method may underestimate the power spectrum, the power loss can be corrected by using  a transfer function (e.g.Cunnington et al. 2023).Machine learning may also help to compensate for the power loss (e.g.Makinen et al. 2021;Zhou et al. 2023).On the other hand, however, the foreground may be not that smooth in frequency space, because of the polarization leakage (e.g.Jelić et al. 2010;Gao et al. 2023).For this reason, in this paper, we only involve the instrumental noise and cosmic variance when making the forecast for constraints, as optimistic estimates.
We leave the impact of the foreground to further explicit investigations.
Another effect is the foreground wedge.Because of the chromatic instrument response, the foreground leakage contaminates the cosmological signal with Fourier modes in a wedge in the ( ⊥ ,  ∥ ) space defined as (Datta et al. 2010;Morales et al. 2012;Liu et al. 2014a,b;Jensen et al. 2016 and references therein) where  ⊥ and  ∥ are Fourier wavenumbers perpendicular and parallel to the line of sight respectively. FoV is the field of view angular size.  = / 0 and   is the transverse comoving distance, where  is the light speed and  0 is the Hubble constant. () = √︁ Ω  (1 + ) 3 + Ω Λ .Although in principle one can carefully perform foreground removal to reduce the foreground contamination in the wedge, the simplest way is just to discard all modes in the wedge.Foreground removal is necessary because the first peak of the VAO wiggles appears at  ∼ 0.05 Mpc −1 , however, the foreground dominates over the cosmological signal at least for  ∥ ≲ 0.1 Mpc −1 , see e.g.Chapman et al. 2016. In Fig. 21 we plot the power spectra and the extracted VAO wiggles of the foreground removed and wedge discarded mock data cubes.Surprisingly and interestingly, we find that if we perform foreground removal first and then discard the modes suspected to be contaminated in the wedge, the power loss discussed in the last paragraph is somewhat corrected.At  ≳ 0.01 Mpc −1 the recovered 21 cm power spectra and the extracted VAO wiggles are quite close to the original ones.
Probably because the power loss in foreground removal is mainly in the modes with small  ∥ .When we discard the modes in the wedge, we also discard many modes with small  ∥ .As a result, they will not bias the recovered power spectra.
For this reason, we believe that in our previous predictions, if we evaluate the influence of the foreground by such a method, the results do not change much.However, we note that this is because we assume the cosmological signal is isotropic.The real 21 cm power spectrum could be anisotropic, discarding the wedge modes may bias the isotropically averaged power spectrum (Jensen et al. 2016).

Figure 1 .
Figure 1.The streaming velocity field (top left), density field (top right), and collapse fraction field in the CDM model (bottom left), and in the axion model with  a = 10 −20 eV (bottom right), at  = 20.

Figure 2 .𝑐
Figure 2. The global spectrum of the 21 cm signal in different dark matter models.

Figure 3 .
Figure 3. Top: The 21 cm power spectra for different dark matter models.Errorbars refer to the sample variance of the simulated signal.The hatched regions around the CDM model are the SKA1-low uncertainties for Ω survey = 20 deg 2 , while filled regions are SKA2low uncertainties for Ω survey = 200 deg 2 .The observation time  obs = 2000 hour.Bottom: The pure wiggles signal extracted from the 21 cm power spectrum by removing a polynomial component.

Figure 4 .
Figure 4.The VAO wiggles extracted from 21 cm power spectra at different redshifts for CDM models with various  * .We mark the 21 cm frequencies on the right -axis.In each panel, the solid white line marks the redshift when x ≈ 1.The empty regions have  k ≥ 10 3 K, where our models are not valid.The hatched regions are beyond the frequency range of the SKA-low band.The vertical dashed line marks  = 0.05 Mpc −1 , the location of the strongest VAO peak.

Figure 5 .
Figure 5. Same to Fig. 4, however here we vary the   .

Figure 8 .
Figure 8.The extracted VAO wiggles in the axion model with various masses.For the purpose of doing comparison between different panels, we set the same colorbar ranges for all panels.

Figure 9 .
Figure 9.The 21 cm VAO wiggles at different redshifts for CDM models with different  esc,LW  LW .Here we use Eq.(25) as the critical mass for Pop III formation in minihalos.Since we need to compare the VAO wiggles of different models, we set the same colorbar ranges for all panels.The three panels correspond to  esc,LW  LW = 0.01 BL05 LW , 0.1 BL05 LW and  BL05 LW respectively.and freeze other astrophysical and cosmological parameters (suppose they can be derived from other observations or theories), we show simple demonstrations of the constraints on parameters  a and  * in Fig. 14, from the suspected SKA1low observations with survey area 20 deg 2 and SKA2-low observations with 200 deg 2 respectively.Δ = 20MHz, and  obs = 2000 hour.For SKA1-low, the main constraints are from 0.01 Mpc −1 ≲  ≲ 0.1 Mpc −1 .At smaller scale the instrumental noise is large, at larger scale the cosmic variance is large.In the bottom panel of Fig.3and in Fig.13, seemingly at smaller  the uncertainties are smaller, however, the relative uncertainties increase rapidly with decreasing , due to the cosmic variance.In such constraints, we do not involve the influence of the foreground, so this is an optimistic forecast.The foreground may bias the measured large-scale power spectrum, however, it can be somewhat corrected for, see the discussions in Sec.4.2.3.There is a degeneracy between the  a and  * in the SKA1-low constraints, since both lower  * and

Figure 10 .
Figure 10.Similar to Fig. 9, however here we use the Eq.(29) as the critical mass for Pop III formation in minihalos.The three panels correspond to  esc,LW  LW =  BL05 LW , 10 BL05 LW and 100 BL05 LW respectively.

Figure 11 .
Figure 11. to Fig. 9, however here we use the virial mass corresponding to circular velocity Eq. (32) as the critical mass for Pop III formation in minihalos.The three panels correspond to  esc,LW  LW = 0.1 BL05 LW ,  BL05 LW and 10 BL05 LW respectively.

Figure 18 .
Figure18.The amplitude of the first peak at  = 0.05 Mpc −1 for different axion massess with various fractions.The meanings of the filled regions are same to Fig.17.Since the panel will be too crowded if we plot all curves together, we use two panels.

Figure 21 .
Figure 21.Left: The 21 cm power spectra of the mock data cubes if we remove the foreground, and if we remove the foreground and discard the modes in the foreground contaminated wedge.Right: The extracted VAO wiggles from the left panels.Top row: bandwidth 20 MHz.Bottom row: bandwidth 40 MHz.