#  $k$-e$\mu$lator: emulating clustering effects of the $k$-essence dark
energy

A. R. Nouri-Zonoz,1 Farbod Hassani,2 Martin Kunz1
1Université de Genève, Département de Physique Théorique and CAP, Genève 4 ,
Switzerland
2Institute of Theoretical Astrophysics, University of Oslo, 0315 Oslo, Norway
E-mail<EMAIL_ADDRESS>0009-0006-6164-8670 0000-0003-2640-4460
0000-0002-3052-7394

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We build an emulator based on the polynomial chaos expansion (PCE) technique
to efficiently model the non-linear effects associated with the clustering of
the $k$-essence dark energy in the effective field theory (EFT) framework.
These effects can be described through a modification of Poisson’s equation,
denoted by the function $\mu(k,z)$, which in general depends on wavenumber $k$
and redshift $z$. To emulate this function, we perform $200$ high-resolution
$N$-body simulations sampled from a seven-dimensional parameter space with the
Latin hypercube method. These simulations are executed using the k-evolution
code on a fixed mesh, containing $1200^{3}$ dark matter particles within a box
size of $400~{}\text{Mpc}/h$. The emulation process has been carried out
within UQLab, a MATLAB-based software specifically dedicated to emulation and
uncertainty quantification tasks. Apart from its role in emulation, the PCE
method also facilitates the measurement of Sobol indices, enabling us to
assess the relative impact of each cosmological parameter on the $\mu$
function. Our results show that the PCE-based emulator efficiently and
accurately reflects the behavior of the $k$-essence dark energy for the
cosmological parameter space defined by $w_{0}c_{s}^{2}\text{CDM}+\sum
m_{\nu}$. Compared against actual simulations, the emulator achieves sub-
percent accuracy up to the wavenumber $k\approx 9.4~{}h\text{Mpc}^{-1}$ for
redshifts $z\leq 3$. Our emulator provides an efficient and reliable tool for
Markov chain Monte Carlo (MCMC) analysis, and its capability to closely mimic
the properties of the $k$-essence dark energy makes it a crucial component in
Bayesian parameter estimations. The code is publicly available at
https://github.com/anourizo/k-emulator.

###### keywords:

cosmology: dark energy – cosmology: cosmological parameters – methods:
numerical – methods: statistical

††pubyear: 2024††pagerange: $k$-e$\mu$lator: emulating clustering effects of
the $k$-essence dark energy–C

##  1 Introduction

Based on the theory of general relativity, the $\Lambda$CDM cosmological
model, which incorporates the collisionless cold dark matter (CDM) and the
cosmological constant $\Lambda$ associated with dark energy, has long been the
standard paradigm in describing the large-scale structure, evolution, and the
late-time accelerating expansion of the Universe. Strengthened by its
remarkable consistency with diverse observational data, from the analysis of
cosmic microwave background (CMB) power spectrum (Aghanim et al., 2016) to the
large-scale galaxy clustering surveys (Alam et al., 2017), it has earned
widespread acceptance within the cosmology community. However, despite the
many triumphs of the model, the characterization of the cosmological constant
as dark energy has encountered multiple challenges. Among them, the
cosmological constant problem (Adler et al., 1995) and the cosmic coincidence
problem (Martin, 2012) stand as notable issues within the $\Lambda$CDM
cosmological framework. Faced with these persistent issues and several
tensions within the $\Lambda$CDM context, particularly the divergent
measurements of the Hubble constant ($H_{0}$) and the structure growth
parameter ($S_{8}$) between Planck satellite data (Aghanim et al., 2020) and
local observations—such as supernova analyses for $H_{0}$ (Riess et al., 2021)
and weak lensing experiments for $S_{8}$ (Hildebrandt et al., 2017, 2020;
Abbott et al., 2018; Asgari et al., 2021)— the investigation and test of
modified gravity and dark energy theories has become increasingly relevant.
These pursuits are now more feasible than ever, thanks to the advancements in
cosmological instrumentation. Such developments have brought us into a crucial
phase of research, allowing us to constrain different dark energy candidates
with an unprecedented accuracy. At the forefront of these instruments, future
and ongoing surveys such as Euclid111https://sci.esa.int/web/euclid (Laureijs
et al., 2011), the Vera C. Rubin observatory222https://rubinobservatory.org/
(Abell et al., 2009), DESI333https://www.desi.lbl.gov/ (Aghamousa et al.,
2016), and the Nancy Grace Roman Space
Telescope444https://science.nasa.gov/mission/roman-space-telescope/ (Akeson et
al., 2019) will allow us to shrink the error bars of DE parameters through the
combined use of galaxy clustering and weak lensing methods. This will be of
particular interest in the $k$-essence dark energy models (Armendariz-Picon et
al., 2000) where the current data offers weak constraints on the speed of
sound (see e.g. Ade et al., 2016).

As the Universe evolves, the growth of structures on small scales leads to
non-linear structure formation in cosmic evolution where linear perturbation
theory is no longer applicable. This complexity has made it challenging to
simply rely on traditional theoretical approaches, especially considering that
even the effective field theory of large scale structures (EFTofLSS) fails
beyond $k\sim 0.6~{}\text{h/Mpc}$ (Carrasco et al., 2014). Within this
landscape, $N$-body simulations have emerged as indispensable tools. In
particular, the gevolution code (Adamek et al., 2016) is the first $N$-body
simulation for cosmic structure formation that is based on general relativity
employing the weak field approximation. This code is structured in such a way
that facilitates the implementation of dark energy and modified gravity
theories. An example is the implementation of the $k$-essence dark energy
model in the form of the k-evolution code (Hassani et al., 2019).

Yet, while $N$-body simulations offer a reliable tool to test different models
of gravity/dark energy, one cannot ignore the significant computational
demands they impose. The high resolution requirements for accurate predictions
makes these simulations extremely resource-intensive, both in terms of
computational power and time. Running multiple simulations for parameter
exploration can quickly become infeasible, especially in the context of
Bayesian best fit parameter analysis, where thousands of model evaluations are
typically necessary. This computational challenge has given rise to an
increasing interest in the development and application of emulators.

Emulators are fast and accurate approximations of full simulations, trained on
a subset of the parameter space. As surrogate models, they can quickly provide
estimates for different parameter sets without bearing the substantial
computational expense, making them ideal for extensive parameter searches.
Their ability to boost the efficiency of likelihood sampling is crucial for
the optimal use of Markov chain Monte Carlo (MCMC) techniques. Examples of
cosmological emulators include EuclidEmulator (Knabenhans et al., 2019, 2021),
Aemulus project (DeRose et al., 2019; McClintock et al., 2019; Zhai et al.,
2019; Storey-Fisher et al., 2024), FrankenEmu (Heitmann et al., 2014),
CosmicEmu (Heitmann et al., 2009, 2010; Lawrence et al., 2010; Heitmann et
al., 2016; Lawrence et al., 2017; Moran et al., 2023), BE-HaPPY (Valcin et
al., 2019), NGenHalofit (Smith & Angulo, 2019), Dark quest emulator
(Nishimichi et al., 2019), CosmoPower (Spurio Mancini et al., 2022) and
CosmoPower-JAX (Piras & Spurio Mancini, 2023).

In this paper we introduce a new emulator for the $k$-essence dark energy
model based on k-evolution simulations, assuming constant values for the
equation of state parameter $w_{0}$ and the speed of sound $c_{s}$. The key
aim of this emulator is to accurately reproduce the clustering effects of the
$k$-essence dark energy for small values of sound speed, which is a distinct
feature of the k-evolution code. As we show later, when the speed of sound is
large, the results obtained from Boltzmann codes remain reliable.

In the construction of the emulator, we use sparse polynomial chaos expansion
(SPCE) method, as described in detail in Blatman (2009). This process is
conducted using UQLab, a specialized software for uncertainty quantification
(UQ) introduced in Marelli & Sudret (2014). The effectiveness of SPCE is not
only highlighted in our work but has also been prominently featured in the
publications of the EuclidEmulator (Knabenhans et al., 2019, 2021).

This paper is organized as follows: In Section 2 we discuss the theoretical
framework of the $k$-essence dark energy and explain its clustering impact on
Poisson’s equation through the $\mu$ function, which is the target of the our
emulator. In Section 3, we outline the main steps of the emulator construction
and discuss each one in detail. Following this, in Section 4, we explain the
fine-tuning procedure of parameters involved in the construction of the
emulator, aiming to minimize emulation errors. Section 5 is dedicated to the
sensitivity analysis conducted using Sobol indices with the purpose of
identifying the key parameters influencing the emulator’s target, in this case
the $\mu$ function. In Section 6, we transfer the knowledge obtained from the
mock emulator (Halofit-based emulator) to the development of a simulation-
based emulator. This section also covers discussions on convergence tests of
simulations and preprocessing of the training data. Finally, in Section 7, we
explain the application of our emulator in obtaining the gravitational
potential power spectrum using the matter power spectrum. The paper concludes
with a summary of our results.

##  2 Theoretical framework

###  2.1 The $k$-essence dark energy model

$k$-essence models are a family of theories based on a scalar field where a
generic non-canonical Lagrangian term, often denoted as $P(\phi,X)$, governs
the scalar field’s dynamics. Here $\phi$ is the $k$-essence scalar field and
$X=-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$ is the kinetic
component. The action for the $k$-essence models is written as

$S=\int d^{4}x\sqrt{-g}\left(P(\phi,X)+\frac{R}{16\pi
G}+\mathcal{L}_{m}\right),$ (1)

where $g$ is the determinant of the metric $g_{\mu\nu}$, $P(\phi,X)$ is the
generic form of the $k$-essence Lagrangian, $R$ is the Ricci scalar and
$\mathcal{L}_{m}$ is the matter Lagrangian. This form of action was first
proposed by Armendariz-Picon et al. (1999) as an alternative to the standard
slow-roll scenario, describing the inflationary phase of the Universe. Later
on, the scheme’s initial adaptation to dark energy was studied in Chiba et al.
(2000) and Armendariz-Picon et al. (2000).
To study inhomogeneities in the Universe it is common to adopt the Friedman-
Lemaître-Robertson-Walker (FLRW) metric in the Poisson gauge which reads

$ds^{2}=a^{2}(\eta)[-e^{2\Psi}d\eta^{2}-2B_{i}dx^{i}d\eta+(e^{-2\Phi}\delta_{ij}+h_{ij})dx^{i}dx^{j}],$
(2)

where $a(\eta)$ is the scale factor, $d\eta=dt/a(t)$ is the differential
element of conformal time, $\Psi$ and $\Phi$ denote the temporal and spatial
scalar perturbations, and $B_{i}$ and $h_{ij}$ represent the transverse vector
($\partial^{i}B_{i}=0$) and transverse-traceless tensor
($\partial^{i}h^{j}_{i}=h^{i}_{i}=0$) perturbations, were we have used the
Einstein summation convention over the repeated indices.

From the action defined in (1), the stress-energy tensor of the $k$-essence
scalar field is calculated as,

$\displaystyle T^{(\phi)}_{\mu\nu}$
$\displaystyle=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}P)}{\delta
g^{\mu\nu}}$ (3)
$\displaystyle=-\frac{2}{\sqrt{-g}}\left(\sqrt{-g}\frac{\partial P}{\partial
X}\frac{\delta X}{\delta g^{\mu\nu}}+\frac{\delta\sqrt{-g}}{\delta
g^{\mu\nu}}P\right)$
$\displaystyle=-\frac{2}{\sqrt{-g}}\left(-\frac{\sqrt{-g}}{2}P_{,X}\partial_{\mu}\phi\partial_{\nu}\phi-\frac{\sqrt{-g}}{2}g_{\mu\nu}P\right)$
$\displaystyle=P_{,X}\partial_{\mu}{\phi}\partial_{\nu}{\phi}+g_{\mu\nu}P.$

By setting

$\rho_{\phi}=2XP_{,X}-P,\qquad P_{\phi}=P,\qquad
u_{\mu}=\frac{\partial_{\mu}\phi}{\sqrt{2X}},$ (4)

it is possible to express the given stress-energy tensor above in the form of
a perfect fluid,
$T_{\mu\nu}=(\rho_{\phi}+P_{\phi})u_{\mu}u_{\nu}+g_{\mu\nu}P_{\phi}$.
Following the expressions in (4), the stress-energy tensor of the $k$-essence
fluid, which lacks anisotropic stress ($\sigma$), can be characterised by its
equation of state $w_{\phi}$ and squared speed of sound $c_{s}^{2}$ (for more
details see Hassani et al., 2019)

$w_{\phi}=\frac{P_{\phi}}{\rho_{\phi}}=\frac{P}{2XP_{,X}-P}\quad,\quad
c_{s}^{2}=\left(\frac{\partial
P}{\partial\rho}\right)\bigg{|}_{\phi}=\frac{P_{,X}}{P_{,X}+2XP_{,XX}},$ (5)

where the comma stands for the partial derivative while the symbol $|_{\phi}$
denotes the derivative when $\phi$ is held constant. Note that within the
context of $k$-essence models, both the equation of state parameter,
$w_{\phi}$, and the speed of sound, $c_{s}$, are time-dependent in general.
However, for the purposes of our analysis, we will treat them as constants.
Allowing the equation of state parameter and the speed of sound to vary with
time leads to an increase in the number of free parameters, which complicates
both the theoretical framework and the numerical analysis. Moreover, from the
theoretical point of view, one can parametrize the time evolution of
$w_{\phi},c_{s}^{2}$ in different ways by considering assumptions for the
function $P(X,\phi)$. This topic is interesting in itself and should be
considered in later investigations. To sidestep these complexities, which
require detailed analysis, we have chosen to defer the exploration of time-
varying parameters to future studies. Generally, when it comes to studying the
behaviour of the $k$-essence dark energy, both $c_{s}^{2}$ and $w_{\phi}$ are
essential parameters to look into (hereafter, we will use $w_{0}$ instead of
$w_{\phi,0}$ to represent today’s value of the equation of state parameter).
It is also worth noting that in our study we use the total stress-energy
tensor which includes contributions from the matter sector and the $k$-essence
field

$T_{\mu\nu}=T_{\mu\nu}^{(m)}+T_{\mu\nu}^{(\phi)}.$ (6)

###  2.2 Impact of dark energy perturbations on Poisson’s equation

We can obtain the gravitational field equations by variying the action (1)
with respect to the metric. In the weak field regime, one equation is the
Hamiltonian constraint, which can be seen as a relativistic Poisson equation
(Hassani et al., 2020),

$\nabla^{2}\Phi=4\pi Ga^{2}\sum_{X}\bar{\rho}_{X}\Delta_{X}+S\,,$ (7)

where

$S\equiv\frac{1}{2}\delta^{ij}\Phi_{,i}\Phi_{,j}-8\pi
Ga^{2}\Phi\sum_{X}\bar{\rho}_{X}\Delta_{X}$ (8)

is called the shortwave correction and is due to non-linear contributions in
the weak field expansion. As shown in Hassani et al. (2020), the shortwave
correction is small and can safely be ignored. $\bar{\rho}_{X}$ and
$\Delta_{X}$ are the background density and the comoving density contrast for
each species including cold dark matter, baryon and $k$-essence dark energy
(as our focus is on the late time stage of the Universe, we will neglect the
contribution from radiation). Switching to Fourier space, contributions to the
$\Phi$ metric perturbation from a clustering dark energy component like
$k$-essence can be captured through a function $\mu(k,z)$,

$-k^{2}\Phi=4\pi
Ga^{2}\mu(k,z)\sum_{X\backslash_{\text{DE}}}\bar{\rho}_{X}\Delta_{X}\,,$ (9)

where the symbol $X\backslash_{\text{DE}}$ means the summation excludes the
contribution from the dark energy component. This formalism is more general
and can capture any deviations from $\Lambda$CDM, it is for example also
applicable to modified-gravity models (Amendola et al., 2018). Generally
speaking, the $\mu$ function tells us how much the Poisson equation is
modified as a result of dark energy perturbations or modification of General
Relativity. Within the $\Lambda$CDM framework where dark energy perturbations
are absent, one would obtain the standard Poisson equation, which implies
$\mu(k,z)=1$ for all scales and redshifts.

As already mentioned, at late times and on small scales the sum in Eq. (9) is
dominated by the contribution from matter so that we can simply set $X=m$.
This suggests defining the $\mu$ function as

$\mu(k,z)=1+\frac{\bar{\rho}_{\text{DE}}(z)\Delta_{\text{DE}}(k,z)}{\bar{\rho}_{m}(z)\Delta_{m}(k,z)}\,.$
(10)

However, in actual observations we will in general be considering power
spectra of the metric and matter density perturbations, and for this purpose
we define a function $\mu^{2}$ through the expression

$\mu^{2}(k,z)=\frac{k^{4}{\langle}\Phi\Phi^{*}{\rangle}}{(4\pi
Ga^{2}\bar{\rho}_{m})^{2}\langle\Delta_{m}\Delta_{m}^{*}\rangle}\,.$ (11)

Naively we might think that we can simply replace $\Delta_{\text{DE}}(k,z)$
and $\Delta_{m}(k,z)$ in Eq. (10) with the square root of the dark energy and
matter power spectra. However, inserting the relativistic Poisson equation
(keeping only matter and dark energy)

$-k^{2}\Phi=4\pi Ga^{2}\left(\bar{\rho}_{\rm DE}\Delta_{\rm
DE}+\bar{\rho}_{m}\Delta_{m}\right)$ (12)

into Eq. (11), we obtain

$\displaystyle\mu^{2}(k,z)$ $\displaystyle=\langle\mu(k,z)\mu(k,z)^{*}\rangle$
(13)
$\displaystyle=\frac{\bar{\rho}^{2}_{\text{DE}}\langle\Delta_{\text{DE}}\Delta_{\text{DE}}^{*}\rangle+\bar{\rho}^{2}_{m}\langle\Delta_{m}\Delta^{*}_{m}\rangle+2\bar{\rho}_{\text{DE}}\bar{\rho}_{m}\langle\Delta_{\text{DE}}\Delta^{*}_{m}\rangle}{\bar{\rho}^{2}_{m}\langle\Delta_{m}\Delta^{*}_{m}\rangle}\,.$

If the dark energy and dark matter are fully correlated such that the
normalized cross-correlation coefficient $f_{\times}$, defined as

$f_{\times}=\frac{\langle\Delta_{\text{DE}}\Delta^{*}_{m}\rangle}{\sqrt{\langle\Delta_{\text{DE}}\Delta_{\text{DE}}^{*}\rangle\langle\Delta_{m}\Delta^{*}_{m}\rangle}}$
(14)

is equal to unity on all scales, then the $\mu$ defined in Eq. (10) can be
recovered by taking the square root of the Eq. (13). In linear perturbation
theory this is indeed the case, but on non-linear scales the situation is less
straightforward and the value of the speed of sound can significantly
influence this correlation. We will discuss this in more details in the
following subsection and in Section 7.

###  2.3 Numerical result from the $k$-evolution code

Figure 1:  Upper panels : Comparison of the $\mu$ function for two cases of
$c_{s}^{2}=10^{-7}$ (left panel) and $c_{s}^{2}=10^{-4}$ (right panel) and
$w_{0}=-0.9$, obtained from CLASS, Halofit and k-evolution at $z=0,1,6$. In
both cases, by moving to higher redshifts and smaller scales, the effect of
$k$-essence perturbations decreases; the former is due to the sub-dominant
contribution of dark energy in the past compared to matter component and the
latter is due to the existence of sound horizons. For $c_{s}^{2}=10^{-4}$ at
$z=0$, the sound horizon wavenumber is positioned at $k_{s}\approx
0.029~{}\text{h/Mpc}$, beyond which, dark energy perturbations start to decay.
As this scale is larger than the matter non-linearity threshold, linear
perturbation theory remains sufficient for the prediction of density
contrasts. Consequently, the results from CLASS, Halofit, and k-evolution
align when transitioning from large to small scales. For the case when
$c_{s}^{2}=10^{-7}$ at $z=0$ , the sound horizon is located at $k_{s}\approx
0.94~{}\text{h/Mpc}$ which coincides with the matter non-linearity scale.
After entering this region the result from CLASS separates from those of
Halofit and k-evolution as the latter two account for the non-linear evolution
of matter (while the dark energy perturbation grows slower compared to
matter), leading them to predict lower values for $\mu$. At smaller scales,
when dark energy undergoes non-linear evolution, a deviation is observed
between the results of k-evolution and Halofit. This difference is attributed
to the non-linear evolution of $\Delta_{\text{DE}}$ in k-evolution, which in
turn yields higher values for the $\mu$ function compared to the Halofit.
Lower panels : Comparison of the $\mu$ function in two scenarios:
$c_{s}^{2}=10^{-7}$ (left panel) and $c_{s}^{2}=10^{-4}$ (right panel), with
$w_{0}=-0.9$. The results are derived from $k$-evolution, taking into account
both the actual cross-correlation $f_{\times}(k,z)\leq 1$ and the assumption
that $f_{\times}(k,z)=1$. For $c_{s}^{2}=10^{-4}$, the agreement of $\mu$
functions across all scales in the two mentioned cases is attributed to dark
energy and matter perturbations being well described by the linear theory and
consequently being fully correlated with each other. In contrast, in the case
of $c_{s}^{2}=10^{-7}$, the transition of the $\mu$ in which
$f_{\times}(k,z)\leq 1$ is suppressed earlier compared to the scenario where
$f_{\times}(k,z)=1$, underscoring the significant influence of cross-
correlation on the behavior of $\mu$ for small speeds of sound. For the
k-evolution data, we ran two simulations, each evolving $3840^{3}$ particles
within boxes of comoving lengths $1280~{}\text{Mpc}/h$ and
$9000~{}\text{Mpc}/h$. Figure 2: Normalized cross-correlation power spectra at
$z=0,5,20$, for the cases of $c_{s}^{2}=10^{-4}$, $c_{s}^{2}=10^{-7}$ and
$c_{s}^{2}=10^{-9}$. For the lower speeds of sound, $c_{s}^{2}=10^{-7}$ and
$c_{s}^{2}=10^{-9}$, where dark energy clusters and has its own dynamics,
there is a noticeable decline in cross-correlation coefficient on small
scales. In contrast, at the larger sound speed of $c_{s}^{2}=10^{-4}$, dark
energy and matter densities maintain a full correlation across almost all
scales.

The k-evolution555https://github.com/FarbodHassani/k-evolution code is an
extension of the relativistic $N$-body code
gevolution666https://github.com/gevolution-code/gevolution-1.2 (Adamek et al.,
2016) in which the effective field theory (EFT) approach to the $k$-essence
dark energy has been implemented, see Hassani et al. (2019) for an in-depth
discussion. A distinguishing feature of k-evolution compared to gevolution is
that it fully couples the $k$-essence equation of motion to the non-linear
evolution of matter. In contrast, gevolution employs the linear solution
provided by the Boltzmann code CLASS (Blas et al., 2011) to compute the dark
energy perturbations.

We illustrate the importance of including non-linearities in Fig. 1 where we
show $\mu(k,z)$ function obtained from Eq. (13) for both cases where we
naively assume $f_{\times}=1$ and more realistic case of $f_{\times}\leq 1$.
We compare $\mu(k,z)$ function at different redshifts using CLASS (linear
perturbation theory only), Halofit (the matter power spectrum is non-linear)
(Takahashi et al., 2012) and k-evolution codes. The cosmological parameters
used in these comparisons are listed in Table 1 and include

  * •

$w_{0}$: the dark energy equation-of-state (EoS)

  * •

$h$: the dimensionless Hubble parameter

  * •

$\Omega_{\text{b}}$: fractional energy density of baryonic matter

  * •

$\Omega_{\text{cdm}}$: fractional energy density of cold dark matter

  * •

$A_{s}$: spectral amplitude

  * •

$n_{s}$: spectral index

We examine these parameters for two specific values of the speed of sound
squared, $c_{s}^{2}=10^{-7}$ and $c_{s}^{2}=10^{-4}$.

Table 1: Cosmological parameters used in the computation of $\mu$ in Fig. 1. $w_{0}$ | h | $\Omega_{\text{b}}$ | $\Omega_{\text{cdm}}$ | $A_{s}$ | $n_{s}$
---|---|---|---|---|---
-0.90 | 0.6755 | 0.0486 | 0.2617 | $2.21\times 10^{-9}$ | 0.9660

As discussed e.g. in Sapone et al. (2009), on large scales, outside the sound
horizon, the dark energy clusters proportionally to the dark matter, with a
relative factor that depends primarily on $w_{0}$. Within the sound horizon,
the dark energy perturbations are suppressed. The overall shape that we expect
for $\mu$ in a $k$-essence cosmology is a constant value on large scales, and
then a transition at the sound horizon, with $\mu$ approaching $1$ for large
values of $k$. This overall shape is clearly visible in Fig. 1. We see that
the plateau on large scales has a value of about $1.043$ today for our case of
$w_{0}=-0.9$. At higher redshift, this value scales with the relative
evolution of the dark energy and dark matter densities, $(1+z)^{3w_{0}}$ in
our model where $w_{0}$ is constant.

In the case of $c_{s}^{2}=10^{-4}$ at $z=0$, the sound horizon is located at
$k_{s}\approx 0.029~{}\text{h/Mpc}$, where we can see a clear cut-off in $\mu$
function. Given that this horizon is significantly larger than the matter non-
linearity scale, $k_{\text{nl}}^{m}\approx 0.1~{}\text{h/Mpc}$ at $z=0$, the
transition in $\mu$ is the same for CLASS, Halofit and k-evolution as both
$\Delta_{\text{DE}},\Delta_{m}$ can be perfectly described by linear
perturbation theory. In the case of $c_{s}^{2}=10^{-7}$ however, the sound
horizon is at $k_{s}=0.94~{}\text{h/Mpc}$ at $z=0$, which falls within the
scale where matter undergoes non-linear evolution. Upon entering the scale of
matter non-linearity, $\Delta_{m}$ starts to increase faster than in the
linear solution, while the dark energy density contrast, $\Delta_{\text{DE}}$,
does not increase as significantly. This leads to a decrease in the $\mu$
function. Consequently, in this region, Halofit and k-evolution yield results
that are lower than CLASS as $\Delta_{m}^{\text{nl}}>\Delta_{m}^{\text{lin}}$.
Both Halofit and k-evolution continue to follow the same pattern until they
approach the phase where dark energy begins its non-linear evolution. In this
phase, the predictions of k-evolution start to deviate from those of Halofit,
indicating higher values for $\mu$. This stems from the k-evolution’s ability
to include $k$-essence clustering as it gets drawn into the dark matter
potential. Such clustering leads to an increase in $\Delta_{\text{DE}}$, and
subsequently to an increase in the $\mu$ function
$(\Delta_{\text{DE}}^{\text{nl}}>\Delta_{\text{DE}}^{\text{lin}})$.

So far, the discussion has been about the upper panels of Fig. 1, with the
underlying assumption that dark energy and dark matter maintain a full
correlation. But as mentioned previously, the value of speed of sound can
affect this correlation. In Fig. 2, we show the cross-correlation coefficient,
$f_{\times}$, for different values of the sound speed, namely
$c_{s}^{2}=10^{-4}$, $c_{s}^{2}=10^{-7}$ and $c_{s}^{2}=10^{-9}$. At large
scales, where both matter and dark energy perturbations evolve linearly, a
complete correlation is maintained between them. This full correlation
continues at smaller scales for $c_{s}^{2}=10^{-4}$, where dark energy does
not cluster strongly and tightly tracks the matter density. This is why in the
lower right panel of Fig. 1, in both cases of $f_{\times}(k,z)\leq 1$ and
$f_{\times}(k,z)=1$, the $\mu$ function is the same across all scales and
redshifts. Conversely, for smaller speeds of sound, $c_{s}^{2}=10^{-7}$ and
$c_{s}^{2}=10^{-9}$, due to the reasons outlined earlier regarding dark
energy’s sound horizon being within the region where matter clusters non-
linearly, dark energy is able to cluster strongly and initiate its own
dynamics to some extent. Due to this, as it can be seen from Fig. 2, the
moralized cross-correlation factor start to decay by moving tho small scales.
Moreover, the same logic clarifies why in the lower left panel of Fig. 1, the
$\mu$ function for the case of $f_{\times}(k,z)\leq 1$ undershoots the one
that assumes $f_{\times}=1$, converging towards $\Lambda$CDM at slightly
larger scales.
For illustrative purposes, in Fig. 3 we have shown 2D snapshots of matter and
dark energy densities obtained from k-evolution with a relatively high
resolution; we simulated the evolution of $1200^{3}$ particles inside a
comoving box with a side length of $100~{}\text{Mpc}/h$, resulting in a
spatial resolution of $0.083~{}\text{Mpc}/h$ . When the speed of sound is
small, it allows dark energy to be influenced more significantly by the
gravitational potentials of matter, leading it to cluster in high density
regions. Conversely, for dark energy with a high speed of sound, the rapid
propagation of its perturbations does not allow it to get trapped by the
matter potential in the same way. This high propagation speed essentially
prevent dark energy from clustering. As a result, regions with dense matter do
not see the same accumulation or influence of dark energy as they do when the
speed of sound is low.

Based on these discussions, we can conclude that the k-evolution results are
primarily important for future tests of $k$-essence dark energy featuring
small values of sound speed, an aspect that we will take into account in the
construction of the emulator. We emulate two versions of $\mu$; 1) when we
assume a full correlation between dark matter and dark energy, i.e.
$f_{\times}=1$ in Eq. (13) and 2) when a realistic correlation between dark
matter and dark energy is considered, meaning $f_{\times}\leq 1$. However, we
will only discuss the results and methods for the version where $f_{\times}=1$
is assumed across all scales and redshifts. This is because we use a similar
approach to develop an emulator for the other version. We revisit the
differences between the two in Section 7, where we explore how to use the
emulator to reconstruct the power spectrum of the gravitational potential
$\Phi$ from a matter power spectrum.

Figure 3: Snapshots displaying matter density (middle panel), dark energy
density having a squared sound speed of $c_{s}^{2}=10^{-7}$ overlaid on matter
density (left panel), and dark energy density with a sound speed squared of
$c_{s}^{2}=10^{-4}$ overlaid on matter density (right panel). For low sound
speeds, dark energy is more prone to the gravitational pull of matter, causing
it to cluster in areas of high density. In contrast, when dark energy has a
higher sound speed, its fluctuations spread quickly and it evades the
confining influence of matter potential. This rapid dispersion prevents dark
energy clustering which leads to less pronounced dark energy presence in dense
matter regions compared to when the sound speed is low.

##  3 Emulator Construction: Principal Stages

In this section, we give an overview of the main steps involved in the
construction and configuration of a PCE-based emulator. It should be
emphasized that the approach we use in our work is significantly influenced by
the methods used in the EuclidEmulator papers (Knabenhans et al., 2019, 2021),
given their demonstrated efficiency.

As previously mentioned, we select $\mu$ function to be our emulation
objective, and we define $\epsilon_{\text{emu}}$ as the general emulation
error

$\epsilon_{\text{emu}}(k,z;{\bm{p}})=\frac{\mu_{\text{emulated}}(k,z;{\bm{p}})-\mu_{\text{simulated}}(k,z;{\bm{p}})}{\mu_{\text{simulated}}(k,z;{\bm{p}})},$
(15)

where $k$ and $z$ stand for wavenumber and redshift, and ${\bm{p}}$ represents
a set of cosmological parameter that are used in the computation of $\mu$
function. Note that in the context of EuclidEmulator, the approach taken
involves emulating the boost to the matter power spectrum as opposed to
directly emulating the matter power spectrum itself. In our case, even though
we have the option to emulate the boost of $\mu$ function by separately
calculating $\mu$ from k-evolution and CLASS, we have chosen to directly
emulate $\mu$. This choice has been validated by the minimal emulation errors
encountered in our case and eliminates the requirement for employing a
Boltzmann solver to obtain $\mu$ afterward.

The main goal is to minimize $\epsilon_{\text{emu}}$, a crucial quantity that
determines the accuracy of the final emulator and guides the selection of
certain fixed hyperparameter values involved in the construction of the
emulator. This, however, needs to be achieved while ensuring that the expenses
of building the emulator remain within reasonable limits.

The key stages in the construction of an emulator using the polynomial chaos
expansion (PCE) approach can be outlined as follow:

  1. I)

Sampling cosmological parameters within prior ranges for the purpose of
training the emulator.

  2. II)

Running simulations to compute the $\mu$ function for each sampled parameter
set, and storing these results in a matrix denoted as $\mathcal{D}$.

  3. III)

Applying dimensionality reduction to the matrix $\mathcal{D}$ using principal
component analysis (PCA).

  4. IV)

Creating surrogate models to represent the eigenvalues of the principal
components as functions of the input cosmological parameters, through the use
of SPCE.

  5. V)

Merging the emulated eigenvalues (surrogate models) with their corresponding
principal components.

In order to carry out a comprehensive investigation with the goal of
determining the optimal number of training samples $(n_{\text{ED}})$777In the
domain of uncertainty quantification, the term experimental design (ED) refers
to the ensemble of training input parameters $(X_{\text{ED}})$ and the
resulting outputs $(Y_{\text{ED}})$, which in our case are generated during
the first and second step of the emulator construction. for the final
simulation-based emulator to achieve a desired accuracy, and also to gain an
in-depth insight into the emulation mechanism—including its functionality,
sensitivity to the parameter space’s dimensionality, and process of
hyperparameter optimization— it is practical to start with the Boltzmann code
rather than using $N$-body simulations. Hence, we begin by constructing mock
emulators, following the outlined five steps, while opting to use the Halofit
extension of CLASS. Subsequently, we will apply our findings from this initial
phase to the final simulation-based emulator using k-evolution.

###  3.1 Step I: Prior ranges and sampling strategy

####  3.1.1 Cosmological priors

We define the emulator on top of a seven dimensional cosmological parameters
space, encompassing {$\Omega_{\text{b}}$, $\Omega_{\text{cdm}}$, $n_{s}$, $h$,
$w_{0}$, $\log c_{s}^{2}$ , $A_{s}$}. To ensure consistency with
EuclidEmulator2, for each parameter we have assumed uniform priors within
ranges that are closely aligned with those in EuclidEmulator2, as shown in
Table 2.

Table 2: Parameter ranges for the emulator. Parameter | Range
---|---
$\Omega_{\text{b}}$ | [0.04, 0.06]
$\Omega_{\text{cdm}}$ | [0.20, 0.34]
$n_{s}$ | [0.92, 1]
$h$ | [0.61, 0.73]
$w_{0}$ | [-1.3, -0.7]
$\log c_{s}^{2}$ | [-10, -5]
$A_{s}$ | [$1.7\times 10^{-9},2.5\times 10^{-9}$]

In the process of training the emulator, each cosmology we sample from the
parameter box assumes a fixed value for $\Omega_{\text{rad}}$, corresponding
to the CMB temperature $T_{\text{CMB}}=2.725$ K. Meanwhile, the dark energy
density parameter $\Omega_{\text{DE}}$, is calculated using the flatness
condition $\Omega_{\text{tot}}=1$, expressed as:

$\Omega_{\text{tot}}=\Omega_{\text{DE}}+\Omega_{\text{cdm}}+\Omega_{\text{b}}+\Omega_{\text{rad}}.$
(16)

Differently to the case of EuclidEmulator2, where a dynamical dark energy
model is adopted using both $w_{0}$ and $w_{a}$, our parameterisation for the
$k$-essence model of dark energy employs $w_{0}$ and the speed of sound,
$c_{s}$. For all cosmological models considered in this study, we set $w_{a}$
to zero, resulting in a more simplified yet robust description of the dark
energy. Furthermore we have limited our analysis to small values of the sound
speed between $c_{s}^{2}=10^{-10}$ and $c_{s}^{2}=10^{-5}$ since as showed in
the previous section larger values yield results that are consistent with
those obtained from Boltzmann code.

####  3.1.2 Sampling strategy

The effectiveness of an emulator depends significantly on the strategy used
for sampling from the parameter space. As a widely adopted and popular option,
we choose latin hypercube sampling (LHS) ( e.g. McKay et al., 1979) for
several reasons:

  1. (i)

LHS ensures that each parameter is uniformly sampled across its entire range.
This uniformity is crucial in high-dimensional spaces where the emulator needs
to capture the wide variability of parameters.

  2. (ii)

Compared to simple random sampling, LHS can provide more extensive coverage of
the parameter space with fewer samples. This efficiency is especially valuable
in context of expensive simulations, where the number of training simulations
is restricted by computational resources.

  3. (iii)

The stratified sampling approach of LHS contributes to better predictive
accuracy of the emulator. It allows for a more representative and diverse set
of simulation runs.

We produce Latin hypercube samples of varying sizes,
$n_{\text{ED}}=\\{10,25,50,100,125,200,300\\}$, to subsequently identify the
optimal sample size that meets our required accuracy standards. It is
important to note that before initiating the LHS process, the total number of
samples has to be specified. LHS inherently leads to non-unique outcomes due
to the random nature of the sampling within each stratum of the parameter
space. To address this non-uniqueness in order to optimize our sampling
strategy, we employ the method proposed in Knabenhans et al. (2019, 2021); For
each predetermined sample size, we generate $10^{4}$ distinct LHS sets; the
optimal set is then chosen based on the criterion of maximizing the minimal
Euclidean distance between points, ensuring a sample distribution that is as
spread out and evenly distributed as possible.

###  3.2 Step II: Computation of the $\mu$ function

Following the sampling of cosmological parameter sets, for each set we run
simulations (Halofit in the case of mock emulator) to calculate the $\mu$
function. The results of this computation are then stored in a dataset denoted
as $\mathcal{D}$. This dataset is represented as matrix with $n_{\text{ED}}$
rows and $n_{k}.n_{z}$ columns

$\displaystyle\mathcal{D}$ $\displaystyle=$

It is important to note that redshift is not included in the cosmological
parameters that are sampled. Instead, the matrix $\mathcal{D}$ stores $\mu$
functions at specific redshifts corresponding to various interval steps of the
simulation888In fact, similar to Knabenhans et al. (2019, 2021), we store the
logarithms of these $\mu$ functions inside the matrix $\mathcal{D}$ as doing
so will result in lower emulation error in the end.. Consequently, for any
randomly chosen cosmological parameter from Table 2, the final emulator
produces the $\mu$ function at redshifts and wavenumbers for which it is
trained. To enable the user to obtain the emulated $\mu$ values at redshifts
not present in the matrix, we implement linear interpolation between the two
nearest columns that encompass the requested redshift. Since for now, our
current focus is solely on studying the emulation error, we compute the
spectra at only five specific redshifts: $\\{0,0.5,1,2,3\\}$, across a range
of $1000$ wavenumbers within
$10^{-2}h\,\text{Mpc}^{-1}<k<10h\,\text{Mpc}^{-1}$.

###  3.3 Step III: Principal component analysis (PCA)

Given the high dimensionality of the matrix $\mathcal{D}$, computational
processing is notably challenging. Furthermore, not all details within the
data, such as simulation noises and other non-physical phenomena, are
desirable for our analysis and are thus intended to be removed. In this
context, techniques for reducing dimensionality, such as PCA, are highly
efficient. An additional compelling reason for employing PCA is the
limitations of our emulation method. A basic PCE is inherently designed to
predict scalar outputs and can not be directly applied to non-scalar
quantities like the $\mu$ function. By using PCA, we can decompose the
covariance matrix of the dataset $\mathcal{D}$ to identify the principal
components and extract their associated eigenvalues. These eigenvalues being
scalar quantities, can then be emulated.

PCA decomposes a dataset into orthogonal principal components, where each
component corresponds to an axis along which the variance of the projected
data is maximized. The first axis, known as the first principal component
($\text{PC}_{1}$), captures the largest proportion of the dataset’s variance.
Subsequently, the second principal component ($\text{PC}_{2}$) accounts for
the maximum variance residual from that captured by $\text{PC}_{1}$, under the
constraint of orthogonality to $\text{PC}_{1}$. This process iterates across
all dimensions of the dataset, allowing for the identification of up to $n$
principal components within an $n$-dimensional space, each orthogonal to the
others and sequentially optimizing the variance accounted for. For additional
information on PCA, see e.g Chapter 8 of Geron (2022).

To determine the principal components of the data matrix $\mathcal{D}$, we
employ a matrix factorization method called the Singular Value Decomposition
(SVD) technique. Moreover, PCA presupposes that the dataset is centered; thus,
prior to applying SVD, the dataset $\mathcal{D}$ is centered by subtracting
the mean of each column from its respective elements. This preprocessing step
results in a centered dataset $\mathcal{D}_{\text{centered}}$, which SVD then
decomposes into three matrices as shown below

$\displaystyle\mathcal{D_{\text{centered}}}$
$\displaystyle=\mathbf{U}~{}\mathbf{\Sigma}~{}\mathbf{V^{T}}$ (17)
$\displaystyle=\sum_{i=1}^{n_{\mathrm{ED}}}\lambda_{i}(\Omega_{b},\Omega_{\text{CDM}},n_{s},h,w_{0},\log
c_{s}^{2},A_{s})\mathrm{PC}_{i}(k,z).$

In this decomposition, $\mathbf{V^{T}}$ contains the principal components. The
singular values in $\mathbf{\Sigma}$, when multiplied by $\mathbf{U}$,
correspond to the weights (or magnitudes) of these principal components,
reflecting their contribution to capturing the variance in
$\mathcal{D_{\text{centered}}}$. While there could potentially be
$n_{k}.n_{z}$ principal components for $\mathcal{D}_{\text{centered}}$, the
summation effectively considers only those up to $n_{\mathrm{ED}}$, as the
eigenvalues $\lambda_{i}$ for $i>n_{\mathrm{ED}}$ vanish.

Once all the principal components are determined, the matrix
$\mathcal{D}_{\text{centered}}$ can be projected into a $d$-dimensional space.
This is achieved by considering only the first $d$ columns of $\mathbf{V}$,
which correspond to the first $d$ principal components. But, rather than
selecting the number of principal components arbitrarily, a preferable
strategy is to initially determine $a_{\text{PCA}}$, the total variance we
wish to retain in the projection, and then accordingly choose the appropriate
number of dimensions, $d$, to which we reduce the dimensionality. Following
this approach, the inverse transformation allows for the reconstruction of
$\mathcal{D_{\text{centered}}}$ back into its original dimensional space, as
shown in Eq. (17). However, due to the loss of information during the
projection, this will not be an exact match to $\mathcal{D_{\text{centered}}}$
(unless all the principal components are considered in the summation). The
reconstruction error can be measured by computing the mean squared distance
between the two datasets, i.e. between $\mathcal{D_{\text{centered}}}$ and its
reconstructed counterpart.

We generate a sequence of $a_{\text{PCA}}$ values, including
$a_{\text{PCA}}=\\{0.9,0.99,\ldots,0.999999\\}$, for the purpose of
determining the most suitable value through a grid search, which will be
explained in detail later.

###  3.4 Step IV: Polynomial chaos expansion (PCE)

Emulating the scalar outputs $\lambda_{i}$, obtained from the PCA in the
previous step, constitutes the essential phase of this work. Since each
$\lambda_{i}$ captures a distinct aspect of variance in the data, it is
crucial to emulate them separately. To do this, we employ PCE method
integrated within the UQLab999https://www.uqlab.com software to provide a
functional approximation of each eigenvalue. PCE achieves this by representing
the model’s spectral characteristics on an appropriately constructed
orthonormal basis of polynomial functions. In our case, the expansion is
expressed as follows

$\lambda_{i}\left(\Omega_{b},\Omega_{\text{CDM}},n_{s},h,w_{0},\log
c_{s}^{2},A_{s}\right)=\sum_{\bm{\alpha}\in\mathbb{N}^{M}}\beta_{\bm{\alpha}}\Psi_{\bm{\alpha}}({\bf
X}).$ (18)

Here, $\bm{\alpha}=\\{\alpha_{1},\alpha_{2},...,\alpha_{M}\\}$ represents a
multi-index, $\Psi_{\bm{\alpha}}$ corresponds to the orthonormal multivariate
polynomial basis, and $\beta_{\bm{\alpha}}$ is its associated coefficient.
Moreover, ${\bf X}=\\{X_{1},X_{2},...,X_{M}\\}^{\text{T}}$ denotes the input
random vector of the cosmological parameters, and $M$ indicates the
dimensionality of the parameter space, which in our case is seven. The
multivariate polynomials $\Psi_{\bm{\alpha}}$, are constructed through the
product of univariate polynomials $\phi_{\bm{\alpha}}^{(j)}(X_{j})$, with the
specific type of each univariate polynomial determined by the distribution of
the respective input variables $X_{i}$. Given the uniform distribution of
input cosmologies in our context, each constructed
$\phi_{\bm{\alpha}}^{(j)}(X_{j})$ is a member of the Legendre polynomial
family,

$\Psi_{\bm{\alpha}}(\bm{X})\equiv\prod_{j=1}^{M}\phi_{\alpha}^{(j)}(X_{j})=\prod_{l=1}^{M}\sqrt{2\alpha_{l}+1}P_{\alpha_{l}}(X_{l}).$
(19)

####  3.4.1 Sparse polynomial chaos expansion: Truncation Strategy

For practical implementation of the Eq. (18), the summation has to be
truncated to fit within the computational limitations. A common strategy to
achieve this involves restricting the number of multi-indices which results in
the reduction of the number of multivariate polynomial bases
$\Psi_{\bm{\alpha}}$ participating in the summation. This restriction is
implemented by constructing a limited index set (also known as the truncation
set) denoted as $\mathcal{A}$, such that $\mathcal{A}\subset\mathbb{N}^{M}$.

A logical initial step in constructing $\mathcal{A}$ would be to filter the
polynomial bases in order to keep only those whose total degree does not
exceed $p$. This leads to the definition of an index set characterized as
follow:

$\mathcal{A}^{M,p}=\\{\bm{\alpha}\in\mathbb{N}^{M},||\bm{\alpha}||_{1}\leq
p\\},$ (20)

where the term $||\bm{\alpha}||_{1}$ represents the length of the vector
$\bm{\alpha}$, and is defined as

$||\bm{\alpha}||_{1}=\sum_{j=1}^{M}\alpha_{j}.$ (21)

However, this truncation scheme solely may not sufficiently reduce the
complexity in practical situations, since the total number of polynomial
terms, denoted by $N$, follows the formula

$N=\begin{pmatrix}p+M\\\ p\end{pmatrix}=\frac{(p+M)!}{p!\,M!}.$ (22)

Consequently, the number of terms grow rapidly with the increase in the number
of input variables, making the process computationally challenging and often
impractical in scenarios involving high-dimensional parameter spaces.

To address this, we turn to strategies rooted in the sparsity-of-effects
principle. This principle is based on the observation that in many complex
systems, only a few interaction terms (low-order interactions) significantly
influence the behaviour of systems. Most of the higher-order interactions or
terms have negligible effects. Therefore, identifying and focusing on these
lower order influential terms can lead to a more efficient and computationally
manageable model. For that, we employ two widely used truncation scheme as
detailed in Blatman (2009); Blatman & Sudret (2011); Marelli et al. (2022a).

One such approach is known as hyperbolic truncation. This strategy modifies
(extends) the standard degree-based truncation scheme to maintain significant
lower-order terms, further penalizing higher-rank indices. Unlike the standard
method, which includes all terms up to a specified degree, hyperbolic
truncation concentrates on terms with a lower combined total degree,
determined by the $q$-norms condition. This corresponds to an index set
defined as

$\mathcal{A}^{M,p,q}=\\{\bm{\alpha}\in\mathcal{A}^{M,p},||\bm{\alpha}||_{q}\leq
p\\},$ (23)

where

$||\bm{\alpha}||_{q}=\left(\sum_{j=1}^{M}{\alpha_{j}}^{q}\right)^{1/q}~{}~{}~{},~{}~{}~{}0<q\leq
1.$ (24)

Conceptually, the hyperbolic condition can be seen as the selection of
polynomial terms that fall under a hyperbolic curve within the space defined
by the multi-index $\bm{\alpha}$ (see Fig. 1 of Marelli et al., 2022a).

Another complementary strategy is the maximal interaction. This technique
limits the highest order of interactions between variables in the constructed
polynomials, rather than just limiting the degree of individual polynomials.
The model can be specifically adjusted to limit the number of non-zero
elements in the multi-index $\bm{\alpha}$ to a maximum value $r$, where $r\leq
M$. This is mathematically represented as

$\mathcal{A}^{M,p,r}=\\{\bm{\alpha}\in\mathcal{A}^{M,p},||\bm{\alpha}||_{0}\leq
r\\},$ (25)

with

$||\bm{\alpha}||_{0}=\sum_{j=1}^{M}1_{\\{\alpha_{j}>0\\}},$ (26)

being defined as the rank of $\bm{\alpha}$. In practical terms, this
translates to allowing only $r$ or fewer univariate polynomials
$\phi_{\alpha}^{(i)}(X_{i})$ to participate in constructing each term of the
multivariate polynomial bases. Essentially, it restricts the model to consider
interactions among at most $r$ input parameters in each $\Psi_{\bm{\alpha}}$.

Integrating both truncation methods, the definitive index set can be
formulated as

$\mathcal{A}^{M,p,q,r}=\\{\bm{\alpha}\in\mathcal{A}^{M,p},||\bm{\alpha}||_{q}\leq
p~{}~{}{\text{and}}~{}~{}||\bm{\alpha}||_{0}\leq r\\}.$ (27)

Thus the SPCE of each considered eigenvalues reads

$\lambda_{i}\approx\sum_{\bm{\alpha}\in\mathcal{A}_{i}^{M,p,q,r}}\beta_{\bm{\alpha}}\Psi_{\bm{\alpha}}({\bf
X}).$ (28)

####  3.4.2 Calculation of the coefficients

In Blatman (2009), several methods for calculating the expansion coefficients
$\beta_{\alpha}$ are proposed. However, UQLab implements only those methods
that involve post-processing of model evaluations. The method we employ here
is the Least Angle Regression (LARS) technique (Efron et al., 2004). This
approach solves the following optimization problem

$\hat{\bm{\beta}}_{i,\alpha}=\underset{\bm{\beta}_{i}\in\mathbb{R}^{\left|\mathcal{A}_{i}^{p,q,r}\right|}}{\operatorname{argmin~{}}}\mathbb{E}\left[\left(\bm{\beta}_{i}^{\top}\Psi_{i}(\bm{X})-\lambda_{i,\text{true}}(\bm{X})\right)^{2}+\gamma\|\bm{\beta}\|_{1}\right].$
(29)

Here, the loss function is defined as the mean squared error between the
emulated and true eigenvalue of each principal component. Additionally, it
includes a regularization term

$\gamma\|\bm{\beta}\|_{1}=\gamma\sum\limits_{\alpha\in\mathcal{A}_{i}^{p,q,r}}|\beta_{i,\alpha}|,$
(30)

which is introduced to penalize high-rank solutions in the minimization
process. Further details on the calculation of PCE coefficients based on LARS
are available in Marelli et al. (2022a); Efron et al. (2004).

####  3.4.3 Basis adaptivity of PCE

Determining the most suitable finite basis for an effective PCE modeling is
often demanding in practical contexts. The basis adaptivity algorithm in PCE
addresses this challenge by selecting the most relevant polynomial terms from
a broad range of candidates. The process involves an iterative approach, where
the degree $p$ and the $q$-norm parameters are adjusted to optimize the
validation error of the emulator. Based on predefined set of $p$ and $q$
values, the algorithm executes a sequence of operations that can be outlined
as follow:

  1. (i)

Starting the process by constructing a basis with initial values $p=p_{0}$ and
$q=q_{0}$.

  2. (ii)

Computing the PCE coefficients and measuring the generalization error of the
model. This measurement involves using the Leave-One-Out (LOO) cross-
validation error, $\varepsilon_{\text{LOO}}$, which is used to prevent
overfitting. It is calculated by constructing $n_{\text{ED}}$ separate
emulated eigenvalues, $\lambda^{\backslash i}_{\text{emu}}$, with each being
built upon a reduced experimental design $X\backslash
x(i)=\\{x(j),j=1,\ldots,n_{\text{ED}},j\neq i\\}$, where $\backslash x(i)$
denotes the exclusion of the $i$-th data point from $X$ . The method
sequentially leaves out each data point, trains the model on the remaining
data (on $n_{\text{ED}}-1$ samples), and then assesses its accuracy by
comparing the model’s predictions at the excluded data point $x(i)$ with the
actual value $\lambda_{\text{true}}$. The LOO error is thus given by

$\varepsilon_{\text{LOO}}=\frac{\displaystyle\sum\limits_{i=1}^{n_{\text{ED}}}\left(\lambda_{\text{true}}(x(i))-\lambda^{\backslash
i}_{\text{emu}}(x(i))\right)^{2}}{\displaystyle\sum\limits_{i=1}^{n_{\text{ED}}}\left(\lambda_{\text{true}}(x(i))-\hat{\mu}_{\lambda}\right)^{2}},$
(31)

where the denominator serves to normalize the error by scaling it against the
variance of the true output. Note that The computation of
$\varepsilon_{\text{LOO}}$ is carried out entirely within the experimental
design, which includes both the training dataset and the corresponding
responses, and does not involve a separate test set.

  3. (iii)

Evaluating the calculated error against a predetermined threshold. If the
error is within this threshold or if the algorithm reaches its maximum
iteration limit, it terminates and selects the PCE model with the minimal
$\varepsilon_{\text{LOO}}$. If not, the algorithm adjusts the parameters by
incrementing $p$ or $q$, and restarts the process from the basis generation
phase.

The default setting in UQLab is such that if $\epsilon_{\text{LOO}}$ does not
decrease over two consecutive iterations, the degree-adaptive scheme stops
increasing the maximum degree. For the case of q-norm adaptivity, when the
polynomial degree is low, the effect of the q-norm is minimal and usually does
not lead to the inclusion of extra basis functions, resulting in the
$\epsilon_{\text{LOO}}$ remaining unchanged. Because of that, in UQLab,
iterations of $q$-norm are only considered relevant if the PCE’s basis size or
the $\epsilon_{\text{LOO}}$ are affected. Subsequently, the increase of
$q$-norm is automatically stopped if there is no decrease in the
$\epsilon_{\text{LOO}}$ over two relevant iterations. In our configuration,
however, we deactivated this early stop mechanism and allowed for a complete
examination of all $q$-norm candidates for every degree. Doing so, we ensured
that the optimal candidate is chosen after all possibilities have been
considered.

###  3.5 Step V: Combining emulated $\lambda_{i}$ and principal components

In the last phase, the emulated $\lambda_{i}$ are merged with the
corresponding principal components that we have taken into account in the SVD
process. It is important to recall that we centered the $\mathcal{D}$ matrix
before conducting SVD by subtracting the mean of each column from the
respective $\mu$ function. Therefore, to formulate the emulated matrix we must
re-add these mean values to the combined set of emulated eigenvalues and
principal components. This leads to the final derivation of the emulated $\mu$
as shown in the following equation

$\mu_{\text{emu}}(k,z;\bm{p})\approx\mu_{\mathrm{PCA}}(k,z)+\sum_{i=1}^{n_{\mathrm{PCA}}}\sum_{\alpha\in\mathcal{A}_{i}^{p,q,r}}\hat{\beta}_{i,\alpha}\Psi_{i}^{\alpha}(\bm{p})\mathrm{PC}_{i}(k,z).$
(32)

##  4 Fine-Tuning emulator’s parameters

In the preceding section, we have seen how the construction of the emulator at
various stages of its assembly is fundamentally influenced by a series of
critical parameters, namely the Latin hypercube sample size ($n_{\text{ED}}$),
the percentage of the dataset’s variance we wish to retain during PCA
($a_{\text{PCA}}$), the maximum total degree of the polynomial bases ($p$),
and the truncation parameters ($q$ and $r$). So, it becomes clear that
choosing the best values for these parameters is crucial for ensuring that the
emulator operates at its highest level of performance. To this end, we employ
a grid search analysis. In a grid search we explore in details a manually
specified subset of these parameters101010Here, there is no need to specify
grid values for $r$; UQLAB automatically determines the optimal $r$ during
polynomial coefficient computation. For instance, if $r$ is found to be $3$,
it implies that polynomials with interactions beyond $3$ do not significantly
contribute to the model, as their coefficients would be negligible or zero..
The objective is to evaluate different combinations of parameter values to
determine which combination yields the best performance in terms of the
emulator’s accuracy and computational efficiency. For this purpose, we have
created the following subset of parameters

$\displaystyle n_{\text{ED}}$ $\displaystyle=\\{10,25,50,100,125,200,300\\},$
(33) $\displaystyle a_{\text{PCA}}$
$\displaystyle=\\{0.9,0.99,0.999,...,0.999999\\},$ $\displaystyle p$
$\displaystyle=\\{3,4,5,...,20\\},$ $\displaystyle q$
$\displaystyle=\\{0.25,0.30,0.35,...,1\\}.$

The fine-tuning of the emulator’s parameters is strategically divided into two
steps. In the first step, the hyperparameters $n_{\text{ED}}$ and
$a_{\text{PCA}}$ are adjusted based on their performance against a separate
test set. Since the largest emulation errors are often found near the edges of
the parameter spaces, we have generated $10^{6}$ random test samples and then
filtered out $15\,000$ samples that fall inside a hypersphere of radius 1,
inscribed in the normalized parameter box, where each parameter is mapped
between $-1$ and $1$. This approach effectively concentrates our validation on
the more reliable central region while minimizing the impact of edge-related
errors. Note that in addition to the cosmological parameters listed in Table
2, the test set that we generate also includes the sum of neutrino masses in
the range $0.0~{}\text{eV}\leq\sum m_{\nu}\leq 0.15~{}\text{eV}$, thus making
the test sets effectively eight-dimensional. It is important to clarify that
the emulator itself is not trained to include the sum of neutrino masses as a
variable. However, as we will discuss later, the inclusion of this parameter
has a minimal impact on the $\mu$ function, indicating that the emulator is
not particularly sensitive to it. By adding this parameter to the test set,
users are given the option to include the sum of neutrino masses in their
analyses for a range of applications. However, it is worth mentioning that the
effect of including the sum of neutrino masses is only introduced through the
change in either $\Omega_{\rm m}$ or $\Omega_{\rm DE}$, as massive neutrinos
contribute to the total energy density with $\Omega_{\rm tot}=1$.

In the second step, the remaining parameters, $p$, $q$, and $r$, are fine-
tuned internally using the LOO error measurement, which, as discussed before,
evaluates the cross-validation error exclusively based on the experimental
design, i.e. the training datasets and their corresponding outputs.

Similar to the approach taken in Knabenhans et al. (2019, 2021), we conducted
an analysis, illustrated in Fig. 4, to examine how the maximum emulation error
varies with the number of training samples and the fractional variance
retained during PCA.

Figure 4: Demonstration of the maximum emulation error (y-axis) as a function
of the number of training samples, $n_{\text{ED}}$ (x-axis), and the retained
fractional variance, $a_{\text{PCA}}$, indicated by different line colors. The
graph shows the dual impact of the training sample size and the variance
retained during PCA on the accuracy of the emulator tested against $15\times
10^{3}$ test samples. For each set of $n_{\text{ED}}$ and $a_{\text{PCA}}$,
the parameters $p$, $q$ and $r$ are selected based on the LOO error. Figure 5:
Emulation error distribution across $1000$ wavenumbers at $z=0$, evaluated on
$15\times 10^{3}$ eight-dimensional test sets, at the 99th percentile (black
area), 95th percentile (red area), and 68th percentile (yellow area). These
percentiles represent the values below which $99\%$, $95\%$, and $68\%$ of
data points, respectively, are included.

In Fig. 4, our analysis shows that even when up to $90\%$ of the total
variance is preserved in PCA, the maximum emulation error does not fall below
$1\%$, despite using as many as $300$ training samples. This outcome,
highlighted by the plateaus in the graph’s curves for different values of
$a_{\text{PCA}}$, implies that a mere increase in training sample size, beyond
a certain threshold, is ineffective in further reducing the maximal emulation
error when the retained fractional variance during PCA is insufficient.
Based on this analysis, we have decided to use $200$ training samples in the
construction of the k-evolution-based emulator, ensuring that
$a_{\text{PCA}}=99.999\%$ of the total variance is preserved through PCA. In
the case of the mock emulator, this amount of $a_{\text{PCA}}$ has
corresponded to $11$ principal components, resulting in $11$ separate PCE for
emulating the associated eigenvalues. Details regarding the k-evolution-based
emulator will be provided in Section 6.

In Fig. 5, the percentile representation of the mock emulator, constructed
with these $11$ principal components is depicted. The accuracy was evaluated
using the previously mentioned $15\,000$ test samples. We also conducted an
in-depth performance analysis of the mock emulator within the mentioned eight-
dimensional parameter space, constructing 28 distinct 2D planes by pairing
different parameters. On each plane, we created a $40\times 40$ grid of test
sample points, resulting in a total of $44\,800$ samples across all the
planes, and then tested the emulator’s performance at each point while keeping
the other parameters, not involved in the respective plane, constant. Doing
so, we could analyze the effects of specific parameter interactions on the
emulator’s accuracy and study the emulator’s performance across various areas
of the parameter spaces. To visualize the results, we generated heat maps for
each of the 28 planes, as shown in Fig. 14.

##  5 Sensitivity analysis based on Sobol indices

To assess the influence of individual input parameters on the eigenvalues
$\lambda_{i}$, we use the Hoeffding-Sobol decomposition technique, as detailed
in Sobol (2001); Hoeffding (1948). This method employs variance-based
sensitivity analysis, where Sobol indices (denoted by $S$) are used to
distribute the variance of $\lambda_{i}$ across each input parameter, and
across combinations of parameters including pairs, triplets, and higher-order
interactions. For the computation of Sobol indices, we adopt a polynomial
expansion approach (PCE-based Sobol indices, as referenced in Sudret (2008);
Marelli et al. (2022b)) instead of relying on traditional Monte Carlo
simulations.

Figure 6: First-order Sobol indices for the first four principal components.

More description of this approach, including its theoretical framework, can be
found in Appendix A.

Fig. 6 shows the result of the first-order Sobol indices ($S_{1}$) for the
first four principal components, based on $200$ simulations. Clearly, $w_{0}$
has the most impact on the output model, as it holds the highest first-order
Sobol index in the first principal component. It is important to recognize
that a significant portion of the total variance in the original data, about
$97.77\%$, is already captured in the projection onto the first principal
component. Consequently, the significance of the first order Sobol indices for
the cosmological parameter diminishes generally as one progresses to higher
principal components. Following the same logic, $c_{s}^{2}$ and
$\Omega_{\text{cdm}}$ are identified as the subsequent influential
cosmological parameters due to their secondary contribution to the second and
third principal components. These can be justified by considering that $w_{0}$
and $\Omega_{\text{cdm}}$ collectively influence the amplitude of the $\mu$
function, whereas $c_{s}^{2}$ affects its cut-off scale, transitioning from
large to small scales.
While Sobol analysis effectively reveals the relative influence of
cosmological parameters on the $\mu$ function, it particularly underscores the
significant influence of the $k$-essence parameters, namely $w_{0}$ and
$c_{s}^{2}$, as somewhat expected. To better investigate the impacts of other
cosmological parameters, beyond those associated with $k$-essence, we keep the
$k$-essence parameters constant, setting $w_{0}=-0.9$ and $c_{s}^{2}=10^{-7}$.
We then proceed by once again conducting a targeted PCE-based Sobol analysis
on the remaining parameters, including the sum of neutrino masses, to
investigate its impact as well. The result of this is presented in Fig. 7. As
can be observed from the analysis, among all the cosmological parameters, the
sum of the neutrino masses exhibits the lowest impact across all the first
four principal components. Due to the negligible influence of the total
neutrino masses, we decided to omit it from the training set. However, users
have the option to enter a value ranging from $0.0~{}\text{eV}$ to
$0.15~{}\text{eV}$ in the emulator’s input for their own purposes (see the
explanation in Section 4). This is why in Section 4, we added $\sum m_{\nu}$
into the test set to demonstrate explicitly that this does not lead to
unacceptable emulation errors.

Figure 7: First-order Sobol indices for the first four principal components.
$w_{0}$ and $c_{s}^{2}$ are maintained as fixed values (and hence not involved
in the PCE) to explore the relative impacts of other parameters.

##  6 Simulation-based emulator

Up to this point, our main focus has been on the configuration and the
construction of the mock emulator using Halofit within the CLASS framework.
The main reason for this choice was its lower computational demand compared to
$N$-body simulations, which made it possible to generate extensive test
datasets—a task practically impossible with $N$-body simulations. In this
section, we intend to apply the insights gained from developing the mock
emulator to construct an emulator based on the k-evolution code.

###  6.1 Convergence tests

Before proceeding with simulations to train the emulator, conducting
convergence tests is imperative. These tests aim to determine the optimal
simulation parameters, such as box size, number of particles, and
corresponding spatial resolution, necessary to achieve an acceptable level of
accuracy. In our case, we are interested in results as accurate as possible
within the $k$ value range of
$10^{-2}h\,\text{Mpc}^{-1}<k<10h\,\text{Mpc}^{-1}$.

Figure 8: Convergence test analysis for the $\mu$ function, with mesh sizes of
$N_{\text{grid}}=\\{380^{3},512^{3},1024^{3},1200^{3},1600^{3},1984^{3}\\}$
where $N_{\text{pcl}}=N_{\text{grid}}$, demonstrates that a grid size of
$N_{\text{grid}}=1200^{3}$ reaches a convergence level of $0.1\%$ within the
wavenumber range of $10^{-2}h\,\text{Mpc}^{-1}$ to $10h\,\text{Mpc}^{-1}$.
Increasing the grid size to $N_{\text{grid}}=1600^{3}$ and
$N_{\text{grid}}=1984^{3}$, results in only marginal improvement in the
precision of convergence.

For the Flagship simulation, we conducted two simulations based on the
cosmological parameters specified in Table 1. The first simulation involved
$3840^{3}$ particles and a box size of $L=1280~{}h^{-1}~{}\text{Mpc}$, leading
to a Nyquist wavenumber of $k_{\text{N}}=9.42~{}h~{}\text{Mpc}^{-1}$. The
second simulation contained $2304^{3}$ particles and a smaller box size of
$L=300~{}h^{-1}~{}\text{Mpc}$, yielding a Nyquist wavenumber of
$k_{\text{N}}=24.12~{}h~{}\text{Mpc}^{-1}$. By linking these two simulation we
managed to provide reliable results on both large and small scales. In our
series of test simulations, we primarily focused on varying the grid size
($N_{\text{grid}}$) while maintaining a consistent box size for the majority
of the tests. The details of these simulations are shown in Table 3. For each
test, the results were truncated at the corresponding Nyquist wavenumber to
prevent sampling errors and preserve the validity of the simulation data. In
Fig. 8 we have shown the result of these convergence tests.

Table 3: Simulation parameters and computational resources

Label | $N_{\text{grid}}$ |  $L$ ($h^{-1}$ Mpc) |  $k_{\text{N}}$ ($h~{}\text{Mpc}^{-1}$) |  Runtime (core hours) | Machine
---|---|---|---|---|---
1FS | $3840^{3}$ | 1280 | 9.42 |  50688 | CSCS
2FS | $2304^{3}$ | 300 | 24.13 |  129024 | CSCS
Sim1 | $380^{3}$ | 400 | 2.98 |  19 | Baobab
Sim2 | $512^{3}$ | 400 | 4.02 |  64 | Baobab
Sim3 | $1024^{3}$ | 400 | 8.04 |  1536 | Baobab
Sim4 | $1200^{3}$ | 400 | 9.42 |  2560 | Baobab
Sim5 | $1600^{3}$ | 400 | 12.56 |  9627 | Baobab
Sim6 | $1984^{3}$ | 620 | 10.05 |  17641 | Baobab

Our analysis reveals that the test simulation, using
$N_{\text{grid}}=N_{\text{pcl}}=1200^{3}$ and a box length of
$L=400h~{}\text{Mpc}^{-1}$, has successfully attained a satisfactory level of
convergence, falling below $0.1\%$ up to $k=9.42h~{}\text{Mpc}^{-1}$. Further
increase in number of particles and the size of the simulation volume appear
to have minimal impact on improving this convergence. Given the satisfactory
results achieved with the settings $N_{\text{grid}}=1200^{3}$ and
$L=400h\,\text{Mpc}^{-1}$, we have decided to select these specific parameters
to generate the training set for the simulation-based emulator.

We should highlight that the acceptable convergence level achieved with this
relatively small number of particles, owes much to the way of calculating the
$\mu$ function by dividing the square root of one power spectrum by another.
This division significantly diminishes the resolution effects and the impact
of cosmic variance, the latter being a consequence of the limited size of the
simulation box that leads to statistical uncertainties on large scales. Such
reductions in these non-physical effects are the main reason behind our
decision to emulate the $\mu$ function instead of directly emulating the
potential power spectrum. The high susceptibility of the potential power
spectrum to such issues implies that maintaining an acceptable level of
accuracy across a wide spectrum of wavenumbers would require a substantially
larger number of particles and greater simulation volume. By opting for the
$\mu$ function instead, we mitigate the necessity for extensive computational
resources.

Figure 9: Demonstration of the impact of SG filter: employing a 14-data point
window and a polynomial of order 3 leads to a notable RMSE reduction,
specifically in the wavenumber range from $k=0.025h\,\text{Mpc}^{-1}$ to
$k=1h\,\text{Mpc}^{-1}$.

To further diminish the impact of cosmic variance on the $\mu$ function, we
use the Savitzky-Golay (SG) filter (Savitzky & Golay, 1964). This filter
essentially reduces signal noise while preserving its characteristics through
the application of polynomial fits within localized segments defined by moving
windows. By using an SG filter with a window length of $14$ data points and a
polynomial order of $3$, i.e. fitting third-order polynomials to $14$ points
at each step, we successfully achieved a notable reduction in noise across our
data, especially at larger scales as depicted in Fig. 9. By doing so, we could
lower down the root mean square error (RMSE), described as

$\text{RMSE}=\sqrt{\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}\left(\frac{\mu_{i}-\hat{\mu}_{i}}{\mu_{i}}\right)^{2}},$
(34)

from $0.064\%$ to $0.043\%$ across the $k$ range of
$0.025h\,\text{Mpc}^{-1}<k<1h\,\text{Mpc}^{-1}$.

###  6.2 From mock to simulation-based emulator

Following the analysis presented in Fig. 4, we proceed with the same set of
$200$ Latin hypercube-sampled cosmologies for the purpose of training the
k-evolution-based emulator. These $200$ simulations are executed on the Baobab
and Yggdrasil clusters111111https://doc.eresearch.unige.ch/hpc/start. As
previously discussed, all simulations maintain the fixed parameters of
$N_{\text{grid}}=1200$ and $L=400h\,\text{Mpc}^{-1}$, and for every cosmology
set, the resulting $\mu$ function is processed through an SG filter with a
window length of $14$ and a polynomial order of $3$. Here, spectra are
calculated at $53$ distinct $z$-values121212The selection of these redshifts
is designed to ensure that the $\mu$ values across different redshifts remain
as close as possible, especially at larger scales, to avoid any significant
gap. Doing so helps in minimizing errors when performing linear interpolation
between neighboring redshifts. within the range of $0<z<3$ and at 1024
distinct $k$-bins. After the truncation of the results beyond the Nyquist
wavenumber, this results in 591 $k$-values ranging from
$10^{-2}h\,\text{Mpc}^{-1}$ to $9.42h\,\text{Mpc}^{-1}$. However, prior to
restoring the data in a matrix format similar to that discussed in Section
3.2, we perform interpolation on the smoothed $\mu$ function using the
available $k$-values. Subsequently, we calculate the $\mu$ values at $1000$
evenly spaced wavenumbers within the mentioned range. Therefore, the final
matrix has $n_{\text{ED}}=200$ rows and $n_{k}.n_{z}=53\,000$ columns. The
remaining stages of the emulator’s construction are the same as explained
before in Section 3. This includes the principal component analysis of the
matrix $\mathcal{D}$, polynomial chaos expansion of the scalar outputs
$\lambda_{i}$, and finally re-merging the emulated eigenvalues with the
corresponding principal components. Since these steps are identical to those
described in the construction of the mock emulator, we will not repeat them in
detail here. Instead, we will focus on presenting the the results of these
analyses.

Table 4: Results of the polynomial chaos expansions for k-evolution
experimental design. The SFB (size of full basis) column indicates the total
number of polynomial terms generated after truncation, using $q$ parameter,
while the SSB (size of sparse basis) column reflects the number of terms with
non-zero coefficients after employing the LARS method. Polynomial chaos output
---
PC |  p |  q |  r |  $\epsilon_{\text{LOO}}$ |  SFB |  SSB
1 |  6 |  0.85 |  4 |  $1.03\times 10^{-5}$ |  533 |  82
2 |  10 |  0.60 |  3 |  $7.39\times 10^{-4}$ |  673 |  73
3 |  8 |  0.60 |  3 |  $3.01\times 10^{-4}$ |  302 |  94
4 |  10 |  0.65 |  4 |  $5.2\times 10^{-3}$ |  960 |  48
5 |  10 |  0.75 |  4 |  $5.80\times 10^{-3}$ |  2304 |  88
6 |  6 |  0.85 |  4 |  $2.55\times 10^{-2}$ |  533 |  53
7 |  12 |  0.55 |  3 |  $4.04\times 10^{-2}$ |  771 |  59
8 |  10 |  0.75 |  4 |  $1.30\times 10^{-2}$ |  2304 |  60
9 |  14 |  0.60 |  4 |  $2.24\times 10^{-2}$ |  2150 |  53

By retaining $99.999\%$ of the total variance of the dataset $\mathcal{D}$
when projecting it onto a lower-dimensional hyperplane during PCA, we identify
$9$ principal components. Additional details and plots are available in
Appendix C. Consequently, we create $9$ PCEs to emulate the weights these
identified components. Each PCE undergoes an optimization process (basis
adaptivity), with the parameters $p$ and $q$ being individually fine-tuned.
The outcomes of this fine-tuning, including the LOO validation error
($\epsilon_{\text{LOO}}$) and the sizes of the full and sparse basis (denoted
as SFB and SSB, respectively) are detailed in Table 4. In the absence of
applying the hyperbolic truncation parameter $q$, the size of the full
polynomial basis for any given degree $p$ is calculated according to Eq. (22).
For instance, considering the polynomial expansion for the first principal
component with $p=6$ and given the input dimensionality of $M=7$, this
calculation yields a basis size of $1716$ terms. Nevertheless, as clearly
demonstrated in Table 4, the intervention of $q$ parameter has significantly
reduced the basis size to $533$ terms. Subsequent use of LARS method for
calculating the coefficients further reduced this number, yielding only $82$
of these terms having non-zero coefficients, with each of these polynomials
depending on at most 4 variables, as indicated by $r=4$. The graphical
representation of non-zero coefficient for all the principal components are
illustrated in Appendix C.

Figure 10: First-order Sobol indices for the first four principal components.
The results are largely consistent with the mock emulator case, with a notable
exception being the prominence of $c_{s}^{2}$ over $w_{0}$ in the second
principal component, reflecting the specific impact of non-linear clustering
of dark energy in k-evolution model.

We also conducted a first order Sobol index analysis for the first four
principal components, with the findings illustrated in Fig. 10. The outcomes
closely match with those from the mock emulator scenario (Fig. 6), except for
the distinct prominence of $c_{s}^{2}$ over $w_{0}$ in the second principal
component. This difference can be attributed to the stronger influence of
$c_{s}^{2}$ within k-evolution framework due to the non-linear clustering of
dark energy.

###  6.3 Emulator’s performance

Due to computational constraints, we confined the evaluation of the emulator’s
performance to $20$ test sets of cosmological parameters. These sets are
randomly selected from a larger collection of $15\,000$ test sets previously
generated during the configuration of the mock emulator. They are then
subjected to k-evolution simulation, employing the same number of particles
and volume size as those used in the training process. In Fig. 11, we show the
average emulation error over the $20$ test sets at each wavenumber for
different redshifts (upper panel), alongside the emulation errors for all $20$
test set cosmologies at $z=0$ (lower panel). From this group, we identified
the cosmology corresponding to the highest emulation error (approximately
$0.08\%$) and compare the simulated and emulated $\mu$ at different redshifts.
The comparative results are depicted in Fig. 12.

In terms of speed, the emulator—implemented within a Python wrapper—takes
approximately $0.2$ seconds ($5.36\times 10^{-5}$ core hours) per evaluation
on a regular laptop, compared to the $2.56\times 10^{3}$ core hours required
for a simulation at the same resolution. This corresponds to a speed-up factor
of nearly fifty million.

Figure 11: Upper panel: Average emulation error across $20$ cosmologies at
different wavenumbers, observed at redshifts $z=0,0.5,1,2$.
Lower panel: Emulation error at $z=0$ for each set of 20 cosmologies (depicted
in grey lines) together with the average emulation error at each $k$ (shown by
the black line) and the standard deviation (represented in green). The maximum
emulation error recorded is $\approx 0.08$ percent.  Figure 12: Comparison
between simulated and emulated values of $\mu$ across various redshifts for
the test set exhibiting the highest emulation error, approximately $0.08\%$ at
$z=0$.

##  7 Reconstructing the potential power spectrum

Figure 13: Reconstruction of the potential power spectrum using Eq. (35)
through the use of the $\mu$ function and the matter power spectrum at $z=0$
(left panel) and $z=0.3$ (right panel). Employing both the simulated and
emulated $\mu$ function in which the suppression in cross-correlation between
dark matter and dark energy at small scales is considered (respectively
denoted as $\mu^{\text{simulated}}_{f_{\times}\leq 1}$ and
$\mu^{\text{emulated}}_{f_{\times}\leq 1}$), results in consistent outcomes
compared to the direct output of the potential power spectrum. This
consistency is observed within the wavenumber range unaffected by cosmic
variance and resolution effects at both $z=0$ and $z=0.3$. At scales affected
by these factor, the results surpass the $1\%$ threshold. Conversely, when
using $\mu^{\text{simulated}}_{f_{\times}=1}$, i.e. assuming a full
correlation between dark matter and dark energy, the deviation from the direct
output reaches a maximum of $\sim 2.5\%$ at $z=0$, while exhibiting almost
$1\%$ maximum deviation at $z=0.3$ due to the less effective clustering of
dark energy.

In this section, we explain the process of determining the potential power
spectrum by using the $\mu$ function together with the matter power spectrum,
showing how accurate estimations of both can enable a precise calculation of
the potential power spectrum.

To derive the potential power spectrum, one may rearrange Eq. (11) such that

$k\Delta^{2}_{\Phi}=\frac{\mu(k,z)^{2}(4\pi
G_{N}a^{2}\bar{\rho}_{m})^{2}\Delta^{2}_{m}}{k^{3}},$ (35)

where $\Delta^{2}_{\Phi}$ and $\Delta^{2}_{m}$ are respectively dimensionless
potential and matter power spectrum given by

$\Delta_{X}^{2}(k)=\frac{k^{3}}{2\pi^{2}}P_{X}(k).$ (36)

As demonstrated in Fig. 13, efforts to recreate the true potential power
spectrum based on the cosmological parameters of Table 1, reveal that using
the CLASS-based $\mu$ function ($\mu_{\text{CLASS}}$) results in a maximum
information loss of more than $5\%$ at $z=0$ compared to a direct calculation
of the true potential power spectrum obtained with k-evolution code. Even
applying a simulated $\mu$ function that disregards the true cross-correlation
between dark energy and dark matter, $\mu^{\text{simulated}}_{f_{\times}=1}$,
a maximum deviation of approximately $2\%$ at the same redshift is observed.
In contrast, consistent reproduction of $\Phi$ power spectrum is only achieved
when employing the $\mu$ function that accounts for the true cross-
correlation, i.e. $\mu^{\text{simulated}}_{f_{\times}\leq 1}$. At $z=0.3$,
both $\mu_{\text{CLASS}}$ and $\mu^{\text{simulated}}_{f_{\times}=1}$
demonstrate less deviations compared to the scenario at $z=0$. This is
because, at higher redshifts, the less pronounced clustering of dark energy
results in a stronger correlation with dark matter.
All the results shown in Fig. 13, including the computation of $\mu$ function
and the direct computation of the potential power spectrum, are derived from
the combination of the two mentioned Flagship simulations in Section 6.1. It
is evident that reproducing the potential power spectrum, using the high-
resolution $\mu^{\text{simulated}}_{f_{\times}\leq 1}$ in one instance and the
emulated $\mu^{\text{emulated}}_{f_{\times}\leq 1}$ in another, has produced
nearly identical results. This supports our choice to opt for the simulation
settings $N_{\text{grid}}=1200^{3}$ and $L=400~{}h^{-1}~{}\text{Mpc}$ in the
convergence test investigation and indicates that an exact predicted value of
$\mu$ function (using our emulator) and the matter power spectrum (using, for
example, EuclidEmulator) will result in a precise reproduction of the
potential power spectrum through the Eq. (35).

The reproduction error using the emulated $\mu$ function remains well below
the $1\%$ limit across the wavenumber range unaffected by cosmic variance and
resolution effects, i.e. from $k=0.02h\,\text{Mpc}^{-1}$ to
$k=6h\,\text{Mpc}^{-1}$. However, it exceeds the $1\%$ limit on small scales
due to resolution constraints and on large scales due to cosmic variance.
These deviations might be attributed to how differently cosmic variance and
resolution effects are incorporated into the potential and dark energy spectra
by different solvers. To further investigate this issue, we conducted
simulations at varying resolutions and calculated the $\mu$ function using Eq.
(11). We then compared these results within their reliable range –not skewed
by cosmic variance and resolution effects– to the $\mu$ output from a single
fixed-resolution simulation obtained from Eq. (13). By observing that the
results from the variable-resolution simulations using Eq. (13) match with the
one obtained from Eq. (13), we then concluded that using both the matter power
spectrum and $\mu$ function obtained from Eq. (10), gives a more accurate
estimation of potential power spectrum compared to the direct output provided
by the simulation.

##  8 Conclusions

We develop an emulator to accurately model the non-linear clustering effects
of the $k$-essence dark energy, using the PCE approach implemented in the
UQLab software. The emulator is built on a seven dimensional parameter space
and is based on k-evolution simulations using $200$ training samples. We use a
mock (Halofit-based) emulator to study the optimal emulator configuration, and
find a maximum emulation error of $0.15\%$ using only $11$ principal
components. This was assessed against the outcome of $15\,000$ sets of
cosmological parameters that were randomly sampled within a hypersphere of
radius 1, inscribed in the normalized parameter space. By applying the
insights gained from the construction of the mock emulator, we find a $0.08\%$
maximum emulation error in the case of the k-evolution-based emulator
(projected onto $9$ principal components) which is assessed against only $20$
test sets due to the limited computational budget. The emulator is executed in
a wall time of $\sim 0.2s$ on a usual laptop which is comparable to the case
of EuclidEmulator1 ($\sim 0.4s$ per evaluation) and EuclidEmulator2 ($\sim
0.3s$ per evaluation) in which the same emulation method is used.

Using the polynomial expansions obtained for the emulator, we perform a
sensitivity analysis based on Sobol indices, which in turn reveals the
dominant influence of $w_{0}$ and $c_{s}^{2}$ on $\mu(k,z)$. Through the same
analysis, we also observed that the sum of neutrino masses have the least
impact on the $\mu$ function among other cosmological parameters.

Our emulator predicts the $\mu$ function that encodes the amount of clustering
in dark energy relative to the dark matter. Since it is based on the ratio of
two power spectra, the $\mu$ function suffers much less from cosmic variance
and resolution effects than the power spectra themselves. This is a crucial
advantage, since it significantly reduces the requirements for particle
numbers and simulation volume for a given accuracy.

One important application of our emulator is the accurate reconstruction of
the lensing potential power spectrum from a matter power spectrum. To obtain
the correct power spectrum on non-linear scales, one needs to account for the
true cross-correlation between dark energy and dark matter in the definition
of $\mu$. Overlooking this aspect leads to a non-negligible deviation from the
true potential power spectrum at low redshifts. This is particularly critical
for the lensing signal, where errors accumulate due to its integrating effect.
For this reason, we built two emulators, one that includes the true cross-
correlations between dark energy and dark matter and one that assumes a full
correlation between the two at all scales and redshifts. We demonstrate
explicitly that the former, i.e. $\mu^{\text{emulated}}_{f_{\times}\leq 1}$ ,
allows to recover the $\Phi$ power spectrum to high accuracy from the matter
power spectrum.

The emulator is available for public use at
https://github.com/anourizo/k-emulator. The repository includes the data for
polynomial bases and their respective coefficients for both variants of the
$\mu$ function, with $f_{\times}=1$ and $f_{\times}\leq 1$. Additionally, it
includes the Python wrapper and a DEMO Jupyter notebook that guide users
through the process of using the emulator to generate emulated $\mu$ function.

## Acknowledgements

ARN and MK acknowledge funding from the Swiss National Science Foundation. FH
acknowledges the Research Council of Norway and the resources provided by
UNINETT Sigma2 – the National Infrastructure for High Performance Computing
and Data Storage in Norway. FH would also like to thank Elif Ilhan, Álvaro de
la Cruz-Dombriz, Peter Dunsby, Roy Maartens and James Adam for their support
and hospitality during the research visit to Cape Town. ARN would like to
thank the Vahabzadeh Foundation for their continuous financial support. This
work is also supported by a grant from the Swiss National Supercomputing
Centre (CSCS) under project ID s1051. Most of the computations were performed
at the University of Geneva using Baobab and Yggdrasil HPC services.

## Data Availability

The emulator code and data is publicly available at
https://github.com/anourizo/k-emulator .

## References

  * Abbott et al. (2018) Abbott T. M. C., et al., 2018, Phys. Rev. D, 98, 043526
  * Abell et al. (2009) Abell P. A., et al., 2009, LSST Science Book, Version 2.0 (arXiv:0912.0201)
  * Adamek et al. (2016) Adamek J., Daverio D., Durrer R., Kunz M., 2016, JCAP, 07, 053
  * Ade et al. (2016) Ade P. A. R., et al., 2016, Astron. Astrophys., 594, A14
  * Adler et al. (1995) Adler R. J., Casey B., Jacob O. C., 1995, Am. J. Phys., 63, 620
  * Aghamousa et al. (2016) Aghamousa A., et al., 2016, The DESI Experiment Part I: Science,Targeting, and Survey Design (arXiv:1611.00036)
  * Aghanim et al. (2016) Aghanim N., et al., 2016, Astron. Astrophys., 594, A11
  * Aghanim et al. (2020) Aghanim N., et al., 2020, Astron. Astrophys., 641, A6
  * Akeson et al. (2019) Akeson R., et al., 2019, The Wide Field Infrared Survey Telescope: 100 Hubbles for the 2020s (arXiv:1902.05569)
  * Alam et al. (2017) Alam S., et al., 2017, Mon. Not. Roy. Astron. Soc., 470, 2617
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Rev. Rel., 21, 2
  * Armendariz-Picon et al. (1999) Armendariz-Picon C., Damour T., Mukhanov V. F., 1999, Phys. Lett. B, 458, 209
  * Armendariz-Picon et al. (2000) Armendariz-Picon C., Mukhanov V. F., Steinhardt P. J., 2000, Phys. Rev. Lett., 85, 4438
  * Asgari et al. (2021) Asgari M., et al., 2021, Astron. Astrophys., 645, A104
  * Blas et al. (2011) Blas D., Lesgourgues J., Tram T., 2011, JCAP, 07, 034
  * Blatman (2009) Blatman G., 2009, PhD thesis, Université Blaise Pascal, Clermont-Ferrand, France
  * Blatman & Sudret (2011) Blatman G., Sudret B., 2011, Journal of Computational Physics, 230, 2345
  * Carrasco et al. (2014) Carrasco J. J. M., Foreman S., Green D., Senatore L., 2014, JCAP, 07, 057
  * Chiba et al. (2000) Chiba T., Okabe T., Yamaguchi M., 2000, Phys. Rev. D, 62, 023511
  * DeRose et al. (2019) DeRose J., et al., 2019, Astrophys. J., 875, 69
  * Efron et al. (2004) Efron B., Hastie T., Johnstone I., Tibshirani R., 2004, The Annals of Statistics, 32, 407
  * Geron (2022) Geron A., 2022, Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems, 2nd edn. O’Reilly Media
  * Hassani et al. (2019) Hassani F., Adamek J., Kunz M., Vernizzi F., 2019, JCAP, 12, 011
  * Hassani et al. (2020) Hassani F., L’Huillier B., Shafieloo A., Kunz M., Adamek J., 2020, JCAP, 04, 039
  * Heitmann et al. (2009) Heitmann K., Higdon D., White M., Habib S., Williams B. J., Wagner C., 2009, Astrophys. J., 705, 156
  * Heitmann et al. (2010) Heitmann K., White M., Wagner C., Habib S., Higdon D., 2010, Astrophys. J., 715, 104
  * Heitmann et al. (2014) Heitmann K., Lawrence E., Kwan J., Habib S., Higdon D., 2014, Astrophys. J., 780, 111
  * Heitmann et al. (2016) Heitmann K., et al., 2016, Astrophys. J., 820, 108
  * Hildebrandt et al. (2017) Hildebrandt H., et al., 2017, Mon. Not. Roy. Astron. Soc., 465, 1454
  * Hildebrandt et al. (2020) Hildebrandt H., et al., 2020, Astron. Astrophys., 633, A69
  * Hoeffding (1948) Hoeffding W., 1948, Ann. Math. Stat
  * Knabenhans et al. (2019) Knabenhans M., et al., 2019, Mon. Not. Roy. Astron. Soc., 484, 5509
  * Knabenhans et al. (2021) Knabenhans M., et al., 2021, Mon. Not. Roy. Astron. Soc., 505, 2840
  * Laureijs et al. (2011) Laureijs R., et al., 2011, Euclid Definition Study Report (arXiv:1110.3193)
  * Lawrence et al. (2010) Lawrence E., Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010, Astrophys. J., 713, 1322
  * Lawrence et al. (2017) Lawrence E., et al., 2017, Astrophys. J., 847, 50
  * Marelli & Sudret (2014) Marelli S., Sudret B., 2014, in Proceedings of the 2nd International Conference on Vulnerability and Risk Analysis and Management (ICVRAM 2014). University of Liverpool, United Kingdom, pp 2554–2563, doi:10.1061/9780784413609.257
  * Marelli et al. (2022a) Marelli S., Lüthen N., Sudret B., 2022a, QLab user manual – Polynomial chaos expansions, Report UQLab-V2.0-104, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2022. https://www.uqlab.com/pce-user-manual
  * Marelli et al. (2022b) Marelli S., Lamas C., Konakli K., Mylonas C., Wiederkehr P., Sudret B., 2022b, UQLab user manual – Sensitivity analysis, Report UQLab-V2.0-106, Chair of Risk, Safety and Uncertainty Quan- tification, ETH Zurich, Switzerland, 2022. https://www.uqlab.com/sensitivity-user-manual
  * Martin (2012) Martin J., 2012, Comptes Rendus. Physique, 13, 566
  * McClintock et al. (2019) McClintock T., et al., 2019, Astrophys. J., 872, 53
  * McKay et al. (1979) McKay M. D., Beckman R. J., Conover W. J., 1979, Technometrics, 21, 239
  * Moran et al. (2023) Moran K. R., et al., 2023, Mon. Not. Roy. Astron. Soc., 520, 3443
  * Nishimichi et al. (2019) Nishimichi T., et al., 2019, Astrophys. J., 884, 29
  * Piras & Spurio Mancini (2023) Piras D., Spurio Mancini A., 2023, The Open Journal of Astrophysics, 6, 20
  * Riess et al. (2021) Riess A. G., Casertano S., Yuan W., Bowers J. B., Macri L., Zinn J. C., Scolnic D., 2021, Astrophys. J. Lett., 908, L6
  * Sapone et al. (2009) Sapone D., Kunz M., Kunz M., 2009, Phys. Rev. D, 80, 083519
  * Savitzky & Golay (1964) Savitzky A., Golay M. J. E., 1964, Analytical Chemistry, 36, 1627
  * Smith & Angulo (2019) Smith R. E., Angulo R. E., 2019, Mon. Not. Roy. Astron. Soc., 486, 1448
  * Sobol (2001) Sobol I. M., 2001, Math Comput Simul
  * Spurio Mancini et al. (2022) Spurio Mancini A., Piras D., Alsing J., Joachimi B., Hobson M. P., 2022, Mon. Not. Roy. Astron. Soc., 511, 1771
  * Storey-Fisher et al. (2024) Storey-Fisher K., Tinker J. L., Zhai Z., DeRose J., Wechsler R. H., Banerjee A., 2024, Astrophys. J., 961, 208
  * Sudret (2008) Sudret B., 2008, Reliability Engineering & System Safety, 93, 964
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, Astrophys. J., 761, 152
  * Valcin et al. (2019) Valcin D., Villaescusa-Navarro F., Verde L., Raccanelli A., 2019, JCAP, 12, 057
  * Zhai et al. (2019) Zhai Z., et al., 2019, Astrophys. J., 874, 95

##  Appendix A Analysis of variance (ANOVA) based on polynomial chaos
expansion

###  A.1 Hoeffding-Sobol decomposition

We define $g({\bm{x}})$ as a square-integrable function representing a
computational model $Y$, where ${\bm{x}}=\\{x_{1},x_{2},\ldots,x_{M}\\}$ is an
$M$-dimensional vector of independent variables within the support
$\mathcal{D}_{\bm{x}}$ and probability density function $f$, given by

$\mathcal{D}_{\bm{x}}:=\\{0\leq x_{i}\leq 1|i\in
1,\ldots,M;M\in\mathbb{N}\\},\quad f({\bm{x}})=1.$ (37)

The goal is to analyze how sensitive $g({\bm{x}})$ is to each input parameter,
their pairs, triplets, and so on. According to Sobol (2001) , the function
$g(\bm{x})$ can be decomposed into components of increasing dimensions,
expressed as follow:

$\displaystyle g({\bm{x}})=g_{0}$
$\displaystyle+\sum_{i=1}^{M}g_{i}\left(x_{i}\right)+\sum_{1\leq i<j\leq
M}g_{ij}\left(x_{i},x_{j}\right)$ $\displaystyle+\sum_{1\leq
i_{1}<\ldots<i_{s}\leq M}g_{i_{1}\ldots
i_{s}}\left(x_{i_{1}},\ldots,x_{i_{s}}\right)+\ldots+g_{1,\ldots,M}\left({\bm{x}}\right).$
(38)

In the above expression $g_{0}$ is a constant term and is equal to the
expected value of $g(\bm{x})$. It is also assumed that the integrals of the
summands $g_{i_{1}\ldots i_{s}}$ with respect to any of their own variables
vanish

$\int_{0}^{1}g_{i_{1},\ldots,i_{s}}\left(x_{i_{1}},\ldots,x_{i_{s}}\right)dx_{i_{k}}=0,~{}1\leq
k\leq s,~{}1\leq i_{1}\leq...i_{s}\leq M,$ (39)

which also yields to the fact that the summands in the expression (38) are
mutually orthogonal in the integration space, and can be calculated in a
recursive approach as proposed in Marelli et al. (2022b):

$\displaystyle g_{i}\left(x_{i}\right)$
$\displaystyle=\int_{0}^{1}\cdots\int_{0}^{1}g(\bm{x})d\bm{x}_{\sim
i}-g_{0}=\mathbb{E}(Y|x_{i})-g_{0},$ $\displaystyle
g_{ij}\left(x_{i},x_{j}\right)$
$\displaystyle=\int_{0}^{1}\cdots\int_{0}^{1}g(\bm{x})d\bm{x}_{\sim(ij)}-g_{0}-g_{i}\left(x_{i}\right)-g_{j}\left(x_{j}\right)$
$\displaystyle=\mathbb{E}(Y|x_{i},x_{j})-g_{0}-g_{i}-g_{j}.$

where $\sim$ stands for the excluded variables, i.e.

$d\bm{x}_{\sim i}=(dx_{1},dx_{2},...,dx_{i-1},dx_{i+1},...,dx_{M-1},dx_{M}).$
(40)

The orthogonality of these components, as indicated by this approach, also
leads to the unique nature of the decomposition (38).

###  A.2 Sobol indices

Building upon the preceding results, the total variance of $Y$ is expressed as

$\text{D}=\text{Var}{[Y]}=\int_{\mathcal{D}_{\bm{x}}}g^{2}(\bm{x})d\bm{x}-g_{0}^{2}.$
(41)

Integrating the square of the expression (38) over $\mathcal{D}_{\bm{x}}$
allows the decomposition of the total variance D into partial variances

$\text{D}=\sum_{i=1}^{M}\text{D}_{i}+\sum_{1\leqslant i<j\leqslant
M}\text{D}_{ij}+\cdots+\text{D}_{1,2,\ldots,M},$ (42)

where the partial variances are defined as

$\text{D}_{i_{1},\ldots,i_{s}}=\int_{0}^{1}\ldots\int_{0}^{1}g_{i_{1},\ldots,i_{s}}^{2}\left(x_{i_{1}},\ldots,x_{i_{s}}\right)dx_{i_{1}}\ldots
dx_{i_{s}}.$ (43)

The Sobol indices are then introduced as

$S_{i_{1},\ldots,i_{s}}=\frac{\text{D}_{i_{1},\ldots,i_{s}}}{\text{D}},$ (44)

and according to the equation (42), they follow the constraint

$\sum_{i=1}^{M}S_{i}+\sum_{1\leqslant i<j\leqslant
M}S_{ij}+\cdots+S_{1,2,\ldots,M}=1.$ (45)

The first-order Sobol indices $S_{i}$ (main effects) characterize the
fractional contribution of each individual input variable to the total
variance, while the second-order Sobol indices $S_{ij}$, where $i\neq j$,
indicate the fractional contribution of each pair of input variables to the
overall output variance. The interpretation extends similarly to $M$-th order
Sobol indices.

###  A.3 PCE-based Sobol indices

Typically, calculating Sobol indices involves conducting Monte Carlo
simulations. Nonetheless, the Eq. (38) contains a total number of terms
represented by

$\sum_{l=1}^{M}\begin{pmatrix}M\\\ l\end{pmatrix}=2^{M}-1.$ (46)

This implies that a minimum of $2^{M}$ Monte Carlo integrals is necessary for
estimating these indices. However, this approach becomes impractical for
models with a large number of dimensions ($M$). To resolve this issue, we can
turn to PCE for calculating Sobol indices, an approach suggested by Sudret
(2008).

Based on the framework proposed in Marelli et al. (2022b), we can develop the
concept of PCE-based sensitivity analysis as follows: introducing the set
$\mathcal{U}=\\{1,...,M\\}$ and its subset $\bm{u}\subset\mathcal{U}$, the
Hoeffding-Sobol decomposition (38) can be rewritten as follow

$g({\bm{x}})=g_{0}+\sum_{\bm{u}\subset\mathcal{U}}g_{\bm{u}}({\bm{x}}_{\bm{u}}).$
(47)

In this equation, $g_{\bm{u}}$ denotes the partial effect of the parameters
within the subvector $\bm{x}_{\bm{u}}=(x_{i})_{i\in\bm{u}}$.

On the other hand, the truncated polynomial expansion, Eq. (28), can be
rearranged to form an expansion in terms of summands of increasing order. This
is done by imposing a new condition on the finite index set
$\mathcal{A}^{M,p,r}_{q}$, such that for a non-empty set
$\bm{v}\subset\mathcal{U}$, $\mathcal{A}_{\bm{v}}$ is defined as:

$\mathcal{A}_{\bm{v}}:=\left\\{{\bm{\alpha}}\in\mathcal{A}^{M,p,r}_{q},\begin{pmatrix}\alpha_{k}>0,\forall
k=1,\ldots,M\mid k\in{\bm{v}}\\\ \alpha_{k}=0,\forall k=1,\ldots,M\mid
k\notin{\bm{v}}\end{pmatrix}\right\\}.$ (48)

Subsequently, the rearranged truncated PCE is represented as,

$Y=\beta_{0}+\sum\limits_{\bm{v}\subset\mathcal{U}}\sum\limits_{\bm{\alpha}\in\mathcal{A}_{\bm{v}}}\beta_{\alpha}\Psi_{\alpha}(\bm{x}).$
(49)

Since the Hoeffding-Sobol decomposition is unique, comparing Eqs. (47) and
(49) results in

$g_{\bm{u}}(\bm{x}_{\bm{u}})=\sum\limits_{\bm{\alpha}\in\mathcal{A}_{\bm{v}}}\beta_{\alpha}\Psi_{\alpha}(\bm{x}).$
(50)

Given the orthonormality of multivariate polynomial basis, both the partial
and total variances can be attributed to the squares of the PCE coefficients

$\displaystyle\text{Var}[Y]$
$\displaystyle=\sum_{\bm{\alpha}\in\mathcal{A}^{M,p,r}_{q}}\beta_{\bm{\alpha}}^{2}$
(51) $\displaystyle\text{Var}[g_{\bm{u}}(\bm{x}_{\bm{u}})]$
$\displaystyle=\sum_{\bm{\alpha}\in\mathcal{A}_{\bm{v}}}\beta_{\bm{\alpha}}^{2},$

Consequently, the corresponding Sobol indices can be represented as follow:

$S_{\bm{v}}=\frac{\text{Var}[g_{\bm{u}}(\bm{x}_{\bm{u}})]}{\text{Var}[Y]}=\frac{\sum\limits_{\bm{\alpha}\in\mathcal{A}_{\bm{v}}}\beta_{\bm{\alpha}}^{2}}{\sum\limits_{\bm{\alpha}\in\mathcal{A}^{M,p,r}_{q}}\beta_{\bm{\alpha}}^{2}}.$
(52)

##  Appendix B Relative emulation error in 2D planes

In Fig. 14 we show the maximum emulation error of the mock emulator across
$28$ different 2D planes, with each plane containing $1\,600$ test samples. We
see that the most significant deviations appear for extreme choices (corners)
of the dark-energy related parameters $w_{0}$ and $c_{s}^{2}$ as well as
$\Omega_{\rm cdm}$ which indirectly controls $\Omega_{\rm DE}$. It should be
mentioned that, even though these are the most significant deviations, they
still correspond to low emulation errors.

Figure 14: Heat maps displaying the emulator’s performance across 28 distinct
2D planes, each containing $1600$ grid test samples.

##  Appendix C Principal components of the response data

In Fig. 15 we plot the mean of the dataset $\mathcal{D}$ alongside the nine
principal components, which are derived by preserving $99.999\%$ of the
variance in $\mathcal{D}$. These components are presented as functions of the
wavenumbers at $z=0$. Fig. 16 illustrates the magnitude of the PCE
coefficients, represented as $\log_{10}(|\beta_{\alpha}|)$, plotted against
their enumeration index ($\alpha$) for all PCE expansions of the weights of
principal components considered in constructing the k-evolution-based
emulator.

Figure 15: Illustration of the mean of response data and eigenvectors obtained
from the covariance matrix of the response data. These principal components,
after being scaled by their corresponding eigenvalues generated through
emulation, are aggregated with the mean. This combination forms the emulated
$\mu$ function.

Figure 16: Figures of the log-scaled spectra of PCE coefficients obtained
through sparse regression using the LARS method for each principal component.

##

##