Compressing the Cosmological Information in One-dimensional Correlations of the Lyman-α Forest

We begin by describing the simulations used in the analysis, which fall into two categories. First, a set of training simulations are used to construct the emulator. Second, a small number of test simulations are run to represent mock P_1D measurements in a variety of different cosmologies, and which are used to test and validate our analysis pipeline. Both the training and most of the test simulations were presented in Pedersen et al. (2021), where we also described and tested the emulation framework. Here we give an overview of the simulations, and refer the reader to Pedersen et al. (2021) for a more detailed description.

The simulations were run in MP-Gadget ¹⁰ (Feng et al. 2018), a TreeSPH code based on Gadget-2 (Springel 2005). All simulation boxes had a size of L = 67.5 Mpc and 768³ gas and CDM particles. The initial conditions were generated with MP-GenIC at z = 99, with Fourier modes that had random phases but fixed initial amplitudes (Angulo & Pontzen 2016; Villaescusa-Navarro et al. 2018; Anderson et al. 2019; Pedersen et al. 2021). In order to further reduce cosmic variance, for each model (both in the training and test sets) we ran a pair of simulations with inverted phases, and each quantity estimated from the simulations is taken as the average of the pair.

We output 11 snapshots, equally spaced in redshift between z = 2 and z = 4.5. To produce mock Lyα forest spectra from each snapshot, we use fake_spectra ¹¹ (Bird 2017) to calculate a two-dimensional (2D) grid of 500² transmission skewers from each snapshot, with a line-of-sight resolution of 0.05 Mpc. The cosmological and astrophysical parameters used in both the training and test simulations are listed in Table 1.

Table 1. Cosmological and Astrophysical Parameters for the Training and Test Simulations

	Training Set	Central	Neutrino	Running
As (× 10−9)	[1.35–2.71]	2.006	2.251	2.114
ns	[0.92–1.02]	0.9676	0.9676	0.9280
αs	0.0	0.0	0.0	0.015
Ωm	0.316	0.316	0.324	0.316
Σmν (eV)	0.0	0.0	0.3	0.0
${{\rm{\Delta }}}_{p}^{2}(z=3)$	[0.25–0.45]	0.35
np (z = 3)	−[2.35–2.25]	−2.30
zrei	[5.5–15]	10.5
HA	[0.5–1.5]	1.0
HS	[0.5–1.5]	1.0

Note. The limits of the Latin hypercube for the training simulations are shown in the left column, where only the primordial power spectrum and astrophysical parameters are varied. The primordial parameters A_s and n_s here are defined at the CMB pivot scale of k = 0.05 Mpc⁻¹. The Central, Neutrino, and Running simulations are constructed such that they have the same small-scale linear matter power spectrum (${{\rm{\Delta }}}_{p}^{2}$ and n_p ) at z = 3. For all simulations, we fix ω_c = 0.12, ω_b = 0.022, and h = 0.67.

Download table as:  ASCIITypeset image

We provide a brief overview of the emulator parameters and framework. We use the LaCE ¹² framework presented in Pedersen et al. (2021), and refer the reader to this reference for a more complete description. LaCE uses a Gaussian process emulator ¹³ to predict P_1D as a function of six parameters: the dimensionless amplitude (${{\rm{\Delta }}}_{p}^{2}$) and slope (n_p ) of the linear power spectrum around a pivot scale of k_p = 0.7 Mpc⁻¹; the mean transmitted flux fraction (or mean flux, $\bar{F}$); a thermal broadening scale defined in comoving units (${\sigma }_{T}^{\mathrm{com}}$), set by the temperature of the gas at mean density; the slope of the temperature–density relation (γ); and the filtering length in inverse comoving units (${k}_{F}^{\mathrm{com}}$), a proxy for gas pressure. Note that while the P_1D is naturally observed in velocity units, the above quantities are all defined in comoving units in the emulator. The motivation for this is so that the simulated P_1D is estimated at a fixed set of wavenumbers for snapshots at all redshifts.

We train the emulator using 30 pairs of simulations, described in Table 1. These simulations explore different thermal and reionization histories by varying z_rei, H_A , and H_s , which control the redshift of reionization and the heating rates of the gas as rescalings around a fiducial model from Haardt & Madau (2012; see Pedersen et al. 2021 for more detail). The simulations have different amplitudes and slopes of the primordial power spectrum (A_s , n_s ), but have the same value for the physical densities of CDM (ω_c = Ω_c h²) and baryons (ω_b = Ω_b h²), the same value of H₀, and do not include massive neutrinos. All 11 snapshots from all 30 models are used simultaneously, for a total of 330 points in the training sample.

In the implementation of the emulator presented in Pedersen et al. (2021), the Lyα P_1D was emulated directly on a grid of comoving wavenumbers. Because of the limited box size of our simulations, the P_1D measurements on large scales are affected by cosmic variance. Given that each simulation was run with the same random seed, there are random noise spikes in the power spectrum that align at the same comoving wavenumbers in each simulation used to train the emulator. When it came to testing the pipeline on simulated mock data with a different background evolution, we found that these noise spikes are interpreted by the emulator as sharp features in the power spectrum, artificially enhancing the emulator sensitivity to changes in cosmology. This is due to the fact that the likelihood evaluation is performed in velocity units, which require a conversion from comoving units using H(z). The pipeline would therefore try and find the H(z) that would align the noise spikes in the observed data with the noise spikes in the training set, with two consequences: when running on mock simulations with the same training seed, the pipeline is artificially sensitive to the conversion between comoving and velocity units, and therefore H(z), due to the presence of these sharp features; when running on mock simulations with a different random seed, the pipeline struggles to return the correct cosmology, as the likelihood maximization is dominated by the incentive to align the different noise features.

In order to remove residual noise in the measurements of P_1D, we fit a fourth-order polynomial ¹⁴ to the logarithm of P_1D as a function of the logarithm of wavenumber, using scales k_∥ < 8 Mpc⁻¹:

$$\begin{eqnarray}&&\mathrm{log}{P}_{1{\rm{D}}}({k}_{\parallel })=\displaystyle \sum _{n=0}^{4}{c}_{n}\,{\left(\mathrm{log}{k}_{\parallel }\right)}^{n}.\end{eqnarray} \tag{ 1 }$$

Instead of predicting directly P_1D the emulator now predicts the five coefficients, c_n , of this polynomial, that can later be used to predict P_1D on all scales. The use of fitting functions, such as polynomials or principal component analysis, is a standard practice to reduce noise in emulators, such as in the "Coyote" emulator (Lawrence et al. 2010), as well as in early analyses of P_1D (McDonald et al. 2005).

The variance on the emulated coefficients (${\sigma }_{{c}_{n}}^{2}$) can be used to obtain an estimate for the variance of the emulated P_1D:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\sigma }_{{P}_{1{\rm{D}}}}}{{P}_{1{\rm{D}}}}\right)}^{2}=\displaystyle \sum _{n=0}^{4}{\sigma }_{{c}_{n}}^{2}{\left(\mathrm{log}{k}_{\parallel }\right)}^{2n}.\end{eqnarray} \tag{ 2 }$$

We discuss the impact of cosmic variance on emulator predictions in Appendix B.

In order to test our analysis pipeline, we have generated three synthetic data sets (or mocks) for models that are not included in the training set of the emulator:

• 1.  
  Central simulation. This is the simplest case, a simulation without massive neutrinos or running, with the same background expansion as was used in all training simulations, and a primordial power (A_s , n_s ) corresponding to the center of the Latin hypercube used to set up the training set.
• 2.  
  Neutrino simulation. A simulation with Σm_ν = 0.3 eV, where the cosmological constant (Λ) has been lowered to compensate the increase in the total matter density. The amplitude of the primordial power is also ∼10% larger to compensate the suppression of power caused by massive neutrinos. In Pedersen et al. (2021) we used this simulation to show that we could recover unbiased predictions in cosmologies with massive neutrinos, even when the emulator was trained exclusively with simulations with massless neutrinos.
• 3.  
  Running simulation. A simulation with the same cosmology as the Central simulation, except that its primordial power spectrum has a nonzero running of α_s = 0.015. The other parameters describing the primordial power (A_s , n_s ) have been modified to compensate the change in running and have the same linear power around the pivot scale used in the emulator (k_p = 0.7 Mpc⁻¹; see Table 1).

We start by running a pair of simulations (with inverted phases) for each of the three test models. From each of their 11 snapshots we measure P_1D, in comoving units, and fit a fourth-order polynomial as described in Section 2.2 above. In order to roughly simulate the statistical power of DESI, we use a rescaled version of the Sloan Digital Sky Survey (SDSS) Data Release 14 (DR14) covariance matrix of Chabanier et al. (2019b), where all elements are divided by 5 to approximately take into account the difference in the number of spectra between SDSS DR14 and DESI. ¹⁵ As is common in P_1D measurements, the band powers presented in Chabanier et al. (2019b) are defined in velocity units. At each redshift we compute H(z)/(1 + z) using the simulation cosmology to translate these into wavenumbers in comoving units.

We use a Gaussian likelihood, naturally decomposed into 11 independent sublikelihoods, one for each snapshot (redshift bin). The covariance matrix is the sum of the data covariance and an extra term describing the uncertainty in the emulator predictions, computed with Equation (2). The typical emulator uncertainty is smaller than 1% for models near the center of our training set, and it only has a minor impact on likelihood evaluations around the best-fit values of our analyses. However, it can be larger than 10% when evaluating the likelihood near the convex hull of our training sample.

The different subsections in Sections 3 and 4 use a different number of free cosmological parameters, including the amplitude (A_s ), slope (n_s ), and running (α_s ) of the primordial power spectrum at the usual CMB pivot scale of k_s = 0.05 Mpc⁻¹; the physical densities of baryons (ω_b = Ω_b h²) and of CDM (ω_c = Ω_c h²); the sum of the neutrino masses (Σm_ν ); the Hubble parameter H₀; and the angular acoustic scale of the CMB (θ_MC).

We use four functions to describe the thermal and ionization history of the IGM: the effective optical depth as a function of redshift $\tau (z)=-\mathrm{log}\bar{F}(z)$, the thermal broadening scale (in kilometers per second) at mean densities ${\sigma }_{T}^{\mathrm{vel}}(z)$, the slope of the temperature–density relation γ(z), and the filtering/pressure scale ${k}_{F}^{\mathrm{vel}}(z)$ (in s km⁻¹). Following Pedersen et al. (2021), we measure each of these functions from the Central simulation, and use two parameters, α_X and β_X , to describe a power-law rescaling for each of the four functions. Here the subscript X refers to one of the four IGM parameters. For instance, the thermal broadening scale, ${\sigma }_{T}^{\mathrm{vel}}(z)$, is parameterized as

$$\begin{eqnarray}&&{\left.\mathrm{ln}{\sigma }_{T}^{\mathrm{vel}}(z)=\mathrm{ln}{\sigma }_{T}^{\mathrm{vel}}(z)\right|}_{\mathrm{cen}}+{a}_{{\sigma }_{T}}+{b}_{{\sigma }_{T}}\mathrm{ln}\displaystyle \frac{1+z}{1+3},\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{T}^{\mathrm{vel}}(z){| }_{\mathrm{cen}}$ is the thermal broadening scale in the Central simulation. Therefore, we use a total of eight nuisance parameters related to IGM physics.

There is no guarantee that this simple parameterization is accurate enough to do an analysis on real data, but it should be flexible enough to test the compression of the likelihood in a realistic setting. As described in Table 2, we use combined priors: each parameter is allowed to vary within a given range of values (top-hat prior) and an additional weak Gaussian prior is applied to all parameters; the actual prior is a product of the two.

Table 2. Priors Used for the Cosmological Parameters (Top), and for the Nuisance Parameters Describing the Thermal and Ionization History of the IGM (Bottom)

Parameter	Range Allowed	Gaussian Prior
As (×10−9)	[1.0–3.2]	${ \mathcal N }(2.1,1.1)$
ns	[0.89–1.05]	${ \mathcal N }(0.965,0.08)$
αs	[−0.8–0.8]	${ \mathcal N }(0.0,0.8)$
ωb	[0.018–0.026]	${ \mathcal N }(0.022,0.004)$
ωc	[0.10–0.14]	${ \mathcal N }(0.12,0.02)$
Σmν (eV)	[0.0–1.0]	${ \mathcal N }(0.0,0.5)$
H0	[50–100]	${ \mathcal N }(67.0,25.0)$
θMC(×10−3)	[9.9–10.9]	${ \mathcal N }(10.4,0.5)$
aτ	[−0.1–0.1]	${ \mathcal N }(0.0,0.05)$
bτ	[−0.2–0.2]	${ \mathcal N }(0.0,0.1)$
${a}_{{\sigma }_{T}}$	[−0.4–0.4]	${ \mathcal N }(0.0,0.2)$
${b}_{{\sigma }_{T}}$	[−0.4–0.4]	${ \mathcal N }(0.0,0.2)$
aγ	[−0.2–0.2]	${ \mathcal N }(0.0,0.1)$
bγ	[−0.4–0.4]	${ \mathcal N }(0.0,0.2)$
akF	[−0.2–0.2]	${ \mathcal N }(0.0,0.1)$
bkF	[−0.4–0.4]	${ \mathcal N }(0.0,0.2)$

Note. All parameters have a limited range of values allowed and a Gaussian prior.

Download table as:  ASCIITypeset image

In the next sections we will discuss the constraints on two derived parameters that are able to capture most of the cosmological information in P_1D: the (dimensionless) amplitude and slope of the linear power spectrum at a pivot point k_⋆ = 0.009 s km⁻¹ and redshift z_⋆ = 3 ¹⁶ :

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\star }^{2}=\displaystyle \frac{{k}_{\star }^{3}{P}_{L}({k}_{\star },{z}_{\star })}{2{\pi }^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\left.{n}_{\star }=\displaystyle \frac{{\rm{d}}\ \mathrm{ln}{P}_{L}(k,z)}{{\rm{d}}\ \mathrm{ln}k}\right|}_{{k}_{\star },{z}_{\star }},\end{eqnarray} \tag{ 5 }$$

where P_L (k, z) is the linear power spectrum in velocity units. It is important to highlight that these parameters are defined in velocity units, since P_1D measurements are also presented in velocity units and parameters defined in comoving units would be model dependent.

Let us finish this section by summarizing the steps needed to make a likelihood evaluation:

• 1.  
  Given a set of cosmological parameters, we use the Boltzman solver camb (Lewis et al. 2000) to make predictions for P_L (z, k) and H(z) at all redshifts and scales.
• 2.  
  For each redshift z_i in our mock P_1D measurement, we compute the value of the amplitude (${{\rm{\Delta }}}_{p}^{2}$) and slope (n_p ) of the linear power, P_L (z_i , k), around the pivot point k_p = 0.7 Mpc⁻¹. These are two of the six parameters that will be passed to the emulator to get a prediction of P_1D at z_i .
• 3.  
  The other four parameters ($\bar{F}$, ${\sigma }_{T}^{\mathrm{com}}$, γ, ${k}_{F}^{\mathrm{com}}$) are computed from the eight nuisance parameters and the four IGM-related functions measured from the Central simulation. For instance, we use Equation (3) to compute the thermal broadening scale (${\sigma }_{T}^{\mathrm{vel}}$) in velocity units at redshift z_i , and the comoving scale passed to the emulator is ${\sigma }_{T}^{\mathrm{com}}={\sigma }_{T}^{\mathrm{vel}}(1+{z}_{i})/H({z}_{i})$.
• 4.  
  For each redshift, we ask the emulator to predict the P_1D corresponding to the six emulator parameters computed above. The emulator prediction is in comoving units, and we use H(z_i ) to translate it to velocity units.
• 5.  
  The emulator also returns an uncertainty associated to the prediction, that we add to the data covariance (after translating the emulator covariance to velocity units).
• 6.  
  We use these ingredients to compute a Gaussian likelihood, and multiply it by the prior probability described above.

We use emcee (Foreman-Mackey et al. 2013) to run Monte Carlo Markov Chains, and we use GetDist (Lewis 2019) to make contour plots with marginalized posteriors.

The red contours in Figure 1 show a simplified version of the analysis where only the primordial power parameters (A_s , n_s ) and the eight IGM parameters are varied. In other words, we use a fixed template ¹⁷ for the linear power, P_L (z, k), and rescale it with these two parameters. This analysis is significantly faster than the standard analysis, since we only need to call camb a single time to compute the transfer function for the fiducial cosmology.

The template analysis can be seen as an analysis with infinitely tight priors on the other cosmological parameters (ω_c , H₀, Σm_ν ). While the constraints on the traditional cosmological parameters (A_s , n_s ) are strongly affected by this change in the priors, the constraints on the compressed parameters (${{\rm{\Delta }}}_{\star }^{2}$, n_⋆; top-right panel) remain the same.

In this particular realization of the analysis, we have used a template computed with the same cosmology that was used to run the simulation. Even in the standard analysis (blue contours in Figure 1) we had to assume a value for the baryon density (ω_b = 0.022). We will use the term fiducial cosmology to refer to the cosmological parameters that are being kept fixed in the analysis. Obviously in a real analysis the true cosmology is not known, and so we next test the effect of changing this fiducial comsmology on our results.

In the top panels of Figure 2 we redo the template analysis when using different fiducial cosmologies, with the wrong CDM density (in red) or the wrong sum of the neutrino masses (in blue). While there is a clear bias on the primordial power parameters (left), the compressed parameters are much less affected by the choice of fiducial cosmology.

Zoom In Zoom Out Reset image size

The bottom panels of the same figure show a template analysis for the three test simulations described in Table 1. In all three analyses we use the Central cosmology as our fiducial cosmology. As can be seen in the bottom-left panel, this results in biased posteriors for the primordial power parameters in the Neutrino and Running simulations (stars identify the true values used in each simulation). However, the marginal posteriors of the compressed parameters are again recovered successfully (bottom-right panel). These marginal posteriors in the bottom-right panel will be used in the next section.

Instead of looking at posteriors of f_⋆ and g_⋆ computed as derived parameters in fits for a particular model, we would like to directly sample these without assuming any cosmological model. For instance, in a ΛCDM universe, without curvature or massive neutrinos, both f_⋆ and g_⋆ would be just a function of Ω_m . However, more exotic models could decouple the linear growth from the expansion of the universe, making these parameters independent.

Therefore, in this appendix we directly sample the four compressed parameters (${{\rm{\Delta }}}_{\star }^{2}$, n_⋆, f_⋆, and g_⋆) and the same eight nuisance parameters used in the main text to model the IGM. We use a uniform prior range of [0.24, 0.47], [−2.352, −2.25], [0.9, 1.0], and [0.9, 1.0] for each parameter, respectively. The details of how we do this are detailed later in Section C.2.

In Figure 7 we show constraints on compressed parameters, after marginalizing over the IGM. We use as a mock data set the Neutrino simulation, and show two sets of constraints. In red, we have fixed f_⋆ and g_⋆ to the values of the fiducial cosmology (f_⋆ = 0.981, g_⋆ = 0.968), whereas in blue they are left as free parameters. The dashed lines show the true values in the mock simulation (f_⋆ = 0.989, g_⋆ = 0.969). We note that f_⋆ is very poorly constrained, implying that the P_1D alone is not highly sensitive to the redshift evolution of the linear power spectrum. This result is consistent with the findings of McDonald et al. (2005), although we confirm that this is still the case when using high-precision data sets. The posterior for g_⋆ is slightly better constrained, although it can only rule out very low values of g_⋆ < 0.9. Additionally, there is very little effect on the posteriors for ${{\rm{\Delta }}}_{\star }^{2}$ and n_⋆ when marginalizing over f_⋆ and g_⋆ when compared to fixing them.

Zoom In Zoom Out Reset image size

Note that the red contours were constructed assuming the wrong background cosmology (wrong values of f_⋆ and g_⋆), but that the constraints on ${{\rm{\Delta }}}_{\star }^{2}$ and n_⋆ are nevertheless unbiased.

Here we describe the procedure for mapping from a set of values for the compressed parameters (${{\rm{\Delta }}}_{\star }^{2}$, n_⋆, f_⋆, and g_⋆) to the 11 pairs of emulator parameters (${{\rm{\Delta }}}_{p}^{2}$, n_p ) the values required to generate theoretical predictions for the P_1D from the emulator, one at each redshift. This is done using a fiducial cosmology, as outlined below. We will use k to refer to the (modulus of the) 3D wavenumbers in comoving coordinates, i.e., in megaparsecs. We will use q to refer to the same wavenumber in velocity units. They are related by

$$\begin{eqnarray}&&k=\displaystyle \frac{{\rm{H}}(z)}{1+z}\,q=M(z)\,q.\end{eqnarray} \tag{ C3 }$$

M(z) will play an important role in this discussion.

We will use P(k) to refer to (3D) power spectra in comoving units, i.e., in units of cubic megaparsecs. We will use Q(q) to refer to (3D) power spectra in velocity units, i.e., with units of (km s⁻¹)³. They are related by

$$\begin{eqnarray}&&Q(z,q)={M}^{3}(z)\,P(z,k=M(z)q).\end{eqnarray} \tag{ C4 }$$

In our code we will use a fiducial cosmology as a reference, and parameterize our models as deviations from that cosmology. We will use either subscripts "₀" or superscripts "⁰" to identify functions for the fiducial cosmology. We address changes to the shape of the power spectrum, as described by ${{\rm{\Delta }}}_{\star }^{2}$, n_⋆ (and α_⋆) first, and then later address changes to the redshift evolution using f_⋆ and g_⋆. We can now define the ratio of the linear power between any model and the fiducial one, at the central redshift z_⋆, and in velocity units:

$$\begin{eqnarray}&&B(q)=\displaystyle \frac{{Q}_{\star }(q)}{{Q}_{\star }^{0}(q)}.\end{eqnarray} \tag{ C5 }$$

This will be another important function, tightly related to the linear power parameters that we will end up using.

We fit a second-order polynomial to the logarithm of the linear power spectrum at z_⋆, in velocity units, around a pivot point q_⋆. By default we use z_⋆ = 3 and q_⋆ = 0.009 s km⁻¹, and we fit the polynomial in a range of wavenumbers defined as q_⋆/2 < q < 2q_⋆: ¹⁹

$$\begin{eqnarray}&&{Q}_{\star }(q)\approx A{\left(\displaystyle \frac{q}{{q}_{\star }}\right)}^{{n}_{\star }+{\alpha }_{\star }/2\mathrm{ln}(q/{q}_{\star })},\end{eqnarray} \tag{ C6 }$$

or, equivalently,

$$\begin{eqnarray}&&\mathrm{ln}{Q}_{\star }(q)\approx \mathrm{ln}A+\left[{n}_{\star }+{\alpha }_{\star }/2\mathrm{ln}(q/{q}_{\star })\right]\mathrm{ln}(q/{q}_{\star }).\end{eqnarray} \tag{ C7 }$$

n_⋆ is the first log-derivative around q_⋆, and α_⋆ is the second log-derivative around the same point. Note that the polynomial fit, however, returns ($\mathrm{ln}A$, n_⋆, α_⋆/2). Finally, we define a dimensionless parameter describing the amplitude, ${{\rm{\Delta }}}_{\star }^{2}\,=A\,{q}_{\star }^{3}/(2{\pi }^{2})$. When reconstructing the linear power spectrum using a fiducial cosmology, we use differences in the shape parameters with respect to the fiducial ones:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}B(q) & = & \mathrm{ln}{Q}_{\star }(q)-\mathrm{ln}{Q}_{\star }^{0}(q)\\ & \approx & \left({{\rm{\Delta }}}_{\star }^{2}-{{\rm{\Delta }}}_{\star \,0}^{2}\right)+\left[\left({n}_{\star }-{n}_{\star }^{0}\right)+\displaystyle \frac{{\alpha }_{\star }-{\alpha }_{\star }^{0}}{2}\mathrm{ln}(q/{q}_{\star })\right]\\ & & \times \,\mathrm{ln}(q/{q}_{\star }).\end{array}\end{eqnarray} \tag{ C8 }$$

We are also concerned with reconstructing the linear power spectrum at redshifts other than z_⋆. We ignore neutrinos for now, and work with just the CDM+baryon power spectrum. In this case, we can use the linear growth factor D(z), defined as

$$\begin{eqnarray}&&P(z,k)={\left[\displaystyle \frac{D(z)}{{D}_{\star }}\right]}^{2}{P}_{\star }(k),\end{eqnarray} \tag{ C9 }$$

where, in general, functions y_⋆ = y(z_⋆). We can write the power spectrum at an arbitrary redshift as a function of the fiducial one:

$$\begin{eqnarray*}\begin{array}{rcl}Q(z,q) & = & {M}^{3}(z)\,P(z,k=M(z)q)\\ & = & {M}^{3}(z)\,{\left[\displaystyle \frac{D(z)}{{D}_{\star }}\right]}^{2}{P}_{\star }(k=M(z)q)\\ & = & {\left[\displaystyle \frac{M(z)}{{M}_{\star }}\right]}^{3}{\left[\displaystyle \frac{D(z)}{{D}_{\star }}\right]}^{2}{Q}_{\star }({q}^{{\prime} }=M{(z)/}_{\star }q)\\ & = & {\left[\displaystyle \frac{M(z)}{{M}_{\star }}\right]}^{3}{\left[\displaystyle \frac{D(z)}{{D}_{\star }}\right]}^{2}B({q}^{{\prime} }=M{(z)/}_{\star }q)\\ & & \times \,\,{Q}_{\star }^{0}({q}^{{\prime} }=M{(z)/}_{\star }q)\\ & = & {\left[\displaystyle \frac{M(z)}{{M}_{\star }}\right]}^{3}{\left[\displaystyle \frac{D(z)}{{D}_{\star }}\right]}^{2}B({q}^{{\prime} }=M{(z)/}_{\star }q){\left[{M}_{\star }^{0}\right]}^{3}\\ & & \times \,{P}_{\star }^{0}(k={M}_{\star }^{0}M(z)/{M}_{\star }q)\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl} & = & {\left[\displaystyle \frac{M(z)}{{M}_{\star }}\right]}^{3}{\left[\displaystyle \frac{D(z)}{{D}_{\star }}\right]}^{2}B({q}^{{\prime} }=M{(z)/}_{\star }q){\left[{M}_{\star }^{0}\right]}^{3}\\ & \times & {\left[\displaystyle \frac{{D}_{\star }^{0}}{{D}_{0}(z)}\right]}^{2}{P}^{0}(z,k={M}_{\star }^{0}M(z)/{M}_{\star }q)\\ & = & {\left[\displaystyle \frac{M(z)}{{M}_{\star }}\right]}^{3}{\left[\displaystyle \frac{D(z)}{{D}_{\star }}\right]}^{2}B({q}^{{\prime} }=M{(z)/}_{\star }q){\left[\displaystyle \frac{{M}_{\star }^{0}}{{M}_{0}(z)}\right]}^{3}\\ & \times & {\left[\displaystyle \frac{{D}_{\star }^{0}}{{D}_{0}(z)}\right]}^{2}{Q}^{0}(z,{q}^{{\prime} }=({M}_{\star }^{0}M(z))/({M}_{\star }{M}_{0}(z))q)\\ & = & {\left[m(z)\right]}^{3}{\left[d(z)\right]}^{2}B({q}^{{\prime} }=m(z){M}_{0}{(z)/}_{\star }^{0}q)\,{Q}^{0}(z,{q}^{{\prime} }=m(z)q),\end{array}\end{eqnarray*}$$

where for convenience we have defined two functions:

$$\begin{eqnarray}&&m(z)=\displaystyle \frac{M(z)}{{M}_{\star }}\displaystyle \frac{{M}_{\star }^{0}}{{M}_{0}(z)}\end{eqnarray} \tag{ C10 }$$

and

$$\begin{eqnarray}&&d(z)=\displaystyle \frac{D(z)}{{D}_{\star }}\displaystyle \frac{{D}_{\star }^{0}}{{D}_{0}(z)},\end{eqnarray} \tag{ C11 }$$

which describe differences in expansion rate and in linear growth, respectively.

Using the definition of g_⋆ in Equation (C2), we approximate m(z) using the difference of g_⋆ between the input and the fiducial cosmology as

$$\begin{eqnarray}&&\mathrm{ln}m(z)\approx \displaystyle \frac{3}{2}\left({g}_{\star }-{g}_{\star }^{0}\right)\mathrm{ln}\left(\displaystyle \frac{1+z}{1+{z}_{\star }}\right),\end{eqnarray} \tag{ C12 }$$

or, equivalently,

$$\begin{eqnarray}&&m(z)\approx {\left(\displaystyle \frac{1+z}{1+{z}_{\star }}\right)}^{3/2({g}_{\star }-{g}_{\star }^{0})}.\end{eqnarray} \tag{ C13 }$$

Similarly, we approximate d(z) using the difference of f_⋆ between the input and the fiducial cosmology as

$$\begin{eqnarray}&&\mathrm{ln}d(z)\approx -\left({f}_{\star }-{f}_{\star }^{0}\right)\mathrm{ln}\left(\displaystyle \frac{1+z}{1+{z}_{\star }}\right),\end{eqnarray} \tag{ C14 }$$

or, equivalently,

$$\begin{eqnarray}&&d(z)\approx {\left(\displaystyle \frac{1+z}{1+{z}_{\star }}\right)}^{-({f}_{\star }-{f}_{\star }^{0})}.\end{eqnarray} \tag{ C15 }$$

With these equations, for a given set of (${{\rm{\Delta }}}_{\star }^{2}$, n_⋆, α_⋆, f_⋆, g_⋆), Q(z, q) can be estimated. We then use the approximation of m(z) to convert the velocity unit power spectrum to a comoving power spectrum, and fit a polynomial over the range k_p /2 < k < 2k to obtain values for ${{\rm{\Delta }}}_{p}^{2}$ and n_p . Note that the emulator returns a P_1D in comoving units. The final step is to convert this into velocity units, once again using the above approximation of m(z). This reconstruction process and the composite approximations have been compared against the true values generated in camb , and we verified that they are accurate to within the percent level across all redshifts and extended model spaces considered in this paper.