Introduction. Discovering the fundamental nature of dark matter (DM) is one of the leading problems in the current era of physics. While the model of cold dark matter (CDM) is successful at explaining the large-scale structure of the Universe [1], it is largely unconstrained on smaller scales (\(\sim\) sub-Mpc). Additionally, the leading particle candidate for CDM, a weakly interacting massive particle, has only had consistent null results from collider and direct detection searches [2]–[4]. These results have motivated considerations for non-minimal models of DM which deviate from the small-scale structure predicted in \(\Lambda\)CDM. A distinctive feature of many such models is the presence of dark acoustic oscillations (DAOs) [5]–[7], which result in an overall suppression and oscillation feature in the Universe’s matter power spectrum (see Fig. 1)

DAOs can arise due to a coupling between DM and dark radiation in the dark sector, completely analogous to the process that causes baryon acoustic oscillations (BAOs) in \(\Lambda\)CDM. DAOs are additionally motivated as they can ameliorate discrepancies between observations and \(\Lambda\)CDM that have emerged in the last few years, the most significant being the Hubble, \(S_{8}\) and CMB-BAO tensions [8]–[13].

One of the most theoretically motivated dark sectors giving rise to DAOs is atomic dark matter (aDM) [5]. An aDM subcomponent can cause deviations at larger [6], [7], [9], [9], [14], [15], smaller [16]–[19] and stellar [20]–[29] scales, and is predicted by many models solving the Little Hierarchy Problem [30]–[36]. However, since DAOs also arise in a wider class of dark sector models [37]–[42], we will focus on DAOs as a generic dark sector signature in this letter.

To this end, it is useful to implement an effective parameterization of DAOs in the linear matter power spectrum. In the effective theory of structure (ETHOS) formalism [43], [44], the DAO transfer function (ratio of the linear matter power spectrum to the \(\Lambda\)CDM limit) is specified by the wavenumber and height of the first DAO peak, assuming a fully interacting dark sector. Further extending the range of theories covered by this parameterization, Ref. [14] allowed the DM that interacts with the dark radiation to be a subset of the total DM abundance. We adopt their model in this work.

Since DAOs imprint themselves on the matter power spectrum at small scales, the Lyman-\(\alpha\) forest, as a tracer of the small-scale (sub-Mpc), high-redshift (\(z \sim 5\)) linear matter power spectrum [45], [46], makes an excellent probe. In particular, high-resolution Lyman-\(\alpha\) forest spectra (which we use here) probe much smaller scales [47] than current cosmic microwave background (CMB) experiments like Planck [1], the Atacama Cosmology Telescope [48]–[50] and the South Pole Telescope [51]. Measurements of the high-redshift galaxy UV luminosity function (UVLF) have placed the strongest bounds on the DAO scale to date [14], albeit given different priors than we use here. The Lyman-\(\alpha\) forest is a powerful probe of even smaller scales than the UV luminosity function, and so we use the Lyman-\(\alpha\) forest in this work to extend our sensitivity to DAO models.

The Lyman-\(\alpha\) forest is a spectral absorption signal arising from neutral hydrogen in the intergalactic medium (IGM), away from the highest-density regions of the cosmic web. This environment makes the Lyman-\(\alpha\) forest a powerful probe of DM behavior, as, at any given scale, it is less affected by non-linearities and galactic feedback processes than galaxy tracers. Previous studies have used the Lyman-\(\alpha\) forest to constrain alternative DM models such as warm dark matter (WDM) [52]–[55], axion DM [56]–[58], interacting DM [59], decaying DM [60], [61], compact object DM [62] and mixed DM models [63], [64]. Further, the Lyman-\(\alpha\) forest appears to be uniquely sensitive to DAOs, as previously pointed out in Ref. [65]. Existing work on the Lyman-\(\alpha\) forest signature of dark matter models with DAOs has focused on a few benchmark models [65], the weak DAO limit [66], [67], or ignored the oscillations, using a WDM-like model [68]. Our work is the first to obtain constraints while fully accounting for the DAO feature as modeled through a suite of hydrodynamical simulations.

However, cosmological hydrodynamical simulations of the IGM are computationally expensive. A direct sampling of the parameter posterior distribution by such simulations is impractical. Accurate comparison of data to theory is enabled by machine learning (ML) emulator models and active learning methods [57], [59], [69]–[71]. The emulator model is trained on a limited set of simulations and is then used to interpolate the Lyman-\(\alpha\) forest flux power spectrum across the model parameter space. Active learning allows informed selection of new training simulation points after each iteration based on the knowledge gained from previous training sets. The emulator solution has been shown to accurately constrain cosmological, DM and IGM parameters using the Lyman-\(\alpha\) forest [55], [56], [59].

In this work, we combine the probabilistic modeling used in previous works with deep learning methods. This hybrid approach, known as deep kernel learning (DKL), combines Gaussian processes (GPs) [72] with neural networks. We perform cosmological hydrodynamical simulations for different cosmological, IGM and DAO transfer function parameters. We train a DKL emulator that produces Lyman-\(\alpha\) forest flux power spectra. We then use this emulator to place constraints in the DAO parameter space by performing a Bayesian inference with Lyman-\(\alpha\) forest data [73] collected from eleven quasar spectra from Keck-HIRES [74] and four from VLT-UVES [75].

Methods. We adopt the effective DAO transfer-function parameterization of Ref. [14]. To minimize the computational cost of hydrodynamical simulations, we train an emulator on simulations selected via the Bayesian optimization active learning procedure of Ref. [70]. The resulting emulator predicts the Lyman-\(\alpha\) forest flux power spectrum across the parameter space and is used to evaluate the likelihood, as part of Markov chain Monte Carlo sampling of the posterior.

DAO transfer function model — We model the DAO feature in the linear matter power spectrum with the same phenomenological model as Ref. [14]. The model has four free parameters: \(f\), \(A\), \(k_{\rm{peak}}\), \(k_{\rm{damp}}\). The suppression of the matter power spectrum is modeled by a mixed WDM and CDM-like transfer function [76], with the depth of the suppression set by the fraction of the dark matter not in the CDM component \(f\). The DAOs are modeled as a Gaussian-damped sinusoid whose first peak is at wavenumber \(k_{\rm{peak}}\), with amplitude \(A\) and the damping envelope set by \(k_{\rm{damp}}\). \(k_{\rm{peak}}\) also sets the wavelength of the oscillations. The bottom panel of Figure 1 shows examples of linear matter power spectra with strong DAOs (\(A\)=1) and a WDM-like model with no DAOs (\(A\)=0).

Since the simulations are expensive, it is preferable to minimize the dimensionality of the parameter space. We fix \(k_{\rm{damp}}= 10 k_{\rm{peak}}\), a conservative approximation for setting constraints, which we justify in Appendix 1, reducing our DAO model to three parameters \(\{f,A,k_{\rm{peak}}\}\).

For the DAO parameters, we assume uniform priors: \(k_{\rm peak}\sim\mathcal{U}(0.2\: h\,\mathrm{Mpc}^{-1},50\, h\,\rm Mpc^{-1})\), \(A\sim\mathcal{U}(0,1.5)\), and \(f\sim\mathcal{U}(0,1)\). The range of \(f\) spans the full physically allowed interval by construction, while the range of \(A\) encompasses values expected in viable atomic dark matter models. The upper bound on \(k_{\rm peak}\) is chosen because the Lyman-\(\alpha\) forest becomes increasingly insensitive to DAO features at smaller scales, preventing meaningful constraints beyond \(k_{\rm peak}\simeq 50\: h\,\rm Mpc^{-1}\).

Cosmological and intergalactic medium parameters —  The matter power spectrum inferred from the Lyman-\(\alpha\) forest is also sensitive to the cosmological model and the thermalization and ionization of the IGM, primarily due to the uncertain nature of reionization. To ensure our constraints are robust, we thus marginalize over cosmological and IGM parameters that capture these uncertainties.

For the cosmological parameters, we vary the spectral index \(n_s\) and amplitude \(A_s\) of the primordial power spectrum, assuming Gaussian priors with means and standard deviations \((0.9649,\,0.004)\) and \((2.101\times10^{-9},\,2.962\times10^{-11})\), respectively, derived from Planck CMB results [1]. The Lyman-\(\alpha\) forest flux power spectrum is not sensitive to other standard cosmological parameters [77], [78]. The transfer function parameters only capture the effects of the DAOs due to the coupled plasmas. We thus also vary the effective number of additional relativistic species \(\Delta N_\mathrm{eff}\) to ensure the effects of free-streaming radiation on the matter power spectrum are considered. We assume a uniform prior on the dark sector temperature ratio \(\xi \sim \mathcal{U}(0, 0.5)\), which relates to \(\Delta N_\mathrm{eff} = \frac{8}{7}(\frac{11}{4})^{4/3}\xi^4 \approx 4.4 \;\xi^4\).

To marginalize over uncertainty in the properties of the IGM, we vary the simulation input parameters \(\{H_A, H_S, z_{\rm{rei}}, T_{\rm{rei}}\}\), which modify baseline ultraviolet background (UVB) rates, taken from Ref. [79]. \(H_A\) is a multiplicative factor, while \(H_S\) applies an overdensity \(\Delta\)-dependent rescaling, such that the new photoheating rates implemented are \(\epsilon_{i} = H_A \times \epsilon_{0,i} \times \Delta^{H_S}\), where \(\epsilon_{0,i}\) are the default photoheating rates, with \(i \in\) [HI, HeI, HeII]. We additionally modify the default photoionization rates by varying the reionization redshift \(z_{\rm{rei}}\) and total heat injection \(T_{\rm{rei}}\) as implemented in Ref. [80]. Thus, for the IGM, we vary \(H_A \in\) [0.05, 2.5], \(H_S \in\) [-1, 0.7], \(T_{\rm{rei}} \in\) [\(1.5\times10^4\), \(4\times10^4\)] K, \(z_{\rm{rei}} \in\) [6, 7.8], with the maximum \(z_{\rm{rei}}\) value set by the default value in the baseline photoionization rates.

While the above parameters are varied as simulation inputs, we use the following output parameters to characterize the IGM in each simulation \(\{T_0(z = z_i),\tilde{\gamma}(z = z_i),u_0(z = z_i),\tau_0(z = z_i)\}\), for each redshift that we consider \(z_i = [4.2, 4.6, 5.0]\). The vast majority of neutral gas in the IGM follows a temperature \(T\)-density relation: \(T(z)=T_0(z)\Delta^{\tilde{\gamma}(z)-1}\). We thus describe the IGM thermal state using \(T_0(z)\), \(\tilde{\gamma}(z)\), and also \(u_0 (z)\), the cumulative thermal energy injected into the IGM per unit mass until redshift \(z\), in units of \({\rm eV}\,m_p^{-1}\). We further vary a redshift-dependent normalization \(\tau_0(z)\) of the effective optical depth, \(\tau_{\rm eff}(z)=\tau_0(z)\tau_{\rm eff}^{\rm fid}(z)\), with \(\tau_{\rm eff}\equiv -\ln\langle F\rangle\) defined by the mean transmitted flux fraction \(F\). For \(\tau_\mathrm{eff}^\mathrm{fid}\), we use the model of Ref. [73]. Following previous work [53], [56], [59], [63], [81] and to disfavor unphysically cold IGMs with an inverted temperature-density relation, Gaussian priors are adopted for the IGM parameters: \(\tau_0(z)\sim\mathcal{N}(1.0,0.05^2)\), \(\gamma(z)\sim\mathcal{N}(1.2,0.2^2)\), and \(T_0(z)\sim\mathcal{N}((8022,7651,8673)\,\mathrm{K},(3000\,\mathrm{K})^2)\) at \(z=(5.0,4.6,4.2)\), respectively. We additionally impose a uniform prior within the convex hull of the IGM simulations to prevent unphysical combinations of \(T_0\) and \(u_0\) and prevent unphysical jumps in \(T_0\) greater than 5000 K and in \(u_0\) greater than \(10\,\mathrm{eV}\,m_\mathrm{p}^{-1}\) from each redshift bin to the next.

Cosmological hydrodynamical simulations —  Our observable is the Lyman-\(\alpha\) forest 1D flux power spectrum at \(z=[4.2, 4.6, 5.0]\), i.e., the two-point line-of-sight correlation in the transmitted flux contrast in Fourier space. In order to model this quantity, we run cosmological hydrodynamical simulations of the IGM using the publicly-available code GIZMO [82]. We first produce \(\Lambda\)CDM-like linear matter power spectra using the Boltzmann code CLASS-aDM [47], which are then modified by the DAO transfer function as described in Appendix 1. The initial conditions of the simulations are generated using this modified matter power spectrum and MUSIC [83] at \(z=99\). We generate separate initial conditions for the dark matter and baryonic components [84]. We then evolve \(512^3\) particles each of dark matter and gas in a periodic \((10\, h^{-1}\) Mpc)\(^3\) box from \(z = 99\) to \(z = 4.2\), saving snapshots of particle data at \(z = [4.2, 4.6, 5.0]\). As we explain in Appendix 2, we do not need to include the hydrodynamical effects of non-minimal DM in our simulations.

To reduce significantly the computational expense of each simulation, but with negligible effect on the flux power spectrum, we implement a simplified star formation criterion (QuickLymanAlpha), following Ref. [85]: gas particles at overdensities \(>\) 1000 and with temperatures \(<\) \(10^5\) K are converted to collisionless star particles. We verify this approach in Appendix 3, where we also confirm insensitivity to the choice of feedback model.

From each particle snapshot, we generate 32000 mock quasar spectra along one axis (with spectral pixel widths \(\Delta v = 1\,\mathrm{km}\,\mathrm{s}^{-1}\)) of the Lyman-\(\alpha\) absorption line and then calculate the 1D (line-of-sight) flux power spectrum using fake_spectra [86]. These 1D flux power spectra are generated at \(z=[4.2, 4.6, 5.0]\) for each parameter point and used as training data for the emulator (below). A total number of 413 simulations were used for training. For each simulation, we vary the mean flux rescaling \(\tau_0\) and produce ten flux power spectra, which increases the size of the training dataset to 4130. We perform tests of numerical convergence with respect to particle number and box size in Appendix 4. We correct for the box size effect by rescaling the flux power spectra. The rescaling ratios used are 0.979 for \(z=5\), 0.959 for \(z=4.6\), and 0.940 for \(z=4.2\).

Deep kernel learning emulator and Bayesian optimization — We emulate the flux power spectrum as a function of the eighteen model parameters \(\boldsymbol{\theta}=\{k_{\rm{peak}},A,f,n_s,A_s,\xi,T_0(z=z_i),\tilde{\gamma}(z=z_i),u_0(z=z_i),\tau_0(z=z_i)\}\), for \(z = [4.2,4.6,5.0]\). The emulator is trained on our Lyman-\(\alpha\) forest simulations (above). Previous work has used Gaussian process emulators for Lyman-\(\alpha\) forest flux power spectra [56], [57], [59], [70], [71], [87], [88]. A GP emulator interpolates simulation outputs by modeling their covariance through a kernel function [72]. Standard kernels assume stationarity, i.e., that correlations depend only on the separation between points in parameter space. However, the response of the Lyman-\(\alpha\) forest flux power spectrum to DAO parameters is highly nonstationary.

When \(f=0\), there are no DAO effects, regardless of \(k_{\rm{peak}}\) and \(A\). As \(f\) increases, \(k_{\rm{peak}}\) and \(A\) begin to influence the flux power spectrum, although their effects are still highly dependent on \(f\). We thus use deep kernel learning (DKL) emulation [89], as implemented in GPyTorch. The DKL emulator combines a neural network (NN) with a GP, using the NN as a feature extractor to map the original parameter space into a latent representation where a stationary kernel can accurately model the covariance structure.

The NN feature extractor is physically informed, guiding the learned latent representation toward physically relevant parameter combinations and improving emulator performance. The feature extractor and GP are trained jointly until convergence, defined by the change in the NN hyperparameters falling below a prescribed threshold. Then, we freeze the NN hyperparameters and fine-tune the GP kernel with the evidence lower bound (ELBO) as the loss function [90]. The DKL emulator allows us to represent the non-stationarity and degeneracies in our parameters while retaining the predictive variance from GP models. Because a full GP scales poorly with the size of the training set, retraining the GP during optimization of the feature extractor would be computationally prohibitive. We therefore employ a variational GP, which represents the latent space covariance through a small set of learnable inducing points. These inducing points act as a compressed representation of the training data and are optimized jointly with the neural network and GP. For a detailed explanation, see Appendix 5.

The initial emulator was constructed using 68 parameter points arranged in a Latin hypercube. From this initial set, we used Bayesian optimization (as previously used by Refs. [56], [59], [70]) to iteratively add simulations to the training set. This approach uses observed data and emulator uncertainty quantification from the GP to decide the optimal construction of the training set. We ran 413 simulations in total, and stopped adding training simulations once the estimated posterior distribution converged with respect to the training set (see Appendix 5 for details).

Posterior distribution sampling — We compare the emulator output to the flux power spectrum measured in Ref. [73], derived from eleven quasar spectra from Keck-HIRES [74] and four from VLT-UVES [75], which contains the smallest scales measured to-date in the Lyman-\(\alpha\) forest (velocity wavenumbers \(k_f < 0.2\,\mathrm{s}\,\mathrm{km}^{-1}\); see Fig. 1). We assume a Gaussian likelihood function which includes emulator uncertainty. We use Markov chain Monte Carlo sampling to estimate the posterior distribution using the emcee sampler [91], declaring convergence once each chain is fifty times the auto-correlation length.

Results and discussion. The lower panel of Fig. 1 shows the ratio of the linear matter power spectra of models with strong DAOs (\(f=1\), \(k_{\rm{peak}}=20\:h\)/Mpc, \(A=1\)) and WDM-like suppression (\(f=1\), \(k_{\rm{peak}}=20\:h\)/Mpc, \(A=0\)) to the corresponding \(\Lambda\)CDM model (\(f=0\)), with all other cosmological and IGM parameters fixed. The matter power spectrum ratio equals the square of the transfer function \(T(k)\) (see Appendix 1), so the first (negative) DAO trough in \(T(k)\) becomes the first DAO peak in the power spectrum, while the second DAO peak is at \(k_{\rm{peak}}\). The middle panel of Fig. 1 shows the imprint of these two models on the flux power spectrum at the three redshifts we consider. As previously seen, a small-scale WDM-like linear matter power spectrum suppression causes a strong small-scale suppression in the flux power spectrum. We demonstrate that the addition of a strong DAO, however, significantly changes the flux power. In particular, the oscillatory feature largely washes out, but it does reinstate flux power at all scales below the initial cut-off. We thus anticipate that the Lyman-\(\alpha\) forest will be sensitive to DAO features beyond simple suppressions and enhancements.

The upper panel in Fig. 1 compares emulated flux power spectra to observations. The dot-dashed lines show the same strong DAO model as above, illustrating that such features at the scales considered are strongly disfavored. The dashed lines show the emulated flux power spectra at the maximum posterior parameter point, while the solid lines show those for the best-fit point. These points are not identical owing to prior volume effects. The fit is in general good, apart from the final wavenumber bin, where the flux power is underestimated. This result is consistent with previous analyses [55], who attributed this feature to either the noise model or a signature of power enhancement beyond \(\Lambda\)CDM. We investigate this feature further in Appendix 7, finding that it does not affect our conclusions since the dominant DAO sensitivity derives from the largest affected scales.

Figure 2 shows the 2D posterior distribution of the DAO parameters \(k_{\rm{peak}}\), \(A\) and \(f\) given observed Lyman-\(\alpha\) forest flux power spectra after marginalizing over all cosmological and IGM parameters, compared to the parameter space previously allowed by the CMB. The \(k_{\rm{peak}}\)-\(A\) panel demonstrates that Lyman-\(\alpha\) forest data are indeed sensitive to the presence of DAOs. The data disfavor strong DAOs (larger \(A\)) when \(k_{\rm{peak}}\) is smaller. The \(k_{\rm{peak}}\) - \(f\) panel demonstrates that the Lyman-\(\alpha\) forest limits the abundance of DM forming DAOs to be \(f \lesssim 0.2\) when \(k_{\rm{peak}}\sim 0.2\,h\,\mathrm{Mpc}^{-1}\), but that this limit relaxes to \(f \lesssim 0.3\) when \(k_{\rm{peak}}\sim 50\,h\,\mathrm{Mpc}^{-1}\). Otherwise, we find that the data are consistent with the \(\Lambda\)CDM limit (\(f = 0\)). The uniform prior on \(k_{\rm{peak}}\) means that we do not sample the other \(\Lambda\)CDM limit when \(k_{\rm{peak}}\to \infty\); however, we find that these data lose sensitivity to DAOs (i.e., the parameter \(A\)) already at the maximum \(k_{\rm{peak}}\) that we consider. We discuss parameter degeneracies between DAO and cosmological and IGM parameters in Appendix 6. The 1D marginalized 95% credible intervals are given in Table 1 in Appendix 6.

Studies of the impact of strong DAOs on the CMB [6], [15], [47] have worked directly in the parameter space of atomic dark matter, yielding constraints that are most clearly expressed in terms of the fraction of dark matter that is atomic, the dark photon temperature, and the dark sound horizon. Translated into the parameters of our model, the most recent constraints with Planck and ACT DR6 data [48] reach down to the \(f\leq0.05\) level for \(k_{\rm{peak}}\approx 1.2\; h/\mathrm{Mpc}\), but lose sensitivity past \(k_{\rm{peak}}\approx 4.5\;h/\mathrm{Mpc}\) (indicated by the red line in Fig. 2). Lyman-\(\alpha\) forest data are sensitive to significantly smaller scales than current CMB data, extending sensitivity to DAOs from \(\sim 1\,h\,\mathrm{Mpc}^{-1}\) to \(\sim 50\,h\,\mathrm{Mpc}^{-1}\). In Ref. [14], UVLF data are used to constrain DAOs, with the same transfer function model that we use. Bounds on \(f\) and \(k_{\rm{peak}}\) are derived, but, unlike this work, the DAO amplitude \(A\) was not constrained. While the use of different priors complicates further quantitative comparison, our analysis appears to provide stronger bounds at the smallest scales.

Our baseline model assumes that the same time-varying UVB rates apply at all locations, but this is a known approximation [92]–[94]. We thus re-perform the inference using a model of a spatially-inhomogeneous reionization calibrated by radiative transfer simulations [95], finding that our results do not change (see Appendix 7). While our results apply to generic models of interacting dark matter that produce DAOs, we additionally verify that they are directly applicable to models of atomic dark matter. In this case, additional dark radiative cooling can affect astrophysical structures [16]–[29], but we find that this extra physics has negligible influence on the Lyman-\(\alpha\) forest (see Appendix 2). This test means that we can directly combine these results with other analyses to probe atomic dark matter, which we will consider in future work.

Conclusion and outlook. In this work, we present the first Lyman-\(\alpha\) forest constraints on DAOs (and indeed any features beyond simple suppression and enhancement) in the linear matter power spectrum using a full statistical analysis that models the flux power spectrum using cosmological hydrodynamical simulations. Given the high cost of simulations and the non-trivial, non-stationary covariance of DAO model parameters, we use a novel combination of deep kernel learning and Bayesian optimization to construct an emulator of the flux power spectra trained on simulations. We use the emulator to perform inference on DAO, cosmological and IGM parameters given the smallest-scale Lyman-\(\alpha\) forest data. We find that the Lyman-\(\alpha\) forest limits the fraction of dark matter forming DAOs \(f \lesssim 0.3\), with stronger bounds \(f \lesssim 0.2\) when \(k_{\rm{peak}}\sim 1\,h\,\mathrm{Mpc}^{-1}\). These results constitute the first small-scale (\(k_{\rm{peak}}\gtrsim 1\,h\,\mathrm{Mpc}^{-1}\)) cosmological constraints on the DAO amplitude using a full forward model from initial conditions to the IGM. Future work will combine our results with other cosmological datasets [14], [15] to probe concrete physics models like atomic dark matter. We further anticipate that our results can apply to other models like millicharged DM [96]–[99], DM with an electric dipole moment [100], DM with massive boson exchange [39], DM interacting with massless sterile neutrinos via a broken dark \(U(1)\) gauge symmetry [101]–[103], DM charged under a non-Abelian gauge symmetry [104], and to searches for inflationary potential features [105].

Acknowledgements↩︎

The authors thank Elisa Ferreira for helpful discussions. This work was enabled by computational resources provided by Compute Canada and the Digital Research Alliance of Canada. The work of ZY, DC and NM was in part supported by Discovery Grants from the Natural Sciences and Engineering Research Council of Canada, the Canada Research Chair program, the Ontario Early Researcher Award, and the University of Toronto McLean Award. KKR is supported by an Ernest Rutherford Fellowship from the UKRI Science and Technology Facilities Council (grant no. ST/Z510191/1). The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. JB acknowledges support from NSF grants PHY-2210533 and PHY-2513893. SR acknowledges support from the Eric and Wendy Schmidt AI in Science Fellowship.

1 Dark acoustic oscillation transfer function↩︎

The linear matter power spectrum in the DAO model is defined by (the square of) a transfer function \(T(k)\) multiplying the \(\Lambda \text{CDM}\) linear matter power spectrum such that \(P_\mathrm{DAO}(k,z) \equiv T^{2}(k)P_{\Lambda \text{CDM}}(k,z)\). The transfer function is the sum of two terms \(T=T_{\alpha\beta\gamma\delta} + T_{\mathrm{osc}}\). The first term accounts for the WDM-like power spectrum suppression at high \(k\) [76], [106], [107]: \[T_{\alpha \beta \gamma \delta}(k) \equiv f(1 + (\alpha k)^\beta)^\gamma + (1-f), \label{model:fabgd}\tag{1}\] where \(\alpha\), \(\beta\) and \(\gamma\) are free parameters that we will calibrate below. \(f\) is the fraction of the total dark matter that forms DAOs; the remainder is cold dark matter. The second term \(T_\mathrm{osc}\) accounts for the damped DAOs. The oscillations are modeled as sinusoidal with frequency \(\omega\) and amplitude \(A\), beginning at wavenumber \(k_{\mathrm{start}}\), peaking at wavenumber \(k_{\rm{peak}}\) and damped at wavenumber \(k_\mathrm{damp}\): \[\begin{align} T_{\mathrm{osc}}(k) \equiv \Theta(k - k_{\mathrm{start}}) \times \left[fA \cos\left(\omega \left(\frac{k}{k_{\rm{peak}}} - 1\right)\right)e^{-\left(\frac{k}{k_{\rm{damp}}}\right)^{2}}\right] \,. \end{align} \label{model:osc}\tag{2}\]

To ensure that our model matches realistic power spectra with strong DAOs, and with future atomic dark matter analyses in mind, we calibrate some of the free parameters by fitting to linear matter power spectra computed using the modified Boltzmann code CLASS-aDM [47]. However, we stress that our parameterization choice nonetheless can be mapped to a wider range of dark sector models with DAOs. We find \(\beta = 4.15\), \(\gamma = -20\), \(\omega=2.083\pi\), \(k_{\mathrm{start}}=0.28k_{\rm{peak}}\). The parameter \(\alpha\) is chosen such that \(T_{\alpha \beta \gamma \delta}(k_{\textrm{start}}) = 0.1 f + (1-f)\), so that the power spectrum suppression matches onto the start of the DAOs: \[\alpha = \frac{1}{k_{\mathrm{start}}}\left(0.1^{1/\gamma}-1\right)^{1/\beta} \,.\]

After this calibration, we have four free parameters \(\{k_{\rm{peak}},k_{\rm{damp}},f,A\}\). To reduce the dimensionality of the parameter space further, we choose to fix \(k_{\rm{damp}}= 10 k_{\rm{peak}}\). This decision is motivated by the fact that a larger damping wavenumber (i.e., a less damped oscillation) leads to more power being returned to the power spectrum. Thus, the flux power spectra will be more CDM-like and the constraints we find on DAOs will be necessarily conservative. Fig. 3 illustrates this fact and also that the flux power spectrum is in any case only weakly sensitive to \(k_{\rm{damp}}\).

2 Tests of dark radiative cooling in the atomic dark matter model↩︎

Dark acoustic oscillations arise in many models of interacting dark matter and dark radiation, but one model of particular interest is atomic dark matter (aDM) [5]. In addition to forming DAOs, aDM can undergo radiative and collisional cooling, which can also affect astrophysical structure [108], [109]. While DAOs occur at much earlier redshifts and are accounted for in CLASS, the cooling effects become important at much later redshifts and require a specialized version of GIZMO [16] to incorporate. In this implementation, aDM is treated hydrodynamically as a gas, similar to baryons in the original version, making it more computationally expensive. For the DAO results we present here to apply to models of aDM, we must verify that aDM cooling has a negligible effect on the Lyman-\(\alpha\) forest.

We hypothesize that since the Lyman-\(\alpha\) forest signal derives from the IGM, and that aDM cooling occurs in the centers of halos, the cooling will have little effect on the flux power spectra. We explicitly test this approximation by considering an aDM parameter point expected to cool significantly, with a large dark matter fraction in aDM (\(f = 0.9\)), a large dark temperature ratio (\(\xi = 0.5\)), and a small \(\beta_{\rm{cool}}\) (\(\sim0.001\))¹ as defined in Ref. [18]. At this parameter point, we run one simulation with the full aDM hydrodynamics, and one CDM+baryons-only simulation but where we initialize the CDM particles to have the same matter power spectrum as the aDM simulation at \(z=99\). Fig. 4 shows that we indeed find no appreciable effect on the flux power spectrum from aDM cooling.

However, while aDM cooling does not have a direct effect on the Lyman-\(\alpha\) forest through the distribution of matter in the IGM, cooling in halos could have a large impact on star formation rates and thus the UV background photons that heat and ionize the IGM. Since we already consider a large range of IGM histories by varying the input simulation parameters, \(\{H_A, H_S,z_{\rm{rei}},T_{\rm{rei}}\}\), we anticipate that this aDM effect would be degenerate with the output IGM parameters we already consider.

3 Tests of QuickLymanAlpha flag and feedback models↩︎

Previous Lyman-\(\alpha\) forest studies have used a simplified star formation criterion to speed up simulations while having negligible impact on the flux power spectra [85]. Usually referred to as the QuickLymanAlpha flag, the criterion converts gas particles at overdensities \(>\) 1000 and with temperatures \(<\) \(10^5\) K into collisionless star particles. Fig. 5 confirms that there is no significant difference between the GIZMO simulation with and without the QuickLymanAlpha implementation (blue line). Further, studies have shown that active galactic nuclei (AGN) feedback can have a marginal effect on the 1D flux power spectra [110]–[112]. To explore this effect, we compare the effect of the FIRE-2 feedback model [113] to the default setting (green line). While the FIRE-2 module results in a modest \(\sim(5-10)\%\) difference at smaller scales, it is within the error bars of the data. Thus, for the sake of computational speed, we opt to use the QuickLymanAlpha flag and neglect using FIRE-2.

4 Numerical convergence tests↩︎

To ensure that we have chosen a sufficient number of simulation particles, we run test simulations with different particle resolutions for a fixed box size. We run this test at a parameter point with strong matter power spectrum suppression, so that any numerical effects from small-scale fragmentation would be most apparent. Fig. 6 shows the convergence tests for particle resolution. The ratio between the flux power spectrum from the \(512^3\) particle simulation to that from the \(1024^3\) particle simulation remains within the data error.

In optically-thin Lyman-\(\alpha\) forest simulations, the finite box size introduces a systematic bias in the predicted flux power spectra due to the absence of long-wavelength modes and their nonlinear coupling to smaller scales [114], [115]. The missing large-scale modes lead to smaller bulk flow velocities and less shock heating in the IGM, thereby yielding systematically colder gas. These effects reduce thermal broadening and Jeans smoothing, enhancing small-scale flux power. The power enhancement is expected to increase monotonically with smaller box size, as seen in Refs. [114], [115]. We confirm this monotonicity with simulations with box sizes of \(5\,h^{-1}\) Mpc, \(10\,h^{-1}\) Mpc, and \(20\,h^{-1}\) Mpc, as shown in Fig. 7. To reduce the effects of sample variance, we average over ten simulations with different random seeds for the \(5\,h^{-1}\) Mpc and \(10\,h^{-1}\) Mpc boxes. We correct for the box size effect by rescaling the flux power spectra of our \(10\,h^{-1}\) Mpc simulations to match the larger boxes, as described in the main text.

5 Deep kernel learning emulator and Bayesian optimization↩︎

Gaussian processes (GPs) have been used previously to emulate the Lyman-\(\alpha\) forest flux power spectrum [56], [57], [59], [70], [71], [87], [88]. A GP model describes a prior distribution over functions \(g(x)\) which can be specified by a mean function \(m(x)\) and a covariance (kernel) function \(\mathrm{ker}(x,x')\). In our context here, \(x\) is the DAO, cosmological and IGM parameter vector and \(g\) is the mapping to the flux power spectrum vector evaluated at fixed values of the velocity wavenumber \(k_\mathrm{f}\) and redshift \(z\). We are free to select different kernel functions to encode desired properties of the model. If we condition this Gaussian process prior on a training dataset of simulations, we obtain a posterior distribution that allows us to make predictions of the flux power spectrum at other parameter points. Given a simulated training dataset of pairs of parameter points and flux power spectra \({X,Y}\) with noise \(\sigma_n\), the covariance matrix becomes

\[K_Y = \mathrm{ker}(X,X)+\sigma_n^2I,\] where \(I\) is the identity matrix. If we want to make a prediction at parameter point \(x'\), the joint prior distribution is \[\left[\begin{array}{c} \mathbf{Y} \\ g(x') \end{array}\right] \sim \mathcal{N}\left(\left[\begin{array}{c} m \\ m' \end{array}\right],\left[\begin{array}{cc} K_Y & \mathrm{ker}(X,x') \\ \mathrm{ker}(X,x')^T & \mathrm{ker}(x', x') \end{array}\right]\right).\] Conditioning on the simulated training data, we get a Gaussian-distributed predictive posterior \(p(g(x')|Y)\) given by the posterior mean \[\mu'=m'+\mathrm{ker}(X,x')^TK_Y^{-1}(Y-m)\] and posterior variance \[\sigma'^2=\mathrm{ker}(x',x')-\mathrm{ker}(X,x')^TK_Y^{-1}\mathrm{ker}(X,x'). \label{eq:GP95variance}\tag{3}\] This posterior is the emulator prediction for the flux power spectrum at a new DAO, cosmological and IGM parameter point.

However, despite its previous success for Lyman-\(\alpha\) forest flux power spectrum emulation, a traditional GP emulator is not sufficient to describe the full complexity of the DAO effects on the Lyman-\(\alpha\) forest flux power spectra, even by considering different covariance kernel function choices. Specifically, the covariance of our parameter space is highly non-stationary. As \(f \rightarrow 0\), the effects of \(k_{\rm{peak}}\) and \(A\) on the flux power vanish, i.e., their covariance changes. There do exist non-stationary covariance kernels such as the Gibbs kernel but we leave detailed comparisons to future work.

To represent fully these behaviors, we turn to deep kernel learning (DKL) [89], where a small neural network works as a feature extractor and maps the original parameter space to a latent parameter space, where we enforce stationarity. We then use a GP to emulate the flux power spectrum given this new latent parameter space. Other emulation strategies use a neural network to map directly from parameters to simulation outputs like the power spectrum [116]–[118]. However, these do not typically provide a global uncertainty quantification (UQ). Bayesian neural networks [119], [120] can combine deep learning emulation with global UQ, but also typically suffer from underdetermined hyperparameters. UQ is important in our case as we will use posterior predictive uncertainties to design the simulation training set by active learning (Bayesian optimization).

Figure 8 shows a flow chart that illustrates the architecture of our DKL emulator. The feature extractor network is at the top. Its architecture is physically informed. The neural network is split into two branches: one for DAO transfer function parameters and the other for cosmological and IGM parameters. In the transfer function branch, we first multiply a soft gate function \(G(f)\) to \(k_{\rm{peak}}\) and \(A\) to enforce the \(\Lambda\)CDM behavior at \(f=0\), where \(k_{\rm{peak}}\) and \(A\) do not change the transfer function:

\[G(f)=1 - e^{-(f/\tau)^p},\] where hyperparameters \(\tau=0.02\) and \(p=3\) are fixed to ensure that the gate function turns on smoothly as \(f\) increases. The network is then followed by two linear layers with widths 16 and 8, respectively, each followed by a Softplus activation. An additional feature-wise linear modulation (FiLM) layer [121] is applied to the latent representation of the transfer function branch, allowing the learned features to be conditioned continuously on the DAO parameter \(f\). The FiLM layer thus applies a transformation \[h' = \gamma(f)h+\beta(f),\] where \(h\) and \(h'\) are the layer before and after modulation, \(\gamma\) and \(\beta\) are functions of \(f\) learned through a small network. They are optimized such that \(h'=h\) when \(f=0\), corresponding to the desired \(\Lambda\)CDM behavior. This modulation thus smoothly encodes the physical behavior of \(f\) so that the total network does not have to “jump" between drastically different behaviors when approaching the \(\Lambda\)CDM limit. The cosmological and IGM parameter branch consists of a single 16-dimensional fully connected layer. The 8-dimensional physics representation and the 16-dimensional nuisance representation are concatenated into a 24-dimensional feature vector. This vector is subsequently transformed through a linear layer followed by a Softplus activation, yielding 10 output latent dimensions, which matches the input dimensionality.²

For the second stage of the DKL (see bottom of Fig. 8), we use a Gaussian process emulator that now takes pairs of parameters in the latent space and corresponding simulated flux power spectra and then outputs predicted flux power spectra at new parameter points with uncertainties quantified. We use a variational Gaussian process, also called a sparse Gaussian process [72], [90], [122]. The main advantage of the variational GP over the traditional GP described above is its scalability with large training sets. In this work, the training set has 4130 points (including accounting for the mean flux rescaling). The variational GP has much smaller memory requirements and faster evaluation time since it avoids inverting the matrix \(K_Y\) for the total training set. This scalability is necessary for DKL, because the memory requirement and evaluation time quickly becomes substantial if each training step requires the evaluation of a deep neural network and an exact GP for the full dataset. Variational GP thus naturally supports batch training, which makes it much more compatible with neural networks.

We now describe the variational GP approach. For a training dataset, let \(X'\) be the set of \(n\) parameter points and \(Y\) be the noisy evaluation of some latent function \(g(X')\). For clarity, as above, \(X'\) is the set of training DAO, cosmological and IGM parameters, now mapped to the new latent space, \(g\) is the mapping to simulated flux power spectra, at fixed values of \(k_f\) and \(z\), and \(Y\) is the set of flux power spectra evaluated at the training points. The noise here refers to sample variance in the simulations. As in the case of the exact GP emulator described above, we introduce a Gaussian process prior to \(g\), specified by a mean function and a covariance kernel, and we want to obtain the predictive posterior. For the variational GP, \(m\) “inducing points” are chosen in the same (latent) space of \(X'\), but at different points than the training set. Let the set of inducing points be \(Z\). The flux power spectrum values \(u \equiv g(Z)\) evaluated at each inducing point are treated as a compact set of variables that summarizes the behavior of the full GP. Unlike the training targets, the inducing values are not directly observed. Instead, their posterior distribution is learned during training. The exact predictive posterior is thus approximated as \[p(g|X',Y)\approx q(g)=\int p(g|u)q(u) \mathop{}\!\mathrm{d}u,\] where \(q(u)\) is the prior distribution for the inducing points. On the covariance level, this means we can approximate the true covariance \(K_{n n}\approx K_{n m} K_{m m}^{-1} K_{m n}\) where \(K_{m m}\) is the covariance matrix on the inducing points, and \(K_{n m}\) is the cross-covariance matrix between training and inducing points. The computation now scales with \(\mathcal{O}\left(nm^2\right)\) instead of \(\mathcal{O}\left(n^3\right)\). In our case, we have 4130 training points and only 128 inducing points.

The inducing points are selected by minimizing the Kullback-Leibler (KL) divergence between the approximate posterior \(q\) using the inducing points and the exact posterior \(p\) using the training points: \[KL(q(g)||p(g|X',Y))=\int \mathop{}\!\mathrm{d}f q(g) \log {\frac{q(g)}{p(g|X',Y)}}. \label{eq:KL}\tag{4}\] However, we want to avoid directly computing the exact posterior \(p\), which is the motivation of using the variational GP. To minimize the KL divergence without such an explicit computation, we use Bayes’ theorem such that \[\log{p(g|X',Y)}=\log{p(Y|g,X')}+\log{p(g)}-\log{p(Y|X')}. \label{eq:bayes}\tag{5}\] It follows by combining Eqs. (4 ) and (5 ) that \[\log{p(Y|X')}-KL(q(g)||p(g|Y))=\int \mathop{}\!\mathrm{d}g q(g)\log{p(Y|g)}-KL(q(g)||p(g)) \equiv \mathrm{ELBO}.\] We define the above to be the evidence lower bound (ELBO), which we maximize during training. From the left hand side, we can see that maximizing ELBO simultaneously maximizes the marginal likelihood \(p(Y|X')\), which improves the fit to the simulated data, and minimizes the KL divergence to the true posterior, which optimizes inducing point selection. The right hand side is directly used to compute ELBO during each training step. In our architechture, the feature extractor and the GP are jointly trained for 150 epochs before we freeze the neural network hyperparameters and we then fine-tune the GP, using the ELBO [90] as the loss function. The model is trained using the Adam optimizer. We illustrate the projection of training and inducing points onto two of the latent parameters in Fig. 8. As expected, the density of inducing points roughly traces the density of training points. The streaks in the latent space projection come from the mean flux rescalings that more densely sample the \(\tau_0\) dimensions.

In summary (see also Fig. 8), our input parameters are fed into a feature extractor that consists of a neural network split into two branches, where the physically informed branch contains a soft gate and a layer of FiLM modulation to express known parameter dependencies. The feature extractor produces latent parameters, which are then fed into a variational GP model. The GP model consists of a posterior distribution approximated using a set of inducing points in the latent parameter space, which are selected during training to maximally represent training data. We use this approximate posterior distribution to predict the flux power spectrum. We iteratively expand the training set using the Bayesian optimization procedure described in Ref. [70]. At each Bayesian optimization step, we re-train the DKL emulator using the procedure described above. In practice, we train three emulators, one for the flux power spectrum at each redshift that we consider.

Figure 9 shows the DAO, cosmological and IGM model parameters’ shift by number of sigmas between each Bayesian optimization training set iteration. We start with 68 simulations sampled by a Latin hypercube in the parameter space (680 training points in total after mean flux rescaling post-processing). Bayesian optimization adds training simulations by a balance of selecting points that have high emulator prediction uncertainty (exploration) and points that have high posterior probability given the true data (exploitation). The simulation batch size that we add at each step varies from 10 to 5, with an exception in the last few fine-tuning batches, which only contain 2-3 simulations each. The final training set contains 413 simulations (4130 training points in total), now clustered in regions of high posterior probability given the flux power spectrum data. We observe that the variation between the posteriors given successive training set iterations is small by the end, indicating an emulator that is stable with respect to the training set design. Indeed, most model parameters have already converged by batch 16, but the DAO parameters and the IGM parameters \(u_0(z=z_i)\) are not stable until later.

Figures 10 and 11 show kernel density estimates of the distributions of the ratio of empirical emulator error to estimated emulator error and data error, respectively. Empirical emulator error is defined as the leave-one-out cross-validation error, i.e., re-training the emulator after leaving out each of the training simulations³ and comparing the flux power prediction \(P_\mathrm{f,pred}\) to the true simulated flux power \(P_\mathrm{f,sim}\). Estimated emulator error is the standard deviation of the posterior predictive distribution of the flux power spectrum given the Gaussian process emulator in the leave-one-out cross-validation setting, i.e., \(\sigma\) as defined in Eq. 3 . Data error is as estimated in Ref. [73]. We find that the emulator error is consistently overestimated with respect to the true error at all wavenumbers and redshifts (the blue distributions are more peaked than the orange), i.e., the emulator is underfitting. After propagation of the emulator uncertainty into the data likelihood function, this will lead to more conservative model parameter constraints. We find that for the low-\(k_\mathrm{f}\) and mid-\(k_\mathrm{f}\) bins, the emulator error is usually much less than the data uncertainty. For the high-\(k_\mathrm{f}\) bins, the uncertainties are comparable, which, again after considering propagation to the likelihood, will lead to more conservative bounds. This result is consistent with previous flux power spectrum emulators [123], where emulator uncertainties contribute non-negligibly to the total error budget. We leave to future work an investigation into strategies to reduce further this part of the error budget, e.g., by using multi-fidelity emulators [124], [125] to distribute training simulations even more optimally.

6 Degeneracies between dark acoustic oscillation and intergalactic medium parameters↩︎

Figure 12 shows 2D marginalized posterior distribution degeneracies between DAO transfer function and IGM parameters at \(z = 4.2\) (equivalent results are found at the other redshifts). The degeneracies are non-trivial since we have DAO parameters that both suppress and enhance the flux power spectrum in a scale-dependent way. As \(A\) increases, this increases small-scale flux power (see Fig. 1). This increase can be compensated by extra IGM heating from the UVB background that suppresses the small-scale power (making the IGM more diffuse), thus leading to the observed degeneracy between \(A\) and \(u_0\). We also find that hotter IGMs (higher \(u_0\)) are preferred when \(k_{\rm{peak}}\) decreases. This correlation occurs since the marginalized posterior is integrated over \(A\) (and all other parameters). The power boosting effect of \(A\) is stronger for lower \(k_{\rm{peak}}\) meaning that the increase in power is preferentially compensated by suppressing the power through extra IGM heating. This degeneracy diminishes as \(k_{\rm{peak}}\) increases, since the DAO peak in the linear transfer function (that causes the boost in flux power) gets pushed to higher \(k\) beyond the sensitivity of the data. Otherwise, the IGM posterior distributions are broadly consistent with the literature [54]–[56], [59], although we caution against a direct comparison to previous results given the extra dark matter parameters that we consider here which opens up new parameter degeneracies. The primordial power spectrum and effective optical depth posteriors are consistent with the prior distributions (see Table 1).

7 Tests of patchy reionization and noise model↩︎

In our simulations, we assume that reionization is a spatially-homogeneous process, i.e., the same (time-dependent) photoionization and photoheating rates are applied at each spatial position in the box. It is known that this is an approximation and that reionization proceeds by expanding ionizing bubbles around sources of ionizing sources, leading to large-scale (\(\sim 40\,\mathrm{Mpc}\)) spatial inhomogeneities during and for a time after reionization in the ionization and temperature fields [94], [95]. The scale of these inhomogeneities is much larger than the scale of the DAOs we test in this work. However, there is a secondary effect where the local mean flux will spatially fluctuate, leading to a \(\sim 10\%\) effect on the smallest scales we probe in the flux power spectrum as seen in radiative transfer simulations [92], [93]. This effect however is not significant compared to current data uncertainty. We explicitly test this effect by correcting the flux power spectra from the emulator according to the model presented in Ref. [95] that was calibrated to radiative transfer simulations. We find that, after applying this correction, there is no statistically significant change in the posterior distribution (see Fig. 13) and that our results are robust to this effect given current statistical uncertainties. We anticipate, however, that the consideration of such effects will become increasingly important as the number of high-quality, high-resolution, high-redshift spectra increases, e.g., from the onset of extremely large telescopes [126].

As commented in the main letter, even the \(\Lambda\)CDM limit of the emulator does not return a good fit to the largest wavenumber bin of the data. This discrepancy was previously noted in Ref. [55]. They argue that this may arise from a mis-modeling of the noise in the data in the original data reduction [73], although this statement cannot be disentangled from any potential simulation systematics. To account for this effect, they perform a test where they add a white noise term to the emulated flux power spectrum with a free amplitude parameter \(a_\mathrm{n}(z = z_i)\) that is varied. A different parameter for each redshift bin is allowed. They find that the posteriors on \(a_\mathrm{n}(z = z_i)\) mildly deviate from the no noise mis-modeling limit \(a_\mathrm{n}(z = z_i)=0\). We repeat this test in this analysis (see Fig. 14). While we find that \(a_\mathrm{n}(z = z_i)\) posteriors peak near unity, they remain consistent with the no mis-modeling limit also. We attribute the difference between the two noise analyses as to the effect of emulator uncertainty reducing the constraining power of the largest wavenumber data bin. Indeed, we do find that the fit to the data is improved with the addition of the \(a_\mathrm{n}\) parameters. Most importantly, we find that the DAO posteriors are very insensitive to the noise model. As discussed elsewhere, most of the constraining power on the DAO parameters comes from the larger affected scales in these data, where the effect of the DAO amplitude \(A\) is strongest.

References↩︎