Fitting functions on the cheap: the relative nonlinear matter power spectrum

TL;DR

This paper introduces a computationally efficient method for constructing accurate nonlinear matter power spectrum fitting functions using the relative spectrum approach, enabling precise predictions with minimal simulations.

Contribution

The authors propose a novel relative spectrum-based fitting method that reduces computational costs and achieves 1% accuracy across diverse cosmologies.

Findings

Achieves 1% accuracy for a wide range of cosmologies up to k≈10/Mpc.

Uses only five inexpensive simulations to construct the fitting function.

Compatible with existing emulators like CosmicEmu for improved predictions.

Abstract

We propose an alternative approach to the construction of fitting functions to the nonlinear matter power spectrum extracted from $N$-body simulations based on the relative matter power spectrum $\delta(k,a)$, defined as the fractional deviation in the absolute matter power spectrum produced by a target cosmology away from a reference $\Lambda$CDM prediction. From the computational perspective, $\delta(k,a)$ is fairly insensitive to the specifics of the simulation settings, and numerical convergence at the 1%-level can be readily achieved without the need for huge computing capacity. Furthermore, $\delta(k,a)$ exhibits several interesting properties that enable a piece-wise construction of the full fitting function, whereby component fitting functions are sought for single-parameter variations and then multiplied together to form the final product. Then, to obtain 1%-accurate absolute power spectrum predictions for any target cosmology only requires that the community as a whole invests in producing one single ultra-precise reference $\Lambda$CDM absolute power spectrum, to be combined with the fitting function to produce the desired result. To illustrate the power of this approach, we have constructed the fitting function RelFit using only five relatively inexpensive $w$CDM simulations (box length $L=256 h^{-1}$Mpc, $N=1024^3$ particles, initialised at $z_i=49$). In a 6-parameter space spanning $\{\omega_m,A_s,n_s,w,\omega_b,h\}$, the output relative power spectra of RelFit are consistent with the predictions of the CosmicEmu emulator to 1% or better for a wide range of cosmologies up to $k\simeq 10$/Mpc. Thus, our approach could provide an inexpensive and democratically accessible route to fulfilling the 1%-level accuracy demands of the upcoming generation of large-scale structure probes, especially in the exploration of "non-standard" or "exotic" cosmologies on nonlinear scales.

Tables

Table 1. Table 1 : Cosmological parameter values of the reference Λ Λ \Lambda CDM model.

Parameter	Symbol	Value
Total physical matter density	${\omega }_m$	0.1422
Physical baryon density	${\omega }_b$	0.0221
Physical neutrino density	${\omega }_{\nu }$	0
Spatial curvature	${\Omega }_k$	0
Effective number of neutrinos	$N_{\mathrm{eff}}$	3.04
Dark energy equation of state parameter	$w$	$-1$
Dimensionless Hubble parameter	$h$	0.673
Primordial scalar fluctuation amplitude at $k_{\mathrm{piv}}=0.05$/Mpc	${10}^9{A}_s$	$2.198$
Scalar spectral index	$n_s$	0.96
Running of the scalar spectral index	$n_{\mathrm{run}}$	0
Tensor-to-scalar ratio	$r$	0
Optical depth to reionisation	$\tau$	0.09

Table 2. Table 2 : Simulations used in section  2 : L 𝐿 L is the simulation box length, N 𝑁 N the number of simulation particles, z i subscript 𝑧 𝑖 z_{i} the initial redshift, r s subscript 𝑟 𝑠 r_{s} the gravitational softening length, and { ω m , A s , n s , w } subscript 𝜔 m subscript 𝐴 𝑠 subscript 𝑛 𝑠 𝑤 \{\omega_{\rm m},A_{s},n_{s},w\} are cosmological model parameters described in table  1 . All except the last two simulations have been performed using Gadget -2/ Camb ; the last two are outputs of Pkdgrav3 / Class .

Run	$L(h^{-1}$Mpc)	$N$	$z_i$	$r_s(h^{-1}$kpc)	${\omega }_m$	${10}^9{A}_s$	$n_s$	$w$
Ref	320	${1282}^3$	49	6	0.1422	2.198	0.96	$-1.00$
Ref$w2$	320	${1282}^3$	49	6	0.1422	2.198	0.96	$-0.85$
Ref2	960	${1282}^3$	49	6	0.1422	2.198	0.96	$-1.00$
Ref2$w2$	960	${1282}^3$	49	6	0.1422	2.198	0.96	$-0.85$
Ref3	1920	${1282}^3$	49	6	0.1422	2.198	0.96	$-1.00$
Ref3$w2$	1920	${1282}^3$	49	6	0.1422	2.198	0.96	$-0.85$
1024Ref-512	512	${1024}^3$	49	6	0.1422	2.198	0.96	$-1.00$
1024$w2$-512	512	${1024}^3$	49	6	0.1422	2.198	0.96	$-0.85$
1024Ref	256	${1024}^3$	49	6	0.1422	2.198	0.96	$-1.00$
1024$w2$	256	${1024}^3$	49	6	0.1422	2.198	0.96	$-0.85$
1024Ref-128	128	${1024}^3$	49	6	0.1422	2.198	0.96	$-1.00$
1024$w2$-128	128	${1024}^3$	49	6	0.1422	2.198	0.96	$-0.85$
768Ref-512	512	${768}^3$	49	6	0.1422	2.198	0.96	$-1.00$
768$w2$-512	512	${768}^3$	49	6	0.1422	2.198	0.96	$-0.85$
768Ref-256	256	${768}^3$	49	6	0.1422	2.198	0.96	$-1.00$
768$w2$-256	256	${768}^3$	49	6	0.1422	2.198	0.96	$-0.85$
768Ref-128	128	${768}^3$	49	6	0.1422	2.198	0.96	$-1.00$
768$w2$-128	128	${768}^3$	49	6	0.1422	2.198	0.96	$-0.85$
512Ref-512	512	${512}^3$	49	6	0.1422	2.198	0.96	$-1.00$
512$w2$-512	512	${512}^3$	49	6	0.1422	2.198	0.96	$-0.85$
512Ref-256	256	${512}^3$	49	6	0.1422	2.198	0.96	$-1.00$
512$w2$-256	256	${512}^3$	49	6	0.1422	2.198	0.96	$-0.85$
512Ref-128	128	${512}^3$	49	6	0.1422	2.198	0.96	$-1.00$
512$w2$-128	128	${512}^3$	49	6	0.1422	2.198	0.96	$-0.85$
1024Ref-256-$z_i$29	256	${1024}^3$	29	6	0.1422	2.198	0.96	$-1.00$
1024$w2$-256-$z_i$29	256	${1024}^3$	29	6	0.1422	2.198	0.96	$-0.85$
1024Ref-256-$r_s$12	256	${1024}^3$	49	12	0.1422	2.198	0.96	$-1.00$
1024$w2$-256-$r_s$12	256	${1024}^3$	49	12	0.1422	2.198	0.96	$-0.85$
pkdgravRef	384	${1024}^3$	49	6	0.1422	2.198	0.96	$-1.00$
pkdgrav$w2$	384	${1024}^3$	49	6	0.1422	2.198	0.96	$-0.85$

Table 3. Table 3 : Simulations discussed in sections  3 , a subset of which—indicated by an asterisk—will be used in section  4 to calibrate our fitting function: L 𝐿 L is the simulation box length, N 𝑁 N the number of simulation particles, z i subscript 𝑧 𝑖 z_{i} the initial redshift, r s subscript 𝑟 𝑠 r_{s} the gravitational softening length, and { ω m , A s , n s , w } subscript 𝜔 m subscript 𝐴 𝑠 subscript 𝑛 𝑠 𝑤 \{\omega_{\rm m},A_{s},n_{s},w\} are cosmological model parameters described in table  1 . All simulations have been performed using Gadget -2/ Camb .

Run	$L(h^{-1}$Mpc)	$N$	$z_i$	$r_s(h^{-1}$kpc)	${\omega }_m$	${10}^9{A}_s$	$n_s$	$w$	Cal
1024Ref	256	${1024}^3$	49	6	0.1422	2.198	0.96	$-1.00$	*
1024$w1$	256	${1024}^3$	49	6	0.1422	2.198	0.96	$-0.92$
1024$w2$	256	${1024}^3$	49	6	0.1422	2.198	0.96	$-0.85$	*
1024$w3$	256	${1024}^3$	49	6	0.1422	2.198	0.96	$-0.75$
1024$w4$	256	${1024}^3$	49	6	0.1422	2.198	0.96	$-1.15$	*
1024${\omega }_{m,l}$	256	${1024}^3$	49	6	0.1381	2.198	0.96	$-1.00$	*
1024${\omega }_{m,h}$	256	${1024}^3$	49	6	0.1361	2.198	0.96	$-1.00$	*
1024${\omega }_{m,l}w1$	256	${1024}^3$	49	6	0.1381	2.198	0.96	$-0.92$
1024${\omega }_{m,h}w1$	256	${1024}^3$	49	6	0.1461	2.198	0.96	$-0.92$
1024${\omega }_{m,l}w2$	256	${1024}^3$	49	6	0.1381	2.198	0.96	$-0.85$
1024${\omega }_{m,h}w2$	256	${1024}^3$	49	6	0.1461	2.198	0.96	$-0.85$
1024$n_{s,l}$	256	${1024}^3$	49	6	0.1422	2.198	0.93	$-1.00$	*
1024$n_{s,h}$	256	${1024}^3$	49	6	0.1422	2.198	0.98	$-1.00$	*
1024$n_{s,l}w2$	256	${1024}^3$	49	6	0.1422	2.198	0.93	$-0.85$
1024$n_{s,h}w2$	256	${1024}^3$	49	6	0.1422	2.198	0.98	$-0.85$
1024$A_{s,l}$	256	${1024}^3$	49	6	0.1422	2.100	0.96	$-1.00$	*
1024$A_{s,h}$	256	${1024}^3$	49	6	0.1422	2.300	0.96	$-1.00$	*
1024$A_{s,l}w2$	256	${1024}^3$	49	6	0.1422	2.100	0.96	$-0.85$
1024$A_{s,h}w2$	256	${1024}^3$	49	6	0.1422	2.300	0.96	$-0.85$
1024${\omega }_{m,l}{A}_{s,l}w2$	256	${1024}^3$	49	6	0.1381	2.100	0.96	$-0.85$
1024${\omega }_{m,h}{A}_{s,h}w2$	256	${1024}^3$	49	6	0.1461	2.300	0.96	$-0.85$

Table 4. Table 4 : Additional simulations of extended cosmological models performed using Gadget -2/ Camb , used in section  4.2 as blind tests (i.e., not calibration) of RelFit : L 𝐿 L is the simulation box length, N 𝑁 N the number of simulation particles, z i subscript 𝑧 𝑖 z_{i} the initial redshift, r s subscript 𝑟 𝑠 r_{s} the gravitational softening length, N eff subscript 𝑁 eff N_{\rm eff} is the effective number of neutrinos, while { w 0 , w a } subscript 𝑤 0 subscript 𝑤 𝑎 \{w_{0},w_{a}\} replace w 𝑤 w to parameterise a possible time dependence of the dark energy equation of state by way of equation ( 4.7 ). All other cosmological parameters not listed here are held at their reference Λ Λ \Lambda CDM values given in table  1 .

Run	$L(h^{-1}$Mpc)	$N$	$z_i$	$r_s(h^{-1}$kpc)	$w_0$	$w_a$	$N_{\mathrm{eff}}$
1024$w2w_{a}1$	256	${1024}^3$	49	6	$-0.85$	$0.3$	3.04
1024$w2w_{a}2$	256	${1024}^3$	49	6	$-0.85$	$0.2$	3.04
1024$w2w_{a}3$	256	${1024}^3$	49	6	$-0.85$	$0.1$	3.04
1024$w2w_{a}4$	256	${1024}^3$	49	6	$-0.85$	$-0.1$	3.04
1024$N_{\mathrm{eff}}3.3$	256	${1024}^3$	49	6	$-1.00$	$0.0$	3.34
1024$N_{\mathrm{eff}}4.0$	256	${1024}^3$	49	6	$-1.00$	$0.0$	4.04

Table 5. Table 5 : Randomly sampled cosmologies on the “surface of 15%-variation” defined by the parameter ranges ( 4.12 ).

Model	${\omega }_m$	$A_s$	$n_s$	$w$	$h$	${\omega }_b$	${\sigma }_8$
M01	0.140092	2.13003	0.914590	$-0.858985$	0.683514	0.022283	0.7725
M02	0.145117	2.19429	0.918589	$-1.13693$	0.615461	0.022034	0.8451
M03	0.140688	2.24912	0.993800	$-0.815920$	0.651725	0.022436	0.8008
M04	0.143777	2.23732	0.915382	$-0.973216$	0.606349	0.021910	0.8206
M05	0.145488	2.35429	0.968273	$-0.849820$	0.615829	0.022298	0.8372
M06	0.145904	2.12986	0.990991	$-1.18144$	0.698259	0.022046	0.8976
M07	0.139331	2.13040	0.911948	$-0.889670$	0.669114	0.022372	0.7734
M08	0.141269	2.00676	0.966808	$-0.879664$	0.638494	0.022906	0.7600
M09	0.148683	2.22116	0.957208	$-0.975085$	0.588572	0.022523	0.8362
M10	0.138605	1.99789	0.941659	$-1.16323$	0.672825	0.021899	0.8125

Equations

$$\delta (\mathbf{\Theta };{\mathbf{\Theta }}_0;k;a)\equiv \frac{P(\mathbf{\Theta };k;a)-P({\mathbf{\Theta }}_0;k;a)}{P({\mathbf{\Theta }}_0;k;a)}$$

$$\mathbf{\Theta }$$

$${\mathbf{\Theta }}_0$$

$$0.1381$$

$$2.1$$

$$0.93$$

$$-1.15$$

$$\gamma (\mathbf{\Theta },{\mathbf{\Theta }}_0;k;a)\equiv \frac{P_{\mathrm{L}}(\mathbf{\Theta };k;a)-P_{\mathrm{L}}({\mathbf{\Theta }}_0;k;a)}{P_{\mathrm{L}}({\mathbf{\Theta }}_0;k;a)},$$

$${\mathbf{\Theta }}_2$$

$${\mathbf{\Theta }}_{1a}$$

$${\mathbf{\Theta }}_{1b}$$

$$1+\delta ({\mathbf{\Theta }}_2,{\mathbf{\Theta }}_0)$$

$$=\left[1+\delta ({\mathbf{\Theta }}_{1b},{\mathbf{\Theta }}_0)\right]\left[1+\delta ({\mathbf{\Theta }}_2,{\mathbf{\Theta }}_{1b})\right],$$

$$\gamma ({\mathbf{\Theta }}_2,{\mathbf{\Theta }}_{1a})$$

$$\gamma ({\mathbf{\Theta }}_2,{\mathbf{\Theta }}_{1b})$$

$$P_{\mathrm{L}}(\mathbf{\Theta }=\{{\theta }_{w,A_s};w,A_s\})=f(w)g(A_s)h({\theta }_{w,A_s}).$$

$$1+\delta ({\mathbf{\Theta }}_2,{\mathbf{\Theta }}_0)\simeq \left[1+\delta ({\mathbf{\Theta }}_{1b},{\mathbf{\Theta }}_0)\right]\left[1+\delta ({\mathbf{\Theta }}_{1a},{\mathbf{\Theta }}_0)\right]$$

$$P_{\mathrm{L}}(\mathbf{\Theta };k;a)=\mathcal{N}(A_s,{\omega }_{\mathrm{m}}){\left(\frac{k}{k_{\mathrm{piv}}}\right)}^{n_s-1}{D}^2(w,{\omega }_{\mathrm{m}},h;a)T^2({\omega }_{\mathrm{m}},{\omega }_{\mathrm{b}};k)$$

$$1+\delta (\mathbf{\Theta },{\mathbf{\Theta }}_0;k;a)$$

$$\simeq \left[1+(\delta /\gamma {)}_{\mathcal{N}}\frac{\Delta \mathcal{N}}{\bar{\mathcal{N}}}\right]\left[1+(\delta /\gamma {)}_{n_s}\frac{\Delta \mathcal{Q}}{\bar{\mathcal{Q}}}\right]\left[1+(\delta /\gamma {)}_D\frac{\Delta D^2}{{\bar{D}}^2}\right]\left[1+(\delta /\gamma {)}_T\frac{\Delta T^2}{{\bar{T}}^2}\right],$$

$$\mathbf{\Theta }$$

$${\mathbf{\Theta }}_0$$

$$1+(\delta /\gamma {)}_{{\omega }_{\mathrm{m}}}{\gamma }_{{\omega }_{\mathrm{m}}}\simeq \left[1+(\delta /\gamma {)}_{\mathcal{N}}\frac{\Delta {\mathcal{N}}_1}{\bar{\mathcal{N}}}\right]\left[1+(\delta /\gamma {)}_D\frac{\Delta D_1^2}{{\bar{D}}^2}\right]\left[1+(\delta /\gamma {)}_T\frac{\Delta T^2}{{\bar{T}}^2}\right],$$

$$\Delta {\mathcal{N}}_1$$

$$\Delta D_1^2$$

$$1+\delta (\mathbf{\Theta },{\mathbf{\Theta }}_0;k;a)\simeq \frac{1+(\delta /\gamma {)}_{\mathcal{N}}\frac{\Delta \mathcal{N}}{\bar{\mathcal{N}}}}{1+(\delta /\gamma {)}_{\mathcal{N}}\frac{\Delta {\mathcal{N}}_1}{\bar{\mathcal{N}}}}\left[1+(\delta /\gamma {)}_{n_s}\frac{\Delta \mathcal{Q}}{\bar{\mathcal{Q}}}\right]\left[1+(\delta /\gamma {)}_{{\omega }_{\mathrm{m}}}{\gamma }_{{\omega }_{\mathrm{m}}}\right]\frac{1+(\delta /\gamma {)}_D\frac{\Delta D^2}{{\bar{D}}^2}}{1+(\delta /\gamma {)}_D\frac{\Delta D_1^2}{{\bar{D}}^2}},$$

$$g^{\prime \prime }+\left[\frac{7}{2}-\frac{3}{2}\frac{w(a)}{1+X(a)}\right]\frac{g^{\prime }}{a}+\frac{3}{2}\frac{1-w(a)}{1+X(a)}\frac{g}{a^2}=0$$

$$X(a)$$

$$=\frac{{\omega }_{\mathrm{m}}}{h^2-{\omega }_{\mathrm{m}}}{a}^{3(w_0+w_a)}\exp \left[3w_a(1-a)\right],$$

$$g(a)=\exp {\int }_0^{a}d\ln a^{\prime }\left[\Omega (a^{\prime }{)}^{\rho }-1\right],$$

$${\gamma }_{n_s}$$

$$(\delta /\gamma {)}_{n_s}=\frac{(\delta /\Gamma {)}_{n_s}}{\ln (k/k_{\mathrm{piv}})},$$

$$1+\delta (\mathbf{\Theta },{\mathbf{\Theta }}_0;k;a)$$

$$\simeq \left[1+(\delta /\gamma {)}_{\mathcal{N}}\frac{\Delta \mathcal{N}}{\bar{\mathcal{N}}}\right]\left[1+(\delta /\Gamma {)}_{n_s}\Gamma \right]\left[1+(\delta /\gamma {)}_D\frac{\Delta D^2}{{\bar{D}}^2}\right]\left[1+(\delta /\gamma {)}_T\frac{\Delta T^2}{{\bar{T}}^2}\right],$$

$$1+\delta (\mathbf{\Theta },{\mathbf{\Theta }}_0;k;a)\simeq \frac{1+(\delta /\gamma {)}_{\mathcal{N}}\frac{\Delta \mathcal{N}}{\bar{\mathcal{N}}}}{1+(\delta /\gamma {)}_{\mathcal{N}}\frac{\Delta {\mathcal{N}}_1}{\bar{\mathcal{N}}}}\left[1+(\delta /\Gamma {)}_{n_s}\Gamma \right]\left[1+(\delta /\gamma {)}_{{\omega }_{\mathrm{m}}}{\gamma }_{{\omega }_{\mathrm{m}}}\right]\frac{1+(\delta /\gamma {)}_D\frac{\Delta D^2}{{\bar{D}}^2}}{1+(\delta /\gamma {)}_D\frac{\Delta D_1^2}{{\bar{D}}^2}}.$$

$$\delta (\mathbf{\Theta },{\mathbf{\Theta }}_0)=\sum _X{\frac{\partial \ln P}{\partial \ln X}|}_{\mathbf{\Theta }={\mathbf{\Theta }}_0}\frac{\Delta X}{\bar{X}},$$

$$y(k,a)\equiv \frac{k}{k_{\mathrm{peak}}(a)},$$

Full text

1 Introduction

The upcoming generation of large-scale structure surveys such as the ESA Euclid mission [1] and the Large Synoptic Survey Telescope (LSST) [2] have the potential to measure cosmological observables at an unprecedented level of precision. In terms of the matter power spectrum, the measurement uncertainty is expected to be at the 1% level down to length scales corresponding to wavenumbers $k\sim O(5)\,h$/Mpc. Such high precisions in turn put heavy demands on theoretical calculations of the observables.

On large scales where perturbations are expected to remain well below $O(1)$, linear perturbation theory can easily satisfy the 1% precision requirement. Likewise, perturbative methods can be extended to higher orders on weakly nonlinear scales ($k\sim 0.05\to 0.1\,h$/Mpc at scale factor $a=1$), and much effort has been devoted recently towards improving the convergence of these computations (see, e.g., [3]). Calculations in the fully nonlinear scales, i.e,. $k\gtrsim O(0.1)\,h$/Mpc at $a=1$, however, belong in the domain of numerical simulations.

However, simulations are inherently computationally expensive, and it is currently not economical to run full simulations for more than a select $O(10\to 100)$ parameter combinations “representative” of a large cosmological parameter space. In fact, achieving the required 1% precision for even one single set of cosmological parameters is a computational challenge that necessitates the use of some of the largest computing facilities in the world [4, 5]. As an example, each cosmology in the Mira–Titan suite of $w$CDM simulations is realised by two high-resolution simulations with 30 billion+ and 60 billion+ particles each, plus 16 lower-resolution 100-million-particle runs [6, 7]. Only a select few researchers in the world have access to the requisite computing power to carry out such calculations en masse.

Currently, in order to explore large parameter spaces with parameter combinations running into $O(10^{5})$—as is required in a typical Markov Chain Monte Carlo parameter estimation analysis—the favoured approach is to employ fitting functions such as Halofit [8, 9] or HMCode [10, 11] that have been calibrated against simulation results. Alternatively, one can interpolate between a set of simulations spanning the parameter spaces of interest, such as the emulator approach of [12, 4, 13, 14]. However, as the accuracy of any fitting or interpolation function is contingent upon there being sufficient calibrators to fairly sample the parameter space and the calibrating simulations themselves having the required level of precision, the burden is again back on the simulations and the same select few research groups that have the computing resources to supply these calculations. Such a strong reliance on computing resources clearly poses severe limitations on the participation of the wider scientific community, especially in the exploration of “exotic” cosmologies such as decaying or interacting dark matter (e.g., [15, 16]), or dark energy perturbations (e.g., [17]) on nonlinear scales.

In this paper we put forward a different approach to constructing fitting functions to the nonlinear matter power spectrum that will alleviate to a large extent the precision burden on the calibrating simulations and potentially democratise the exploration of precision cosmology on nonlinear scales: Instead of the usual practice of fitting or interpolating directly the absolute simulated matter power spectrum $P({\bf\Theta};k;a)$ for a select few cosmological parameter combinations ${\bf\Theta}$, we propose to construct a fitting function to a set of spectra $\delta({\bf\Theta};{\bf\Theta}_{0};k;a)$, defined as

[TABLE]

relative to the absolute matter power spectrum of a reference cosmological model, $P({\bf\Theta}_{0};k;a)$. As we shall demonstrate, there are a number of reasons why fitting the relative power spectra may be superior to fitting their absolute counterparts:

From the computational perspective, relative power spectra can be calculated much more precisely than absolute power spectra from $N$-body simulations using the same box size and number of particles [18]. This is because many systematic uncertainties are multiplicative and affect all simulations in the same way; taking the ratio of two simulation results therefore enable these uncertainties to cancel to a large extent. Indeed, the use of ratios to “get around” systematic uncertainties that may not be completely well understood is a well-known technique used in many areas of physics, e.g., collider phenomenology [19], precision cosmology [20], and neutrino physics [21, 22].

An immediate corollary of this observation is that a nominal accuracy goal can be achieved at much a lower computational cost using relative power spectrum simulations than their absolute counterparts. Once a fitting function to $\delta=\delta({\bf\Theta};{\bf\Theta}_{0};k;a)$ is available as a function of the underlying cosmology, to obtain an accurate estimation of a target $P({\bf\Theta};k;a)$ for any parameter combination requires only that we perform one single ultra-high precision simulation of the reference cosmological model to establish $P({\bf\Theta}_{0};k;a)$ and then combine this result with the fitting function. In this regard, our proposal parallels the “halo model reaction” approach of [23, 24] and the CosmicEFT approach of [25], wherein the equivalent of $\delta$ is computed using semi-analytical methods such as the halo model and effective field theory. 2. 2.

The present generation of linear cosmological probes, e.g., measurements of the cosmic microwave background (CMB) temperature and polarisation anisotropies by the Planck mission [26, 27], already constrains cosmology to the extent that variations in the absolute power spectra are typically $\lesssim 10\%$. This means that any fitting function to $\delta$ need only be calibrated to at most $\sim 10$%-precision in order to reproduce a target $P({\bf\Theta};k;a)$ with $\lesssim 1$%-level error (assuming, of course, that an ultra-precise reference $P({\bf\Theta}_{0};k;a)$ is available), and the smaller $\delta$ is the laxer the calibration precision requirement. This is a trivial demand in comparison with the 1% calibration precision required of fitting functions designed to directly reproduce $P({\bf\Theta};k;a)$. 3. 3.

Since typically $|\delta|\lesssim 0.1$, it is strongly suggestive that the relative matter power spectrum $\delta$ may be computable perturbatively from similarly small deviations in the linear power spectrum away from the reference cosmology. Indeed, we find that $\delta$ can be related to relative changes in, e.g., the linear growth function, the primordial power spectrum, etc., in a remarkably cosmology-independent way. This attractive feature enables a multiplicative construction of the full fitting function, whereby component fitting functions are sought for variations of cosmological model parameters (or their proxies such as the linear growth function) one at a time and the full fitting function pasted together via a simple multiplication of the components.

The paper is organised as follows. We begin in section 2 with a discussion of the convergence of the absolute and the relative matter power spectrum, using cosmologies with a non-canonical dark energy equation of state parameter as an example. Section 3 examines the properties of the relative power spectrum under single- and multi-parameter variations, through which we motivate a strategy for the construction of a fitting function for $\delta$. We propose specific functional forms for the fitting function in section 4 , which we then calibrate against $N$-body simulations to produce RelFit. Comparisons of the predictions of RelFit and other approaches are presented in the same section. Section 5 contains our conclusions. Throughout the work we use as our reference cosmology ${\bf\Theta}_{0}$ a $\Lambda$CDM model with parameter values given in table 1, roughly comparable to the best-fit of the Planck 2015 CMB data [26]. Where confusion is unlikely to arise, we shall sometimes omit writing out the dependences of the absolute and relative matter power spectra on $k$ and/or $a$.

2 Numerical convergence of the absolute and the relative spectrum

Many factors may influence the numerical convergence of a simulation result. Chief amongst these are the simulation box size and the number of particles employed to sample the cosmological fluid (i.e., cold dark matter in a $\Lambda$CDM-type cosmology) phase space. Other important factors include the redshift at which a simulation is initialised, and the gravitational softening length adopted in the simulation to prevent spurious relaxation. In this section we examine the extent to which numerical convergence of the absolute and relative power spectra depends on these factors, using a series of $N$-body simulations performed with the Gadget-2 code [28]. The specifics of each simulation are summarised in table 2.

For each simulation we employ initial conditions generated via the Zel’dovich approximation from linear transfer functions outputted by Camb [29]. We include baryons in the computation of the linear transfer function required for initial condition generation, but do not distinguish baryons from cold dark matter in the actual simulations. The latter is certainly an oversimplification in precision calculations of an absolute power spectrum, but can be expected to be a reasonable approximation in the case of a relative power spectrum.

2.1 Box size and number of particles

It is well known that numerical convergence of the absolute matter power spectrum requires simulations in large boxes with many particles. If the box size is too small sample (cosmic) variance becomes a serious issue. Increasing the box size however requires that we also up the number of particles in order to suppress shot noise on small scales. These issues have been discussed in detail in, e.g., a series of papers related to the Coyote simulations (e.g., [4]) and more recently in [5]. The general conclusion is that to achieve an absolute power spectrum calculation at the 1% level of precision requires box lengths exceeding $L=500\,h^{-1}$Mpc and particle numbers of order $N=10000^{3}$.

To illustrate the lack of convergence of the absolute matter power spectrum $P({\bf\Theta};k;a)$, we show in figure 1 $P({\bf\Theta}_{0};k;a)$ constructed from various reference $\Lambda$CDM simulations using different box sizes and particle numbers (but keeping for now the initialisation redshift and softening length fixed at $z_{i}=49$ and $r_{s}=6\,h^{-1}$kpc respectively) summarised in table 2. These are normalised to a benchmark power spectrum $P_{\rm bench}({\bf\Theta}_{0};k;a)$, constructed from amalgamating the power spectra extracted from three “high-quality” runs—Ref3 at $k<1\,h$/Mpc, Ref2 in the range $1\,h$/Mpc$<k<3\,h$/Mpc, and Ref at $k>3\,h$/Mpc. Clearly, no single simulation is able to converge to the benchmark at better than 10% across the entire $k$-range, and convergence worsens as the scale factor $a$ approaches unity.

In contrast, the relative change in the matter power spectrum between two cosmologies with different parameter values, $\delta=\delta({\bf\Theta},{\bf\Theta}_{0};k;a)$ as defined in equation (1.1), is much less susceptible to sample variance, provided the two simulations used to construct $\delta$ have been run under identical conditions and initialised with identical phases in the density field. We emphasise that these requirements of identical simulations settings are crucial, as it is precisely this sameness that ensures two simulations suffer largely the same systematic effects that eventually cancel out when forming a ratio, leaving a $\delta$ that is ultimately relatively insensitive to the simulation settings. A similar observation has also been made in [18].

This relative insensitivity to the simulation settings also means that numerical convergence in $\delta$ can be achieved using much smaller boxes and hence smaller numbers of simulation particles than in the case of the absolute power spectrum. Figure 2 illustrates this point by way of the relative change in power $\delta({\bf\Theta},{\bf\Theta}_{0};k;a)$ between two cosmological models specified respectively by the parameter values

[TABLE]

where $w$ denotes the dark energy equation of state parameter, and $\theta_{w}=\bar{\theta}_{w}$ stipulates that all other model parameters besides $w$ are to be held fixed at their reference values $\bar{\theta}_{w}$ given in table 1. As in figure 1, the relative power spectra here have been constructed from the simulations of table 2 using different combinations of box sizes and particle numbers, and for clarity we have subtracted away the benchmark relative power spectrum $\delta_{\rm bench}$ constructed from the “high-quality” Ref$w$2, Ref2$w$2, Ref3$w$2, Ref, Ref2, and Ref3 runs of table 2.

Clearly, independently of the number of simulation particles employed, sample variance dominates when the box size is too small, but becomes manageable once the box length reaches $L=256\,h^{-1}$Mpc. In terms of particle numbers, we find $N=1024^{3}$ to be sufficient to eliminate to a large extent shot noise in boxes of side length $L\geq 256\,h^{-1}$Mpc, enabling numerical convergence at the $0.01$ level down to wavenumbers close to the Nyquist frequency at $a=0.85$ and better than $0.005$ at $a=0.50$; even at $a=1$, convergence at the (not unacceptable) 0.02 level is possible for a large range of wavenumbers. Importantly, these conclusions are independent of the choice of initial phases, as demonstrated in figure 3, where we have re-simulated the relative power spectrum of the two cosmologies of equation (2.1) using four different sets of initials seeds for the setting $L=256\,h^{-1}$Mpc and $N=1024^{3}$, and plotted their deviations from the ensemble average $\langle\delta\rangle$.

Note that the alternative choice of $L=512\,h^{-1}$Mpc and $N=1024^{3}$ could even enable the attainment of 0.01 numerical convergence at $a=1$, as shown in figure 2. The downside, however, is that such a setting yields power spectrum predictions only up to $k=5\,h$/Mpc, and to achieve a better resolution in $L=512\,h^{-1}$Mpc boxes would require a computing capacity beyond our current means. Henceforth, we shall adopt the setting $L=256\,h^{-1}$Mpc and $N=1024^{3}$, a fair compromise between computing power and the accuracy demands of future large-scale structure probes,111A scale factor $a=0.85$ corresponds to a redshift $z=0.176$, reasonably low relative to the median redshift $z_{\rm m}=0.8\to 0.9$ of the Euclid and LSST galaxy redshift surveys [1, 2]. Similarly, while cosmic shear is in principle sensitive to the matter distribution at $z=0$, in practice the lensing weights are dominated by structures at roughly half the source-to-observer comoving distance; for a shear tomographic bin at $z=0.5\to 1.0$, such as used in the Euclid parameter sensitivity forecast [1], the weight peaks at $z\sim 0.3$. This encourages us to think that 0.01 numerical convergence of the matter power spectrum down to $z=0.176$ may suffice. and restrict our attention to $a\leq 0.85$.

2.2 Initial redshift and gravitational softening

We consider also the sensitivity of the absolute and relative matter power spectrum to the simulation initial redshift $z_{i}$ and gravitational softening length $r_{s}$, and vary these simulation parameters from the default $z_{i}=49$ and $r_{s}=6\,h^{-1}$kpc to $z_{i}=29$ and $r_{s}=12\,h^{-1}$kpc respectively. The results at $a=0.85,0.50$ are shown in figure 4.

Evidently, changing the gravitational softening length has no discernible effect on the relative power spectrum at either $a=0.85$ or $a=0.50$, and alters the absolute power spectrum only at the percent level at $k=10\,h$/Mpc. On the other hand, with initial conditions set by the Zel’dovich approximation, both initialisation redshifts tested are clearly too low to achieve reasonable accuracy for the absolute power spectrum because of long-lived transients (although the problem of transients can be avoided by adopting 2LPT initial conditions [30]). The relative power spectrum, however, appears to be largely insensitive to $z_{i}$ within the 0.01 accuracy requirement.

Of course the case of varying only $w$ away from its reference $\Lambda$CDM value is particularly benevolent in the sense that even for $w=-0.85$ the evolution history of the density perturbations at $z\gg 1$ is essentially identical to the $w=-1$ case. This means that any transient excited as a result of the initialisation procedure must be identical in both cases, and cancel out exactly when we form the relative power spectrum.

2.3 Code comparison

Lastly, we compare the stability of the relative matter power spectrum predictions with respect to the numerical codes used to generate the results. This we achieve by performing a new set of simulations using the Pkdgrav3 $N$-body code [31] with initial conditions generated from the linear transfer function output of Class [32] using the method outlined in [33]. (Recall that our main suite of simulation results have been generated with Gadget-2 and Camb.) The simulation specifics are summarised in the last two lines of table 2.

Figure 5 shows the relative matter power spectrum between the $w=-0.85$ and the $w=-1$ cosmologies at $a=0.85,0.50$ computed in this manner relative to the Gadget-2/Camb benchmark $\delta_{\rm bench}$. Clearly, the discrepancy between the two different code outputs is well within our 0.01 accuracy tolerance and comparable to the deviations one would expect from initialising the simulations with different sets of random phases (see figure 3). We therefore conclude that the relative power spectrum is largely insensitive to $N$-body code from which it is generated or to the linear Boltzmann solver that provides the initial conditions.

3 Properties of the relative power spectrum

Having established the advantage of the relative matter power spectrum $\delta$ over its absolute counterpart in terms of numerical convergence, we now examine its properties more closely, in order to devise a fitting strategy and eventually a functional form that can directly fit $\delta$.

To this end we have performed a suite of simulations summarised in table 3 using Gadget-2/Camb, varying the parameter values of $\omega_{\rm m}$, $A_{s}$, $n_{s}$, and $w$ away from their reference $\Lambda$CDM values one at a time as well as in combination in the ranges

[TABLE]

In terms of measurement uncertainties, our choice of $\omega_{\rm m}$ values spans a range comparable to $9.2$ times the standard deviation inferred from the 2018 Planck+external data combination222For $\omega_{\rm m}$, $A_{s}$, and $n_{s}$, this means the 2018 Planck TT+TE+EE+lowE+lensing+BAO combination [27], while for $w$ the set includes also SNe. in a vanilla 6-parameter $\Lambda$CDM fit ($\sigma(\omega_{\rm m})=0.00087$) [27]; $6.6$ times for the primordial fluctuation amplitude ($\sigma(10^{9}A_{s})=0.030$); $13$ for the spectral index ($\sigma(n_{s})=0.0038$); and $12$ for a time-independent dark energy equation of state parameter ($\sigma(w)=0.032$). While these ranges differ between parameters in terms of the number of standard deviations, the effects the parameter variations produce on the nonlinear matter power spectrum are of very similar magnitudes—typically no more than $20$% at $a=0.85$, as shown in figures 6 to 8.

Two interesting properties of $\delta({\bf\Theta};{\bf\Theta}_{0};k;a)$ can be discerned from our simulation set: an approximate universality and multiplicability. We discuss these properties in detail below, and propose how they can be jointly exploited as a strategy for constructing a fitting function to any general $\delta({\bf\Theta},{\bf\Theta}_{0};k;a)$ in a multivariate parameter space.

3.1 Approximate universality: varying one parameter at a time

Consider figure 6. The solid lines in the top panels show the relative matter power spectra $\delta=\delta({\bf\Theta},{\bf\Theta}_{0};k;a)$ at $a=0.85,0.50,0.30$ of four target cosmological models described by ${\bf\Theta}=\{\theta_{w}=\bar{\theta}_{w};w=-0.73,-0.85,-0.92,-1.15\}$, where all non-$w$ model parameters $\theta_{w}$ are held fixed at their reference values $\bar{\theta}_{w}$, relative to the canonical reference $\Lambda$CDM model ${\bf\Theta}_{0}=\{\theta_{w}=\bar{\theta}_{w};w=-1\}$ of table 1. The middle panels show the same cosmological models in a similar construct $\gamma=\gamma({\bf\Theta},{\bf\Theta}_{0};k;a)$, defined as

[TABLE]

i.e., akin to $\delta({\bf\Theta},{\bf\Theta}_{0};k;a)$, but with the target and reference absolute power spectra $P({\bf\Theta})$ and $P({\bf\Theta}_{0})$ replaced with their linear counterparts $P_{\rm L}({\bf\Theta})$ and $P_{\rm L}({\bf\Theta}_{0})$ outputted by Camb [29]. The bottom panels show the ratios $\delta/\gamma$.

An immediately notable feature in figure 6 is that despite their differences in $\delta$ and $\gamma$, at each scale factor $a$ and over a wide range of wavenumbers $k$, all four target cosmologies return a functional form for the ratio $\delta/\gamma$ that is quantitatively remarkably independent of the chosen value of $w$; the function tends to unity in the linear regime, peaks at an $a$-dependent $k_{\rm peak}$, and drops off to zero at large $k$ values. At $k\lesssim 4\ h$/Mpc the agreement between models is always better than 10%. This “approximate universality” of $\delta/\gamma$ likewise holds for a reference $w$ value different from the canonical choice of $-1$ in ${\bf\Theta}_{0}$, as demonstrated by the dashed lines in figure 6 (which feature $w=-0.92$ in ${\bf\Theta}_{0}$), provided of course that we choose the same reference $w$ for both $P({\bf\Theta}_{0})$ and $P_{\rm L}({\bf\Theta}_{0})$ in the construction of $\delta/\gamma$.

Approximate universality in $\delta/\gamma$ extends also to the case in which we employ a set of the non-$w$ cosmological parameters $\theta_{w}$ different from $\bar{\theta}_{w}$, again on the understanding that whatever values we choose for $\theta_{w}$ in the construction of $\delta/\gamma$ are held constant across the four target and reference absolute power spectra, $P({\bf\Theta}),P({\bf\Theta}_{0}),P_{\rm L}({\bf\Theta})$, and $P_{\rm L}({\bf\Theta}_{0})$. This is illustrated in figure 7 by the solid lines, representing $\delta,\gamma$, and $\delta/\gamma$ constructed from a selection of target and reference cosmologies from the simulations of table 3, where ${\bf\Theta}=\{\theta_{w}=\bar{\theta}_{w};w=-0.85\}$, ${\bf\Theta}_{0}=\{\theta_{w}=\bar{\theta}_{w};w=-1\}$, and $\theta_{w}\neq\bar{\theta}_{w}$. In the same figure, the cosmological models represented by the dashed lines feature in addition a non-canonical reference $w$ value in ${\bf\Theta}_{0}$ (in this instance, $w=-0.92$); again, their respective $\delta/\gamma$ conforms to the same approximately universal form already observed amongst the solid lines as well as in figure 6.

So far we have discussed the approximate universality of $\delta/\gamma$ exclusively in the context wherein the target and reference cosmologies, ${\bf\Theta}$ and ${\bf\Theta}_{0}$, differ only by their $w$ parameter value. To further test the hypothesis of a $\delta/\gamma$ approximate universality under variation of any one cosmological parameter besides $w$, we show in figure 8 $\delta,\gamma$, and $\delta/\gamma$ for three families of relative power spectra at $a=0.85$ described by

1. 1a.

Solid: ${\bf\Theta}=\{\theta_{A_{s}};10^{9}A_{s}=2.100,2.300\}$, ${\bf\Theta}_{0}=\{\theta_{A_{s}};10^{9}A_{s}=2.198\}$; 2. b.

Dashed: ${\bf\Theta}=\{\theta_{A_{s}};10^{9}A_{s}=2.300\}$, ${\bf\Theta}_{0}=\{\theta_{A_{s}};10^{9}A_{s}=2.100\}$; 3. 2a.

Solid: ${\bf\Theta}=\{\theta_{\omega_{\rm m}};\omega_{\rm m}=0.1381,0.1461\}$, ${\bf\Theta}_{0}=\{\theta_{\omega_{\rm m}};\omega_{\rm m}=0.1422\}$; 4. b.

Dashed: ${\bf\Theta}=\{\theta_{\omega_{\rm m}};\omega_{\rm m}=0.1461\}$, ${\bf\Theta}_{0}=\{\theta_{\omega_{\rm m}};\omega_{\rm m}=0.1381\}$; 5. 3a.

Solid: ${\bf\Theta}=\{\theta_{n_{s}};n_{s}=0.93,0.98\}$, ${\bf\Theta}_{0}=\{\theta_{n_{s}};n_{s}=0.96\}$; 6. b.

Dashed: ${\bf\Theta}=\{\theta_{n_{s}};n_{s}=0.98\}$, ${\bf\Theta}_{0}=\{\theta_{n_{s}};n_{s}=0.93\}$.

See also figure 9 for $\delta/\gamma$ of these models at $a=0.50,0.30$. Here, the convention $\theta_{X}$ again denotes all model parameters other than $X$, and we consider both $\theta_{X}=\bar{\theta}_{X}$ and $\theta_{X}\neq\bar{\theta}_{X}$ selected from the simulations of table 3. Again, the close similarity of $\delta/\gamma$ within each family is unmistakable. In the case of variations in $A_{s}$ and $\omega_{\rm m}$, we see that the $a$-dependent locations of the peaks are similar to $k_{\rm peak}$ previously identified for variations in $w$.

Note that in the case of variation of $\omega_{\rm m}$, $\delta/\gamma$ exhibits prominent oscillations at $k\lesssim 1\ h$/Mpc. Oscillations arise in the first place from a phase difference in the baryon acoustic oscillations between cosmologies with different matter densities, and can already be seen in both $\delta$ and $\gamma$. Nonlinear evolution additionally alters the amplitudes and phases of these oscillations, so that a residual survives in $\delta/\gamma$.

A final remark concerns the singularities in $\delta/\gamma$ under variation of $n_{s}$ observed in figures 8 and 9. These are artefacts following from our choice of pivot scale $k_{\rm piv}=0.05$/Mpc for the primordial power spectrum ${\cal P_{R}}(k)=A_{s}(k/k_{\rm piv})^{n_{s}-1}$. In fact, a singularity will arise in $\delta/\gamma$ whenever the linear power spectra of the target and reference cosmologies cross over. A judicious choice of $k_{\rm piv}$, e.g., $k_{\rm piv}=0.002$/Mpc, would have confined such cross-overs to scales outside of the range of interest and facilitated the task of finding a fitting function. However, as we shall discuss in section 3.3, rather than re-running simulations with a different $k_{\rm piv}$, it transpires that for power-law primordial power spectra the remedy is very simple.

Then, to summarise section 3.1, for a family of relative matter power spectra described by the target and reference cosmological model parameters ${\bf\Theta}=\{\theta_{X};X\}$ and ${\bf\Theta}_{0}=\{\theta_{X};\bar{X}\}$, the ratio of the relative (nonlinear) power spectrum to the relative linear power spectrum, $\delta/\gamma$, is, at each scale factor $a$ and over a wide range of wavenumbers $k$, largely independent of the values of $\theta_{X}$, $X$, and $\bar{X}$. We shall denote this approximately universal ratio $(\delta/\gamma)_{X}$.

3.2 Multiplicability: varying two or more parameters at a time

Consider now three target cosmological models specified respectively by the parameters

[TABLE]

where $\theta_{X,Y}$ denotes all model parameters besides $X$ and $Y$, $\bar{\theta}_{X,Y}$ their reference values in table 1, and our canonical reference model is again defined by ${\bf\Theta}_{0}=\{\theta_{w}=\bar{\theta}_{w};w=-1\}$.

From the definition (1.1) it is easy to establish that the three target cosmologies must have relative power spectra $\delta$ satisfying at all times the general relations

[TABLE]

irrespective of our exact choice of model parameter values. For the particular target cosmologies (3.3) under consideration, the corresponding relative linear power spectra $\gamma$ also happen to obey

[TABLE]

because of the especially simple and, importantly, separable effects variations of $w$ and $A_{s}$ induce on the absolute linear power spectrum, in the sense that $P_{\rm L}({\bf\Theta})$ is a separable function of $w$, $A_{s}$, and $\theta_{w,A_{s}}$:

[TABLE]

It then follows straightforwardly from the approximate universality of $\delta/\gamma$ discussed in section 3.1 that $\delta({\bf\Theta}_{2},{\bf\Theta}_{1a})\simeq\delta({\bf\Theta}_{1b},{\bf\Theta}_{0})$ and $\delta({\bf\Theta}_{2},{\bf\Theta}_{1b})\simeq\delta({\bf\Theta}_{1a},{\bf\Theta}_{0})$, and hence

[TABLE]

as an approximation to equation (3.4).

The top panels of figure 10 demonstrate the remarkable correspondence between the exact $\delta({\bf\Theta}_{2},{\bf\Theta}_{0})$ and its approximation constructed from $\delta({\bf\Theta}_{1b},{\bf\Theta}_{0})$ and $\delta({\bf\Theta}_{1a},{\bf\Theta}_{0})$ via equation (3.7) for $w=-0.85$ and $10^{9}A_{s}=2.100,2.300$; at all scale factors and for the entire range of wavenumbers under consideration, the approximation is able to reproduce the exact relative matter power spectrum to $0.01$ or better. The bottom panels provide a second example of this excellent correspondence for the target cosmologies ${\bf\Theta}_{2}=\{\theta_{w,n_{s}}=\bar{\theta}_{w,n_{s}};w=-0.85,n_{s}=0.93,0.98\}$, ${\bf\Theta}_{1a}=\{\theta_{w}=\bar{\theta}_{w};w=-085\}$, and ${\bf\Theta}_{1b}=\{\theta_{n_{s}}=\bar{\theta}_{n_{s}};n_{s}=0.93,0.98\}$ (for which the equivalents of equations (3.5), (3.6), and hence (3.7) also hold).

Naturally, alternatively to equation (3.7), the approximate universality of $\delta/\gamma$ under variation of one parameter means that we could also have approximated 1+$\delta({\bf\Theta}_{2},{\bf\Theta}_{0})$ using instead $\left[1+\delta({\bf\Theta}_{2},{\bf\Theta}_{1b})\right]\left[1+\delta({\bf\Theta}_{2},{\bf\Theta}_{1a})\right]$—or, indeed, any other combination of two relative power spectra in which we vary only one parameter at a time and whose linear counterparts equate to the relations (3.5)—with similarly good although not identical results to figure 10. The essence of equation (3.7), however, lies in its suggestion that the multiplicative nature of the relative power spectrum and the approximate-universal form $(\delta/\gamma)_{X}$ under variation of $X$ may be jointly exploited as a relatively simple strategy for constructing a fitting function to any general $\delta({\bf\Theta},{\bf\Theta}_{0})$ in a multivariate parameter space.

Furthermore, the condition of separability (3.6) implies that the natural division of cosmological models into families (for the purpose of finding the approximate-universal forms $(\delta/\gamma)_{X}$) is not in terms of the model parameters per se, but rather their linear “proxies”—the linear transfer function $T$, the linear growth function $D$, etc.—that naturally cast the absolute linear power spectrum in power-law $w$CDM-type cosmologies into a separable function:

[TABLE]

following the textbook convention of [34], where ${\cal N}\equiv A_{s}/\omega_{\rm m}^{2}$ is the overall normalisation of the linear matter power spectrum up to an irrelevant multiplicative constant.

Then, for the parameter variations represented by the independent parameters of equation (3.8), it follows from the same reasoning of approximate universality and multiplicability that a multivariate relative power spectrum may be most conveniently approximated by

[TABLE]

where $\Delta X\equiv X-\bar{X}$ with the revised understanding that $X$ may be a model parameter or a linear proxy, $(\delta/\gamma)_{X}$ is the approximate-universal form of $\delta/\gamma$ under variation of $X$ alone, ${\cal Q}\equiv(k/k_{\rm piv})^{n_{s}-1}$, and

[TABLE]

with $\bar{\cal N}\equiv{\cal N}(\bar{A}_{s},\bar{\omega}_{\rm m})$, $\bar{D}\equiv D(\bar{w},\bar{\omega}_{\rm m},\bar{h};a)$, and $\bar{T}\equiv T(\bar{\omega}_{\rm m},\bar{\omega}_{\rm b};k)$, specify the target and the reference cosmology respectively.

3.3 Further remarks

Equation (LABEL:eq:form) serves as a starting point for the construction of a fitting function of the relative power spectrum; Three more remarks are in order before we proceed.

Remark 1: Fitting functions

The salient feature of equation (LABEL:eq:form) is that the cosmological dependence of the relative power spectrum has been largely subsumed by the linear quantities $\Delta X/\bar{X}$. Thus, the task of finding a full fitting function for $\delta({\bf\Theta},{\bf\Theta}_{0};k;a)$ boils down at the most elementary level to writing down a cosmology-independent functional form in terms of the wavenumber $k$ and scale factor $a$ for each of the four familial approximate-universal forms $(\delta/\gamma)_{X}$. For fixed values of $a$ this is a trivial exercise. A more useful endeavour would be to model the approximate-universal forms’ dependence on the scale factor $a$, to be pursued in section 4.

In a more sophisticated model one could of course also incorporate the small, cosmology-dependent deviations from the approximate-universal forms that inevitably creep in at large wavenumbers. We do not however see this as a necessary step at this stage: the one-parameter approximate-universal forms $(\delta/\gamma)_{X}$ are in the worst case $10\to 20$% “off” at $k\gtrsim 4\ h$/Mpc (see figures 6 to 9), while the linear deviations $\Delta X/\bar{X}$ are typically $O(0.1)$. Thus, barring an unfortunate add-up of errors, we can be confident that $\delta$ can be reproduced to $\pm 0.01\to 0.02$ up to $k\sim 10\ h$/Mpc.

Remark 2: Varying the matter density

Equation (LABEL:eq:form) is amenable to further algebraic manipulation, a property that is especially useful in those cases where a cosmological model parameter controls more than one linear proxy. The case in point is the physical matter density $\omega_{\rm m}$, the only parameter that controls the linear transfer function $T$ in the cosmologies under consideration. Because $\omega_{\rm m}$ affects also the linear growth function $D$ and the normalisation ${\cal N}$, it is a priori not possible to establish the approximate-universal form $(\delta/\gamma)_{T}$ directly from a set of $N$-body simulations such as detailed in table 3 that uses $\omega_{\rm m}$ as a base parameter.

However, equation (LABEL:eq:form) permits us to write

[TABLE]

where

[TABLE]

and $\gamma_{\omega_{\rm m}}\equiv\gamma({\bf\Theta}=\{\bar{\theta}_{\omega_{\rm m}};\omega_{\rm m}\},{\bf\Theta}_{0}=\{\bar{\theta}_{\omega_{\rm m}};\bar{\omega}_{\rm m}\})$ denotes the relative linear matter power spectrum under variations in $\omega_{\rm m}$ alone. Then, solving for $(\delta/\gamma)_{T}$ and substituting back into equation (LABEL:eq:form) itself yields an alternative form

[TABLE]

which has the desirable feature that $(\delta/\gamma)_{\omega_{\rm m}}$ can be determined directly from simulations. Indeed, for the cosmologies of table 3, equation (3.13) may be the more convenient albeit less general fitting form than equation (LABEL:eq:form).

Figure 11 shows the exact relative matter power spectrum from the simultaneous variation of $\{\omega_{\rm m},w,A_{s}\}$ and its approximate form constructed from single-parameter variations via equation (3.13), for two target cosmologies $\{\omega_{\rm m}=0.1461,w=-0.85,10^{9}A_{s}=2.300\}$ and $\{\omega_{\rm m}=0.1381,w=-0.85,10^{9}A_{s}=2.100\}$. The agreement is excellent: typically much better than 0.01, and in the worst case $\sim 0.02$ at large $k$ values.

Note that to calculate the linear growth functions $D$ and $D_{1}$ we have solved numerically the differential equation [35]

[TABLE]

with the initial conditions $g(a_{\rm ini})=1$, $g^{\prime}(a_{\rm ini})=0$, and $a_{\rm ini}=10^{-3}$. Here, $g\equiv D/a$, the prime $(\cdots)^{\prime}$ denotes a derivative with respect to the scale factor $a$, $w(a)$ is the dark energy equation of state parameter which may be time-dependent, and

[TABLE]

where the second equality applies in the case $w(a)=w_{0}+w_{a}(1-a)$ [36, 37]. It is usually understood that the solution to equation (3.14) can be approximated to high accuracy by the integral [38]

[TABLE]

where $\Omega(a)$ is the reduced matter density at $a$, and $\rho=0.55+0.05[1+w(z=1)]$ was originally proposed in [38]. Indeed, we have checked that for even time-dependent equations of state, the approximate formula (3.16) is able to reproduce numerical solutions to roughly 1 part in $10^{4}$, sufficient to approximate the growth function differences $\Delta D^{2}/\bar{D}^{2}$ to $O(0.001)$ accuracy for the models tested. Nonetheless, we prefer to err on the side of caution and work directly with the differential equation (3.14).

Remark 3: Varying the scalar spectral index

As already pointed out in section 3.1, because of our choice of pivot scale $k_{\rm piv}=0.05$/Mpc, variation of the scalar spectral index $n_{s}$ introduces a singularity in $(\delta/\gamma)_{n_{s}}$ in the $k$ range of interest. This singularity can be easily removed by recognising that the relative linear power spectrum $\gamma$ under variation of $n_{s}$ alone, $\gamma_{n_{s}}\equiv\Delta{\cal Q}/\bar{\cal Q}$, can be recast as

[TABLE]

where $\Delta n_{s}\equiv n_{s}-\bar{n}_{s}$, and $\Gamma(k)\equiv\Delta n_{s}\sum_{i=0}^{\infty}\left[\Delta n_{s}\ln(k/k_{\rm piv})\right]^{i}/(i+1)!$ is always finite and, at leading order, equal to $\Delta n_{s}$. It then follows that the corresponding approximate-universal form is equivalently

[TABLE]

and hence $(\delta/\gamma)_{n_{s}}\,\Delta{\cal Q}/\bar{\cal Q}=(\delta/\Gamma)_{n_{s}}\Gamma$, where the ratio $(\delta/\Gamma)_{n_{s}}$ must also be approximately universal for all variations of $n_{s}$, albeit better-behaved than the original $(\delta/\gamma)_{n_{s}}$.

Then, applying this understanding to equation (LABEL:eq:form) and its restricted form (3.13), we find respectively

[TABLE]

and

[TABLE]

Our fitting function for the relative matter power spectrum, RelFit, will be based upon these expressions; we shall determine the functional forms for $(\delta/\Gamma)_{n_{s}}$ and $(\delta/\gamma)_{X}$ in section 4.

4 RelFit fitting functions

That the ratio $\delta/\gamma$ should take on an essentially cosmology-independent form under variation of one cosmological model parameter or proxy is perhaps not very surprising upon scrutiny. As the top and middle panels of figures 6 to 8 demonstrate, the current generation of observations on linear scales already constrains cosmology to the extent that $\delta,\gamma,\Delta X/\bar{X}\ll 1$. Such tight constraints imply that perturbing $P({\bf\Theta})$ in $X$ around a reference model ${\bf\Theta}_{0}$ will always yield to leading order in $\Delta X/\bar{X}$ a linear dependence of $\delta$ on $\Delta X/\bar{X}$, i.e.,

[TABLE]

regardless of the exact functional dependence of $P({\bf\Theta};k;a)$ on $X$.

Furthermore, while the functional derivatives $\left.\partial\ln P/\partial\ln X\right|_{{\bf\Theta}={\bf\Theta}_{0}}$ depend in principle on our choice of expansion point ${\bf\Theta}_{0}$, the correction incurred by choosing a different ${\bf\Theta}_{0}$ must be $\ll{\cal O}(1)$ if the new expansion point remains within the observationally allowed range. Thus, in this restricted sense the derivatives $\left.\partial\ln P/\partial\ln X\right|_{{\bf\Theta}={\bf\Theta}_{0}}$ are essentially “approximately universal”, and we identify them with the approximate-universal forms $(\delta/\gamma)_{X}$ defined in section 3.1. Then, to first order in small $\Delta X/\bar{X}$, equations (LABEL:eq:form) and (4.1) are the same.

Identifying the approximate-universal forms $(\delta/\gamma)_{X}$ with finite-difference estimates of the functional derivatives of $P({\bf\Theta};k;a)$ immediately suggests that a reasonable approximation of their functional forms can be established using as few as two simulations per family $X$, where $\Delta X/\bar{X}$ should be chosen to be as close to zero as is permitted by the precision limitations. Then, the full $w$CDM fitting function can in principle be constructed with as few as five simulations in total. Given however that we have already at our disposal a set of some 20 simulations, we opt instead to compute the derivatives based on double-sided estimation, which ups the number of required simulations to nine.

In finding functional forms for $(\delta/\gamma)_{X}$ we adopt a strictly empirical approach and simply match rational functions to our simulated spectra, irrespective of their limiting behaviours on very small scales. This also means that extrapolating RelFit to outside the calibration $k$-region may return nonsensical results. We note however that our simulated $(\delta/\gamma)_{X}$, especially $X={\cal N},D$, appear to exhibit small-scale behaviours consistent with the stable clustering ansatz [39]. This suggests that a physically motivated fitting function could be constructed from, e.g., recasting the well-known Peacock–Dodds fitting formula [40] as a fitting function directly for $(\delta/\gamma)_{X}$. We refer the interested reader to appendix A for details.

4.1 Functional forms for $(\delta/\gamma)_{X}$

Following the findings of section 3, we choose as the independent variable in our fitting functions

[TABLE]

where $k_{\rm peak}$ specifies the locations of the peak features in $(\delta/\gamma)_{{\cal N},\omega_{\rm m},D}$. Interpolating our simulation outputs at $a=0.85,0.7,0.5,0.3$, we find $k_{\rm peak}$ to be well described by

[TABLE]

with an $a$-dependent $k_{\sigma}\equiv 1/x$ defined by the condition

[TABLE]

evaluated for the reference cosmological model. For the reference $\Lambda$CDM cosmology of table 1, $k_{\sigma}/(h/{\rm Mpc})\simeq 0.844,1.11,2.06,7.92$ at $a=0.85,0.7,0.5,0.3$ respectively.

4.1.1 $X={\cal N},\omega_{\rm m},D,n_{s}$

We use a subset of the $N$-body simulation results of table 3 to calibrate $(\delta/\gamma)_{X}$ in the wavenumber range $k=0.05\to 10\ h$/Mpc at $a=0.85,0.7,0.50,0.30$. Specifically, we use relative matter power spectra formed from the pairs:

• •

$X={\cal N}$: {1024$A_{s,l}$, 1024Ref}, {1024$A_{s,h}$, 1024Ref};

• •

$X=\omega_{\rm m}$: {1024$\omega_{{\rm m},l}$, 1024Ref}, {1024$\omega_{{\rm m},h}$, 1024Ref};

• •

$X=D$: {1024$w2$, 1024Ref}, {1024$w4$, 1024Ref};

• •

$X=n_{s}$: {1024$n_{s,l}$, 1024Ref}, {1024$n_{s,h}$, 1024Ref}.

At each scale factor $a$, we construct for each pair the corresponding ratio $\delta/\gamma$ and combine them to form a mean $(\delta/\gamma)_{X}$ for each family $X$ weighted by the inverse of the linear relative power spectrum, $|\gamma|^{-1}$, at that scale factor.

We fit each weighted mean $(\delta/\gamma)_{X}$ using rational functions of quadratic polynomials in $\log_{10}y$, where the fitting coefficients are themselves functions of the scale factor $a$. In all cases, $(\delta/\gamma)_{X}$ must converge to the predictions of linear theory at $k\to 0$, a condition we explicitly enforce in all of our fitting functions by tuning down the rational functions with a $1-e^{-y}$ factor. Specifically, for variations in $X={\cal N},\omega_{\rm m},D$, we use the functional form

[TABLE]

while variations in $n_{s}$ are well described by

[TABLE]

In all cases $X={\cal N},\omega_{\rm m},D,n_{s}$, the coefficients $b_{0,1,2}^{X}=b_{0,1,2}^{X}(a)$ and $c_{1,2}^{X}=c_{0,1,2}^{X}(a)$ are polynomials of the scale factor $a$ alone given in appendix B, and we remind the reader again that no attempts have been made to model the $k\to\infty$ behaviours of the fitting functions.

Figure 12 shows the predictions of the restricted form of RelFit, $\delta_{\rm fit}$, based on equation (3.20), against the relative matter power spectra, $\delta_{\rm sim}$, constructed from the simulations of table 3; figure 13 shows the corresponding fitting errors formed from their differences. The fit is across the board excellent. At $a=0.85$ and for the whole range of wavenumbers explored, no individual error exceeds 0.01 in magnitude for the eight calibration models, or exceeds 0.025 for the remaining 12 models not used in the calibration of RelFit. The fit improves as we move to smaller scale factors: at $a=0.30$, the fitting error is always well below 0.01 for the entire $k$-range.

The reasoning behind RelFit together with the parameter dependence of the linear matter power spectrum in $w$CDM cosmologies, equation (3.8), also suggests that varying the dimensionless Hubble parameter $h$ should produce an effect on the nonlinear matter power spectrum identical to varying the linear growth function $D$. Likewise, RelFit in its present form imposes no restriction on the time-dependence of the dark energy equation of state parameter. These scenarios will be explored further in sections 4.3 and 4.2 respectively.

4.1.2 $X=T$

While none of the pairs of simulations in table 3 models explicitly a variation in the linear transfer function $T$ alone, following the arguments of section 3.3 it is possible to construct a fitting function for $(\delta/\gamma)_{T}$ using a combination of our set of 1024$\omega_{\rm m}$ simulations and the fitting functions derived in section 4.1.1. With $(\delta/\gamma)_{T}$ available, a more general form of RelFit based on equation (LABEL:eq:formns) could be achieved, potentially widening the applicability of the fitting function also to target cosmologies involving variations in the physical baryon density $\omega_{{\rm b}}$ (section 4.3) as well as in the effective number of neutrinos $N_{\rm eff}$ (section 4.2).

Recall that varying $\omega_{\rm m}$ changes simultaneously the normalisation ${\cal N}$, the linear transfer function $T$, and the linear growth function $D$. Then, beginning with the relative nonlinear matter power spectra formed from the pairs {1024$\omega_{{\rm m},l}$, 1024Ref} and {1024$\omega_{{\rm m},h}$, 1024Ref}, a simple procedure based on equation (3.11) can be used to recover $(\delta/\gamma)_{T}$ in each case:

Compute the variation in $D$ due to the change in $\omega_{\rm m}$ alone, and use it in RelFit to calculate the corresponding nonlinear variation; 2. 2.

Repeat the above for the nonlinear variation in ${\cal N}$ due to $\omega_{\rm m}$; 3. 3.

Remove the $D$ and ${\cal N}$ contributions of steps 1 and 2 from the simulated relative nonlinear matter power spectrum via equation (3.11) to form the relative nonlinear power spectrum under variations in the linear transfer function $T$ alone; 4. 4.

Divide the relative nonlinear power spectrum of step 3 through by its linear counterpart to form $(\delta/\gamma)_{T}$.

We have repeated this process for the two pairs of relative power spectra, formed a weighted average as described in section 4.1.1, and fit it using a rational function of the form (4.5). The resulting fitting coefficients $b_{0,1,2}^{T}=b_{0,1,2}^{T}(a)$ and $c_{1,2}^{T}=c_{0,1,2}^{T}(a)$ can be found in appendix B.

Figure 14 shows the fitting errors of the general form of RelFit, equation (LABEL:eq:formns), for the eight 1024$\omega_{{\rm m}}$XXX simulations of table 3 relative to 1024Ref.333Recall that relative matter power spectra formed from pairs of simulations without variations in $\omega_{\rm m}$ are not affected by the choice between the general and the restricted form of RelFit. Again, we see that the fit is across the board excellent, and at $a=0.85,0.70,0.50$, comparable to that of the restricted form (figure 13). At $a=0.30$, however, the general form of RelFit appears to systematically underestimate the simulated power spectra at $k\gtrsim 1\,h$/Mpc by some 0.005 to 0.01. This is likely an artefact of the admittedly convoluted method with which we have extracted $(\delta/\gamma)_{T}$ in this section, and can potentially be improved with calibrations against dedicated simulations in which only the linear transfer function is varied. We shall defer this exercise to a later publication. Suffice it to say here that figure 14 demonstrates the robustness of the general strategy of fitting function construction proposed in this work.

4.2 Application to extended models

The form of RelFit—phrased in terms of variations in the linear transfer function, linear growth function, etc.—suggests that its applicability extends beyond the $w$CDM cosmologies we have used to calibrate its free parameters. In order to test this possibility, we have performed an additional set of simulations using Gadget-2/Camb, detailed in table 4, that go beyond $w$CDM in two different ways: (i) a time-dependent dark energy equation of state parameter $w(a)$, which at the linear level affects only the growth function, and (ii) a linear transfer function modified by a non-canonical effective number of neutrinos $N_{\rm eff}$.

(i) Time-dependent dark energy equation of state

Dynamical dark energy models such as quintessence typically predict effective equations of state for the dark energy component that change with time (e.g., [41]). The exact time dependence varies from model to model. Here, we use for simplicity a time dependence parameterised by [36, 37]

[TABLE]

where we fix $w_{0}=-0.85$, but allow $w_{a}$ to vary in the interval $w_{a}\in[-0.1,0.3]$ in our simulations. Current cosmological measurements do not provide strong constraints on the time dependence of $w(a)$, and the models represented by our choices of $w_{a}$ values, while spanning a parameter range comparable to only about 1.5 times the standard deviation inferred from the 2018 Planck+SNe+BAO data [27], do in fact deviate strongly from the reference $\Lambda$CDM cosmology in their matter power spectrum predictions.

Extending RelFit to include a time-dependent dark energy equation of state parameter simply requires that we redefine the linear growth function variations $\Delta D^{2}$ and $\Delta D_{1}^{2}$ that appear in equations (LABEL:eq:formns) and (3.20) as

[TABLE]

where $\bar{w}(a)$ denotes the reference $\Lambda$CDM choices of $\bar{w}_{0}=-1$ and $\bar{w}_{a}=0$ in equation (4.7).

(ii) Non-canonical effective number of neutrinos

Any light thermal particle species that decouples while relativistic will behave in the cosmological context essentially like a standard-model neutrino, and contribute to the non-photon radiation energy density, conventionally parameterised as the effective number of thermalised neutrinos $N_{\rm eff}$. Well-known examples of such particle species include sterile neutrinos and axions (e.g., [42, 43]).

Phenomenologically, increasing $N_{\rm eff}$ alone shifts the epoch of matter–radiation equality to a lower redshift according to [44],

[TABLE]

where $\omega_{\gamma}$ is the present-day photon energy density. For the linear matter power spectrum, changes in $z_{\rm eq}$ are manifested primarily as a shift in the location of the turning point $k_{\rm eq}$ according to

[TABLE]

which, within the structure of RelFit, is captured by a variation in the linear transfer function. Then, incorporating $N_{\rm eff}$ into RelFit simply requires that we use the general form (LABEL:eq:formns) of the fitting function together with

[TABLE]

where the linear transfer function $T$ is now a function of three cosmological parameters.

The 2018 Planck+external data combination currently constrains $N_{\rm eff}$ most tightly to $N_{\rm eff}=2.99^{+0.34}_{-0.33}$ (95% C.I.) [27]; our two choices of $N_{\rm eff}=3.34$ and $N_{\rm eff}=4.04$ in table 4 therefore represent respectively a $2\sigma$ and a $20\sigma$ variation away from the 2018 Planck best-fit.

Figure 15 shows the predictions of RelFit—as calibrated originally in section 4.1—against the relative matter power spectra constructed from the simulations of table 4, together with the corresponding fitting errors. Again, the differences between the predictions of RelFit and the simulated relative power spectra up to $k\sim 1\,h$/Mpc generally do not exceed about 0.01; in the case of 1024$N_{\rm eff}4.0$, the large fluctuations around zero seen at $k\sim 0.1\to 1\,h$/Mpc are a consequence of nonlinear corrections to the baryon acoustic oscillations, which in principle can be modelled approximately using a suppression factor (as has been implemented in, e.g., HMCode [10], but not in RelFit).

Beyond $k\gtrsim 1\,h$/Mpc the fitting errors tend to increase, although for most $w(a)$ and $N_{\rm eff}$ cosmologies tested here the RelFit predictions still fall within 0.02 of the simulation results. The only exception is the case of 1024$w2w_{a}1$, where at $k\gtrsim 4\,h$/Mpc the deviation is up to 0.03. We note however that the particular $w(a)$ cosmology represented by this simulation is fairly far away from the $\Lambda$CDM reference cosmological model in terms of the deviation of its linear matter power spectrum from the reference case ($\gtrsim 15$% at $a\geq 0.70$). Given the “perturbative” nature of RelFit, it is perhaps not surprising that its simple linear prescription should break down at large wavenumbers.

We conclude section 4.2 with the emphasis that none of the simulations of table 4 has been used to calibrate RelFit. In particular, the $(\delta/\gamma)_{T}$ fitting function that forms the basis of the RelFit predictions in the two $N_{\rm eff}$ scenarios has been extracted from a combination of target cosmology simulations that have nothing to do with varying $N_{\rm eff}$ at face value. That RelFit is still capable of predicting to $0.01\to 0.02$ the relative power spectra of these target cosmologies speaks again for the general soundness of our strategy.

4.3 Comparison with CosmicEmu, Halofit, and HMCode

4.3.1 Single-parameter variations

The essence of RelFit is a set of first-order logarithmic functional derivatives of the nonlinear matter power spectrum $P({\bf\Theta};k;a)$ with respect to variations in the linear matter power spectrum $P_{\rm L}({\bf\Theta};k;a)$ evaluated at the reference cosmology ${\bf\Theta}={\bf\Theta}_{0}$. Predicting a target nonlinear $P({\bf\Theta};k;a)$ relative to the reference $P({\bf\Theta}_{0};k;a)$ simply consists in multiplying these derivatives with the relevant variations in the linear $P_{\rm L}({\bf\Theta};k;a)$ away from the reference $P_{\rm L}({\bf\Theta}_{0};k;a)$. One immediately concludes that the smaller the linear variations a target cosmology produces, the higher the fidelity of RelFit in predicting its nonlinear variations.

We take as a formal assessment of “smallness” the fractional variation in the linear matter power spectrum at the “peak” wavenumber $k_{\rm peak}$, defined in equation (4.3), of the reference cosmology. Then, at $a=0.85$ and for single-parameter variations, a maximum 10% (15%) variation corresponds to target cosmological parameter values falling in the region

[TABLE]

where the equivalent $\sigma_{8}$ range assumes all parameters but $A_{s}$ held fixed at their reference values, and we have included in this list the (as-yet-unexplored) Hubble parameter $h$ and physical baryon density $\omega_{{\rm b}}$. Where applicable the parameter region (4.12) is larger than that of equation (3.1) used to establish RelFit, while the $\omega_{\rm b}$ range, representing $2\to-5$% variations in the linear matter power spectrum, has been chosen so as to stay within the confines of the Mira–Titan simulations [6, 7].444In the same vein, $10^{9}A_{s}=2.495$, or equivalently, $\sigma_{8}=0.9$, represents only a 13.4% variation from the reference cosmology, and $w=-1.3$ only 11.4% at $a=0.85$. Simple power counting then suggests that the output of RelFit in the region (4.12) should be accurate to $\lesssim 0.01$ ($\lesssim 0.02$).

The left panel of figure 16 compares the output of RelFit in the parameter region (4.12) at the calibration scale factors $a=0.85,0.50$, with the predictions of the CosmicEmu emulator trained on the Mira–Titan simulations [6, 7]. For comparable cosmological parameters, CosmicEmu interpolates in the parameter region

[TABLE]

which, with the exception of $\omega_{\rm m}$ and $\sigma_{8}$, is marginally larger than the “15%-variation” parameter region defined in equation (4.12). As can be seen, the agreement between RelFit and CosmicEmu in the region (4.12) is remarkable: with few exceptions, the two sets of predictions agree to 0.01 (0.02) or better across the whole wavenumber range tested.

The same comparison at the “off-calibration” scale factors $a=0.60,0.40$ is shown in the right panel of figure 16, which serves to test the $a$-dependence of the fitting coefficients presented in appendix B. At $a=0.60$ the agreement between RelFit and CosmicEmu is as good as or at marginally worse than the “on-calibration” comparisons discussed above. The $a=0.40$ results are likewise concordant for variations in $\omega_{\rm m}$, $n_{s}$, $w$, $h$, and $\omega_{\rm b}$ across the whole $k$-range, but appear to diverge at $k\sim 2\ h$/Mpc by as much as 4% for variations in $A_{s}$. This may be an error of interpolation in RelFit consequent to a sparsely sampled $a$-space—recall that we have calibrated RelFit at only four instances ($a=0.85,0.70,0.50,0.30$). Interestingly, however, while CosmicEmu uses eight samples in a similar timeframe ($a=1.0,0.91,0.81,0.70,0.60,0.50,0.38,0.33$), the particular instance of $a=0.40$ is likewise off-calibration. To pin down the exact source of discrepancy would require new simulations, which we defer to a later publication.

Lastly, while it may be tempting to interpret figure 16 as an accuracy test of RelFit, it must be kept in mind that CosmicEmu itself has a claimed error margin of 4% on the absolute power spectrum [7]. Likewise, relative power spectra formed from its output are in some cases—particularly when the target and reference cosmologies are far apart—demonstrably erroneous by up to 2% as $k\to 0$, due to convergence to linear perturbation theory not having been explicitly enforced in the emulation process (in contrast to the calibration of RelFit, which does respect convergence to linear theory). Nonetheless, it is encouraging that agreement to $0.01\to 0.02$ or better can be achieved in a fairly broad parameter region, especially given that RelFit and CosmicEmu have been calibrated against completely independent simulations generated from two different $N$-body codes.

4.3.2 Multi-parameter variations

Next we test RelFit against CosmicEmu in the full 6-parameter space of equation (4.12). To do so we draw 10 sets of six random numbers on a 5-sphere of unit radius centred on the origin, where each axis represents a cosmological parameter direction. These random numbers are then rescaled according to the parameter ranges of equation (4.12), assuming that the reference $\Lambda$CDM cosmology sits at the centre of the sphere. Table 5 shows the 10 target cosmologies sampled in this manner. The sampling procedure ensures that all 10 target cosmologies reside on the “surface of 15% variation” away from the reference, where we generically expect the fitting error of RelFit to be the largest—up to $\sim 0.02$ by the arguments of section 4.3.1; by the same token we can expect the errors of RelFit to be smaller than $\sim 0.02$ for those cosmologies contained within the surface.

Figure 17 compares the output of RelFit for the 10 target cosmologies of table 5 at $a=0.85,0.50$ with the predictions of CosmicEmu. As with the single-parameter comparisons of section 4.3.1, the consistency between the two sets of predictions is remarkable: for the most part the output of RelFit is within 0.01 of the CosmicEmu predictions for the whole range of wavenumbers tested, and offers a clearly better concordance than can be achieved with HMCode (2016 version) [11] and especially Halofit [8] as updated in [9] for the same set of target cosmologies, also shown in figure 17.

The “worst-performing” cosmology of the RelFit set has a maximum deviation from the CosmicEmu predictions of 0.03 at $k\sim 1.5\,\bar{h}$/Mpc and $a=0.50$. Interestingly, the same cosmology also exhibits the largest deviation from CosmicEmu under corrections with both HMCode and Halofit. Bearing in mind that CosmicEmu has a claimed accuracy of 4% [7], these deviations could suggest that the inaccuracy lies with CosmicEmu itself rather than with the three fitting functions. It is likewise intriguing that except at $k\lesssim 0.5\,\bar{h}$/Mpc, the agreement between RelFit and CosmicEmu does not improve with a decreasing scale factor $a$, in contrast to the fidelity of RelFit to simulation results, which, as shown in figures 13 and 14, does improve significantly from $a=0.85$ to $a=0.30$ across the board. Further investigation of these oddities is however beyond the scope of the present work.

We conclude our study with a comparison of an alternative calibration of RelFit—against only five simulations,555The simulations used are 1024$A_{s,l}$, 1024$w4$, 1024$\omega_{{\rm m},l}$, 1024$n_{s,{\rm h}}$, 1024Ref of table 3. the minimum required to map out the four approximate-universal forms $(\delta/\gamma)_{X}$—to CosmicEmu, using again the 10 target cosmologies of table 5. This comparison is shown in the bottom panels of figure 17 as “RelFit (5 simulations)”. Evidently, this even “cheaper” version of RelFit agrees with CosmicEmu almost as well as the default version (calibrated against nine simulations), with only marginal deteriorations (and possibly a hint of systematic bias) in the agreement at $k\gtrsim 1\,\bar{h}$/Mpc and still outperforming both HMCode and Halofit. While we do not advocate this alternative calibration because of potential biases introduced by the one-sided derivative estimates (see section 4.1), this exercise serves to illustrate succinctly the power of the RelFit method, and supports our thesis that an accurate fitting function to the relative nonlinear matter power spectrum can indeed be obtained very cheaply.

5 Conclusions

The central message of this work is twofold: (i) The relative matter power spectrum, defined as the fractional deviation in the absolute matter power spectrum produced by a target cosmology away from a reference $\Lambda$CDM prediction, is fairly insensitive to the specifics of a simulation and can be computed to 1%-level accuracy at a much lower computational cost than can the absolute matter power spectrum itself. (ii) Within the $w$CDM class of models tested, the relative nonlinear power spectrum has the interesting property that when divided through by its linear counterpart under single-parameter variations, the result exhibits an approximate universality for each class of variations at the onset of nonlinearity. Exploiting this and the property of multiplicability of the relative power spectrum, it is possible to construct full fitting functions to any cosmology in the vicinity of $\Lambda$CDM in a piece-wise manner, whereby component fitting functions are sought for single-parameter variations and then multiplied together to form the full fitting function.

Point 1 offers an advantage in the exploration of the nonlinear matter power spectrum under variation of cosmology, in that once an ultra-precise reference absolute matter power spectrum has been computed, variations away from the reference cosmology can be investigated at a relatively low cost, enabling a larger swath of parameter space to be explored or a particular parameter region of interest to be more densely sampled. Point 2 enables independent, piece-wise studies of cosmological models on nonlinear scales, by which we mean a fitting function for variations in, e.g., the primordial curvature power spectrum can be constructed independently of that for variations in, e.g., the dark energy equation of state. Both have particular implications for the investigation of “non-standard” or “exotic” cosmologies: Because computational costs have been significantly reduced, the task of exploring exotic model parameter spaces is now possible for a much wider section of the scientific community. Computing the nonlinear matter power spectrum at 1%-level accuracy can be made a far more egalitarian exercise than is currently feasible with conventional methods.

As an illustration of the approach, we have used nine relatively inexpensive $w$CDM simulations (box length $L=256\,h^{-1}$Mpc and $N=1024^{3}$ particles, initialised at $z_{i}=49$) spanning the parameter directions $\{\omega_{\rm m},A_{s},n_{s},w\}$ to construct the fitting function RelFit that is able to reproduce to $0.01\to 0.02$ accuracy or better the relative nonlinear matter power spectra of 20-odd $w$CDM cosmologies at $0.85\geq a\geq 0.30$ up to $k\simeq 10\,h$/Mpc. RelFit is likewise applicable—without modification and to the same accuracy—to cosmologies in which $w(a)=w_{0}+w_{a}(1-a)$ parameterises a time-dependent dark energy equation of state, and where $N_{\rm eff}$ may deviate from the canonical $N_{\rm eff}=3.04$.

Testing RelFit against the output of the CosmicEmu emulator trained on the Mira–Titan simulations [6, 7], we find again consistency at better than $0.01\to 0.02$ in a large region of the 6-parameter space $\{\omega_{\rm m},A_{s},n_{s},,w,\omega_{b},h\}$, despite RelFit not having been calibrated against the same simulations—or any high-quality simulation for that matter. For the set of 10 randomly selected cosmologies examined, the ability of RelFit to replicate the CosmicEmu predictions surpasses that of both Halofit [8] (with updates [9]) and HMCode (2016 version) [11]. The same success can be reproduced even with only five calibrating simulations, although for reasons of minimising potential systematic biases, the nine-simulation calibration is preferable—this version of RelFit is summarised in appendix B.

To conclude, the relative matter power spectrum is an inexpensive and democratically accessible route to fulfilling the 1%-level accuracy demands of the forthcoming generation of large-scale structure probes. Our prototype fitting function RelFit for $w(a)$CDM+$N_{\rm eff}$ cosmologies, which takes the linear matter power spectrum as an input, can be readily implemented in publicly available linear Boltzmann codes such as Camb [29] and Class [32] together with, e.g., an output of CosmicEmu as a placeholder for the ultra-precise reference absolute power spectrum yet to come. In the future we shall extend the approach to cosmologies including massive neutrinos, as well as more “exotic” scenarios such as decaying dark matter, interacting dark matter, and dark energy perturbations.

Acknowledgments

Y3W thanks Amol Upadhye for useful discussions about CosmicEmu. STH is supported by a grant from the Villum Foundation. Y3W is supported in part by the Australian Research Council’s Discovery Project (project DP170102382) and Future Fellowship (project FT180100031) funding schemes.

Appendix A Connection to the stable clustering ansatz

The Peacock–Dodds (PD) fitting formula renders the dimensionless nonlinear matter power spectrum, $\Delta_{\rm NL}^{2}({\bf\Theta};k;a)\equiv k^{3}P({\bf\Theta};k;a)/(2\pi^{2})$, in terms of a simple function of the dimensionless linear power spectrum, $\Delta_{\rm L}^{2}({\bf\Theta};k;a)\equiv k^{3}P_{\rm L}({\bf\Theta};k;a)/(2\pi^{2})$, and the $k$-independent linear growth function $D$ [40]:

[TABLE]

Here, the nonlinear and linear wavenumbers, $k_{\rm NL}$ and $k_{\rm L}$, satisfy a scaling relation

[TABLE]

following from the stable clustering ansatz [39, 40]. The dependences of $f$ on $\Delta_{\rm L}^{2}$ and $D$ on strongly nonlinear scales are likewise stipulated by stable clustering. The general form of $f$, however, must be determined from and calibrated against simulations, which introduces an additional, empirical dependence on the effective linear spectral index

[TABLE]

in the fitting coefficients of the PD formula that we shall ignore for now.

Then, perturbing around ${\bf\Theta}_{0}$ it is straightforward to establish that

[TABLE]

to leading order in $\delta,\gamma,\Delta D/\bar{D}$, with expansion coefficients

[TABLE]

Note that with the exception of $U$, all quantities on the RHS of equation (A.4) are technically functions of $k_{\rm L}$, and must be mapped to $k_{\rm NL}$ via equation (A.2) before they can be meaningfully interpreted.

Figure 18 shows the expansion coefficients as functions of $k_{\rm NL}$ for the canonical reference $\Lambda$CDM model at $a=0.85,0.70,0.50,0.30$. We have used the (cosmology- and simulation-dependent) PD fitting formula [40] to make this plot.666Note that in evaluating the expansion coefficients in figure 18, we have applied the PD formula to a “dewiggled” linear power spectrum, in order to show more clearly the coefficients’ broad-band characteristics. However, a number of cosmology- and simulation-independent features can still be identified:

At any given epoch, the function $R$ generally evolves from unity in the linear regime to a peak around $k_{\rm peak}$ (see equation (4.3)) at which structures begin to acquire nonlinearities. It then drops to a constant asymptotic value $R_{\infty}$; here, the stable clustering ansatz stipulates that

[TABLE]

on very nonlinear scales, leading to $R_{\infty}=3/2$. 2. 2.

In the linear regime the function $f$ approaches $\Delta^{2}_{\rm L}$ and does not depend on the growth factor $D$; in this limit $V$ asymptotes to zero. In the nonlinear limit the assumption of stable clustering and hence equation (A.8) lead to $V_{\infty}=-3$. The transition between the two regimes again takes place around $k_{\rm peak}$. Note that $V$ is always negative in the PD fitting formula. 3. 3.

Barring baryon acoustic oscillations, the effective linear spectral index $n(k_{\rm L})$ is always positive for all cosmologies consistent with current observations; it takes on a value $\sim 4$ at wavenumbers $k\lesssim k_{\rm eq}$ defined in equation (4.10), and decreases monotonically beyond $k\sim k_{\rm eq}$ to zero on strongly nonlinear scales. 4. 4.

The $U$ term traces its origin to the cosmology-dependent mapping of the nonlinear wavenumber $k_{\rm NL}$ to the linear wavenumber $k_{\rm L}$ via equation (A.2). For $\Lambda$CDM-type cosmologies, however, we expect $U$ to increase monotonically from zero in the linear regime to $U_{\infty}=1/3$ in the strongly nonlinear regime.

With these considerations in mind it is straightforward to establish the qualitative behaviours of some of the approximate-universal forms $(\delta/\gamma)_{X}$. For $X={\cal N}$, e.g., from varying only the primordial fluctuation amplitude $A_{s}$ between the reference and target cosmologies, the relative linear power spectrum $\gamma$ is $k$-independent and $\Delta D/\bar{D}=0$. Then, we find from equation (A.4)

[TABLE]

where we can immediately conclude from the first approximate equality that the relative nonlinear power spectrum $\delta$ at $k_{\rm NL}$ is always larger in magnitude than the corresponding $\gamma$ at $k_{\rm L}$, i.e., $(\delta/\gamma)_{\cal N}\gtrsim 1$. In the stable clustering limit we expect $(\delta/\gamma)_{\cal N}$ to tend to the expression in the second line of equation (A.9), which for our reference $\Lambda$CDM cosmology evaluates numerically to $(\delta/\gamma)_{\cal N}\simeq 1.2$ at $k\sim 10\,h$/Mpc, a prediction that appears to be borne out by our simulation results at $a=0.85,0.50$ shown in figures 8 and 9.

In the case of $X=D$, e.g., from varying the dark energy equation of state parameter $w$, we have again a $k$-independent $\gamma=2\,\Delta D/\bar{D}$, such that equation (A.4) becomes

[TABLE]

Recall that on linear scales $V$ evaluates to zero. Thus, we expect $(\delta/\gamma)_{D}$ to evolve on mildly nonlinear scales approximately like $(\delta/\gamma)_{\cal N}\gtrsim 1$. As we move further into the nonlinear regime, however, $V(g)$ becomes increasingly negative; when the condition $V<2(1-R+R\,U\,n)$ is satisfied, $(\delta/\gamma)_{D}$ flips from $\gtrsim 1$ to $\lesssim 1$, also observed in figures 6 and 7. In the stable clustering limit, $(\delta/\gamma)_{D}$ should vanish exactly (to all orders): our simulation results in figures 6 and 7 appear to suggest this behaviour at large wavenumbers $k$. The simulations do not however have the necessary dynamical range to fully resolve this limit.

Appendix B RelFit fitting coefficients

The fitting function RelFit for the relative nonlinear matter power spectrum in $w$CDM cosmologies comes in a general form,

[TABLE]

and a restricted form,

[TABLE]

the latter of which applies to a restricted set of $w$CDM parameters.

Here, the target and reference $\Lambda$CDM cosmologies are specified respectively by

[TABLE]

where ${\cal N}\equiv A_{s}/\omega_{\rm m}^{2}$ is the overall normalisation of the linear matter power spectrum, $D$ is the linear growth function, $T$ the linear transfer function, and $\Delta X\equiv X-\bar{X}$ denotes the variation in $X={\cal N},D,T$ between the target and reference cosmologies. The function

[TABLE]

specifies the variation in the shape of the primordial curvature power spectrum, taken to be of a power-law form, with $\Delta n_{s}\equiv n_{s}-\bar{n}_{s}$ for a $k$-independent $n_{s}$. In the case of the restricted form of RelFit, which applies if the only parameter varied in the linear transfer function is the physical matter density $\omega_{\rm m}$, we require also the auxiliary definitions

[TABLE]

where the last entry denotes the relative linear power spectrum in which only the physical matter density $\omega_{\rm m}$ is varied away from its reference $\Lambda$CDM value.

The crux of RelFit are the approximate-universal forms $(\delta/\gamma)_{X}$ and $(\delta/\Gamma)_{n_{s}}$. For $X={\cal N},\omega_{\rm m},D,T$, these are given by

[TABLE]

while for $X=n_{s}$ we have

[TABLE]

with the $\Gamma(k)$ series truncated at $i=2$. Here, the independent variable is

[TABLE]

where $k_{\sigma}\equiv 1/x$ is defined by the condition

[TABLE]

The corresponding $a$-dependent coefficients, calibrated against nine simulations at output scale factors $a=0.85,0.70,0.50,0.30$ are as follows:

• •

$X={\cal N}$:

[TABLE]

• •

$X=\omega_{\rm m}$:

[TABLE]

• •

$X=D$:

[TABLE]

• •

$X=T$:

[TABLE]

• •

$X=n_{s}$:

[TABLE]