Probing Minihalos in the Late Dark Ages through Velocity Acoustic Oscillations in the 21-cm Power Spectrum

Introduction

The Dark Ages offer an ideal opportunity for probing the primordial fluctuations and dark matter (DM) properties, as the cosmological information has not yet been contaminated by radiation from astrophysical luminous objects. The redshifted 21 cm line provides the only known direct tomographic probe of the neutral hydrogen during the Dark Ages (see reviews e.g. [38, 89, 9, 34]. In addition to this, there are other potential observables, include the secondary CMB anisotropies from scattering by moving primordial molecular clouds ([93] (R. Maoli et al. 1994)), the CMB spectrum distortions arising from the transition lines of the first molecules ([55]), as well as the high-energy radiation of DM annihilation or decay ([14]). These observables provide indirect or redshift-integrated constraints on the matter distribution during the Dark Ages; whereas the 21 cm signal can directly map their three-dimensional distribution and redshift evolution. The 21 cm signal from the Dark Ages can probe small-scale matter density fluctuations ([60, 57, 39]), primor-

Corresponding author: Bin Yue

[email protected] dial non-Gaussianity ([19, 20, 96, 33, 8]), primordial magnetic fields ([109, 76]), DM particle properties ([97]), DM decay/annihilation or interactions ([112, 28, 117, 68, 82]), primordial black holes ([44, 118] (Z. Zhang et al. 2026)), and cosmological models ([67, 77]).

In the absence of radiation from luminous astrophysical objects, the 21 cm spin temperature is primarily determined by collisional coupling, whereas the gas kinetic temperature is governed by adiabatic cooling, Compton heating, and shock heating. Previous investigations have pointed out that, at $z \gtrsim10$ a significant fraction of the gas can be heated to above the CMB temperature by shock heating ([43, 36, 51] (S. R. Furlanetto et al. 2006)). M. Kuhlen et al. [54] investigated the spin temperature and brightness temperature of the neutral hydrogen in the late Dark Ages ($z \sim20$) by hydrodynamical simulations and found that collisions are only efficient in minihalos and filaments. Shock heating generates emission signal in these dense regions while the remaining regions still appear in absorption. The spatial contrast in temperature that is created between the over- and under-dense environments provides a favorable opportunity to probe the distribution of neutral hydrogen during this epoch. While Y. Xu et al. (2021) found that such effects can reduce the 21 cm global absorption by up to $\sim17\%$ in the absence of luminous astrophysical objects. H. Park et al. [79] and S. Sikder et al. [102] proposed that the clumping effect that arises from small-scale density fluctuations can change the 21 cm global spectrum and large-scale power spectrum, thus opening a new window into small-scale structures and DM properties through the global signal or the large-scale fluctuations.

Forming from primordial fluctuations and growing continuously during the Dark Ages, the first nonlinear DM structures act as tracers of the primordial fluctuations field. However, such structures are difficult to detect directly, unless the DM particles can annihilate into cosmic rays or photons ([4, 26]). Once the mass of a DM halo exceeds the Jeans mass ([3]), it can accrete gas from the intergalactic medium (IGM). The accreted gas in such a minihalo is finally heated to virial temperature which is generally much higher than the CMB temperature, then produces the 21 cm signal in emission ([48] (S. R. Furlanetto & S. P. Oh 2006)). In contrast, in the Dark Ages the IGM temperature is lower than the CMB, so its 21 cm signal is always in absorption.

I. T. Iliev et al. [48]; P. R. Shapiro et al. [98]; S. R. Furlanetto & S. P. Oh (2006) studied the 21 cm emission from minihalos in the Dark Ages. They found that a single minihalo can reach a high brightness temperature ($\sim$ several hundred mK); however, direct detection is severely limited by its compact angular scale and the blending of multiple minihalos within a single beam. Thus, for an interferometer, the received beam-flux is the cumulative contribution from many minihalos. The fluctuations level of the beam-flux would then be much lower, just several mK, but still detectable for realistic telescopes. Their works provide the theoretical basis for using the 21 cm signal to probe the first nonlinear structures. Nevertheless, in the standard $\Lambda$CDM cosmology, the 21 cm power spectrum of minihalos is generally much weaker than the IGM ([119, 96] (S. R. Furlanetto & S. P. Oh 2006)), because most of the gas resides in the IGM, just a tiny fraction collapses into minihalos ([17]). Therefore, it is difficult to extract the contribution of minihalos from the total 21 cm signal received.

The relative streaming motion between the DM and baryon is generated before recombination and is still supersonic in the Dark Ages. Y. Ali-Haïmoud et al. [1] investigated its effects on the 21 cm signal from the Dark Ages and found that at large scales, the signal is enhanced by several to tens of percent; at small scales, the signal is suppressed. Interestingly, at the extremely small scales, the signal is boosted instead. They pointed out that such effects allow us to probe small-scale matter power spectrum by measuring the large-scale signal.

This relic relative motion can prevent gas from being accreted into smaller minihalos. As a result, the gas content of minihalos is modulated by the relative velocities. As the spatial correlation of the relic DM-baryon relative streaming velocities has oscillatory features ([111, 110, 22, 1, 62, 29]), all physical processes relevant to the gas content in minihalos must also have similar oscillatory features. These are the velocity acoustic oscillations (VAO) features.

During the Cosmic Dawn, the formation of Pop III stars is modulated by the DM-baryon relative motion ([108, 45, 75, 30, 113, 56, 92]). The 21 cm signal of the IGM influenced by the Ly$\alpha$ and X-ray radiation from Pop III stars therefore exhibits oscillations. The earliest studies of this oscillatory imprint in the 21 cm signal were presented by N. Dalal et al. [22] and E. Visbal et al. [113], subsequently followed by M. McQuinn & R. M. O’Leary [62]; A. Fialkov et al. [31]; J. Mirocha et al. [65]; J. B. Muñoz [70]; S. C. Hotinli et al. [46]; D. Sarkar & E. D. Kovetz [91]; H. Shimabukuro et al. [101]; X. Zhang et al. [121]; Z. Zhang et al. [122]; X. Du et al. (2026) and so on, in a variety of astrophysical and cosmological settings. These studies showed that the large-scale acoustic structures of the DM-baryon relative velocity field can generate prominent oscillations on the 21 cm power spectrum in the Cosmic Dawn.

In the Dark Ages, both minihalos and the IGM have the 21 cm signal. For minihalos the signal is in emission, while for the IGM the signal is in absorption before it is heated to above the CMB. However, it is hard to separate minihalos and the IGM in the received total signal, because their signals have similar large-scale power spectrum shape and the IGM always dominates over the total signal. The DM-baryon relative motion provides us a chance to detect the contributions from minihalos in the received total signal. Since the gas content of minihalos is reduced by the relative motion (S. Naoz et al. [75]), their 21 cm signal must also have the oscillatory features. Such features remain largely unexplored in literature. If the 21 cm signal of the IGM is not or at most just very slightly influenced by this relative motion, then the VAO features provide us with unique information about minihalos through the 21 cm signal.

This is particularly important for the late Dark Ages ($z \sim30 - 20$, sometimes may extend to $z \sim15$ if the IGM has not yet been influenced by radiation from luminous objects) when minihalos are still small enough and more sensitive to the relative velocities, while the 21 cm signal of the IGM becomes weak because the collision-coupling becomes less efficient (R. Barkana & A. Loeb [3]; S. R. Furlanetto et al. 2006; J. R. Pritchard & A. Loeb [89]). Therefore, the late Dark Ages is a golden epoch for probing minihalos via the VAO features.

Regarding the IGM, in the presence of Pop III stars whose formation is modulated by the relative motion, its 21 cm signal would be influenced by the relative motion as well, through the Ly$\alpha$ and X-ray radiation from these stars (e.g. N. Dalal et al. [22]; E. Visbal et al. [113]). However, in the Dark Ages when Pop III stars have not yet formed, it remains unclear whether the 21 cm signal of the IGM also has oscillatory features. It relies on whether the shock heating—the dominant mechanism capable of changing the gas kinetic temperature—works efficiently in the IGM, because the shock heating

is modulated by relative motion (R. M. O’Leary & M. McQuinn [78]). Many papers (N. Y. Gnedin & P. A. Shaver [43]; M. McQuinn & R. M. O’Leary [62]; Y. Xu et al. 2021; J.-Y. Jia et al. [51]; H. Park et al. [79]) found that the shock heating can change the gas temperature significantly, but they did not analyze the IGM and halos separately. S. R. Furlanetto & A. Loeb [36] pointed out that even at $z \sim20$, a significant fraction of the IGM can be heated to $\gtrsim10^4$ K by shocks. If this is true, then the 21 cm signal of the IGM would also exhibit oscillatory features. However their work was based on analytical model, therefore probably not accurate for this problem. P. R. Shapiro et al. [98] also analyzed the 21 cm signal of the IGM and minihalos separately in their hydrodynamical simulations. Their results show that at $z \sim20$ the 21 cm signal of the IGM is basically in absorption, and is unlikely to be heated significantly. Neither of these two works incorporated the DM-baryon relative motion. Moreover, R. M. O’Leary & M. McQuinn [78] simulated the relative motion-induced shock heating and found that it is confined to the immediate vicinities of halos, whereas the IGM in a large volume is not significantly heated. Based on this, M. McQuinn & R. M. O’Leary [62] found that the temperature fluctuations raised by the relative motion-induced shock heating should be $\lesssim10\%$. Although not explicitly illustrated in the paper, it can be inferred from their results that in the Dark Ages, the IGM is unlikely to exhibit prominent VAO features in the 21 cm power spectrum. In addition to calculating the VAO signals from minihalos in the late Dark Ages, in this paper we will also evaluate the 21 cm signal from the diffuse IGM in the presence of DM-baryon relative motion, and aim to quantify whether the IGM will contaminate the VAO features of minihalos, via hybrid of hydrodynamic simulations and semi-numerical simulations.

The redshifted 21 cm signal from the Dark Ages falls in a very low-frequency band, where terrestrial ionospheric absorption/distortion and radio-frequency interference are severe. The lunar far-side is therefore one of the best locations for such observations. A number of lunar-based or lunar-orbit projects/concepts have been proposed, including DAPPER (J. O. Burns et al. [10]), LuSEE-Night (S. D. Bale et al. [2]), DSL (also known as Hongmeng) (X. Chen et al. [15]), FARSIDE (J. O. Burns et al. [11]), LARAF (X. Chen et al. [16]), FarView (R. S. Polidan et al. [88]; J. O. Burns et al. [12]), CoDEX (L. V. E. Koopmans et al. [53]), and DEX (C. D. Brinkerink et al. [6]). H. T. J. Bevins [5] gave precise predictions of the magnitude and frequency range of the sky-averaged 21 cm signal based on existing observational constraints on cosmology. J. C. Pober & W. Smith [84] quantified how foreground wedge avoidance would degrade the sensitivity of the Dark Ages 21 cm power spectrum observations. Y. Wen et al. [114] investigated the feasibility of designing a lunar far-side interferometer array for probing the small-scale 21 cm power spectrum and constraining cosmological parameters. C. Zhang et al. (2026) developed a wedge-aware Fisher framework to estimate the influence of foreground wedge-avoidance on constraining the primordial non-Gaussianity.

In this paper, we investigate the VAO features on the 21 cm power spectrum of minihalos in the late Dark Ages, and estimate their detectability for future lunar-based low-frequency arrays. The paper is organized as follows. In § 2, we introduce the 21 cm brightness temperature of individual minihalos. In § 3, we calculate the VAO wiggles on the 21 cm power spectrum of minihalos in the late Dark Ages. In § 4, we estimate the detectability of this signal using the lunar-based low-frequency array. In § 5, we investigate whether the 21 cm signal of the IGM will also be influenced by the DM-baryon relative streaming velocities and whether it will contaminate the signal from minihalos. In § 6 we discuss uncertainties relevant to the lightcone effect and the foreground. Finally in § 7 we give our summary. Throughout this paper, we adopt the Planck 2018 cosmology parameters for ΛCDM model, $\Omega_{\mathrm{m}} = 0.3153$, $\Omega_{\Lambda} = 0.6847$, $\Omega_{\mathrm{b}} = 0.0493$, $h = 0.6736$, $\sigma_{8} = 0.8111$, and $n_{\mathrm{s}} = 0.9649$ ([83]).

The 21 cm signal of an individual minihalo

In the late Dark Ages, the gas temperature of a minihalo is generally much larger than the CMB temperature, so at the 21 cm frequency it acts as a bright spot. I. T. Iliev et al. [48]; ); P. R. Shapiro et al. [98] has developed a model for calculating the brightness temperature of such a point source by assuming that its gas profile follows the truncated isothermal sphere (TIS) model ([99, 47]). Given the gas density profile $\rho_{\mathrm{gas}}(r)$ and temperature profile $T_{\mathrm{K}}(r)$, the spin temperature profile is ([35])

where $T_{\mathrm{CMB}}$ is the CMB temperature at redshift $z$; $x_{\mathrm{c}}$ is the collisional coefficient which depends on gas density and temperature ([124, 35]); $x_{\alpha}$ is the coefficient of hydrogen-Lyα scattering (the Wouthuysen-Field effect, [115, 32]).

Throughout this paper, we limit our investigation to the Dark Ages where there is no radiation field other than the CMB, so we ignore the $x_{\alpha}$ term in Eq. (1). Moreover, since gas in the virialized halo is dense (the overdensity $\delta\gtrsim$ several tens) and optically thick to Lyα photons, even in the presence of Lyα radiation, the coupling is still dominated by collisional processes.

In the TIS model, $T_{\mathrm{K}}$ is independent of $r$ and is set to be the virial temperature. Following P. R. Shapiro et al. [99]; ); I. T. Iliev & P. R. Shapiro [47], the virial temperature of a halo with mass $M_{\mathrm{h}}$ is

where $\mu$ is the mean molecular weight, which we set to $\mu= 1.22$ for the neutral gas in minihalos, and $h$ is the dimensionless Hubble constant, and

in which $\Omega_{\mathrm{m}}$ is the matter density parameter at present. The quantity

is the nonlinear virial overdensity with respect to the critical density at the collapse redshift ([7]) and

in which $\Omega_{\Lambda}$ is the dark energy density parameter today.

In the TIS model, the gas density profile of a halo is

where $\Omega_{\mathrm{b}}$ is the baryon density parameter, and

is the TIS normalization density. Following the standard TIS model, we identify the halo radius with the truncation radius,

and define the TIS core radius by

with $\zeta_t \simeq30.04$ (P. R. Shapiro et al. [99]; I. T. Iliev & P. R. Shapiro [47]). The spin-temperature profile $T_{\mathrm{S}}(r)$ is obtained by substituting the above $T_{\mathrm{vir}}$ and $\rho_{\mathrm{gas}}$ into Eq. (1).

The brightness temperature at the rest-frame 21 cm transition frequency $\nu_0 = 1420.4057$ MHz for a line of sight with projected distance $r_{\mathrm{p}}$ from the minihalo center is

where the geometry relation between $r'$ and $s'$ is given in Figure 1. The optical depth $\tau'_{\nu_0}$ is the integration of the absorption coefficient $\kappa_{\nu_0}$ from $s'$ to $2\sqrt{R_{\mathrm{h}}^2-r_{\mathrm{p}}^2}$. $\tau_{\nu_0}$ in the first term of the right hand side of Eq. (10) is the optical depth for $s'=0$. The absorption coefficient $\kappa_{\nu_0}(r)$ is

where $A_{10} = 2.85 \times10^{-15}\ \mathrm{s}^{-1}$ is the spontaneous decay rate of the 21 cm transition, $T_\star= h_{\mathrm{P}}\nu_0/k_{\mathrm{B}} = 0.068$ K, and $\phi(\nu_0)$ is the center of the normalized line profile. Following the thermal-Doppler approximation used by I. T. Iliev et al. [48], we take

and

for which $\sigma_{\mathrm{V}}$ is the one-dimensional velocity dispersion of a TIS halo, $c$ is the light speed. It is calculated from the TIS solution as

where $G$ is the gravitational constant, $H(z)$ is the Hubble parameter, and

Here $R_{\mathrm{TIS}} = 530$ and $\alpha_{\mathrm{p}} = 3.68$ are the dimensionless TIS parameters adopted in our calculation.

Figure 1 illustrates the geometry for Eq. (10). The face-averaged brightness temperature of this halo is

The differential brightness temperature in contrast to the CMB is

In Figure 2 we plot the $\langle\delta T_{\mathrm{b},\nu_0}\rangle_{\mathrm{h}}$ as a function of $M_{\mathrm{h}}$ at $z = 25$, $20$ and $15$ respectively. The brightness temperature drops rapidly below $\sim10^4\ M_\odot$.

The VAO features in the 21 cm power spectrum of minihalos

X. Zhang et al. [121] modified the original 21cmFAST (A. Mesinger et al. [63]; S. Murray et al. [72]; Y. Qin et al. [90]; J. Park et al. [80]; J. B. Muñoz et al. [71]; J. E. Davies et al. [24]) to involve the this code to generate the cosmological density field and the DM-baryon relative streaming velocity field. We set box length $L_{\mathrm{box}} = 1200$ Mpc and divide it into $400^3$ cells. Each cell has a size $x_{\mathrm{cell}} = 3$ Mpc, which is consistent with the coherence scale of the DM-baryon relative streaming velocities (D. Tseliakhovich & C. Hirata [111]; A. Fialkov et al. [30]; J. B. Muñoz [70]). Figure 3 shows slices of the density and velocity fields at $z = 20$.

In the presence of DM-baryon relative streaming velocity, the gas fraction in smaller minihalos is reduced. S. Naoz et al. [75] has investigated this using numerical simulations and has provided a fitting formula for the gas fraction as a function of halo mass and relative velocity $v_{\mathrm{db}}$,

where the shape parameter for the profile $\alpha= 1$, as given by N. Y. Gnedin [42]. The pivotal quantity in this formalism is the characteristic mass, $M_{\mathrm{c}}(v_{\mathrm{db}}, z)$, which represents the mass scale at which the gas content is reduced to half of the cosmic mean fraction.

We fit the characteristic mass $M_{\mathrm{c}}$ reported by S. Naoz et al. [75] as

where

$v_{\mathrm{ref}}^{\mathrm{N13}}$ is actually the rms of the relative streaming velocities in S. Naoz et al. [75].

Parameters $a_0$, $a_1$ and $a_2$ are independent of $z$ and are approximated by,

The S. Naoz et al. [75] simulations used the WMAP 5-year cosmology parameters (E. Komatsu et al. [52]), while we use the Planck 2018 cosmology parameters. In principle there are inconsistencies. However, we believe this will just have minor impacts on our results. We check that in the Planck 2018 cosmology, the rms of the relative streaming velocities $\sigma_{v_{\mathrm{db}}} = 28.85\ {\rm km\ s}^{-1}$ at the recombination era, just $\sim4\%$ different from the WMAP 5-year cosmology which is 30 ${\rm km\ s}^{-1}$. In our calculations, we use the $M_{\mathrm{c}}$ in Eq. (19) as a sharp cutoff. Minihalos above this characteristic mass are not influenced by $v_{\mathrm{db}}$; while minihalos below this mass contain no gas. In such a treatment, the detailed dependence of the $f_{\mathrm{g}}$ shape on cosmology is irrelevant to our calculations. The oscillatory features in the 21 cm signal are generated by the spatial variations of $M_{\mathrm{c}}$ that arise from the spatial variations of $v_{\mathrm{db}}$. We test that, if we vary the $v_{\mathrm{db}}$ by $\sim4\%$ globally, the amplitude of the minihalos VAO wiggles would just change by $\sim3.6\%$ at $z = 20$.

For a cell with density $\delta_{\mathrm{cell}}$ and DM-baryon relative streaming velocity $v_{\mathrm{db}}^{\mathrm{cell}}$, we calculate the 21 cm signal from minihalos in the cell

where $M_0^{\mathrm{cell}} = \rho_{\mathrm{m}}x_{\mathrm{cell}}^3$ is the mean mass of one cell, $\rho_{\mathrm{m}}$ is the mean matter density of the Universe, and $A(M_{\mathrm{h}},z) = \pi R_{\mathrm{h}}^2$ is the geometric cross section of a minihalo.

For simplicity we set the minimum mass of minihalos sustaining gas in the presence of $v_{\mathrm{db}}$ as $M_{\min}=M_{\mathrm{c}}$. We set $M_{\max}$ as the virial mass corresponding to $T_{\mathrm{vir}}=10^{4}\ \mathrm{K}$. The conditional halo mass function $\mathrm{d}n/\mathrm{d}M_{\mathrm{h}}$ is given by the Extended Press–Schechter (EPS) scenario H. J. Mo & S. D. M. White (1996),

in which $\sigma(M_{\mathrm{h}})$ is the rms of the density field fluctuations, smoothed at mass scale $M_{\mathrm{h}}$, and

The overdensity

where $\delta_{\mathrm{c}}(z)$ and $\delta_{\mathrm{L}}$ are linear extrapolations of critical overdensity for halo collapse and $\delta_{\mathrm{cell}}$ to present-day H. J. Mo & S. D. M. White (1996) respectively.

Since the volume-averaged EPS abundance does not exactly reproduce the Sheth–Tormen (ST) mass function R. K. Sheth & G. Tormen [100] for the mean Universe, we renormalize the conditional EPS mass function in Eq. (23) by a factor defined as

where the numerator is the ST halo mass function, and the denominator is the average of Eq. (23) taken over all cells in the cosmological box. This prescription is quite popular (e.g. A. Mesinger et al. [63]), it preserves the cell-to-cell modulation predicted by EPS, while ensuring that the box-averaged halo abundance is normalized to the ST mass function.

The 21 cm brightness temperature of the IGM in each cell is J. R. Pritchard & A. Loeb [89]

where the spin temperature of the IGM, $T_{\mathrm{S}}^{\mathrm{IGM}}$, is calculated from Eq. (1), by replacing the halo density and temperature with those of the IGM in that cell. The IGM overdensity satisfies $1+\delta_{\mathrm{IGM}}^{\mathrm{cell}}=(1+\delta_{\mathrm{cell}})(1-f_{\mathrm{coll}}(\delta_{\mathrm{cell}},v_{\mathrm{db}}^{\mathrm{cell}}))$, where $f_{\mathrm{coll}}$ is the collapse fraction of minihalos above $M_{\mathrm{c}}$ in that cell. The IGM temperature in the absence of astrophysical radiation is simply $T_{\mathrm{K}}^{\mathrm{IGM}}=T_{0}(1+\delta_{\mathrm{IGM}}^{\mathrm{cell}})^{2/3}$, where $T_{0}$ provided by RECFAST S. Seager et al. [94, 95].

Figure 4 exhibits slices of the 21 cm fields at $z=20$ for minihalos and the IGM, respectively. For the three-dimensional boxes, at $z=20$, the mean signal is $0.23\ \mathrm{mK}$ for minihalos, and $-0.35\ \mathrm{mK}$ for the IGM. While at $z=15$, these two mean signals are $0.71\ \mathrm{mK}$ and $-0.07\ \mathrm{mK}$ respectively. The power spectra of these respective 21 cm signal and the total signal fields are shown in the top panel of Figure 5. We note that at $z=20$ the power spectrum of the total signal is even smaller than that of the individual components, i.e. minihalos or the IGM. This is because the total power spectrum includes a cross-correlation term between the minihalo emission and the IGM absorption fields. This term arises because minihalos form preferentially in overdense regions where the IGM density is also higher. Moreover, the 21 cm absorption of the IGM is also stronger in these overdense regions, where the collision coupling works more efficiently. Given that $\delta T_{\mathrm{b,MH}}^{\mathrm{cell}}>0$ (emission) and $\delta T_{\mathrm{b,IGM}}^{\mathrm{cell}}<0$ (absorption) in the late Dark Ages, this cross-correlation is expected to be negative, thereby partially suppressing the total power spectrum amplitude. This is actually good for detecting the VAO wiggles. The VAO wiggles are primarily imprinted by the modulation of the minihalo abundance by the large-scale relative velocity field, whereas the IGM response to these velocities is negligible. Consequently, the shape and phase of the VAO features remain dominated by the minihalo contribution, but the overall baseline amplitude is suppressed by the cross-correlation term. Therefore the VAO wiggles become more apparent. This ensures that the detectability forecasts presented in § 4, which focus on the significance of the wiggle features, remain robust. Same to J. B. Muñoz [70], we fit the 21 cm power spectrum by fifth-order polynomial, then remove the fitted smooth component to get the net VAO wiggles. They are plotted in the bottom panel of Figure 5.

The detectability

We then estimate the detectability of the net VAO wiggles in the 21 cm power spectrum. We adopt the same interferometric lunar-array configuration as the first line of the second item in Tab. 2 of W. Smith & J. C. Pober [104]: each sub-array is composed of 3 × 3 cross-shaped dipoles; then, 4 × 4 sub-arrays constitute a group; 2 × 2 groups, in turn, form one cluster; finally, the entire array consists of 12 × 12 clusters and has 82944 dipoles in total, spread over ∼ 3.5 km. The layout of the array is similar to Fig. 1 of W. Smith & J. C. Pober [104], except that the 4 × 4 dipoles in each sub-array (denoted by red squares in that figure) are replaced with 3 × 3 dipoles. This configuration has signal-to-noise ra tio (S/N) larger than the original FarView design [88, 12].

We quantify the detectability of the net VAO wiggles by summing the S/N over the relevant $k$ bins,

where $\Delta^2_{\mathrm{wiggles}}(k_i)$ is the net VAO wiggles amplitude in the $i$th $k$ bin. The uncertainties $\sigma_{\Delta^2}$ combine the 21cmSense forecast on observational uncertainties (include thermal noise and cosmic variance) with the uncertainties of the theoretical signal in quadrature,

The top panel of Figure 6 shows the forecast observational uncertainties on the net VAO wiggles at $z=15$ using this array configuration. The uncertainties are calculated using the package 21cmSense [73, 85, 86 A. R.]

[81]), for an integration time of 24 lunar nights¹ which are equal to 8064 hr. We adopt the optimistic foreground-subtraction model in 21cmSense, which removes only the modes inside the primary-beam foreground wedge [87]. The wiggles are obviously detectable, with $(S/N)_{\mathrm{tot}}\approx4.7$. At $z=20$, however, the signal is too weak to detect with this array configuration. Thus, we extend our array to 22$\times$22 clusters and show the predicted observational uncertainties for integration time of 60 lunar nights (20160 hr) in the bottom panel of Figure 6. With such an array configuration and integration time, the net VAO wiggles are detectable as well, with $(S/N)_{\mathrm{tot}}\approx7.4$.

The VAO wiggles in cosmology beyond the standard $\Lambda$CDM model

In the standard $\Lambda$CDM model, at $z=20$, the first peak (at $k=0.05\ {\rm Mpc}^{-1}$) of the net VAO wiggles from minihalos has amplitude $\sim3.1 \times10^{-4}\ {\rm mK}^{2}$. For an interferometer to measure this scale at $z=20$, a baseline of $\sim400$ m is required, which necessitates the construction of a huge array. Moreover, since minihalos are not luminous, to achieve a statistically significant detection, it requires the dense baseline distribution formed by orders of $10^{4}$–$10^{5}$ antennas, as we have just demonstrated.

There are, however, many beyond standard $\Lambda$CDM models in which the amplitude of small-scale matter power spectrum is much larger than in the standard $\Lambda$CDM model, but still consistent with the CMB observations. Examples include a running spectral index for primordial power spectrum (T. Minoda et al. [64]); a broad enhancement followed by a DM microphysical cutoff (D. Inman & K. Kohri [49]); a bump or tilt on dwarf galaxy scales (O. R. Skorikov et al. [103]); and a localized Gaussian form power excess with a possible cutoff at still smaller scales (E. O. Nadler et al. [74]). M. Cielo et al. [18] proposed a model for which, compared with the standard $\Lambda$CDM model, the primordial scalar power spectrum is enhanced on small scales by sudden phase transitions during inflation. There is a characteristic wavenumber $k_{\rm trans}$, around which the power spectrum first increases as $\propto k^{6}$; reaching a peak; then drops gradually. It forms a “bump” on the power spectrum as exhibited in Fig. 1 of H. Zhu et al. [123]. This model has the advantage that, only the halos in a certain mass range is boosted, see Fig. 2 of H. Zhu et al. [123]. By choosing appropriate $k_{\rm trans}$, the number of minihalos is boosted, but massive halos that have the potential to host galaxies are unchanged or just minorly changed. H. Zhu et al. [123] found that in this model, the constraints on reionization history are not broken because the boost on ionizing photons sinks (minihalos) balances the boost on ionizing photons sources. However, the morphologies of the ionization and the 21 cm fields change a lot, resulting in different bubble size distribution and 21 cm power spectrum.

In all these models, because the matter power spectrum is only enhanced in a certain $k$ range, compared with the standard $\Lambda$CDM model, minihalos could be orders of magnitude more numerous, while the number of massive halos capable of hosting luminous objects remains almost unaffected. The resulting enhancement of the 21 cm signal from minihalos provides an opportunity to constrain these non-standard small-scale power spectra through their VAO features. For the quantitative illustration below, we adopt the sudden-transition model of M. Cielo et al. [18].

In Figure 7, we plot the net VAO wiggles in the 21 cm power spectrum from minihalos at $z=20$ with a transition wavenumber of $k_{\rm trans}=100$, 200 and 1000 ${\rm Mpc}^{-1}$ respectively, along with predicted uncertainties for different array configuration and/or integration time. Other parameters are set the same as H. Zhu et al. [123]. The figure exhibits that in such a model, the VAO wiggles become much larger. For $k_{\rm trans}=100\ {\rm Mpc}^{-1}$, at $z=20$, the VAO signal improves by roughly an order of magnitude than in the $\Lambda$CDM model. When the transition wavenumber is pushed to 1000 ${\rm Mpc}^{-1}$, the signal is already on the order of several ${\rm mK}^{2}$. As a result, the wiggles can be detected with arrays consisting of fewer clusters of antennas. Under the $\Lambda$CDM model, $22 \times22$ clusters and 60 lunar nights (20160 hr) integration time are required for detection at $z=20$. In this model, $12 \times12$ clusters would suffice for the more conservative case of $k_{\rm trans}=100\ {\rm Mpc}^{-1}$. When the signal is greatly enhanced in models of a larger $k_{\rm trans}=1000\ {\rm Mpc}^{-1}$, only $4 \times4$ clusters with 6 lunar nights (2016 hr) of integration are needed—required observational resources are greatly reduced.

The contamination from the IGM

The IGM produces 21 cm signal whose absolute amplitude and power spectrum could be much larger than those of the minihalos. If the IGM temperature is also influenced by the DM-baryon relative motion, then its power spectrum would also exhibit oscillatory features and may contaminate the VAO features of minihalos. However, if its temperature is not or just slightly influenced by the relative motion, then its power spectrum will always be smooth and will not be confused with the VAO features of minihalos. In the above calculations we have ignored this possibility. In this section, we will investigate whether the IGM also has oscillatory features and whether it would contaminate the minihalso by using a hybrid of high-resolution small-box hydrodynamic simulations and large-scale semi-numerical simulations.

The small-box high-resolution hydrodynamical simulations

We modify the Gadget-4 code (V. Springel et al. [107]; V. Springel [105]; V. Springel et al. [106]) to include cooling and heating (Y. Xu et al. 2021). We perform hydrodynamical simulations with a 3 Mpc box size and $2 \times1024^3$ particles. The mass resolutions of DM and gas particles are approximately $8.4 \times10^2\ M_\odot$ and $1.6 \times10^2\ M_\odot$, respectively. Such a hydrodynamical simulation box mimics one cell in our cosmological box with $L_{\mathrm{box}} = 1200$ Mpc. We then add different values of bulk DM-baryon relative streaming velocity $v_{\mathrm{db}}^{\mathrm{cell}}$ to the hydrodynamical simulation box to study its impacts on the IGM. We use the Friends-of-Friends (FoF) algorithm (M. Davis et al. [25]) with linking length $b_{\mathrm{lin}} = 0.2$ to find all DM particles in halos. Then, we remove all gas particles within the virial radius of each halo before calculating the net contribution to the mean 21 cm brightness temperature from the IGM.

Figure 8 shows the 21 cm signal slices in our hydrodynamical simulation box, for bulk velocities of zero and $\sim3\sigma_{v_{\mathrm{db}}}$ respectively. Results for minihalos and the IGM, and their sum, are displayed separately. The figure clearly indicates that, the bulk velocity changes the 21 cm signal of the minihalos significantly, while just very slightly influences that of the IGM. We note that even in the IGM only panel, there are still some bright emission spots. We check that they are hot and dense gas particles in the filaments outside virial radius of any minihalos. Although we present the 21 cm signal from minihalos in Figure 8, we also check that for these minihalos, the gas temperature drops to roughly half of the virial temperature at the virial radius; and the inner gas density is several factors lower than the TIS model. Therefore the 21 cm signal of individual minihalos in our simulations is smaller than the TIS model. It is not clear whether this is because our simulations suffer the resolution problem (many halos just have several tens to several hundreds gas particles), or because the TIS model is too optimistic for the 21 cm emission. Nevertheless, we still adopt the TIS model to calculate the 21 cm emission of minihalos. The minihalos panel in Figure 8 is just used to intuitively show the DM-baryon relative motion effects.

The difference between the mean 21 cm signal of the IGM in a 3 Mpc hydrodynamical simulation box with and without the bulk relative streaming velocity,

as a function of $v_{\mathrm{db}}^{\mathrm{cell}}$, is plotted in Figure 9 at $z = 20$. The magnitude of $\Delta[\delta T_{\mathrm{b,IGM}}^{\mathrm{cell}}]$ appears to increase with $v_{\mathrm{db}}^{\mathrm{cell}}$, but by less than $\sim0.01$ mK when the bulk velocity $\lesssim3\sigma_{v_{\mathrm{db}}}$. To further investigate the positions of the gas particles contributing to $\Delta[\delta T_{\mathrm{b,IGM}}^{\mathrm{cell}}]$, we also remove the gas particles within $2R_{\mathrm{vir}}$ and $3R_{\mathrm{vir}}$ respectively. The resulting $\Delta[\delta T_{\mathrm{b,IGM}}^{\mathrm{cell}}]$ is shown in Figure 9 as well. We see that in this case, the magnitude becomes $\lesssim0.002$ mK. It means that, most of the shock-heated IGM gas particles are actually located in the vicinity of halos. This is consistent with the discoveries in R. M. O’Leary & M. McQuinn [78]. By inspecting the simulation outputs, we also find that some shock-heated gas particles are located in filaments even beyond $3R_{\mathrm{vir}}$ from the center of FoF halos. This indicates that the separation between the immediate halo environment and the diffuse IGM is not perfectly sharp, but it does not affect our conclusion that, for the 21 cm signal, the IGM response to the bulk relative velocity is negligible compared with the minihalos signal.

The hybrid method and large-scale VAO features of the IGM

Although the small-box high-resolution hydrodynamical simulations reveal that the 21 cm signal of the IGM is quite less sensitive to the relative motion, they cannot directly provide the estimate of the VAO wiggles signal for the IGM, because the wiggles are large-scale features. On the other hand, however, it is quite impractical to run hydrodynamical simulations with box size as large as our cosmological box, while at the same time can resolve the smallest minihalos. A popular method for solving this problem in literature is the hybrid method (e.g. E. Visbal et al. [113]; A. Fialkov et al. [31]; M. McQuinn & R. M. O’Leary [62]; A. Fialkov [29]; H. Park et al. [79]; S. Sikder et al. [102]). Using many small-box hydrodynamical simulations with same box size equal to the cell size of the large-scale cosmological box, however with different bulk relative velocities and background overdensities, to construct a table of 21 cm signals as a function of velocity and overdensity. Then, for each cell in the large-scale box, its 21 cm signal is interpolated based on this table.

In our work, using the hybrid method we should add a velocity- and overdensity-dependent correction term

to the 21 cm field generated from Eq. (27), to account for the relative motion effects on the IGM’s 21 cm signal.

However, all our small-box hydrodynamical simulations have zero background overdensity, in principle we do not know how their 21 cm signal depends on the background overdensity. We analyze the cross-correlation between our large-scale density field and relative velocity field, and find that, the cross-correlation coefficient is at $\mathcal{O}(10^{-4})$ level. So these two fields are almost independent of each other (C. Cain et al. [13]). We therefore reasonably hypothesize that, the behavior of $\Delta[\delta T_{\mathrm{b,IGM}}^{\mathrm{cell}}(v_{\mathrm{db}}^{\mathrm{cell}},\delta_{\mathrm{cell}})]$ as a function of background overdensity $\delta_{\mathrm{cell}}$, is independent of relative motion. That is

We try three hypotheses: (1) $\mathcal{F}=1$; (2) $\mathcal{F}=1+\delta_{\mathrm{IGM}}^{\mathrm{cell}}$; and (3) $\mathcal{F}=\delta T_{\mathrm{b,IGM}}(\delta_{\mathrm{IGM}}^{\mathrm{cell}})/\delta T_{\mathrm{b,IGM}}(0)$, for which $\delta T_{\mathrm{b,IGM}}$ is the 21 cm brightness temperature of the IGM in the absence of relative motion, expressed in Eq. (27).

For all of these three hypotheses, we find that the oscillatory features of the IGM arising from the relative motion are almost negligible compared with minihalos. In Figure 10 we plot the net wiggles in the 21 cm power spectrum of the IGM at $z=20$, in the absence (top left panel) and presence (top right panel for hypothesis 1, bottom left panel for hypothesis 2 and bottom right panel for hypothesis 3) of relative motion correction respectively. For comparison we also plot the VAO wiggles from the minihalos (same to the bottom panel in Figure 6). We find that, even in the absence of DM-baryon relative motion correction, the 21 cm power spectrum of the IGM has its own wiggles. The locations of these wiggles are different from the VAO wiggles, they are actually generated by the baryon acoustic oscillations (BAO) in the density field [121, 21]. Such BAO wiggles may confuse the VAO signal, fortunately, their amplitude is obviously smaller than the VAO wiggles of minihalos. Moreover, even after adding the bulk motion-induced temperature correction to the IGM, the wiggles in its 21 cm power spectrum change very little. It means that the effects of the DM-baryon relative motion on the IGM is negligible compared with the IGM’s BAO wiggles and the VAO wiggles of minihalos. To further confirm this, we also plot the power spectra of the correction term for the three hypotheses respectively in Figure 10. Indeed, they are negligible compared with the VAO wiggles of minihalos. We conclude that the contamination from the IGM can be safely ignored in detecting the VAO features of minihalos in the late Dark Ages.

We clarify that, when comparing the VAO features with the BAO features, the situation is different for the Dark Ages and the Cosmic Dawn, and different for minihalos and the IGM.

In the Dark Ages, astrophysical luminous objects have not yet formed, so the amplitudes of the VAO and BAO wiggles do not depend on astrophysical parameters. The 21 cm signal from minihalos has both the VAO wiggles and BAO wiggles. This is because the BAO component is inherited from acoustic features in the density field, while the abundance of minihalos is modulated by their environment densities. In contrast, the VAO component considered here arises because the coherent DM-baryon relative velocity field modulates gas content of low-mass minihalos. Its amplitude is therefore sensitive to the velocity-dependent mechanism that prevent gas from being accreted into minihalos. That is to say, the amplitude of the VAO wiggles is sensitive to the lower limit of minihalos; while the amplitude of the BAO wiggles depends on the cumulative gas content in all minihalos. So the shape of minihalo mass function is crucial in determining their relative amplitude, which of course depends on cosmology parameters. Regarding the IGM in the Dark Ages, according to our hypotheses, its VAO wiggles are much weaker than the BAO wiggles. The amplitude of BAO wiggles relies on the primordial matter power spectrum and the growth of density perturbations, which is further dependent on cosmology parameters.

In the Cosmic Dawn, the IGM has both the VAO and BAO wiggles, and their relative amplitude evolves along with the formation of first stars [121, 21]. At the beginning of the Cosmic Dawn, a small amount of the first stars just form, the IGM only exhibits BAO wiggles and its VAO wiggles are negligible, like in the Dark Ages. As the first stars form in large numbers and the IGM is significantly influenced by their Ly$\alpha$ and X-ray radiation, the VAO wiggles become much higher than the BAO wiggles. Later on, when the Ly$\alpha$ scattering and X-ray heating become saturated, the VAO wiggles disappear but the BAO wiggles always exist because it is directly from the coherent structures in the density field. The formation of first stars is regulated by both astrophysical parameters and cosmological parameters, so in the Cosmic Dawn, both the VAO and BAO wiggles and their relative amplitude depend on these parameters.

In this work, we limit our investigations to the Dark Ages. We clarify that our results about the amplitudes of the VAO and BAO wiggles, and their comparison, apply only to the cosmology and minihalo model adopted in this paper. Although the characteristic acoustic shape of the VAO template is comparatively robust [70], the relative BAO and VAO amplitude varies with model parameters.

Discussion

The VAO features from minihalos evolve rapidly and the lightcone effect must be significant. For example, at $z=30$, the amplitude of the first peak at $k=0.05\ \mathrm{Mpc}^{-1}$ is $9.7 \times10^{-7}\ \mathrm{mK}^{2}$; at $z=20$ it is $3.1 \times10^{-4}\ \mathrm{mK}^{2}$; and at $z=15$ it is $1.4 \times10^{-3}\ \mathrm{mK}^{2}$. We have ignored the lightcone effect and may overestimate the (S/N)tot. In such a case, the multi-frequency angular power spectrum is more appropriate than the spherically-averaged 3D power spectrum (X. Du et al. 2026).

Another challenge is the foreground radiation from our Milky Way and extragalactic sources. It can be $\sim5$ orders of magnitude higher than the 21 cm signal from the Dark Ages. Although it is spectrally smooth and many foreground removal algorithms take advantages of this feature, the chromatic response of the interferometer mixes angular foreground structures into frequency structures, contaminating the signal in a wedge-like region in the $(k_{\perp}, k_{\parallel})$ space, where $k_{\perp}/k_{\parallel}$ is the wavenumber perpendicular/parallel to the line of sight (e.g. A. Datta et al. [23]; M. F. Morales et al. [69]; A. R. Parsons et al. [81]; J. C. Pober et al. [86]; A. Liu et al. [58]; H. Jensen et al. [50]; S. Gagnon-Hartman et al. [40]). The wedge border is set either by the horizon (horizon-limit) or by the primary beam (primary beam-limit). Modes outside this region are relatively cleaner. The most straightforward method is to only use these clean modes, and this is the foreground-avoidance. However, the impact of the wedge becomes more severe at higher redshift because the slope of its upper boundary increases with redshift, reducing the number of accessible modes outside the wedge. Therefore, the foreground contamination presents a greater challenge for the Dark Ages than the epoch of reionization. J. C. Pober & W. Smith [84] quantified this effect for the Dark Ages experiments and found that, using the horizon-limit wedge, foreground-avoidance typically reduces the power spectrum detection significance by about an order of magnitude. For their fiducial array, at $z = 30$, the nominal no-foreground detection significance for the 21 cm power spectrum decreases from $10\sigma$ to $\sim1.1\sigma$ after the modes in the foreground wedge are discarded. Similarly, P. Bull et al. [8] found that, excluding the modes in the worst case foreground wedge degrades the constraints on the non-Gaussianity in the bispectrum by a factor of $\gtrsim5 - 100$.

S. Ghosh et al. [41] performed the first end-to-end simulation pipeline for the Dark-ages EXplorer. Their mock data includes the 21 cm signal from the Dark Ages, the foreground, thermal noise, and the lunar topography. They found that: using the polynomial fitting and spectral tapering, the 21 cm signal from $z \sim42$ can be recovered, although the antenna-position and lunartopography errors can bias the recovered power spectrum by $\sim10 - 30\%$; but the signal from $z \sim148$ is still difficult to access. Their findings boost our confidence in detecting the 21 cm signal from the late Dark Ages ($z \sim20$).

Since the VAO wiggles are large-scale features, foreground-removal is necessary ([121]). We note that in S. Ghosh et al. [41], even for the $z \sim42$ signal, the modes with $k_{\parallel} \lesssim0.2\ \mathrm{Mpc}^{-1}$ are still contaminated by foreground. In this case, the multi-frequency angular power spectrum that only counts the correlation from modes perpendicular to the line of sight could be useful (X. Du et al. 2026), and probably more advanced foreground-removal algorithms are required.

In this paper, we adopt the foreground-avoidance strategy. However, our primary beam is the synthetic beam of the sub-array, which is much smaller than a dipole. Therefore our foreground-wedge is actually much smaller than the horizon-limit and the dipole beam-limit. Our results are actually close to the perfect foreground-removal case, so the forecasts presented in this paper are considered as the optimistic case.

Summary

In this paper, we investigated the VAO features on the 21 cm power spectrum of minihalos in the late Dark Ages ($z \sim20$) and their detectability by low-frequency array at the farside of lunar surface. During this epoch, if the IGM has never been influenced by the radiation from Pop III stars, its 21 cm signal is still in absorption, but much weaker due to the less efficient collisional coupling compared with higher redshifts, while the signal from minihalos is in emission. At $z \gtrsim30$, the fluctuations amplitude of the IGM’s 21 cm signal is much larger than that of minihalos, but they become comparable at $z \sim30 - 20$. Generally, it is not easy to isolate the contribution from minihalos in the detected total 21 cm signal. However, the relic supersonic DM-baryon relative streaming velocities suppress the gas accretion onto minihalos and modulate their 21 cm signal, finally leaving the VAO features in the power spectrum. On the other hand, the 21 cm signal from the IGM is just slightly influenced by such relative velocities and exhibit negligible VAO features compared with that from minihalos. This provides us with a chance to probe minihalos in the late Dark Ages. We calculated the VAO wiggles and estimated the observational uncertainties, we found that:

1. In the standard $\Lambda$CDM cosmology, the first peak of the net VAO wiggles at $k = 0.05\ \mathrm{Mpc}^{-1}$ has amplitude $\sim1.4 \times10^{-3}\ \mathrm{mK}^{2}$ at $z \sim15$, and $\sim3.1 \times10^{-4}\ \mathrm{mK}^{2}$ at $z \sim20$, respectively.

2. If the IGM has not yet been influenced by radiation from Pop III stars, an array with 82944 dipoles, like the one proposed in W. Smith & J. C. Pober [104], can detect the VAO features from minihalos at $z \sim15$ with 24 lunar nights (8064 hr). If the array can be extended by a factor of $\sim3$, then it can detect the VAO features from $z \sim20$ with 60 lunar nights (20160 hr).

3. In some cosmology models, the small-scale matter power spectrum is enhanced compared with the standard $\Lambda$CDM model. In such models, minihalos can produce much stronger VAO wiggles that are detectable for a smaller array, providing us with a new chance for constraining these models.

Our investigation provides a new strategy for probing the small-scale nonlinear structures in the late Dark Ages by the lunar-based low-frequency array on the farside of the Moon in the future.

Acknowledgments

We thank the anonymous referee very much for the constructive suggestions that help to improve this manuscript. This work is supported by the NSFC International (Regional) Cooperation and Exchange Project No. 12361141814, the National SKA Program of China Nos. 2020SKA0110402 & 2020SKA0110401, the China’s Space Origin Exploration Program Nos. GJ1110405 & GJ11010401, the National Key R&D Program of China No. 2022YFF0504300, and the Project Supported by the Specialized Research Fund for State Key Laboratory of Radio Astronomy and Technology. The work was carried out at National Supercomputer Center in Tianjin, and the calculations were performed on Tianhe new generation supercomputer.